数字信号处理课程设计DFT的快速算法快速傅里叶变换FFT.doc

1、燕 山 大 学 课 程 设 计 说 明 书目录前言第一章 离散傅里叶变换DFT31.1 DFT定义31.2 DFT的快速算法快速傅里叶变换(FFT) 3第二章 基2 DIT-FFT算法42.1 按时域抽取的基2 DIT-FFT算法4第三章 基于C语言设计16点基2DIT-FFT程序及运行结果 63.1 按时间抽取的基2DIT-FFT程序 63.2 程序运行结果 15第四章 课程设计的总结 17参考文献 17 前 言信号(signal)是一种物理体现，或是传递信息的函数。而信息是信号的具体内容。模拟信号(analog signal)：指时间连续、幅度连续的信号。数字信号(digital sign

2、al)：时间和幅度上都是离散(量化)的信号。 数字信号可用一序列的数表示，而每个数又可表示为二制码的形式，适合计算机处理。数字信号处理是将信号以数字方式表示并处理的理论和技术。数字信号处理与模拟信号处理是信号处理的子集。 其目的是对真实世界的连续模拟信号进行测量或滤波。因此在进行数字信号处理之前需要将信号从模拟域转换到数字域，这通常通过模数转换器实现。而数字信号处理的输出经常也要变换到模拟域，这是通过数模转换器实现的。数字信号处理的核心算法是离散傅立叶变换(DFT)，是DFT使信号在数字域和频域都实现了离散化，从而可以用通用计算机处理离散信号。而使数字信号处理从理论走向实用的是快速傅立叶变换(

3、FFT)，FFT的出现大大减少了DFT的运算量，使实时的数字信号处理成为可能、极大促进了该学科的发展。随着大规模集成电路以及数字计算机的飞速发展，加之从60年代末以来数字信号处理理论和技术的成熟和完善，用数字方法来处理信号，即数字信号处理，已逐渐取代模拟信号处理。 随着信息时代、数字世界的到来，数字信号处理已成为一门极其重要的学科和技术领域。第一章 离散傅里叶变换DFT离散傅立叶变换（DFT）实现了信号首次在频域表示的离散化，使得频域也能够用计算机进行处理。并且这种DFT变换可以有多种实用的快速算法。使信号处理在时、频域的处理和转换均可离散化和快速化，因而具有重要的理论意义和应用价值。1.1

4、DFT定义设序列x(n)长度为M，定义x(n)的N点DFT为式中，N称为离散傅里叶变换区间长度，要求N M。为书写简单，令 ，因此通常将N点DFT表示为1.2 DFT的快速算法快速傅里叶变换(FFT)DFT使计算机在频域处理信号成为可能，但是当N很大时，直接计算N点DFT的计算量非常大。快速傅里叶变换（FFT，Fast Fourier Transform）可使实现DFT的运算量下降几个数量级，从而使数字信号处理的速度大大提高，工程应用成为可能。人们已经研究出多种FFT算法，它们的复杂度和运算效率各不相同。快速傅里叶变换就是不断地将长序列的DFT分解为短序列的DFT，并利用 的周期性和对称性及其

5、一些特殊值来减少DFT运算量的快速算法。FFT算法分类：1.时间抽选法 DIT：Decimation-In-Time2.频率抽选法 DIF：Decimation-In-Frequency时间域抽取：基2时间抽取(Decimation in time)DIT-FFT算法频率域抽取：基2频率抽取(Decimation in frequency)DIF-FFT算法第二章 基2 DIT-FFT算法2.1 按时间抽取的基2 DIT-FFT算法：1、按时间抽取的基2 DIT-FFT算法原理先设序列点数为N=2M，M为整数。如果不满足这个条件，可以人为地加上若干零值点，使之达到这一要求。这种N为2的整数幂的

6、FFT称基-2 FFT。设输入序列长度为N=2M (M为正整数) ，将该序列按时间顺序的奇偶分解为越来越短的子序列，称为按时间抽取(DIT )的FFT算法。 序列x(n)的N点DFT为 将上面的和式按n的奇偶性分解为 令x1(l)=x(2l),x2(l)=x(2l+1)，因为 ，所以上式可写成上式说明，按n的奇偶性将x(n)分解为两个N/2长的序列x1(l)和x2(l)，则N点DFT可分解为两个N/2点DFT来计算。用X1(k)和X2(k)分别表示x1(l)和x2(l)的N/2点DFT，即有上述公式，及X1（k）、X2(k)的隐含周期性得到：这样，就将N点DFT的计算分解为计算两个N/2点离散

