1228.DFT快速算法的研究

1、 目录摘要.1第一部分 fft的回顾2第一节 dft与fft的过渡.2第二节 按时间抽取以2为基的算法(fft-dit2)（库利-图基算法）.2第三节 按频率抽取以2为基的算法(fft-dit2)（桑德-图基算法）.6第二部分 混合基与分裂基的研究.7第一节 混合基的相关知识7第二节 分裂基的相关知识9第三部分 fft算法改进.12第四部分 fft实例运行结果及比较14第五部分 总结与感想.18参考文献 .19摘要快速傅立叶变换fft是离散傅立叶变换dft的一种快速算法,我们知道有限长序列的主要特点是其频域也可以离散化成有限长序列，即能进行dft变换。dft广泛应用于信号的频谱分析、滤波器的设

2、计以及系统的分析、设计和实现中。而在相当长时间中由于n值大时，计算量大。无法达到实时处理的要求。1965年,库利(cooley)和图基(tukey)首先提出fft算法.对于n点dft,仅需(n/2)log2n 次复数乘法运算。大大提高了计算效率。后来又有桑德快速算法以及一些算法改进。fft是dft的一种快速计算机算法，无任何近似，无精度损失。本论文主要介绍基本fft算法与实现，以及混合基、分裂基等一些先进算法，最后附有方案改进和算法程序设计。第一部分 fft的算法回顾第一节 dft与fft的过渡 n点有限长序列x(n)的dft为二者仅有一符号之差和一个常熟因子1/n。由dft可知完成整个dft

3、运算共需n2次复数乘法和n(n-1)次复数加法。由复数知识可知，此运算4n2次实数乘法和 2n(n-1)次实数加法运算。因此dft 正比于n2，当n很大时，其运算量很可观。为此我们对此加以改进。 在dft中我们知道wnnk具有一些特定性质：可得出由此可对dft中的一些项合并，将其分解为短序列dft，fft就是由此而得出的。 第二节 按时间抽取以2为基的算法(fft-dit2)（库利-图基算法）设序列点数n=2l，lz。如条件不满足可人为加上若干个0，n=2l序列n的奇偶分两组 于是得经过wnnk的性质转换使x（k）前后分离如下：对x(k)的前一半（k=0,1, ,n/2-1）,由dft定义:对

4、x(k)的后一半（k=n/2, n/2+1, ,n-1）由此我们只需求0n/21区间的，就可求出x（k）。其算法流图如下图一。 x1（0）n/2 点 dft x(0) x（0） x1（1） wn0x(2) x（1） x1（2） wn1 x(4) x（2） x1（3） wn2x(6) x（3） y1（0） wn4 wn3n/2 点 dftx(1) x（4） y1（1） wn5x(3) x（5） y1（2） wn6x(5) x（6） y1（3） wn7x(7) x（7） 图一 n点离散傅立叶变换的计算化为两个n/2点离散傅立叶变换的计算（n=8）由以上流图可见，其fft的运算量大大减少。直接计算d

5、ft与fft算法的计算量比值见下表1。表1nn2n22414.041644.0864125.416256328.03210248012.864409619221.41281638444836.625665536102464.0由流图可知，输入为顺序，而输出为倒位序，其完成是通过变址运算来实现的。下表2是n=8时的顺序二进制数以及相应的倒位序二进制数。 表2自然数二进制数倒位序二进制数倒位序顺序0000000010011004201001023011110641000011510110156110011371111117通过上面的推导可以看出，n点dft可以分解为两个n/2点的dft，每个n/2

6、点的dft又可以分解为两个n/4点的dft。依此类推，当n为2的整数次幂时（n=2m），由于每分解一次降低一阶幂次，所以通过m次的分解，最后全部成为一系列2点dft运算。以上就是按时间抽取的快速傅立叶变换（fft）算法。序列x（k）的离散傅立叶变换的区别在于将wn改变为wn-1，并多了一个除仪n的运算。因为wn 和wn-1对于推导按时间抽取的快速傅立叶变换算法并无实质性区别，因为可将fft和快速傅立叶反变换（ifft）算法合并在同一个程序中。如下语句（c语言）： 1. 子函数语句void srfft(x，y，n，sign)2. 形参说明 x双精度实型一维数组，长度为n。开始时存放要变换数据的实

