快速傅立葉變換(FFT)的計算機實現(xiàn)

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

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

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

4、f1*t+pi/180*w1)+x2*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;%原信號描點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);xl

5、im(0 N-1);grid;title(采樣信號);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進行采樣FFT2(1,1,1,90,2,2,18

6、0,3,3,180,32)FFT2(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ù)時除卻k=0和N/2外（N為奇數(shù)時，除卻k=0外），圖形呈中心對稱分布。理論上,由于復指數(shù)的周期性，長度為N（假設N為偶數(shù)，N為奇數(shù)時類似）的的DFT頻譜分析，可把k=N/2,N/2+1,N-1看作是負頻

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

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

9、FFT2(2,1,1,90,2,2,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處相位值不變。 鍵入：%改變第三個余弦成分的初相位（由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ā)生相應的變化，即恰合頻譜的物

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

11、793;int M,N;/*定義復數(shù)*/struct Complex_double real;double img;*a,*b;/*定義2的冪計算*/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用二進制表示，然后碼位倒讀*/i-=(int)(i/x_2(j)*x_2(j);return(ii);/*定義復數(shù)乘法*/struct Complex_ Multi(struct Complex_ a,struct Complex_ b)struct Complex_ c;c.real=a.rea

12、l*b.real-a.img*b.img;c.img=a.real*b.img+a.img*b.real;return(c);/*定義復數(shù)冪運算*/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);/*定義復數(shù)加法*/struct Complex_ Add(struct Complex_ a,struct Complex_ b)struct

