关于DFT与FFT运算速度的比较.docx
- 文档编号:5876860
- 上传时间:2023-01-01
- 格式:DOCX
- 页数:17
- 大小:420.69KB
关于DFT与FFT运算速度的比较.docx
《关于DFT与FFT运算速度的比较.docx》由会员分享,可在线阅读,更多相关《关于DFT与FFT运算速度的比较.docx(17页珍藏版)》请在冰豆网上搜索。
关于DFT与FFT运算速度的比较
离散傅里叶变换
赵德华
(哈尔滨工业大学电子与信息工程学院0705201班)
摘要:
本文简要介绍了DFT的原理及其FFT实现算法。
以实例验证DFT与FFT的等效性,并在VC6.0环境下对DFT与FFT的计算时间做了系统比较,验证理论得出的关系。
关键词:
DFT,FFT,计算速度,C程序
DiscreteFourierTransform
ZhaoDehua
(ElectronicsandInformationDepartment,HarbinInstituteofTechnology)
Abstract:
ThispaperbrieflyintroducestheprincipleofDFTandFFT.ItverifiestheequivalentofDFTandFFTwithsevelexamples.ItalsocomplestheelapsedtimeofDFTandFFTinVC6.0toverifythetheory.
Keywords:
DFT,FFT,calculatingspeed,Cprogram
目录
0.引言…………………………………………………………………………………………………………………………………………………………1
1.DFT与FFT的定义………………………………………………………………………………………………………………………………..1
2.DFT与FFT等价性验证………………………………………………………………………………………………………………………..2
3.DFT与FFT运算速度比较…………………………………………………………………………………………………………………….3
4.总结………………………………………………………………………………………………………………………………………………………..6
参考文献……………………………………………………………………………………………………………………………………………………6
0.引言
离散傅里叶变换(DiscreteFourierTransform),是连续傅里叶变换在时域和频域上都离散的形式,将时域信号的采样变换为在离散时间傅里叶变换(DTFT)频域的采样。
在形式上,变换两端(时域和频域上)的序列是有限长的,而实际上这两组序列都应当被认为是离散周期信号的主值序列。
即使对有限长的离散信号作DFT,也应当将其看作经过周期延拓成为周期信号再作变换。
离散傅里叶变化的出现解决了信号离散化的问题,从而使其在数字滤波、功率谱分析、仿真、系统分析、通信方面得以应用。
快速傅里叶变换(FastFourierTransform),是离散傅里叶变换的快速算法,它是根据离散傅里叶变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。
由于它应用的理论基础仍是离散傅里叶变换,所以它对离散傅里叶变换的理论并没有新的发现,但是对于在计算机系统或者说数字系统中应用离散傅立叶变换,可以说是进了一大步。
快速傅里叶变化最早于1965年由Cooly和Tukey提出,这使离散傅里叶变化运算次数由N2减少为Nlog2N次,使得离散傅里叶变换应用于实际变成现实。
目前,快速傅里叶变化技术已广泛应用于各个领域,成为数字信号处理技术的一个重要组成部分。
离散傅里叶变化除了本文介绍的基2快速傅里叶变化算法外,他们都不同程度上减少了运算次数,如戈泽尔(Goertzel)算法、Chirp-Z变化(CZT)算法、WinogradFastTransformAlogrithm(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称基2FFT。
设输入序列长度为N=2M(M为正整数),将该序列按时间顺序的奇偶分解为越来越短的子序列,称为按时间抽取(DIT)的FFT算法。
下面给出8点基
2FFT的运算流图:
图1-1
下面对两种变化的运算量进行比较。
设x(n)为N项的复数序列,由DFT变换,任一X(m)的计算都需要N次复数乘法和N-1次复数加法,而一次复数乘法等于四次实数乘法和两次实数加法,一次复数加法等于两次实数加法,即使把一次复数乘法和一次复数加法定义成一次“运算”(四次实数乘法和四次实数加法),那么求出N项复数序列的X(m),即N点DFT变换大约就需要N2次运算。
当N=1024点甚至更多的时候,需要N2=1048576次运算,在FFT中,利用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%,点数越多,运算量的节约就越大,这就是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点DFT需进行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次试验反映出的时间差是可以忽略的。
由此可见,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-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/s
TFFT/s
mean(TDFT)/s
mean(TFFT)/s
TDFT/TFFT
理论比值
2
10000
0.032
0.031
3.2e-6
3.1e-6
1.03
2
4
10000
0.156
0.063
1.56e-5
6.3e-6
2.48
2
8
1000
0.062
0.016
6.2e-5
1.6e-5
3.875
2.67
16
1000
0.25
0.047
2.5e-4
4.7e-5
5.32
4
32
1000
1.016
0.109
1.016e-3
1.09e-4
9.32
6.4
64
100
0.406
0.016
4.06e-3
1.6e-4
25.38
10.67
128
100
1.703
0.063
1.703e-2
6.3e-4
27.03
18.29
256
100
6.609
0.141
6.609e-2
1.41e-3
46.87
32
512
100
26.218
0.313
2.6218e-1
3.13e-3
83.76
56.89
1024
10
10.531
0.063
1.0531
6.3e-3
167.16
102.4
2048
10
42.172
0.156
4.2172
1.56e-2
270.33
186.18
4096
10
167.047
0.344
1.67047e+1
3.44e-2
485.60
341.33
8192
1
67.141
0.078
6.7141e+1
7.8e-2
860.78
630.15
16384
1
267.485
0.156
2.67485e+2
1.56e-1
1714.65
1170.29
表3-1
在Matlab中对以上数据绘图如下:
图3-18
理论分析中已经指出,N点DFT运算次数将按N2规律增长,而N点FFT以Nlog2N规律增长,两者运算时间比值为N/log2N。
当N较小时,两者运算时间尚可比拟。
但随着N的增大,FFT相对DFT减少的运算量相对值就越大。
在N=1024时,运算次数已相差两个数量级,而当N=16384时,运算次数更是相差达三个数量级。
故信号序列点数越多,FFT相对DFT越有优势。
这一点在上述试验中得到了充分证实。
当N<256时,一次DFT运算时间还较小,而当N>1024后一次DFT运算时间就已经超过1秒。
事实上,在实际工程应用时不可能只进行一次DFT变化,因此即便在短序列分析中FFT的优势也是明显的。
上述实验中,DFT与FFT运算时间比与理论值有一定偏差。
其原因有以下两点:
一,程序与理论推导不同,程序难免在基于理论的同时存在一定的冗余,而两者冗余度的不同将直接影响其运算时间比;二,本次试验的实验环境是VC6.0,其在计时方面本身就劣于汇编等机器语言。
尽管如此,在数量级层面试验反映的结果还是和理论分析一致的。
此外,由于N=2M,M为整数,DFT的运算时间以N2规律增长,即本实验DFT运算时间应以4为倍数增长,这一点在实验数据中有很好的体现。
4.总结
由上面的分析我们可以看出,快速傅里叶变化是离散傅里叶变化的一种高速算法,它的理论基础仍然是离散傅里叶变换,而并无任何创新。
但是由于快速傅里叶变化充分利用了离散傅里叶变换的奇、偶、虚、实等特性,并对离散傅立叶变换的算法进行改进,使得离散傅里叶变化运算次数由N2减少为Nlog2N次,使得离散傅里叶变换应用于实际变成现实。
参考文献
[1]AlanV.OppenheimAlanS.WillskywithS.HamidNawab.Signals&Systems.Prentice-Hall,1997;
[2]AlanV.OppenheimRonaldW.ShaferwithJhonR.Buck.Discrete-TimeSignalProcessing.Prentice-Hall,1999;
[3]赵淑清.郑薇.随机信号分析.哈尔滨工业大学出版社.1999;
附录
程序
#include
#include
#include
#defineN256
#defineM_PI3.14159
#defineEXP2.7182818288
voidFFT(float*xr,float*xi,float*yr,float*yi,intn,intinv);
voidDFT(intn,float*x,intm,float*yr,float*yi);
voidSignalRect(intn,intm,float*yr,float*yi);
voidSignalCos(intn,intf,float*yr,float*yi);
voidtra(float*x,float*y);
main()
{
/*以不同程序实现不同功能,具体见下文*/
}
voidSignalRect(intn,intm,float*yr,float*yi)
{
inti;
for(i=0;i yr[i]=yi[0]=0; for(i=0;i yr[i]=1; } voidSignalCos(intn,intf,float*yr,float*yi) { inti; for(i=0;i yr[i]=yi[0]=0; for(i=0;i yr[i]=cos(2*M_PI*f*i/n); } voidDFT(intn,float*x,intm,float*yr,float*yi) { intj,k; for(k=0;k { for(j=0;j { yr[k]=yr[k]+x[j]*cos(2*M_PI*k*j/m); yi[k]=yi[k]+x[j]*sin(2*M_PI*k*j/m); } } return; } voidFFT(float*xr,float*xi,float*yr,float*yi,intn,intinv) { inti,j,a,b,k,m; intep,arg,mt,s0,s1; floatsign,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;i { yr[i]=xr[i]/n; yi[i]=xi[i]/n; } } elsesign=-1.0; for(i=0;i { if(i { tra(&yr[i],&yr[j]); tra(&yi[i],&yi[j]); } 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;i { s[i]=sign*sin(ph*i); c[i]=cos(ph*i); } a=2; b=1; for(mt=1;mt<=ep;mt++) { s0=n/a; s1=0; for(k=0;k { i=k; while(i { arg=i+b; if(k==0) { pr=yr[arg]; pi=yi[arg]; } else { pr=yr[arg]*c[s1]-yi[arg]*s[s1]; pi=yr[arg]*s[s1]+yi[arg]*c[s1]; } yr[arg]=yr[i]-pr; yi[arg]=yi[i]-pi; yr[i]=yr[i]+pr; yi[i]=yi[i]+pi; i=i+a; } s1=s1+s0; } a=2*a; b=b*2; } free(c); free(s); } voidtra(float*x,float*y) { floatt; t=(*x); (*x)=(*y); (*y)=t; } 验证DFT与FFT等价性主函数: main() { floatyr[N]={0.0},yi[N]={0.0},xr[N]={0.0},xi[N]={0.0}; inti; SignalCos(N,16,xr,xi);/*产生64点cos(pi*i/2)离散信号*/ /*SignalRect(64,16,xr,xi);产生矩形脉冲*/ DFT(N,xr,N,yr,yi); printf("theoutput(byDFT,thefirst5)is"); for(i=0;i<5;i++) printf("%f+%fj",yr[i],yi[i]); printf("\n"); FFT(xr,xi,yr,yi,N,0); printf("theoutput(byFFT,thefirst5)is"); for(i=0;i<5;i++) printf("%f+%fj",yr[i],yi[i]); printf("\n"); getch(); } 验证DFT与FFT运算时间仅与N有关主程序: main() { floatyr[N]={0.0},yi[N]={0.0},xr[N]={0.0},xi[N]={0.0}; inti=100; clock_tstart,end; SignalCos(N,3,xr,xi);/*SignalCos(N,3,xr,xi)orSignalRect(64,16,xr,xi);*/ start=1*clock(); while(i--) { FFT(xr,xi,yr,yi,N,0); } end=1*clock(); printf("ThetimeforSignalCosbyFFT(N=%d,100times)was: %fseconds\n",N,(double)(end-start)/CLK_TCK); i=100; start=1*clock(); while(i--) { DFT(N,xr,N,yr,yi); } end=1*clock(); printf("ThetimeforSignalCosbyDFT(N=%d,100times)was: %fseconds\n",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("ThetimeforSignalRectbyFFT(N=%d,100times)was: %fseconds\n",N,(double)(end-start)/CLK_TCK); i=100; start=1*clock(); while(i--) { DFT(N,xr,N,yr,yi); } end=1*clock(); printf("ThetimeforSignalRectbyDFT(N=%d,100times)was: %fseconds\n",N,(double)(end-start)/CLK_TCK); getch(); } 对不同点数正弦信号DFT与FFT运算时间比较主程序(需对N,i赋值进行相应调整): main() { floatyr[N]={0.0},yi[N]={0.0},xr[N]={0.0},xi[N]={0.0}; inti=10000; clock_tstart,end; SignalCos(N,3,xr,xi);/*SignalCos(N,3,xr,xi)orSignalRect(64,16,xr,xi);*/ start=1*clock(); while(i--) { FFT(xr,xi,yr,yi,N,0); } end=1*clock(); printf("ThetimeforSignalCosbyFFT(N=%d,10000times)was: %fseconds\n",N,(double)(end-start)/CLK_TCK); i=10000; start=1*clock(); while(i--) { DFT(N,xr,N,yr,yi); } end=1*clock(); printf("ThetimeforSignalCosbyDFT(N=%d,10000times)was: %fseconds\n",N,(double)(end-start)/CLK_TCK); getch(); }
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 关于 DFT FFT 运算 速度 比较