快速傅里葉變換

1、實(shí)驗(yàn).快速傅立葉變換一、實(shí)驗(yàn)?zāi)康? .學(xué)習(xí)和掌握快速傅立葉變換（FFT）的實(shí)現(xiàn)過(guò)程和編程技術(shù)2. 運(yùn)用FFT分析正弦信號(hào)的頻譜3. 測(cè)試FFT的運(yùn)算時(shí)間，比較 FFT與DFT的運(yùn)算速度，獲得對(duì) FFT “快速”的感性認(rèn)識(shí)。4. 鍛煉和提高數(shù)字信號(hào)處理的程序設(shè)計(jì)和調(diào)試能力。二、實(shí)驗(yàn)原理與方法FFT并不是與DFT不同的另一種變換，而是為了減少DFT運(yùn)算次數(shù)的一種快速算法。它是對(duì) DFT變換式進(jìn)行一次次分解，使其成為若干小點(diǎn)數(shù)的組合，從而減少運(yùn)算量o常用的FFT是以2為基數(shù)的，其長(zhǎng)度 N=2%它的效率高，程序簡(jiǎn)單，使用非常方便，當(dāng)要變換的序列長(zhǎng)度不等于2的整數(shù)次方時(shí)， 為了使用以2為基數(shù)的FFT,可

2、以用末尾補(bǔ)零的方法，使其長(zhǎng)度延長(zhǎng)至 2的整數(shù)次方。本實(shí)驗(yàn)運(yùn)用時(shí)間抽取基2 FFT,其原理、信號(hào)流圖和運(yùn)算過(guò)程可參閱課堂筆記、教材和其它教科書(shū)。FFT的實(shí)現(xiàn)要比DFT復(fù)雜，通常采用三個(gè)嵌套循環(huán)來(lái)實(shí)現(xiàn)。最外面的循環(huán)是分級(jí)循環(huán)，N=2M點(diǎn)的FFT分為M級(jí)計(jì)算。中間一層是分組循環(huán)，一級(jí)內(nèi)蝶形系數(shù)Wk相同的蝶形構(gòu)成一組。最內(nèi)層為蝶形計(jì)算的循環(huán)，一組內(nèi)不同輸入數(shù)據(jù)的蝶形逐個(gè)計(jì)算。編程時(shí)要注意各級(jí)蝶形組之間的間隔、組內(nèi)蝶形之間間隔、以及系數(shù)Wk變化。有不少書(shū)中有 FFT程序可供參考，但要注意理解和弄懂，不可一味照搬，要盡量自己去編。三、實(shí)驗(yàn)內(nèi)容1. 設(shè)計(jì)說(shuō)明編制時(shí)間抽取基 2 FFT程序計(jì)算前面 DFT程序

3、分析過(guò)的正弦信號(hào)，與 DFT計(jì)算的結(jié)果進(jìn)行比較，以驗(yàn)證所編程序的正確性。其中一個(gè)用于FFT為復(fù)數(shù)運(yùn)算，若用實(shí)數(shù)運(yùn)算來(lái)實(shí)現(xiàn)，需要設(shè)置兩個(gè)數(shù)組來(lái)存放輸入輸出數(shù)據(jù)。存放數(shù)據(jù)的實(shí)部，另一個(gè)存放數(shù)據(jù)的虛部。正弦輸入信號(hào)為實(shí)數(shù)，則令其虛部為零。同樣，碟形運(yùn)算也要化成實(shí)數(shù)來(lái)進(jìn)行，分別算出實(shí)部和虛部。當(dāng)然也可以直接用復(fù)數(shù)數(shù)組和語(yǔ)句實(shí)現(xiàn)。程序設(shè)計(jì)的難點(diǎn)和重點(diǎn)在于要合理安排和正確設(shè)置碟形運(yùn)算的分級(jí)循環(huán)、級(jí)內(nèi)分組循環(huán)和組內(nèi)分碟形循環(huán)。要注意乘法系數(shù)時(shí)的變化以及各循環(huán)變量的變化和調(diào)整。程序中，在進(jìn)行 FFT計(jì)算之前，先要將輸入數(shù)據(jù)排列成二進(jìn)制倒序的形式，這可通過(guò)將有關(guān)數(shù) 據(jù)相互換位來(lái)實(shí)現(xiàn)。正弦抽樣信號(hào)的產(chǎn)生和輸出頻譜

