FFT算法DSP实现(DIT基2FFT运算流图)

1、0J国科技疋乂任线hnp /www paper edu cnFFT算法的研究与DSP实现王靖琰=王靖斌$|中国科学院上海应用物理研究所，上海(201800)'山东农业大学信息工程与技术学院，泰安(271018)E-mail: ivcamail ustc edu cn0J国科技疋乂任线hnp /www paper edu cn摘本文介绍了 DIT基2FFT、DIF基2PFT,分裂基FFT和混合基FFT算法的原理,并在TMS320C6713B浮点型DSP上实现了算法，最后对各种算法进行了对比.关键词：时域抽取基2FFT,频域抽取基2FFT,分裂基FFT,混合基FFT, TMS320C671

2、3B1-引言傅里叶变换是一种将信号从时域变换到频域的变换形式。是声学、图像、电估和信号处 理等领域中一种垂要的分析工具叭其中,离散傅里叶变换(DFT)更是数字信号处理领域不 町缺少的工几么一，特别是快速傅里叶变换(FFT丿的出现使得DFT在实际应用屮得到了广泛 的应用。FFT算法就是不断地把氏序列的DFT逐次分解为几个短序列的DFT,并利用的周期 性和対称性來减少DFT的运算次数这一基本原理來加快运算速度的。-般高基数FFT算法耍比低棊数FFT算法运算杲耍小，计算速度要快，但需要貝冇更强 大的数据并行处理能力的器件支持且编程难度也増犬了。本文分别在浮点熨 DSPTMS320C67I3B I：实

3、现了DIT、DIF某2FFT算法，分裂基FFT算法和混合某FFT算法及其 逆变换，并对语音信号进行变换和逆变换的实验，対比分析J'各种算法的时间复朵度与准确 度。2.各种FFT算法原理介绍与比较现冇的FFT算法包括针对N=2L点的基2FFT、甚4FFT、分裂旱FFT以及针对N不等J'2L 的混合基FFT。其中毎种算法乂可以分为时域抽取算法(DIT)和频域抽取算法(DIF)o F 面分别介绍它们各H的特点并对它们进行比较，导出各H的使用场合.2.1 DIT 和 IDF 基 2 FFT 算法对F N=2I-点DFT可按照如卜方法做时间抽取：JJ丨JX (灯 £ x(2r)

4、W + X x(2r +1)比严从=士 x(2r)W + 必士如 1)町(1)r=0r=0r»0Vr=0VNN这样就把一个N点DFT转换为个基本碟形运算和2个一点DFT运算，如时域抽取22® 2FFT 一级分解框图所示：x(2r)x(2i+l)NX(k+号)图L时域抽収基2FFT级分解松图对框图屮每个N/2点DFT继续进行同样的分解，直到只仃最某本的2点DFT为止.就 得到了时域抽取基2FFTo由时域抽取，造成了输入数据的二进制倒序。N Nk T如果对J： N=2I-点DFT按照如卜方法做频域抽取：=£ v(/) + Wjx(n + ) WtJAn竺，X 伙)=&

5、#163; x(n)V；k + £ x(n += _n=0n=0Z/ts0fdjW/T = (_/ ,分别另 k=2i和 2r+l, J是得到x)£H-0N x(n)-x(n + -) Wm)+m+号)吩X(2“l)=£/i"0nkN_0J国科技疋乂任线hnp /www paper edu cn0J国科技疋乂任线hnp /www paper edu cn这样就把-个N点DFT转换为号个基本碟形运算和2个号点DBT运算，如频域抽取某2FFT 级分解框图所示:x(n)x(n+y)X(2i)X(2i+1)图2频域抽取基2FFT 级分解松图对框图屮每个N/2点D

6、FT继续进行同样的分解.直到只白最某本的2点DFT为止，就得到 了频域抽取基2FFT。由J：频域抽取，造成了输出数据的二进制倒序。以上两种算法就是Coolev和Tukev在1965年提出的FFT算法，这两种算法的计算彊是一N. Mlog；, Nmi02 n样的，部使DFT的乘法计算由N2次降为2次，加法次数由N(N-1)次降为皿°禺" ID国科技疋乂任线hnp /wvw paper edu cn次。DIT与DIF的实质区别在J：先乘旋转肉子再进行碟形运算(DIT丿还是先进彳j：碟形运算再 乘旋转因子(DIF) 3o后來又提出的基4算法、分裂墓算法等都是针对N等2的整数次邪的

