搜索
华中科技大学信号与系统课程设论文快速傅里叶变换（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.5FFT算法时间抽取算法的基本思想是把时间间隔细分，细分后的时间间隔内包含的点数较少。的计算可以分成两部分，首先对符号进行简化，令,复数是单位1的N次开放，即：贓熱俣阃歲匱阊邺镓騷。假设N1，这样。根据的表示方法，N点DFT和DFT反变换为：现在令N是一个偶整数，以便N/2是一个整数。已知信号x[n],令a[n]为x[n]的偶数次项，b[n]为x[n]的级数次项。令和分别代表和的（N/2）点DFT，即：坛摶乡囂忏蒌鍥铃氈淚。令代表的N点DFT，可有：由上两式计算需要次乘法。为了看清这一点，首先注意到计算需要次乘法，和一样多，在上两式中计算需要N/2次乘法，所以，总的乘法次数为，比次乘法少次，因此，当N非常大时，可以大大减少乘法次数。蜡變黲癟報伥铉锚鈰赘。如果N/2是偶数，和又可以分别表示为两部分，进而重复上面的过程。如果q为正整数，这个递减的过程可以重复进行，知道信号只剩下一个非零值。買鲷鴯譖昙膚遙闫撷凄。下图给出了N=8时FFT算法的流程图（图1-1），在最左侧输入的是已知信号，需要注意信号输入的顺序，该顺序可以通过“倒位序”的方法来确定。假设，已知整数n的范围是0~N-1，时间标号n可以用q位二进制数来表示，把这q位二进制数进行倒位序排列，记得FFT算法每一行的输入信号。表1-1给出了图1-1中FFT算法的输入信号值的顺序。綾镝鯛駕櫬鹕踪韦辚糴。图1-1N=8时的FFT算法流程图表1-1N=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]存的是倒位序。IJ的时候需要变序，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’=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.1时域函数为，取样点数N=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，时域信号函数file，此时，maxharmanicorder=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}尧侧閆繭絳闕绚勵蜆贅。时域信号函数为时结果为：时域信号函数为时结果为：时域信号函数为时结果为：幅度谱（从左至右高次谐波分量有增加）：由上图可知，随着谐波分量的增加，幅度谱的个别取样点幅值有变化，由于谐波分量增量很小，因此幅值变化也很小，同时，整体波形变化也不大。识饒鎂錕缢灩筧嚌俨淒。由上图分析可知，随着谐波分量的增加，相位谱完全没有变