快速傅立葉變換(FFT)的計(jì)算機(jī)實(shí)現(xiàn)

1、信號與系統(tǒng)課程設(shè)計(jì)FFT的計(jì)算機(jī)實(shí)現(xiàn)快速傅里葉變換（FFT）的計(jì)算機(jī)實(shí)現(xiàn)賴智鵬華中科技大學(xué)電氣與電子工程學(xué)院0809班UEmail: 摘要：本文是信號與系統(tǒng)課程的課程設(shè)計(jì)，旨在熟悉FFT的計(jì)算過程，結(jié)合DFT物理意義和實(shí)驗(yàn)結(jié)果加深對傅立葉變換的理解。文章首先用MATLAB對一個(gè)簡單信號進(jìn)行FFT仿真，得出頻譜圖；其次完成了FFT的C語言實(shí)現(xiàn)，結(jié)合MATLAB作圖及數(shù)據(jù)處理功能得出了C實(shí)現(xiàn)下的FFT結(jié)果；最后，討論分析實(shí)驗(yàn)結(jié)果。關(guān)鍵詞：DFT、基-2按時(shí)間抽取FFT算法、MATLAB、C、頻譜、物理意義1. 算法描述1) DFT的運(yùn)算量2) 減少運(yùn)算的方法：v 化長序列為短序列。如將長度為N的序

2、列分解為兩個(gè)長度為N/2的序列v 利用的性質(zhì)（注：本文中的C程序未用到此性質(zhì)） 3) C程序采用基-2按時(shí)間抽取的FFT算法設(shè)輸入序列長度為 (M為正整數(shù))，將該序列按時(shí)間順序的奇偶分解為越來越短的子序列，稱為基2按時(shí)間抽取的FFT算法,也稱為Coolkey-Tukey算法。若N不滿足條件，則人為地加上若干零值，使。2. 實(shí)驗(yàn)設(shè)計(jì)與步驟為簡單起見,同時(shí)不失一般性,本實(shí)驗(yàn)采用三個(gè)余弦成分和一個(gè)支流偏置成分疊加所得的信號作為信號源。3. MATLAB的FFT仿真和C的FFT實(shí)現(xiàn)1) 關(guān)于頻譜的理論分析：2) MATLAB的FFT仿真v MATLAB程序代碼function output_args

3、= FFT2( x0,x1,f1,w1,x2,f2,w2,x3,f3,w3,fs )%這是一個(gè)自定義函數(shù)，輸入被采樣的信號的參數(shù)和采樣頻率，然后輸出原信號波形、采樣信號序列、采樣序列幅度譜和相位譜。函數(shù)規(guī)定信號由一個(gè)直流成分（大小為x0），三個(gè)余弦成分（各自頻率分別為f1,f2,f3,幅值為x1,x2,x3,初相位為w1,w2,w3）,采樣頻率為fs，采樣時(shí)間為1秒。x0,x1,f1,w1,x2,f2,w2,x3,f3,w3,fs作為參數(shù)輸入t=0:1/fs:1-1/fs;%定義采樣時(shí)刻N(yùn)=length(t);%采樣序列長度s=x0+x1*cos(2*pi*f1*t+pi/180*w1)+x2

4、*cos(2*pi*f2*t+pi/180*w2)+x3*cos(2*pi*f3*t+pi/180*w3);%采樣信號y=fft(s);%快速傅立葉變換tt=0:1/10000:1;%原信號描點(diǎn)ss=x0+x1*cos(2*pi*f1*tt+pi/180*w1)+x2*cos(2*pi*f2*tt+pi/180*w2)+x3*cos(2*pi*f3*tt+pi/180*w3);%原信號figure;plot(tt,ss);grid;title(原始信號);%原信號波形figure;subplot(3,1,1);plot(0:N-1,s(1:N),-o);xlim(0 N-1);grid;tit

