离散傅里叶变换课件

1、第3章 离散傅里叶变换(DFT),3.1 离散傅里叶变换的定义及物理意义 3.2 离散傅里叶变换的基本性质 3.3 频率域采样 3.4 DFT的应用举例 习题与上机题,傅里叶变换和Z变换是数字信号处理中常用的重要数学变换。对于有限长序列，还有一种更为重要的数学变换，即本章要讨论的离散傅里叶变换(Discrete Fourier Transform，DFT)。DFT之所以更为重要，是因为其实质是有限长序列傅里叶变换的有限点离散采样，从而实现了频域离散化，使数字信号处理可以在频域采用数值运算的方法进行，这样就大大增加了数字信号处理的灵活性。更重要的是，DFT有多种快速算法，统称为快速傅里叶变换(F

2、ast 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)为 式

3、中， ，N称为DFT变换区间长度，NM。通常称(3.1.1)式和(3.1.2)式为离散傅里叶变换对。为了叙述简洁，常常用DFTx(n)N和IDFTX(k)N分别表示N点离散傅里叶变换和N点离散傅里叶逆变换。下面证明IDFTX(k)的唯一性。,(3.1.2),把(3.1.1)式代入(3.1.2)式，有 由于 ,所以，在变换区间上满足下式： IDFTX(k)N=x(n) 0nN-1 由此可见，(3.1.2)式定义的离散傅里叶逆变换是唯一的。 【例3.1.1】 x(n)=R4(n), 求x(n)的4点和8点DFT。 解 设变换区间N=4，则,设变换区间N=8，则 由此例可见，x(n)的离散傅里叶变换

4、结果与变换区间长度N的取值有关。对DFT与Z变换和傅里叶变换的关系及DFT的物理意义进行讨论后，上述问题就会得到解释。,3.1.2 DFT与傅里叶变换和Z变换的关系 设序列x(n)的长度为M，其Z变换和N(NM)点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(ej)在区间0, 2上的N点等间隔采样。这就是DFT的物理意义。由此可见，DFT的变换区间长度N不同，表示对X(ej)在区间0, 2上的采样间隔和采样点数不同，所以

