第4章快速计算离散傅里叶变换.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共需多少时间？ 10242×50 ×10-6 1023×1024×10×10-6=63(s) 2 Cooley和Tukey发明的方法计算DFT需：复乘(N/2)log2N次，复加Nlog2N次。用来计算上面的DFT共需多少时间？ 512×10× 50 ×10-6 1024×10×10×10-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)2×2+(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？,4.2.3 时域抽取法快速傅里叶变换的运算量 从分解的级来看 每级需复乘N/2次，？复加N次；？ M=log2N级需复乘N/2×M次，？复加N×M次。？ 对于复乘1次需50s，复加1次需10s的计算机，现在做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,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弧度,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*b f=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)+j·x2(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)+j·x2(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), 0n7 x14(n)=x(n)-x(n+8)-jx(n+4)-x(n+12)Wn16, 0n3 x24(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形流图(图中省去箭头),第二次分解： 先对图4.4.3中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), 0n3 y14(n)=x2(n)-x2(n+4)-x2(n+2)x(n+6)Wn8,n=0,1 y24(n)=x2(n)-x2(n+4)+jx2(n+2)x2(n+6)W3n8,n=0,1,图4.4.4 图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,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('按键对应双音频信号的频谱');,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),'“'),Problems,1 已知两个序列为x1=0.95nR64(n)和x2=sin(0.5n)R48(n)，它们的卷积是y=x1*x2。求直接计算方法和快速圏积计算方法的运算量。 2 已知数字振荡器的采样频率是2500Hz，它能输出频率为150Hz，375Hz，620Hz或850Hz的正弦波信号。当接收到该振荡器传来的信号时，采用DFT检测出信号的频率。请确定DFT的最小点数N 。,