5、le(采樣信號);subplot(3,1,2);plot(0:N-1,abs(y(1:N),-o);xlim(0 N-1);grid;xlabel(k);ylabel(幅度);title(理想采樣信號的幅度譜);subplot(3,1,3);plot(0:N-1,angle(y(1:N),-o);grid;xlabel(k);ylabel(相位);axis(0 N-1 -pi pi);title(理想采樣信號的相位譜);end;v 改變采樣頻率 依次鍵入：%其中FFT2()是自定義函數(shù),對信號在1秒內(nèi)以頻率fs進(jìn)行采樣FFT2(1,1,1,90,2,2,180,3,3,180,32)FFT2(

6、1,1,1,90,2,2,180,3,3,180,16)FFT2(1,1,1,90,2,2,180,3,3,180,8)FFT2(1,1,1,90,2,2,180,3,3,180,4)原始信號如圖：圖1 圖2、fs =32Hz圖3、fs =16Hz圖4、fs =8Hz圖5、fs =4Hz注意到：1、幅度譜，除卻k=0外，圖形呈軸對稱分布；2、相位譜，N為偶數(shù)時(shí)除卻k=0和N/2外（N為奇數(shù)時(shí)，除卻k=0外），圖形呈中心對稱分布。理論上,由于復(fù)指數(shù)的周期性，長度為N（假設(shè)N為偶數(shù)，N為奇數(shù)時(shí)類似）的的DFT頻譜分析，可把k=N/2,N/2+1,N-1看作是負(fù)頻率-N/2,-N/2+1,-N/2+

7、2,-2,-1,特別的，對于實(shí)序列而言，正負(fù)頻率成分的幅值必須相等，而初相位必須相反。下面取fs=8Hz和4Hz部分計(jì)算結(jié)果分析，其他情況可類似處理?？煽吹剑篺s=8Hz時(shí)，由頻譜結(jié)合物理意義計(jì)算的得到的余弦成分的幅值和初相位與原信號一致；fs=4Hz時(shí)，由頻譜結(jié)合物理意義計(jì)算的得到的余弦成分的幅值和初相位與原信號不同。對此，可根據(jù)抽樣定理解釋。因?yàn)閒1=1Hz，f2=2Hz， f3=3Hz，所以連續(xù)時(shí)間信號的截至頻率fm=f3=3Hz，所以抽樣頻率fs=4Hz時(shí)，fs %改變信號的直流成分（由1變?yōu)?），原信號和頻譜如下FFT2(2,1,1,90,2,2,180,3,3,180,32)圖9

8、圖10 fs=32Hz變化趨勢：幅度譜中k=0對應(yīng)的幅度變?yōu)?N=64(原為32)，其余處無明顯變化，相位譜中，k=0,1,2,3對應(yīng)的相位值都不變。 鍵入： %改變第三個(gè)余弦成分的幅值（由3改為4），原信號圖和頻譜圖如下FFT2(2,1,1,90,2,2,180,4,3,180,32)圖11 圖12 fs=32Hz變化趨勢：幅度譜中k=3對應(yīng)的幅度變?yōu)?4（4*32/2=64,原來這個(gè)值為3*32/2=48），其余處無明顯變化，相位譜中k=0,1,2,3對應(yīng)的相位值都未發(fā)生變化。 鍵入： %改變第三個(gè)余弦成分的頻率（由3Hz改為4Hz），原信號圖和頻譜如下FFT2(2,1,1,90,2,2,

9、180,4,4,180,32)圖13 圖14 fs=32Hz變化趨勢：幅度譜中，k=4處幅度為N/2*4=64,k=3處幅度變?yōu)?，k=0,1,2處幅度不變；相位譜中，k=4處相位變?yōu)?80，k=0,1,2,3處相位值不變。 鍵入：%改變第三個(gè)余弦成分的初相位（由180改為270），原信號和頻譜如下圖FFT2(2,1,1,90,2,2,180,4,4,270,32)圖15 圖16 fs=32Hz變化趨勢：幅度譜波形無明顯變化；相位譜中，k=4處相位值變?yōu)?90度，即270度。由此可見，改變信號的波幅、頻率和相位，幅度譜和相位譜將發(fā)生相應(yīng)的變化，即恰合頻譜的物理意義。3) C語言實(shí)現(xiàn)FFTv C