5、DFT的变换结果不同。上例中， x(n)=R4(n)，DFT变换区间长度N分别取8、16时，X(ej)和X(k)的幅频特性曲线图如图3.1.1所示。由此容易得到x(n)=R4(n)的4点DFT为X(k)=DFTx(n)4=4(k)，这一特殊的结果在下面将得到进一步解释。 ,图3.1.1 R4(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的周期序列

6、 都可以看做长度为N的有限长序列x(n)的周期延拓序列，而x(n)则是 的一个周期，即,为了以后叙述简洁，当N大于等于序列x(n)的长度时，将(3.1.5)式用如下形式表示： (3.1.7) 式中x(n) N表示x(n)以N为周期的周期延拓序列，(n)N表示模N对n求余，即如果 n=MN+n1 0n1N1, M为整数 则 (n)N=n1 例如， , 则有 所得结果符合图3.1.2(a)和(b)所示的周期延拓规律。,图3.1.2 x(n)及其周期延拓序列,现在解释DFTR4(n)4=4(k)。根据DFT第二种物理解释可知，DFTR4(n)4表示R4(n)以4为周期的周期延拓序列R4(n)4的频谱

7、特性，因为R4(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

8、.2】 设x(n)=R4(n)，X(ej)=FTx(n)。分别计算X(ej)在频率区间0，2上的16点和32点等间隔采样，并绘制X(ej)采样的幅频特性图和相频特性图。 解 由DFT与傅里叶变换的关系知道，X(ej)在频率区间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 %

9、以下为绘图部分（省略，程序集中有） 程序运行结果如图3.1.3所示。,图3.1.3 程序ep312.m 运行结果,3.2 离散傅里叶变换的基本性质 3.2.1 线性性质 如果x1(n)和x2(n)是两个有限长序列，长度分别为N1和N2，且 y(n)=ax1(n)+bx2(n) 式中，a、b为常数，取N=maxN1, N2, 则y(n)的N点DFT为 Y(k)=DFTy(n)N=aX1(k)+bX2(k) 0kN1 (3.2.1) 其中X1(k)和X2(k)分别为x1(n)和x2(n)的N点DFT。 ,显然，y(n)是长度为N的有限长序列。观察图3.2.1可见，循环移位的实质是将x(n)左移m位

10、，而移出主值区(0nN-1)的序列值又依次从右侧进入主值区。“循环移位”就是由此得名的。 由循环移位的定义可知，对同一序列x(n)和相同的位移m，当延拓周期N不同时，y(n)=x(n+m)NRn(n)则不同。请读者画出N = M=6，m=2时，x(n)的循环移位序列y(n)波形图。,图3.2.1 x(n)及其循环移位过程,2 时域循环移位定理 设x(n)是长度为M(MN)的有限长序列，y(n)为x(n)的循环移位，即 则 (3.2.3) 其中 ,证明 令n+m=n，则有 由于上式中求和项 以N为周期，因此对其在任一周期上的求和结果相同。将上式的求和区间改在主值区，则得 ,3 频域循环移位定理

11、如果 X(k)=DFTx(n)N 0kN-1 Y(k)=X(k+l)NRN(k) 则 (3.2.4) (3.2.4)式的证明方法与时域循环移位定理类似，直接对Y(k)=X(k+l)NRN(k)进行IDFT即得证。,3.2.3 循环卷积定理 时域循环卷积定理是DFT中最重要的定理，具有很强的实用性。已知系统输入和系统的单位脉冲响应，计算系统的输出，以及FIR滤波器用FFT实现等，都是基于该定理的。下面首先介绍循环卷积的概念和计算循环卷积的方法，然后介绍循环卷积定理。 1 两个有限长序列的循环卷积 设序列h(n)和x(n)的长度分别为N和M。h(n)与x(n)的L点循环卷积定义为 （3.2.5）,

12、式中，L称为循环卷积区间长度，LmaxN，M。上式显然与第1章介绍的线性卷积不同，为了区别线性卷积，用 表示循环卷积，用 表示L点循环卷积，即yc(n)=h(n) x(n)。观察(3.2.5)式，x(nm)L是以L为周期的周期信号，n和m的变化区间均是0, L-1，因此直接计算该式比较麻烦。计算机中采用矩阵相乘或快速傅里叶变换(FFT)的方法计算循环卷积。下面介绍用矩阵计算循环卷积的公式。,当n = 0, 1, 2, , L1时，由x(n)形成的序列为: x(0), x(1), , x(L1)。令n=0, m=0, 1, , L1，由式（3.2.5）中x(n-m)L形成x(n)的循环倒相序列为

13、 与序列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, L1的序列值再从左边移进。 再令n = 2, m = 0, 1, , L-1，此时得到的序列又是上面的序列向右循环移1位。依次类推，当n和m均从0变化到L-1时，得到一个关于x(nm)L的矩阵如下：,（3.2.6）,上面矩阵称为x(n)的L点“循环卷积矩阵

14、”，其特点是: （1） 第1行是序列x(0), x(1), , x(L1)的循环倒相序列。注意，如果x(n)的长度MN。若仍选取LNM1，以L为循环卷积区间，并用上述快速卷积法计算线性卷积，则要求对短序列补很多零点，而且长序列必须全部输入后才能进行快速计算。因此要求存储容量大，运算时间长，并使处理延时很大，不能实现实时处理。,况且在某些应用场合，序列长度不定或者认为是无限长，如电话系统中的语音信号和地震检测信号等。显然，在要求实时处理时，直接套用上述方法是不行的。解决这个问题的方法是将长序列分段计算，这种分段处理方法有重叠相加法和重叠保留法两种。下面只介绍重叠相加法，重叠保留法作为本章习题题2

15、1，留给读者讨论。 设序列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)的线性卷积时，可先计算分段线性卷积yk(n)=h(n)*xk(n)，然后把分段卷积结果叠加起来即可，如图3.4.3所示。每一分段卷积yk(n)的长度为NM1，因此相邻分段卷积yk(n)与yk1(n)有N1个点重叠，必须把重叠部分的yk(n)与yk1(n)相加，才能得到正确的卷积序列y(n)。,显然, 可用图3.4.1所示的快速卷积法计算分段卷积yk(n)， 其中

16、L=NM1。由图3.4.3可以看出，当第二个分段卷积y1(n)计算完后，叠加重叠点便可得输出序列y(n)的前2M个值；同样道理，分段卷积yi(n)计算完后，就可得到y(n)第i段的M个序列值。因此，这种方法不要求大的存储容量，且运算量和延时也大大减少，最大延时TDmax=2MTs+To，Ts是系统采样间隔，To是计算1个分段卷积所需时间，一般要求ToNM-1时，y(n)应补L-N个零点，而h(n)应从M-L到M-1区间上截取或按上述区间-N+1 M-1截取后在-N+1点前面补L-(N+M-1)个零点后，以L为周期进行周期延拓。,综上所述，可归纳出具体计算步骤如下： (1) 形成hL(n)序列：

17、 (2) (3),(4) Y(k)=DFTy(n) 0kL-1 (5) 计算Y(k)Y(k)； (6) V(k)=IDFTY(m)H(m) 0kL-1 (7) 与标准DFT(FFT)算法相比较，ChirpZ变换有以下特点： (1) 输入序列长度N和输出序列长度不需要相等，且二者均可为素数。,(2) 分析频率点zk的起始点z0及相邻两点的夹角 是任意的(即频率分辨率是任意的)，因此可从任意频率上开始，对输入数据进行窄带高分辨率的谱分析。 (3) 谱分析路径可以是螺旋形的。 (4) 当 时， zk均匀分布在单位圆上，此时ChirpCD*2Z变换就是序列的DFT。因此可以说，DFT是ChirpZ变换

18、的特例。 总之，ChirpZ变换用作谱分析时，具有灵活、适应性强和运算效率高等优点。,4 用DFT进行谱分析的误差问题 DFT(实际中用FFT计算)可用来对连续信号和数字信号进行谱分析。在实际分析过程中，要对连续信号采样和截断，有些非时限数据序列也要截断，由此可能引起分析误差。下面分别对可能产生误差的三种现象进行讨论。,(1) 混叠现象。对连续信号进行谱分析时，首先要对其采样，变成时域离散信号后才能用DFT(FFT)进行谱分析。采样速率Fs必须满足采样定理，否则会在=(对应模拟频率f=Fs/2)附近发生频谱混叠现象。这时用DFT分析的结果必然在f=Fs/2附近产生较大误差。因此，理论上必须满足

19、Fs2fc(fc为连续信号的最高频率)。对Fs确定的情况，一般在采样前进行预滤波，滤除高于折叠频率Fs/2的频率成分，以免发生频率混叠现象。 ,(2) 栅栏效应。我们知道，N点DFT是在频率区间0，2上对时域离散信号的频谱进行N点等间隔采样，而采样点之间的频谱函数是看不到的。这就好像从N个栅栏缝隙中观看信号的频谱情况，仅得到N个缝隙中看到的频谱函数值。因此称这种现象为栅栏效应。由于栅栏效应，有可能漏掉(挡住)大的频谱分量。,为了把原来被“栅栏”挡住的频谱分量检测出来，对有限长序列，可以在原序列尾部补零；对无限长序列，可以增大截取长度及DFT变换区间长度，从而使频域采样间隔变小，增加频域采样点数

20、和采样点位置，使原来漏掉的某些频谱分量被检测出来。对连续信号的谱分析，只要采样速率Fs足够高，且采样点数满足频率分辨率要求(见(3.4.14)式)，就可以认为DFT后所得离散谱的包络近似代表原信号的频谱。,(3) 截断效应。实际中遇到的序列x(n)可能是无限长的，用DFT对其进行谱分析时，必须将其截短，形成有限长序列y(n)=x(n)w(n)，w(n)称为窗函数，长度为N。w(n)=RN(n), 称为矩形窗函数。根据傅里叶变换的频域卷积定理，有,其中 对矩形窗数w(n)=RN(n)，有 幅度谱Wg()曲线如图3.4.12所示(Wg()以2为周期，只画低频部分)。图中，|2/N的部分称为主瓣，其

21、余部分称为旁瓣。,图3.4.12 矩形窗的幅度谱,例如，x(n)=cos(0n)，0=/4, 其频谱为 x(n)的频谱X(ej)如图3.4.13(a)所示。将x(n)截断后，y(n)=x(n)RN（n)的幅频曲线如图3.4.13(b)所示。,图3.4.13 x(n)=cos(0n)加矩形窗前、后的幅频特性,由上述可见，截断后序列的频谱Y(ej)与原序列频谱X(ej)必然有差别，这种差别对谱分析的影响主要表现在如下两个方面： (1) 泄露。 由图3.4.13(b)可知，原来序列x(n)的频谱是离散谱线，经截断后，使原来的离散谱线向附近展宽，通常称这种展宽为泄露。显然，泄露使频谱变模糊，使谱分辨率

