快速傅里叶变换(FFT)的计算机实现信号与系统课程设计Word

华中科技大学信号与系统课程设论文快速傅里叶变换（FFT）的计算机实现学院：班级：学号：姓名：指导老师：二〇一三年八月摘要用C语言编程完成对输入波形的时域采样的FFT变换以及频域分析，同时用DFT 变换来验证FFT变换结果的正确性。

时域信号的输入有两种方式：一种是依次输入时域信号（仅支持多个正弦及余弦信号之和的形式）各谐波分量的幅值和角频率值，另一种是直接输入时域信号的采样值。

然后进行DFT变换和FFT变换，两者结果应该是一样的。

DFT变换的实现直接脱胎于定义，FFT变换采用的是基2时间抽取算法，用倒位序算法和蝶形算法实现。

关键词：傅里叶变换；DFT； FFT；倒位序1背景知识1.1 为什么要进行傅里叶变换在进行科学研究的过程中，很重要的一点就是为一个系列的问题找到一个通法，从而为实际的应用打下基础。

在信号分析这个方向，如果直接对时域信号进行分析，那么我们将发现，很难找到一种固定的方法来进行分析，这是因为信号对时间t的函数形式存在无数种，我们是无法将这些函数以及这些函数的各种形式的组合都统一起来进行研究的。

而进行傅里叶变换之后，我们就能很好达到这个目的——用一种方法来研究所有的信号（这些信号也需要满足一定的条件，但范围是非常广的）。

那么为什么傅里叶变换可以达到这样的目的呢？对于一个时域信号，我们习惯从时间的角度进行理解，也就是以时间为自变量，以信号值为因变量，一个信号是该信号在所有的时间点的值的叠加，每个信号分量在时间域上只占据了一个点，要想得到这个信号的所有信息，我们需要知道这个信号在时间轴上每一点的值，缺一不可。

傅里叶变换之后，依然是一个叠加的问题，只不过由时间角度的叠加变成了频率角度的叠加，这时每一个信号分量都覆盖了一个时间域上的整个区间，并且每一个信号分量都是周期性的（三角函数的周期性），这样，只需要知道每个信号分量的幅值、频率、相位就能在时间轴上得到它的所有信息，而所有信号分量的叠加就得到了原来的信号。

并且我们并非需要将所有信号分量都叠加起来，这是因为傅里叶变换之后，随着信号分量频率值的上升，信号分量幅值呈一个较快的下降趋势，在精度允许的范围内，我们并不需要去考虑频率值超出某个范围的那部分信号分量，我们所考虑的那个频率区间的信号分量的叠加已经能够很好的还原出原信号，而这个频率区间只占据了频率轴很少的部分，从而简化了后面的分析过程。

同时，若原信号是周期性的，那该信号在频率轴上将只占据有限个点，分析难度更是大大减小。

1.2傅里叶级数1.2.1频率分量与频率成分对于时域信号来说，频率含量就是信号被分成的正弦波簇所确定的所有频率分量。

例如，有（1-1）式中，N rad/s正弦波的相位。

根据式（1-1）其1-1）定义的信号完全由频率，振幅和相信的相位所确定。

由式（1-1）定义的信号特征，可以通过对组成该信号的正弦频率，振幅，相位来研究。

1.2.2周期信号的三角傅里叶级数式：在式（1-2）rad/s）在式（1-21-6）确定：表达式（1-2k次谐波项为的重要特性。

由式（1-2）给定的三角傅里叶级数可以在写成余弦函数和相位的形式：式中：并且：1.2.3复指数级数由式（1-2）给定的三角傅里叶级数可以表示成如下的指数形式：在式（1-10rad/s），且式（1-10应该注意到，式（1-11）给定的也能在任何整周期上的积分计算出来，例如：1.3傅里叶变换及其直角和极坐标表达式一般意义上，（1-13）式中的积分收敛（存在），是“表现良好”的和绝对可积的，那么积分就是收敛的。

绝对可积是指下式成立：“表现良好”就是信号在任何有限的时间段内存在有限的间断点和最大、最小值。

除了脉冲信号，我们感兴趣的大部分信号都是“表现良好”的。

所有实际信号（即物理上可以产生的信号）都是“表现良好”的，且满足（1-14）式。

数。

因为复数既可以表示成直角坐标形式，也可以表示成极坐标形式，所以，也可以表示成这两种形式。

下面将定义这两种形式。

则的直角坐标表达式为：式中，是的模，是的辐角。