10、代碼/*此代碼用于fft計(jì)算，采用基-2按時(shí)間抽取的FFT算法Decimation-in-Time(DIT)(Coolkey-Tukey)。為方便起見，同時(shí)不失一般性，把含有一個(gè)直流成分和三個(gè)正弦成分的信號作為被采樣信號，信號的輸入由函數(shù)input()完成，若想改變信號波形只需改變input()函數(shù)代碼。代碼分別定義了復(fù)數(shù)結(jié)構(gòu)體、復(fù)數(shù)運(yùn)算和碼位倒讀函數(shù)。計(jì)算結(jié)果分別在命令窗口和文件e:keshe.txt中輸出，keshe.txt文件數(shù)據(jù)可用于matlab作圖分析*/#include#include/*圓周率*/double pi=3.9793;int M,N;/*定義復(fù)數(shù)*/struct Co

11、mplex_double real;double img;*a,*b;/*定義2的冪計(jì)算*/int x_2(int a)int i,r;r=1;for(i=0;i=0;j-)ii+=(int)(i/x_2(j)*x_2(M-1-j);/*先把i用二進(jìn)制表示，然后碼位倒讀*/i-=(int)(i/x_2(j)*x_2(j);return(ii);/*定義復(fù)數(shù)乘法*/struct Complex_ Multi(struct Complex_ a,struct Complex_ b)struct Complex_ c;c.real=a.real*b.real-a.img*b.img;c.img=a.

12、real*b.img+a.img*b.real;return(c);/*定義復(fù)數(shù)冪運(yùn)算*/struct Complex_ W_K(struct Complex_ W,int k)int i;struct Complex_ x,y;x=W;if(k=0)x.real=1;x.img=0;return(x);elsefor(i=1;ik;i+)y=Multi(x,W);x=y;return(x);/*定義復(fù)數(shù)加法*/struct Complex_ Add(struct Complex_ a,struct Complex_ b)struct Complex_ c;c.real=a.real+b.re

13、al;c.img=a.img+b.img;return(c);/*離散序列的輸入*/void input()int NN,i;float f1,f2,f3,a0,a1,a2,a3,w1,w2,w3;M=0;/*在一秒鐘內(nèi)，采點(diǎn)數(shù)*/printf(采點(diǎn)數(shù):);scanf(%d,&NN);printf(n);/*信號的直流成分*/printf(信號的直流成分:n);scanf(%f,&a0);/*信號有三個(gè)余弦成分*/*余弦成分1*/printf(第1個(gè)正弦成分的頻率:n);scanf(%f,&f1);printf(幅值:n);scanf(%f,&a1);printf(初相位：n);scanf(%

14、f,&w1);/*余弦成分2*/printf(第2個(gè)正弦成分的頻率:n);scanf(%f,&f2);printf(幅值:n);scanf(%f,&a2);printf(初相位：n);scanf(%f,&w2);/*余弦成分3*/printf(第3個(gè)正弦成分的頻率:n);scanf(%f,&f3);printf(幅值:n);scanf(%f,&a3);printf(初相位：n);scanf(%f,&w3);for(i=0;x_2(i)NN;i+)M=i+1;N=x_2(M);/*N為不小于NN的最小的2的冪*/a=(struct Complex_*)calloc(N,sizeof(struct

15、 Complex_);/*動(dòng)態(tài)開辟N個(gè)單位*/b=(struct Complex_*)calloc(N,sizeof(struct Complex_);/*動(dòng)態(tài)開辟N個(gè)單位*/*把采樣信號用碼位倒讀的方法存入計(jì)算機(jī)中*/for(i=0;iNN;i+)areverse(i).real=a0+a1*cos(2*pi*f1*i/N+w1/180*pi)+a2*cos(2*pi*f2*i/N+w2/180*pi)+a3*cos(2*pi*f3*i/N+w3/180*pi);/*其中reverse()為碼位倒讀函數(shù)*/areverse(i).img=0;/*動(dòng)態(tài)空間空閑處，補(bǔ)零*/for(i=NN;iN

16、;i+)areverse(i).real=0;areverse(i).img=0;/*主函數(shù)*/void main()int i,j,k,m,n,q,r;FILE *fp;/*文件指針，指向存儲fft結(jié)果的文件*/struct Complex_ d,W,Wk,x,y;j=0;input();m=1;n=N;for(i=0;iM;i+)/*M級運(yùn)算*/n=n/2;m=m*2;q=x_2(i+1);/*q=2(i+1)*/W.real=cos(-2*pi/q);/*旋轉(zhuǎn)因子實(shí)部*/W.img=sin(-2*pi/q);/*旋轉(zhuǎn)因子虛部*/ for(j=0;jn;j+)/*n群運(yùn)算*/ for(k=

17、0;km/2;k+)Wk=W_K(W,k); bj*m+k=Add(aj*m+k,Multi(Wk,aj*m+k+m/2);/*Add()為復(fù)數(shù)相加函數(shù)，Multi()為復(fù)數(shù)相乘函數(shù)*/ for(k=m/2;km;k+) Wk=W_K(W,k); bj*m+k=Add(aj*m+k-m/2,Multi(Wk,aj*m+k); /*Add()為復(fù)數(shù)相加函數(shù)，Multi()為復(fù)數(shù)相乘函數(shù)*/ for(r=0;rN;r+) ar=br; /*將結(jié)果寫入文件e:keshe.txt*/if(fp=fopen(e:keshe.txt,w+)=NULL)printf(Cannot open file str

18、ike any key exit!);getch();exit(1);for(n=0;nN;n+)fprintf(fp,%16.14f %16.14fn,an.real,an.img);for(n=0;nN;n+)printf(%16.14f+%16.14fin,an.real,an.img);free(a);/*釋放內(nèi)存空間*/free(b);/*釋放內(nèi)存空間*/v 改變采樣頻率實(shí)驗(yàn)數(shù)據(jù)將由keshe.txt調(diào)入到MATLAB中，做出頻譜圖。 參數(shù)設(shè)置1圖17 fs=32Hz(C實(shí)現(xiàn))圖18 fs=32Hz（MATLAB仿真）注意到兩圖的幅度譜圖，很相似，而相位譜卻相差很大。用C實(shí)現(xiàn)的FFT

19、得到的相位譜也似乎沒有什么對稱性，但可以看到兩個(gè)相位譜在k=0,1,2,3和31,30,29處的相位值是基本一致的，而其他點(diǎn)處，相位值差異則很大。（表1，列出了具體的數(shù)值計(jì)算結(jié)果）表1kFFT（C）FFT（MATLAB）FFT結(jié)果相位值FFT結(jié)果相位值032.00 + 0.000i032.00 + 0.000i010.000 + 16.00i1.2.792e-15 - 16.00i1.2-32.00 + 1.000e-14i3.-32.00 - 9.138e-15i3.3-48.00 + 3.000e-14i3.-48.00 - 4.557e-14i3.4-1.000e-14 + 0.000i

20、3.-1.118e-14 - 9.265e-16i3.0605650.000 + 0.000i0-3.056e-15 - 1.025e-15i2.61.000e-14 + 1.000e-14i0.7.106e-15 - 4.069e-15i0.71.000e-14 + 1.000e-14i0.1.516e-15 - 5.824e-15i1.8-1.000e-14 + 1.000e-14i2.-1.015e-14 - 7.100e-15i2.5535990.000 - 1.000e-14i-1.57081.516e-15 + 7.352e-15i-1.41815100.000 + 0.000i0

21、-2.031e-15 - 2.031e-15i2.110.000 + 0.000i01.516e-14 - 5.075e-15i0.4475212-1.000e-14 + 0.000i3.-9.124e-15 - 9.265e-16i3.130.000 + 0.000i03.050e-15 - 1.019e-15i0.14-1.000e-14 + 0.000i3.3.050e-15 + 0.000i015-2.000e-14 + 1.000e-14i2.-1.397e-14 - 1.016e-14i2.161.000e-14 + 0.000i0-3.050e-15 + 0.000i3.17-2

22、.000e-14 - 1.000e-14i-2.67795-1.397e-14 + 1.016e-14i-2.512818-1.000e-14 + 0.000i3.3.050e-15 + 0.000i019-1.000e-14 + 0.000i3.3.050e-15 + 1.019e-15i-0.3587720-1.000e-14 + 0.000i3.-9.124e-15 + 9.265e-16i-3.04494211.000e-14 + 0.000i01.516e-14 + 5.075e-15i-0.44752220.000 + 0.000i0-2.031e-15 + 2.031e-15i-

23、2.35619230.000 + 1.000e-14i1.1.516e-15 - 7.352e-15i1.24-1.000e-14 - 1.000e-14i-2.35619-1.015e-14 + 7.100e-15i-2.55359250.000 - 1.000e-14i-1.57081.516e-15 + 5.824e-15i-1.36116261.000e-14 - 1.000e-14i-0.78547.106e-15 + 4.069e-15i-0.574327-1.000e-14 + 0.000i3.-3.056e-15 + 1.025e-15i-2.7233728-1.000e-14

24、 + 0.000i3.-1.118e-14 + 9.265e-16i-3.0605629-48.00 - 5.000e-14i-3.14159-48.00 + 4.557e-14i-3.1415930-32.00 - 2.000e-14i-3.14159-32.00 + 9.138e-15i-3.14159311.000e-14 - 16.00i-1.57082.792e-15 + 16.00i-1.5708結(jié)合圖17、18和表1，可以發(fā)現(xiàn)兩個(gè)相位譜在k=0,1,2,3和31,30,29處一致，而在其他點(diǎn)相差很大，注意到這些點(diǎn)的幅值很?。ɡ碚撋蠎?yīng)該為0，實(shí)際上卻由于計(jì)算機(jī)的誤差，而使之非0），

25、可認(rèn)為為0，可忽略這些點(diǎn)，而主要考慮主頻率成分。事實(shí)上，結(jié)合圖表可發(fā)現(xiàn)，C程序所得的結(jié)果也并非完全無對稱性的，并且可注意到所有不對稱點(diǎn)的相位值接近圓周率或0。這種差異是由于計(jì)算機(jī)的誤差導(dǎo)致的，在一般情況下，這不會造成太大的影響，而在相位處于邊界值（和0）時(shí)，這種誤差將計(jì)算結(jié)果發(fā)生很大變化。注：在先前編寫的C程序中，把復(fù)數(shù)實(shí)部和虛部的數(shù)據(jù)類型設(shè)置為float時(shí)，相位譜的不對稱更加嚴(yán)重，甚至C和MATLAB仿真的兩個(gè)相位譜在k= 31,30,29都不一致，此處實(shí)驗(yàn)結(jié)果是把數(shù)據(jù)類型改為double類型后（即提高精度后）所得的，仍存在這種不對稱，理論上可通過不斷提高精度進(jìn)行改進(jìn)，但只要是某頻率成分初相

26、位為或0，就存在這種嚴(yán)重不對稱的風(fēng)險(xiǎn)。 參數(shù)設(shè)置2圖19 fs=16Hz(C實(shí)現(xiàn)) 圖20 fs=16Hz（MATLAB仿真） 可發(fā)現(xiàn)C實(shí)現(xiàn)下的FFT幅度譜圖隨采樣頻率變化的趨勢與MATLAB仿真實(shí)驗(yàn)類似，而由于參數(shù)1實(shí)驗(yàn)相同的原因，C實(shí)現(xiàn)的FFT相位譜并不對稱。 參數(shù)設(shè)置3圖21 fs=8Hz(C實(shí)現(xiàn)) 圖22 fs=8Hz（MATLAB仿真）類似的，可發(fā)現(xiàn)幅度譜圖隨采樣頻率變化的趨勢與MATLAB仿真實(shí)驗(yàn)類似。而C和MATLAB實(shí)現(xiàn)下，k=6時(shí)的相位值一個(gè)為3.14，一個(gè)為-3.14，具體原因同上。而和相差，理論上初相位為是無差異的，但是體現(xiàn)在相位譜上就是造成不對稱。 參數(shù)設(shè)置4圖23 f

27、s=4Hz(C實(shí)現(xiàn))圖24 fs=4Hz（MATLAB仿真）同MATLAB仿真結(jié)果一樣出現(xiàn)混疊現(xiàn)象。 參數(shù)設(shè)置5圖25 fs=14Hz(C實(shí)現(xiàn)) 圖26 fs=14Hz（MATLAB仿真）注意到幅度譜和相位譜都有很大區(qū)別，其中相位譜的區(qū)別是由兩種不同的計(jì)算機(jī)語言規(guī)則造成的（正如之前提到的），而幅度譜的區(qū)別是由算法不同造成的， C程序是采用補(bǔ)零的方法把序列長度湊成，然后把長度為N的序列輸出為結(jié)果，這樣做的后果是改變了進(jìn)行DFT的序列， 得到的自然不是所求的結(jié)果，所以這個(gè)C程序，只能針對序列長度滿足的情況。v 改變被采樣信號波形 參數(shù)設(shè)置6圖27 fs=32Hz(C實(shí)現(xiàn))圖28 fs=32Hz(M

28、ATLAB仿真)二者幅度譜相同，相位譜只在中間部分有差異（此差異對實(shí)際分析影響很小，且也無實(shí)際意義，可忽略），而在k=0,1,2,3和31,30,29處兩幅圖的情況一致，故驗(yàn)證了之前關(guān)于出現(xiàn)圖17和圖18二者相位譜如此大差異原因的猜想。 參數(shù)設(shè)置7圖29 fs=32Hz(C實(shí)現(xiàn))圖30 fs=32Hz(MATLAB仿真)二者幅度譜圖相同，且較上一參數(shù)時(shí)發(fā)生的變化趨勢與MATLAB仿真試驗(yàn)中一樣；相位譜圖在k=0,1,2,3和31,30,29處相位值相同。 參數(shù)設(shè)置8圖31 fs=32Hz(C實(shí)現(xiàn))圖32 fs=32Hz(MATLAB仿真)二者幅度譜圖相同，且較上一參數(shù)時(shí)發(fā)生的變化趨勢與MATL

29、AB仿真試驗(yàn)中一樣；相位譜圖在k=0,1,2,4和31,30,28處相位值相同。至此，考察了序列長度和信號諧波成分對MATLAB的FFT仿真和C實(shí)現(xiàn)的FFT的影響，并結(jié)合了DFT的物理意義，對實(shí)驗(yàn)現(xiàn)象做出了解釋。對比MATLAB和C實(shí)現(xiàn)下的FFT結(jié)果的差異，做出了自己的解釋。4. 課程設(shè)計(jì)心得與自我評價(jià)這個(gè)課程設(shè)計(jì)的主要部分是在暑假完成的，課程設(shè)計(jì)是關(guān)于FFT的計(jì)算機(jī)實(shí)現(xiàn)，要求掌握FFT算法，并分別用MATLAB進(jìn)行FFT仿真和編寫實(shí)現(xiàn)FFT的C代碼。記得之前學(xué)習(xí)電路理論、復(fù)變函數(shù)、數(shù)理方程和信號系統(tǒng)時(shí)都遇到過把時(shí)域上的問題轉(zhuǎn)化為頻域上的問題解決的情況，電路中主要是Laplace變換的實(shí)際應(yīng)用

30、，而復(fù)變函數(shù)里面主要講的是各種變換的技術(shù)新性問題，數(shù)理方程里面研究的是一些數(shù)學(xué)方程的物理模型，而信號與系統(tǒng)課程則引入了離散的概念。首先是時(shí)域上的離散，由對連續(xù)時(shí)間進(jìn)行信號采樣，即x(t)xn，把連續(xù)時(shí)間信號的各種性質(zhì)在離散時(shí)間信號做推廣，如CTFTDTFT，Laplace變換Z-變換；其次是頻域上的離散，對連續(xù)頻率進(jìn)行信號采樣，即DTFTDFT，而FFT則是處理DFT的一種快速算法。課程設(shè)計(jì)完成的過程中遇到許多問題：首先是MATLAB的仿真。這時(shí)候問題不大，直接利用現(xiàn)成的內(nèi)置函數(shù)。其次是FFT的C實(shí)現(xiàn)。此處在代碼的編寫和程序的調(diào)試耗費(fèi)了較多的時(shí)間，主要問題在于對算法的理解不夠以及對C語言的一些問題未完全搞懂。后來在C實(shí)現(xiàn)的FFT結(jié)果與MATLAB的仿真結(jié)果對比時(shí)，遇到較多問題：1、初相位為180度（或-180度）或0的頻率成分在C實(shí)現(xiàn)的頻