2017-10-18 8 views
1

この件に関するお手伝いをお待ちしております。私は最近、ノイズが含まれている場合、離散フーリエ変換のためのParsevalの定理を解くことを試みてきました。私のコードはthis codeに基づいています。Parsevalの定理は正弦波+ノイズのFFTのために保持されませんか?

ノイズが含まれていない場合と同様に、周波数領域のパワーの合計は、負の周波数をカットしているため、時間領域のパワーの半分になります。

しかし、時間領域信号に多くのノイズが加えられると、信号+ノイズのフーリエ変換の総電力は、信号+ノイズの合計電力の半分よりもずっと少なくなります。次のように

私のコードは次のとおりです。

import numpy as np 
import numpy.fft as nf 
import matplotlib.pyplot as plt 

def findingdifference(randomvalues): 

    n    = int(1e7) #number of points 
    tmax   = 40e-3 #measurement time 
    f1    = 30e6 #beat frequency 

    t    = np.linspace(-tmax,tmax,num=n) #define time axis 
    dt    = t[1]-t[0] #time spacing 

    gt    = np.sin(2*np.pi*f1*t)+randomvalues #make a sin + noise 

    fftfreq   = nf.fftfreq(n,dt) #defining frequency (x) axis 
    hkk    = nf.fft(gt) # fourier transform of sinusoid + noise 
    hkn    = nf.fft(randomvalues) #fourier transform of just noise 

    fftfreq   = fftfreq[fftfreq>0] #only taking positive frequencies 
    hkk    = hkk[fftfreq>0] 
    hkn    = hkn[fftfreq>0] 

    timedomain_p = sum(abs(gt)**2.0)*dt #parseval's theorem for time 
    freqdomain_p = sum(abs(hkk)**2.0)*dt/n # parseval's therom for frequency 
    difference  = (timedomain_p-freqdomain_p)/timedomain_p*100 #percentage diff 

    tdomain_pn = sum(abs(randomvalues)**2.0)*dt #parseval's for time 
    fdomain_pn = sum(abs(hkn)**2.0)*dt/n # parseval's for frequency 
    difference_n = (tdomain_pn-fdomain_pn)/tdomain_pn*100 #percent diff 

    return difference,difference_n 

def definingvalues(max_amp,length): 

    noise_amplitude  = np.linspace(0,max_amp,length) #defining noise amplitude 
    difference   = np.zeros((2,len(noise_amplitude))) 
    randomvals   = np.random.random(int(1e7)) #defining noise 

    for i in range(len(noise_amplitude)): 
     difference[:,i] = (findingdifference(noise_amplitude[i]*randomvals)) 

    return noise_amplitude,difference 

def figure(max_amp,length): 

    noise_amplitude,difference = definingvalues(max_amp,length) 

    plt.figure() 
    plt.plot(noise_amplitude,difference[0,:],color='red') 
    plt.plot(noise_amplitude,difference[1,:],color='blue') 
    plt.xlabel('Noise_Variable') 
    plt.ylabel(r'Difference in $\%$') 
    plt.show() 

    return 
figure(max_amp=3,length=21) 

私の最後のグラフは、このfigureようになります。これを解決するときに私は何か間違っていますか?ノイズが増えてこの傾向が生じるという物理的な理由はありますか?完全に正弦波ではない信号に対してフーリエ変換を行うこととは関係がありますか?私がこれをやっている理由は、私が実際のデータを持っている非常にノイズの多い正弦波信号を理解することです。

答えて

1

電力を計算するためにスペクトル全体(正の負)を使用すると、Parsevalの定理が一般的に成立します。

不一致の理由は、DC(f = 0)コンポーネントで、これはやや特殊な扱いです。

まず、DCコンポーネントはどこから来ていますか? np.random.randomを使用して、0と1の間のランダムな値を生成します。平均的に、信号を0.5 * noise_amplitudeだけ上げると、多くの電力が必要になります。このパワーは時間領域で正しく計算されます。

しかし、周波数領域では、f = 0に対応する単一のFFTビンのみが存在します。他のすべての周波数の電力は2つのビンに分配され、1つのビンにはDC電力のみが含まれます。

ノイズをスケーリングすると、DC電源が追加されます。負の周波数を除去することによって信号パワーの半分を除去しますが、ほとんどのノイズパワーは完全に使用されるDCコンポーネントにあります。

あなたは、いくつかのオプションがあります:

  1. 力を計算するために、すべての周波数を使用します。直流成分を含まない
  2. 使用ノイズ:2番目はOKかもしれない、私はオプション1で行くだろうhkk[fftfreq==0] /= np.sqrt(2)

とI:randomvals = np.random.random(int(1e7)) - 0.5

  • 直流電力の半分を除去することにより、電力計算を「修正します」

    fftfreq   = fftfreq[fftfreq>0] #only taking positive frequencies 
    hkk    = hkk[fftfreq>0] 
    hkn    = hkn[fftfreq>0] 
    

    これは本当に意味がありません。本当に3.

    最後に、コードにマイナーな問題があるのはお勧めしません。それをより良いものに変更するには、

    hkk    = hkk[fftfreq>=0] 
    hkn    = hkn[fftfreq>=0] 
    

    に変更するか、オプション1で完全に削除してください。

  • +0

    ありがとうございました!これは今動作します。私は第一と第二の選択肢を試してみました。 –

    関連する問題