快速傅里葉變換

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

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

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

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

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

6、行實驗一中的 DFT程序，記錄運行時間，比較它與FFT的速度差距。四、C語言參考程序(1)用FFT分析正弦信號頻譜的程序#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;/* 計算FFT的級數(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)進行碟形計算*/n1=1;for(l=1;l<=m;l+)/*分級循環(huán) */(n1=2*n1;n2=n1/2;e=3n2;for(k=0;k<n2;k+) /*分組循環(huán) */(c=cos(k*e); /*計算 W */s=-sin(k*e);(for(i=k;i<n;i+=n1) /*碟形循環(huán)*/j=i+n2;tr=c*xj-s*yj; /*計算碟形 */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)用子程序計算FFT */fft(x,y,N);/*計算FFT的幅度*/for(k=0;k<N;k+)wk=sqr

10、t(xk*xk+yk*yk);/*畫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)測量FFT運行時間的程序FFT程序中以下給出的是用于測量FFT運算時間的有關(guān)程序語句，應(yīng)將它們插入前面main()函數(shù)內(nèi)適當(dāng)?shù)奈恢谩nsigned int start,end;float time;int i; /*提示：這兩行放在主程序開頭數(shù)據(jù)類型說明部分/*提示：以下部分替代主程序中的fft(x,y,N)start=biostime(0,0);

13、/*for (n=0;n<2000;n+) /*讀取開始時間*/循環(huán)計算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é)束時間*/time=(end-start)/2000.0*55; /* biostime()讀取的一個時間單位約為55ms */printf("It took %f msn", time); /*屏幕輸出FFT的運算時間*/* draw(N); */ /*不運行作圖語句*/(3)測量DFT運行時間的程序DFT程序中以下給出的是用于測量DFT運算時間的有關(guān)程序語句，應(yīng)將它們插入實驗main()函數(shù)內(nèi)適當(dāng)?shù)奈恢谩nsigned int start,end;float time;int j;start=biostime(0,0); /*for(j=0;j<50;j+)(/*讀取開始時間*/循環(huán)計算 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é)束時間*/t