7、部，最后存放变换结果的实部。 y双精度实型一维数组，长度为n。开始时存放要变换数据的虚部，最后存放变换结果的虚部。 n整型变量。数据长度必须是2 的整数次幂，即n2m。sign整型变量。当sign时，子函数fft( )计算离散傅立叶正变换（dft）；当sign时，子函数计算离散傅立叶反变换（dft）。3、子函数程序（文件名：srfftc）#includemath.hvoid fft(x,y,n,sign);int n, sign;double x ,y ; int i,j,k,l,m ,n1,n2,;double c,c1,e,s,s1t,tr;ti;for(j=1,i=1;i16;i+) m

8、=i; j=2*j; if(j=n)break; n1=n-1;for(j=0,i=0;in;i+) if(ij) tr=xj;ti=yj;xj=xi;yj=yi;xi=tr;yi=ti; k=n/2; while(k(j+1)j=j-k; k=k/2; j=j+k; n1=1;for(l=1;l=m;l+)n1=2*n1; n2=n1/2; e=3n2; c=1.0; s=0.0; c1=cos(e); s1=-sign*sin(e); for(j=0;jn2;j+)for(i=j;in;i+=n1) k=i+n2; tr=c*xk-s*xk; xk=xi-tr;

9、yk=yi-ti; xi=xi+tr; yi+yi+ti; t=c;c=c*c1-s*s1;s=t*s1+s*c1; if(sign=-1) for(i=0;in;i+) xi/=n; yi/=n; 第三节 按频率抽取以2为基的算法(fft-dit2)（桑德-图基算法）算法原理：对时域有限长序列x(n)，设n（可补零）。1)按频率抽取以2为基的fft算法(fft-dif2)是以频域奇偶分组方式将x(k)不断二分为一系列短序列，减小了dft的基数。2)在每个短序列的计算中，利用对称性和周期性，以蝶形运算来合并相同项，从而降低计算量，提高了计算速度。由的性质可知按k的奇数和偶数可将x(k)分为两部

10、分。x(2k)x(2k)这样，我们就把一个n点dft分解为两个n/2点的dft了。其分解过程如下图二所示。x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)x(0)x(2)x(4)x(6)x(1)x(3)x(5)x(7)vn/2点dftn/2点dft 图二按频率抽取以2为基的fft算法(fft-dif2)与按时间抽取以2为基的fft算法的特点和程序基本类似，这里不在介绍。第二部分 混合基与分裂基的研究第一节 混合基的相关知识混合基算法：混合基算法主要应用于当序列值n为复合数时的快速傅立叶变换。它是将复合数n分解成一系列因子乘积的整数情况的fft算法。当n=ab,计算dft如下：

11、（1）其中n=ab=ba将十进制数n表示为b进制数，n0末位，n1进位，得n=n1b+n0, n1=0,1,n-1, n0=0,1,b-1 （2）同理将频率变量k表示为a进制，k1为进位，k0为末位。得k=k1a+k0, k1=0,1,b-1, k0=0,1,a-1 （3）将（2）,（3）代入（1）式中，得 由以上过程中得 ， k00，1，a-1 （4） （5） k10，1，b-1 （6） （7）以上是n为复合数的fft算法得全过程，从中我们归纳出它的一般计算步骤如下：a) 将n分为两个因子，例如n=ab,将序列变换为：x(n)=x(bn1+n0 )=x(n1,n0) n1=0，1，a-1 n

12、00，1，b-1；b) 利用式（4）求c) 利用式（5）求d) 利用式（6）求e) 利用式（7）求整序，的傅立叶变换以n428为例，得其混合基fft算法流程图，如图三所示：x(n)=x(n1,n0) x1(k0,n0) x2(k0,k1)=x(k) x1(k1,k0)=x(k) 0 00 00 00 0 00 0 1 01 01 wn0 01 4 01 1 -12 10 -j10 10 1 02 2 -1 3 11 j 11 wn1 11 5 03 3 -14 20 -1 20 20 2 10 4 1 -1 5 21 21 wn2 21 6 11 5 j -1 -16 30 -1 30 30

