快速傅里叶变换

1、实验.快速傅立叶变换一、实验目的1 .学习和掌握快速傅立叶变换（FFT）的实现过程和编程技术2. 运用FFT分析正弦信号的频谱3. 测试FFT的运算时间，比较 FFT与DFT的运算速度，获得对 FFT “快速”的感性认识。4. 锻炼和提高数字信号处理的程序设计和调试能力。二、实验原理与方法FFT并不是与DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对 DFT变换式进行一次次分解，使其成为若干小点数的组合，从而减少运算量o常用的FFT是以2为基数的，其长度 N=2%它的效率高，程序简单，使用非常方便，当要变换的序列长度不等于2的整数次方时， 为了使用以2为基数的FFT,可

2、以用末尾补零的方法，使其长度延长至 2的整数次方。本实验运用时间抽取基2 FFT,其原理、信号流图和运算过程可参阅课堂笔记、教材和其它教科书。FFT的实现要比DFT复杂，通常采用三个嵌套循环来实现。最外面的循环是分级循环，N=2M点的FFT分为M级计算。中间一层是分组循环，一级内蝶形系数Wk相同的蝶形构成一组。最内层为蝶形计算的循环，一组内不同输入数据的蝶形逐个计算。编程时要注意各级蝶形组之间的间隔、组内蝶形之间间隔、以及系数Wk变化。有不少书中有 FFT程序可供参考，但要注意理解和弄懂，不可一味照搬，要尽量自己去编。三、实验内容1. 设计说明编制时间抽取基 2 FFT程序计算前面 DFT程序

3、分析过的正弦信号，与 DFT计算的结果进行比较，以验证所编程序的正确性。其中一个用于FFT为复数运算，若用实数运算来实现，需要设置两个数组来存放输入输出数据。存放数据的实部，另一个存放数据的虚部。正弦输入信号为实数，则令其虚部为零。同样，碟形运算也要化成实数来进行，分别算出实部和虚部。当然也可以直接用复数数组和语句实现。程序设计的难点和重点在于要合理安排和正确设置碟形运算的分级循环、级内分组循环和组内分碟形循环。要注意乘法系数时的变化以及各循环变量的变化和调整。程序中，在进行 FFT计算之前，先要将输入数据排列成二进制倒序的形式，这可通过将有关数 据相互换位来实现。正弦抽样信号的产生和输出频谱

4、图的绘制与前面的DFT程序相同。2. 实验步骤(1) 将编制的程序输入计算机，调试、运行正确。即输入实验一中的12组数据，看FFT输出谱线的结果是否与前面DFT的结果一致。(2) 在不同的FFT长度下(可分别令 N=64, 128, 256, 512, 1024)，运行 FFT程序，观察、比较和记录运行时间。由于现在计算机的运算速度很快，计算一次FFT的时间很短(在ms数量级)，不便观察和记录。可在程序中添置循环语句，让其重复计算FFT许多次(可重复 100次、1000次或更多，视情况而定)，使总的运算时间达到数秒或数十秒钟，从而可以通过人工来观察和记录。此时注意，为了防止 FFT多次循环计算

5、，使输入输出数据越来越大而发生溢出，可令输入数据全部为零。这样，FFT的运算量保持不变，但输入输出始终为零，不会因多次循环而发生溢出。在观测FFT的运算时间时，作图部分的语句不要运行(一方面，它不应该计入FFT或DFT的运算时间。另一方面，当N较大时，原有程序中的作图语句会出错)，可将其注释掉。测量FFT运行时间的一个更好、更精确的方法是在程序中调用读取计算机时钟的函数。TURBOC中提供多种读取计算机时间的方法和函数，一种比较简单的方法是用biostime() 函数,后面给出了相应的参考程序。(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的级数 M */for(j=1,i=1;i<n;i+)m=i;j=2*j;if(j=n) break;/*改变输入数据的顺序*/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;)/*三重嵌套循环进行碟形计算*/n1=1;for(l=1;l<=m;l+)/*分级循环 */(n1=2*n1;n2=n1/2;e=3n2;for(k=0;k<n2;k+) /*分组循环 */(c=cos(k*e); /*计算 W */s=-sin(k*e);(for(i=k;i<n;i+=n1) /*碟形循环*/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);/*产生正弦输入数据*/for (n=0;n<N;n+)xn=sin(2*PI*f*n*T);yn=0;/*调用子程序计算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运算时间的有关程序语句，应将它们插入前面main()函数内适当的位置。unsigned int start,end;float time;int i; /*提示：这两行放在主程序开头数据类型说明部分/*提示：以下部分替代主程序中的fft(x,y,N)start=biostime(0,0);

13、/*for (n=0;n<2000;n+) /*读取开始时间*/循环计算FFT 2000次*/fft(x,y,N);for (i=0;i<N;n+)/*清除x,y的值，防止下一次的FFT数据溢出(xi=0;yi=0;)end=biostime(0,0); /*读取结束时间*/time=(end-start)/2000.0*55; /* biostime()读取的一个时间单位约为55ms */printf("It took %f msn", time); /*屏幕输出FFT的运算时间*/* draw(N); */ /*不运行作图语句*/(3)测量DFT运行时间的程序DFT程序中以下给出的是用于测量DFT运算时间的有关程序语句，应将它们插入实验main()函数内适当的位置。unsigned int start,end;float time;int j;start=biostime(0,0); /*for(j=0;j<50;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); /*读取结束时间*/t