7、 算法，他们都右DIT和DIFZ分，但是其本质区别与基2FFT足一样的。2.2分裂基FFT算法1984年提出的分裂啟spgradix丿算法同时使用康2和取4算法，被认为是II前对J'N-2L 齐类算法中最为理想的一种。j'N=2L点、DFT,取写式的DIF的偶序号输出项，即叫X(2r)=£/i=0x()+W +号)W,对k的奇序号项用基4算法,即X(4r+l)=X/=0X(4r+3)=士n=0(g) “ + 如)-7( + %) 5 + >比咋(x(“) - x(n + 斗)+ j(x(n + y) - (x(n +24匕面三式构成了分裂基算法的L形算法结构,如

8、图3所示:兀)4x(n+y)x(n+y)4X(2i)X(4i 十 1)X(4i十 3)图3频域抽取分裂丛FFT 级分解椎图対框图屮每个N/2点和N/4点DFT继续进行同样的分解.宜到只有最基木的2点DFT 和4点DFT为止，就得到了频域抽取分裂基FFT。由频域抽取，造成了输出数据的二进 制倒序。事实上分裂基的到序算法与垄二的是样的，可以使用同样的程序。分裂基的每个L形仃两个复数乘法运算，所以N点分裂丛FFT乘法个数是L形个数的-#-ID国科技疋乂任线hnp /wvw paper edu cn2倍，而加I法个数与某2FFT相同叫 第j级中的L形个数匚与第J-1级的关系如卜式所示:厶=普_号八1&

9、#39;2，厶_1全部的复数乘法运算次数为-X, N'j)3.vz z ii2?cm = y厶亍+丫一歹厶=3川1昭2川_§+(-1)厶6(6)(7)-5-ID国科技疋乂任线hnp /wvw paper edu cn-#-ID国科技疋乂任线hnp /wvw paper edu cn与基2FFT算法的复数乘法相比，N呃 N前而的系数由+变到+，仅这-项就可以使 复数乘法次数卜降33%。2.3混合基FFT算法上而讨论的是对J； N-2l的FFT,而若不满足N-2l,可以将號补一些零点使N増长 到最临近的一个21数值在使用上面的算法；但如果要求准确的N点DFT,而N又是素数， 只能

10、采用直接DFT算法或CZT方法：若N是一个复合数，即它可以分解成一些因子的乘积, 则可以使用混合基FFT算法現设DFT点数N =血是一个复介数，将n用卜面公式表达:-#-ID国科技疋乂任线hnp /wvw paper edu cn-#-ID国科技疋乂任线hnp /wvw paper edu cn这样可以把泯®分为乙个/；点的序列。令N =叨,将频率变帚k用卜式衷达:则DFT公式可以化为:X(k、,ko)=、:/io=OWj=O(9)(10)j是我们得到了 n =応点、的混介甚fft r/；z： n先将d表示为按乙进制表示的形式=也+ 0，即将分为乙个斤点的序列，然后做个斤点的DFT得

11、到X伙o，o)：将 X(A0,/?0 )乘以旋转因子呛人后，再做片个乙点的DFT得到兀伙。虫)，最后进行乙进制 整序,得到X(kl9k0)=X伙)。N =仁点混介基算法的框图如图 4 所示。-#-0J国科技疋乂任线hnp /wvw paper edu cn图4频域抽取混介荃FFT 级分解框图这是N = W点、频域抽収混合基FFT -级分解的算法，将£和人进彳J：进一步的分解，直 到所有因子都是素数就得到了混合基的FFT算法。可以证明，当N = /匚时，采用混合某算法所需要的总乘法次数为：3 "&/；)+1 (11)2.4各种FFT算法的比较1965年Coolv与Tu

