快速计算离散傅里叶变换.ppt

第4章 快速计算离散傅里叶变换,4.1 引言4.2 基2FFT算法4.3 进一步减少运算量的措施4.4 分裂基FFT算法4.5 离散哈特莱变换(DHT),4.1 引言,与序列的傅里叶变换相比，离散傅里叶变换有实用价值。但是按定义直接计算DFT有实用价值吗？只有一些。因为这种算法的计算数度太慢了。特别是与后人发明的算法相比，它的慢更显突出。1965年，J.W.Cooley 和 J.W.Tukey在Mathematics of Computation上发表了An algorithm for the machine calculation of complex Fourier series。它极大的提高了计算离散傅里叶变换的速度。,从定义来看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.2 基2(radix2)FFT算法,4.2.1 直接计算DFT的特点及减少运算量的方法 直接计算N个采样值的DFT 需要有N2次复数乘法和N(N-1)次复数加法。如果把N分成几小段，降低DFT的规模，是不是可以大幅度地减少乘法和加法的运算次数？还有，WNkn具有对称性和周期性，是不是可以巧妙地利用？,例如，当N=8时，从形式上看，W8kn共有64个值。但从图来看，Wkn实际上只有W0W7这8个值是独立的；而且，其中有一半是对称的。科学家Cooley和Tukey正是巧妙地利用这些特性加快了DFT的运算速度。周期性：对称性：,4.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)22+(N/2)N2/2次，复加(N/2)(N/2-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.2 8点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),(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.3 8点DFT的第二次时域抽取分解图,3次分解DFT，长度为N/23，8点DITFFT运算流图需要几次分解DFT，才会使DFT变为1点的DFT？,时域抽取法快速傅里叶变换的运算量从分解的级来看每级需复乘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（倍）！,例如：分析序列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.2.4 DITFFT的运算规律及编程思想 1.运算规律 原位计算从蝶形来看这种运算的好处；有M级从每次分解DFT次数和DFT变短的规律来看；旋转因子，L指第几级，J是序号，从后往前看；各级蝶形的点距，从后往前看。,2.编程思想 循环1 一级一级地计算蝶形，给出每个蝶的两点距离2L-1；循环2一种一种蝶形地计算，给出旋转因子 的指数J，每级有2L-1种不同的蝶；循环3 同一种蝶里一个一个蝶形地计算，给出同一种蝶形里各蝶形的间隔距离2L。看图说明,3.程序框图,图4.2.6 DITFFT运算程序框图,4.倒序的意思 因为DIT-FFT是对x(n)的序列按偶奇不断地分解，使得N=2M的序号按2倍不断地变短；造成了在蝶形运算时的输入信号排列顺序与原来的顺序不一样。所以倒序就是从序号的2进制的低位向高位不断地把0（代表偶数）和1（代表奇数）分开。,图4.2.7 N=23时的倒序图,表4.2.1 顺序和倒序二进制数对照表,图4.2.8 倒序规律,图4.2.9 倒序程序框图,习题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.2.5 频域抽取法基2FFT基本原理 设序列x(n)的长度为N=2M，M为自然数。(1)缩短DFT，将x(n)按前后对半分开。其DFT可表示为如下形式：,（2）重组DFT，按DFT的定义重新组合变短的DFT。将X(k)分解成偶数组与奇数组，变成N/2点的DFT。当k取偶数时当k取奇数时 该运算结构中方括号部份可以用蝶形图表示,图4.2.10 DIFFFT蝶形运算流图符号 采用蝶形符号可以表示N=8 点的DFT运算，下面是经过1次分解的DFT的示意图。注意：上半部份有4点，用“”的公式做；下半部份有4点，用“”的公式做。,图4.2.11 N=8的 DIFFFT一次分解运算流图,图4.2.12 N=8的 DIFFFT二次分解运算流图,图4.2.13 N=8的 DIFFFT运算流图,图4.2.14 DITFFT的一种变形运算流图,图4.2.15 DITFFT的一种变形运算流图,4.2.6 IDFT的快速算法 方法1：按照FFT的方法编造IDFT的快速算法 程序。那么将产生时域抽取法和频域抽取法两种。好处是？坏处是？方法2：利用FFT的程序计算IDFT。将FFT中的WNkn改为 WN-kn，并且输出乘上1/N。比较DFT和IDFT的运 算公式就明白：这么做产生哪两种方法？好处是？坏处是？,图4.2.16 是DITIFFT运算流图。它是从哪种FFT方法转变过来的？为什么称它为DITIFFT运算流图？,有时，为了防止运算过程中发生溢出，将1/N分配到每一级的蝶形运算中。下图是DITIFFT 防止溢出的运算流图这种做法的乘法次数比前面的增加 次。,方法3：先对X(k)取共轭，然后直接利用FFT程序计算 x*(n)，最后输出再取共轭和乘上1/N。怎么知道呢？根据是，对某序列两次取共轭等于没有取共轭。好处是？坏处是？,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 用一个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.4 分裂基FFT算法,从理论上讲，用较大的基数可以进一部减少DFT的运算次数，但要使程序变得更复杂为代价。分裂基FFT算法原理是这样：它和基2FFT的基本原理大体相同；不同的是分裂基FFT的做法是把N点长的时序分成4段，再按基2FFT的频域抽取法的组合方法把 这4段变成1个(N/2-1)长的DFT和 2个(N/4-1)长的DFT。,(4.4.2),(4.4.3),令,则(4.4.2)式可写成如下更简明的形式：,(4.4.4),图4.4.1 分裂基第一次分解L形流图,例如，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),0k7 X(4k+1)=DFTx14(n),0k3 X(4k+3)=DFTx24(n),0k3,图4.4.2 分裂基FFT算法L形排列示意图与结构示意图(a)分裂基FFT算法L形排列示意图；(b)分裂基FFT算法运算流图结构示意图,图4.4.3 16点分裂基第一次分解L形流图(图中省去箭头),第二次分解：先对图中N/2点DFT进行分解。令X1(l)=X(2l)，则有 X1(2l)=DFTy2(n),0l3 X1(4l+1)=DFTy14(n),0l1 X1(4l+3)=DFTy24(n),0l1,其中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.4.4 图中N/2点DFT的分解L形流图,图4.4.5 4点分裂基L形运算流图,图4.4.6 16点分裂基FFT运算流图,4.5 离散哈特莱变换(DHT),4.5.1 离散哈特莱变换定义 设x(n),n=0,1,N-1，为一实序列，其DHT定义为,式中，cas()=cos+sin。其逆变换(IDHT)为,(4.5.3),4.5.2 DHT与DFT之间的关系x(n)的DFT可表示为同理，x(n)的DHT可表示为XH(k)=ReX(k)-ImX(k),4.5.3 DHT的主要优点(1)DHT是实值变换，在对实信号或实数据进行谱分析时避免了复数运算，从而提高了运算效率，相应的硬件也更简单、更经济；(2)DHT的正、反变换(除因子1/N外)具有相同的形式，因而，实现DHT的硬件或软件既能进行DHT，也能进行IDHT；(3)DHT与DFT间的关系简单，容易实现两种谱之间的相互转换。,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(按键对应双音频信号的频谱);,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),),Problems,1 已知两个序列为x1=0.95nR64(n)和x2=sin(0.5n)R48(n)，它们的卷积是y=x1*x2。求直接计算方法和快速圏积计算方法的运算量。2 已知数字振荡器的采样频率是2500Hz，它能输出频率为150Hz，375Hz，620Hz或850Hz的正弦波信号。当接收到该振荡器传来的信号时，采用DFT检测出信号的频率。请确定DFT的最小点数N。,