數(shù)字信號處理課程設(shè)計(jì)-DFT的快速算法——快速傅里葉變換(FFT)

1、目錄前言第一章 離散傅里葉變換DFT31.1 DFT定義31.2 DFT的快速算法快速傅里葉變換(FFT) 3第二章 基2 DIT-FFT算法42.1 按時域抽取的基2 DIT-FFT算法4第三章 基于C語言設(shè)計(jì)16點(diǎn)基2DIT-FFT程序及運(yùn)行結(jié)果 63.1 按時間抽取的基2DIT-FFT程序 63.2 程序運(yùn)行結(jié)果 15第四章 課程設(shè)計(jì)的總結(jié) 17參考文獻(xiàn) 17 前 言信號(signal)是一種物理表達(dá)，或是傳遞信息的函數(shù)。而信息是信號的具體內(nèi)容。模擬信號(analog signal)：指時間連續(xù)、幅度連續(xù)的信號。數(shù)字信號(digital signal)：時間和幅度上都是離散(量化)的信號

2、。 數(shù)字信號可用一序列的數(shù)表示，而每個數(shù)又可表示為二制碼的形式，適合計(jì)算機(jī)處理。數(shù)字信號處理是將信號以數(shù)字方式表示并處理的理論和技術(shù)。數(shù)字信號處理與模擬信號處理是信號處理的子集。 其目的是對真實(shí)世界的連續(xù)模擬信號進(jìn)行測量或?yàn)V波。因此在進(jìn)行數(shù)字信號處理之前需要將信號從模擬域轉(zhuǎn)換到數(shù)字域，這通常通過模數(shù)轉(zhuǎn)換器實(shí)現(xiàn)。而數(shù)字信號處理的輸出經(jīng)常也要變換到模擬域，這是通過數(shù)模轉(zhuǎn)換器實(shí)現(xiàn)的。數(shù)字信號處理的核心算法是離散傅立葉變換(DFT)，是DFT使信號在數(shù)字域和頻域都實(shí)現(xiàn)了離散化，從而可以用通用計(jì)算機(jī)處理離散信號。而使數(shù)字信號處理從理論走向?qū)嵱玫氖强焖俑盗⑷~變換(FFT)，F(xiàn)FT的出現(xiàn)大大減少了DFT的運(yùn)

3、算量，使實(shí)時的數(shù)字信號處理成為可能、極大促進(jìn)了該學(xué)科的開展。隨著大規(guī)模集成電路以及數(shù)字計(jì)算機(jī)的飛速開展，加之從60年代末以來數(shù)字信號處理理論和技術(shù)的成熟和完善，用數(shù)字方法來處理信號，即數(shù)字信號處理，已逐漸取代模擬信號處理。 隨著信息時代、數(shù)字世界的到來，數(shù)字信號處理已成為一門極其重要的學(xué)科和技術(shù)領(lǐng)域。第一章 離散傅里葉變換DFT離散傅立葉變換DFT實(shí)現(xiàn)了信號首次在頻域表示的離散化，使得頻域也能夠用計(jì)算機(jī)進(jìn)行處理。并且這種DFT變換可以有多種實(shí)用的快速算法。使信號處理在時、頻域的處理和轉(zhuǎn)換均可離散化和快速化，因而具有重要的理論意義和應(yīng)用價值。1.1 DFT定義設(shè)序列x(n)長度為M，定義x(n)

4、的N點(diǎn)DFT為式中，N稱為離散傅里葉變換區(qū)間長度，要求N M。為書寫簡單，令 ，因此通常將N點(diǎn)DFT表示為1.2 DFT的快速算法快速傅里葉變換(FFT)DFT使計(jì)算機(jī)在頻域處理信號成為可能，但是當(dāng)N很大時，直接計(jì)算N點(diǎn)DFT的計(jì)算量非常大。快速傅里葉變換FFT，F(xiàn)ast Fourier Transform可使實(shí)現(xiàn)DFT的運(yùn)算量下降幾個數(shù)量級，從而使數(shù)字信號處理的速度大大提高，工程應(yīng)用成為可能。人們已經(jīng)研究出多種FFT算法，它們的復(fù)雜度和運(yùn)算效率各不相同?？焖俑道锶~變換就是不斷地將長序列的DFT分解為短序列的DFT，并利用 的周期性和對稱性及其一些特殊值來減少DFT運(yùn)算量的快速算法。FFT算法

