快速计算离散傅里叶变换

,第二级,第三级,第四级,第五级,*,第4章 快速傅里叶变换(FFT),第4章 快速计算离散傅里叶变换,4.1 引言,4.2 基2FFT算法,4.3 进一步减少运算量的措施,4.4 分裂基FFT算法,4.5 离散哈特莱变换(DHT),1,4.1 引言,与序列的傅里叶变换相比，离散傅里叶变换有实用价值。但是按定义直接计算DFT有实用价值吗？只有一些。因为这种算法的计算数度太慢了。特别是与后人发明的算法相比，它的慢更显突出。,1965年，J. W. Cooley 和 J. W. Tukey在Mathematics of Computation上发表了An algorithm for the machine calculation of complex Fourier series。它极大的提高了计算离散傅里叶变换的速度。,2,从定义来看N点长的DFT的运算量。,1 直接计算DFT需：复乘N,2,次，复加(N-1)N次。因为 1个k需复乘N次，复加(N-1)次。,对于复乘1次需50s，复加1次需10s的计算机，用直接法做N=1024点长的DFT共需多少时间？,1024,2,50 10,-6,102310241010,-6,=63(s),2 Cooley和Tukey发明的方法计算DFT需：复乘(N/2)log,2,N次，复加Nlog,2,N次。用来计算上面的DFT共需多少时间？,51210 50 10,-6,1024101010,-6,=0.36(s),3,4.2 基2(radix2)FFT算法,4.2.1 直接计算DFT的特点及减少运算量的方法,直接计算N个采样值的DFT 需要有N,2,次复数乘法和N(N-1)次复数加法。,如果把N分成几小段，降低DFT的规模，是不是可以大幅度地减少乘法和加法的运算次数？,还有，W,N,kn,具有对称性和周期性，是不是可以巧妙地利用？,4,例如，当N=8时，从形式上看，W,8,kn,共有64个值。但从图来看， W,kn,实际上只有W,0,W,7,这8个值是独立的；而且，其中有一半是对称的。 科学家Cooley和Tukey正是巧妙地利用这些特性加快了DFT的运算速度,。,周期性：,对称性：,Im,W,8,6,j,W,8,5,W,8,7,-1 1 Re,W,8,4,W,8,0,W,8,3,W,8,1,W,8,2,-j,0,5,4.2.2 时域抽取法基2FFT基本原理,设序列x(n)的长度N=2,M,，M为自然数。,(1) 缩短DFT，把x(n)按n的奇偶顺序分成两半。,则x(n)的DFT为,6,(2) 重组DFT，按DFT的定义重新组合变短的DFT。变短后的DFT中X,1,(k)和X,2,(k)分别为x,1,(r)和x,2,(r)的N/2点DFT，周期为N/2；对称性W,N,k+N/2,= W,N,k,。X(k)又可表示为,经过这两步骤处理后，1个N点的DFT就变成了2个N/2点的DFT。运算量变成：,复乘(N/2),2,2+(N/2)N,2,/,2次，,复加(N/2) (N/2-1),2 +(N/2) 2=N,2,/,2次。,比原来多了还是少了？,(4.2.7),(4.2.8),7,将式(4.2.7)和式(4.2.8)用流图符号表示，称为蝶形运算符号。,采用蝶形符号可以表示N=8 点的DFT运算，下面是经过1次分解的DFT的示意图。,注意：上半部份有4点，用“”的公式做；,下半部份有4点，用“”的公式做。,8,图4.2.2 8点DFT的一次时域抽取分解图,9,2次分解x(n)的DFT:,(1) 缩短x,1,(r)和x,2,(r)的DFT，与第一次分解相同，将x,1,(r)按奇偶分解成两个N/2,2,长的子序列x,3,(l)和x,4,(l)，即,则x,1,(r)的DFT为,(4.2.9),10,(2)重组DFT，按DFT的定义重新组合变短的DFT。变短后的DFT中X,3,(k)和X,4,(k)分别为x,3,(l)和x,4,(l)的N/4点DFT，周期为N/2,2,；对称性W,N/2,k+N/4,= W,N/2,k,。X,1,(k)也可表示为,用同样的方法可以计算出,如果是8点的DFT，经两次分解，DFT的长度是多少,？,有几个这种长度的DFT,？,11,图4.2.3 8点DFT的第二次时域抽取分解图,12,3,次分解DFT， ，长度为N/2,3,，,8点DITFFT运算流图,需要几次分解DFT，才会使DFT变为1点的DFT,？,13,4.2.3,时域抽取法快速傅里叶变换的运算量,从分解的级来看,每级需复乘N/2次，,？,复加N次；,？,M=log,2,N级需复乘N/2M次，,？,复加NM次。,？,对于复乘1次需50s，复加1次需10s的计算机，现在做N=1024点的DFT运算。,按定义直接运算需要,1024,2,50 10,-6,102310241010,-6,=63(s),按DIT-FFT运算需,要,51210 50 10,-6,1024101010,-6,=0.36(s),它们的速度相差630.36=175 （倍）！,14,例如：分析序列x(n)=sin(1.8n)+cos(1.8n)的频谱。,clear,close all %用两种方法计算DFT,n=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弧度,15,4.2.4 DITFFT的运算规律及编程思想,1. 运算规律,原位计算从蝶形来看这种运算的好处；,有M级从每次分解DFT次数和DFT变短的规律来看；,旋转因子 ，L指第几级，J是序号，从后往前看；,各级蝶形的点距 ，从后往前看。,16,2. 编程思想,循环1,一级一级地计算蝶形，给出每个蝶的两点距离2,L-1,；,循环2,一种一种蝶形地计算，给出旋转因子 的指数J，每级有2,L-1,种不同的蝶；,循环3 ,同一种蝶里一个一个蝶形地计算，给出同一种蝶形里各蝶形的间隔距离2,L,。,看图说明,17,3. 程序框图,图4.2.6 DITFFT运算程序框图,18,4,. 倒序的意思,因为DIT-FFT是对x(n)的序列按偶奇不断地分解，使得N=2,M,的序号按2倍不断地变短；造成了在蝶形运算时的输入信号排列顺序与原来的顺序不一样。,所以倒序就是从序号的2进制的低位向高位不断地把0（代表偶数）和1（代表奇数）分开。,图4.2.7 N=2,3,时的倒序图,19,表4.2.1 顺序和倒序二进制数对照表,20,图4.2.8 倒序规律,21,图4.2.9 倒序程序框图,22,习题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*b,f=N/T(3)/2,23,4.2.5 频域抽取法基2FFT基本原理,设序列x(n)的长度为N=2,M,，M为自然数。,(1) 缩短DFT，将x(n)按前后对半分开。其DFT可表示为如下形式：,24,（2）重组DFT，按DFT的定义重新组合变短的DFT。将X(k)分解成偶数组与奇数组，变成N/2点的DFT。,当k取偶数时,当,k取奇数时,该运算结构中方括号部份可以用蝶形图表示,25,图4.2.10 DIFFFT蝶形运算流图符号,采用蝶形符号可以表示N=8 点的DFT运算，下面是经过1次分解的DFT的示意图。,注意：上半部份有4点，用“”的公式做；,下半部份有4点，用“”的公式做。,26,图4.2.11 N=8的 DIFFFT一次分解运算流图,27,图4.2.12 N=8的 DIFFFT二次分解运算流图,28,图4.2.13 N=8的 DIFFFT运算流图,29,图4.2.14 DITFFT的一种变形运算流图,30,图4.2.15 DITFFT的一种变形运算流图,31,4.2.6 IDFT的快速算法,方法1： 按照FFT的方法编造IDFT的快速算法 程序。,那么将产生时域抽取法和频域抽取法两种。,好处是？ 坏处是？,方法2：,利用FFT的程序计算IDFT。将FFT中的W,N,kn,改为,W,N,-kn,，并且输出乘上1/N。,比,较DFT和IDFT的运,算公式就明白：,这么做产生哪两种方法？好处是？坏处是？,32,图4.2.16 是DITIFFT运算流图 。,它是从哪种FFT方法转变过来的？,为什么称它为DITIFFT运算流图 ？,33,有时，为了防止运算过程中发生溢出，将1/N分配到每一级的蝶形运算中。下图是DITIFFT 防止溢出的运算流图,这种做法的乘法次数比前面的增加 次。,34,方法3：先对,X(k),取共轭，然后直接利用FFT程序计算,x,*,(n)，最后输出再取共轭和乘上1/N。,怎么知道呢？,根据是，对某序列两次取共轭等于没有取共轭。,好处是？坏处是？,35,4.3 实序列的FFT算法,1 直接利用FFT程序。,好处是？坏处是？,2 编一个专用于实序列的FFT程序。,好处是？坏处是？,3 用一个FFT程序算两个实序列的FFT。,根据是DFT的共轭对称性：,若 x(n)=x,1,(n)+jx,2,(n)，,则 X,1,(k)=X(k)+X,*,(N-k)/2,X,2,(k)= -jX(k)-X,*,(N-k)/2,好处是？坏处是？,36,4 用一个N/2点的FFT程序算一个N点实序列的FFT 。,将N点实序列x(n)分成偶数点和奇数点两个序列x,1,(r)和,x,2,(r) ，然后造出一个N/2点的复序列，,y(r)= x,1,(r)+jx,2,(r) ，r=0, 1, , (N/2-1),对y(r)利用N/2点的FFT程序可以算出,Y(k)=DFTy(r)，k=0, 1, , (N/2-1),这时，利用方法3得,X,1,(k)=Y(k)+Y,*,(N/2-k)/2,X,2,(k)= -jY(k)-Y,*,(N/2-k)/2,最后，利用时域抽取法快速傅里叶变换的原理得,X(k)=X,1,(k) +W,N,k,X,2,(k)，k=0, 1, , (N-1),好处是？坏处是？,37,4.4 分裂基FFT算法,从理论上讲，用较大的基数可以进一部减少DFT的运,算次数，但要使程序变得更复杂为代价。,分裂基FFT算法原理是这样：,它和基2FFT的基本原理大体相同；,不同的是分裂基FFT的做法是把N点长的时序分成4段，,再按基2FFT的频域抽取法的组合方法把 这4段变成1个(N/2-1)长的DFT和 2个(N/4-1)长的DFT。,38,(4.4.2),39,(4.4.3),令,40,则(4.4.2)式可写成如下更简明的形式：,(4.4.4),41,图4.4.1 分裂基第一次分解L形流图,42,例如，N=16，第一次抽选分解时，由式(4.4.3)得,x,2,(n)=x(n)+x(n+8), 0n7,x,14,(n)=x(n)-x(n+8)-jx(n+4)-x(n+12)W,n,16,0n3,x,24,(n)=x(n)-x(n+8)+jx(n+4)-x(n+12)W,3n,16,0n3,把上式代入式(4.4.4)，可得,X(2k)=DFTx,2,(n), 0k7,X(4k+1)=DFTx,14,(n), 0k3,X(4k+3)=DFTx,24,(n), 0k3,43,图4.4.2 分裂基FFT算法L形排列示意图与结构示意图,(a)分裂基FFT算法L形排列示意图；,(b)分裂基FFT算法运算流图结构示意图,44,图4.4.3 16点分裂基第一次分解L形流图(图中省去箭头),45,第二次分解：,先对图4.4.3中N/2点DFT进行分解。,令X,1,(l)=X(2l)，则有,X,1,(2l)=DFTy2(n), 0l3,X,1,(4l+1)=DFTy,1,4,(n), 0l1,X,1,(4l+3)=DFTy,2,4,(n), 0l1,46,其中,y,2,(n)=x,2,(n)+x,2,(n+4), 0n3,y,1,4,(n)=x,2,(n)-x,2,(n+4)-x,2,(n+2)x(n+6)W,n,8,n=0,1,y,2,4,(n)=x,2,(n)-x,2,(n+4)+jx,2,(n+2)x,2,(n+6)W,3,n,8,n=0,1,47,图4.4.4 图4.4.4中N/2点DFT的分解L形流图,48,图4.4.5 4点分裂基L形运算流图,49,图4.4.6 16点分裂基FFT运算流图,50,4.5 离散哈特莱变换(DHT),4.5.1 离散哈特莱变换定义,设x(n),n=0,1,N-1，为一实序列，其DHT定义为,式中，cas()=cos+sin。其逆变换(IDHT)为,(4.5.3),51,4.5.2 DHT与DFT之间的关系,x(n)的DFT可表示为,同理，x(n)的DHT可表示为,X,H,(k)=ReX(k)-ImX(k),52,4.5.3 DHT的主要优点,(1)DHT是实值变换，在对实信号或实数据进行谱分析时避免了复数运算，从而提高了运算效率，相应的硬件也更简单、更经济；,(2)DHT的正、反变换(除因子1/N外)具有相同的形式，因而，实现DHT的硬件或软件既能进行DHT，也能进行IDHT；,(3)DHT与DFT间的关系简单，容易实现两种谱之间的相互转换。,53,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,end,end,f1=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(按键对应双音频信号的频谱);,54,k=18,20,22,24,31,34,38,42;,val=zeros(1,8);,for m=1:8;,Fx(m)=X(k(m)+1);,end,figure(2);val=abs(Fx);,stem(k,val);grid;xlabel(k);,limit=80; ylabel(|X(k)|);,for s=5:8;if val(s)limit,break,end,end,for r=1:4;if val(r)limit,break,end,end,disp(按键符号是,setstr(tm(r,s-4),),55,Problems,1 已知两个序列为x,1,=0.95,n,R,64,(n)和x,2,=sin(0.5n)R,48,(n)，它们的卷积是y=x,1,*x,2,。求直接计算方法和快速圏积计算方法的运算量。,2 已知数字振荡器的采样频率是2500Hz，它能输出频率为150Hz，375Hz，620Hz或850Hz的正弦波信号。当接收到该振荡器传来的信号时，采用DFT检测出信号的频率。请确定DFT的最小点数N 。,56,