离散傅里叶变换

第,3,章,离散傅里叶变换,(DFT),第,3,章 离散傅里叶变换,(DFT),3.1,离散傅里叶变换的定义及物理意义,3.2,离散傅里叶变换的基本性质,3.3,频率域采样,3.4 DFT,的应用举例,习题与上机题,傅里叶变换和,Z,变换是数字信号处理中常用的重要数学变换。对于有限长序列，还有一种更为重要的数学变换，即本章要讨论的离散傅里叶变换,(Discrete Fourier Transform,，,DFT),。,DFT,之所以更为重要，是因为其实质是有限长序列傅里叶变换的有限点离散采样，从而实现了频域离散化，使数字信号处理可以在频域采用数值运算的方法进行，这样就大大增加了数字信号处理的灵活性。更重要的是，,DFT,有多种快速算法，统称为快速傅里叶变换,(Fast Fourier Transform,，,FFT),，从而使信号的实时处理和设备的简化得以实现。,因此，时域离散系统的研究与应用在许多方面取代了传统的连续时间系统。所以说，,DFT,不仅在理论上有重要意义，而且在各种信号的处理中亦起着核心作用。,本章主要讨论,DFT,的定义、物理意义、基本性质以及频域采样和,DFT,的应用举例等内容。,3.1,离散傅里叶变换的定义及物理意义,3.1.1 DFT,的定义,设,x,(,n,),是一个长度为,M,的有限长序列，则定义,x,(,n,),的,N,点离散傅里叶变换为,(3.1.1),X,(,k,),的离散傅里叶逆变换,(Inverse Discrete Fourier Transform, IDFT),为,式中，,，,N,称为,DFT,变换区间长度，,N,M,。通常称,(3.1.1),式和,(3.1.2),式为离散傅里叶变换对。为了叙述简洁，常常用,DFT,x,(,n,),N,和,IDFT,X,(,k,),N,分别表示,N,点离散傅里叶变换和,N,点离散傅里叶逆变换。下面证明,IDFT,X,(,k,),的唯一性。,(3.1.2),把,(3.1.1),式代入,(3.1.2),式，有,由于,所以，在变换区间上满足下式：,IDFT,X,(,k,),N,=,x,(,n,) 0,n,N-,1,由此可见，,(3.1.2),式定义的离散傅里叶逆变换是唯一的。,【,例,3.1.1,】,x,(,n,)=,R,4,(,n,),求,x,(,n,),的,4,点和,8,点,DFT,。,解,设变换区间,N,=4,，则,设变换区间,N,=8,，则,由此例可见，,x,(,n,),的离散傅里叶变换结果与变换区间长度,N,的取值有关。对,DFT,与,Z,变换和傅里叶变换的关系及,DFT,的物理意义进行讨论后，上述问题就会得到解释。,3.1.2 DFT,与傅里叶变换和,Z,变换的关系,设序列,x,(,n,),的长度为,M,，其,Z,变换和,N,(,N,M,),点,DFT,分别为,比较上面二式可得关系式,(3.1.3),或,(3.1.4),(3.1.3),式表明序列,x,(,n,),的,N,点,DFT,是,x,(,n,),的,Z,变换在单位圆上的,N,点等间隔采样。,(3.1.4),式则说明,X,(,k,),为,x,(,n,),的傅里叶变换,X,(e,j,),在区间,0, 2,上的,N,点等间隔采样。这就是,DFT,的物理意义。由此可见，,DFT,的变换区间长度,N,不同，表示对,X,(e,j,),在区间,0, 2,上的采样间隔和采样点数不同，所以,DFT,的变换结果不同。上例中，,x,(,n,)=,R,4,(,n,),，,DFT,变换区间长度,N,分别取,8,、,16,时，,X,(e,j,),和,X,(,k,),的幅频特性曲线图如图,3.1.1,所示。由此容易得到,x,(,n,)=,R,4,(,n,),的,4,点,DFT,为,X,(,k,)=DFT,x,(,n,),4,=4,(,k,),，这一特殊的结果在下面将得到进一步解释。,图,3.1.1,R,4,(,n,),的,FT,和,DFT,的幅度特性关系,3.1.3 DFT,的隐含周期性,前面定义的,DFT,变换对中，,x,(,n,),与,X,(,k,),均为有限长序列，但由于,的周期性，使,(3.1.1),和,(3.1.2),式中的,X,(,k,),隐含周期性，且周期均为,N,。对任意整数,m,，总有,所以,(3.1.1),式中，,X,(,k,),满足：,实际上，任何周期为,N,的周期序列都可以看做长度为,N,的有限长序列,x,(,n,),的周期延拓序列，而,x,(,n,),则是,的一个周期，即,上述关系如图,3.1.2,（,a,）和（,b,）所示。一般称周期序列中从,n,=0,到,N,1,的第一个周期为的主值区间，而主值区间上的序列称为的主值序列。因此,x,(,n,),与的上述关系可叙述为：是,x,(,n,),的周期延拓序列，,x,(,n,),是的主值序列。,(3.1.5),(3.1.6),为了以后叙述简洁，当,N,大于等于序列,x,(,n,),的长度时，将,(3.1.5),式用如下形式表示：,(3.1.7),式中,x,(,n,),N,表示,x,(,n,),以,N,为周期的周期延拓序列，,(,n,),N,表示模,N,对,n,求余，即如果,n,=,MN,+,n,1,0,n,1,N,1,M,为整数,则,(,n,),N,=,n,1,例如，,则有,所得结果符合图,3.1.2(a),和,(b),所示的周期延拓规律。,图,3.1.2,x,(,n,),及其周期延拓序列,应当说明，若,x,(,n,),实际长度为,M,，延拓周期为,N,，则当,N,M,时，,(3.1.5),式仍表示以,N,为周期的周期序列，但,(3.1.6),和,(3.1.7),式仅对,N,M,时成立。图,3.1.2(a),中,x,(,n,),实际长度,M,=6,，当延拓周期,N,=4,时，如图,3.1.2(c),所示。,如果,x,(,n,),的长度为,M,，且，,N,M,，则可写出的离散傅里叶级数表示式,(3.1.8),(3.1.9),式中,即,X,(,k,),为的主值序列。将,(3.1.8),和,(3.1.9),式与,DFT,的定义,(3.1.1),和,(3.1.2),式相比较可知，有限长序列,x,(,n,),的,N,点离散傅里叶变换,X,(,k,),正好是,x,(,n,),的周期延拓序列,x,(,n,),N,的离散傅里叶级数系数的主值序列，即。后面要讨论的频域采样理论将会加深对这一关系的理解。我们知道，周期延拓序列频谱完全由其离散傅里叶级数系数确定，因此，,X,(,k,),实质上是,x,(,n,),的周期延拓序列,x,(,n,),N,的频谱特性，这就是,N,点,DFT,的第二种物理解释（物理意义）。,(3.1.10),现在解释,DFT,R,4,(,n,),4,=4,(,k,),。根据,DFT,第二种物理解释可知，,DFT,R,4,(,n,),4,表示,R,4,(,n,),以,4,为周期的周期延拓序列,R,4,(,n,),4,的频谱特性，因为,R,4,(,n,),4,是一个直流序列，只有直流成分（即零频率成分）。,3.1.4,用,MATLAB,计算序列的,DFT,MATLAB,提供了用快速傅里叶变换算法,FFT(,算法见第,4,章介绍,),计算,DFT,的函数,fft,，其调用格式如下,:,Xk = fft (xn, N),；,调用参数,xn,为被变换的时域序列向量，,N,是,DFT,变换区间长度，当,N,大于,xn,的长度时，,fft,函数自动在,xn,后面补零。函数返回,xn,的,N,点,DFT,变换结果向量,Xk,。当,N,小于,xn,的长度时，,fft,函数计算,xn,的前面,N,个元素构成的,N,长序列的,N,点,DFT,，忽略,xn,后面的元素。 ,Ifft,函数计算,IDFT,，其调用格式与,fft,函数相同，可参考,help,文件。,【,例,3.1.2,】,设,x,(,n,)=,R,4,(,n,),，,X,(e,j,)=FT,x,(,n,),。分别计算,X,(e,j,),在频率区间,0,，,2,上的,16,点和,32,点等间隔采样，并绘制,X,(e,j,),采样的幅频特性图和相频特性图。,解,由,DFT,与傅里叶变换的关系知道，,X,(e,j,),在频率区间,0,，,2,上的,16,点和,32,点等间隔采样，分别是,x,(,n,),的,16,点和,32,点,DFT,。调用,fft,函数求解本例的程序,ep312.m,如下,:,%,例,3.1.2,程序,ep312.m,% DFT,的,MATLB,计算,xn=1 1 1 1; %,输入时域序列向量,xn=R4(n),Xk16=fft(xn, 16); %,计算,xn,的,16,点,DFT,Xk32=fft(xn, 32); %,计算,xn,的,32,点,DFT,%,以下为绘图部分（省略，程序集中有）,程序运行结果如图,3.1.3,所示。,图,3.1.3,程序,ep312.m,运行结果,3.2,离散傅里叶变换的基本性质,3.2.1,线性性质,如果,x,1,(,n,),和,x,2,(,n,),是两个有限长序列，长度分别为,N,1,和,N,2,，且,y,(,n,)=,ax,1,(,n,)+,bx,2,(,n,),式中，,a,、,b,为常数，取,N,=max,N,1,N,2,则,y,(,n,),的,N,点,DFT,为,Y,(,k,)=DFT,y,(,n,),N,=,aX,1,(,k,)+,bX,2,(,k,) 0,k,N,1,(3.2.1),其中,X,1,(,k,),和,X,2,(,k,),分别为,x,1,(,n,),和,x,2,(,n,),的,N,点,DFT,。 ,3.2.2,循环移位性质,1,序列的循环移位,设,x,(,n,),为有限长序列，长度为,M,，,M,N,，则,x,(,n,),的循环移位定义为,y,(,n,)=,x,(,n,+,m,),N,R,N,(,n,) (3.2.2),(3.2.2),式表明，将,x,(,n,),以,N,为周期进行周期延拓得到,，再将左移,m,得到，最后取的主值序列则得到有限长序列,x,(,n,),的循环移位序列,y,(,n,),。,M,=6,N,=8,m,=2,时，,x,(,n,),及其循环移位过程如图,3.2.1,所示。,显然，,y,(,n,),是长度为,N,的有限长序列。观察图,3.2.1,可见，循环移位的实质是将,x,(,n,),左移,m,位，而移出主值区,(0,n,N-,1),的序列值又依次从右侧进入主值区。“循环移位”就是由此得名的。,由循环移位的定义可知，对同一序列,x,(,n,),和相同的位移,m,，当延拓周期,N,不同时，,y,(,n,)=,x,(,n,+,m,),N,R,n,(,n,),则不同。请读者画出,N,=,M,=6,，,m,=2,时，,x,(,n,),的循环移位序列,y,(,n,),波形图。,图,3.2.1,x,(,n,),及其循环移位过程,2, 时域循环移位定理,设,x,(,n,),是长度为,M,(,M,N,),的有限长序列，,y,(,n,),为,x,(,n,),的循环移位，即,则,(3.2.3),其中,证明,令,n,+,m,=,n,，则有,由于上式中求和项,以,N,为周期，因此对其在任一周期上的求和结果相同。将上式的求和区间改在主值区，则得,3, 频域循环移位定理,如果,X,(,k,)=DFT,x,(,n,),N,0,k,N-,1,Y,(,k,)=,X,(,k,+l),N,R,N,(,k,),则,(3.2.4),(3.2.4),式的证明方法与时域循环移位定理类似，直接对,Y,(,k,)=,X,(,k,+l),N,R,N,(,k,),进行,IDFT,即得证。,3.2.3,循环卷积定理,时域循环卷积定理是,DFT,中最重要的定理，具有很强的实用性。已知系统输入和系统的单位脉冲响应，计算系统的输出，以及,FIR,滤波器用,FFT,实现等，都是基于该定理的。下面首先介绍循环卷积的概念和计算循环卷积的方法，然后介绍循环卷积定理。,1, 两个有限长序列的循环卷积,设序列,h,(,n,),和,x,(,n,),的长度分别为,N,和,M,。,h,(,n,),与,x,(,n,),的,L,点循环卷积定义为,（,3.2.5,）,式中，,L,称为循环卷积区间长度，,L,max,N,，,M,。上式显然与第,1,章介绍的线性卷积不同，为了区别线性卷积，用 表示循环卷积，用表示,L,点循环卷积，即,y,c,(,n,)=,h,(,n,),x,(,n,),。观察,(3.2.5),式，,x,(,n,m,),L,是以,L,为周期的周期信号，,n,和,m,的变化区间均是,0,L-,1,，因此直接计算该式比较麻烦。计算机中采用矩阵相乘或快速傅里叶变换,(FFT),的方法计算循环卷积。下面介绍用矩阵计算循环卷积的公式。,当,n,= 0, 1, 2,L,1,时，由,x,(,n,),形成的序列为,: ,x,(0),x,(1),x,(,L,1),。令,n,=0,m,=0, 1,L,1,，由式（,3.2.5,）中,x,(,n,-,m,),L,形成,x,(,n,),的循环倒相序列为,与序列,x,(,n,),进行对比，相当于将第一个序列值,x,(0),不动，将后面的序列反转,180,再放在,x,(0),的后面。这样形成的序列称为,x,(,n,),的循环倒相序列。,令,n,= 1,m,= 0, 1,L,-1,，由式（,3.2.5,）中,x,(,n,-,m,),L,形成的序列为,观察上式等号右端序列，它相当于,x,(,n,),的循环倒相序列向右循环移一位，即向右移,1,位，移出区间,0,L,1,的序列值再从左边移进。,再令,n,= 2,m,= 0, 1,L-,1,，此时得到的序列又是上面的序列向右循环移,1,位。依次类推，当,n,和,m,均从,0,变化到,L,-1,时，得到一个关于,x,(,n,m,),L,的矩阵如下：,（,3.2.6,）,上面矩阵称为,x,(,n,),的,L,点“循环卷积矩阵”，其特点是,:,（,1,） 第,1,行是序列,x,(0),x,(1),x,(,L,1),的循环倒相序列。,注意，如果,x,(,n,),的长度,M,L,，则需要在,x,(,n,),末尾补,L,M,个零后，再形成第一行的循环倒相序列。,（,2,） 第,1,行以后的各行均是前一行向右循环移,1,位形成的。,（,3,） 矩阵的各主对角线上的序列值均相等。,有了上面介绍的循环卷积矩阵，就可以写出式,(3.2.5),的矩阵形式如下,:,（,3.2.7,）,按照上式，可以在计算机上用矩阵相乘的方法计算两个序列的循环卷积，这里关键是先形成循环卷积矩阵。上式中,如果,h,(,n,),的长度,N,L,，则需要在,h,(,n,),末尾补,L-N,个零。,【,例,3.2.1,】,计算下面给出的两个长度为,4,的序列,h,(,n,),与,x,(,n,),的,4,点和,8,点循环卷积。,解,按照式（,3.2.21,）写出,h,(,n,),与,x,(,n,),的,4,点循环卷积矩阵形式为,h,(,n,),与,x,(,n,),的,8,点循环卷积矩阵形式为,h,(,n,),和,x,(,n,),及其,4,点和,8,点循环卷积结果分别如图,3.2.2(a),、,(b),、,(c),和,(d),所示。请读者计算验证本例的,8,点循环卷积结果等于,h,(,n,),与,x,(,n,),的线性卷积结果。后面将证明，当循环卷积区间长度,L,大于等于,y,(,n,) =,h,(,n,)*,x,(,n,),的长度时，循环卷积结果就等于线性卷积。,图,3.2.2,序列及其循环卷积波形,2.,环卷积定理,有限长序列,x,1,(,n,),和,x,2,(,n,),的长度分别为,N,1,和,N,2,，,N,= max,N,1,N,2,，,x,1,(,n,),和,x,2,(,n,),的,N,点循环卷积为,则,x,(,n,),的,N,点,DFT,为,其中,N,（,3.2.9,）,（,3.2.8,）,证明,直接对,(3.2.8),式两边进行,DFT,，则有,令,n-m,=,n,，则有,因为上式中,是以,N,为周期的，所以对其在任一个周期上求和的结果不变。因此,由于，因此,即循环卷积亦满足交换律。,作为习题请读者证明以下频域循环卷积定理,:,如果,x,(,n,)=,x,1,(,n,),x,2,(,n,),，则,(3.2.10a),N,或,(3.2.10b),式中,相对频域循环卷积定理，称,(3.2.9),式为时域循环卷积定理。,N,3.2.4,复共轭序列的,DFT,设,x,*(,n,),是,x,(,n,),的复共轭序列，长度为,N,，,X,(,k,)=DFT,x,(,n,),N,，则,且,X(N,)=X(0),。,证明,根据,DFT,的唯一性，只要证明,(3.2.11),式右边等于左边即可。,(3.2.11),又由,X,(,k,),的隐含周期性，有,X,(,N,)=,X,(0),用同样的方法可以证明,(3.2.12),3.2.5 DFT,的共轭对称性,1, 有限长共轭对称序列和共轭反对称序列,为了区别于傅里叶变换中所定义的共轭对称（或共轭反对称）序列，下面用,x,ep,(,n,),和,x,op,(,n,),分别表示有限长共轭对称序列和共轭反对称序列，则二者满足如下关系式：,(3.2.13a),(3.2.13b),当,N,为偶数时，将上式中的,n,换成,N,/2,n,，可得到：,上式更清楚地说明了有限长序列共轭对称序列是关于,n,=,N,/2,点对称。容易证明，如同任何实函数都可以分解成偶对称分量和奇对称分量一样，任何有限长序列,x,(,n,),都可以表示成其共轭对称分量和共轭反对称分量之和，即,(3.2.14),将上式中的,n,换成,N-n,，并取复共轭，再将,(3.2.13a),式和,(3.2.13b),式代入，得到：,(3.2.15),(3.2.14),式分别加减,(3.2.15),式，可得,(3.2.16a),(3.2.16b),2,DFT,的共轭对称性,(1),如果将,x,(,n,),表示为,x,(,n,)=,x,r,(,n,)+j,x,i,(,n,) (3.2.17),其中,那么，由,(3.2.11),式和,(3.2.16a),式可得,(3.2.18),由,(3.2.11),式和,(3.2.16b),式可得,(3.2.19),由,DFT,的线性性质即可得,(3.2.20),其中，,X,op,(,k,)=DFT,x,r,(,n,),是,X,(,k,),的共轭对称分量，,X,op,(,k,)=DFTj,x,i,(,n,),是,X,(,k,),的共轭反对称分量。,(2),如果将,x,(,n,),表示为,(3.2.21),其中，,是,x,(,n,),的共轭对称分量，,是,x,(,n,),的共轭反对称分量,那么，由,(3.2.12),式可得,因此,(3.2.22),其中,综上所述，可总结出,DFT,的共轭对称性质：如果序列,x,(,n,),的,DFT,为,X,(,k,),，则,x,(,n,),的实部和虚部（包括,j,）的,DFT,分别为,X,(,k,),的共轭对称分量和共轭反对称分量；而,x,(,n,),的共轭对称分量和共轭反对称分量的,DFT,分别为,X,(,k,),的实部和虚部乘以,j,。,另外，请读者根据上述共轭对称性证明有限长实序列,DFT,的共轭对称性（见本章习题题,7,）。,设,x,(,n,),是长度为,N,的实序列，且,X,(,k,)=DFT,x,(,n,),N,，则,X,(,k,),满足如下对称性：,(1),X,(,k,),共轭对称，即,X,(,k,)=,X,*,(,N,-,k,),k,=0, 1,N,-1 (3.2.23),(2),如果,x,(,n,),是偶对称序列，即,x,(,n,)=,x,(,N-n,),，则,X,(,k,),实偶对称，即,X,(,k,)=,X,(,N-k,) (3.2.24),(3),如果是奇对称序列，即,x,(,n,)=-,x,(,N-n,),，则,X,(,k,),纯虚奇对称，即,X,(,k,)=-,X,(,N-k,) (3.2.25),实际中经常需要对实序列进行,DFT,，利用上述对称性质，可减少,DFT,的运算量，提高运算效率。例如，计算实序列的,N,点,DFT,时，当,N,=,偶数时，只需计算,X,(,k,),的前面,N,/2+1,点，而,N,=,奇数时，只需计算,X,(,k,),的前面,(,N,+1)/2,点，其他点按照,(3.2.23),式即可求得。例如，,X,(,N-,1)=,X,*,(1),X,(,N-,2)=,X,*,(2),这样可以减少近一半运算量。,【,例,3.2.2,】,利用,DFT,的共轭对称性，设计一种高效算法，通过计算一个,N,点,DFT,，就可以计算出两个实序列,x,1,(,n,),和,x,2,(,n,),的,N,点,DFT,。,解,构造新序列,x,(,n,)=,x,1,(,n,)+j,x,2,(,n,),，对,x,(,n,),进行,DFT,，得到：,由,(3.2.17),、,(3.2.18),和,(3.2.19),式得到：,所以，由,X,(,k,),可以求得两个实序列,x,1,(,n,),和,x,2,(,n,),的,N,点,DFT,：,3.3,频 率 域 采 样,时域采样定理告诉我们，在一定条件下，可以由时域离散采样信号恢复原来的连续信号。那么能不能也由频域离散采样恢复原来的信号（或原连续频率函数）？其条件是什么？内插公式又是什么形式？本节就上述问题进行讨论。,设任意序列,x,(,n,),的,Z,变换为,且,X,(,z,),的收敛域包含单位圆,(,即,x,(,n,),存在傅里叶变换,),。在单位圆上对,X,(,z,),等间隔采样,N,点,得到：,显然，,(3.3.1),式表示在区间,0, 2,上对,x,(,n,),的傅里叶变换,X,(e,j,),的,N,点等间隔采样。将,X,(,k,),看做长度为,N,的有限长序列,x,N,(,n,),的,DFT,，即,下面推导序列,x,N,(,n,),与原序列,x,(,n,),之间的关系，并导出频域采样定理。,(3.3.1),由,DFT,与,DFS,的关系可知，,X,(,k,),是,x,N,(,n,),以,N,为周期的周期延拓序列的离散傅里叶级数系数的主值序列，即,将,(3.3.1),式代入上式得,式中,-,=,-,=,-,-,=,-,=,-,=,=,m,N,k,n,m,k,N,N,k,m,kn,N,km,N,W,N,m,x,W,W,m,x,N,n,x,1,0,),(,1,0,1,),(,),(,1,),(,因此,所以,(3.3.3),(3.3.2),式,(3.3.3),说明，,X,(,z,),在单位圆上的,N,点等间隔采样,X,(,k,),的,N,点,IDFT,是原序列,x,(,n,),以,N,为周期的周期延拓序列的主值序列。综上所述，可以总结出频域采样定理：,如果序列,x,(,n,),的长度为,M,，则只有当频域采样点数,N,M,时，才有,x,N,(,n,)=IDFT,X,(,k,)=,x,(,n,),即可由频域采样,X,(,k,),恢复原序列,x,(,n,),，否则产生时域混叠现象。,满足频域采样定理时，频域采样序列,X,(,k,),的,N,点,IDFT,是原序列,x,(,n,),，所以必然可以由,X,(,k,),恢复,X,(,z,),和,X,(e,j,),。下面推导用频域采样,X,(,k,),表示,X,(,z,),和,X,(e,j,),的内插公式和内插函数。设序列,x,(,n,),长度为,M,，在频域,0, 2,上等间隔采样,N,点，,N,M,，则有,因为满足频域采样定理，所以式中,将上式代入,X,(,z,),的表示式中，得到,:,(3.3.4a),式中，因此,(3.3.4b),令,(3.3.5),则,(3.3.6),式,(3.3.6),称为用,X,(,k,),表示,X,(,z,),的内插公式，,k,(,z,),称为内插函数。将,z,=e,j,代入,(3.3.4a),式，并进行化简，可得,(3.3.7),(3.3.8),在数字滤波器的结构与设计中，我们将会看到，频域采样理论及有关公式可提供一种有用的滤波器结构和滤波器设计途径，,(3.3.7),式有助于分析,FIR,滤波器频率采样设计法的逼近性能。,【,例,3.3.1,】,长度为,26,的三角形序列,x,(,n,),如图,3.3.1(a),所示。编写,MATLAB,程序验证频域采样理论。,解,解题思想： 先计算,x,(,n,),的,32,点,DFT,，得到其频谱函数,X,(e,j,),在频率区间,0,2,上等间隔,32,点采样,X,32,(,k,),，再对,X,32,(,k,),隔点抽取，得到,X,(e,j,),在频率区间,0,2,上等间隔,16,点采样,X,16,(,k,),。最后分别对,X,16,(,k,),和,X,32,(,k,),求,IDFT,得到：,绘制,x,16,(,n,),和,x,32,(,n,),波形图验证频域采样理论。,MATLAB,求解程序,ep331.m,如下：,%,数字信号处理,(,第三版）,第,3,章例,3.3.1,程序,ep331.,%,频域采样理论验证,M=26; N=32; n=0:M; ,xa=0:M/2; xb= ceil(M/2)-1:-1:0; xn=,xa, xb,;,%,产生,M,长三角波序列,x(n),Xk=fft(xn, 512);,%512,点,FFT,x(n),X32k=fft(xn, 32);,%32,点,FFT,x(n),x32n=ifft(X32k);,%32,点,IFFT,X32(k),得到,x32(n),X16k=X32k(1:2:N); %,隔点抽取,X32k,得到,X16(k),x16n=ifft(X16k, N/2); %16,点,IFFT,X16(k),得到,x16(n),以下绘图部分省略。,程序运行结果如图,3.3.1,所示。图,3.3.1(a),和,(b),分别为,X,(e,j,),和,x,(,n,),的波形；图,3.3.1(c),和,(d),分别为,X,(e,j,),的,16,点采样,|,X,16,(,k,)|,和,x,16,(,n,)=IDFT,X,16,(,k,),16,波形图；图,3.3.1(e),和,(f),分别为,X,(e,j,),的,32,点采样,|,X,32,(,k,)|,和,x,32,(,n,)=IDFT,X,32,(,k,),32,波形图；由于实序列的,DFT,满足共轭对称性，因此频域图仅画出,0,，,上的幅频特性波形。本例中,x,(,n,),的长度,M,=26,。从图中可以看出，当采样点数,N,=16,M,时，无时域混叠失真，,x,32,(,n,)=IDFT,X,32,(,k,),=,x,(,n,),。,图,3.3.1,频域采样定理验证,3.4 DFT,的应用举例,3.4.1,用,DFT,计算线性卷积,用,DFT,计算循环卷积很简单。设,h,(,n,),和,x,(,n,),的长度分别为,N,和,M,其,L,点循环卷积为,且,L,则由,DFT,的时域循环卷积定理有,由此可见，循环卷积既可以在时域直接计算，也可以按照图,3.4.1,所示的计算框图在频域计算。由于,DFT,有快速算法，当,L,很大时，在频域计算循环卷积的速度快得多，因而常用,DFT(FFT),计算循环卷积。,图,3.4.1,用,DFT,计算循环卷积的原理框图,在实际应用中，为了分析时域离散线性时不变系统或者对序列进行滤波处理等，需要计算两个序列的线性卷积。与计算循环卷积一样，为了提高运算速度，也希望用,DFT(FFT),计算线性卷积。而,DFT,只能直接用来计算循环卷积，因此，下面先导出线性卷积和循环卷积之间的关系以及循环卷积与线性卷积相等的条件，最后得出用图,3.4.1,线性卷积的条件。,假设,h,(,n,),和,x,(,n,),都是有限长序列，长度分别是,N,和,M,。它们的线性卷积和循环卷积分别表示如下：,其中,所以,(3.4.1),(3.4.2),L,对照,(3.4.1),式可以看出，上式中,即,(3.4.3),(3.4.3),式说明，,y,c,(,n,),等于,y,l,(,n,),以,L,为周期的周期延拓序列的主值序列。我们知道，,y,l,(,n,),长度为,N,M,1,，因此只有当循环卷积长度,L,N,M,1,时，,y,l,(,n,),以,L,为周期进行周期延拓时才无时域混叠现象。此时取其主值序列显然满足,y,c,(,n,)=,y,l,(,n,),。由此证明了循环卷积等于线性卷积的条件是,L,N,M,1,。图,3.4.2,中画出了,h,(,n,),、,x,(,n,),、,h,(,n,)*,x,(,n,),以及,L,分别取,6,、,8,、,10,时,h,(,n,),L,x,(,n,),的波形。由于,h,(,n,),长度,N,=4,，,x,(,n,),长度,M,=5,，,N,M,1=8,，因此只有,L,8,时，,h,(,n,),L,x,(,n,),波形才与,h,(,n,)*,x,(,n,),相同。,图,3.4.2,线性卷积与循环卷积波形图,综上所述，取,L,N,M,1,，则可按照如图,3.4.1,所示的计算框图用,DFT(FFT),计算线性卷积。其中,DFT,和,IDFT,通常用快速算法,(FFT),来实现，故常称其为快速卷积。,实际上，经常遇到两个序列的长度相差很大的情况，例如,M,N,。若仍选取,L,N,M,1,，以,L,为循环卷积区间，并用上述快速卷积法计算线性卷积，则要求对短序列补很多零点，而且长序列必须全部输入后才能进行快速计算。因此要求存储容量大，运算时间长，并使处理延时很大，不能实现实时处理。,况且在某些应用场合，序列长度不定或者认为是无限长，如电话系统中的语音信号和地震检测信号等。显然，在要求实时处理时，直接套用上述方法是不行的。解决这个问题的方法是将长序列分段计算，这种分段处理方法有重叠相加法和重叠保留法两种。下面只介绍重叠相加法，重叠保留法作为本章习题题,21,，留给读者讨论。,设序列,h,(,n,),长度为,N,，,x,(,n,),为无限长序列。将,x,(,n,),等长分段，每段长度取,M,，则,（,3.4.4a,）,于是，,h,(,n,),与,x,(,n,),的线性卷积可表示为,（,3.4.4b,）,式中,(3.4.4b),式说明，计算,h,(,n,),与,x,(,n,),的线性卷积时，可先计算分段线性卷积,y,k,(,n,)=,h,(,n,),*,x,k,(,n,),，然后把分段卷积结果叠加起来即可，如图,3.4.3,所示。每一分段卷积,y,k,(,n,),的长度为,N,M,1,，因此相邻分段卷积,y,k,(,n,),与,y,k,1,(,n,),有,N,1,个点重叠，必须把重叠部分的,y,k,(,n,),与,y,k,1,(,n,),相加，才能得到正确的卷积序列,y,(,n,),。,显然,可用图,3.4.1,所示的快速卷积法计算分段卷积,y,k,(,n,),， 其中,L,=,N,M,1,。由图,3.4.3,可以看出，当第二个分段卷积,y,1,(,n,),计算完后，叠加重叠点便可得输出序列,y,(,n,),的前,2,M,个值；同样道理，分段卷积,y,i,(,n,),计算完后，就可得到,y,(,n,),第,i,段的,M,个序列值。因此，这种方法不要求大的存储容量，且运算量和延时也大大减少，最大延时,T,Dmax,=2,MT,s,+,T,o,，,T,s,是系统采样间隔，,T,o,是计算,1,个分段卷积所需时间，一般要求,T,o,MT,s,。这样，就实现了边输入边计算边输出，如果计算机的运算速度快，可以实现实时处理。,图,3.4.3,用重叠相加法计算,线性卷积时域关系示意图,用,DFT,计算分段卷积,y,k,(,n,),的方法如下：,(1),i,=0,；,L,=,N,M,1,；计算并保存,H,(,k,)=DFT,h,(,n,),L,;,(2),读入,x,k,(,n,)=,x,(,n,),R,M,(,n,kM,),，构造变换区间,0,，,L,1,上的序列，实际中就是将,x,i,(,n,),的,M,个值存放在长度为,M,的数组中,并计算,(3),；,(4),，,n,= 0,1,2,L,1,；,(5),计算：,(6),i,=,i,1,，返回,(2),。,应当说明，一般,x,(,n,),是因果序列，假设初始条件,y,1,(,n,)=0,。,MATLAB,信号处理工具箱中提供了一个函数,fftfilt,，该函数用重叠相加法实现线性卷积的计算。调用格式为：,y=fftfilt(h, x,，,M),。式中, h,是系统单位脉冲响应向量；,x,是输入序列向量；,y,是系统的输出序列向量（,h,与,x,的卷积结果）；,M,是由用户选择的输入序列,x,的分段长度，缺省,M,时，默认输入序列,x,的分段长度,M=512,。,【,例,3.4.1,】,假设,h,(,n,)=,R,5,(,n,),，,x,(,n,)=,cos(,n,/10),+cos(2,n,/5),u,(,n,),，用重叠相加法计算,y,(,n,)=,h,(,n,)*,x,(,n,),，并画出,h,(,n,),、,x,(,n,),和,y,(,n,),的波形。,解,h,(,n,),的长度为,N,=5,对,x,(,n,),进行分段，每段长度为,M,=10,。计算,h,(,n,),和,x,(,n,),的线性卷积的,MATLAB,程序如下：,%,例,3.4.1,重叠相加法的,MATLAB,实现程序：,ep341.m,Lx=41; N=5; M=10; %Lx,为信号序列,x(n),长度,hn=ones(1, N); hn1=,hn zeros(1, Lx-N),;,%,产生,h(n),，其后补零是为了绘图好看,n=0:L-1; ,xn=cos(pi*n/10)+cos(2*pi*n/5); %,产生,x(n),的,Lx,个样值,yn=fftfilt(hn, xn, M); %,调用,fftfilt,用重叠相加法计算卷积,%=,%,以下为绘图部分,省略,运行程序画出,h,(,n,),、,x,(,n,),和,y,(,n,),的波形如图,3.4.4,所示。请读者从理论上证明,y,(,n,),的稳态波形是单一频率的正弦波。,运行绘图程序,fig345.m,可以得到用重叠相加法求解本例题的,x,k,(,n,), y,k,(,n,),和,y(,n,)=y,0,(,n,)+y,1,(,n,)+y,2,(,n,)+y,3,(,n,),如图,3.4.5,所示。,图,3.4.4,例,3.4.1,的求解程序运行结果,图,3.4.5,重叠相加法时域波形,3.4.2,用,DFT,对信号进行谱分析,1, 用,DFT,对连续信号进行谱分析,工程实际中，经常遇到连续信号,x,a,(,t,),，其频谱函数,X,a,(j,),也是连续函数。为了利用,DFT,对,x,a,(,t,),进行频谱分析，先对,x,a,(,t,),进行时域采样，得到,x,(,n,)=,x,a,(,nT,),，再对,x,(,n,),进行,DFT,，得到的,X,(,k,),则是,x,(,n,),的傅里叶变换,X,(e,j,),在频率区间,0,，,2,上的,N,点等间隔采样。这里,x,(,n,),和,X,(,k,),均为有限长序列。,然而，,由傅里叶变换理论知道，若信号持续时间有限长，则其频谱无限宽；若信号的频谱有限宽，则其持续时间必然为无限长,。所以严格地讲，持续时间有限的带限信号是不存在的。因此，按采样定理采样时，上述两种情况下的采样序列,x,(,n,)=,x,a,(,nT,),均应为无限长，不满足,DFT,的变换条件。实际上对频谱很宽的信号，为防止时域采样后产生频谱混叠失真，可用预滤波器滤除幅度较小的高频成分，使连续信号的带宽小于折叠频率。,对于持续时间很长的信号，采样点数太多,以致无法存储和计算，只好截取有限点进行,DFT,。由上述可见，用,DFT,对连续信号进行频谱分析必然是近似的，其近似程度与信号带宽、采样频率和截取长度有关。实际上从工程角度看，滤除幅度很小的高频成分和截去幅度很小的部分时间信号是允许的。因此，,在下面分析中，假设,x,a,(,t,),是经过预滤波和截取处理的有限长带限信号。,设连续信号,x,a,(,t,),持续时间为,T,p,，最高频率为,f,c,，如图,3.4.6(a),所示。,x,a,(,t,),的傅里叶变换为,对,x,a,(,t,),以采样间隔,T,(2,f,c,),-1,(,即,f,s,=1/,T,2,f,c,),采样得,x,(,n,)=,x,a,(,nT,),。设共采样,N,点，并对,X,a,(,jf,),作零阶近似,(,t,=,nT,dt,=,T,),得,显然，,X,a,(,jf,),仍是,f,的连续周期函数，,x,(,n,),和,X,(,jf,),如图,3.4.5(b),所示。 对,X,(,jf,),在区间,0,F,s,上等间隔采样,N,点， 采样间隔为,F,， 如图,3.4.5(c),所示。 参数,F,s,、,T,p,、,N,和,F,满足如下关系式：,由于,NT=T,p,， 所以,(3.4.5),(3.4.6),将,f=kF,和式,(3.4.5),代入,X,(,jf,),中可得,X,a,(,jf,),的采样,0,k,N,-1,令,则,(3.4.8),(3.4.7),同理,由,可以推出,图,3.4.6,用,DFT,分析连续信号谱的原理示意图,(3.4.7),式说明，可以通过对连续信号采样并进行,DFT,再乘以,T,，近似得到模拟信号频谱的周期延拓函数在第一个周期,0,f,s,上的,N,点等间隔采样 ，如图,3.4.6,所示。对满足假设的持续时间有限的带限信号，在满足时域采样定理时， 包含了模拟信号频谱的全部信息（,k,=0, 1, 2,N,/2,表示正频率频谱采样,;,k,=,N,/2+1,，,N,/2+2,，,，,N,1,表示负频率频谱采样）。,所以，上述分析方法不丢失信息，即可由,X,(,k,),恢复,X,a,(j,),或,x,a,(,t,),，但直接由分析结果,X,(,k,),看不到,X,a,(j,),的全部频谱特性，而只能看到,N,个离散采样点的谱线，这就是所谓的栅栏效应。对实信号，其频谱函数具有共轭对称性，所以分析正频率频谱就足够了。不存在频谱混叠失真时，正频率,0,F,s,/2,频谱采样为,（,3.4.12,）,值得注意，如果,x,a,(,t,),持续时间无限长，上述分析中要进行截断处理，所以会产生所谓的截断效应，从而使谱分析产生误差。本节最后将讨论上述误差问题产生的原因及改进措施。 ,下面举例说明截断效应。理想低通滤波器的单位冲激响应,h,a,(,t,),及其频响函数,H,a,(,f,),如图,3.4.7(a),、,(b),所示（图,3.4.7(a),中只画出,h,a,(,t,),所截取的一段）。图中，,现在用,DFT,来分析,h,a,(,t,),的频率响应特性。由于,h,a,(,t,),的持续时间为无穷长，因此要截取一段,T,p,，假设,T,p,=8 s,，采样间隔,T,=0.25 s(,即采样频率,F,s,=4 Hz),，采用点数,N,=,T,p,/,T,=32,；频域采样间隔,F,=1/,T,p,=0.125 Hz,；由于,h,a,(,t,),为实信号，因此仅取正频率,0,f,s,/2,频谱采样：,其中,H,a,(,kF,),如图,3.4.7(c),中黑点所示。由图可见，低频部分近似理想低通频响特性，而高频误差较大，且整个频响都有波动。这些误差就是由于对,h,a,(,t,),截断所产生的，所以通常称之为截断效应。为减少这种截断误差，可适当加长,T,p,，增加采样点数,N,或用窗函数处理后再进行,DFT,。有关窗函数的内容将在,FIR,数字滤波器设计中详细叙述。,图,3.4.7,用,DFT,计算理想低通滤波器的频响曲线,在对连续信号进行谱分析时，主要关心两个问题，这就是谱分析范围和频率分辨率。谱分析范围为,0,F,s,/2,，直接受采样频率,F,s,的限制。为了不产生频率混叠失真，通常要求信号的最高频率,f,c,F,s,/2,。频率分辨率用频率采样间隔,F,描述，,F,表示谱分析中能够分辨的两个频谱分量的最小间隔。显然，,F,越小，谱分析就越接近,X,a,(j,f,),，所以,F,较小时，我们称频率分辨率较高。下面讨论用,DFT,对连续信号谱分析的参数选择原则。,在已知信号的最高频率,f,c,(,即谱分析范围,),时，为了避免频率混叠现象，要求采样速率,F,s,满足下式：,(3.4.13),按照,(3.4.11),式，谱分辨率,F,=,F,s,/,N,，如果保持采样点数,N,不变，要提高频谱分辨率,(,减小,F,),，就必须降低采样频率，采样频率的降低会引起谱分析范围变窄和频谱混叠失真。如维持,F,s,不变，为提高频率分辨率可以增加采样点数,N,，因为,，,只有增加对信号的观察时间,T,p,，才能增加,N,。,T,p,和,N,可以按照下面两式进行选择：,(3.4.14),(3.4.15),【,例,3.4.2,】,对实信号进行谱分析，要求谱分辨率,F,10 Hz,，信号最高频率,f,c,=2.5 kHz,， 试确定最小记录时间,T,p min,，最大的采样间隔,T,max,，最少的采样点数,N,min,。如果,f,c,不变，要求谱分辨率提高,1,倍，最少的采样点数和最小的记录时间是多少？,解,因此,T,p min,=0.1 s,。因为要求,F,s,2,f,c,，所以,为使用,DFT,的快速算法,FFT,，希望,N,符合,2,的整数幂，为此选用,N,=512,点。为使频率分辨率提高,1,倍，即,F,=5 Hz,，要求：,用快速算法,FFT,计算时，选用,N,=1024,点。,上面分析了为提高谱分辨率，又保持谱分析范围不变，必须增长记录时间,T,p,，增加采样点数。应当注意，这种提高谱分辨率的条件是必须满足时域采样定理，甚至采样速率,F,s,取得更高。,2, 用,DFT,对序列进行谱分析,我们知道单位圆上的,Z,变换就是序列的傅里叶变换，即,X,(e,j,),是,的连续周期函数。如果对序列,x,(,n,),进行,N,点,DFT,得到,X,(,k,),，则,X,(,k,),是在区间,0,2,上对,X,(e,j,),的,N,点等间隔采样，频谱分辨率就是采样间隔,2/,N,。因此序列的傅里叶变换可利用,DFT(,即,FFT),来计算。,对周期为,N,的周期序列，由,(2.3.10),式知道，其频谱函数为,其中,由于以,N,为周期，因而,X,(e,j,),也是以,2,为周期的离散谱，每个周期有,N,条谱线，第,k,条谱线位于,=(2/,N,),k,处，代表的,k,次谐波分量。而且，谱线的相对大小与成正比。由此可见，周期序列的频谱结构可用其离散傅里叶级数系数表示。由,DFT,的隐含周期性知道，截取的主值序列,，并进行,N,点,DFT,，得到：,所以可用,X,(,k,),表示的频谱结构。