5、分類：1.時間抽選法 DIT：Decimation-In-Time2.頻率抽選法 DIF：Decimation-In-Frequency時間域抽?。夯?時間抽取(Decimation in time)DIT-FFT算法頻率域抽?。夯?頻率抽取(Decimation in frequency)DIF-FFT算法第二章 基2 DIT-FFT算法2.1 按時間抽取的基2 DIT-FFT算法：1、按時間抽取的基2 DIT-FFT算法原理先設(shè)序列點(diǎn)數(shù)為N=2M，M為整數(shù)。如果不滿足這個條件，可以人為地加上假設(shè)干零值點(diǎn)，使之到達(dá)這一要求。這種N為2的整數(shù)冪的FFT稱基-2 FFT。設(shè)輸入序列長度為N=2M

6、 (M為正整數(shù)) ，將該序列按時間順序的奇偶分解為越來越短的子序列，稱為按時間抽取(DIT )的FFT算法。 序列x(n)的N點(diǎn)DFT為 將上面的和式按n的奇偶性分解為 令x1(l)=x(2l),x2(l)=x(2l+1)，因?yàn)?，所以上式可寫成上式說明，按n的奇偶性將x(n)分解為兩個N/2長的序列x1(l)和x2(l)，那么N點(diǎn)DFT可分解為兩個N/2點(diǎn)DFT來計(jì)算。用X1(k)和X2(k)分別表示x1(l)和x2(l)的N/2點(diǎn)DFT，即有上述公式，及X1k、X2(k)的隱含周期性得到：這樣，就將N點(diǎn)DFT的計(jì)算分解為計(jì)算兩個N/2點(diǎn)離散傅里葉變換X1(k)和X2(k)，再計(jì)算上式。蝶形

7、圖：2、 按時間抽選的FFT算法的特點(diǎn)：1 原位運(yùn)算由圖4.2.4可以看出，DIT-FFT的運(yùn)算過程很有規(guī)律。N=2M點(diǎn)的FFT共進(jìn)行M級運(yùn)算，每級由N/2個蝶形運(yùn)算組成。同一級中，每個蝶形的兩個輸入數(shù)據(jù)只對計(jì)算本蝶形有用，而且每個蝶形的輸入、輸出數(shù)據(jù)結(jié)點(diǎn)又同在一條水平線上，這就意味著計(jì)算完一個蝶形后，所得輸出數(shù)據(jù)可立即存入原輸入數(shù)據(jù)所占用的存儲單元。這樣，經(jīng)過M級運(yùn)算后，原來存放輸入序列數(shù)據(jù)的N個存儲單元中便依次存放 的N個值。這種利用同一存儲單元存儲蝶形計(jì)算輸入、輸出數(shù)據(jù)的方法稱為原位計(jì)算。原位計(jì)算可節(jié)省大量內(nèi)存，從而使設(shè)備本錢降低。2倒序規(guī)律按原位計(jì)算時，F(xiàn)FT的輸出 是按正常

8、順序排列在存儲單元中，但是這時輸入x(n)卻不是、按自然順序存儲的，而是按x(0)，x(4)，x(7)的順序存入存儲單元，看起來好象是“混亂無序的，實(shí)際上是有規(guī)律的，我們稱之為倒序。第三章 基于C語言設(shè)計(jì)16點(diǎn)基2DFT-FFT程序及運(yùn)行結(jié)果設(shè)計(jì)參數(shù)及要求：本次試驗(yàn)要求利用C語言設(shè)計(jì)16點(diǎn)基2DIT-FFT程序，并對如下程序進(jìn)行分析，給出輸出結(jié)果并畫出信號的幅值圖和相位譜。n012345678910X(n)01.222.232.85 2.982.601.760.62-0.62-1.76-2.601112131415-2.98-2.85-2.23-1.2203.1 按時間抽取的基2DFT-FFT

9、程序：/*時間抽選基2 FFT算法C語言實(shí)現(xiàn)*/#include <stdio.h>#include <math.h>#include <stdlib.h>#include <conio.h>#define N 1000#include<graphics.h>/*定義復(fù)數(shù)類型*/typedef struct double real; double img;complex;complex xN, *W; /*輸入序列,變換核*/int size_x=0; /*輸入序列的大小，在本程序中僅限2的次冪*/double PI; /*圓周率*/

