实序列FFT算法的C语言实现.docx
- 文档编号:8410710
- 上传时间:2023-01-31
- 格式:DOCX
- 页数:36
- 大小:260.32KB
实序列FFT算法的C语言实现.docx
《实序列FFT算法的C语言实现.docx》由会员分享,可在线阅读,更多相关《实序列FFT算法的C语言实现.docx(36页珍藏版)》请在冰豆网上搜索。
实序列FFT算法的C语言实现
实序列FFT算法的C语言实现
学生:
XX指导教师:
XX
内容摘要:
DFT和IDFT是数字信号分析与处理中的一种重要运算和变换,但直接根据定义计算DFT时,运算量大,不能实时得到计算结果。
特别是在实际应用中,N都取得比较大,此时,由于乘法和加法的次数都近似与N的平方成正比,这将极大增加DFT计算所需时间。
为此,提出了许多DFT和IDFT的快速算法,称为快速傅里叶变换(FFT)和快速傅里叶反变换(IFFT)。
本文较为系统地阐述了快速傅里叶变换的算法原理然后用MATLAB实现了快速傅里叶变换。
论文首先首要介绍了FT与DFT的定义、DFT与FFT的关系,然后重点介绍基2时域抽取FFT算法以及其原理和运算流图,再应用C语言实现了实序列的FFT。
最后在Matlab软件上进行仿真,仿真结果验证了设计的正确性。
关键词:
傅里叶变换快速傅立叶变换Matlab仿真
RealizationofFFTalgorithmforrealsequence
WithCprogram
Abstract:
DFTandIDFTareimportanttransformandprocessingindigitalsignalprocessing.However,therearelargeamountofcomputationbydirectlycalculatingaccordingtothedefinitionofDFT.Especiallyinthepracticalapplication,Nisbigger,atthistime,becausethetimeofmultiplicationandadditionareapproximatelyproportionaltothesquareofN,whichwillgreatlyincreasethecalculationtimeneededforDFT.Therefore,manyDFTandIDFTfastalgorithmareraised,whicharecalledFFTandIFFT.
InthispaperrelativelysystematicallyelaboratedthefastFouriertransformalgorithmprincipleanduseMATLABsoftwaretorealizethefastFouriertransform.ThepaperfirstintroducesthedefinitionofFTandDFT,therelationshipbetweenDFTandFFT,andthenmainlyintroducesDIT-FFT,includingitsprincipleandoperationflowdiagram,andfinallyusedClanguagetorealizetherealsequenceFFT.ThedesignsaresimulatedinMatlabsoftware,theresultsofthesimulationconfirmtheexactnessofthedesign.
Keywords:
FouriertransformationfastFouriertransformationMatlabsimulation
实序列FFT算法的C语言实现
前言
在实际的数字系统中,DFT是一种得到了广泛的应用的、重要的信号处理手段,但它的运算效率非常低。
随着DFT输入的点数增加到数百或数千,DFT需要的运算量变得非常大。
快速傅里叶变换(FFT)可使实现DFT的运算量下降几个数量级,从而使数字信号处理的速度大大提高[1]。
FFT是离散傅立叶变换(DFT)的快速算法,可以将一个信号变换到频域。
有些信号在时域上是不易看出有什么特征的,但是如果变换到频域之后,就很容易了。
FFT可以将一个信号的频谱提取出来,这在频谱分析方面也是经常用的。
实际中需要做快速傅里叶变换的多为实序列数据,而其变换算法都是以复数序列作为输入。
本设计利用频域的性质,将实序列数据变为复数序列,再进行FFT,以提高效率。
本设计就是要求在熟悉DFT及FFT算法基本原理的基础上,编制C语言程序实现实数序列的FFT算法。
1序列的FT和DFT
序列的傅里叶变换(FourierTransform,FT)和离散傅里叶变换(DFT)都是对序列的频域描述,它们揭示了序列由那些分量构成,各分量的幅度和相位大小。
1.1序列的FT
序列x(n)的傅里叶变换又称为离散时间傅里叶变换(DTFT),其定义为
(1.1-1)
式中,称为数字角频率。
如果已知x(n)的傅里叶变换X(ej),则可下式求得其时域表达式
(1.1-2)
上式称为序列的傅里叶反变换(IFT)。
序列的傅里叶变换是的连续函数,且一般情况下为复变函数,并可表示为
其中,
和
分别称为幅度谱和相位谱。
此外,根据
的周期性可知,序列的频谱及其幅度谱和相位谱都是关于以2为周期的周期函数。
1.2序列的DFT
序列的FT变换
为的连续函数,不便于用计算机程序进行辅助分析和计算。
为了便于用计算机辅助求解和分析,提出了离散傅里叶变换(DFT)。
1.2.1DFT的定义和计算
长度为M的有限长序列x(n),n=0~M-1,其N点DFT和IDFT(离散傅里叶反变换)分别定义为
,k=0~N-1(1.2.1-1)
,n=0~N-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=0~N/2-1的点,即可根据共轭对称性,由式(1.2.1-5)求出另外N/2个点,即k=N/2~N-1的点,合起来得到完整的N点DFT。
在式(1.2.1-5)中分别令k=0和N/2,又可得到
以上两式说明,在实序列的N点DFT中,当k=0和N/2时一定为实数。
2FFT算法
DFT和IDFT是数字信号分析与处理中的一种重要运算和变换,但直接根据定义计算DFT时,运算量大,不能实时得到计算结果。
特别是在实际应用中,N都取得比较大,此时,由于乘法和加法的次数都近似与N的平方成正比,这将极大增加DFT计算所需时间。
为此,提出了许多DFT和IDFT的快速算法,称为快速傅里叶变换(FFT)和快速傅里叶反变换(IFFT)[2]。
FFT和IFFT算法的基本思想是根据旋转因子WNm的对称性,将N点DFT分解几个较短的DFT,从而减少乘法次数,缩短计算时间。
根据这一思想,出现了许多DFT和IDFT的快速算法。
下面结合本设计重点介绍基2时域抽取FFT算法。
2.1基2时域抽取FFT算法
时域抽取FFT算法(DIT-FFT)的基本思想是在时域中将N点序列x(n)中的各点按序号的奇偶进行抽取分组,得到两个N/2点序列。
在计算出两个N/2点序列的DFT后,通过蝶形运算得到原序列的N点DFT[3]。
2.1.1基本原理
假设序列x(n)的长度为N=2M,其中M为正整数。
根据序号n的奇偶将x(n)分解为两个N/2点的子序列,即
则
令
分别为序列x1(r)和x2(r)的N/2点DFT,则得到
,k=0~N/2-1(2.1.1-1)
由于X1(k)和X2(k)都以N/2点为周期,且
因此式(2.1.1-1)可表示为
,k=0~N/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)。
图2.1.1-1蝶形运算符号
重复上述分解过程,继续将序列x1(n)和x2(n)按序号的奇偶进行抽取分解,直到最终分解为N/2个两点序列。
通过每次抽取分解,都将序列和DFT的点数缩短一半,同时转换为若干蝶形运算。
2.1.2DIT-FFT算法的运算流图
根据上述分解过程,容易知道,对N点DFT,只需通过M=log2N次抽取,就可以将原来的N点DFT分解为N个1点DFT和M级蝶形运算,而1点DFT就是序列本身。
以
,从而得到完整的蝶形运算流图如图2.1.2-1所示。
图2.1.2-18点DIT-FFT运算流图
2.1.3DIT-FFT算法的运算量和存储量
根据上述分解过程,容易知道,对N点DFT,只需通过M=log2N次抽取,就可以将根据上述DIT-FFT的基本原理,通过时域抽取,将N点DFT转换为M级蝶形运算,每级包括N/2个蝶形运算,因此N点FFT共包括MN/2个蝶形运算。
每个蝶形运算需要一次乘法运算和两次加法运算,则N点FFT的运算量为:
乘法:
加法:
而DFT直接算法共需N2次乘法和N(N-1)次加法。
当N>>1时,采用FFT算法将大大减少运算次数,提高运算速度。
由图2.1.2-1所示运算流图可知,DIT-FFT的运算过程很有规律。
在同一级的N/2个蝶形运算中,每个蝶形运算的两个输入输出数据只对该蝶形有用,而且都位于同一条水平线上。
因此,每个蝶形的两个输出数据可以只占用输入数据的存储单位,而不要另外分配存储单元。
这种存储计算数据的方法称为原位计算。
原位计算可以极大地节省存储单元。
对N点FFT算法,只需要N个存储单元。
运算一开始,这N个单元用于存储原始序列。
计算过程中,用于存储各蝶形运算的输出数据。
运算完毕后,这N个单元中存储的就是输出N点FFT[4]。
2.2实序列的FFT算法
实际系统中处理的序列都是实数序列。
如果直接按照前述的FFT算法进行计算,就是将其视为虚部都为0的复数序列,这就增加了存储量,降低了运算效率。
处理该问题的方法有很多。
例如,通过一个N点FFT同时计算出两个实数序列的N点FFT;根据基2时域抽取FFT算法的基本思想,用一个N/2点FFT计算出一个实序列的N点FFT[5]。
本设计采用另外一种方法,以减少运算量和存储量。
已经知道,实序列的mDFT都具有共轭对称性,或者说,其实部都是偶对称的、虚部都是奇对称的。
利用这种共轭对称性,计算时可以只计算出序列的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)一定是实数,因此计算过程中只需为X(0)和X(4)分别分配一个存储单元。
3实序列FFT算法的C语言实现
根据设计要求,利用VS2010软件编制了C语言程序,实现实序列的基2时域抽取FFT算法。
3.1VS2010简介
VS2010提供了两类容器,以便有效地管理开发工作所需的项,如引用、数据连接、文件夹和文件。
这两类容器分别叫做解决方案和项目。
此外,VS2010还提供解决方案文件夹,用于将相关的项目组织成项目组,然后对这些项目组执行操作。
作为查看和管理这些容器及其关联项的界面。
为了使集成开发环境(IDE)能够应用它的各种工具、设计器、模板和设置,VisualStudio2010(简称VS2010)实现了概念上的容器(称为解决方案和项目)。
另外,VisualStudio还提供了解决方案文件夹,用于将相关的项目组织成组,然后对这些项目组执行操作。
解决方案管理VisualStudio配置、生成和部署相关项目集的方式。
解决方案可以只包含一个项目,也可以包含由开发团队联合生成的多个项目。
复杂的应用程序可能需要多个解决方案。
项目包含一组源文件以及相关的元数据,如组件参考和生成说明。
生成项目时通常会生成一个或多个输出文件。
根据实际需要,项目可以简单,也可以复杂。
一个简单的项目可能由一个窗体或HTML文档、源代码文件和一个项目文件组成。
更加复杂的项目可能由这些项以及数据库脚本、存储过程和对现有XMLWebservices的引用组成。
创建新项目时,VS2010会自动生成一个解决方案。
然后,可以根据需要将其他项目添加到该解决方案中。
也可以创建不包含项目的空白解决方案,从而使用VS2010编辑器和设计器修改独立的文件。
因为每个项目或解决方案都由一个目录及其内容组成,所以,可以在Windows资源管理器中移动、复制或删除解决方案和项目。
VS2010将解决方案的定义存储在两个文件中,即.sln和.suo。
这两个文件相当于早期版本中的工作区文件(.dsw)。
解决方案定义文件(.sln)存储定义解决方案的元数据,每当解决方案活动时,都使用构建该解决方案并设置其属性时存储在.suo文件中的元数据来自定义IDE。
3.1.1新建项目
为了调试本设计所要求的程序,首先启动VS2010,然后依次选中文件菜单中的新建、项目命令,打开如图3.1.1-1所示的新建项目对话框。
图3.1.1-1新建项目对话框
在该对话框中,选中左边VisualC++下面的Win32,在右边选中Win32控制台应用程序。
然后在名称框中输入项目名称。
这里设定项目名称为rfft,并设定项目存放的位置为E盘。
注意,一旦设定项目名称后,解决方案的名称默认与项目名称相同。
单击确定按钮,在后面出现的个对话框中都选默认值,在应用程序设置对话框中确保选中“空项目”,之后,即可创建项目rfft,同时创建一个同名的解决方案。
现在打开E盘,可以看到一个rfft文件夹,其中有一个rfft.sln文件,代表解决方案,同时还有一个rfft子文件,其中有一个文件rfft.vcxproj,代表rfft项目。
3.1.2新建文件
为了编辑源程序代码,在VS2010的文件菜单中选中新建、文件命令,在弹出的新建文件对话框中选中VisualC++、C++文件。
单击打开按钮后,可以看到在VS2010的工作窗口中出现一个新文件“源1.cpp”。
选中菜单中的另存为命令,将其重新保存为C语言源程序文件,这里命名为rfft.c,用于新建实序列FFT算法子程序。
注意在另存为对话框中一定要选中C源程序选项,并找到文件夹rfft,将文件保存到该文件中。
重复上述步骤,再新建一个C语言源程序文件,命名为mrfft.c,该文件是主程序文件,用于调用rfft.c文件中定义的rfft函数,计算实序列的N点DFT。
分别在mrfft.c和rfft.c文件中编辑实序列FFT算法的子程序和主程序,之后,为了将文件添加到项目中,在VS2010工作窗口左边的解决方案资源管理器中右键单击项目名称rfft,在弹出的菜单中依次选中添加、现有项命令,找到刚才创建的文件mrfft.c,单击添加按钮,将其添加到项目和解决方案中。
将主程序文件添加到解决方案后,通过生成菜单中的解决方案命令,对程序进行编译连接。
如果没有错误,则生成解决方案,并在rfft文件夹下自动生成一个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点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顺序和倒序对照表
顺序
倒序
十进制
二进制
十进制
二进制
0
000
0
000
1
001
4
100
2
010
2
010
3
011
6
110
4
100
1
001
5
101
5
101
6
110
3
011
7
111
7
111
由表3.2.1-1可知,自然顺序序号加1,是在其二进制数最低位加1,并逢2向高位进1。
而倒序序号是在其对应的二进制数最高位加1,逢2向低位进位。
对于N=2M,M位二进制数各位的权值从左向右依次为N/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所示。
根据上述框图编制的实现倒序的程序代码如下:
图3.2.1-1倒序程序框图
n1=n-1;
for(j=0,i=0;i {if(i {xtmp=x[j]; x[j]=x[i]; x[i]=xtmp; } k=n/2; while(k<(j+1)) {j=j-k; k=k/2; } j=j+k; } 代码中,i和j分别为自然序号和倒序序号,设置i从0到N-1,共循环N次。 在每次循环中,当i 代码中之后的while循环语句用于实现倒序值j的计算,以便下次循环实现相应一对数据的交换。 3.2.2蝶形运算 将输入序列倒序后,进行各级蝶形运算。 程序中首先用如下语句: for(j=1,i=1;i<16;i++) {m=i; j=2*j; if(j==n)break; } 求得N点FFT蝶形运算的级数。 之后,进行第一级蝶形运算,再用k控制共循环m-1次,依次进行后面各级蝶形运算。 由FFT算法的运算流图可知,对第一级中的每个蝶形运算,旋转因子都为 =1。 再考虑到输入序列各点都是实数,因此将序列中相邻每两输入数据直接相加减,即得到第一级各蝶形的输出。 相应的代码如下: for(i=0;i {xtmp=x[i]; x[i]=xtmp+x[i+1]; x[i+1]=xtmp-x[i+1]; } 在后面各级蝶形运算中,由运算流图可知,有些蝶形所需的旋转因子仍然为 =1。 对这些旋转因子对应的蝶形运算,同样化简为将其两个输入数据直接相加减,以得到该蝶形的两个输出。 实现这一运算的主要代码如下: xtmp=x[i]; x[i]=xtmp+x[i+n2]; x[i+n2]=xtmp-x[i+n2]; 其中,n2=2k-1,k=2~m。 例如,对第二级蝶形,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)-x(6) 对第三级,n1=8,因此上述语句只执行1次,实现的运算是: x(0)=x(0)+x(4);x(4)=x(0)-x(4) 程序中还有如下语句: x[i+n2+n4]=-x[i+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点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*x[a3]+wi*x[a4]; t2=wi*x[a3]-wr*x[a4]; 执行时,a3=5,a4=7;wr和wi分别为 的实部和虚部的负,即wr=C1,wi=-S1。 因此,如果令x[a4]=-x(7),则t1和t2就分别等于式(3.2.2-3)中两个中括号中的项。 这也就解释了为什么在执行第二级运算时需要将x(7)取负。 程序中后面有如下两条语句: x[a4]=x[a2]-t2; x[a1]=x[a1]+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)取负,再由这两条语句计算,得到的输出刚好为最后N点FFT结果中X (1)的实部和虚部,并且X (1)的实部和虚部分别存放在序号为1和7的单元。 根据同样的分析,程序中另外两条语句 x[a3]=-x[a2]-t2; x[a2]=x[a1]-t1; 用于计算最后结果中的X(3),并且将其实部和虚部分别存放在序
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 序列 FFT 算法 语言 实现