快速傅里叶变换

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

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

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

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

5、可以通过人工来观察和记录。此时注意，为了防止FFT多次循环计算，使输入输出数据越来越大而发生溢出，可令输入数据全部为零。这样，FFT的运算量保持不变，但输入输出始终为零，不会因多次循环而发生溢出。在观测FFT的运算时间时，作图部分的语句不要运行（一方面，它不应该计入FFT或DFT的运算时间。另一方面，当N较大时，原有程序中的作图语句会出错），可将其注释掉。      测量FFT运行时间的一个更好、更精确的方法是在程序中调用读取计算机时钟的函数。TURBO C中提供多种读取计算机时间的方法和函数，一种比较简单的方法是用biostime()函数，后

6、面给出了相应的参考程序。（3）   用上述类似的方法，在不同的长度下（N=64，128，256，512，1024），运行实验一中的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

7、;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;

8、      xj=xi;      yj=yi;      xi=tr;      yi=ti;        k=n/2;   while(k<(j+1)             j=j-k; 

9、0;    k=k/2;         j=j+k;    /* 三重嵌套循环进行碟形计算 */   n1=1;   for(l=1;l<=m;l+) /* 分级循环 */         n1=2*n1;     n2=n1/2;     for(k=0;k<n2

10、;k+)  /* 分组循环 */             c=cos(k*e);  /* 计算WNk */       s=-sin(k*e);                     for(i=k;i<n;i+

11、=n1) /* 碟形循环 */         j=i+n2;     tr=c*xj-s*yj;  /* 计算碟形 */     ti=c*yj+s*xj;     xj=xi-tr;     yj=yi-ti;     xi=xi+tr;     yi=yi+ti;      

12、0;                   /* 主程序部分 */  #include<stdio.h>  #include<math.h>  #include<conio.h>  #include<graphics.h>  #define PI 3.1415926  float x1024,y1024,w1024;  void draw(i

13、nt);  void axis(int,int);   main()         int N,f,n,k;    float T; /* 键盘输入f、N、T */printf("the frequency of the sine wave f=");    scanf("%d",&f);    printf("the number of samp

14、les N=");    scanf("%d",&N);    printf("the time step of samples T=");    scanf("%f",&T); /* 产生正弦输入数据 */    for (n=0;n<N;n+)           

15、0; xn=sin(2*PI*f*n*T);       yn=0;      /* 调用子程序计算FFT */    fft(x,y,N);/* 计算FFT的幅度 */    for(k=0;k<N;k+)             wk=sqrt(xk*xk+yk*yk);   &

16、#160;  /* 画FFT的幅度谱 */          draw(N);      getch();      closegraph();   /* 作图子程序 */  void draw(int N)         int k,driver,mode,step,x,y;   &#

17、160; driver=VGA;     mode=VGAHI;     step=512/N;     x=64;y=400;     registerbgidriver(driver);     initgraph(&driver,&mode,"c:tc");     axis(x,y);  

18、60;  outtextxy(x+step+546),y+8,"k");     outtextxy(x+step-8),y-290,"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

19、 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;        &

20、#160;   void axis(int x,int y)              line(10,y,630,y);        line(610,(y-4),630,y);        line(610,(y+4),630,y);        l

21、ine(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

22、start,end;float time;int i;    /*提示：这两行放在主程序开头数据类型说明部分        /*提示：以下部分替代主程序中的fft(x,y,N)start=biostime(0,0);       /* 读取开始时间 */    for (n=0;n<2000;n+)    /* 循环计算FFT 2000次 */   

23、  fft(x,y,N);for (i=0;i<N;n+)      /*清除x,y的值，防止下一次的FFT数据溢出             xi=0;       yi=0;              end=biostime(0,0); &

24、#160;   /* 读取结束时间 */    time=(end-start)/2000.0*55;    /* biostime()读取的一个时间单位约为55ms */    printf("It took %f msn", time);    /* 屏幕输出FFT的运算时间 */      /*      draw(N); */&#

25、160;   /* 不运行作图语句 */（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;      &