12、kev捉出基2FFT算法以后，新的算法不断涌现，主耍仃两个方向： 一是针对N等于2的整数次幕的算法,如分裂基算法,另一种是针对N不等于2的整数次 耶的算法，如混介基础算法。分裂卑算法的优点是时间复杂度小，速度快，但是对J；N不等2的整数次幕的FFT 只能采取补零的方法，如果要求精确的N点FFT则无能为力，灵活性较差：而混介甚FFT 在保证时间复杂度较小的条件卜，可以对任意可以分裂为素数因子的N点信号进行FFT, 灵活性好。将这些算法的特点具列一表，进行对比如卜：衣1并种FFT榇法的特点DFT基 2FFT分裂基FFT混合基FFT时间复杂度差良优良灵活性好差好3. FFT算法在DSP上的实现本文算

13、法在DSP上的实现車点是程序开发，由I：人最运用裔级语言，减少了程序开发 的周期。开发匸作主要完成在TMS320C6713B上的频率抽取卑4FFT程序及其逆变换程序, 以及对语音的FFT和IFFT变换与算法的评价。3.1 TMS320C6713B 浮点 DSP为了满足算法的高将度以及信匕的较丿、动态范H；l的耍求，则需耍浮点DSP來实现本算 法。这里选取了 TI浮点DSP中性价比较高的C6713B作为整个系统核心。TMS320C6713B是一款高性能的32位浮点DSP,适用专业音频信匕处理。该芯片的 内部结构是在TMS320C67XX的茲础上加以改进制成的，其可工作在225MHz的主频上， 拥

14、仃8个独立功能单元，并釆用了先进的VLIW体系结构及流水线技术；独立L1/L2存储 器结构，多种外设接I I方便扩展；它采用3.3V的I/O电斥和1.8V的内核电床供电方式，并 貝有两级cache缓存结构。除此之外，它还有行速度快和蒂度高两个主要特点，完全能够满 足算法的耍求。本文采用的TDS6713EVM丿I：发板是闻亭公司仏新研制的高速语音信号（采集）处理平 台，可作为专用语抒信号编解码处理测试平台，也町用J：各种对数据桔度令特殊要求的浮点 数字信兮处理场合。TDS6713EVM系统结构图如廉5所示。IVI 5 TDS67L3EVM系统结构图4.2DIF 与 DIT 基 2FFT 的 DS

15、P 实现为了实现DIT与DIF基2FFT,肓先观察它们的的运算流图，进行对比研究，得出一般 规律。8点DIT和DIF的基2FFT运算流图如图6所示。通过对比观察以卜两图町以发现，他们都是分级运算的，共L = log2A/级。每级仃 N/2个某本碟形，而H是分组的。DIT算法的第一级有一个组，往后逐级二分宜到有N/2个 组，毎组所乘的旋转因子是同样的，并I!毎组的的毎个蝶形输入数据Z间的间隔是逐级减小 到2； DIF算法的第一级有N/2个组，往后逐级：分宜到只有一个组，毎组所乘的旋转因子 也是同样的，并II.毎组的的每个蝶形输入数据Z间的幼隔是逐级増加到N/2的。DIT的输入 数据是按二进制倒序