10、int z;float wN,pN;char strN;/*初始化變換核*/void initW() int i; W=(complex *)malloc(sizeof(complex) * size_x); for(i=0;i<size_x;i+) Wi.real=cos(2*PI/size_x*i); Wi.img=-1*sin(2*PI/size_x*i); /*加法*/void add(complex a,complex b,complex *c) c->real=a.real+b.real; c->img=a.img+b.img; /*乘法*/void mul(co

11、mplex a,complex b,complex *c) c->real=a.real*b.real - a.img*b.img; c->img=a.real*b.img + a.img*b.real; /*減法*/void sub(complex a,complex b,complex *c) c->real=a.real-b.real; c->img=a.img-b.img; /*除法*/void divi(complex a,complex b,complex *c) c->real=( a.real*b.real+a.img*b.img )/( b.re

12、al*b.real+b.img*b.img); c->img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img); /*變址計(jì)算，將x(n)碼位倒置*/void change() complex temp; unsigned short i=0,j=0,k=0; double t; for(i=0;i<size_x;i+) k=i;j=0; t=(log(size_x)/log(2); while( (t-)>0 ) j=j<<1; j|=(k & 1); k=k>>1; if(j&

13、gt;i) temp=xi; xi=xj; xj=temp; /*快速傅里葉變換*/void fft() int i=0,j=0,k=0,l=0; complex up,down,product; change(); for(i=0;i< log(size_x)/log(2) ;i+) /*一級蝶形運(yùn)算*/ l=1<<i; for(j=0;j<size_x;j+= 2*l ) /*一組蝶形運(yùn)算*/ for(k=0;k<l;k+) /*一個蝶形運(yùn)算*/ mul(xj+k+l,Wsize_x*k/2/l,&product); add(xj+k,product,

14、&up); sub(xj+k,product,&down); xj+k=up; xj+k+l=down; /*快速傅里葉逆變換*/void ifft() int i=0,j=0,k=0,l=size_x; complex up,down; for(i=0;i< (int)( log(size_x)/log(2) );i+) /*一級蝶形運(yùn)算*/ l/=2; for(j=0;j<size_x;j+= 2*l ) /*一組蝶形運(yùn)算*/ for(k=0;k<l;k+) /*一個蝶形運(yùn)算*/ add(xj+k,xj+k+l,&up); up.real/=2;u

15、p.img/=2; sub(xj+k,xj+k+l,&down); down.real/=2;down.img/=2; divi(down,Wsize_x*k/2/l,&down); xj+k=up; xj+k+l=down; change();/*輸出傅里葉變換的結(jié)果和幅值結(jié)果*/void output() int i; clrscr() ; printf("The result are as followsnn"); for(i=0;i<size_x;i+) printf("%.4f",xi.real); if(xi.img&g

16、t;=0.0001)printf("+%.4fjn",xi.img); else if(fabs(xi.img)<0.0001)printf("n"); else printf("%.4fjn",xi.img); for( z=0;z<16;z+) wz=sqrt(xz.real)*(xz.real)+(xz.img)*(xz.img); printf("nnThe Amplitudes are as followsnn"); for(z=0;z<size_x;z+) printf("%

17、.4fn",wz); /*輸出傅里葉變換相位的結(jié)果*/void output2() int i; setbkcolor(0); for(i=0;i<size_x;i+) for( z=0;z<16;z+) pz=atan(xz.img)/(xz.real); printf("nnnThe Phase are as followsnn"); for(z=0;z<size_x;z+) setcolor(1); printf("%.4fn",pz); void draw(int n) /* 作幅值圖 */ int k,driver,

18、mode,step,x,y; int i,j,xstep=512/n,ystep=55; char buffer10,bufferaN; step=512/n; driver=VGA; mode=VGAHI; x=64; y=400; initgraph(&driver,&mode,""); setbkcolor(15); /* 繪制坐標(biāo)及刻度 */ setcolor(1); moveto(200,30); outtext("The Figure of Amplitude is as follow!"); line(x,80,x,400)

19、; /* 畫縱坐標(biāo)*/ line(x-4,100,x,80); line(x,80,x+4,100); line(40,y,600,y); /*畫橫坐標(biāo)*/ line(600,y,580,y+4); line(600,y,580,y-4); j=400; settextjustify(CENTER_TEXT,CENTER_TEXT); /*縱坐標(biāo)刻度*/ for(i=0;i<=25;i+=5) line(64,j,54,j); itoa(i,buffer,10); outtextxy(30,j,buffer); j-=ystep; setcolor(1); outtextxy(x,70,

20、"Amplitude");/*縱坐標(biāo)標(biāo)注*/ j=64; settextjustify(CENTER_TEXT,TOP_TEXT); /*畫橫坐標(biāo)刻度*/ for(i=0;i<=15;i+) line(j,400,j,410); itoa(i,buffer,10); outtextxy(j,415,buffer); j+=xstep; outtextxy(600,y+10,"xn");/*橫坐標(biāo)標(biāo)注*/ /* 繪制頻譜圖 */ for(k=0;k<n;k+) setcolor(4); setlinestyle(0,0,3); line(x,y

21、,x,(y-(wk)*11); x+=step; sprintf(buffera,"%.4f",wk); if(wk>=0) outtextxy(x-step+2,(y-(wk)*12)-10,buffera); else outtextxy(x-step+2,(y-(wk)*12)+10,buffera); getch(); cleardevice(); void draw2(int n) /* 作相位圖 */ int k,driver,mode,step,x,y; float kd; int i,j,xstep=512/n,ystep=30; char buffe

22、r10,bufferpN,bufferkdN; step=512/n; driver=VGA; mode=VGAHI; x=64; y=250; initgraph(&driver,&mode,""); setbkcolor(15); /* 繪制坐標(biāo)及刻度 */ setcolor(1); moveto(200,30); outtext("The Figure of phase is as follow!"); line(x,50,x,450); /* 畫縱坐標(biāo)*/ line(x-4,70,x,50); line(x,50,x+4,70);

23、 line(40,y,600,y); /*畫橫坐標(biāo)*/ line(600,y,580,y+4); line(600,y,580,y-4); j=250; settextjustify(CENTER_TEXT,CENTER_TEXT); /*縱坐標(biāo)刻度*/ outtextxy(15,j,"0"); for(kd=0.5;kd<=3.0;kd+=0.5) line(64,j-ystep,54,j-ystep); sprintf(bufferkd,"%.4f",kd); outtextxy(30,j-ystep,bufferkd); j-=ystep;

24、j=250; for(kd=0.5;kd<=3.0;kd+=0.5) line(64,j+ystep,54,j+ystep); sprintf(bufferkd,"%.4f",-kd); outtextxy(30,j+ystep,bufferkd); j+=ystep; setcolor(1); outtextxy(x,40,"Phase");/*縱坐標(biāo)標(biāo)注*/ j=64; settextjustify(CENTER_TEXT,TOP_TEXT); /*畫橫坐標(biāo)刻度*/ for(i=0;i<=15;i+) line(j,240,j,250);

25、 itoa(i,buffer,10); outtextxy(j+10,265,buffer); j+=xstep; outtextxy(600,y+10,"xn");/*橫坐標(biāo)標(biāo)注*/ /* 繪制相位圖 */ for(k=0;k<n;k+) setcolor(4); setlinestyle(0,0,3); line(x,y,x,(y-(pk)*60); x+=step; sprintf(bufferp,"%.4f",pk); if(pk>=0) outtextxy(x-step+2,(y-(pk)*60)-10,bufferp); else

26、 outtextxy(x-step+2,(y-(pk)*60)+5,bufferp); int main() int i,method; void fft(); /*快速傅里葉變換*/ void ifft(); void initW(); /*初始化變換核*/ void change(); /*變址*/ void add(complex ,complex ,complex *); /*復(fù)數(shù)加法*/ void mul(complex ,complex ,complex *); /*復(fù)數(shù)乘法*/ void sub(complex ,complex ,complex *); /*復(fù)數(shù)減法*/ voi

27、d divi(complex ,complex ,complex *);/*復(fù)數(shù)除法*/ void output(); void output2(); /*輸出結(jié)果*/ system("cls"); PI=atan(1)*4; printf("Please input the size of x:n"); scanf("%d",&size_x); printf("Please input the data in xN:n"); for(i=0;i<size_x;i+) scanf("%lf%

28、lf",&xi.real,&xi.img); initW(); printf("Use FFT(0) or IFFT(1)?n"); scanf("%d",&method); if(method=0) fft(); else ifft(); output(); getch();/* 畫頻譜圖 */ draw(size_x); output2(); getch(); draw2(size_x); getch(); closegraph();3.2 程序運(yùn)行結(jié)果：Please input the size of X:16Please input the data in XN：xn = 0 1.2200 2.2300 2.8500 2.9800 2.6000 1.7600 0.62