《实序列FFT算法的C语言实现.doc》由会员分享，可在线阅读，更多相关《实序列FFT算法的C语言实现.doc（41页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上实序列FFT算法的C语言实现学生：XX 指导教师：XX内容摘要：DFT和IDFT是数字信号分析与处理中的一种重要运算和变换，但直接根据定义计算DFT时，运算量大，不能实时得到计算结果。特别是在实际应用中，N都取得比较大，此时，由于乘法和加法的次数都近似与N的平方成正比，这将极大增加DFT计算所需时间。为此，提出了许多DFT和IDFT的快速算法，称为快速傅里叶变换（FFT）和快速傅里叶反变换（IFFT）。本文较为系统地阐述了快速傅里叶变换的算法原理然后用MATLAB实现了快速傅里叶变换。论文首先首要介绍了FT与DFT的定义、DFT与FFT的关系, 然后重点介绍基2时域抽
2、取FFT算法以及其原理和运算流图，再应用C语言实现了实序列的FFT。最后在Matlab软件上进行仿真，仿真结果验证了设计的正确性。关键词： 快速傅立叶变换 Matlab 仿真Realization of FFT algorithm for real sequenceWith C programAbstract：DFT and IDFT are important transform and processing in digital signal processing. However, there are large amount of computation by directly cal
3、culating according to the definition of DFT. Especially in the practical application, N is bigger, at this time, because the time of multiplication and addition are approximately proportional to the square of N, which will greatly increase the calculation time needed for DFT. Therefore, many DFT and
4、 IDFT fast algorithm are raised, which are called FFT and IFFT.In this paper relatively systematically elaborated the fast Fourier transform algorithm principle and use MATLAB software to realize the fast Fourier transform. The paper first introduces the definition of FT and DFT,the relationship bet
5、ween DFT and FFT, and then mainly introduces DIT-FFT ,including its principle and operation flow diagram, and finally used C language to realize the real sequence FFT.The designs are simulated in Matlab software, the results of the simulation confirm the exactness of the design.Keywords：Fourier tran
6、sformation fast Fourier transformation Matlab 目 录专心-专注-专业实序列FFT算法的C语言实现前言在实际的数字系统中，DFT是一种得到了广泛的应用的、重要的信号处理手段，但它的运算效率非常低。随着DFT输入的点数增加到数百或数千，DFT需要的运算量变得非常大。快速傅里叶变换（FFT）可使实现DFT的运算量下降几个数量级，从而使数字信号处理的速度大大提高1。FFT是离散傅立叶变换（DFT）的快速算法，可以将一个信号变换到频域。有些信号在时域上是不易看出有什么特征的，但是如果变换到频域之后，就很容易了。FFT可以将一个信号的频谱提取出来，这在频谱分析
7、方面也是经常用的。实际中需要做快速傅里叶变换的多为实序列数据，而其变换算法都是以复数序列作为输入。本设计利用频域的性质，将实序列数据变为复数序列，再进行FFT,以提高效率。本设计就是要求在熟悉DFT及FFT算法基本原理的基础上，编制C语言程序实现实数序列的FFT算法。1 序列的FT和DFT序列的傅里叶变换（Fourier Transform，FT）和离散傅里叶变换（DFT）都是对序列的频域描述，它们揭示了序列由那些分量构成，各分量的幅度和相位大小。1.1 序列的FT序列x(n)的傅里叶变换又称为离散时间傅里叶变换（DTFT），其定义为 （1.1-1）式中，w 称为数字角频率。如果已知x(n)的
8、傅里叶变换X(ejw)，则可下式求得其时域表达式 （1.1-2）上式称为序列的傅里叶反变换（IFT）。序列的傅里叶变换是w 的连续函数，且一般情况下为复变函数，并可表示为其中，和分别称为幅度谱和相位谱。此外，根据的周期性可知，序列的频谱及其幅度谱和相位谱都是关于w 以2p为周期的周期函数。1.2 序列的DFT序列的FT变换为w 的连续函数，不便于用计算机程序进行辅助分析和计算。为了便于用计算机辅助求解和分析，提出了离散傅里叶变换（DFT）。1.2.1 DFT的定义和计算长度为M的有限长序列x(n)，n=0M-1，其N点DFT和IDFT（离散傅里叶反变换）分别定义为，k=0N-1 （1.2.1-
9、1），n=0N-1 （1.2.1-2）式中，称为旋转因子。借助于矩阵，可以将以上两式所示序列的N点DFT和IDFT运算用矩阵表示为 （1.2.1-3） （1.2.1-4）1.2.2 实序列的DFT 实际系统中所处理的信号大多数是实序列，对这些实序列x(n)，根据式（1.2.1-1）得到则因为x(n)为实序列，则。此外，有因此或者 （1.2.1-5）上式说明，所有实数序列的DFT都具有共轭对称性。因此，对实序列，只需要计算出其N点DFT中的前面N/2个点，即k=0N/2-1的点，即可根据共轭对称性，由式（1.2.1-5）求出另外N/2个点，即k=N/2N-1的点，合起来得到完整的N点DFT。在式
10、（1.2.1-5）中分别令k=0和N/2，又可得到以上两式说明，在实序列的N点DFT中，当k=0和N/2时一定为实数。2 FFT算法DFT和IDFT是数字信号分析与处理中的一种重要运算和变换，但直接根据定义计算DFT时，运算量大，不能实时得到计算结果。特别是在实际应用中，N都取得比较大，此时，由于乘法和加法的次数都近似与N的平方成正比，这将极大增加DFT计算所需时间。为此，提出了许多DFT和IDFT的快速算法，称为快速傅里叶变换（FFT）和快速傅里叶反变换（IFFT）2。FFT和IFFT算法的基本思想是根据旋转因子WNm的对称性，将N点DFT分解几个较短的DFT，从而减少乘法次数，缩短计算时间
11、。根据这一思想，出现了许多DFT和IDFT的快速算法。下面结合本设计重点介绍基2时域抽取FFT算法。2.1 基2时域抽取FFT算法时域抽取FFT算法（DIT-FFT）的基本思想是在时域中将N点序列x(n)中的各点按序号的奇偶进行抽取分组，得到两个N/2点序列。在计算出两个N/2点序列的DFT后，通过蝶形运算得到原序列的N点DFT3。2.1.1 基本原理假设序列x(n)的长度为N=2M，其中M为正整数。根据序号n的奇偶将x(n)分解为两个N/2点的子序列，即则令分别为序列x1(r)和x2(r)的N/2点DFT，则得到，k=0N/2-1 （2.1.1-1）由于X1(k)和X2(k)都以N/2点为周
12、期，且因此式（2.1.1-1）可表示为 ，k=0N/2-1 （2.1.1-2）上式代表的运算称为蝶形运算，可用图2.1.1-1所示的流图符号表示。例如，令上式中的k=1，则得到一个蝶形运算为通过上述过程，就将序列x(n)的N点DFT分解为两个N/2点的DFT，在得到X1(k)和X2(k)以后，通过N/2个蝶形运算即可得到X(k)。ABCA+BCA-BC图2.1.1-1 蝶形运算符号重复上述分解过程，继续将序列x1(n)和x2(n)按序号的奇偶进行抽取分解，直到最终分解为N/2个两点序列。通过每次抽取分解，都将序列和DFT的点数缩短一半，同时转换为若干蝶形运算。2.1.2 DIT-FFT算法的运
13、算流图根据上述分解过程，容易知道，对N点DFT，只需通过M=log2N次抽取，就可以将原来的N点DFT分解为N个1点DFT和M级蝶形运算，而1点DFT就是序列本身。以，从而得到完整的蝶形运算流图如图2.1.2-1所示。x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(0)X(1)X(2)X(3)X(4)X(5)X(6)X(7)图2.1.2-1 8点DIT-FFT运算流图2.1.3 DIT-FFT算法的运算量和存储量根据上述分解过程，容易知道，对N点DFT，只需通过M=log2N次抽取，就可以将根据上述DIT-FFT的基本原理，通过时域抽取，将N点DFT转换为M级蝶形运算，每级
14、包括N/2个蝶形运算，因此N点FFT共包括MN/2个蝶形运算。每个蝶形运算需要一次乘法运算和两次加法运算，则N点FFT的运算量为：乘法：加法：而DFT直接算法共需N2次乘法和N(N-1)次加法。当N1时，采用FFT算法将大大减少运算次数，提高运算速度。由图2.1.2-1所示运算流图可知，DIT-FFT的运算过程很有规律。在同一级的N/2个蝶形运算中，每个蝶形运算的两个输入输出数据只对该蝶形有用，而且都位于同一条水平线上。因此，每个蝶形的两个输出数据可以只占用输入数据的存储单位，而不要另外分配存储单元。这种存储计算数据的方法称为原位计算。原位计算可以极大地节省存储单元。对N点FFT算法，只需要N
15、个存储单元。运算一开始，这N个单元用于存储原始序列。计算过程中，用于存储各蝶形运算的输出数据。运算完毕后，这N个单元中存储的就是输出N点FFT4。2.2 实序列的FFT算法实际系统中处理的序列都是实数序列。如果直接按照前述的FFT算法进行计算，就是将其视为虚部都为0的复数序列，这就增加了存储量，降低了运算效率。处理该问题的方法有很多。例如，通过一个N点FFT同时计算出两个实数序列的N点FFT；根据基2时域抽取FFT算法的基本思想，用一个N/2点FFT计算出一个实序列的N点FFT5。本设计采用另外一种方法，以减少运算量和存储量。已经知道，实序列的mDFT都具有共轭对称性，或者说，其实部都是偶对称
16、的、虚部都是奇对称的。利用这种共轭对称性，计算时可以只计算出序列的N点DFT中前面N/2+1个点，即X(0)X(N/2)。在得到这些点后，利用共轭对称性，通过简单的替换就可得到完整的N点DFT。采用这种方法计算时，不需要存储后面N/2-1个点，即X(N/2+1)X(N-1)。此外，实序列的N点DFT中，X(0)和X(N/2)一定为实数，因此计算过程中也不需要存储其虚部。这就极大地减少了存储量。以实序列的8点DFT为例，根据共轭对称性有：因此，只需计算出X(0)X(4)，即可根据以上各式得到X(5)X(7)，进而得到完整的8点DFT。此外，对实数序列的8点DFT，X(0)和X(4)一定是实数，因
17、此计算过程中只需为X(0)和X(4)分别分配一个存储单元。3 实序列FFT算法的C语言实现根据设计要求，利用VS2010软件编制了C语言程序，实现实序列的基2时域抽取FFT算法。3.1 VS2010简介VS2010提供了两类容器，以便有效地管理开发工作所需的项，如引用、数据连接、文件夹和文件。这两类容器分别叫做解决方案和项目。此外，VS2010还提供解决方案文件夹，用于将相关的项目组织成项目组，然后对这些项目组执行操作。 作为查看和管理这些容器及其关联项的界面。 为了使集成开发环境 (IDE) 能够应用它的各种工具、设计器、模板和设置，Visual Studio 2010（简称VS2010）实
18、现了概念上的容器（称为解决方案和项目）。 另外，Visual Studio 还提供了解决方案文件夹，用于将相关的项目组织成组，然后对这些项目组执行操作。 解决方案管理 Visual Studio 配置、生成和部署相关项目集的方式。解决方案可以只包含一个项目，也可以包含由开发团队联合生成的多个项目。复杂的应用程序可能需要多个解决方案。项目包含一组源文件以及相关的元数据，如组件参考和生成说明。 生成项目时通常会生成一个或多个输出文件。根据实际需要，项目可以简单，也可以复杂。一个简单的项目可能由一个窗体或 HTML 文档、源代码文件和一个项目文件组成。更加复杂的项目可能由这些项以及数据库脚本、存储过
19、程和对现有 XML Web services 的引用组成。 创建新项目时，VS2010会自动生成一个解决方案。然后，可以根据需要将其他项目添加到该解决方案中。也可以创建不包含项目的空白解决方案，从而使用 VS2010编辑器和设计器修改独立的文件。因为每个项目或解决方案都由一个目录及其内容组成，所以，可以在 Windows 资源管理器中移动、复制或删除解决方案和项目。VS2010将解决方案的定义存储在两个文件中，即.sln 和 .suo。这两个文件相当于早期版本中的工作区文件 (.dsw)。解决方案定义文件 (.sln) 存储定义解决方案的元数据，每当解决方案活动时，都使用构建该解决方案并设置其
20、属性时存储在 .suo 文件中的元数据来自定义 IDE。3.1.1 新建项目为了调试本设计所要求的程序，首先启动VS2010，然后依次选中文件菜单中的新建、项目命令，打开如图3.1.1-1所示的新建项目对话框。图3.1.1-1 新建项目对话框在该对话框中，选中左边Visual C+下面的Win32，在右边选中Win32控制台应用程序。然后在名称框中输入项目名称。这里设定项目名称为rfft，并设定项目存放的位置为E盘。注意，一旦设定项目名称后，解决方案的名称默认与项目名称相同。单击确定按钮，在后面出现的个对话框中都选默认值，在应用程序设置对话框中确保选中“空项目”，之后，即可创建项目rfft，同
21、时创建一个同名的解决方案。现在打开E盘，可以看到一个rfft文件夹，其中有一个rfft.sln文件，代表解决方案，同时还有一个rfft子文件，其中有一个文件rfft.vcxproj，代表rfft项目。3.1.2 新建文件为了编辑源程序代码，在VS2010的文件菜单中选中新建、文件命令，在弹出的新建文件对话框中选中Visual C+、C+文件。单击打开按钮后，可以看到在VS2010的工作窗口中出现一个新文件“源1.cpp”。选中菜单中的另存为命令，将其重新保存为C语言源程序文件，这里命名为rfft.c，用于新建实序列FFT算法子程序。注意在另存为对话框中一定要选中C源程序选项，并找到文件夹rff
22、t，将文件保存到该文件中。重复上述步骤，再新建一个C语言源程序文件，命名为mrfft.c，该文件是主程序文件，用于调用rfft.c文件中定义的rfft函数，计算实序列的N点DFT。分别在mrfft.c和rfft.c文件中编辑实序列FFT算法的子程序和主程序，之后，为了将文件添加到项目中，在VS2010工作窗口左边的解决方案资源管理器中右键单击项目名称rfft，在弹出的菜单中依次选中添加、现有项命令，找到刚才创建的文件mrfft.c，单击添加按钮，将其添加到项目和解决方案中。将主程序文件添加到解决方案后，通过生成菜单中的解决方案命令，对程序进行编译连接。如果没有错误，则生成解决方案，并在rfft
23、文件夹下自动生成一个Debug文件夹，可以看到其中有一个看到其中有一个rfft.exe文件。此时，VS2010窗口中的解决方案资源管理器如图3.1.2-1所示。项目rfft下面的mrfft.c文件是手动加入的，而外部依赖项下面的所有文件（包括子程序文件rfft.c）是通过编译连接自动加入的。图3.1.2-1 解决方案资源管理器3.2 实序列FFT算法子程序根据上述原理编制的实序列FFT算法子程序参见附录1。程序中，首先将输入序列进行倒序排列，然后进行各级蝶形运算。3.2.1 倒序由FFT的运算流图可知，在进行FFT运算之前，首先需要将按自然顺序送入的输入序列按照一定的顺序排好。例如，对N=8点
24、FFT，原始输入序列x(0)、x(1)、x(2)、x(7)，计算时的顺序应为x(0)、x(4)、x(2)、x(7)。这一操作称为倒序。以8点FFT为例，输入序列序号的自然顺序和倒序如表3.2.1-1所示。表3.2.1-1 顺序和倒序对照表顺 序倒 序十进制二进制十进制二进制0000000010014100201020103011611041001001510151016110301171117111由表3.2.1-1可知，自然顺序序号加1，是在其二进制数最低位加1，并逢2向高位进1。而倒序序号是在其对应的二进制数最高位加1，逢2向低位进位。对于N=2M，M位二进制数各位的权值从左向右依次为N/
25、2、N/4、2、1。因此，在最高位加1，相当于加N/2。如果最高位为0，即序号小于N/2，则直接加N/2得到下一个倒叙序号；否则，先将最高位变为0，然后次高位加1，相当于加N/4。在次高位加1时，同样需要判断该位是0还是1，再用同样的方法进行加1运算。依此类推，直到完成最高位加1，逢2向低位进位的运算。形成倒序序号后，将输入序列x(n)重新按倒序排列。实现倒序的程序框图如图3.2.1-1所示。根据上述框图编制的实现倒序的程序代码如下：i=0:N-1n1N-1ij?x(i)和x(j)交换kN/2kj+1?jj-kkk/2jj+kYNNY图3.2.1-1 倒序程序框图n1=n-1;for(j=0,
26、i=0;in1;i+)if(ij)xtmp=xj;xj=xi;xi=xtmp;k=n/2;while(k(j+1)j=j-k;k=k/2;j=j+k;代码中，i和j分别为自然序号和倒序序号，设置i从0到N-1，共循环N次。在每次循环中，当ij时，将x(i)和x(j)交换，否则不交换。代码中之后的while循环语句用于实现倒序值j的计算，以便下次循环实现相应一对数据的交换。3.2.2 蝶形运算将输入序列倒序后，进行各级蝶形运算。程序中首先用如下语句：for(j=1,i=1;i16;i+)m=i;j=2*j;if(j=n) break;求得N点FFT蝶形运算的级数。之后，进行第一级蝶形运算，再用k
27、控制共循环m-1次，依次进行后面各级蝶形运算。由FFT算法的运算流图可知，对第一级中的每个蝶形运算，旋转因子都为=1。再考虑到输入序列各点都是实数，因此将序列中相邻每两输入数据直接相加减，即得到第一级各蝶形的输出。相应的代码如下： for(i=0;in;i+=2)xtmp=xi;xi=xtmp+xi+1;xi+1=xtmp-xi+1; 在后面各级蝶形运算中，由运算流图可知，有些蝶形所需的旋转因子仍然为=1。对这些旋转因子对应的蝶形运算，同样化简为将其两个输入数据直接相加减，以得到该蝶形的两个输出。实现这一运算的主要代码如下：xtmp=xi;xi=xtmp+xi+n2;xi+n2=xtmp-xi
28、+n2;其中，n2=2k-1，k=2m。例如，对第二级蝶形，k=2，则n2=2。因此利用以上三条语句将第二级蝶形中对应的各蝶形运算的输入输出直接加减，结果再放回原单元。对第三级蝶形，k=3，则n2=4。利用以上三条语句将第三级蝶形中对应的各蝶形运算的输入输出直接加减，结果再放回原单元。对8点FFT，第二级有两个这样的蝶形，第三级有一个这样的蝶形。程序中用i和n1控制，其中n1=2k，i每次循环递增n1。对第二级，n1=4，则每次循环i从0到7递增4，因此上述三条语句共执行2次，实现的运算是：x(0)=x(0)+x(2);x(2)=x(0)-x(2)x(4)=x(4)+x(6);x(6)=x(4
29、)-x(6)对第三级，n1=8，因此上述语句只执行1次，实现的运算是：x(0)=x(0)+x(4);x(4)=x(0)-x(4)程序中还有如下语句：xi+n2+n4=-xi+n2+n4;该语句的作用是将蝶形运算中的某些数据直接取负。例如，对8点FFT，在进行第二级运算时，将x(3)和x(7)取负；在进行第三级运算时，将x(6)取负。后面将解释为什么要这样做。程序中后面的语句用于执行旋转因子不为时对应的其它蝶形运算。用j控制从1（n4-1），共循环n4-1次，其中n4=2k-2。对第二级，n4=1，因此该循环一次也不执行。对第三级，n4=2，因此该循环执行1次。为了解释该循环中各语句的功能，以8
30、点FFT为例，推导第三级蝶形运算的输出结果。以旋转因子WN1对应的蝶形为例，该蝶形中的一个输出为 （3.2.2-1）其中下标2和3分别代表第2和第3级蝶形的输出。根据运算流图，第2级的输出中，有将以上两式代入式（3.2.2-1）得到其中，则 （3.2.2-2）上式中，第一级蝶形的输出X1(1)、X1(3)、X1(5)和X1(7)都是由原序列通过直接加减得到的，因此都为实数。由此得到第三级的输出中，X3(1)的实部和虚部分别为 （3.2.2-3）式中，C1和S1分别为的实部和虚部。程序中的如下两条语句：t1=wr*xa3+wi*xa4;t2=wi*xa3-wr*xa4;执行时，a3=5，a4=7
31、；wr和wi分别为的实部和虚部的负，即wr=C1，wi=-S1。因此，如果令xa4=-x(7)，则t1和t2就分别等于式（3.2.2-3）中两个中括号中的项。这也就解释了为什么在执行第二级运算时需要将x(7)取负。程序中后面有如下两条语句：xa4=xa2-t2;xa1=xa1+t1;执行时，a4=7，a1=1，a2=3，则这两条语句实现的运算分别为 （3.2.2-4）将以上两式与式（3.2.2-3）比较可知，如果令X2(3)=-X1(3), X2(1)=X1(1)，则X3(1)和X3(7)就分别等于X3r(1)和X3i(1)。因此在第二级运算时，将X(3)取负，再由这两条语句计算，得到的输出刚
32、好为最后N点FFT结果中X(1)的实部和虚部，并且X(1)的实部和虚部分别存放在序号为1和7的单元。根据同样的分析，程序中另外两条语句xa3=-xa2-t2;xa2=xa1-t1;用于计算最后结果中的X(3)，并且将其实部和虚部分别存放在序号为3和5的单元。对8点FFT，通过循环执行上述4条语句1次，即得到结果中的X(1)和X(3)，再结合前面的语句，至此就得到了8点FFT中的X(0)、X(1)、X(3)和X(4)。对8点FFT中的X(2)和X(6)，根据运算流图得到由于在第二级蝶形运算时，X2(2)和X2(6)对应蝶形中的旋转因子都为WN0，因此X2(2)和X2(6)都为实数，则根据上式得到
33、因此，在第三级中，并不需要计算X(2)，只需将原序号为6的单元取负，即可得到X(2)的虚部，并存回序号为6的单元，而X(2)的实部通过第二级计算，存放在序号为2的单元。总结上述分析，各结果存放的单元序号如表3.2.2-1所示。表中数据的下标r和i分别代表实部和虚部。表3.2.2-1 计算结果存储顺序单元序号存放结果单元序号存放结果0X(0)4X(4)1Xr(1)5Xi(3)2Xr(2)6Xi(2)3Xr(3)7Xi(1)3.3 实序列FFT算法主程序调用上述FFT算法子程序实现实序列N点FFT运算的主程序代码参见附录2。下面对该程序代码做详细说明。主程序的流程比较简单，首先从文件in_data
34、.dat读取原始序列及其长度，然后调用子程序rfft计算该序列的DFT。计算得到DFT一方面直接在屏幕上显示出来，另一方面输出到文件out_data.dat。3.3.1 原始序列的产生和读取需要读取的原始序列首先通过在Matlab中执行seq_in.m文件产生，并存入到文件in_data.dat中。附录3给出了seq_in.m程序代码。代码中首先设置序列的长度N，然后以数组的形式直接设置输入序列x(n)，之后调用Matlab中的fopen、fprintf等函数将其写入到文件in_data.dat中。在调用FFT的主程序中，首先从in_data.dat文件中读取序列长度n。如果序列长度n大于DF
35、T的点数N，则令n=N，然后从文件中读取序列的前面n点，并存入数组x(n)，之后在屏幕上显示出原始序列。3.3.2 计算结果的显示和输出读取得到序列后，调用rfft函数计算其N点DFT。如果序列长度nN，则先通过补零将序列的长度扩展到N，然后再调用rfft函数计算FFT。表3.3.2-1为rfft函数返回计算结果的存放格式。对实序列的8点DFT，根据共轭对称性有：即因此，调用rfft函数计算得到的结果，必须经过整理。主程序中的如下语句：x0=x0;y0=0;for(i=1;iN;i+)if(i=N/2)xi=xi;yi=0;else if(iN/2) xi=xi;yi=xN-i;elsexi=
36、xN-i;yi=-yN-i; 就是实现这一功能的。整理后N点DFT的实部和虚部依次存入数组x和y中。其中X(0)和X(N/2)一定为实数，对应的虚部y(0)和y(N/2)存入0。之后根据iN/2，分别求得DFT各点的实部和虚部。整理后的计算结果首先直接在屏幕上显示，然后调用fprintf函数将其输出到文件out_data.dat。为便于在Matlab中读取，将DFT的点数N作为第一个数据也存入该文件中。3.4 运行结果分析程序运行后，为观察计算结果，可以有两种方法。第一种方法，直接在屏幕上显示；第二种方法，在Matlab中通过执行seq_out.m文件，读取计算得到N点DFT，然后利用Matl
37、ab提供的绘图函数画出N点DFT的波形。3.4.1 计算结果数据分析首先在Seq_in.m中设置实数序列x(n)= 1,1,1,1，设置序列的长度N=4。然后再MATLAB中运算该程序，将序列及其长度写入文件in_data.dat。然后启动VS2010，打开mrfft.c主程序，在程序开始的define语句处设置DFT的点数N=8，并在最后一个大括号处设置一个断点。运行程序后，屏幕上的显示如图3.4.1-1所示。图3.4.1-1 4点序列的8点FFT及IFFT在MATLAB中调用FFT函数计算得到的结果为4，1-2.4142i，0，1-0.4142i，0，1+0.4142i，0,1+2.414
38、2i将上述数据与图3.4.1-1比较可知，除了有一定的精度误差外，计算结果完全正确， 3.4.2 N点DFT波形分析在Seq_out.m文件中，有如下语句：for n=1:N z(n)=x(2*n)+i*x(2*n+1)endstem(0:N-1,abs(z),.)这些语句用于从Out_data.dat文件中读取计算结果，存入数组z中。MATLAB的数组中，每个数据可以是复数。然后，用stem语句和abs函数绘制序列DFT模的波形。在MATALB中执行Seq_out.m文件，读取计算结果，并绘制上述4点长门序列的8点DFT波形如图3.4.2-1所示。图3.4.2-1 4点长门序列的8点DFT波
39、形将seq_in.m程序中产生序列的相关语句修改为：N=64;x=cos(pi/4*n)+0.4*cos(3*pi/8*n)用以产生一个周期序列，该序列包含两个余弦序列分量，其周期分别为2p/(p/4)=8和16，因此整个序列的周期为16。设置序列的长度为64，刚好4个周期。修改主程序mrfft.c开始的define语句，将DFT点数N设为等于序列的周期。图3.4.2-1给出了该周期序列的64点DFT的波形。对上述结果分析如下。余弦序列的傅里叶变换为在w=02p范围内可以表示为由于序列的N点DFT是其FT在w=02p范围内的N点等间隔采样，则得到图3.4.2-1 周期序列的64点DFT对上述周
40、期序列中的两个余弦分量，其频率w0分别为p/4和3p/8，则由上式分别得到其64点DFT为再根据DFT的线性性得到原周期序列的64点DFT为这说明，该周期序列的64点DFT分别在k=8、12、52和56处有一个点，高度分别等于32、12.8、12.8和32。由图3.4.2-1可知，计算结果与理论分析结果完全一致。4 结束语 FFT是离散傅里叶（DFT）的快速算法，它利用旋转因子的周期性和对称性等性质，将DFT的计算逐次分解成较短序列的DFT，计算较短序列的DFT，再组合成原序列，从而减少了运算量，使运算更加高效。利用FFT频域特性，通过将输入的实序列转换为复数序列再进行FFT变换，提出一种FF
41、T有效处理实序列的算法并用C语言实现。本文重点介绍基2时域抽取FFT算法，介绍了DIT-FFT的基本原理以及其运算流图，并通过C语言实现以及MATLAB仿真，得出了计算结果与理论分析结果完全一致。本论文是在XX老师的悉心指导下完成的。向老师渊博的专业知识，严谨的治学态度，精益求精的工作作风，诲人不倦的高尚师德，、宽以待人的崇高风范，朴实无华、平易近人的对我影响深远。本论文从选题到完成，每一步都是在向老师的指导下完成的，倾注了大量的心血。从这次论文的完成我也学到了不少的知识，不仅使我树立了远大的学术目标、掌握了基本的，学会了独立自主的学习和查询参考文献，还使我明白了许多待人接物与为人处世的道理。
42、在此，谨向向老师表示崇高的敬意和衷心的感谢！附 录：附录1：实序列FFT算法子程序#include math.h#include stdio.hvoid rfft(x,n)int n;float x;int i,j,k,m,a1,a2,a3,a4,n1,n2,n4;float wnp,wnp0,wr,wi,xtmp,t1,t2;/倒序n1=n-1;for(j=0,i=0;in1;i+)if(ij)xtmp=xj;xj=xi;xi=xtmp;k=n/2;while(k(j+1)j=j-k;k=k/2;j=j+k;/求蝶形运算级数for(j=1,i=1;i16;i+)m=i;j=2*j;if(j=
43、n) break;/蝶形运算for(i=0;in;i+=2)xtmp=xi;xi=xtmp+xi+1;xi+1=xtmp-xi+1;n2=1;for(k=2;k=m;k+)n4=n2;n2=2*n4;n1=2*n2;wnp0=6./n1;for(i=0;in;i+=n1)xtmp=xi;xi=xtmp+xi+n2;xi+n2=xtmp-xi+n2;xi+n2+n4=-xi+n2+n4;wnp=wnp0;for(j=1;j=(n4-1);j+)a1=i+j; a2=i-j+n2; a3=i+j+n2; a4=i-j+n1; wr=cos(wnp);wi=sin(wnp);t1=wr*xa3+wi*xa4;t2=wi*xa3-wr*xa4;xa4=xa2-t2;xa3=-xa2-t2;xa2=xa1-t1;xa1=xa1+t1;wnp=wnp+wnp0;附录2：主程序#includemath.h#includerfft.c#includestdio.h#define N 64main()