16、排列的，而DIF的输出数据是按照二进制倒序排列的。ID国科技疋乂任线hup /www paper edu cnx(0)M4)x(2)x(6)Ml)x(5)M3)x(7)-1X(0)X(l)XXG)XX(5)X(6)X(7)x(0)x(l)x(2)x(3)M4)x(5)x(6)x(7)(a)8点DIT基2FFT流图X(0)XXX(6)X(l)X(5)XX(7)(b丿8点DIF基2FFT流图图6 8点DIT和DIF的基2FFT运知流图根据以1：总结出的规律，我们发现编写DIT和DIF基2FFT程序的基本步骤相同，即先 分级，然后按照每级中的组进行碟形预算和旋转因子的相乘，不同Z处在与每级中的组中蝶

17、 形运算与旋转因子相乘运算的控制，以及倒序在运算前还是运算后。这里给出DIT卑2FFT 的流程图：定义DIT某2FFT算法程序如卜：void fftdit(double *x,double *yunt n),泄义DIF基2FFT算法程序如卜：void fftdif(double *x,double11),对/ DIT算法程序中的运算级控制采用以卜循坏：for(l= 1 J<=m,l+)而DIF则采用：for(L=m,L>=IX)对于该级中每组的元素数,DIT采用nl-2*nl的方式进行计算,DIF则采用nl-nl/2. 这些编程上的不同之处都显示了 DIT和DIF各自的特点。采用C

18、卄语言编制以上程序，并在CCS上编译卜锻到DSP开发板上就可以实现甚2FFT 的DSP实现了。ID国科技疋乂任线http /www paper edu cn4.2分裂基FFT算法的DSP实现首先我们给出16点DIF分裂某FFT运算流图如图8所示。AAIA-4图8 16点DIF分裂堆FFT运舁流图由此图可以看出，分裂基FFT算法结构同基2FFT算法结构相似，适用J' N=2L的场合, 并右L级运算实现。运算流图的输入顺序，输出为二进制倒序。当然，由J：L形的特殊结 构，便秤序实现比基2PFT或基4FFT算法稍复杂起叭整个程序需耍需L-1次分级运算，其程序流程图如图9所示。定义DIF分裂基

19、FFT函数如卜：void srfftdouble *x,doublen),这就实现了分裂某FFT算法，用C卄编制程序，在CCS环境卜编译卜载到TMS320C6713 中就可以了。图9 DIF分裂准FFT程序流程图4.3混合基FFT算法的DSP实现实现混合基FFTfM序较为复杂，首先耍对运算长度N进行因子分解.分解为N=巨-#-0J国科技疋乂任线hnp /wvw paper edu cn-#-0J国科技疋乂任线hnp /wvw paper edu cn1-3rL：然后分为L级进彳j：蝶形运算，毎级蝶形运算分为塁个£,i = 0,l,厶-1点DFT运算；対J： DIT算法每级蝶形前耍与旋

20、转因子相乘.而对J：DIF算法每级蝶形后耍与旋转 因子柑乘；对J'' DIT算法输入数据要先进行混m寸J-DIF算法最后要对输出数据进彳J： 混序。图10 DIT混合里FFT程庁流程图其体编制程序时碰到的难点首先是对N进行分解时，如何确定分解的内子。针对这一 问题，笔者认为要遵循以卜原则：（1）分解出來的£丿=0,1,上-1耍尽最有萬效的石点DFT运算：（2）1；,/= 0,1,的积要等PH,但己经编好高效£点DFT运算程序的人之积与 N不一定相等，但耍使差不人，以使剩余的未编制的DFT级数尽帚少。（3）要使级数L尽起少，以使混序的运算复杂度尽杲小。根据以上

21、原则，结合输入点数的规律，笔音在程序屮首先尝试把N顺序分解为2, 3, 4.5, 8, 10之积，然后把N与这些因子之商作为最后一个因子，可以证明这个商为整数。需LL- 耍特别说明的是，対J：16点的DFT,分解为4x4点DFT的效率j;2x8点DFT.所以因 子分解完毕后，如果最后有因子2与8相邻,则用两个4替换。程序中实现如下:if (factorsj = 2)L = J，while (U>0丿 && (fhctoi:sL != 8)if (i>0)factorsfj = 4,factorsfi 4,定义因子分解的函数如K：void factonze(mt n,

22、 int *nFact, inr fact),丙子分解Z后.为了控制毎级旋转I大I子乘法和蝶形运算，需要产牛:两组控制因子:sofar(i)= 口 actual® j= 12 丄冃ienjam(i)= , j= 12 丄IIactual(j)此中，actual是第1级的棊因子，sofar是第1级前己经运算完毕的因子的乘积.而reimin 为第还未运算的因子的乘积。这一功能通过以卜换数实现：void transTableSetup(int sofar, mt actuaL int rerminf,int *nFact, int *nPoints),由本算法是DIT算法，需要对输入信号进

23、行混序，以便让卜面的变换可以在原址进 行，绘终的结果是按顺序排列的。对输入信号进行混序也是按照分级的思想进行的，即先按 照进行混序，然后在対每组E 点的数据进行按照7)，依次类推直到最后一级。混序通过以下函数实现：int fact, int remain, double xRe, double xlin, double yRe, double vlm),然后就是进行每级的旋转因子乘法和蝶形运算了。每级的蝶形和旋转因子都不同，而本 程序是通过twiddleTransf 个程序实现的，所以需要传入每级的控制参数，也即是 transTableSetup两数的运算结果.首先通过lintTrig怖数初试