22、降低。从图3.4.13可以看出，频谱泄露程度与窗函数幅度谱的主瓣宽度直接相关，在第7章将证明，在所有的窗函数中，矩形窗的主瓣是最窄的，但其旁瓣的幅度也最大。,(2) 谱间干扰。 在主谱线两边形成很多旁瓣，引起不同频率分量间的干扰(简称谱间干扰)，特别是强信号谱的旁瓣可能湮没弱信号的主谱线，或者把强信号谱的旁瓣误认为是另一频率的信号的谱线，从而造成假信号，这样就会使谱分析产生较大偏差。 ,由于上述两种影响是由对信号截断引起的，因此称之为截断效应。由图3.4.12可以看出，增加N可使Wg()的主瓣变窄，减小泄露，提高频率分辨率，但旁瓣的相对幅度并不减小。为了减小谱间干扰，应用其它形状的窗函数w(n

23、)代替矩形窗(窗函数将在FIR数字滤波其设计中介绍)。但在N一定时，旁瓣幅度越小的窗函数，其主瓣就越宽。所以，在DFT变换区间(即截取长度)N一定时，只能以降低谱分析分辨率为代价，换取谱间干扰的减小。 通过进一步学习数字信号处理的功率谱估计等现代谱估计内容可知，减小截断效应的最好方法是用近代谱估计的方法。但谱估计只适用于不需要相位信息的谱分析场合。 ,最后要说明的是，栅栏效应与频率分辨率是不同的两个概念。如果截取长度为N的一段数据序列，则可以在其后面补N个零，再进行2N点DFT，使栅栏宽度减半，从而减轻了栅栏效应。但是这种截短后补零的方法不能提高频率分辨率。因为截短已经使频谱变模糊，补零后仅使

