0. 預(yù)備知識

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

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

FFT應(yīng)用時需要主意的點包括：

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

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

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 總結(jié)

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

示例程序運行結(jié)果如下：

1.2.1 分析頻譜，右下角為最終結(jié)果

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重新組合信號的結(jié)果

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。