一、c語言編寫的FFT程序

c語言程序

正弦波表達式為： s(t) = 0.6 sin( 2π 50t ) 和s(t) = 0.6 sin( 2π 500t )
頻率為8000Hz,近似為8192
FFT程序來源于徐士良老師算法書中
1.FFT.c

#include "stdio.h" #include "math.h" #include "G:\C code\FFTDuiBi\FFTDuiBi\kfft.c"//注意路徑 #include <time.h>#define PI 3.1415926535main() {int i, j;double pr[8192], pi[8192], fr[8192], fi[8192], t[8192];clock_t begin, end;double cost1, cost2, persent;for (i = 0; i <8192; i++) //生成輸入信號{t[i] = i*0.001;pr[i] = 0.6*sin(2 * PI * 500 * i) + 0.6*sin(2 * PI * 50 * i); pi[i] = 0.0; //0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i)}begin = clock(); //begin=time(NULL); //開始記錄時間kfft(pr, pi, 8192, 13, fr, fi); //調用FFT函數end = clock(); //結束記錄時間cost1 = (double)(end - begin ) / CLOCKS_PER_SEC;//cost1 = (double)(end - begin);for (i = 0; i<8192; i++){printf("%d\t%lf\n", i, pr[i]); //輸出結果}printf("花費時間為：%lfs",cost1);}

2.kfft.c

#include "math.h" void kfft(pr, pi, n, k, fr, fi) int n, k; double pr[], pi[], fr[], fi[]; {int it, m, is, i, j, nv, l0;double p, q, s, vr, vi, poddr, poddi;for (it = 0; it <= n - 1; it++) //將pr的實部和虛部循環賦值給fr[]和fi[]{m = it;is = 0;for (i = 0; i <= k - 1; i++){j = m / 2;is = 2 * is + (m - 2 * j);m = j;}fr[it] = pr[is];fi[it] = pi[is];}pr[0] = 1.0;pi[0] = 0.0;p = 6.283185306 / (1.0*n);pr[1] = cos(p); //將w=e^-j2pi/n用歐拉公式表示pi[1] = -sin(p);for (i = 2; i <= n - 1; i++) //計算pr[]{p = pr[i - 1] * pr[1];q = pi[i - 1] * pi[1];s = (pr[i - 1] + pi[i - 1])*(pr[1] + pi[1]);pr[i] = p - q; pi[i] = s - p - q;}for (it = 0; it <= n - 2; it = it + 2){vr = fr[it];vi = fi[it];fr[it] = vr + fr[it + 1];fi[it] = vi + fi[it + 1];fr[it + 1] = vr - fr[it + 1];fi[it + 1] = vi - fi[it + 1];}m = n / 2;nv = 2;for (l0 = k - 2; l0 >= 0; l0--) //蝴蝶操作{m = m / 2;nv = 2 * nv;for (it = 0; it <= (m - 1)*nv; it = it + nv)for (j = 0; j <= (nv / 2) - 1; j++){p = pr[m*j] * fr[it + j + nv / 2];q = pi[m*j] * fi[it + j + nv / 2];s = pr[m*j] + pi[m*j];s = s*(fr[it + j + nv / 2] + fi[it + j + nv / 2]);poddr = p - q;poddi = s - p - q;fr[it + j + nv / 2] = fr[it + j] - poddr;fi[it + j + nv / 2] = fi[it + j] - poddi;fr[it + j] = fr[it + j] + poddr;fi[it + j] = fi[it + j] + poddi;}}for (i = 0; i<n; i++){pr[i] = sqrt(fr[i] * fr[i] + fi[i] * fi[i]); //幅度的計算}return; }

編譯結果并繪圖

在TCC下進行編譯，并用gnuplot繪圖

plot "<FFT.exe" w l

二、python編寫

python程序

import numpy as np#導入一個數據處理模塊 import timeimport matplotlib.pyplot as plt#導入一個繪圖模塊# 依據快速傅里葉算法得到信號的頻域 def test_fft():sampling_rate = 8192 # 采樣率fft_size = 8192 # FFT取樣長度t = np.arange(0, 8.192, 1.0 / sampling_rate)#np.arange(起點，終點，間隔)產生8.192s長的取樣時間x=0.6*np.sin(2*np.pi*500*t)+0.6*np.sin(2*np.pi*50*t)# 兩個正弦波疊加，500HZ和50HZ# N點FFT進行精確頻譜分析的要求是N個取樣點包含整數個取樣對象的波形。# 因此N點FFT能夠完美計算頻譜對取樣對象的要求是n*Fs/N（n*采樣頻率/FFT長度），# 因此對8KHZ和512點而言，完美采樣對象的周期最小要求是8000/512=15.625HZ,# 所以156.25的n為10,234.375的n為15。t1=time.perf_counter()xs = x[:fft_size]# 從波形數據中取樣fft_size個點進行運算xf = np.fft.rfft(xs) / fft_size # 返回fft_size/2+1 個頻率t2 = time.perf_counter()#利用np.fft.rfft()進行FFT計算，rfft()是為了更方便對實數信號進行變換，# 由公式可知 / fft_size為了正確顯示波形能量# rfft函數的返回值是N/2+1個復數，分別表示從0(Hz)到sampling_rate/2(Hz)的分。# 于是可以通過下面的np.linspace計算出返回值中每個下標對應的真正的頻率：freqs = np.linspace(0, sampling_rate*10, fft_size/2+1 ) # 表示頻率#freqs = np.linspace(0, sampling_rate/2 , fft_size/2 + 1) # 表示頻率xfp = 20 * np.log10(np.clip(np.abs(xf), 1e-20, 1e100))#xfp = np.abs(xf) # 代表信號的幅值，即振幅# 最后我們計算每個頻率分量的幅值，并通過 20*np.log10()將其轉換為以db單位的值。# 為了防止0幅值的成分造成log10無法計算，我們調用np.clip對xf的幅值進行上下限處理plt.figure(figsize=(8, 4))plt.subplot(211)plt.plot(t[:fft_size], xs)plt.xlabel(u"時間(秒)", fontproperties='FangSong')plt.title(u"500Hz和50Hz的波形和頻譜", fontproperties='FangSong')plt.subplot(212)plt.plot(freqs, xfp)plt.xlabel(u"頻率(Hz)", fontproperties='FangSong')#字體FangSongplt.ylabel(u'幅值', fontproperties='FangSong')plt.subplots_adjust(hspace=0.4)'''subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=None)有六個可選參數來控制子圖布局。值均為0~1之間。其中left、bottom、right、top圍成的區域就是子圖的區域。wspace、hspace分別表示子圖之間左右、上下的間距。實際的默認值由matplotlibrc文件控制的。'''plt.show()print(f"python的FFT花費的時間為：{t2-t1}")test_fft()

結果顯示：

花費時間對比

python運行FFT花費時間：

python運行FFT花費時間為：3.1e-6

c運行FFT花費時間：

c運行FFT花費時間為：0.002000s

總結

以上是生活随笔為你收集整理的频谱分析：c和python对比FFT的效率并画出幅度谱的全部內容，希望文章能夠幫你解決所遇到的問題。

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