第4章 快速计算离散傅里叶变换

null第4章 快速计算离散傅里叶变换 第4章 快速计算离散傅里叶变换 4.1 引言 4.2 基2FFT算法 4.3 进一步减少运算量的措施 4.4 分裂基FFT算法 4.5 离散哈特莱变换(DHT) 4.1 引言 4.1 引言 与序列的傅里叶变换相比，离散傅里叶变换有实用价值。但是按定义直接计算DFT有实用价值吗？只有一些。因为这种算法的计算数度太慢了。特别是与后人发明的算法相比，它的慢更显突出。 1965年，J. W. Cooley 和 J. W. Tukey在Mathematics of Computation上发 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 了An algorithm for the machine calculation of complex Fourier series。它极大的提高了计算离散傅里叶变换的速度。null 从定义来看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) 2 Cooley和Tukey发明的 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 计算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 基2(radix2)FFT算法 4.2.1 直接计算DFT的特点及减少运算量的方法 直接计算N个采样值的DFT 需要有N2次复数乘法和N(N-1)次复数加法。 如果把N分成几小段，降低DFT的规模，是不是可以大幅度地减少乘法和加法的运算次数？ 还有，WNkn具有对称性和周期性，是不是可以巧妙地利用？ null 例如，当N=8时，从形式上看，W8kn共有64个值。但从图来看， Wkn实际上只有W0~W7这8个值是独立的；而且，其中有一半是对称的。 科学家Cooley和Tukey正是巧妙地利用这些特性加快了DFT的运算速度。 周期性： 对称性：null 4.2.2 时域抽取法基2FFT基本原理 设序列x(n)的长度N=2M，M为自然数。 (1) 缩短DFT，把x(n)按n的奇偶顺序分成两半。 则x(n)的DFT为null (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)  null将式(4.2.7)和式(4.2.8)用流图符号表示，称为蝶形运算符号。 采用蝶形符号可以表示N=8 点的DFT运算，下面是经过1次分解的DFT的示意图。 注意：上半部份有4点，用“＋”的公式做； 下半部份有4点，用“－”的公式做。null图4.2.2 8点DFT的一次时域抽取分解图null2次分解x(n)的DFT: (1) 缩短x1(r)和x2(r)的DFT，与第一次分解相同，将x1(r)按奇偶分解成两个N/22长的子序列x3(l)和x4(l)，即 则x1(r)的DFT为(4.2.9) null(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？ null图4.2.3 8点DFT的第二次时域抽取分解图null3次分解DFT，… ，长度为N/23，… 8点DIT―FFT运算流图 需要几次分解DFT，才会使DFT变为1点的DFT？null 4.2.3 时域抽取法快速傅里叶变换的运算量 从分解的级来看—— ①每级需复乘N/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 （倍）！ null例如： 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 序列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弧度null 4.2.4 DIT―FFT的运算规律及编程思想 1. 运算规律 原位计算——从蝶形来看这种运算的好处； 有M级——从每次分解DFT次数和DFT变短的规律来看； 旋转因子 ，L指第几级，J是序号，从后往前看； 各级蝶形的点距 ，从后往前看。 null2. 编程思想 循环1—— 一级一级地计算蝶形，给出每个蝶的两点距离2L-1； 循环2—— 一种一种蝶形地计算，给出旋转因子 的指数J，每级有2L-1种不同的蝶； 循环3 —— 同一种蝶里一个一个蝶形地计算，给出同一种蝶形里各蝶形的间隔距离2L。 看图说明null 3. 程序框图 图4.2.6 DIT―FFT运算程序框图 null 4. 倒序的意思 因为DIT-FFT是对x(n)的序列按偶奇不断地分解，使得N=2M的序号按2倍不断地变短；造成了在蝶形运算时的输入信号排列顺序与原来的顺序不一样。 所以倒序就是从序号的2进制的低位向高位不断地把0（代表偶数）和1（代表奇数）分开。图4.2.7 N=23时的倒序图null表4.2.1 顺序和倒序二进制数对照表 null 图4.2.8 倒序规律 null 图4.2.9 倒序程序框图 习 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 1和2的解习题1和2的解clear; N=1024; A=[N^2,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)/2null 4.2.5 频域抽取法基2FFT基本原理 设序列x(n)的长度为N=2M，M为自然数。 (1) 缩短DFT，将x(n)按前后对半分开。其DFT可表示为如下形式：null（2）重组DFT，按DFT的定义重新组合变短的DFT。将X(k)分解成偶数组与奇数组，变成N/2点的DFT。 当k取偶数时 当k取奇数时 该运算结构中方括号部份可以用蝶形图表示null 图4.2.10 DIF―FFT蝶形运算流图符号 采用蝶形符号可以表示N=8 点的DFT运算，下面是经过1次分解的DFT的示意图。 注意：上半部份有4点，用“＋”的公式做； 下半部份有4点，用“－”的公式做。null图4.2.11 N=8的 DIF―FFT一次分解运算流图null图4.2.12 N=8的 DIF―FFT二次分解运算流图null图4.2.13 N=8的 DIF―FFT运算流图null图4.2.14 DIT―FFT的一种变形运算流图null图4.2.15 DIT―FFT的一种变形运算流图null 4.2.6 IDFT的快速算法 方法1： 按照FFT的方法编造IDFT的快速算法 程序。 那么将产生时域抽取法和频域抽取法两种。 好处是？ 坏处是？ 方法2：利用FFT的程序计算IDFT。将FFT中的WNkn改为 WN-kn，并且输出乘上1/N。比较DFT和IDFT的运 算公式就明白： 这么做产生哪两种方法？好处是？坏处是？null图4.2.16 是DIT―IFFT运算流图 。 它是从哪种FFT方法转变过来的？ 为什么称它为DIT―IFFT运算流图 ？null 有时，为了防止运算过程中发生溢出，将1/N分配到每一级的蝶形运算中。下图是DIT―IFFT 防止溢出的运算流图 这种做法的乘法次数比前面的增加 次。null方法3：先对X(k)取共轭，然后直接利用FFT程序计算 x*(n)，最后输出再取共轭和乘上1/N。 怎么知道呢？ 根据是，对某序列两次取共轭等于没有取共轭。 好处是？坏处是？4.3 实序列的FFT算法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)= -j[X(k)-X*(N-k)]/2 好处是？坏处是？ null 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)=DFT[y(r)]，k=0, 1, … , (N/2-1) 这时，利用方法3得 X1(k)=[Y(k)+Y*(N/2-k)]/2 X2(k)= -j[Y(k)-Y*(N/2-k)]/2 最后，利用时域抽取法快速傅里叶变换的原理得 X(k)=X1(k) +WNkX2(k)，k=0, 1, … , (N-1) 好处是？坏处是？ 4.4 分裂基FFT算法 4.4 分裂基FFT算法 从理论上讲，用较大的基数可以进一部减少DFT的运 算次数，但要使程序变得更复杂为代价。 分裂基FFT算法原理是这样： 它和基2FFT的基本原理大体相同； 不同的是分裂基FFT的做法是把N点长的时序分成4段， 再按基2FFT的频域抽取法的组合方法把 这4段变成1个(N/2-1)长的DFT和 2个(N/4-1)长的DFT。null(4.4.2) null(4.4.3) 令 null 则(4.4.2)式可写成如下更简明的形式： (4.4.4) null图4.4.1 分裂基第一次分解L形流图 null例如，N=16，第一次抽选分解时，由式(4.4.3)得  x2(n)=x(n)+x(n+8), 0≤n≤7 x14(n)={［x(n)-x(n+8)］-j［x(n+4)-x(n+12)］}Wn16, 0≤n≤3 x24(n)={［x(n)-x(n+8)］+j［x(n+4)-x(n+12)］}W3n16, 0≤n≤3  把上式代入式(4.4.4)，可得 X(2k)=DFT［x2(n)］, 0≤k≤7 X(4k+1)=DFT［x14(n)］, 0≤k≤3 X(4k+3)=DFT［x24(n)］, 0≤k≤3 null图4.4.2 分裂基FFT算法L形排列示意图与结构示意图 (a)分裂基FFT算法L形排列示意图； (b)分裂基FFT算法运算流图结构示意图 null图4.4.3 16点分裂基第一次分解L形流图(图中省去箭头) null 第二次分解： 先对图4.4.3中N/2点DFT进行分解。 令X1(l)=X(2l)，则有 X1(2l)=DFT［y2(n)］, 0≤l≤3 X1(4l+1)=DFT［y14(n)］, 0≤l≤1 X1(4l+3)=DFT［y24(n)］, 0≤l≤1 null 其中 y2(n)=x2(n)+x2(n+4), 0≤n≤3 y14(n)={［x2(n)-x2(n+4)］-［x2(n+2)x(n+6)］}Wn8,n=0,1 y24(n)={［x2(n)-x2(n+4)］+j［x2(n+2)x2(n+6)］}W3n8,n=0,1 null图4.4.4 图4.4.4中N/2点DFT的分解L形流图 null 图4.4.5 4点分裂基L形运算流图 null 图4.4.6 16点分裂基FFT运算流图 4.5 离散哈特莱变换(DHT) 4.5 离散哈特莱变换(DHT) 4.5.1 离散哈特莱变换定义 设x(n),n=0,1,…,N-1，为一实序列，其DHT定义为式中，cas(α)=cosα+sinα。其逆变换(IDHT)为 (4.5.3) null 4.5.2 DHT与DFT之间的关系 x(n)的DFT可表示为 同理，x(n)的DHT可表示为 XH(k)=Re[X(k)]-Im[X(k)] null 4.5.3 DHT的主要优点 (1)DHT是实值变换，在对实信号或实数据进行谱分析时避免了复数运算，从而提高了运算效率，相应的硬件也更简单、更经济； (2)DHT的正、反变换(除因子1/N外)具有相同的形式，因而，实现DHT的硬件或软件既能进行DHT，也能进行IDHT； (3)DHT与DFT间的关系简单，容易实现两种谱之间的相互转换。 nullclf;%双音频信号的检测 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('按键对应双音频信号的频谱');nullk=[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)),'"'])ProblemsProblems1 已知两个序列为x1=0.95nR64(n)和x2=sin(0.5n)R48(n)，它们的卷积是y=x1*x2。求直接计算方法和快速圏积计算方法的运算量。 2 已知数字振荡器的采样频率是2500Hz，它能输出频率为150Hz，375Hz，620Hz或850Hz的正弦波信号。当接收到该振荡器传来的信号时，采用DFT检测出信号的频率。请确定DFT的最小点数N 。