7、傅里叶变换X1(k)和X2(k)，再计算上式。蝶形图：2、 按时间抽选的FFT算法的特点：（1）原位运算由图4.2.4可以看出，DIT-FFT的运算过程很有规律。N=2M点的FFT共进行M级运算，每级由N/2个蝶形运算组成。同一级中，每个蝶形的两个输入数据只对计算本蝶形有用，而且每个蝶形的输入、输出数据结点又同在一条水平线上，这就意味着计算完一个蝶形后，所得输出数据可立即存入原输入数据所占用的存储单元。这样，经过M级运算后，原来存放输入序列数据的N个存储单元中便依次存放 的N个值。这种利用同一存储单元存储蝶形计算输入、输出数据的方法称为原位计算。原位计算可节省大量内存，从而使设备成本降低。（2

8、倒序规律按原位计算时，FFT的输出 是按正常顺序排列在存储单元中，但是这时输入x(n)却不是、按自然顺序存储的，而是按x(0)，x(4)，x(7)的顺序存入存储单元，看起来好象是“混乱无序”的，实际上是有规律的，我们称之为倒序。第三章 基于C语言设计16点基2DFT-FFT程序及运行结果设计参数及要求：本次试验要求利用C语言设计16点基2DIT-FFT程序，并对如下程序进行分析，给出输出结果并画出信号的幅值图和相位谱。n012345678910X(n)01.222.232.85 2.982.601.760.62-0.62-1.76-2.601112131415-2.98-2.85-2.23-

9、1.2203.1 按时间抽取的基2DFT-FFT程序：/*时间抽选基2 FFT算法C语言实现*/#include #include #include #include #define N 1000#include/*定义复数类型*/typedef struct double real; double img;complex;complex xN, *W; /*输入序列,变换核*/int size_x=0; /*输入序列的大小，在本程序中仅限2的次幂*/double PI; /*圆周率*/int z;float wN,pN;char strN;/*初始化变换核*/void initW() int

10、 i; W=(complex *)malloc(sizeof(complex) * size_x); for(i=0;ireal=a.real+b.real; c-img=a.img+b.img; /*乘法*/void mul(complex 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.

11、img; /*除法*/void divi(complex a,complex b,complex *c) c-real=( a.real*b.real+a.img*b.img )/( b.real*b.real+b.img*b.img); c-img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img); /*变址计算，将x(n)码位倒置*/void change() complex temp; unsigned short i=0,j=0,k=0; double t; for(i=0;i0 ) j=j1; if(ji) temp=x

12、i; 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+) /*一级蝶形运算*/ l=1i; for(j=0;jsize_x;j+= 2*l ) /*一组蝶形运算*/ for(k=0;kl;k+) /*一个蝶形运算*/ mul(xj+k+l,Wsize_x*k/2/l,&product); add(xj+k,product,&up); sub(xj+k,product,&down); xj+k=

13、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+) /*一级蝶形运算*/ l/=2; for(j=0;jsize_x;j+= 2*l ) /*一组蝶形运算*/ for(k=0;kl;k+) /*一个蝶形运算*/ add(xj+k,xj+k+l,&up); up.real/=2;up.img/=2; sub(xj+k,xj+k+l,&down); down.real/=2;down.img/=2;

14、divi(down,Wsize_x*k/2/l,&down); xj+k=up; xj+k+l=down; change();/*输出傅里叶变换的结果和幅值结果*/void output() int i; clrscr() ; printf(The result are as followsnn); for(i=0;i=0.0001)printf(+%.4fjn,xi.img); else if(fabs(xi.img)0.0001)printf(n); else printf(%.4fjn,xi.img); for( z=0;z16;z+) wz=sqrt(xz.real)*(xz.real

15、)+(xz.img)*(xz.img); printf(nnThe Amplitudes are as followsnn); for(z=0;zsize_x;z+) printf(%.4fn,wz); /*输出傅里叶变换相位的结果*/void output2() int i; setbkcolor(0); for(i=0;isize_x;i+) for( z=0;z16;z+) pz=atan(xz.img)/(xz.real); printf(nnnThe Phase are as followsnn); for(z=0;zsize_x;z+) setcolor(1); printf(%.

16、4fn,pz); void draw(int n) /* 作幅值图 */ int k,driver,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); /* 绘制坐标及刻度 */ setcolor(1); moveto(200,30); outtext(The Figure of Amplitude is as follow!)

17、 line(x,80,x,400); /* 画纵坐标*/ line(x-4,100,x,80); line(x,80,x+4,100); line(40,y,600,y); /*画横坐标*/ line(600,y,580,y+4); line(600,y,580,y-4); j=400; settextjustify(CENTER_TEXT,CENTER_TEXT); /*纵坐标刻度*/ for(i=0;i=25;i+=5) line(64,j,54,j); itoa(i,buffer,10); outtextxy(30,j,buffer); j-=ystep; setcolor(1); o