13、3 12 6 -j7 31 31 wn3 31 7 13 7 -1 图三 n=4x2=8的fft运算流程图上面讨论的是将n分解为两个素数的情况，当n为高组合素数是，同样可按以上步骤进行计算，如nam，mcd，根据以上流程同样可以得到其dft。以此类推可以到任意基的dft，课堂中学的基2fft混合基运算里的一个特例。第二节 分裂基的相关知识1. 分裂基概述：自从基2快速算法出现后，人们再不断探求更快的算法。1984年，法国的p.dohamel和h.hollmann将基2分解和基4分解糅合在一起，提出了分裂基fft算法。该算法的基本想法是对偶序号使用基2fft算法，对奇序号使用基4算法。其运算量比

14、前诉几种算法都有所减少，运算流程图却于基2ff很接近。它在fft算法中是一种实用的高效率新算法。2. 方法简介序列x(n)的离散傅立叶变换为将x(k)按序列号k的奇偶分成两组。当k为偶数时，进行基2频率抽取分解，x(k)可表示为x(2k)当k为奇数时，进行基4频率抽取分解，x(k)可表示为 由此可见，一个n点的dft可以分解为一个n2点的dft和两个n4点的dft。这种分解既有基2的部分，又有基4的部分，因此称为分裂基分解。上面的n2点dft又可分解为一个n8点的dft和两个n8点的dft，而两个n4点的dft也可以分解为一个n8点的dft和两个n16点的dft。依次类推，直至分解到最后一级为

15、止。这就是按频率抽取的分裂基快速傅立叶变换算法。3. 分裂基算法的优点：a) 在已知的n2m的各种算法中，所需的乘法数是最少的，并接近理论上的最小值；b) 分裂基算法有着和基4、基2算法一样的规则结构，可以同址运算，这在用ic芯片来实现这些算法时特别重要；c) 若把基2，基4和分裂基算法中无关紧要的旋转因子都考虑进去，那么三者所需的计算量是一样的。分裂基算法的特点时合理安排了算法结构，使无关紧要的旋转因子最大程度地减少。4. 分裂基c语言程序： 一。程序说明：1 子函数语句void srfft(x,y,n)2,形参说明 x双精度实型一维数组，长度为n。开始时存放要变换数据的实部，最后存放变换结

16、果的实部。 y双精度实型一维数组，长度为n。开始时存放要变换数据的虚部，最后存放变换结果的虚部。 n整型变量。数据长度必须是2 的整数次幂，即n2m。二。子函数程序（文件名：srfft。c）#include”math.h”void strfft(x,y,z)int n;double x,y;int i,j,k,m,i1,i2,i3,n1,n2,n4,id,is;double a,b,c,e,a3,r1,r2,r3,r4;double c1,c3,s1,s2,s3,cc1,cc3,ss1,ss3;for(j=1,i=1;i16;i+)m=i; j=2*j; if(j=n)break; n2=2*

17、n; for(k=1;km;k+) n2=n2/2;n4=n2/4;e=6.28318530718/n2;a=0;for(j=0;jn4;j+)a3=3*a;cc1=cos(a);ss1=sin(a);cc3=cos(a3);ss3=sin(a3);a=(j+1)*e;is=j;id=2*n2;do for(i=is;i(n-1);i=i+id) i1=i+n4;i2=i1+n4;i3=i2+n4;r1=xi-xi2;xi=xi+xi2;r2=xi1-xi3;xi1=xi1+xi3;s1=yi-yi2;yi=yi+yi2;s2=yi1-yi3;yi1=yi1+yi3;s3=r1-s2;r1=r

