快速计算离散傅里叶变换

《快速计算离散傅里叶变换》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、第4章快速计算离散傅里叶变换4.1引言4.2基2FFT算法4.3进一步减少运算量的措施4.4分裂基FFT算法4.5离散哈特莱变换(DHT)4.1引言与序列的傅里叶变换相比，离散傅里叶变换有实用价值。但是按定义直接计算DFT有实用价值吗？只有一些。因为这种算法的计算数度太慢了。特别是与后人发明的算法相比，它的慢更显突出。1965年，J.W.Cooley和J.W.Tukey在MathematicsofComputation上发表了AnalgorithmforthemachinecalculationofcomplexFourierseries。它

2、极大的提高了计算离散傅里叶变换的速度。从定义来看N点长的DFT的运算量。1直接计算DFT需：复乘N2次，复加(N-1)N次。因为1个k需复乘N次，复加(N-1)次。对于复乘1次需50μs，复加1次需10μs的计算机，用直接法做N=1024点长的DFT共需多少时间？10242×50×10-6＋1023×1024×10×10-6=63(s)2Cooley和Tukey发明的方法计算DFT需：复乘(N/2)log2N次，复加Nlog2N次。用来计算上面的DFT共需多少时间？512×10×50×10-6＋1024×10×10×10-6=0.36(s)4

3、.2基2(radix2)FFT算法4.2.1直接计算DFT的特点及减少运算量的方法直接计算N个采样值的DFT需要有N2次复数乘法和N(N-1)次复数加法。如果把N分成几小段，降低DFT的规模，是不是可以大幅度地减少乘法和加法的运算次数？还有，WNkn具有对称性和周期性，是不是可以巧妙地利用？例如，当N=8时，从形式上看，W8kn共有64个值。但从图来看，Wkn实际上只有W0~W7这8个值是独立的；而且，其中有一半是对称的。科学家Cooley和Tukey正是巧妙地利用这些特性加快了DFT的运算速度。周期性：对称性：ImW86jW85W87-11

4、ReW84W80W83W81W82-j04.2.2时域抽取法基2FFT基本原理设序列x(n)的长度N=2M，M为自然数。(1)缩短DFT，把x(n)按n的奇偶顺序分成两半。则x(n)的DFT为(2)重组DFT，按DFT的定义重新组合变短的DFT。∵变短后的DFT中X1(k)和X2(k)分别为x1(r)和x2(r)的N/2点DFT，周期为N/2；∵对称性WNk+N/2=－WNk。∴X(k)又可表示为经过这两步骤处理后，1个N点的DFT就变成了2个N/2点的DFT。运算量变成：复乘(N/2)2×2+(N/2)≈N2/2次，复加(N/2)×(N/2

5、-1)×2+(N/2)×2=N2/2次。比原来多了还是少了？(4.2.7)(4.2.8)将式(4.2.7)和式(4.2.8)用流图符号表示，称为蝶形运算符号。采用蝶形符号可以表示N=8点的DFT运算，下面是经过1次分解的DFT的示意图。注意：上半部份有4点，用“＋”的公式做；下半部份有4点，用“－”的公式做。图4.2.28点DFT的一次时域抽取分解图2次分解x(n)的DFT:(1)缩短x1(r)和x2(r)的DFT，与第一次分解相同，将x1(r)按奇偶分解成两个N/22长的子序列x3(l)和x4(l)，即则x1(r)的DFT为(4.2.9)(

6、2)重组DFT，按DFT的定义重新组合变短的DFT。∵变短后的DFT中X3(k)和X4(k)分别为x3(l)和x4(l)的N/4点DFT，周期为N/22；∵对称性WN/2k+N/4=－WN/2k。∴X1(k)也可表示为用同样的方法可以计算出如果是8点的DFT，经两次分解，DFT的长度是多少？有几个这种长度的DFT？图4.2.38点DFT的第二次时域抽取分解图3次分解DFT，…，长度为N/23，…8点DIT―FFT运算流图需要几次分解DFT，才会使DFT变为1点的DFT？4.2.3时域抽取法快速傅里叶变换的运算量从分解的级来看——①每级需复乘N

7、/2次，？复加N次；？②M=log2N级需复乘N/2×M次，？复加N×M次。？对于复乘1次需50μs，复加1次需10μs的计算机，现在做N=1024点的DFT运算。按定义直接运算需要10242×50×10-6＋1023×1024×10×10-6=63(s)按DIT-FFT运算需要512×10×50×10-6＋1024×10×10×10-6=0.36(s)它们的速度相差63÷0.36=175（倍）！例如：分析序列x(n)=sin(1.8n)+cos(1.8n)的频谱。clear,closeall%用两种方法计算DFTn=0:1023;w=1.8

8、;x=sin(w*n)+cos(w*n);subplot(2,1,1),stem(n,x,'.');%axis([250,350,-1.5,1.5])w=linsp