第3章.离散傅里叶变换及其快速算法与应用

1、n离散傅里叶变换（DFT)的定义nDFT与DTFT和Z变换的关系nDFT的数字计算nDFT的性质nDFT 的快速算法FFTnDFT(FFT)的应用第第3 3章章 离散傅里叶变换及快速算法与应用离散傅里叶变换及快速算法与应用 傅里叶变换对的数据如果连续或无限，则不适合在计算机上运算。从数值计算的角度出发，我们感兴趣的是时域及频域都离散且有限的情况，这就是我们本章要研究的离散傅里叶变换。离散傅里叶变换(Discrete Fourier Transform，简称，简称DFT)是数字信号处理中非常重要的一种数学变换，除了在理论上相当重要之外，因其存在有效的快速算法，故在各种数字信号处理的算法中起核心作

2、用。引言引言3.1.1 离散傅里叶变换的定义离散傅里叶变换的定义10210(0 -1)(:( ) ( )( )0,1,11( )( )( )0,1,1NknNnjNnNknNkx nM nMx nNDFTX kDFT x nx n WkNWeNMx nIDFT X kX k WnkNN）是长为 （的有限长序列，则定义）的 点其中为旋转因子，离散傅里叶逆变换IDFT：3.1.2 DFT与与DTFT和和Z变换的关系变换的关系210101202( ),0,1,-1( ) ( )( )() ( )( )( ) ( )( )()|( )|20,1,1,jkNMnnMjj nnMjknNNnjkz eNx

3、 nM nMX zZT x nx n zX eDTFT x nx n eX kDFT x nx n eX eX zkNkN长表示频率：结论：结论： DFT是有限长序列x(n)的傅里叶变换在区间0，2上的N点等间隔采样；DFT是对有限长序列x(n) 进行Z变换后在单位圆上进行等间隔采样的采样值 。3.1.2 DFT与与DTFT和和Z变换的关系变换的关系kN223.1.4 DFT的数字计算的数字计算（直接定义）（直接定义）functionX=dft(x,N)dw=2*pi/N;n=0:N-1;k=0:N-1;M=length(x);if NM for l=M+1:N x(l)=0; endendX

4、=x*(exp(-j*dw).(n*k);DFT计算子函数：IDFT计算子函数：3.1.4 DFT的数字计算的数字计算（直接计算）（直接计算）functionx=idft(X,N)dw=2*pi/N;n=0:N-1;k=0:N-1;x=X*(exp(j*dw).(n*k)/N;3.1.4 DFT的数字计算的数字计算（直接计算）（直接计算）41 ( )( ),( )48.x nR nx nDFT例 求的 点和 点x=1 1 1 1 ;n=0:3;k=-1500:1500;dw=(pi/500);X,w=dift(x,n,dw,k);magX=abs(X);figure(1)plot(w/pi,m

5、agX);N=4;X=dft(x,N);MX=abs(X);figure(2)stem(MX);x1=idft(X,N);figure(3)stem(real(x1)n1）隐含的周期性n 2）循环移位性质3.2 DFT的的性质性质( )( ).knNWX kx nN因的周期性，和隐含周期性，且周期均为11( ),( )()( )( ),011,( ):mod( ,)NNNNx nMy nx nmR nnNnnMNnnNMnnmatlabn N长为则其循环移位定义：表示模 对 求余：为整数循环移位的实现子程序：function y=cirshift(x,m,N)x=x zeros(1,N-len

6、gth(x);n=0:N-1;n=mod(n-m,N);y=x(n+1);152( )10 (0.8) ,010,(6)nx nnx n例 已知求n=0:10;x=10*(0.8).n;y=cirshift(x,6,15);n=0:14;x=x,zeros(1,4);subplot(2,1,1)stem(n,x);subplot(2,1,2)stem(n,y);024681012140510 xn024681012140246810 x(n-6)3.2 DFT的的性质性质n2）循环移位性质时域循环移位定理频率循环移位定理( )()( )( )( )NNkmNy nx nmR nY kWX k则

7、( )()( )( )( )NNnlNY kXklR ky nW x n则n3）循环卷积定理1.两个有限长序列的循环卷积3.2 DFT的的性质性质10( )( ). ( )( )( )( ) () ( )= ( )( )max,LcLLmh nx nNM h nx nLy nh m x nmR nh nx nLN M设与的长度分别为 和与的的循环卷积：循环卷积的计算：function y=circonv(x1,x2,L)x1=x1 zeros(1,L-length(x1);x2=x2 zeros(1,L-length(x2);n=0:L-1;ix2=x2(mod(-n,L)+1);for k=

8、0:L-1 y1=cirshift(ix2,k,L); yc(k+1)=x1*y1;endy=yc;123( )1,2,2 ,( )1,2,3,4 ,04.x nx nn例 已知求其4点、5点和6点循环卷积。x1=1 2 2;x2=1 2 3 4;y4=circonv(x1,x2,4);y5=circonv(x1,x2,5);y6=circonv(x1,x2,6);y4 = 15 12 9 14y5 = 9 4 9 14 14y6 = 1 4 9 14 14 8n时域循环卷积定理n频域循环卷积定理3.2 DFT的的性质性质1212122121( )( )max,( )( )( ),( )( )

9、( )( )x nx nNNNN Nx nx nx nx nNDFTX kXkX k和的长度分别为和，则的 点为：1212( )( )( ),1( )( )( )x nx n x nX kX kXkN如果则循环卷积定理验证12124( )1,2,2,0 ,( )1,2,3,4 ,04.4( )( )( )( ),x nx nny nX kXkY k例 已知求其 点循环卷积、并验证循环卷积定理。x1=1 2 2 0;x2=1 2 3 4;N=4;y=circonv(x1,x2,N); X1=dft(x1,N);X2=dft(x2,N);Y=dft(y,N); yi=idft(X1.*X2,N);

10、%yi=y?nDFT的对称性0 -1,/ 2NN因DFT定义区间为这里的对称是指关于的对称性。( )( ) ( ) ,( )(- )-|( )| |(- )|-arg( )-arg(- )-Nx nX kDFT x nXkX N kX kX N kX kX N k若是实序列，则共轭对称模复角15( )1,2,3,4 ,010.x nn例 其DFT:X1 = 10 , -2 + 2i , -2 , -2 - 2i3.3 频率域采样频率域采样( ),(),0 2(),)( )( )( )()( ),( )( )jjNNiNx nMDTFTX eX eNX kX kx nNxnIDFT X kx n

11、iN R nNM xnx n如果序列的长度为其为在对等间隔抽样 个点得到抽样序列（能由（重建吗？ 应满足何条件？110( )( )11( )( )1NNkkNNMX kX zzX zX kNW z内插公式：时，可由从建例 3.3.1M=26;n=0:M;xa=0:M/2;xb=ceil(M/2)-1:-1:0;xn=xa,xb;N=27;k=0:N-1;dw=2*pi/N;X,w=dift(xn,n,dw,k);figure(1)subplot(1,2,1)stem(n,xn);xlabel(xn);x1=idft(X,N);subplot(1,2,2)stem(k,real(x1)xlabe

12、l(x1);010203002468101214xn010203002468101214x1N=27010203002468101214xn05101502468101214x1N=16010203002468101214xn01020304002468101214x1N=32 直接直接DFT存在的问题存在的问题结论：结论：直接DFT运算共需N N2 2次复数乘法次复数乘法, N(N-1)N(N-1)次复数加法次复数加法。由于一次复乘需4次实乘和2次实加，一次复加需2次实加， 因此直接DFT运算共需4N2次实乘及2N(2N-1)次实加。当N很大时运算量非常大，如果要进行实时信号处理就必须设法减

13、少运算量。 1,.,1 , 0 )()(DFT)(10NkWnxnxkXNnknN4.DFT 的快速算法的快速算法FFT 减少减少DFT运算量的方法运算量的方法思路：思路：将长序列分解为短序列，利用DFT系数的对称性和周期性，分解计算再合并处理，可大大减少运算量。 方法：方法：根据分解与合并方法的不同，有多种快速算法。基本算法可分为两类： 按时间抽取法（Decimation in Time，缩写为DIT） 按频率抽取法（Decimation in Frequency，缩写为DIF）时间抽取法：时间抽取法：将序列按时间顺序的奇偶分解为越来越短的子序列进行计算的快速算法，也称Cooley-Tuke

14、y算法。频率抽取法：频率抽取法：将序列按时间顺序的前后分解为越来越短的子序列进行计算的快速算法， 也称Sande-Tukey算法。 FFT( Fast Fourier Transform )并不是一种新的傅里叶变换。FFT是计算离散傅里叶变换（DFT）的快速算法。重点以按时间抽取的基2FFT(DIT-FFT)算法为例，介绍FFT算法的基本原理，归纳FFT算法的运算规律，总结FFT算法的编程框图，列出FFT算法的实现程序，从而建立起一种从算法理论到程序实现全过程的完整概念。 另外，简单介绍进一步减少运算量的措施和其它快速算法。 快速傅里叶变换快速傅里叶变换将x(n)按 n 的奇偶分为两组：)()

15、()()() 12()2()()()()(211202212021120)12(12021210kXWkXWrxWWrxWrxWrxWnxWnxWnxkXkNNrrkNkNNrrkNNrkrNNrrkNNnnkNNnnkNNnnkN为奇为偶 基基2DIT-FFT2DIT-FFT算法的算法的算法原理算法原理 DIT蝶形运算流图符号)()()(21kXWkXkXkN)()()2()2()2(212)2(1kXWkXNkXWNkXNkXkNNkN说明：式中得到的X(k)只有N/2个点，实际上X(k)有N个点，要用X1(k)和X2(k)表示全部的X(k)值，还要应用系数的周期性和对称性。 基基2DIT

16、-FFT2DIT-FFT算法的算法的运算流程运算流程 按时间抽取将一个N点DFT分解为两个N/2点DFT（N=8）8点DFT二次时域抽取分解运算流图N=8 点DIT-FFT运算流图例：计算x(n)=0,1,2,3,4,5,6,7的DFT值X(k) 基基2DIT-FFT2DIT-FFT算法的算法的运算规律运算规律 蝶形运算：蝶形运算：蝶形运算是分级进行的；每级的蝶形运算可以按旋转因子的指数大小排序进行；如果指数大小一样则可从上往下依次蝶算。蝶形运算的规律是编程的基础 。原位计算：原位计算：当数据输入到存储器以后，每一级运算的结果仍可储存在同一组存储器中，直到最后输出，中间无需其它存储器，这叫原位

17、计算。 序列倒序序列倒序 ：运算前要先对输入的序列进行位序颠倒 。顺序数（I） 顺序二进制（IB） 倒序二进制（JB） 倒序数 （J） 0123456700000101001110010111011100010001011000110101111104261537规律：规律：倒序二进制数是顺序二进制数的倒置。 基基2DIT-FFT2DIT-FFT算法的算法的程序框图程序框图 基基2DIT-FFT2DIT-FFT算法的实现算法的实现程序程序 clc;close all;clear;format compact;%倒序xn=0,1,2,3,4,5,6,7;%可取任意序列M=nextpow2(len

18、gth(xn), N=2M,%将数据输入到存储地址A=xn,zeros(1,N-length(xn); disp(输入到各存储单元的数据:),disp(A);%数据倒序操作J=0;%给倒序数赋初值for I=0:N-1;%按顺序交换数据和算倒序数 if I=K; J=J-K;K=K/2; end J=J+K;enddisp(倒序后各存储单元的数据:),disp(A);%分级按序依次进行蝶形运算WN=exp(-j*2*pi/N);%计算常量for L=1:M;%分级计算 disp(运算级次:),disp(L); B=2(L-1); for R=0:B-1;%各级按序计算 P=2(M-L)*R;

19、for K=R:2L:N-2;%每序依次计算 T=A(K+1)+A(K+B+1)*WNP; A(K+B+1)=A(K+1)-A(K+B+1)*WNP; A(K+1)=T; end end disp(本级运算后各存储单元的数据:),disp(A); enddisp(输出各存储单元的数据:),Xk=A,disp(调用fft函数运算的结果:),fftxn=fft(xn,N), 基基2DIT-FFT2DIT-FFT算法的运算量分析算法的运算量分析 对于任何一个2的整数次幂 N=2M 点的DFT运算，总可以通过M 次分解，直到DFT 就是其本身的1点序列。这样的M次分解，构成从x(n)到X(k)的M 级

20、运算。从流程图可看到，每级运算都由N/2个蝶形运算构成。因此每级运算都需要 (N/2) 次复乘和N次复加(每个结作加、减各一次)。故总共需要的运算： 复乘 (N/2)M=(N/2)log2N 复加 NM=N log2N 实际运算量与这个数字稍有出入，但按时间抽取法所需的计算量，大概与N log2N 成正比，而直接运算的运算量与N2成正比，当N较大时，FFT算法优势明显。 IDFT IDFT的快速算法的快速算法方法方法1 1：只要把DFT运算中的蝶形系数W nkN 改为W -nkN 。结果再除N，则FFT运算程序可直接进行IDFT运算。FFT算法可用于IDFT运算，简称为IFFT。方法方法2 2

21、：不改FFT程序，直接用它作IFFT(1)将X（k）取共轭（虚部乘以-1） (2)然后直接作FFT运算(3)最后对FFT结果取共轭并除N得x(n) ( )X k 求共轭( )XkFFT 求( )FFT Xk( )FFT XkN 除以( )x n 求共轭1010)(1)()()(NknkNNnnkNWkXNnxWnxkX)(1)(1)()(1)(1010kXDFTNWkXNnxWkXNnxNknkNNknkNn3.4 DFT 的应用举例DFT的广泛应用是在其快速算法FFT出现之后，主要应用在卷积和相关的具体计算，及用DFT近似计算信号的傅里叶变换。本节介绍DFT的两个重要应用：一是用DFT计算卷积；二是用DFT对连续信号或序列进行谱分析。 3.4.1 用DFT计算线性卷积补零至 L求FFT)()()(nhnxny补零至 L求FFT)()()(kHkXkY