13、Complex_ c;c.real=a.real+b.real;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)，采點數(shù)*/printf(采點數(shù):);scanf(%d,&NN);printf(n);/*信號的直流成分*/printf(信號的直流成分:n);scanf(%f,&a0);/*信號有三個余弦成分*/*余弦成分1*/printf(第1個正弦成分的頻率:n);scanf(%f,&f1);printf(幅值:n);scanf(

14、%f,&a1);printf(初相位：n);scanf(%f,&w1);/*余弦成分2*/printf(第2個正弦成分的頻率:n);scanf(%f,&f2);printf(幅值:n);scanf(%f,&a2);printf(初相位：n);scanf(%f,&w2);/*余弦成分3*/printf(第3個正弦成分的頻率: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 Com

15、plex_*)calloc(N,sizeof(struct Complex_);/*動態(tài)開辟N個單位*/b=(struct Complex_*)calloc(N,sizeof(struct Complex_);/*動態(tài)開辟N個單位*/*把采樣信號用碼位倒讀的方法存入計算機中*/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).im

16、g=0;/*動態(tài)空間空閑處，補零*/for(i=NN;iN;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級運算*/n=n/2;m=m*2;q=x_2(i+1);/*q=2(i+1)*/W.real=cos(-2*pi/q);/*旋轉(zhuǎn)因子實部*/W.img=sin(-2*pi/q);/*旋轉(zhuǎn)因子虛部*/

17、for(j=0;jn;j+)/*n群運算*/ for(k=0;km/2;k+)Wk=W_K(W,k); bj*m+k=Add(aj*m+k,Multi(Wk,aj*m+k+m/2);/*Add()為復數(shù)相加函數(shù)，Multi()為復數(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()為復數(shù)相加函數(shù)，Multi()為復數(shù)相乘函數(shù)*/ for(r=0;rN;r+) ar=br; /*將結(jié)果寫入文件e:keshe.txt*/if(fp=fopen(e:keshe.txt,w+)=NUL

18、L)printf(Cannot open file strike 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ù)據(jù)將由keshe.txt調(diào)入到MATLAB中，做出頻譜圖。 參數(shù)設置1圖17 fs=32Hz(C實現(xiàn))圖18 fs=32Hz（MATLAB仿真）注意到兩

19、圖的幅度譜圖，很相似，而相位譜卻相差很大。用C實現(xiàn)的FFT得到的相位譜也似乎沒有什么對稱性，但可以看到兩個相位譜在k=0,1,2,3和31,30,29處的相位值是基本一致的，而其他點處，相位值差異則很大。（表1，列出了具體的數(shù)值計算結(jié)果）表1kFFT（C）FFT（MATLAB）FFT結(jié)果相位值FFT結(jié)果相位值032.0000000000000 + 0.00000000000000i032.0000000000000 + 0.00000000000000i010.00000000000000 + 16.0000000000000i1.5707962.99760216648792e-15 - 16

20、.0000000000000i1.5707962-32.0000000000000 + 1.00000000000000e-14i3.141593-32.0000000000000 - 9.76996261670138e-15i3.1415933-48.0000000000000 + 3.00000000000000e-14i3.141593-48.0000000000000 - 4.04121180963557e-14i3.1415934-1.00000000000000e-14 + 0.00000000000000i3.141593-1.16001965116118e-14 - 9.420

21、55475210265e-16i3.0605650.00000000000000 + 0.00000000000000i0-3.99680288865056e-15 - 1.77635683940025e-15i2.72336861.00000000000000e-14 + 1.00000000000000e-14i0.7853987.54951656745106e-15 - 4.88498130835069e-15i0.57430571.00000000000000e-14 + 1.00000000000000e-14i0.7853981.11022302462516e-15 - 5.218

22、04821573824e-15i1.3611568-1.00000000000000e-14 + 1.00000000000000e-14i2.356194-1.06581410364015e-14 - 7.10542735760100e-15i2.5535990.00000000000000 - 1.00000000000000e-14i-1.57081.11022302462516e-15 + 7.21644966006352e-15i-1.41815100.00000000000000 + 0.00000000000000i0-2.22044604925031e-15 - 2.22044

23、604925031e-15i2.356194110.00000000000000 + 0.00000000000000i01.11022302462516e-14 - 5.32907051820075e-15i0.4475212-1.00000000000000e-14 + 0.00000000000000i3.141593-9.71608556119124e-15 - 9.42055475210265e-16i3.044936130.00000000000000 + 0.00000000000000i03.55271367880050e-15 - 1.33226762955019e-15i0

24、.35877114-1.00000000000000e-14 + 0.00000000000000i3.1415933.55271367880050e-15 + 0.00000000000000i015-2.00000000000000e-14 + 1.00000000000000e-14i2.677945-1.58761892521397e-14 - 1.15463194561016e-14i2.512796161.00000000000000e-14 + 0.00000000000000i0-3.55271367880050e-15 + 0.00000000000000i3.1415931

25、7-2.00000000000000e-14 - 1.00000000000000e-14i-2.67795-1.58761892521397e-14 + 1.15463194561016e-14i-2.512818-1.00000000000000e-14 + 0.00000000000000i3.1415933.55271367880050e-15 + 0.00000000000000i019-1.00000000000000e-14 + 0.00000000000000i3.1415933.55271367880050e-15 + 1.33226762955019e-15i-0.3587

26、720-1.00000000000000e-14 + 0.00000000000000i3.141593-9.71608556119124e-15 + 9.42055475210265e-16i-3.04494211.00000000000000e-14 + 0.00000000000000i01.11022302462516e-14 + 5.32907051820075e-15i-0.44752220.00000000000000 + 0.00000000000000i0-2.22044604925031e-15 + 2.22044604925031e-15i-2.35619230.0000

27、0000000000 + 1.00000000000000e-14i1.5707961.11022302462516e-15 - 7.21644966006352e-15i1.41814724-1.00000000000000e-14 - 1.00000000000000e-14i-2.35619-1.06581410364015e-14 + 7.10542735760100e-15i-2.55359250.00000000000000 - 1.00000000000000e-14i-1.57081.11022302462516e-15 + 5.21804821573824e-15i-1.36

28、116261.00000000000000e-14 - 1.00000000000000e-14i-0.78547.54951656745106e-15 + 4.88498130835069e-15i-0.574327-1.00000000000000e-14 + 0.00000000000000i3.141593-3.99680288865056e-15 + 1.77635683940025e-15i-2.7233728-1.00000000000000e-14 + 0.00000000000000i3.141593-1.16001965116118e-14 + 9.420554752102

29、65e-16i-3.0605629-48.0000000000000 - 5.00000000000000e-14i-3.14159-48.0000000000000 + 4.04121180963557e-14i-3.1415930-32.0000000000000 - 2.00000000000000e-14i-3.14159-32.0000000000000 + 9.76996261670138e-15i-3.14159311.00000000000000e-14 - 16.0000000000000i-1.57082.99760216648792e-15 + 16.0000000000

30、000i-1.5708結(jié)合圖17、18和表1，可以發(fā)現(xiàn)兩個相位譜在k=0,1,2,3和31,30,29處一致，而在其他點相差很大，注意到這些點的幅值很?。ɡ碚撋蠎摓?，實際上卻由于計算機的誤差，而使之非0），可認為為0，可忽略這些點，而主要考慮主頻率成分。事實上，結(jié)合圖表可發(fā)現(xiàn)，C程序所得的結(jié)果也并非完全無對稱性的，并且可注意到所有不對稱點的相位值接近圓周率或0。這種差異是由于計算機的誤差導致的，在一般情況下，這不會造成太大的影響，而在相位處于邊界值（和0）時，這種誤差將計算結(jié)果發(fā)生很大變化。注：在先前編寫的C程序中，把復數(shù)實部和虛部的數(shù)據(jù)類型設置為float時，相位譜的不對稱更加嚴重，甚至

31、C和MATLAB仿真的兩個相位譜在k= 31,30,29都不一致，此處實驗結(jié)果是把數(shù)據(jù)類型改為double類型后（即提高精度后）所得的，仍存在這種不對稱，理論上可通過不斷提高精度進行改進，但只要是某頻率成分初相位為或0，就存在這種嚴重不對稱的風險。 參數(shù)設置2圖19 fs=16Hz(C實現(xiàn)) 圖20 fs=16Hz（MATLAB仿真） 可發(fā)現(xiàn)C實現(xiàn)下的FFT幅度譜圖隨采樣頻率變化的趨勢與MATLAB仿真實驗類似，而由于參數(shù)1實驗相同的原因，C實現(xiàn)的FFT相位譜并不對稱。 參數(shù)設置3圖21 fs=8Hz(C實現(xiàn)) 圖22 fs=8Hz（MATLAB仿真）類似的，可發(fā)現(xiàn)幅度譜圖隨采樣頻率變化的趨勢

32、與MATLAB仿真實驗類似。而C和MATLAB實現(xiàn)下，k=6時的相位值一個為3.14，一個為-3.14，具體原因同上。而和相差，理論上初相位為是無差異的，但是體現(xiàn)在相位譜上就是造成不對稱。 參數(shù)設置4圖23 fs=4Hz(C實現(xiàn))圖24 fs=4Hz（MATLAB仿真）同MATLAB仿真結(jié)果一樣出現(xiàn)混疊現(xiàn)象。 參數(shù)設置5圖25 fs=14Hz(C實現(xiàn)) 圖26 fs=14Hz（MATLAB仿真）注意到幅度譜和相位譜都有很大區(qū)別，其中相位譜的區(qū)別是由兩種不同的計算機語言規(guī)則造成的（正如之前提到的），而幅度譜的區(qū)別是由算法不同造成的， C程序是采用補零的方法把序列長度湊成，然后把長度為N的序列輸出

33、為結(jié)果，這樣做的后果是改變了進行DFT的序列， 得到的自然不是所求的結(jié)果，所以這個C程序，只能針對序列長度滿足的情況。v 改變被采樣信號波形 參數(shù)設置6圖27 fs=32Hz(C實現(xiàn))圖28 fs=32Hz(MATLAB仿真)二者幅度譜相同，相位譜只在中間部分有差異（此差異對實際分析影響很小，且也無實際意義，可忽略），而在k=0,1,2,3和31,30,29處兩幅圖的情況一致，故驗證了之前關(guān)于出現(xiàn)圖17和圖18二者相位譜如此大差異原因的猜想。 參數(shù)設置7圖29 fs=32Hz(C實現(xiàn))圖30 fs=32Hz(MATLAB仿真)二者幅度譜圖相同，且較上一參數(shù)時發(fā)生的變化趨勢與MATLAB仿真試驗

34、中一樣；相位譜圖在k=0,1,2,3和31,30,29處相位值相同。 參數(shù)設置8圖31 fs=32Hz(C實現(xiàn))圖32 fs=32Hz(MATLAB仿真)二者幅度譜圖相同，且較上一參數(shù)時發(fā)生的變化趨勢與MATLAB仿真試驗中一樣；相位譜圖在k=0,1,2,4和31,30,28處相位值相同。至此，考察了序列長度和信號諧波成分對MATLAB的FFT仿真和C實現(xiàn)的FFT的影響，并結(jié)合了DFT的物理意義，對實驗現(xiàn)象做出了解釋。對比MATLAB和C實現(xiàn)下的FFT結(jié)果的差異，做出了自己的解釋。4. 課程設計心得與自我評價這個課程設計的主要部分是在暑假完成的，課程設計是關(guān)于FFT的計算機實現(xiàn)，要求掌握FFT

35、算法，并分別用MATLAB進行FFT仿真和編寫實現(xiàn)FFT的C代碼。記得之前學習電路理論、復變函數(shù)、數(shù)理方程和信號系統(tǒng)時都遇到過把時域上的問題轉(zhuǎn)化為頻域上的問題解決的情況，電路中主要是Laplace變換的實際應用，而復變函數(shù)里面主要講的是各種變換的技術(shù)新性問題，數(shù)理方程里面研究的是一些數(shù)學方程的物理模型，而信號與系統(tǒng)課程則引入了離散的概念。首先是時域上的離散，由對連續(xù)時間進行信號采樣，即x(t)xn，把連續(xù)時間信號的各種性質(zhì)在離散時間信號做推廣，如CTFTDTFT，Laplace變換Z-變換；其次是頻域上的離散，對連續(xù)頻率進行信號采樣，即DTFTDFT，而FFT則是處理DFT的一種快速算法。課程設計完成的過程中遇到許多問題：首先是M