快速计算离散傅里叶变换课件

第第4章章 快速傅里叶变换快速傅里叶变换(FFT)第第4章章 快速计算离散傅里叶变换快速计算离散傅里叶变换 4.1 引言引言4.2 基基2FFT算法算法4.3 进一步减少运算量的措施进一步减少运算量的措施4.4 分裂基分裂基FFT算法算法4.5 离散哈特莱变换离散哈特莱变换(DHT)第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.1 引言引言 与与序序列列的的傅傅里里叶叶变变换换相相比比，离离散散傅傅里里叶叶变变换换有有实实用用价价值值。但但是是按按定定义义直直接接计计算算DFT有有实实用用价价值值吗吗？只只有有一一些些。因因为为这这种种算算法法的的计计算算数数度度太太慢慢了了。特特别别是是与与后人发明的算法相比，它的慢更显突出。后人发明的算法相比，它的慢更显突出。1965年年，J.W.Cooley 和和 J.W.Tukey在在Mathematics of Computation上上发发表表了了An algorithm for the machine calculation of complex Fourier series。它极大的提高了计算离散傅里叶变换的速度。它极大的提高了计算离散傅里叶变换的速度。第第4章章 快速傅里叶变换快速傅里叶变换(FFT)从定义来看从定义来看N点长的点长的DFT的运算量。的运算量。1 直接计算直接计算DFT需：复乘需：复乘N2次，复加次，复加(N-1)N次。因次。因为为 1个个k需复乘需复乘N次，复加次，复加(N-1)次。次。对于复乘对于复乘1次需次需50s，复加，复加1次需次需10s的计算机，用直接的计算机，用直接法做法做N=1024点长的点长的DFT共需多少时间？共需多少时间？1024250 10-6 102310241010-6=63(s)2 Cooley和和Tukey发明的方法计算发明的方法计算DFT需：复乘需：复乘(N/2)log2N次，复加次，复加Nlog2N次。用来计算上面的次。用来计算上面的DFT共需多少时间？共需多少时间？51210 50 10-6 1024101010-6=0.36(s)第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.2 基基2(radix2)FFT算法算法 4.2.1 直接计算直接计算DFT的特点及减少运算量的方法的特点及减少运算量的方法 直接计算直接计算N个采样值的个采样值的DFT 需要有需要有N2次复数乘法和次复数乘法和N(N-1)次复数加法。次复数加法。如果把如果把N分成几小段，降低分成几小段，降低DFT的规模，是不是可以大的规模，是不是可以大幅度地减少乘法和加法的运算次数？幅度地减少乘法和加法的运算次数？还有，还有，WNkn具有对称性和周期性，是不是可以巧妙地具有对称性和周期性，是不是可以巧妙地利用？利用？第第4章章 快速傅里叶变换快速傅里叶变换(FFT)例如，当例如，当N=8时，从形式上看，时，从形式上看，W8kn共有共有64个值。但从个值。但从图来看，图来看，Wkn实际上只有实际上只有W0W7这这8个值是独立的；而且，个值是独立的；而且，其中有一半是对称的。其中有一半是对称的。科学家科学家Cooley和和Tukey正是巧妙地正是巧妙地利用这些特性加快了利用这些特性加快了DFT的运算速度的运算速度。周期性：周期性：对称性：对称性：Im W86 j W85 W87 -1 1 Re W84 W80 W83 W81 W82 -j0第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.2.2 时域抽取法基时域抽取法基2FFT基本原理基本原理 设序列设序列x(n)的长度的长度N=2M，M为为自然数。自然数。(1)缩短缩短DFT，把把x(n)按按n的奇偶顺序分成两半。的奇偶顺序分成两半。则则x(n)的的DFT为为第第4章章 快速傅里叶变换快速傅里叶变换(FFT)(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)22+(N/2)N2/2次，次，复加复加(N/2)(N/2-1)2+(N/2)2=N2/2次。次。比原来多了还是少了？比原来多了还是少了？(4.2.7)(4.2.8)第第4章章 快速傅里叶变换快速傅里叶变换(FFT)将式将式(4.2.7)和式和式(4.2.8)用流图符号表示，称为蝶形运用流图符号表示，称为蝶形运算符号。算符号。采用蝶形符号可以表示采用蝶形符号可以表示N=8 点的点的DFT运算，下面是经过运算，下面是经过1次分解的次分解的DFT的示意图。的示意图。注意：上半部份有注意：上半部份有4点，用点，用“”的公式做；的公式做；下半部份有下半部份有4点，用点，用“”的公式做。的公式做。第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图图4.2.2 8点点DFT的一次时域抽取分解图的一次时域抽取分解图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)2次分解次分解x(n)的的DFT:(1)缩缩短短x1(r)和和x2(r)的的DFT，与与第第一一次次分分解解相相同同，将将x1(r)按按奇奇偶偶分分解解成成两两个个N/22长长的的子子序序列列x3(l)和和x4(l)，即即则则x1(r)的的DFT为为(4.2.9)第第4章章 快速傅里叶变换快速傅里叶变换(FFT)(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章章 快速傅里叶变换快速傅里叶变换(FFT)图图4.2.3 8点点DFT的第二次时域抽取分解图的第二次时域抽取分解图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)3次分解次分解DFT，长度为，长度为N/23，8点点DITFFT运算流图运算流图需要几次分解需要几次分解DFT，才会使，才会使DFT变为变为1点的点的DFT？第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.2.3时域抽取法快速傅里叶变换的运算量时域抽取法快速傅里叶变换的运算量从分解的级来看从分解的级来看每级需复乘每级需复乘N/2次，次，？复加复加N次；次；？M=log2N级需复乘级需复乘N/2M次，次，？复加复加NM次。次。？对于复乘对于复乘1次需次需50s，复加，复加1次需次需10s的计算机，现在做的计算机，现在做N=1024点的点的DFT运算。运算。按定义直接运算需要按定义直接运算需要 1024250 10-6 102310241010-6=63(s)按按DIT-FFT运算需运算需要要 51210 50 10-6 1024101010-6=0.36(s)它们的速度相差它们的速度相差630.36=175（倍）！（倍）！第第4章章 快速傅里叶变换快速傅里叶变换(FFT)例如：分析序列例如：分析序列x(n)=sin(1.8n)+cos(1.8n)的频谱。的频谱。clear,close all%用两种方法计算用两种方法计算DFTn=0:1023;w=1.8;x=sin(w*n)+cos(w*n);subplot(2,1,1),stem(n,x,.);%axis(250,350,-1.5,1.5)w=linspace(0,2*pi,1024);tic;X1=x*exp(-j*n*w);toc;%时间约时间约1.36秒，复加秒，复加0.2微秒微秒tic;X2=fft(x);toc;%时间约时间约0秒秒subplot(2,1,2),plot(n,abs(X1),.,n,abs(X2),r);%axis(250,350,0,800);%算出角频率算出角频率1.798弧度弧度第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.2.4 DITFFT的运算规律及编程思想的运算规律及编程思想 1.运算规律运算规律 原位计算原位计算从蝶形来看这种运算的好处；从蝶形来看这种运算的好处；有有M级级从每次分解从每次分解DFT次数和次数和DFT变短的规律来看；变短的规律来看；旋转因子旋转因子 ，L指第几级，指第几级，J是序号，从后往前看；是序号，从后往前看；各级蝶形的点距各级蝶形的点距 ，从后往前看。，从后往前看。第第4章章 快速傅里叶变换快速傅里叶变换(FFT)2.编程思想编程思想循环循环1 一级一级地计算蝶形，给出每个蝶的两点距离一级一级地计算蝶形，给出每个蝶的两点距离2L-1；循环循环2一一种种一一种种蝶蝶形形地地计计算算，给给出出旋旋转转因因子子 的的指指数数J，每每级级有有2L-1种不同的蝶；种不同的蝶；循环循环3 同同一一种种蝶蝶里里一一个个一一个个蝶蝶形形地地计计算算，给给出出同同一一种种蝶蝶形形里里各各蝶形的间隔距离蝶形的间隔距离2L。看图说明看图说明第第4章章 快速傅里叶变换快速傅里叶变换(FFT)3.程序框图程序框图图图4.2.6 DITFFT运算程序框图运算程序框图 第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.倒序的意思倒序的意思因为因为DIT-FFT是对是对x(n)的序列按偶奇不断地分解，使的序列按偶奇不断地分解，使得得N=2M的序号按的序号按2倍不断地变短；造成了在蝶形运算时倍不断地变短；造成了在蝶形运算时的输入信号排列顺序与原来的顺序不一样。的输入信号排列顺序与原来的顺序不一样。所以倒序就是从序号的所以倒序就是从序号的2进制的低位向高位不断地把进制的低位向高位不断地把0（代表偶数）和（代表偶数）和1（代表奇数）分开。（代表奇数）分开。图图4.2.7 N=23时的倒序图时的倒序图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)表表4.2.1 顺序和倒序二进制数对照表顺序和倒序二进制数对照表 第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图4.2.8倒序规律第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图4.2.9倒序程序框图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)习题习题1和和2的解的解clear;N=1024;A=N2,N*(N-1);N/2*log2(N),N*log2(N);N*log2(N)+N,2*N*log2(N)b=5e-6,1e-6;T=A*bf=N/T(3)/2第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.2.5 频域抽取法基频域抽取法基2FFT基本原理基本原理 设序列设序列x(n)的长度为的长度为N=2M，M为为自然数。自然数。(1)缩短缩短DFT，将，将x(n)按前后对半分开按前后对半分开。其其DFT可表示为可表示为如下形式：如下形式：第第4章章 快速傅里叶变换快速傅里叶变换(FFT)（2）重组）重组DFT，按，按DFT的定义重新组合变短的的定义重新组合变短的DFT。将。将X(k)分解成偶数组与奇数组，变成分解成偶数组与奇数组，变成N/2点的点的DFT。当当k取偶数时取偶数时当当k取奇数时取奇数时 该运算结构中方括号部份可以用蝶形图表示该运算结构中方括号部份可以用蝶形图表示第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图图4.2.10 DIFFFT蝶形运算流图符号蝶形运算流图符号 采用蝶形符号可以表示采用蝶形符号可以表示N=8 点的点的DFT运算，下面是经过运算，下面是经过1次分解的次分解的DFT的示意图。的示意图。注意：上半部份有注意：上半部份有4点，用点，用“”的公式做；的公式做；下半部份有下半部份有4点，用点，用“”的公式做。的公式做。第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图图4.2.11 N=8的的 DIFFFT一次分解运算流图一次分解运算流图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图图4.2.12 N=8的的 DIFFFT二次分解运算流图二次分解运算流图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图图4.2.13 N=8的的 DIFFFT运算流图运算流图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图图4.2.14 DITFFT的一种变形运算流图的一种变形运算流图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图图4.2.15 DITFFT的的一一种种变变形形运运算流图算流图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.2.6 IDFT的快速算法的快速算法方法方法1：按照按照FFT的方法编造的方法编造IDFT的快速算法的快速算法 程序。程序。那么将产生时域抽取法和频域抽取法两种。那么将产生时域抽取法和频域抽取法两种。好处是？好处是？坏处是？坏处是？方方法法2：利用利用FFT的程序计算的程序计算IDFT。将。将FFT中的中的WNkn改为改为 WN-kn，并且输出乘上，并且输出乘上1/N。比比较较DFT和和IDFT的运的运 算公式就明白：算公式就明白：这么做产生哪两种方法？好处是？坏处是？这么做产生哪两种方法？好处是？坏处是？第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图图4.2.16 是是DITIFFT运算流图运算流图。它是从哪种它是从哪种FFT方法转变过来的？方法转变过来的？为什么称它为为什么称它为DITIFFT运算流图运算流图？第第4章章 快速傅里叶变换快速傅里叶变换(FFT)有时，为了防止运算过程中发生溢出，将有时，为了防止运算过程中发生溢出，将1/N分配分配到每一级的蝶形运算中。下图是到每一级的蝶形运算中。下图是DITIFFT 防止溢出防止溢出的运算流图的运算流图这种做法的乘法次数比前面的增加这种做法的乘法次数比前面的增加 次。次。第第4章章 快速傅里叶变换快速傅里叶变换(FFT)方方法法3：先对：先对X(k)取共轭，取共轭，然后直接然后直接利用利用FFT程序计算程序计算 x*(n)，最后输出再取共轭和乘上，最后输出再取共轭和乘上1/N。怎么知道呢？怎么知道呢？根据是，对某序列两次取共轭等于没有取共轭。根据是，对某序列两次取共轭等于没有取共轭。好处是？坏处是？好处是？坏处是？第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.3 实序列的实序列的FFT算法算法1 直接利用直接利用FFT程序。程序。好处是？坏处是？好处是？坏处是？2 编一个专用于实序列的编一个专用于实序列的FFT程序。程序。好处是？坏处是？好处是？坏处是？3 用一个用一个FFT程序算两个实序列的程序算两个实序列的FFT。根据是根据是DFT的共轭对称性：的共轭对称性：若若 x(n)=x1(n)+jx2(n)，则则 X1(k)=X(k)+X*(N-k)/2 X2(k)=-jX(k)-X*(N-k)/2 好处是？坏处是？好处是？坏处是？第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4 用一个用一个N/2点的点的FFT程序算一个程序算一个N点实序列的点实序列的FFT。将将N点实序列点实序列x(n)分成偶数点和奇数点两个序列分成偶数点和奇数点两个序列x1(r)和和 x2(r)，然后造出一个，然后造出一个N/2点的复序列，点的复序列，y(r)=x1(r)+jx2(r)，r=0,1,(N/2-1)对对y(r)利用利用N/2点的点的FFT程序可以算出程序可以算出 Y(k)=DFTy(r)，k=0,1,(N/2-1)这时，利用方法这时，利用方法3得得 X1(k)=Y(k)+Y*(N/2-k)/2 X2(k)=-jY(k)-Y*(N/2-k)/2 最后，利用时域抽取法快速傅里叶变换的原理得最后，利用时域抽取法快速傅里叶变换的原理得 X(k)=X1(k)+WNkX2(k)，k=0,1,(N-1)好处是？坏处是？好处是？坏处是？第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.4 分裂基分裂基FFT算法算法 从理论上讲，用较大的基数可以进一部减少从理论上讲，用较大的基数可以进一部减少DFT的运的运算次数，但要使程序变得更复杂为代价。算次数，但要使程序变得更复杂为代价。分裂基分裂基FFT算法原理是这样：算法原理是这样：它和基它和基2FFT的基本原理大体相同；的基本原理大体相同；不同的是分裂基不同的是分裂基FFT的做法是把的做法是把N点长的时序分成点长的时序分成4段，段，再按基再按基2FFT的频域抽取法的组合方法把的频域抽取法的组合方法把 这这4段变成段变成1个个(N/2-1)长的长的DFT和和 2个个(N/4-1)长的长的DFT。第第4章章 快速傅里叶变换快速傅里叶变换(FFT)(4.4.2)第第4章章 快速傅里叶变换快速傅里叶变换(FFT)(4.4.3)令第第4章章 快速傅里叶变换快速傅里叶变换(FFT)则(4.4.2)式可写成如下更简明的形式：(4.4.4)第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图4.4.1分裂基第一次分解L形流图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)例如，N=16，第一次抽选分解时，由式(4.4.3)得x2(n)=x(n)+x(n+8),0n7x14(n)=x(n)-x(n+8)-jx(n+4)-x(n+12)Wn16,0n3x24(n)=x(n)-x(n+8)+jx(n+4)-x(n+12)W3n16,0n3把上式代入式(4.4.4)，可得X(2k)=DFTx2(n),0k7X(4k+1)=DFTx14(n),0k3X(4k+3)=DFTx24(n),0k3第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图4.4.2分裂基FFT算法L形排列示意图与结构示意图(a)分裂基FFT算法L形排列示意图；(b)分裂基FFT算法运算流图结构示意图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图4.4.316点分裂基第一次分解L形流图(图中省去箭头)第第4章章 快速傅里叶变换快速傅里叶变换(FFT)第二次分解：先对图4.4.3中N/2点DFT进行分解。令X1(l)=X(2l)，则有X1(2l)=DFTy2(n),0l3X1(4l+1)=DFTy14(n),0l1X1(4l+3)=DFTy24(n),0l1第第4章章 快速傅里叶变换快速傅里叶变换(FFT)其中y2(n)=x2(n)+x2(n+4),0n3y14(n)=x2(n)-x2(n+4)-x2(n+2)x(n+6)Wn8,n=0,1y24(n)=x2(n)-x2(n+4)+jx2(n+2)x2(n+6)W3n8,n=0,1第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图4.4.4图4.4.4中N/2点DFT的分解L形流图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图4.4.54点分裂基L形运算流图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)图4.4.616点分裂基FFT运算流图第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.5 离散哈特莱变换离散哈特莱变换(DHT)4.5.1 离散哈特莱变换定义离散哈特莱变换定义 设设x(n),n=0,1,N-1，为一实序列，其为一实序列，其DHT定义为定义为式中，cas()=cos+sin。其逆变换(IDHT)为(4.5.3)第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.5.2 DHT与与DFT之间的关系之间的关系x(n)的的DFT可表示为可表示为同理，同理，x(n)的的DHT可表示为可表示为XH(k)=ReX(k)-ImX(k)第第4章章 快速傅里叶变换快速傅里叶变换(FFT)4.5.3 DHT的主要优点的主要优点 (1)DHT是是实实值值变变换换，在在对对实实信信号号或或实实数数据据进进行行谱谱分分析析时时避避免免了了复复数数运运算算，从从而而提提高高了了运运算算效效率率，相相应应的硬件也更简单、更经济；的硬件也更简单、更经济；(2)DHT的的正正、反反变变换换(除除因因子子1/N外外)具具有有相相同同的的形形式式，因因而而，实实现现DHT的的硬硬件件或或软软件件既既能能进进行行DHT，也也能进行能进行IDHT；(3)DHT与与DFT间的关系简单，容易实现两种谱之间的关系简单，容易实现两种谱之间的相互转换。间的相互转换。第第4章章 快速傅里叶变换快速傅里叶变换(FFT)clf;%双音频信号的检测双音频信号的检测d=input(type in the telephone digit=,s);symbol=abs(d);tm=49,50,51,65;52,53,54,66;55,56,57,67;42,48,35,68;for p=1:4;for q=1:4;if tm(p,q)=abs(d);break,end end if tm(p,q)=abs(d);break,endendf1=697,770,852,941;f2=1209,1336,1477,1633;figure(1);n=0:204;x=sin(2*pi*n*f1(p)/8000)+sin(2*pi*n*f2(q)/8000);subplot(2,1,1);stem(n,x,.);xlabel(n);X=fft(x);ylabel(双音频信号双音频信号);subplot(2,1,2);stem(n,abs(X),.);xlabel(k,w=2*pi*k/205(弧度弧度),f=8000*k/205(Hz);ylabel(按键对应双音频信号的频谱按键对应双音频信号的频谱);第第4章章 快速傅里叶变换快速傅里叶变换(FFT)k=18,20,22,24,31,34,38,42;val=zeros(1,8);for m=1:8;Fx(m)=X(k(m)+1);endfigure(2);val=abs(Fx);stem(k,val);grid;xlabel(k);limit=80;ylabel(|X(k)|);for s=5:8;if val(s)limit,break,endendfor r=1:4;if val(r)limit,break,endenddisp(按键符号是按键符号是,setstr(tm(r,s-4),)第第4章章 快速傅里叶变换快速傅里叶变换(FFT)Problems1 已知两个序列为已知两个序列为x1=0.95nR64(n)和和x2=sin(0.5n)R48(n)，它们的卷积是它们的卷积是y=x1*x2。求直接计算方法和快速圏积计。求直接计算方法和快速圏积计算方法的运算量。算方法的运算量。2 已知数字振荡器的采样频率是已知数字振荡器的采样频率是2500Hz，它能输出频率，它能输出频率为为150Hz，375Hz，620Hz或或850Hz的正弦波信号。当的正弦波信号。当接收到该振荡器传来的信号时，采用接收到该振荡器传来的信号时，采用DFT检测出信号检测出信号的频率。请确定的频率。请确定DFT的最小点数的最小点数N。