FFT頻譜分析原理與python的實現

本篇內容介紹了“FFT頻譜分析原理與python的實現”的有關知識，在實際案例的操作過程中，不少人都會遇到這樣的困境，接下來就讓小編帶領大家學習一下如何處理這些情況吧！希望大家仔細閱讀，能夠學有所成！

點擊上面"腦機接口社區"關注我們

頻譜分析

下面是一組用于描述和解釋信號屬性的常用量(matlab的常見形式，python中的常見形式也類似)：

x:  采樣的數據；

n=length(x):   樣本數量；

fs:   采樣頻率(每單位時間或空間的樣本數)(單位常用:赫茲Hz)；

dt=1/fs   ：每樣本的時間或空間增量(如果是時間上的增量，則又稱：采樣間隔或采樣步長，單位常用:s)；

t=(0:n-1)/fs ： 數據的時間或空間范圍；

y=fft(x) ： 數據的離散傅里葉變換(DFT)；

abs(y) ：DFT的振幅；

(abs(y).^2)/n ：DFT的冪；

fs/n  ： 頻率增量；

f=(0:n-1) * (fs/n) ： 頻率范圍；

fs/2  ：Nyquist頻率(頻率范圍的中點)；

頻譜分析是一種將復噪聲號分解為較簡單信號的技術。真實世界中的信號可能由多種簡單信號疊加而成。找出一個信號在不同頻率下的信息（可能是幅度、功率、強度或相位等）的作法就是頻譜分析。

采樣定理：采樣頻率要大于信號頻率的兩倍。

N個采樣點經過FFT變換后得到N個點的以復數形式記錄的FFT結果。

假設采樣頻率為Fs，采樣點數為N。那么FFT運算的結果就是N個復數（或N個點），每一個復數就對應著一個頻率值以及該頻率信號的幅值和相位。

第一個點對應的頻率為0Hz（即直流分量），最后一個點N的下一個點對應采樣頻率Fs。其中任意一個采樣點n所代表的信號頻率：

這表明，頻譜分析得到的信號頻率最大為 (N-1)*Fs/N,對頻率的分辨能力是Fs/N。采樣頻率和采樣時間制約著通過FFT運算能分析得到的信號頻率上限，同時也限定了分析得到的信號頻率的分辨率。

每一個復數的模值對應該點所對應的頻率值的幅度特性，具體的定量關系如下：

假設信號由以下周期的原始信號疊加而成：

那么，在經過FFT分析后得到的第一個點的模值是A1的N倍，而且只有在FFT結果點對應的頻率在ω2,ω3時，其模值才明顯放大，在其他頻率點，模值接近于0。在這些模值明顯放大的點中，除第一個點之外的其它點模值是相應信號幅值的N/2倍。

每個復數的相位就是在該頻率值下信號的相位：φ2，φ3。

FFT結果有對稱性，通常我們只是用前半部分的結果，也就是小于采樣頻率一半的結果。同時也只有采樣頻率一半以內、具有一定幅值的信號頻率才是真正的信號頻率。

下面就用python案例進行說明

案例1

import numpy as npimport pylab as plimport math# 采樣頻率fs=1048# 采樣步長t = [x/1048.0 for x in range(1048)]"""設計的采樣值假設信號y由4個周期信號疊加所得,如下所示"""y = [ 3.0 * np.cos(2.0 * np.pi * 50 * t0 - np.pi * 30/180)     + 1.5 * np.cos(2.0 * np.pi * 75 * t0 + np.pi * 90/180)     +  1.0 * np.cos(2.0 * np.pi * 150 * t0 + np.pi * 120/180)     +  2.0 * np.cos(2.0 * np.pi * 220 * t0 + np.pi * 30/180)     for t0 in t ]pl.plot(t,y)pl.xlabel('time(s)')pl.title("original signal")pl.show()

