关于DFT与FFT运算速度的比较

1、离散傅里叶变换赵德华（哈尔滨工业大学 电子与信息工程学院 0705201班）摘要：本文简要介绍了DFT的原理及其FFT实现算法。以实例验证DFT与FFT的等效性，并在VC6.0环境下对DFT与FFT的计算时间做了系统比较，验证理论得出的关系。关键词：DFT，FFT，计算速度，C程序Discrete Fourier TransformZhao Dehua(Electronics and Information Department, Harbin Institute of Technology)Abstract：This paper briefly introduces the principl

2、e of DFT and FFT.It verifies the equivalent of DFT and FFT with sevel examples.It also comples the elapsed time of DFT and FFT in VC6.0 to verify the theory.Key words: DFT, FFT, calculating speed , C program 目录0.引言11.DFT与FFT的定义.12.DFT与FFT等价性验证.23.DFT与FFT运算速度比较.34.总结.6参考文献60.引言离散傅里叶变换（Discrete Fourie

3、r Transform），是连续傅里叶变换在时域和频域上都离散的形式，将时域信号的采样变换为在离散时间傅里叶变换（DTFT）频域的采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作经过周期延拓成为周期信号再作变换。离散傅里叶变化的出现解决了信号离散化的问题，从而使其在数字滤波、功率谱分析、仿真、系统分析、通信方面得以应用。快速傅里叶变换（Fast Fourier Transform），是离散傅里叶变换的快速算法，它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改

4、进获得的。由于它应用的理论基础仍是离散傅里叶变换，所以它对离散傅里叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。快速傅里叶变化最早于1965年由Cooly和Tukey提出，这使离散傅里叶变化运算次数由N2减少为Nlog2N次，使得离散傅里叶变换应用于实际变成现实。目前，快速傅里叶变化技术已广泛应用于各个领域，成为数字信号处理技术的一个重要组成部分。离散傅里叶变化除了本文介绍的基2快速傅里叶变化算法外，他们都不同程度上减少了运算次数，如戈泽尔（Goertzel）算法、Chirp-Z变化（CZT）算法、Winograd Fast Trans

5、form Alogrithm（WFTA）等。本文将仅介绍基2快速傅里叶变化。1. DFT与FFT的定义DFT的定义：设是连续函数的个抽样值，这N个点的宽度为N的DFT为：。与之对应的反变换，即IDFT的定义：设是连续频率函数的个抽样值， 这N个点的宽度为N的IDFT为：。其中称为N点DFT的变换核函数，称为N点IDFT的变换核函数。它们互为共轭。如果引入，正逆变换的核函数分别可以表示为和。DFT可以表示为：，IDFT可以表示为：。用FFT实现DFT：先设序列点数为N=2M，M为整数。如果不满足这个条件，可以人为地加上若干零值点，使之达到这一要求。这种N为2的整数幂的FFT称基2 FFT。设输入

6、序列长度为N=2M (M为正整数) ，将该序列按时间顺序的奇偶分解为越来越短的子序列，称为按时间抽取(DIT )的FFT算法。下面给出8点基2 FFT的运算流图：图1-1 下面对两种变化的运算量进行比较。设x(n)为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FF

7、T中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成2*（N/2）2=N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是

8、FFT的优越性。2. DFT与FFT等价性验证如上面论述的，FFT只是DFT的快速实现算法，它仍是以DFT为理论基础的，因此他们两者完全等价。这里将在VC6.0下编制DFT与FFT 的程序，用其对两组离散序列进行变化，并输出变换结果。矩形脉冲信号，序列长64，其中前16点为1，其余为0。其DFT与FFT前5点输出如下图：图2-1 正弦脉冲信号，序列长61，数字角频率为PI/2。其DFT与FFT前5点输出如下图：图2-2比较可以发现用DFT和FFT对同一序列变化得到的结果是完全一样的。这说明DFT与FFT是完全等价的，而区别仅在于其具体算法。3. DFT与FFT运算速度比较如上面介绍的，N点DF