4、圖的繪制與前面的DFT程序相同。2. 實(shí)驗(yàn)步驟(1) 將編制的程序輸入計(jì)算機(jī)，調(diào)試、運(yùn)行正確。即輸入實(shí)驗(yàn)一中的12組數(shù)據(jù)，看FFT輸出譜線的結(jié)果是否與前面DFT的結(jié)果一致。(2) 在不同的FFT長(zhǎng)度下(可分別令 N=64, 128, 256, 512, 1024)，運(yùn)行 FFT程序，觀察、比較和記錄運(yùn)行時(shí)間。由于現(xiàn)在計(jì)算機(jī)的運(yùn)算速度很快，計(jì)算一次FFT的時(shí)間很短(在ms數(shù)量級(jí))，不便觀察和記錄?？稍诔绦蛑刑碇醚h(huán)語(yǔ)句，讓其重復(fù)計(jì)算FFT許多次(可重復(fù) 100次、1000次或更多，視情況而定)，使總的運(yùn)算時(shí)間達(dá)到數(shù)秒或數(shù)十秒鐘，從而可以通過(guò)人工來(lái)觀察和記錄。此時(shí)注意，為了防止 FFT多次循環(huán)計(jì)算

5、，使輸入輸出數(shù)據(jù)越來(lái)越大而發(fā)生溢出，可令輸入數(shù)據(jù)全部為零。這樣，F(xiàn)FT的運(yùn)算量保持不變，但輸入輸出始終為零，不會(huì)因多次循環(huán)而發(fā)生溢出。在觀測(cè)FFT的運(yùn)算時(shí)間時(shí)，作圖部分的語(yǔ)句不要運(yùn)行(一方面，它不應(yīng)該計(jì)入FFT或DFT的運(yùn)算時(shí)間。另一方面，當(dāng)N較大時(shí)，原有程序中的作圖語(yǔ)句會(huì)出錯(cuò))，可將其注釋掉。測(cè)量FFT運(yùn)行時(shí)間的一個(gè)更好、更精確的方法是在程序中調(diào)用讀取計(jì)算機(jī)時(shí)鐘的函數(shù)。TURBOC中提供多種讀取計(jì)算機(jī)時(shí)間的方法和函數(shù)，一種比較簡(jiǎn)單的方法是用biostime() 函數(shù),后面給出了相應(yīng)的參考程序。(3) 用上述類似的方法，在不同的長(zhǎng)度下(N=64, 128, 256 , 512 ,1024),運(yùn)

6、行實(shí)驗(yàn)一中的 DFT程序，記錄運(yùn)行時(shí)間，比較它與FFT的速度差距。四、C語(yǔ)言參考程序(1)用FFT分析正弦信號(hào)頻譜的程序#include<math.h>void fft(x,y,n) /* FFT作為子程序 */int n;float x1024,y1024;int i,j,k,l,m,n1,n2;float c,s,e,tr,ti;/* 計(jì)算FFT的級(jí)數(shù) M */for(j=1,i=1;i<n;i+)m=i;j=2*j;if(j=n) break;/*改變輸入數(shù)據(jù)的順序*/n1=n-1;for(j=0,i=0;i<n1;i+)if(i<j)tr=xj;ti=yj