利用下面的关系，可以由直角坐标表达式换成极坐标表达式：1.4 离散时间信号的傅里叶分析1.4.1离散时间傅里叶变换DTFT他的离散时间傅里叶变换（DTFT ）定义为：一般地，由上式定义的般意义的DTFTDTFT如果是时限的离散时间信号，显然，上式中的和是有限值，类似这样的信号都存在一般意义下的DTFT 。

-DTFT坐标其中：:其中：1.4.2离散傅里叶变换DTFT实现DTFT ，必须在频率上离散化，从而引出了下面定义的离散傅里叶变换的概念。

DFT义为：由式（1-16）可见，k角坐标形式。

极坐标形式是：直角坐标形式是：其中：1.5 FFT 算法时间抽取算法的基本思想是把时间间隔细分，计算可以分成两部分，首先对符号进行简化，1的N 次开放，即：假设N>1N 点DFT 和DFT 反变换为：现在令N 是一个偶整数，以便N/2是一个整数。

已知信号x[n],令a[n]为x[n]的偶数次项，b[n]为x[n]N/2）点DFT ，即：N 点DFT ，可有：次乘法。

为了看清这一点，首先注意到计算kN/2次乘法，所以，总N 非常大时，可以大大减少乘法次数。

如果N/2q为正整数，这个递减的过程可以重复进行，知道信号只剩下一个非零值。

下图给出了N=8时FFT算法的流程图（图1-1）n的范围是0~N-1，时间标号n可以用q位二进制数来表示，把这q位二进制数进行倒位序排列，记得FFT1-1给出了图1-1中FFT算法的输入信号值的顺序。

图1-1 N=8时的FFT算法流程图表1-1 N=8时的倒位序排列时间（n）二进制代码倒位序代码排序0000000x[0]1001100x[4]2010010x[2]3011110x[6]4100001x[1]5101101x[5]6110011x[3]7111111x[7]1.5.1倒位序算法在进行FFT运算时，第一步要实现的就是倒位序排列。

假设使用A[I]存的是顺序位序，而B[J]存的是倒位序。

I<J的时候需要变序，I>J的时候就不用。

不然就相当于作了两次变序，又变回去了。

从表1-1可知，按自然顺序排列的二进制数，其下面的数总是比其上面的数大1，即下面的数是上面的数在最低位加1并向高位进位而得到的。

而倒位序二进制数的下面的数是上面的数在最高位加1并由高位向低位进位而得到。

由数组的性质可知，I、J都是从0开始，若已知某个倒位序J，要求下一个倒位序数，则应先判断J的最高位是否为0，这可与flag=N/2相比较。

如果flag>J，则J的最高位为0，只要把该位变为1（J与flag=N/2相加即可），就得到下一个倒位序数；如果flag<=J，则J的最高位为1，可将最高位变为0（J与flag=N/2相减即可）。

然后还需要判断次高位，这可与flag’=N\4相比较，若次高位为0，则需将它变为1（加N\4即可）其他位不变，既得到下一个倒位序数；若次高位是1，则需将它也变为0。

依此类推即可得到最后的倒位序排列。

2c语言实现见附录3利用c程序进行频谱分析3.1直接输入离散信号取样值进行DFT和FFT3.1.1取样点N=4，时域信号取样值xn[]={1，2，2，1}结果见下图可见FFT和DFT的结果基本上是一样的，仅在小数点第六位才有一点差别，见XK[3]。

3.1.2取样点N=5，时域信号取样值xn[]={1，2，2，4，5}结果见下图当N=时，DFT仍然可以进行，因为DFT的代码“翻译”自DFT的定义而来，而FFT的代码是“翻译”自倒位序算法和蝶形算法，这两种算法对取样点数有要求，N必须是以2为底的正指数。

同时，此c程序还对取样点最大值有要求，不得超过128.3.2输入正弦谐波分量信息，让计算机进行取样3.2.1N=4、8、16此时，，最大谐波次数maxharmanicorder=5，谐波分量系数为a[1]=9,a[3]=3,a[5]=1。

N=4时结果见下图：N=8时结果见下图：N=16时，结果见下图（仅出示FFT）幅度谱图（从左至右N=4、6、8）：幅度谱分析：由上图可知，随着取样点数的增加，边界点的幅值急剧上升，不过图像的整体趋势不变相位谱（从左至右N=4、6、8）：由上图可知，随着取样点数的增加，相位谱的变化趋势没有发生变化，只是将原本的取样点变得更密集了而已。