18、uttextxy(x,70,Amplitude);/*纵坐标标注*/ j=64; settextjustify(CENTER_TEXT,TOP_TEXT); /*画横坐标刻度*/ 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);/*横坐标标注*/ /* 绘制频谱图 */ for(k=0;k=0) outtextxy(x-step+2,(y-(wk)*12)-10,buffera); else outtextxy(x-s

19、tep+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 buffer10,bufferpN,bufferkdN; step=512/n; driver=VGA; mode=VGAHI; x=64; y=250; initgraph(&driver,&mode,); setbkcolor(15); /* 绘制坐标及刻度 */ setcolor

20、1); moveto(200,30); outtext(The Figure of phase is as follow!); line(x,50,x,450); /* 画纵坐标*/ line(x-4,70,x,50); line(x,50,x+4,70); line(40,y,600,y); /*画横坐标*/ line(600,y,580,y+4); line(600,y,580,y-4); j=250; settextjustify(CENTER_TEXT,CENTER_TEXT); /*纵坐标刻度*/ outtextxy(15,j,0); for(kd=0.5;kd=3.0;kd+=0

21、5) line(64,j-ystep,54,j-ystep); sprintf(bufferkd,%.4f,kd); outtextxy(30,j-ystep,bufferkd); j-=ystep; 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);/*纵坐标标注*/ j=64; settextjustify

22、CENTER_TEXT,TOP_TEXT); /*画横坐标刻度*/ for(i=0;i=15;i+) line(j,240,j,250); itoa(i,buffer,10); outtextxy(j+10,265,buffer); j+=xstep; outtextxy(600,y+10,xn);/*横坐标标注*/ /* 绘制相位图 */ for(k=0;k=0) outtextxy(x-step+2,(y-(pk)*60)-10,bufferp); else outtextxy(x-step+2,(y-(pk)*60)+5,bufferp); int main() int i,metho

23、d; void fft(); /*快速傅里叶变换*/ void ifft(); void initW(); /*初始化变换核*/ void change(); /*变址*/ void add(complex ,complex ,complex *); /*复数加法*/ void mul(complex ,complex ,complex *); /*复数乘法*/ void sub(complex ,complex ,complex *); /*复数减法*/ void divi(complex ,complex ,complex *);/*复数除法*/ void output(); void o

24、utput2(); /*输出结果*/ 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;isize_x;i+) scanf(%lf%lf,&xi.real,&xi.img); initW(); printf(Use FFT(0) or IFFT(1)?n); scanf(%d,&method); if(method=0) fft(); else ifft(); output(); ge

25、tch();/* 画频谱图 */ draw(size_x); output2(); getch(); draw2(size_x); getch(); closegraph();3.2 程序运行结果：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.6200 -0.6200 -1.7600 -2.6000 -2.9800 -2.8500 -2.2300 -1.2200 004.49118667428684 - 22.57872013

26、41516i-0.863137084989848 + 2.08379725676967i-0.708894177324789 + 1.06093511090429i-0.660000000000000 + 0.660000000000000i-0.667400913065437 + 0.445943033038161i-0.636862915010152 + 0.263797256769669i-0.634891583896613 + 0.126287787982296i-0.640000000000000-0.634891583896613 - 0.126287787982296i-0.63

27、6862915010152 - 0.263797256769669i-0.667400913065437 - 0.445943033038161i-0.660000000000000 - 0.660000000000000i-0.708894177324789 - 1.06093511090429i-0.863137084989848 - 2.08379725676967i4.49118667428684 + 22.5787201341516i序列X(n):傅里叶变换后：第四章 课程设计的总结通过课程设计，熟悉了DFT的性质。加深了对信号频谱的概念及性质的理解。高在做课程设计前，预习课本DFT及FFT及基2DIT-FFT的有关概念，加深了理解，用C语言编程来实现，进一步熟悉了 C语言。 经过这次课程设计，我收获颇多，我不但对2DIT-FFT知识有了更深入的认识也掌握了实际应用的技巧，对其有了更加感性的认识，而且自己运用MATLAB编程能力也有明显的提高。参考文献：1 谢平，王娜，林洪彬 信号处理原理及应用 机械工业出版社 2 MATLAB在数字信号处理中的应用 薛年喜 编著 清华大学出版社3 MATLAB 程序设计及其在信号处理中的应用 聂祥飞 王海宝 谭泽富 主编 西安交通大学出版社 共 17 页 第17页