7、;xj=xi;yj=yi;xi=tr;yi=ti；)k=n/2;while(k<(j+1)(j=j-k;k=k/2;)j=j+k;)/*三重嵌套循環(huán)進(jìn)行碟形計(jì)算*/n1=1;for(l=1;l<=m;l+)/*分級(jí)循環(huán) */(n1=2*n1;n2=n1/2;e=3n2;for(k=0;k<n2;k+) /*分組循環(huán) */(c=cos(k*e); /*計(jì)算 W */s=-sin(k*e);(for(i=k;i<n;i+=n1) /*碟形循環(huán)*/j=i+n2;tr=c*xj-s*yj; /*計(jì)算碟形 */ti=c*yj+s*xj;xj=xi-tr;

8、yj=yi-ti;xi=xi+tr;yi=yi+ti;)/*主程序部分*/#include<stdio.h>#include<math.h>#include<conio.h>#include<graphics.h>#define PI 3.1415926float x1024,y1024,w1024;void draw(int);void axis(int,int);main()int N,f,n,k;float T;/*鍵盤輸入f、N、T */printf("the frequency of the sine wave f="

9、;);scanf("%d",&f);printf("the number of samples N=");scanf("%d",&N);printf("the time step of samples T=");scanf("%f",&T);/*產(chǎn)生正弦輸入數(shù)據(jù)*/for (n=0;n<N;n+)xn=sin(2*PI*f*n*T);yn=0;/*調(diào)用子程序計(jì)算FFT */fft(x,y,N);/*計(jì)算FFT的幅度*/for(k=0;k<N;k+)wk=sqr

10、t(xk*xk+yk*yk);/*畫(huà)FFT的幅度譜*/draw(N);getch();closegraph();(/*作圖子程序*/void draw(int N)(int k,driver,mode,step,x,y;driver=VGA;mode=VGAHI;step=512/N;x=64;y=400;registerbgidriver(driver);initgraph(&driver,&mode,"c:tc");axis(x,y);outtextxy(x+step+546),y+8,"k”);outtextxy(x+step-8),y-29

11、0,"x(k)");outtextxy(x+step-8),y-234,"32");outtextxy(x+step-8),y-121,"16");outtextxy(x+step-19),y+8,"0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1");for(k=0;k<N;k+)(line(x,y,x,(y-floor(wk*7);x+=step;)void axis(int x,int y)line(10,y,630,y);l

12、ine(610,(y-4),630,y);line(610,(y+4),630,y);line(x,(y+20),x,(y-280);line(x-4),(y-260),64,120);line(x+4),(y-260),64,120);)(2)測(cè)量FFT運(yùn)行時(shí)間的程序FFT程序中以下給出的是用于測(cè)量FFT運(yùn)算時(shí)間的有關(guān)程序語(yǔ)句，應(yīng)將它們插入前面main()函數(shù)內(nèi)適當(dāng)?shù)奈恢?。unsigned int start,end;float time;int i; /*提示：這兩行放在主程序開(kāi)頭數(shù)據(jù)類型說(shuō)明部分/*提示：以下部分替代主程序中的fft(x,y,N)start=biostime(0,0);

13、/*for (n=0;n<2000;n+) /*讀取開(kāi)始時(shí)間*/循環(huán)計(jì)算FFT 2000次*/fft(x,y,N);for (i=0;i<N;n+)/*清除x,y的值，防止下一次的FFT數(shù)據(jù)溢出(xi=0;yi=0;)end=biostime(0,0); /*讀取結(jié)束時(shí)間*/time=(end-start)/2000.0*55; /* biostime()讀取的一個(gè)時(shí)間單位約為55ms */printf("It took %f msn", time); /*屏幕輸出FFT的運(yùn)算時(shí)間*/* draw(N); */ /*不運(yùn)行作圖語(yǔ)句*/(3)測(cè)量DFT運(yùn)行時(shí)間的程序DFT程序中以下給出的是用于測(cè)量DFT運(yùn)算時(shí)間的有關(guān)程序語(yǔ)句，應(yīng)將它們插入實(shí)驗(yàn)main()函數(shù)內(nèi)適當(dāng)?shù)奈恢?。unsigned int start,end;float time;int j;start=biostime(0,0); /*for(j=0;j<50;j+)(/*讀取開(kāi)始時(shí)間*/循環(huán)計(jì)算 DFT 50次*/for(k=0;k<N;k+)(r=i=0.0;for(n=0;n<N;n+)(r=r+xn*cos(c*n*k);i=i+xn*sin(c*n*k);)yk=r;wk=i;)end=biostime(0,0); /*讀取結(jié)束時(shí)間*/t