3.2.2取样点数N=8，时域信号函数filemaxharmanicorder=5，a[1，3，5]={9，3，1}，此时，maxharmanicorder=7，a[1，3，5，7]={9，3，1，0.5}，此时，maxharmanicorder=9，a[1，3，5，7，9]={9，3，1，0.5，0.25}幅度谱（从左至右高次谐波分量有增加）：由上图可知，随着谐波分量的增加，幅度谱的个别取样点幅值有变化，由于谐波分量增量很小，因此幅值变化也很小，同时，整体波形变化也不大。

由上图分析可知，随着谐波分量的增加，相位谱完全没有变化。

这是因为增加的谐波分量的频率值时基波的整数倍。

4总结在此次的课程设计中，我是用C语言和MS Visual Studio 2010完成C语言编程部分，Matlab 完成绘图部分，Word完成文字部分。

完成这份课设让我收获良多。

从课程的角度来说，加深了对傅里叶变换和快速傅里叶变换的理解，学会了从频域的角度来理解信号，让我体会到，换一个角度来解决相同的问题，是可以得到更简单的方法的，不过前提是变换和逆变换必须是等价的。

从课程之外的角度来说，这是我第一次用C语言来解决一个问题，这个过程可以说是对C语言的一个重新学习，Matlab同样如此，学过之后就放下了，也很快就忘记了，这次虽然用到的不多，但是也让我更加体会到实践运用才是学习的最好方法；同时，完成课设报告的过程也让我学习了论文的格式和排版。

