2017-03-08 81 views
1
我遇到了麻煩頻譜出傅里葉變換

...我有一些數據: dataFFT功率譜愁楚

,我平均爲中心,並且似乎並不有太多多的趨勢...

我繪製傅里葉變換它:

fourier transform

而我得到的東西,是不是很好....

這裏是我的代碼:

def fourier_spectrum(X, sample_freq=1): 
    ps = np.abs(np.fft.fft(X))**2 
    freqs = np.fft.fftfreq(X.size, sample_freq) 
    idx = np.argsort(freqs) 

    plt.plot(freqs[idx], ps[idx]) 

由於從here採取代碼調整。

這似乎對一些幼稚的正弦波數據的工作:

fourier_spectrum(np.sin(2*np.pi*np.linspace(-10,10,400)), 20./400) 

sin spectrum

所以我的問題是:我期待一個非零幾乎隨處可見光譜,我是什麼做錯了?如果我沒有做錯什麼,我的數據的哪些特徵導致了這種情況?另外,如果我沒有做錯任何事情,並且由於某些原因,fft僅僅不適合我的數據,那麼我應該如何從我的數據中提取重要頻率?

+0

只是不要繪製fft的0 bin,如果有太多「DC」組件。同樣通過眼球,數據看起來在非常低的頻率上有很強的頻率分量 - 通過查看前50個分檔可能看起來更「合理」 – f5r5e5d

+0

在FFT之前減去信號的平均值 - 在「DC 「(即,FFT的倉位0)。另外/替代地,繪製FFT輸出的「log10(abs(...))」(即dB或分貝)以更好地理解發生了什麼。 –

+0

另請參見統計譜估計器,例如[periodogram](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.periodogram.html)和'scipy.signal'中的相關方法。 –

回答

0

事實證明,我只是不明白頻譜中的x軸單位,即Hz。因爲我的樣本間隔大約在一秒左右,而且我的時間在一天左右,所以在我的頻譜中真正可見的唯一單位是〜1/s(邊緣處)到約1/m(近中間)以及比這更長的週期的任何東西都與0無法區分。我的誤解源於this教程中的圖表,他們在那裏進行轉換,以便x軸單位在時間上,而不是反時間。我重寫了我的frequency_spectrum繪圖功能,以在結果圖上做適當的「縮放」...

def fourier_spectrum(X, sample_spacing_in_s=1, min_period_in_s=5): 
    ''' 
     X: is our data 
     sample_spacing_in_s: is the time spacing between samples 
     min_period_in_s: is the minimum period we want to show up in our 
      graph... this is handy because if our sample spacing is 
      small compared to the periods in our data, then our spikes 
      will all cluster near 0 (the infinite period) and we can't 
      see them. E.g. if you want to see periods on the order of 
      days, set min_period_in_s=5*60*60 #5 hours 
    ''' 
    ps = np.abs(np.fft.fft(X))**2 
    freqs = np.fft.fftfreq(X.size, sample_spacing_in_s) 
    idx = np.argsort(freqs) 
    plt.plot(freqs[idx], ps[idx]) 
    plt.xlim(-1./min_period_in_s,1./min_period_in_s) # the x-axis is in Hz