"""現在對上述信號y在0-1秒時間內進行頻譜分析，本案例中采樣頻率為1048Hz,即單位時間內采樣點數為1048"""# 采樣點數N=len(t)# 采樣頻率fs=1048.0# 分辨率df = fs/(N-1)# 構建頻率數組f = [df*n for n in range(0,N)]Y = np.fft.fft(y)*2/N  #*2/N 反映了FFT變換的結果與實際信號幅值之間的關系absY = [np.abs(x) for x in Y]      #求傅里葉變換結果的模pl.plot(f,absY)pl.xlabel('freq(Hz)')pl.title("fft")pl.show()

案例2

from scipy.fftpack import fft, fftshift, ifftfrom scipy.fftpack import fftfreqimport numpy as npimport matplotlib.pyplot as plt%matplotlib inline"""t_s:采樣周期t_start:起始時間t_end:結束時間"""t_s = 0.01t_start = 0.5t_end = 5t = np.arange(t_start, t_end, t_s)f0 = 5f1 = 20# 繪制圖表plt.figure(figsize=(10, 12))# 構建原始信號序列y = 1.5*np.sin(2*np.pi*f0*t) + 3*np.sin(2*np.pi*20*t) + np.random.randn(t.size)ax=plt.subplot(511)ax.set_title('original signal')plt.tight_layout()plt.plot(y)"""FFT(Fast Fourier Transformation)快速傅里葉變換"""Y = fft(y)ax=plt.subplot(512)ax.set_title('fft transform')plt.plot(np.abs(Y))"""Y = fftshift(X) 通過將零頻分量移動到數組中心，重新排列傅里葉變換 X。"""shift_Y = fftshift(Y)ax=plt.subplot(513)ax.set_title('shift fft transform')plt.plot(np.abs(shift_Y))"""得到正頻率部分"""pos_Y_from_fft = Y[:Y.size//2]ax=plt.subplot(514)ax.set_title('fft transform')plt.tight_layout()plt.plot(np.abs(pos_Y_from_fft))"""直接截取 shift fft結果的前半部分"""pos_Y_from_shift = shift_Y[shift_Y.size//2:]ax=plt.subplot(515)ax.set_title('shift fft cut')plt.plot(np.abs(pos_Y_from_shift))plt.show()

眾所周知，傅里葉變換的快速算法FFT可以用來對信號的頻域特征進行分析，然而，FFT僅能用于平穩信號的分析，對于非平穩信號，則需要采用短時傅里葉變換（STFT）進行分析。

對于非平穩信號，短時傅里葉變換所采用的策略是在信號上面加窗，一般是hamming窗，當然也可以是其它類型的窗函數，加窗之后的信號被分割為一組短長度子序列，子序列可以近似的看為平穩序列，可以用傅里葉變換的方式去進行分析，這也是短時傅里葉變換的精髓所在。

用STFT進行腦電信號分析一般有兩種思路，第一種是用STFT來分離EEG信號的波段，從而求得每一個波段的能量作為特征（alpha、beta、theta、gamma、deta）。第二種是利用STFT計算功率譜密度作為特征，功率譜密度（PSD）特征可以針對整個信號子序列也可以針對子序列中特定的波段來計算。這兩種思路中，第二種思路用的比較廣，下面對其進行說明。

matlab中進行STFT的函數為spectrogram，計算功率譜密度（PSD）時使用如下格式：

[S,F,T,P]=spectrogram(x,window,noverlap,nfft,fs)

其中，S為輸入信號x的短時傅里葉變換，F為頻率向量，T為時間向量，P為功率譜密度矩陣，x為輸入信號，window為時間窗，noverlap為overlap的點數，如果為0就是沒有overlap，nfft為DFT的點數，fs為采樣頻率。其中，F向量的維度和P的行數一致，可以根據F向量來選取特定波段的PSD，還可以將alpha、beta、theta、gamma、deta這幾個波段分別分為幾個窄波段，提取窄帶PSD。

“FFT頻譜分析原理與python的實現”的內容就介紹到這里了，感謝大家的閱讀。如果想了解更多行業相關的知識可以關注億速云網站，小編將為大家輸出更多高質量的實用文章！