24、化旋转因子，然肩按照每级的 规律进行乘法。对每级，进行蝶形运算的基不同，通过如卜程序实现：switch(radix)case 2: /l-FFT/break,case 3: /*3-FFT/break,case 4: /*4-FFT/break,case 5:/*5-FFT/break,case 8:break,case 10:break,default:fft_odd(mdLX力 break,值得注意的是，首先尝试分解的6个因子不是独立无关的，比如2与4 Z枳为8, 2与 5 Z积为10,事实上，在程序屮8和10点DFT是通过两个4和5点DFT实现的，这样的 效率高丁直接的8和10点DFTo

25、定义每级的旋转因子乘法和蝶形运算甫数如卜：void twiddleTransfinr sofarRadix, int radix, int remainRadix,double yRe, double ylm力对每级调用此函数即可以实现混合基FFT T：for (count-1, count<-nFactor, count+) nMddleTransf(sofarRadLxcount, actualRadixcount, remainRad ix count.vRe, ylin),这样就实现了混介基FFT篦法，定义混介廉FFT算法函数如2void mLxfftint n, double

26、xRe, doubledouble yRe, double ylm人4.4 DSP上FFT算法程序的优化结合DSP处理器与FFT算法的特点，综合考虑编程复杂性，本文对FFT算法程序进彳J： 了优化，使其效率有所提高。对J基2FFT和混介基FFT,旨先都需耍根据输入的点数N求解级数L=log2N.垠直接 的求法圧调用标准函数log(J函数，通过Lrg*og(N/log求L,但是调用nnth.h中川明的 函数是比较低效率的。由J：N是2的耶，所以可以采用以下方法高效实现求L：for(L= 1 ,(N=N/2'=1 丄+),将以匕两中方法分别在TMS320C6713E匕实现，采用CCS中的P

27、rofiler的程序断点时 钟统计功能,发现调用log。函数求L需要975个时钟，而本文的方法仅需要198个时钟, 可见其高效性。正如上面所说，调用imth.h屮声明的标准库函数是比较低效率，我们可以像上面那样 用比较高效的方法替代它來实现同样的功能來提高效率，但是仃些场介不得不调用它们，这 时候我们可以使用TI公司为用户提供的TMS320C67xFastRTS Ubrary屮的数学函数來捉高 效率。TMS320C67X FastRTS Libniry中的隨数都是TI公司用汇编语言专门针对处理器的特 点高效实现的，下面我通过cos运算來说明。如果我们在程序中包<t<math.h&g

28、t;头文件，使用标准函数库的double cos ( double x ),进 行一次运算需要211个时钟：如果将fastns67x.h头文件包含进来,使用其中声明,在 fastrts67x.lib库文件屮定义的double cosdp( doublex丿函数，只需耍158个时钟，口J见使用库 的优势。爭实上，TI公司为一般的数学运算、数字信号处理和图像处理算法都提供了非常 完善的库，我们在开发屮应该尽最使用以提高开发效率和程序的质駁。-15-ID国科技疋乂任线hnp /wvw paper edu cn5.实验结果以及性能分析我们已经讨论了某2FFT的DIT和DIF算法、分裂皐算法和混介某算法

29、各自的特点并 在DSP上实现了各算法，卜而我们通过DSP上的语音信兮的FFT实验來貝体地对比并种算 法。忖先我们对信勺长度为1024点的语音宿号进行FFT算法和逆变换，对比其时间复杂度 和准确性。从WAV文件中读取1024点语音数据进行变换。1024点语音信号波形如图11所 示。图II L024点语音信号波形利用DIT基2FFT程序对其进行FFT变换，得到的频谱的实、虚部如图12所示:频谱实部波形(b丿频谱虚部波形图13 M.ATLAB1024点语音信弓FFT频谱波形山图町以看出，本程序实现的FFT算法与MATLAB屮的FFT算法输出儿乎一样，这 就验证了本程序的正确性。同时我们也可以看出，i

30、in输入信号是实信兮，造成了频谱实部 NN关J：对称，频谱虚部关J：反対称，这也从-个方而验证了本程序FFT的正确性国。22利用DTT与FFT的关系实现IFFT程序，对频谱信号进行反变换得到语音信号的原始 波形如图14：图141024点语音信号频tniFFT得到的信号波形可以看到，IFFT后得到的信号波形基本与原始信兮一致，从虚部波形上町以看到，何 一些数磺级极小的谋差存在。我们借用通信原理屮的信噪比的概念，并对针对程序的实现对 其他进改造，用以啟化误差的人小佗 设为原始信兮，Y0为FFT后再IFFT得到的倍 号,我们定义SNIV如下:-L5-0J国科技疋乂任线hnp /wvw paper e

31、du cnSNR'=101og| £|y(")-x(/?)|：(13)Ln-0_托单位为分贝(dB丿。対以上1024点信号求SNR245.3592667583791 (dB)。现在我们对不同点数的信号用同样的分析手段进行对比实验。对2点的信号的FFT 貫接调用惭数：対J：非21点的信号，町以宜接调用混合基FFT,具他算法要先补冬到产点, 再调用FFTWi数。实验结果如表2、3所示。表2对不同点数各种FFT知法的运算虽(运算时钟数为单位)点数算法DIT 基 2FFTDIF 基 2FFT分裂基FFT混合基FFT128110994124051949301251092562

32、43137273330203335267000512559008628716453738574287102412182021376184983545129342510011099912405194930916122002431372733242033302318525005590086287164537385822271000123885113761849835441210904衣3対不同点数乞种FFT知法的误至(dB为爪位)点数 算法。1丁基 2FFTDIF 基 2FFT分裂基FFT混介基FFT128-255.5931940857301-251.6398783469884-234.65534

33、9585399565.26011347269841256-257.2115811610556-254.3612832750606-240.733581797111165.39990983698357512-260.9455955425411-251.067344385581-237.079081010208458.643437291787581024-245.3592667583791-240.0836989613196243.033374274061264.95967405994304100-262.9619790458805-259.8242773915911-250.4182470673

34、95262.57913977797035?00-260.1647R61R2612) 254.664064915681.243,182175614716267.18258128621514500-270.8009610444112-255.8232732633446235.201664147528258.318630445567851000259.158724085711-254.5790570309925-243.667816935925962.26979055514795从上而衣2可以看出，并种算法的运算砒都随着点数的增加而增加。当N-21时，无论 点数多还是少，分裂展FFT算法运算吊都优皋

35、2FFT,是算法复杂度最小的算法,这充分 显示了分裂某FFT算法在N=2L时在算法复杂度上的优势，爭实上，分裂幕算法是迄今为止 最仅仅理论复杂度最优的算法。而为N不等-2时时候，非混介某的FFT算法必须要补零 到2l点才能进行FFT,而H.不能进行精确长度FFT运算，而混合基则在FFT运算的长度匕 显示出很人的灵活性，由表中我们可以看出，混合基在保持与基2FFT差不多的运算灵的询 提卜，町以高效第实现100、200、500、1000点的FFT运算；|何在100点的FFT运算中， 其运算量甚至少F分裂基算法。把表2的数据用图表示如图15所示,从中可以明显地看出, 分裂基的运算最明显少其他算法。510X4 2 0 8 6 4 2 1 A 1400500600700运算长度（点）8009001000 1100200300图15对不同点数备种FFT n法的运舁厳图（运外时钟数为巾位表3中我们看出，基2FFT和分裂基FFT有较好的准确度，而混介基FFT则误差较人, 这是由于混合基的运算较为复朵造成的。6 总结以上我们完成了甚2FFT、分裂某FFT和混介展FFT算法的DSP实现，并利用语音信 号进行了实验。根据实验数据，我们总结出，对J' N=2l而分裂甚是最高效的算法，而对J： N不等于离点的FFT,混合基有其算法长度灵活的优势。参考文献1J1史贤勇陈