用 python 實(shí)現(xiàn)FFT，繪制頻譜圖

關(guān)鍵詞 ：fft , scipy 庫(kù), fftshift ，單邊譜，雙邊譜，頻譜泄露

目錄

• 用 python 實(shí)現(xiàn)FFT，繪制頻譜圖
• • 前言
  • 代碼
  • 結(jié)果
  • 總結(jié)和討論

---

前言

之前都是在matlab上實(shí)現(xiàn)FFT，現(xiàn)在因?yàn)樾枰?#xff0c;在python上進(jìn)行實(shí)現(xiàn)，在此做一個(gè)記錄。

---

代碼

直接上代碼

import numpy as np from scipy.fftpack import fft,fftshift import matplotlib.pyplot as pltN = 1024 # 采樣點(diǎn)數(shù) sample_freq=120 # 采樣頻率 120 Hz, 大于兩倍的最高頻率 sample_interval=1/sample_freq # 采樣間隔 signal_len=N*sample_interval # 信號(hào)長(zhǎng)度 t=np.arange(0,signal_len,sample_interval)signal = 5 + 2 * np.sin(2 * np.pi * 20 * t) + 3 * np.sin(2 * np.pi * 30 * t) + 4 * np.sin(2 * np.pi * 40 * t) # 采集的信號(hào)fft_data = fft(signal)# 在python的計(jì)算方式中，fft結(jié)果的直接取模和真實(shí)信號(hào)的幅值不一樣。 # 對(duì)于非直流量的頻率，直接取模幅值會(huì)擴(kuò)大N/2倍， 所以需要除了N乘以2。 # 對(duì)于直流量的頻率(0Hz)，直接取模幅值會(huì)擴(kuò)大N倍，所以需要除了N。 fft_amp0 = np.array(np.abs(fft_data)/N*2) # 用于計(jì)算雙邊譜 fft_amp0[0]=0.5*fft_amp0[0] N_2 = int(N/2) fft_amp1 = fft_amp0[0:N_2] # 單邊譜 fft_amp0_shift = fftshift(fft_amp0) # 使用fftshift將信號(hào)的零頻移動(dòng)到中間# 計(jì)算頻譜的頻率軸 list0 = np.array(range(0, N)) list1 = np.array(range(0, int(N/2))) list0_shift = np.array(range(0, N)) freq0 = sample_freq*list0/N # 雙邊譜的頻率軸 freq1 = sample_freq*list1/N # 單邊譜的頻率軸 freq0_shift=sample_freq*list0_shift/N-sample_freq/2 # 零頻移動(dòng)后的頻率軸# 繪制結(jié)果 plt.figure() # 原信號(hào) plt.subplot(221) plt.plot(t, signal) plt.title(' Original signal') plt.xlabel('t (s)') plt.ylabel(' Amplitude ') # 雙邊譜 plt.subplot(222) plt.plot(freq0, fft_amp0) plt.title(' spectrum two-sided') plt.ylim(0, 6) plt.xlabel('frequency (Hz)') plt.ylabel(' Amplitude ') # 單邊譜 plt.subplot(223) plt.plot(freq1, fft_amp1) plt.title(' spectrum single-sided') plt.ylim(0, 6) plt.xlabel('frequency (Hz)') plt.ylabel(' Amplitude ') # 移動(dòng)零頻后的雙邊譜 plt.subplot(224) plt.plot(freq0_shift, fft_amp0_shift) plt.title(' spectrum two-sided shifted') plt.xlabel('frequency (Hz)') plt.ylabel(' Amplitude ') plt.ylim(0, 6)plt.show()

結(jié)果

分別繪制了原信號(hào)，雙邊譜，單邊譜和移動(dòng)零頻后的譜

總結(jié)和討論

概要：設(shè)定的原信號(hào)的頻率分量分別有 直流量(0Hz，幅值為5)，20Hz (幅值為2)，30Hz(幅值為3)，40Hz(幅值為4)，采樣點(diǎn)數(shù)N為1024，采樣頻率Fs為120Hz

計(jì)算可知, 頻率分辨率為$\frac{120}{1024}$, 信號(hào)長(zhǎng)度L=$\frac{N}{Fs}=\frac{1024}{120}$

關(guān)于頻譜泄露的考慮:
20Hz ，30Hz，40Hz 分量對(duì)應(yīng)的周期為$\frac{1}{20}$(s),$\frac{1}{30}$(s),$\frac{1}{40}$(s),信號(hào)長(zhǎng)度L($\frac{1024}{120}$)能被30Hz對(duì)應(yīng)的周期($\frac{1}{30}$)整除，因此，頻譜圖上30Hz對(duì)應(yīng)的幅值都是準(zhǔn)確的，但是信號(hào)長(zhǎng)度L不能被20Hz和40Hz對(duì)應(yīng)的周期整除，出現(xiàn)了頻譜泄露，因此幅值則略有偏差。從圖上也可以看到，20Hz和40Hz對(duì)應(yīng)的峰的底部有一定的寬度。

----------------------- The END ----------------------------------------

總結(jié)

以上是生活随笔為你收集整理的用 python 实现FFT，绘制频谱图的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

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