一、c語(yǔ)言編寫(xiě)的FFT程序

c語(yǔ)言程序

正弦波表達(dá)式為： s(t) = 0.6 sin( 2π 50t ) 和s(t) = 0.6 sin( 2π 500t )
頻率為8000Hz,近似為8192
FFT程序來(lái)源于徐士良老師算法書(shū)中
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++) //生成輸入信號(hào){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); //開(kāi)始記錄時(shí)間kfft(pr, pi, 8192, 13, fr, fi); //調(diào)用FFT函數(shù)end = clock(); //結(jié)束記錄時(shí)間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]); //輸出結(jié)果}printf("花費(fèi)時(shí)間為：%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的實(shí)部和虛部循環(huán)賦值給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++) //計(jì)算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]); //幅度的計(jì)算}return; }

編譯結(jié)果并繪圖

在TCC下進(jìn)行編譯，并用gnuplot繪圖

plot "<FFT.exe" w l

二、python編寫(xiě)

python程序

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

結(jié)果顯示：

花費(fèi)時(shí)間對(duì)比

python運(yùn)行FFT花費(fèi)時(shí)間：

python運(yùn)行FFT花費(fèi)時(shí)間為：3.1e-6

c運(yùn)行FFT花費(fèi)時(shí)間：

c運(yùn)行FFT花費(fèi)時(shí)間為：0.002000s

總結(jié)

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

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