快速傅里叶变换（FFT）的计算机实现_信号与系统课程设计论文1.docx

华中科技大学信号与系统课程设论文快速傅里叶变换（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)在式（1-2）中，、和都是实数，是基波频率（rad/s），且=/，是周期。系数和可由下面的公式计算出来：, (1-3), (1-4)应该注意的是上面给定的和可以在任何周期内的积分计算出来。例如 , (1-5)在式（1-2）中，项是常数，或者是的直流成分，由式（1-6）确定： (1-6)表达式（1-2）叫做周期信号的三角傅里叶级数。的一次谐波项为，二次谐波项为，第k次谐波项为。注意，组成的谐波频率为整数倍。这是周期信号的重要特性。由式（1-2）给定的三角傅里叶级数可以在写成余弦函数和相位的形式： , (1-7)式中： (1-8)并且： (1-9)1.2.3复指数级数由式（1-2）给定的三角傅里叶级数可以表示成如下的指数形式： (1-10)在式（1-10）中，是实数，对于，一般是复数。是基波的频率（rad/s），且，是基本周期。式（1-10）的复指数级数的系数可由如下公式计算得出： (1-11)应该注意到，式（1-11）给定的也能在任何整周期上的积分计算出来，例如： (1-12)1.3 傅里叶变换及其直角和极坐标表达式给定一个信号，其傅里叶变换被定成频率函数： (1-13)式中，是连续频率分量。一般意义上，（1-13）式中的积分收敛（存在），我们则说信号有傅里叶变换。如果是“表现良好”的和绝对可积的，那么积分就是收敛的。绝对可积是指下式成立： (1-14)“表现良好”就是信号在任何有限的时间段内存在有限的间断点和最大、最小值。除了脉冲信号，我们感兴趣的大部分信号都是“表现良好”的。所有实际信号（即物理上可以产生的信号）都是“表现良好”的，且满足（1-14）式。如前所述，是实变量的复值函数，换句话说，若给定，则通常是复数。因为复数既可以表示成直角坐标形式，也可以表示成极坐标形式，所以，傅里叶变换也可以表示成这两种形式。下面将定义这两种形式。由欧拉公式，可以写成如下形式：令和是的实值函数，定义如下：则的直角坐标表达式为：的极坐标表达式，由下式给出： (1-15)式中，是的模，是的辐角。利用下面的关系，可以由直角坐标表达式换成极坐标表达式：1.4 离散时间信号的傅里叶分析1.4.1离散时间傅里叶变换dtft已知一个离散时间信号,他的离散时间傅里叶变换（dtft）定义为： (1-16)一般地，由上式定义的dtft是实变量的复值函数。对于所有的实值，如果上式中的双边无限求和收敛，则称离散时间信号有一般意义的dtft。有一般意义的dtft的充分条件是绝对可和，即：如果是时限的离散时间信号，显然，上式中的和是有限值，类似这样的信号都存在一般意义下的dtft。-已知离散时间信号，其dtft为，由于一般是复数，可以用直角坐标或极坐标的形式来表示。利用欧拉公式得到的的直角坐标形式：其中：的极坐标形式是:其中：1.4.2离散傅里叶变换是一个离散时间信号，其dtft为，由于是连续变量的函数，除非可以表示为封闭形式，否则不能保存在数字计算及的存储器中。为了能在数字计算及上实现dtft，必须在频率上离散化，从而引出了下面定义的离散傅里叶变换的概念。假设对于所有的整数，离散时间信号等于零，是一个固定的正整数。整数可能很大，例如。的点离散傅里叶变换（dft）定义为： (1-16)由式（1-16）可见，dft是离散变量k的函数。可表示成极坐标形式或直角坐标形式。极坐标形式是： 直角坐标形式是：其中：1.5 fft算法时间抽取算法的基本思想是把时间间隔细分，细分后的时间间隔内包含的点数较少。的计算可以分成两部分，首先对符号进行简化，令,复数是单位1的n次开放，即： 假设n1，这样。根据的表示方法，n点dft和dft反变换为：现在令n是一个偶整数，以便n/2是一个整数。已知信号xn,令an为xn的偶数次项，bn为xn的级数次项。令和分别代表和的（n/2）点dft，即：令代表的n点dft，可有：由上两式计算需要次乘法。为了看清这一点，首先注意到计算需要次乘法，和一样多，在上两式中计算需要n/2次乘法，所以，总的乘法次数为，比次乘法少次，因此，当n非常大时，可以大大减少乘法次数。如果n/2是偶数，和又可以分别表示为两部分，进而重复上面的过程。如果q为正整数，这个递减的过程可以重复进行，知道信号只剩下一个非零值。下图给出了n=8时fft算法的流程图（图1-1），在最左侧输入的是已知信号，需要注意信号输入的顺序，该顺序可以通过“倒位序”的方法来确定。假设，已知整数n的范围是0n-1，时间标号n可以用q位二进制数来表示，把这q位二进制数进行倒位序排列，记得fft算法每一行的输入信号。表1-1给出了图1-1中fft算法的输入信号值的顺序。图1-1 n=8时的fft算法流程图表1-1 n=8时的倒位序排列时间（n）二进制代码倒位序代码排序0000000x01001100x42010010x23011110x64100001x15101101x56110011x37111111x71.5.1 倒位序算法在进行fft运算时，第一步要实现的就是倒位序排列。假设使用ai存的是顺序位序，而bj存的是倒位序。ij的时候就不用。不然就相当于作了两次变序，又变回去了。从表1-1可知，按自然顺序排列的二进制数，其下面的数总是比其上面的数大1，即下面的数是上面的数在最低位加1并向高位进位而得到的。而倒位序二进制数的下面的数是上面的数在最高位加1并由高位向低位进位而得到。由数组的性质可知，i、j都是从0开始，若已知某个倒位序j，要求下一个倒位序数，则应先判断j的最高位是否为0，这可与flag=n/2相比较。如果flagj，则j的最高位为0，只要把该位变为1（j与flag=n/2相加即可），就得到下一个倒位序数；如果flag=j，则j的最高位为1， 可将最高位变为0（j与flag=n/2相减即可）。然后还需要判断次高位，这可与flag=n4相比较，若次高位为0，则需将它变为1（加n4即可）其他位不变， 既得到下一个倒位序数；若次高位是1，则需将它也变为0。依此类推即可得到最后的倒位序排列。2 c语言实现见附录3 利用c程序进行频谱分析3.1直接输入离散信号取样值进行dft和fft3.1.1取样点n=4，时域信号取样值xn=1，2，2，1结果见下图可见fft和dft的结果基本上是一样的，仅在小数点第六位才有一点差别，见xk3。3.1.2取样点n=5，时域信号取样值xn=1，2，2，4，5结果见下图当n=时，dft仍然可以进行，因为dft的代码“翻译”自dft的定义而来，而fft的代码是“翻译”自倒位序算法和蝶形算法，这两种算法对取样点数有要求，n必须是以2为底的正指数。同时，此c程序还对取样点最大值有要求，不得超过128.3.2输入正弦谐波分量信息，让计算机进行取样3.2.1时域函数为，取样点数n=4、8、16此时， ，最大谐波次数maxharmanicorder=5，谐波分量系数为a1=9,a3=3,a5=1。n=4时结果见下图：n=8时结果见下图：n=16时，结果见下图（仅出示fft）幅度谱图（从左至右n=4、6、8）：幅度谱分析：由上图可知，随着取样点数的增加，边界点的幅值急剧上升，不过图像的整体趋势不变相位谱（从左至右n=4、6、8）：由上图可知，随着取样点数的增加，相位谱的变化趋势没有发生变化，只是将原本的取样点变得更密集了而已。3.2.2取样点数n=8，时域信号函数file，此时，maxharmanicorder=5，a1，3，5=9，3，1，此时，maxharmanicorder=7，a1，3，5，7=9，3，1，0.5，此时，maxharmanicorder=9，a1，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#include#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 xnmaxnum;/xnn存放时域信号采样值struct complexnum xkmaxnum;/xkn存放频域信号采样值struct xkstruct xk1maxnum,xk2maxnum;/xk1n存放dft变换后的频域信号采样值，xk2n存放fft变换后的频域信号采样值/double w;/基波角频率/double amaxnum;/谐波分量系数initfunction();/初始化函数，在定义出进行功能介绍while(1)printf(开始运行:n);printf(1、输入时域信号取样点数n nn=);/步骤1：输入时域信号取样点数nscanf(%d,&n);printf(n2、输入时域信号n点离散取样值，存于数组xnnn);/步骤2：依次输入时域信号n点离散取样值，存于数组xnnprintf(选项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;in;i+)printf(xn%d=,i);scanf(%lf,&(xni);break;case 2:timedomainsignalsample(n,xn);break;printf(n3、进行dft变换并显示变换结果n );/步骤3：进行dft变换并显示变换结果dftfunction(n,xn,xk);/dft函数，将时域信号采样数组xnn进行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函数，将时域信号采样数组xnn进行fft，在定义处进行功能介绍complexnumtoxkstruct(n,xk,xk2);/将复数结构体xk转化成复频域信号取样点xk1display(n,xk2);/显示函数，在定义处进行功能介绍printf();/格式控制语句stardivision();printf(n);stardivision();break;default:printf(错误：无法进行fft变换!n原因：时域信号取样点数n不是以2为底的正指数，或者n大于128nn);stardivision();/格式控制语句printf(n);stardivision();void stardivision()/格式控制函数：输出一整排星号作为分隔符printf(*n);void initfunction()/初始化函数stardivision();/输出一排星号，以下为学生信息printf( 班级：1106班n);printf( 学号：u201111932n);printf( 姓名：曾超n);stardivision();stardivision();/本程序解说语句，程序总共分为4个步骤，然后为重复printf(运行此程序后：1、输入时域信号取样点数n nn 2、依次输入时域信号n点离散取样值，存于数组xnnnn 3、进行dft变换并显示变换结果nn 4、进行fft变换并显示变换结果nn 5、重复之前步骤n);stardivision();stardivision();printf(注意事项：1、欲进行dft，n必须为正整数；nn2、欲进行fft变换，n必须是为2为底的正指数数；nn3、欲直接退出，n必须为-1n);/注意事项stardivision();stardivision();/*时域信号取样函数功能：交互输入时域谐波信号，输出时域取样值输入：取样点数n，取样值存储数组指针xn输出：取样值*/void timedomainsignalsample(int n,double *xn)/*变量说明：w：基波角频率t:基波周期t：时间变量maxharmanicorder：波形所含的最大谐波次数，必须为奇数例如：对于y=a1sin(wt)+a3sin(3wt)+a5sin(5wt)+a11sin(11wt)，maxharmanicorder=11amaxharmanic：存储谐波分量系数的数组，仅取其奇数项y:t变量的因变量，暂存取样值*/double y=0,amaxnum,w=0,t=0,t=0;int maxharmanicorder=0,i,j;printf(n 下面请初始化时域信号n表达式举例：y=a1sin(wt)+a3sin(3wt)+a5sin(5wt)+a11sin(11wt)，请依据提示输入n);printf(1、请输入基波角频率w（rad/s）nw=);scanf(%lf,&w);t=2*pi/w;printf(2、请输入最大谐波次数maxharmanicorder,例如上面举例中的最大谐波次数maxharmanicorder=11nmaxharmanicorder=);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,&(ai);for(j=0;jn;j+)t=t/n*j;for(i=1;i=maxharmanicorder;i+=2)y+=ai*sin(i)*w*t);xnj=y;/*显示函数输入：时域信号取样点数n，频域信号取样数组指针输出：以两种形式在命令行显示出频域信号取样点的值，分别是“实部+虚部”“绝对值*角度值”*/void display(int n,struct xkstruct *xk)int i;for(i=0;in;i+)printf(xk%d = %f + i(%f)n,i,xki.real,xki.image);/以“实部+虚部”形式显示频域信号取样点printf(xk%d = (%f)exp(%fpii):n,i,xki.absolutevalue,xki.radianmeasureangle);/以“绝对值*角度值”形式显示频域信号取样点/*格式转换函数功能：将复数结构体数组转化成复频域信号取样点结构体数组输入：时域信号取样点数n，复数结构体指针cnum，频域信号取样点指针xk*/void complexnumtoxkstruct(int n,struct complexnum *cnum,struct xkstruct* xk)int i;for(i=0;in;i+)xki.real=cnumi.real;xki.image=cnumi.image;xki.absolutevalue=sqrt(xki.real*xki.real+xki.image*xki.image);xki.phaseangle=atan(xki.image/xki.real);xki.radianmeasureangle=xki.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，离散时域信号值存储在xndftndft，离散频域信号值存储在xkdftndft*/*/void dftfunction(int ndft,double *xndft,struct complexnum *xkdft)int k,n;/辅助变量double *xndfttemp=xndft;/为了不改变原指针的值，定义了暂时指针struct complexnum *xkdfttemp=xkdft;for(k=0;kndft;k+)(*xkdfttemp).real=0;/初始化为0(*xkdfttemp).image=0;for(n=0;nndft;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，离散时域信号值存储在xnn，离散频域信号值存储在xkn*/*/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,xktempmaxnum;int flag2,flag1,ntemp,level,step,length,lenght1,ip;int i,j,k;for(i=0;infft;i+)xktempi.real=xnffti;xktempi.imag