1、- 0 -FFT 程序的设计一、设计内容1.1 设计要求用 C 程序平台设计快速傅里叶变换的程序，既要有设计的理论内容，也要有相应的流程图，及源代码。1.2 设计目的1.掌握了解快速傅里叶变换程序的设计。2.学习 C 语言平台的使用方法及基本功能。3.加深对快速傅里叶变换的理解，熟悉快速傅里叶变换子程序调用。4.熟悉应用快速傅里叶变换实现两个序列线性卷积的方法二、设计思路2.1 设计原理2.1.1 FFT 算法原理FFT 并不是另一种变换傅里叶变换，而是 DFT 的一种快速算法。它通过对 DFT 变换式的一次次分解，将长序列的 DFT 分解为一系列短序列 DFT 的组合，从而减少运算量。常用的
2、 FFT 是以 2 为基数的，其长度 N=2L。当要变换的序列长度不等于 2 的整数次方时，是为了使用以 2 为基数的 FFT，可以用末位补零的方法，使其长度延长至 2 的整数次方。- 1 -2.1.2 FFT 的运算量对于一个 N 长的序列，直接计算 DFT 需要 N2次复数乘及 N(N-1)=N2复数加。按时间抽选法 FFT，当 N=2L时，共有 L 级蝶形，每级都由 N/2 个蝶形运算组成，每个蝶形有一次复数乘、二次复数加，总的运算量为（N/2）L 次复数乘，NL 次复数加。2.1.3 用 FFT 计算线性卷积用 FFT 可以实现两个序列的圆周卷积。在一定的条件下，可以使圆周卷积等于线性
3、卷积。一般情况，设两个序列的长度分别为 N1和 N2，要使圆周卷积等于线性卷积的充分必要条件是 FFT 的长度 N N 1+N2-1 对于长度不足 N 的两个序列，分别将他们补零延长到N。2.2 C 程序设计法及程序流程图FFT 算法是离散傅立叶变换的快速算法，其结果应该与离散傅立叶变换所算出的结果一致。图 1 是采样点数为 8 的 FFT 算法的一种形式，本程序也是以这种形式设计的。x(0) X(0)x(4) X(1)x(2) X(2)x(6) X(3)x(1) X(4)x(5) X(5)x(3) X(6)x(7) X(7)图 1 FFT 算法形式程序设计的基本思路是在程序开始时刻要求输入采
4、样点数，如果采样点数不是 2 的整数次方（不包括 0 次方） ，则要求输入采样点的数值，并根据采样点数分配响应的数组x0(1)-1x0(2)x0(3)x0(4)x0(5)x0(6)x0(7)-1-1-1x1(1)x1(2)x1(3)x1(4)x1(5)x1(6)x1(7)08W208W2-1-1-1-1x2(1)x2(2)x2(3)x2(4)x2(5)x2(6)x2(7)08W183-1-1-1-1x3(1)x3(2)x3(3)x3(4)x3(5)x3(6)x3(7)x0(0) m=0 x1(0) m=1 x2(0) m=2 x3(0)- 2 -大小，计算迭代次数。然后对采样点进行逆二进制排序
5、，再按上图所示的算法进行计算，程序流程图如图 2：图 2 流程图本程序在 VC+6.0 的开发环境下创建 Win32 Console Application 工程对程序进行设计，下面分别对程序设计中复数类的应用，判断和求迭代次数，逆二进制排序，蝶形运算进行具体说明。2.3 快速傅里叶变换程序设计各类的应用2.3.1 复数类的应用C+语言本身并不包含复数数据类型，但在经过的不断扩展后，如今的STL（Standard Template Library）中已经包含了相应的复数类，但在尝试后发现其功能仍不能满足本程序设计的要求，所以仍需要在其基础上编写函数。程序开始输入采样点数采样点数是否为 2 的不
6、为 0 的整数次方？根据输入点数分配数组大小按迭代，分组，碟形运算的顺序进行运算输出计算结果程序结束NY- 3 -在 FFT 程序设计中，复数类主要被用来计算两复数的加法和乘法以及旋转因子 ，其中KNW，式中 N=2 的 i 次方，i 代表第 i 次迭代，而 k 代表第 k 次蝶形运算，由于lNjW2STL 中的 Complex 类不能进行带参数的复数运算，所以本程序编写了带参数的复数运算cW（int cm,int sign2m） ，用于计算 ，设计的基本思路是化复数乘法和除法运算ikj2为几个复数基本单元的加法运算和除法运算，其中加法运算的次数由函数输入参量 int cm 决定，复数除法运算
7、次数由 int sign2m 决定。基本加法单元为复数类对象 wmt(0,-1*2*pi)，基本除法类对象为 pm(2,0)。下面是 cW（int cm,int sign2m）的代码：complex cW(int cm,int sign2m)int s;complex Wm(0,0),pm(2,0),Wmt(0,-1*2*pi);for(s=0;stemp;inputa=temp;for(a=0;a#include #include using namespace std;#define pi 4*atan(1)int iternumb(int numb)int iternumb=0;if(n
8、umb=0)|(numb=1)cout cW(int cm,int sign2m)int s;complex Wm(0,0),pm(2,0),Wmt(0,-1*2*pi);for(s=0;s tempA1,tempA2,W;coutnum;iternum=iternumb(num);double* input=(double*)malloc(sizeof(double)*num);complex* A=(complex*)malloc(sizeof(complex)*num);int* inverse=(int*)malloc(sizeof(int)*iternum);for(a=0;atem
9、p;inputa=temp;for(a=0;aint fft(complex * a,int l)const double pai=3.141592653589793complex u,w,t;unsigned n=1,nv2,nm1,k,le,lei,ip;unsigned i,j,m;double tmp;n1;nm1=n-1;- 12 -i=0;for(i=0;i,nm1,i+)if(i,j)t=ajaj=ai;ai=t; k=nv2;while(k=1 j+=kle=1for(m=1:m=1:m+)lei=lele=1;- 13 -u=complex(1,0);tmp=pai/lei;
10、w=complex(cos(tmp),-sin(tmp);for(j=0;jlei;j+)for(i=j;in;i+=le)ip=i+lei;t=aip*u;aip=ai-t;ai+=t;u*=w;return 0;2.5 程序运行结果输入采样点数为 4，数据分别为 1,3,5,6.- 14 -图 3 C 程序结果由图 3 知所得结果为 X0=15，X1=-4+3i，X2=-3，X3=-4-3i在 Matlab 中进行检验：图 4 Matlab 验证结果由图 4 结果与所设计程序一致。输入采样点数为 8 分别为 2,3,4,5,7,9,10,11. 详见图 5：- 15 -图 5 C 程序结果
11、在 Matlab 中检验，如图 6：图 6 Matlab 验证结果与所设计程序仍然一致。结论：通过对程序的设计以及对运行结果的验证可以说明所设计的程序运行结果正确，完成了对采样点进行 FFT 变换的要求。（注明：由于程序二中 complex 头文件不是微软公司的，因此无法用 VC 运行。 ）- 16 -三、心得体会学习数字信号处理已经有一个学期了，现在和大家聊一下我学习过程中的一些经历和想法吧。刚接触时，和很多同学都一样，觉得很难，很深奥，不易懂。可是其中却蕴含着许多技巧，渐渐地，对其产生了兴趣从此不能自拔。很多个日夜就这样陪伴着它度过了。期间也遇到过非常多的问题，也一度被这些问题所困惑等到回
12、过头来，看到自己曾经走过的路，唏嘘不已。经常混迹于论坛里，也看到了很多学成者发的帖子，对自己的帮助很大，也交到了一些关于这方面有所长的朋友，经常上网交流。路学习过来的过程中，帮助最大之一无疑来自于网络了。很多时候，通过网络，我们都可以获取到所需要的学习资料。但是，随着我们学习的深入，我们会慢慢发现，网络提供的东西是有限度的，好像大部分的资料都差不多，或者说是适合大部分的初学者所需，而当我们想更进一步提高时，却发现能够获取到的资料越来越少，这也需要我们的帮助之二老师，通过老师的帮助，我们的知识面得到了扩展与提升。我们这里所做的是有关快速傅里叶变换的，它只是数字信号处理中的一小部分，还有很多诸如此类的知识都是令我们受益匪浅。之前没想到 C 语言居然可以跟数字信号处理联系在一起，现在都得以发挥了它们的功效，完美的诠释了我们所学的内容，这也给我们启示，在以后的学习道路中，很多东西都可以联系起来，很多看似平凡的东西都有它独到的用途，使它发光发亮。四、参考文献【1】丁玉美. 数字信号处理. 电子科技大学出版社【2】周辉. 数字信号处理及其 MATLAB 实现. 中国林业出版社【3】程佩青. 数字信号处理教程. 清华大学出版社 【4】余成波. 数字信号处理及 MATLAB 实现. 清华大学出版社【5】张明照. 应用 MATLAB 实现信号分析和处理. 科技出版社