附录 C语言代码#include<stdio.h>#include<math.h>#define maxnum 128 //宏定义最大数：用以初始化数组，且限制了时域信号取样点数须小于128#define PI 3.1415926 //宏定义圆周率struct XKstruct//频域信号采样点结构体{double real; //实部double image; //虚部double absolutevalue; //绝对值double phaseangle; //相位角double radianmeasureangle; //用圆周率表示的角度};struct complexnum//复数结构体{double real; //实部double image; //虚部};struct complexnum complexmul(struct complexnum ,struct complexnum ); //函数声明语句，在定义处进行功能介绍void ComplexnumToXKStruct(int,struct complexnum*,struct XKstruct*);void DFTfunction(int ,double*,struct complexnum*);void FFTfunction(int ,double *,struct complexnum *);void xndistribute(int N,double *xn,double *an,double *bn);void initfunction();void display(int ,struct XKstruct*);void stardivision();void timedomainsignalsample(int ,double *);void main(){int i,N; //i为辅助变量，N为时域信号取样点数double xn[maxnum]; //xn[N]存放时域信号采样值struct complexnum XK[maxnum]; //XK[N]存放频域信号采样值struct XKstruct XK1[maxnum],XK2[maxnum]; //XK1[N]存放DFT变换后的频域信号采样值，XK2[N]存放FFT变换后的频域信号采样值//double w; //基波角频率//double a[maxnum]; //谐波分量系数initfunction(); //初始化函数，在定义出进行功能介绍while(1){printf("开始运行:\n");printf("1、输入时域信号取样点数N \n N="); //步骤1：输入时域信号取样点数Nscanf("%d",&N);printf("\n2、输入时域信号N点离散取样值，存于数组xn[N]\n"); //步骤2：依次输入时域信号N点离散取样值，存于数组xn[N]printf(" 选项1：直接输入离散取样值\n 选项2：输入正弦谐波分量信息，让计算机进行取样\n 选择（1/2）");scanf("%d",&i);while((i!=1)&&(i!=2)){printf("请重新选择：（1/2）");scanf("%d",&i);}switch(i){case 1:for(i=0;i<N;i++){printf(" xn[%d]=",i);scanf("%lf",&(xn[i]));}break;case 2:timedomainsignalsample(N,xn);break;}printf("\n3、进行DFT变换并显示变换结果\n "); //步骤3：进行DFT变换并显示变换结果DFTfunction(N,xn,XK); //DFT函数，将时域信号采样数组xn[N]进行DFT，在定义处进行功能介绍ComplexnumToXKStruct(N,XK,XK1); //将复数结构体XK转化成复频域信号取样点XK1display(N,XK1); //显示函数，在定义处进行功能介绍printf("\n4、进行FFT变换并显示变换结果\n "); //步骤4：进行FFT变换并显示变换结果switch(N) //根据N的值判定能否进行FFT，仅当N是以2为底的正指数数时才可进行，否则直接跳出{case -1: return;case 0:case 2:case 4:case 8:case 16:case 32:case 64:case 128:{FFTfunction(N,xn,XK); //FFT函数，将时域信号采样数组xn[N]进行FFT，在定义处进行功能介绍ComplexnumToXKStruct(N,XK,XK2); //将复数结构体XK转化成复频域信号取样点XK1display(N,XK2); //显示函数，在定义处进行功能介绍printf(""); //格式控制语句stardivision();printf("\n");stardivision();}break;default:{printf(" 错误：无法进行FFT变换!!\n 原因：时域信号取样点数N不是以2为底的正指数，或者N大于128\n\n");stardivision(); //格式控制语句printf("\n");stardivision();}}}}void stardivision() //格式控制函数：输出一整排星号作为分隔符{printf("******************************************************************** ***********\n");}void initfunction() //初始化函数{stardivision(); //输出一排星号，以下为学生信息printf(" 班级：1106班\n");printf(" 学号：U201111932\n");printf(" 姓名：曾超\n");stardivision();stardivision();//本程序解说语句，程序总共分为4个步骤，然后为重复printf("运行此程序后：\1、输入时域信号取样点数N \n\n \2、依次输入时域信号N点离散取样值，存于数组xn[N]\n\n \3、进行DFT变换并显示变换结果\n\n \4、进行FFT变换并显示变换结果\n\n \5、重复之前步骤\n");stardivision();stardivision();printf("注意事项： 1、欲进行DFT，N必须为正整数；\n\n\2、欲进行FFT变换，N必须是为2为底的正指数数；\n\n\3、欲直接退出，N必须为-1\n"); //注意事项stardivision();stardivision();}/*********************************************时域信号取样函数功能：交互输入时域谐波信号，输出时域取样值输入：取样点数N，取样值存储数组指针xn输出：取样值**********************************************/void timedomainsignalsample(int N,double *xn){/*变量说明：w：基波角频率T:基波周期t：时间变量maxharmanicorder：波形所含的最大谐波次数，必须为奇数例如：对于y=a[1]sin(wt)+a[3]sin(3wt)+a[5]sin(5wt)+……+a[11]sin(11wt)，maxharmanicorder=11a[maxharmanic]：存储谐波分量系数的数组，仅取其奇数项y:t变量的因变量，暂存取样值*/double y=0,a[maxnum],w=0,T=0,t=0;int maxharmanicorder=0,i,j;printf("\n 下面请初始化时域信号\n 表达式举例：y=a[1]sin(wt)+a[3]sin(3wt)+a[5]sin(5wt)……+a[11]sin(11wt)，请依据提示输入\n");printf(" 1、请输入基波角频率w（rad/s）\n w=");scanf("%lf",&w);T=2*PI/w;printf(" 2、请输入最大谐波次数maxharmanicorder,例如上面举例中的最大谐波次数maxharmanicorder=11\n maxharmanicorder=");scanf("%d",&maxharmanicorder);while((maxharmanicorder<=0)||(maxharmanicorder%2)==0){printf(" 错误：最大谐波次数应为正奇数，请重新输入\nmaxharmanicorder=");scanf("%d",&maxharmanicorder);}printf(" 3、请依次输入谐波分量系数\n");for(i=1;i<=maxharmanicorder;i+=2){printf(" %d次谐波分量系数a[%d]=",i,i);scanf("%lf",&(a[i]));}for(j=0;j<N;j++){t=T/N*j;for(i=1;i<=maxharmanicorder;i+=2){y+=a[i]*sin((i)*w*t);}xn[j]=y;}}/****************************************************显示函数输入：时域信号取样点数N，频域信号取样数组指针输出：以两种形式在命令行显示出频域信号取样点的值，分别是“实部+虚部”“绝对值*角度值”****************************************************/void display(int N,struct XKstruct *XK){int i;for(i=0;i<N;i++){printf(" XK[%d] = %f + i(%f)\n",i,XK[i].real,XK[i].image); //以“实部+虚部”形式显示频域信号取样点printf(" XK[%d] =(%f)exp(%fPIi):\n",i,XK[i].absolutevalue,XK[i].radianmeasureangle); //以“绝对值*角度值”形式显示频域信号取样点}/**********************************************************************格式转换函数功能：将复数结构体数组转化成复频域信号取样点结构体数组输入：时域信号取样点数N，复数结构体指针cnum，频域信号取样点指针XK**********************************************************************/void ComplexnumToXKStruct(int N,struct complexnum *cnum,struct XKstruct* XK){int i;for(i=0;i<N;i++){XK[i].real=cnum[i].real;XK[i].image=cnum[i].image;XK[i].absolutevalue=sqrt(XK[i].real*XK[i].real+XK[i].image*XK[i].image);XK[i].phaseangle=atan(XK[i].image/XK[i].real);XK[i].radianmeasureangle=XK[i].phaseangle/PI;}}/****************************************************复数乘法函数输入：两个复数结构体，不分先后返回值：两个复数的乘积的复数结构体****************************************************/struct complexnum complexmul(struct complexnum mul1,struct complexnum mul2){struct complexnum result;result.image=mul1.image*mul2.real+mul1.real*mul2.image;result.real=mul1.real*mul2.real-mul1.image*mul2.image;return result;}/****************************************************************************** ****//*DFT函数输入：时域信号取样点数NDFT，时域信号取样数组指针xnDFT，频域信号取样数组指针XKDFT功能：进行DFT，离散时域信号值存储在xnDFT[NDFT]，离散频域信号值存储在XKDFT[NDFT]*//****************************************************************************** *****/void DFTfunction(int NDFT,double *xnDFT,struct complexnum *XKDFT)int k,n; //辅助变量double *xnDFTtemp=xnDFT; //为了不改变原指针的值，定义了暂时指针struct complexnum *XKDFTtemp=XKDFT;for(k=0;k<NDFT;k++){(*XKDFTtemp).real=0; //初始化为0(*XKDFTtemp).image=0;for(n=0;n<NDFT;n++) //变换公式参见报告DFT定义{(*XKDFTtemp).real+=(*xnDFTtemp)*cos(2*PI*k*n/NDFT);(*XKDFTtemp).image+=-(*xnDFTtemp)*sin(2*PI*k*n/NDFT);xnDFTtemp++;}xnDFTtemp=xnDFT;XKDFTtemp++;}}/****************************************************************************** ****//*FFT函数输入：时域信号取样点数NFFT，时域信号取样数组指针xnFFT，频域信号取样数组指针XKFFT功能：进行DFT，离散时域信号值存储在xn[N]，离散频域信号值存储在XK[N]*//****************************************************************************** *****/void FFTfunction(int NFFT,double *xnFFT,struct complexnum *XKFFT){/*变量说明：coefficient指蝶形结运算系数coquotient指运算系数商，当前系数与之前一个系数的商complextemp为辅助变量，临时存储XKtemp用XKtemp：为不改变XKFFT指针值所设置的替代数组flag1、flag2：倒位序算法的变址标志Ntemp：为不改变NFFT值所设置的替代数i、j、k:辅助变量*/struct complexnum coefficient,coquotient,complextemp,XKtemp[maxnum];int flag2,flag1,Ntemp,level,step,length,lenght1,ip;int i,j,k;for(i=0;i<NFFT;i++){XKtemp[i].real=xnFFT[i];XKtemp[i].image=0;}//变址运算，即把自然顺序变成倒位序flag2=NFFT/2;flag1=NFFT-1;for(i=0,j=0;i<flag1;i++){if(i<j) //如果i<j,即进行变址{complextemp=XKtemp[j];XKtemp[j]=XKtemp[i];XKtemp[i]=complextemp;}k=flag2; //求j的下一个倒位序while(k<=j) //如果k<=j,表示j的最高位为1{j=j-k; //把最高位变成0k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0}j=j+k; //把0改为1}//使用蝶形运算完成FFT运算Ntemp=NFFT;for(level=1;(Ntemp=Ntemp/2)!=1;level++); //计算level的值，即计算蝶形总级数level=log（2）Nfor(step=1;step<=level;step++) // 控制蝶形结级数{ //step表示第step级蝶形length=2<<(step-1); //length蝶形结距离，即第m级蝶形的蝶形结相距length点lenght1=length/2; //同一蝶形结中参加运算的两点的距离coefficient.real=1; //coefficient为蝶形结运算系数，初始值为1coefficient.image=0;coquotient.real=cos(PI/lenght1); //coquotient为系数商，即当前系数与前一个系数的商coquotient.image=-sin(PI/lenght1);for(j=0;j<=lenght1-1;j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结{for(i=j;i<=NFFT-1;i=i+length) //控制同一蝶形结运算，即计算系数相同蝶形结{ip=i+lenght1; //i，ip分别表示参加蝶形运算的两个节点complextemp=complexmul(XKtemp[ip],coefficient); //蝶形运算，详见报告中定义XKtemp[ip].real=XKtemp[i].real-complextemp.real;XKtemp[ip].image=XKtemp[i].image-complextemp.image;XKtemp[i].real=XKtemp[i].real+complextemp.real;XKtemp[i].image=XKtemp[i].image+complextemp.image;}coefficient=complexmul(coefficient,coquotient); //改变系数，进行下一个蝶形运算}}for(i=0;i<NFFT;i++) //将XKtemp暂存的频域取样信号转存在XKFFT数组中{XKFFT[i]=XKtemp[i];}（本资料素材和资料部分来自网络，仅供参考。