9、T需进行N2次运算，而N（N=2M）点FFT只需要Nlog2N次运算。因此DFT的运算时间与FFT的运算时间比应为：N/log2N。这里将对上述结论在数量级层次进行验证。为了严谨，这里先验证DFT与FFT运算时间仅与N有关，而与具体信号类型无关。这里分别对128点、256点、512点矩形脉冲信号和正弦脉冲信号进行DFT与FFT变换，运算时间如图：图3-1图3-2图3-3可以看出，对不同类型信号进行DFT与FFT变换，只要其点数相同，其运算时间也是大致相同的。考虑到C程序计时并不如汇编等机器语言准确，且100次重复计算会对每次运算时间之差进行累加，上面3次试验反映出的时间差是可以忽略的。由此可见

10、，DFT与FFT运算时间仅与运算点数N有关，而与具体信号类型无关。当信号点数较多时，对不同类型信号进行DFT与FFT变化将花费大量时间，而上面已经验证DFT与FFT运算时间仅与N有关，而与具体信号类型无关。故下面将仅对同一类型不同点数信号的DFT与FFT的运算时间的比较，而这也具有一般代表性。下面将仅对不同点数正弦信号做DFT与FFT变化，并进行运算时间比较。由于C程序中计时受机器滴答时间限制，故对点数较小信号进行较多次数DFT与FFT变换。而当点数较多时，DFT运算时间将按N2规律增长，故对点数较多信号进行较少次数DFT与FFT变换。各次DFT与FFT运算时间如下诸图：图3-4图3-5图3-

11、6图3-7图3-8图3-9图3-10图3-11图3-12图3-13图3-14图3-15图3-16图3-17下面对以上数据进行整理分析。 （注：e+n表示*10+n）N变换次数TDFT/sTFFT/smean(TDFT)/smean(TFFT)/sTDFT/TFFT理论比值2100000.0320.0313.2e-63.1e-61.0324100000.1560.0631.56e-56.3e-62.482810000.0620.0166.2e-51.6e-53.8752.671610000.250.0472.5e-44.7e-55.3243210001.0160.1091.016e-31.09e

12、-49.326.4641000.4060.0164.06e-31.6e-425.3810.671281001.7030.0631.703e-26.3e-427.0318.292561006.6090.1416.609e-21.41e-346.873251210026.2180.3132.6218e-13.13e-383.7656.8910241010.5310.0631.05316.3e-3167.16102.420481042.1720.1564.21721.56e-2270.331860470.3441.67047e+13.44e-2485.60341.33819

13、2167.1410.0786.7141e+17.8e-2860.786304850.1562.67485e+21.56e-11714.651170.29表3-1在Matlab中对以上数据绘图如下：图3-18 理论分析中已经指出，N点DFT运算次数将按N2规律增长，而N点FFT以Nlog2N规律增长，两者运算时间比值为N/log2N。当N较小时，两者运算时间尚可比拟。但随着N的增大，FFT相对DFT减少的运算量相对值就越大。在N=1024时，运算次数已相差两个数量级，而当N=16384时，运算次数更是相差达三个数量级。故信号序列点数越多，FFT相对DFT越有优势。这一点

14、在上述试验中得到了充分证实。当N1024后一次DFT运算时间就已经超过1秒。事实上，在实际工程应用时不可能只进行一次DFT变化，因此即便在短序列分析中FFT的优势也是明显的。上述实验中，DFT与FFT运算时间比与理论值有一定偏差。其原因有以下两点：一，程序与理论推导不同，程序难免在基于理论的同时存在一定的冗余，而两者冗余度的不同将直接影响其运算时间比；二，本次试验的实验环境是VC6.0，其在计时方面本身就劣于汇编等机器语言。尽管如此，在数量级层面试验反映的结果还是和理论分析一致的。此外，由于N=2M，M为整数，DFT的运算时间以N2规律增长，即本实验DFT运算时间应以4为倍数增长，这一点在实验

15、数据中有很好的体现。4. 总结由上面的分析我们可以看出，快速傅里叶变化是离散傅里叶变化的一种高速算法，它的理论基础仍然是离散傅里叶变换，而并无任何创新。但是由于快速傅里叶变化充分利用了离散傅里叶变换的奇、偶、虚、实等特性，并对离散傅立叶变换的算法进行改进，使得离散傅里叶变化运算次数由N2减少为Nlog2N次，使得离散傅里叶变换应用于实际变成现实。参考文献1Alan V.Oppenheim Alan S.Willsky with S. Hamid Nawab. Signals & Systems. Prentice-Hall,1997;2Alan V.Oppenheim Ronald W.Sha

