0. 預備知識

快速傅里葉變換旨在解決離散傅里葉變換DFT計算量大效率低的問題。當我們想要抑制噪聲提取出某段信號中的有效信息時，如系統模型辨識或者是使用高精度力傳感器測量人體腕部寸關尺脈搏信號這類應用，應該如何設計采樣流程？

• 首先，應當考慮采樣頻率$f_s$的問題，根據香農采樣定理，采樣頻率應大于等于目標信號頻率$f$最高頻段的2倍，工程中通常取2.56到4倍的頻率。采樣頻率可以直接配置傳感器的采樣觸發信號，對于采樣頻率固定的設備，如普通家用攝像頭，則需要根據應用選擇設備型號。采樣頻率最好是2的冪次。
• 其次，采樣時間的問題，在確定采樣頻率后等同于確定采樣點數量$N$。采樣點數量越多，則FFT變換后生成的頻譜的頻段間隔越小，即分辨力越高。采樣數據點數量最好是采樣頻率的倍數關系。

FFT應用時需要主意的點包括：

• FFT輸出的頻譜是雙邊對稱的，關于第$N/2$個數據點左右對稱，每一個數據點是一個復數 $a+b\mathbf{j}$，這些數據點只有幅值和相位是具有實際意義的。
• FFT輸入$N$個數據點，輸出$N$個頻段的復數，可以通過這些復數計算出對應的相位和未歸一化的幅值。對幅值部分進行歸一化時需要主意的是，直流分量需要除以$N$，剩余分量需要除以$N/2$。
• 由于FFT是雙邊對稱的，因此頻率和幅值部分應該取前半段，頻段點為：$f=(n?1)f_s/N$， 需要注意，因為FFT輸出的結果是雙邊的，只有前半段的頻譜是有意義的，因此這里$n$的取值范圍是$n=1,...,N/2$

當我們理解fft的過程，就可以調整不同頻段的幅值，重新組合出一個過濾后的信號。

1. 代碼示例

1.1 代碼

借用這位老哥的案例： https://blog.csdn.net/qq_27825451/article/details/88553441?spm=1001.2014.3001.5506 。
我們重新寫一份FFT的代碼，分析信號的頻譜：

''' Author: Dianye Huang Date: 2023-01-14 10:26:47 LastEditors: Dianye Huang LastEditTime: 2023-01-14 14:37:02 Description: '''import numpy as np from scipy.fftpack import fft, ifft import matplotlib.pyplot as pltclass myFFT:def __init__(self):passdef get_spetrum(self, data, fs, flag_plt=False):N = len(data)fft_y = fft(data) # get complex numberabs_y = np.abs(fft_y) # get magnitudeang_y = np.angle(fft_y) # get phasenrm_y = abs_y/(N/2) # get normailzed magnitude (A0/N, A1.../(N/2))nrm_y[0] = nrm_y[0]/2 nmh_y = nrm_y[:int(N/2)] # half normalized magnitudeagh_y = ang_y[:int(N/2)] # half normalized anglespes = np.arange(int(N/2))*fs/N # specturm axis if flag_plt:plt.figure()x = np.arange(N)plt.subplot(2,3,1)plt.plot(x/fs, data)plt.title('raw signal')plt.xlabel('Time (s)')plt.subplot(2,3,4)plt.plot(x[:50]/fs, data[:50])plt.title('partial raw signal')plt.xlabel('Time (s)')plt.subplot(2,3,2)plt.plot(x, abs_y)plt.title('magnitudes')plt.xlabel('Sample index')plt.subplot(2,3,5)plt.plot(x, ang_y)plt.title('angles')plt.xlabel('Sample index')plt.subplot(2,3,3)plt.plot(x, nrm_y)plt.title('nomalized magnitude')plt.xlabel('Sample index')plt.subplot(2,3,6)plt.plot(spes, nmh_y)plt.title('half nomalized magnitude')plt.xlabel('Frequency (Hz)')plt.show()return nmh_y, agh_y, spesdef get_fft(self, data, fs):N = len(data)fft_y = fft(data)tmp_arr = np.arange(0, fs/2, fs/N)freq_y = np.hstack((tmp_arr, np.flip(tmp_arr, axis=0)))return fft_y, freq_ydef create_demo_signal(self, N=2800, fs=1400):# create a signal along with sample ratet = np.arange(0, N/fs, 1/fs)y = 3 + 7*np.sin(2*np.pi*200*t) + 5*np.sin(2*np.pi*400*t) + 3*np.sin(2*np.pi*600*t)noise_arr = np.random.normal(0, 2, N)return t, y+noise_arr, fsif __name__ == '__main__':mfft=myFFT()t, data, fs = mfft.create_demo_signal(N=4096, fs=2048)# mags, angs, spes = mfft.get_spetrum(data, fs, flag_plt=True) fft_y, freq_y = mfft.get_fft(data, fs)fft_y[freq_y>450] = 0sig = ifft(fft_y).realplt.figure()plt.plot(t[:50], data[:50])plt.plot(t[:50], sig[:50])plt.show()

1.2 總結

• $fft$和$ifft$的求解包在scipy.fftpack中
• 復數的幅值和相位可以使用numpy包中的np.abs()和np.angle()函數

示例程序運行結果如下：

1.2.1 分析頻譜，右下角為最終結果

if __name__ == '__main__':mfft=myFFT()t, data, fs = mfft.create_demo_signal(N=4096, fs=2048)mags, angs, spes = mfft.get_spetrum(data, fs, flag_plt=True)

1.2.2 簡單的濾波，除去大于450Hz的分量利用ifft重新組合信號的結果

if __name__ == '__main__':mfft=myFFT()t, data, fs = mfft.create_demo_signal(N=4096, fs=2048)fft_y, freq_y = mfft.get_fft(data, fs)fft_y[freq_y>450] = 0sig = ifft(fft_y)plt.figure()plt.plot(t[:50], data[:50])plt.plot(t[:50], sig[:50])plt.show()

以上，
Dianye
2023.01.14

總結

以上是生活随笔為你收集整理的快速傅里叶变换FFT和逆变换的python编程的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。