18、1+s2;s2=r2-s1;r2=r2+s1;xi2=r1*cc1-s2*ss1;yi2=-s2*cc1-r1*ss1;xi3=s3*cc3+r2*ss3;yi3=r2*cc3-s3*ss3; is=2*id-n2+j;id=4*id; while(is(n-1);is=0;id=4;dofor(i=is;in;i=i+id) i1=i+1;r1=xi;r2=yi;xi=r1+xi1;yi=r2+yi1;xi1=r1-xi1;yi1=r2-yi1; is=2*id-2; id=4*id;while(is(n-1);n1=n-1;for(j=0,i=0;in1;i+)if(ij) r1=xj;s

19、1=yj; xj =xi; yj=yi; xi=r1; yi=s1;k=n/2;while(k(j+1) j=j-k;k=k/2;j=j+k;)第三部分 fft算法改进自1965年重新发现快速傅立叶变换算法以来，寻求更快的算法的努力从来未停止过。从理论上讲，用更大的基数可以进一步的缩减实数运算的次数，但随之而来的程序的复杂程度和长度都必须付出代价，甚至于得不偿失。除了基8算法尚可实用外，基数取得更高是没有多大意义的。我们知道，dft的计算同周期褶积的计算有着密切的关系。有了fft算法，就可以快速地计算周期褶积。但是也可以直接寻求快速计算周期褶积地算法，甚至也可以用某些快速计算周期褶积地方法反过

20、来用于计算dft。在七十年代中期，上诉两方面都有所进展。前者是所谓数论变换(ntt)，用有限整数环取代了复数域，跟dft一样可以将周期褶积变换为乘积。数论变换可以避免运算误差，这在dft地三角函数地有关运算是不可避免的，但变换的有关参数受到数论上的要求而严格的限制。策二种是所谓的winograd傅立叶变换算法，它将较长点数的dft分解为许多短点数的dft。每种短点数的dft用短周期褶积来算出来，各个短周期褶积的计算则用多项式理论设计出具有尽可能少的乘法次数的算法。虽然在不断的实践中，这两种算法都不理想，但也促进了人们不断进行研究探索！以下是最新的离散傅立叶变换的改良方法：sift:离散傅立叶变

21、换asic/fpga设计的新方法傅立叶变换是数字时代的核心，在数字信号处理中扮演着重要的角色。自傅立叶变换理论提出以来，其计算方法一直是研究的重要课题。离散傅立叶变换是中更为可行和行之有效的傅立叶变换的实际应用。目前已有多种离散傅立叶变换的计算方法。如goertzel的递归算法,bluestein chrip z变换算法，rader算法，winograd 算法，hartly变换等。1965年后，cooley及tukey提出的fft算法成为傅立叶变换计算方法的主流。fftw(fastest fourier transform in the west)是fft最新的算法之一。fft提出之后，傅立叶

22、变换技术的进步更多地来自于微电子工艺的改进。快速傅立叶变换-fft 离散傅立叶变换的基本公式为： 其中 fft的核心是蝶型单元，其设计也是 fft实现的重点，在fpga实现中，是较为复杂的部分，其性能也决定了整个设计的性能。fft需要对样本或输出系数次序做特别的处理。fft使基本的dft计算量由2n2减少到2nlog(n)，计算量大大降低，成为傅立叶变换的关键因素。我们可以注意到，fft算法中，采样点对计算结果是彼此相关的，因此，也需要采集全部采样周期内的样本数据后，才开始运算。fft的计算过程是批处理过程，基本的流程是收集数据、开始计算、输出结果。典型的数据处理流程如图1所示。 fft算法的

23、fpga/asic实现，有很多的研究文献可供参考。 采样积分傅立叶变换 siftsift(sample integrated fourier transform)是以傅立叶变换基本算法为基础的新的离散傅立叶变换实现方法。sift方法中，当采样点对系数的贡献可以独立计算，计算完成后，不需要保存这些样本数据。数据处理是透明型的。由此得到的一个明显的好处在于，在fpga/asic设计中，可以实现无延迟傅立叶变换，全部采样数据输入完即可输出系数。sift处理数据的流程如图2所示.第四部分 fft实例运行结果及比较一、基2fft算法程序运行结果设输入序列x(i)为 x(i)=qi,i=0,1,n-1其离

24、散傅立叶变换为x(k)= (1-qn)/(1-qwk), k=0,1,2,31这里w=e-j2p/n.选取q=0.9+j0.3,n=32，计算离散傅立叶变换x=(k)和离散傅立叶反变换x(i)，并和理论值比较。主函数程序（文件名：fft.m）:#include #include main() inti,j,n; double a1,a2,c,c1,c2,d1,d2,q1,q2,w,w1,w2; double x32,y32,a32,b32; n=32;a1=0.9;a2=0.3;x0=1.0;y0=0.0;for(i=1;in;i+) xi=a1*xi-1-a2*yi-1; yi=a2*xi-

25、1+a1*yi-1; printf(n complex input sequencen); for(i=0;in/2;i+) for(j=0;j2;j+) printf(%10.7f+j%10.7f,x2*i+j,y2*i+j); printf(n)q1=xn-1;q2=yn-2;for(i=1;in;i+)w=6.28318530718/n*i; w1=cos(w); w2=-sin(w); c1=1.0-a1*w1+a2*w2; c2=a1*w2+a2*w1; c=c1*c1+c2*c2; d1=1.0-a1*q1+a2*q2; d2=a1*q2+a2*q1; ai=(c1*d1+c28d

26、2)/c; bi=(c2*d1-c1*d2)c; printf(n theoretical dftn); for(i=0;in/2;i+) for(j=0;j2;j+) printf(%10.7f+j%10.7f,x2*i+j,y2*i+j); printf(n) fft(x,y,n,l);printf(n discrete fourier transformn); for(i=0;in/2;i+) for(j=0;j2;j+) printf(%10.7f+j%10.7f,x2*i+j,y2*i+j); printf(n) fft(x,y,n,-l);printf(n inverse disc

27、rete fourier transformn); for(i=0;in/2;i+) for(j=0;j2;j+) printf(%10.7f+j%10.7f,x2*i+j,y2*i+j); printf(n) 运行结果为：输入序列x(i)为1.0000000 +j 0.0000000 0.9000000 +j 0.30000000.7200000 +j 0.5400000 0.4860000 +j 0.70200000.2268000 +j 0.7776000 -0.0291600 +j 0.7678800- 0.2566080 +j 0.6823440 -0.4356504 +j -0.5

28、371272- 0.5532235 +j 0.3527194 -0.6037170 +j 0.1514804- 0.5887894 +j - 0.0447828 -0.5164756 +j -0.2169413- 0.3997457 +j - 0.3501899 -0.2547141 +j -0.4350946- 0.0987144 +j - 0.4679994 0.0515569 +j - 0.4508137 0.1816453 +j - 0.3902653 0.2805604 +j - 0.2967452 0.3415279 +j - 0.1829025 0.3622459 +j - 0.

29、0621539 0.3446674 +j 0.0527352 0.2943801 +j 0.1508619 0.2196835 +j 0.2240898 0.1304882 +j 0.2675859 0.0371637 +j 0.2799738 -0.0505448 +j 0.2631255- 0.1244280 +j 0.2216495 -0.1784800 +j 0.1621561- 0.2092789 +j 0.0923965 -0.2160699 +j 0.0.203732- 0.2005749 +j - 0.0464851 -0.1665719 +j - 0.1020091x(i)离

30、散傅立叶变换x(k)(k=0,1,2,31)的理论值为 0.6939728 +j 3.4997157 2.7922679 +j 8.0504557 9.4029646 +j -9.1350136 1.8664455 +j -3.8338328 1.131827 +j -2.2341573 0.9047939 +j -1.5346293 0.7995572 +j -1.1396094 0.7396056 +j -0.8823144 0.7008616 +j -0.6985654 0.6735758 +j -0.5584785 0.6531094 +j -0.4462450 0.6369913 +

31、j -0.3526891 0.6237884 +j -0.2720860 0.6126129 +j -0.2006419 0.6028833 +j -0.1357032 0.5942004 +j -0.0753137 0.5862775 +j -0.0179492 0.5788996 +j 0.0376517 0.5718985 +j 0.0926070 0.5651358 +j 0.1479833 0.5584921 +j 0.2048808 0.5518591 +j 0.2645222 0.5451336 +j 0.3283649 0.5382144 +j 0.3982571 0.5310

32、015 +j 0.4766778 0.5234037 +j 0.5671323 0.5153625 +j 0.6748500 0.5069258 +j 0.8081015 0.4984670 +j 0.9809063 0.4913894 +j 1.2192074 0.4907322 +j 1.5770820 0.5173540 +j 2.1888329傅立叶变换x(k)(k=0,1,2,31)为0.6939728 +j 3.4997157 2.7922679 +j 8.0504557 9.4029646 +j -9.1350136 1.8664455 +j -3.8338328 1.13182

33、7 +j -2.2341573 0.9047939 +j -1.5346293 0.7995572 +j -1.1396094 0.7396056 +j -0.8823144 0.7008616 +j -0.6985654 0.6735758 +j -0.5584785 0.6531094 +j -0.4462450 0.6369913 +j -0.3526891 0.6237884 +j -0.2720860 0.6126129 +j -0.2006419 0.6028833 +j -0.1357032 0.5942004 +j -0.0753137 0.5862775 +j -0.0179

34、492 0.5788996 +j 0.0376517 0.5718985 +j 0.0926070 0.5651358 +j 0.1479833 0.5584921 +j 0.2048808 0.5518591 +j 0.2645222 0.5451336 +j 0.3283649 0.5382144 +j 0.3982571 0.5310015 +j 0.4766778 0.5234037 +j 0.5671323 0.5153625 +j 0.6748500 0.5069258 +j 0.8081015 0.4984670 +j 0.9809063 0.4913894 +j 1.21920

35、74 0.4907322 +j 1.5770820 0.5173540 +j 2.1888329傅立叶反变换x(i)(i=0,1,2,31)为1.0000000 +j 0.0000000 0.9000000 +j 0.30000000.7200000 +j 0.5400000 0.4860000 +j 0.70200000.2268000 +j 0.7776000 -0.0291600 +j 0.7678800- 0.2566080 +j 0.6823440 -0.4356504 +j -0.5371272- 0.5532235 +j 0.3527194 -0.6037170 +j 0.151

36、4804- 0.5887894 +j - 0.0447828 -0.5164756 +j -0.2169413- 0.3997457 +j - 0.3501899 -0.2547141 +j -0.4350946- 0.0987144 +j - 0.4679994 0.0515569 +j - 0.4508137 0.1816453 +j - 0.3902653 0.2805604 +j - 0.2967452 0.3415279 +j - 0.1829025 0.3622459 +j - 0.0621539 0.3446674 +j 0.0527352 0.2943801 +j 0.1508

37、619 0.2196835 +j 0.2240898 0.1304882 +j 0.2675859 0.0371637 +j 0.2799738 -0.0505448 +j 0.2631255- 0.1244280 +j 0.2216495 -0.1784800 +j 0.1621561- 0.2092789 +j 0.0923965 -0.2160699 +j 0.0.203732- 0.2005749 +j - 0.0464851 -0.1665719 +j - 0.1020091二、分裂基算法程序运行设输入序列x(i)为 x(i)=qi,i=0,1,n-1其离散傅立叶变换为x(k)= (

38、1-qn)/(1-qwk),k=0,1,2,31这里w=e-j.选取q=0.9+j0.3, n=32，用分裂基快速算法计算离散傅立叶变换x=(k)。 主函数程序（文件名：srfft.m）:#include #include main() int i,j,n; double a1,a2,x32,y32; n=32;a1=0.9;a2=0.3;x0=1.0;y0=0.0;for(i=1;in;i+) xi=a1*xi-1-a2*yi-1; yi=a2*xi-1+a1*yi-1; srfft(x,y,n); printf(n discrete rourier transformn); for(i=0

39、;in/2;i+) for(j=0;j2;j+) printf(%10.7f+j%10.7f,x2*i+j,y2*i+j); printf(n) 运行结果为：傅立叶变换x(k)(k=0,1,2,31)为0.6939728 +j 3.4997157 2.7922679 +j 8.0504557 9.4029646 +j -9.1350136 1.8664455 +j -3.8338328 1.131827 +j -2.2341573 0.9047939 +j -1.5346293 0.7995572 +j -1.1396094 0.7396056 +j -0.8823144 0.7008616 +j -0.6985654 0.6735758 +j -0.5584785 0.6531094 +j -0.4462450 0.6369913 +j -0.3526891 0.6237884 +j -0.2720860 0.612612