16、fer with Jhon R.Buck. Discrete-Time Signal Processing. Prentice-Hall,1999;3赵淑清. 郑薇. 随机信号分析. 哈尔滨工业大学出版社. 1999;附录程序#include#include#include#define N 256#define M_PI 3.14159#define EXP 2.7182818288 void FFT(float *xr,float *xi,float *yr,float *yi,int n,int inv);void DFT(int n,float *x,int m,float *yr,f

17、loat *yi);void SignalRect(int n,int m, float *yr,float *yi);void SignalCos(int n,int f, float *yr,float *yi);void tra(float *x,float *y);main()/*以不同程序实现不同功能，具体见下文*/void SignalRect(int n,int m, float *yr,float *yi) int i; for(i=0;in;i+) yri=yi0=0; for(i=0;im;i+) yri=1; void SignalCos(int n,int f, flo

18、at *yr,float *yi) int i; for(i=0;in;i+) yri=yi0=0; for(i=0;in;i+) yri=cos(2*M_PI*f*i/n); void DFT(int n,float *x,int m,float *yr,float *yi) int j,k; for(k=0;km;k+) for(j=0;jn;j+) yrk=yrk+xj*cos(2*M_PI*k*j/m); yik=yik+xj*sin(2*M_PI*k*j/m); return; void FFT(float *xr,float *xi,float *yr,float *yi,int

19、n,int inv) int i,j,a,b,k,m; int ep,arg,mt,s0,s1; float sign,pr,pi,ph; float *c,*s; c=(float *)calloc(n,sizeof(float); if(c=NULL)exit(1); s=(float *)calloc(n,sizeof(float); if(s=NULL)exit(1); j=0; if(inv=0) sign=1.0; for(i=0;in;i+) yri=xri/n; yii=xii/n; else sign=-1.0; for(i=0;in-1;i+) if(ij) tra(&yr

20、i,&yrj); tra(&yii,&yij); k=n/2; while(k=j) j=j-k; k=k/2; j=j+k; ep=0; i=n; while(i!=1) ep=ep+1; i=i/2; ph=2*M_PI/n; for(i=0;in;i+) si=sign*sin(ph*i); ci=cos(ph*i); a=2; b=1; for(mt=1;mt=ep;mt+) s0=n/a; s1=0; for(k=0;kb;k+) i=k; while(in) arg=i+b; if(k=0) pr=yrarg; pi=yiarg; else pr=yrarg*cs1-yiarg*s

21、s1; pi=yrarg*ss1+yiarg*cs1; yrarg=yri-pr; yiarg=yii-pi; yri=yri+pr; yii=yii+pi; i=i+a; s1=s1+s0; a=2*a; b=b*2; free(c); free(s);void tra(float *x,float *y)float t;t=(*x);(*x)=(*y);(*y)=t;验证DFT与FFT等价性主函数：main()float yrN=0.0,yiN=0.0,xrN=0.0,xiN=0.0;int i;SignalCos(N,16,xr,xi); /*产生64点cos(pi*i/2)离散信号*/

22、* SignalRect(64,16,xr,xi);产生矩形脉冲*/DFT(N,xr,N,yr,yi);printf(the out put (by DFT,the first 5) is );for(i=0;i5;i+)printf(%f+ %fj,yri,yii);printf(n);FFT(xr,xi,yr,yi,N,0);printf(the out put (by FFT,the first 5) is );for(i=0;i5;i+)printf(%f+ %fj,yri,yii);printf(n);getch();验证DFT与FFT运算时间仅与N有关主程序：main()float

23、 yrN=0.0,yiN=0.0,xrN=0.0,xiN=0.0;int i=100;clock_t start, end;SignalCos(N,3,xr,xi); /* SignalCos(N,3,xr,xi) or SignalRect(64,16,xr,xi); */start=1*clock();while(i-)FFT(xr,xi,yr,yi,N,0);end=1*clock();printf(The time for SignalCos by FFT(N=%d,100 times) was: %f secondsn, N,(double)(end - start)/CLK_TCK);i=100;start=1*clock();while(i-)DFT(N,xr,N,yr,yi);end=1*clock();printf(The time for SignalCos by DFT(N=%d,100 times) was: %f secondsn, N,(double)(end - start)/CLK_TCK);i=100;SignalRect(64,16,xr,xi);start=1*clock();while(i-)FFT(xr,xi,yr,yi,N,0);end=1*clock();printf(The time for SignalRect by