FFT是DFT的高效算法,能夠?qū)r域信號轉(zhuǎn)化到頻域上,下面記錄下一段用python實現(xiàn)的FFT代碼。
# encoding=utf-8import numpy as npimport pylab as pl # 導(dǎo)入和matplotlib同時安裝的作圖庫pylabsampling_rate = 8000 # 采樣頻率8000Hzfft_size = 512 # 采樣點(diǎn)512,就是說以8000Hz的速度采512個點(diǎn),我們獲得的數(shù)據(jù)只有這512個點(diǎn)的對應(yīng)時刻和此時的信號值。t = np.linspace(0, 1, sampling_rate) # 截取一段時間,截取是任意的,這里取了0~1秒的一段時間。x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t) # 輸入信號序列,人工生成了一段信號序列,范圍在0~1秒xs = x[:fft_size] # 由上所述,我們只采樣了512個點(diǎn),所以我們只獲得了前512個點(diǎn)的數(shù)據(jù)xf = np.fft.rfft(xs)/fft_size # 調(diào)用np.fft的函數(shù)rfft(用于實值信號fft),產(chǎn)生長度為fft_size/2+1的一個復(fù)數(shù)向量,分別表示從0Hz~4000Hz的部分,這里之所以是4000Hz是因為Nyquist定理,采樣頻率8000Hz,則能恢復(fù)帶寬為4000Hz的信號。最后/fft_size是為了正確顯示波形能量freqs = np.linspace(0, sampling_rate//2, fft_size//2 + 1) # 由上可知,我們得到了數(shù)據(jù),現(xiàn)在產(chǎn)生0~4000Hz的頻率向量,方便作圖xfp = 20*np.log10(np.clip(np.abs(xf), 1e-20, 1e1000)) # 防止幅值為0,先利用clip剪裁幅度,再化成分貝pl.figure(figsize=(8, 4)) # 生成畫布pl.subplot(211) # 生成子圖,211的意思是將畫布分成兩行一列,自己居上面。pl.plot(t[:fft_size], xs) # 對真實波形繪圖pl.xlabel(u"time(s)")pl.title(u"The Wave and Spectrum of 156.25Hz and 234.375Hz")pl.subplot(212) # 同理pl.plot(freqs, xfp) # 對頻率和幅值作圖,xlabel是頻率Hz,ylabel是dBpl.xlabel(u"Hz")pl.subplots_adjust(hspace=0.4) # 調(diào)節(jié)繪圖參數(shù)pl.show()
代碼進(jìn)行了詳細(xì)標(biāo)注。有一個小細(xì)節(jié)是FFT對于取樣時間有要求。N點(diǎn)FFT進(jìn)行精確頻譜分析的要求是N個取樣點(diǎn)包含整數(shù)個取樣對象的波形。因此N點(diǎn)FFT能夠完美計算頻譜,對取樣對象的要求是n*Fs/N(n*采樣頻率/FFT長度)在本例中Fs = 8000Hz,N=512 base_freq=15.625Hz 所以本例中給出了頻率為156.25Hz(n=10)和234.375Hz(n=15)做例子。
效果如下:
以上就是本文的全部內(nèi)容,希望對大家的學(xué)習(xí)有所幫助,也希望大家多多支持武林網(wǎng)之家。
新聞熱點(diǎn)
疑難解答
圖片精選