24、采样间隔变小，但得到的频谱采样的包络仍是已经变模糊的频谱，所以频率分辨率没有提高。因此，要提高频率分辨率，就必须对原始信号截取的长度加长（对模拟信号，就是增加采样时间Tp的长度）。,习题与上机题 1 计算以下序列的N点DFT，在变换区间0nN1内，序列定义为 (1) x(n)=1 (2) x(n)=(n) (3) x(n)=(nn0) 0n0N (4) x(n)=Rm(n) 0mN (5),(6) (7) (8) x(n)=sin(0n)RN(n) (9) x(n)=cos(0n)RN(N) (10) x(n)=nRN(n),2 已知下列X(k)，求x(n)=IDFTX(k)： (1) (2)

25、 其中，m为正整数，0mN/2。,3 已知长度为N=10的两个有限长序列： 做图表示x1(n)、x2(n)、x1(n)与x2(n)的10点和20点循环卷积。,4 证明DFT的对称定理，即假设 X(k)=DFTx(n)，证明 DFTX(n)=Nx(Nk) 5 如果X(k)=DFTx(n)，证明DFT的初值定理 6 设x(n)的长度为N，且 X(k)=DFTx(n) 0kN1 令 h(n)=x(n)NRmN(n) m为自然数 H(k)=DFTh(n)mN 0kmN1 求H(k)与X(k)的关系式。,10 证明离散相关定理。若 则 11 证明离散帕塞瓦尔定理。若X(k)=DFTx(n)，则,12 已

26、知f(n)=x(n)+jy(n)，x(n)与y(n)均为长度为N的实序列。设 F(k)=DFTf(n)N 0kN1 (1) (2) F(k)=1+jN 试求X(k)=DFTx(n)N，Y(k)=DFTy(n)N以及x(n)和y(n)。,13 已知序列x(n)=anu(n)，0a1, 对x(n)的Z变换X(z)在单位圆上等间隔采样N点，采样序列为 求有限长序列IDFTX(k)N。 14 两个有限长序列x(n)和y(n)的零值区间为 x(n)=0 n0, 8n y(n)=0 n0, 20n 对每个序列作20点DFT，即,X(k)=DFTx(n) k=0, 1, , 19 Y(k)=DFTy(n)

27、k=0, 1, , 19 如果 F(k)=X(k)Y(k) k = 0，1，19 f(n)=IDFTF(k) k = 0，1，19 试问在哪些点上f(n)与x(n)*y(n)值相等, 为什么？,15 已知实序列x(n)的8点DFT的前5个值为0.25, 0.125j0.3018, 0, 0.125j0.0518, 0。 （1） 求X(k)的其余3点的值； （2） ，求 （3） x2(n)=x(n)ejn/4，求x2(k)=DFTx2(n)8。,16 x(n)、x1(n)和x2(n)分别如题16图(a)、(b)和(c)所示，已知X(k)=DFTx(n)8。 求 X1(k)=DFTx1(n)8 和

28、 X2(k)=DFTx2(n)8 注: 用X(k)表示X1(k)和X2(k)。 17 设x(n)是长度为N的因果序列，且 试确定Y(k)与X(ej)的关系式。,题16图,18 用微处理机对实数序列作谱分析，要求谱分辨率F50 Hz，信号最高频率为1kHz， 试确定以下各参数： (1) 最小记录时间Tp min； (2) 最大取样间隔Tmax； (3) 最少采样点数Nmin； (4) 在频带宽度不变的情况下，使频率分辨率提高1倍（即F缩小一半）的N值。,19 已知调幅信号的载波频率fc=1 kHz，调制信号频率fm=100 Hz，用FFT对其进行谱分析，试求： (1) 最小记录时间Tp; (2) 最低采样频率fs; (3) 最少采样点数N。 20 在下列说法中选择正确的结论。线性调频Z变换可以用来计算一个有限长序列h(n)在z平面实轴上诸点z