第3章 离散傅里叶变换

1、第第3章章 离散傅立叶变换及其应用离散傅立叶变换及其应用 3.1引言 第2章讲到周期序列的傅里叶级数展开（DFS），它 揭示了离散的周期序列的频谱也是离散，但由于它们 是周期序列，无始无终，还是不能直接使用计算机进 行计算。如果把有限长序列（假设N点）进行周期延 拓。那么就可以借用DFS进行计算，并将最后结果截 取其一个周期；基于这样的思想提出了离散傅里叶变 换（DFT）。并于1965年由库利（Cooley）和图基 （Tukey）提出高效率计算DFT的快速方法，使得 DFT的各种应用突飞猛进，尤其是在卷积计算、调制、 解调、离散余弦变换等数字信号处理算法中。 3.2 离散傅立叶变换 3.2.1

2、 DFT定义 假定一个周期序列 ，它是由长为N点的有限长序列x(n) 经周期延拓而成，即 ( ) 01 ( )=, 0 x nnN x n % 其他 ( )=() r x nx nrN % 周期序列的主值区间主值区间与主值序列主值序列： 对于周期序列 ，定义： 第一个周期 n=0N-1，为 的“主值区间主值区间” 主值区间上的序列为主值序列为主值序列 x(n)。 x(n)与 的关系可描述为： )( nx )( nx )( nx )( )( )()( 主值序列的是 的周期延拓是 nxnx nxnx )()()()( )( )()( nRnxnRnxnx nxnx NNN N RN（n）为矩形序列

3、。 符号(n)N 是余数运算表达式余数运算表达式，表示 n 对 N 求余数。 例： 是周期为 N=8 的序列，求 n=11 和 n=-2 对 N的余数。 )( nx 6)2(68) 1(2 3)11(38111 8 8 n n )6()2( ),3()11( xxxx 频域上的主值区间与主值序列： 周期序列 的离散傅里叶级数 也是一 个周期序列，也可给它定义一个主值区间 ，以及主值序列 X(k)。 )( nx)( kX 10Nk N N kXkX kRkXkX )()( )()( )( 10)( )( )( 1 0 NkWnxnxDFSkX N n kn 10)( 1 )( )( 1 0 Nn

4、WkX N kXIDFSnx N n kn 再看周期序列的离散傅里叶级数变换（DFS）公式： 这两个公式的求和都只限于主值区间（0N-1）， 它们完全适用于主值序列 x(n) 与 X(k) ，因而我们可 得到一个新的定义 有限长序列离散傅里叶变换有限长序列离散傅里叶变换定义。 长度为N的有限长序列 x(n) ，其离散傅里叶变换 X(k) 仍是一个长度为N 的有限长序列，它们的关系 为： 10)( 1 )()( 10)()()( 1 0 1 0 NnWkX N kXIDFTnx NkWnxnxDFTkX N k kn N N n kn N x(n) 与 X(k) 是一个有限长序列离散傅里叶变换

5、对，已知 x(n) 就能唯一地确定 X(k) ，同样已知 X(k) 也就唯一地确定 x(n) 。 x(n) 与 X(k) 都是长度为 N 的序列（复序列） 都有N个独立值，因而具有等量的信息。 可以看出，k=0，1，2，N-1，相当于频率从直 流DC到(N-1)fs的N个频率等分点。这些频率点上的 余弦序列和正弦序列我们称之为频率单元，或分析 频点。也就是说，输入时域序列x(n)与频率单元做 序列点积运算而得到频谱的实部和虚部，即该频率 点所分解到的复系数X(k)。 111 22 000 ( )( )( )cos()( )sin()1 NNN nk NNN nnn X kx nWx nnkjx

6、 nnkkN = ;0 DFT具有隐含周期性，同样地，可以证明IDFT中，x(n)也是 隐含N点周期的。 11 () 00 ()( )( )( ) NN kmN nkn NN nn X kmNx n Wx n WX k DFT表示成矩阵形式，令 x=x(0)，x(1) ，x(2) ，x(N-1)T构成时域序 列的列矩阵。 X=X(0)，X(1) ，X(2) ，X(N-1) T构成频域 序列的列矩阵。 N X=W x 11 21(1) 22 22(1) N (1)2(1)(1) (1) 1111 1 1W 1 N NNN N NNN NNNN NNN WWW WWW WWW 1 N 1 x=WX

7、 N 【例3.2.1】某复合正弦信号 用Matlab计算出16个采样数据。 t=0:1/16000:15/16000; % 时间增量1/16ms xt=sin(2*pi*2000*t)+0.5*sin(2*pi*6000*t-3*pi/4); figure(1);stem(t,xt);xlabel(n);ylabel(x(n); ( )sin(22000 )0.5sin(260003/ 4)x ttt 构造1616的变换矩阵WN，并计算出频谱X(k)。 n=0:15;k=0:15; %两个行向量 WN=exp(-j*2*pi/16).(n*k); %构造变换矩阵 X=xt*WN;Xa=abs(

8、X);Xb=(angle(X)*180/pi; %弧度换成角度。 figure(2); subplot(2,1,1);stem(k,Xa);xlabel(k);ylabel(X(k); %幅度谱 subplot(2,1,2);stem(k,Xb);xlabel(k);ylabel(k); %相位谱 图中k=0和k=16是一样的，对应的是DC或fs(16KHz)； k=2和k=6对应的2KHz和6KHz的分量，且幅度也是成2倍关系。 但有问题：它们的幅值是8与4，对应相位变成了-90和135； 为何呢？ 这就是著名的DFT辅助效应，前者称为DFT的“计算增益”，增 益值0.5N, （负频率部分还

9、有0.5N），显然与数据点数有关， 数据越长，效应越大。这也是IDFT公式中除以N的原因。 后者叫DFT的“附加相位”，是个-90固定值。如6KHz分量的 相位本来是-135，计算出来却是135，这是因为-135- 90=-225，但在习惯的180相位主值表示方式中，- 225等价于135。 至于出现的k=10和k=14的高频分量，那是因为采样带来的镜像谐 波频谱。 幅度图中的细实线是从2KHz和6KHz复合正弦中截取一段2个周期 长的信号的连续频谱，DFT只看到了2个主瓣的最高点，扩散 的连续频谱的其他内容恰恰都躲过了DFT的观察点，即分析频 率单元都落在频谱的过零点上。 DFT性质：性质：

10、 假设有限长序列x1(n)和x2(n)，长度分别为N1和N2，取 N=maxN1, N2 （补0），它们的N点DFT分别为： X1(k)=DFTx1(n) X2(k)=DFTx2(n) (1) 线性线性 y(n)=ax1(n)+bx2(n) Y(k)=DFTax1(n)+bx2(n)=aX(k)+bX2(k) ，a，b为任 意常数 (2) 循环移位循环移位 有限长序列x(n)的循环移位循环移位定义为： y(n)=x(n+m)NRN(n) 含义： 1) x(n+m)N 表示 x(n) 的周期延拓序列的移位： )( )(mnxmnx N 2) x(n+m)NRN(n) 表示对移位的周期序列 x(n

11、+m)N 取主值序列，所以y(n)仍然是一个长度为N的有限长 序列。y（n）实际上可看作序列 x（n）排列在一个 N等分圆周上，并向左旋转向左旋转 m 位位。 循环移位 线性移位和循环移位操作比较 利用周期序列的移位特性： )( )(kXwmnxDFSmnxDFS mk NN ( )()( ) NN DFT y nDFT x nmRn )()()( )()( kXwkRmnxDFS nRmnxDFT mk NN N 序列循环移位后的DFT为： 同样，对于频域有限长序列X(k)的循环移 位，有如下反变换特性： 2 ()( )( )( ) jnr nr N NNN IDFT X krR kW x

12、nex n （3）共轭对称性）共轭对称性 设 x*（n）为 x（n）的共轭复数序列，则 DFTx*（n）=X*(N-k) 证： * 1 0 * )()( N n nk N WnxnxDFT * 1 0 )( N n nk N Wnx kn N kn N njkn N nN N j kn N nN N nkN N WWeWeWWW 2 2 )( )()( )()()( * * * 1 0 )(* kRkNX kNXWnxnxDFT NN N n nkN N 说明： 当k=0时，应为X*(N-0)=X*(0)，因为按定义X（k） 只有N个值，即0kN-1，而XN已超出主值区间， 但一般已习惯于把X

13、（k）认为是分布在N等分的圆周 上，它的末点就是它的起始点，即XN= X0，因此 仍采用习惯表示式 DFTx*（n）=X*（N-k） 以下在所有对称特性讨论中，XN均应理解为 XN=X0， 同样，x(N)=x(0)。 复序列的实部与虚部的DFT变换 以 Rex(n)和 Imx(n)表示序列x（n）的实部 与虚部 * 1 Re( ) ( ( )( ) 2 1 ( )()( ) 2 e DFTx nDFTx nxn X kXNkXk * 1 Im( ) ( ( )( ) 2 1 ( )()( ) 2 o DFT jx nDFTx nxn X kXNkXk Xe（k）和X0（K）表示实部与虚部序列的

14、DFT 由 Xe(k)可得： )()()(kXkXkX oe * )()( 2 1 )(kNNXkNXkNX e )()( 2 1 * kXkNX )()( * kNXkX ee 因此，Xe(k)具有共轭对称性，称为X(k)的共轭偶对 称分量。 显然， 用同样的方法可得到 X0(k)= - X*0(N-k) 即Xo(k)具有共轭反对称特性，称其为X(k)的共轭奇对 称分量。 对于纯实数序列 x（n），即x（n）=Rex(n)，X （k）只有共轭偶对称部分，即：X（k）=Xe（k）， 表明实数序列的DFT满足共轭对称性，利用这一特性 ，只要知道一半数目的X（k），就可得到另一半的 X （k），这

15、一特点在DFT运算中可以加以利用，以提高 运算效率。 根据DFT的对偶特性，分别以xe（n）及x0（n）表示序 列x（n）的循环共轭偶部与循环共轭奇部： )()( 2 1 )( )()( 2 1 )( * * nNxnxnjx nNxnxnx o e 同样应从循环意义上理解 x（N-0）=x（0）。 可证明: DFTxe（n）=ReX（k） DFTx0（n）=jImX（k） (4)选频特性选频特性 对复指数函数 进行采样得复序列 x（n） 0nN-1 其中r为整数。当0=2/N时，x（n）=ej2nr/N，其离散傅里 叶变换为 ( ) o jrn x ne 1 2/2/ 0 () N jn r

16、Njn kN n Xkee 2() 2()/ 1 01 jrk jrkN Nkr e kre ( ) o jrt xte 可见，当输入频率为r0时，变换X（K）的N个值中只 有 X（r）=N，其余皆为零，如果输入信号为若干个不 同频率的信号的组合，经离散傅里叶变换后，不同的k 上，X（k）将有一一对应的输出，因此，离散傅里叶 变换算法实质上对频率具有选择性频率具有选择性。 当0 2k/N 时? 下面分析 现在把【例3.2.1】中的2KHz频率分量改成2.3KHz ，且为了看得更清楚，去掉6KHz分量。即 x(t)=1.sin(2.2300.t)信号，经过T=1/16ms采样，用 N=16个数据

17、计算出的频谱幅度结果，如图3.2.7。 虚线是2.3KHz信号被截取后的连续频谱图。 虚线是2.3KHz信号被截取后的连续频谱图 一个无限长时域信号被截断后，将造成单一频率信 号的能量（频谱幅度平方），泄漏到附近所有频 率区域上。这称为频谱泄漏。 再看该X(k)所对应的采样数据x(n) ，已经不是原序列了。 如何才能避免泄漏呢？ 增大N，来满足条件。 比如：N=160个的，那么做DFT时，其频率分析 点间隔是fs/N=16KHz/160=0.1KHz，第23点就恰 好准确地观察到2.3KHz信号最高幅度。 （5）DFT与与Z变换变换的关系的关系 有限长序列可以进行z变换 1 0 )()()(

18、N n n znxnxZzX )()()()( 1 0 kXnxDFTwnxzX N n nk N wz k N k N wz zXkX )()( k N j k N ewz 2 是z平面单位圆上N等分后的第k点。 k N 2 图 DFT与z变换 N jk eXkX )( 1）X(k)也就是z变换在单位圆上等间隔的采样值。 2）X（k）也可看作是对序列付氏变换X（ej）的 采样，采样间隔为： N=2/N。 即 结论： 采样定律采样定律告诉我们，一个频带有限的信号，可以 对它进行时域采样而不丢失任何信息； DFT变换进一步告诉我们，对于时间有限的信号 （有限长序列），也可以对其进行频域采样频域采

19、样，而不 丢失任何信息，这正反应了傅立叶变换中时域、频域 的对称关系。 时域上的采样，使我们能够采用数字技术来处理这些 时域上的信号（序列）， DFT在频域也离散化，因此使得在频域采用数字技术 处理成为可能。 3.2.3频率域采样 时域采样间隔T，fs=1/T，为防止频谱混叠发生，fs应满 足时域采样定理。 同样，频域里连续频谱被采样成等间隔离散频率点，即 彼此呈谐波关系，而使得时域对应表现为周期化。 设时域序列点数为M，频率采样间隔F=2/M，那么对应 时域周期化的周期大小为ts=1/F，为防止时域混叠发 生，ts应足够大，大于信号长度。 设时域序列点数为Ns=ts/T，则只要M大于等于Ns

20、，那 就不会发生时域序列的混叠。这就是频率域采样定理 。 【例3.2.2】频率域取样的例子 一个序列的连续频谱在一周期里等间隔取样了 32个频率数据 经过IDFT逆变换后得到对应的时域序列： x(n)=2，-1，-3，-5，-2，2.1204e-016，1，2，4，5 ，7，9，8，6，3，1，1，2，1.2204e-016，0， 1.304e-016，0，0，0，0，0，0，0，0，0，0，0 。如图3.2.13，注意x(0)=1，x(1)=-2，x(16)=1， x(17)=2，从x(18)起都为0，说明x(n)是有限长N=18 点的时间序列。 如果改为等间隔取样16个频率数据。 对应的I

21、DFT后得到的时间序列x1(n)为16点， x1(n)=3，1，-3，-5，-2，2.1204e-016，1，2，4， 5，7，9，8，6，3，1 对比发现：x1(0)=3，x1(1)=1，与x(n) 不相同，其 他相同。因为序列x(n)实际长度18，只进行16点 的频域采样，将会发生时域周期化后的混叠。即 x1(0)=3=x(0)+x(-2)= x(0)+x(16)=2+1 ， x1(1)=1=x(1)+x(-1)= x(1)+x(17)=-1+2 。 其他14个数据不受影响。 如果是频谱取18个样点，那么将会刚好获得x(n)， 这是个频率取样密度的临界值。 变量周期分辨率 2 N 2 f、

22、 ss f、 N fs k N 频率 采样 模拟域 数字域 利用循环卷积和共轭对称特性，可证明DFT形式下的 Parseval定律： 1 0 1 0 * )()( 1 )()( N n N k kYkX N nynx 1 0 1 0 22 | )(| 1 | )(| N n N k kX N nx DFT形式下的形式下的Parseval定理定理 当y（n）= x（n）时，即为有限长序列的能量： 3.2.4循环卷积定理 1 0 )()()()()( N m NN nRmnymxkFIDFTnf 1 0 )()()()()( N m NN nRmnxmykFIDFTnf 若 F（k）=X（k）Y（

23、k） 证：这个卷积可看作是周期序列 卷积 后再取其主值序列。将F（k）周期延拓，得： )( )( nynx与 )( )( )( kYkXkF 1 0 1 0 )()()( )( )( N m NN N m mnymxmnymxnf )()()()()( )( 1 0 nRmnymxnRnfnf N N m NN 1 0 )()()()( N m NN nRmnymxnf 则根据DFS的周期卷积公式： 因0mN-1时，x(m)N=x(m)，因此 这一卷积过程与周期卷积比较，过程是一样的，只是 这里只取结果的主值序列，由于卷积过程只在主值区 间0mN-1内进行，所以 实际上就是 y （m）的循环移

24、位，称为“循环卷积循环卷积”，习惯上常用 符号“”表示循环卷积，以区别于线性卷积。 )()( )()()( )()()( )()( 1 0 1 0 nxny nRmnxmy nRmnymx nynx N m NN N m NN N mny)( 1）由有限长序列 x(n)、y(n) 构造周期序列 循环卷积过程: )( )( nynx与 2）计算周期卷积 1 0 )( )( )( N m mnymxnf 3）卷积结果取主值 )()( )(nRnfnf N 1 0 ( )( ) 1 ( ) ()( ) N NN l F kDFT f n X l Y k lR k N 1 0 )()()( 1 N l

25、 NN kRlkXlY N 同样，若 f(n)=x(n)y(n)，则 【例3.2.3】有两个长度都为6点的序列x(n)和y(n)，其 频谱分别记X(k)和Y(k)。验证DFT循环卷积性质。 x(n)=-2，5，-1，3，4，7和y(n)=1，2，7，3 ，4，6。 解：（1）把y(n)周期化； （2）对y(0)=1处左右翻转； （3）计算周期卷积 （4）取主值序列：f(n)= 75 ，68，50，82， 52，41。 Matlab程序如下： x=-2，5，-1，3，4，7; y=1，2，7，3，4，6; N=6； %序列循环卷积长度N m=0:1:N-1; y=y(mod(-m, N)+1);

26、 %对每个序号m求模6的值。即左右 翻转y序列。 A=zeros(N, N); %构造一个66的全0方阵。 for n=1:1:N A(n, : )=cirshftt(y, n-1, N); %对某个n，y序列循环移n- 1位后，对应放在A的第n行。 end f=x*A; %进行乘加运算，得到结果，f=75 ，68，50， 82，52，41。 figure(1); stem(f); %绘制序列杆图。 %N点循环移位函数cirshftt function w=cirshftt(s,m,N) %s是序列，N是其长度，m是 移位点数。 n=0:1:N-1; %得到序号0，1，2，3，4，5，.，N1

27、 q=mod(n-m,N); %根据位移量m值 w=s(q+1); %将循环移m位后的序列放函数出口w中。 %方法2: 计算x(n)和y(n)的DFT频谱序列X(k)和Y(k)。 X=fft(x);% Y=fft(y); % F=X.*Y; % 频域相乘 f1=ifft(F); %查看Workspace 有 f1=75 ，68，50，82，52 ，41，它确实和前面计算的一样！ figure(2); stem(f1); 有限长序列的线性卷积与循环卷积有限长序列的线性卷积与循环卷积 实际问题是求解线性卷积实际问题是求解线性卷积，如信号 x（n）通过系统 h（n），其输出就是线性卷积 y（n）=x

28、（n）*h（n）。 循环卷积比起线性卷积，在运算速度上有很大的优循环卷积比起线性卷积，在运算速度上有很大的优 越性越性，它可以采用快速傅里叶变换（FFT）技术， 问题提出问题提出: 能不能利用计算循环卷积的方法计算线性卷积? 对于上述 x（n）与h（n）的线性卷积，如果 x（n ）、h（n）为有限长序列，则实质上是研究在什么在什么 条件下条件下能用循环卷积代替而不产生失真。 有限长序列的线性卷积： 假定 x（n）为有限长序列，长度为N， y（n）为有限长序列，长度为M， 它们的线性卷积 m mnymxnynxnf)()()(*)()( 因 x（m）的非零区间： 0mN-1， y(n-m)的非零

29、区间： 0n-mM-1， 这两个不等式相加，得： 0nN+M-2， 在这区间以外不是x（m）=0，就是y（n-m）=0， 因而f（n）=0。因此，f（n）是一个长度为）是一个长度为N+M-1的的 有限长序列有限长序列。 循环卷积：循环卷积： 重新构造两个有限长序列 x（n）、y（n），长度 均为 L maxN,M (通过补充的零值)。 为了分析 x（n）与y（n）的循环卷积，先将x（n）， y（n）的周期延拓，得： r q rLnyny qLnxnx )()( )()( 11 00 ( )() ()() () LL L mm fnx m y nmx m y nm % % r r L m L m

30、r rLnf mrLnymx mrLnymx )( )()( )()( 1 0 1 0 结论： x（n）、y（n）周期延拓后的周期卷积，是x（n）、 y（n）线性卷积的周期延拓，周期为L。 它们的周期卷积序列为： 循环卷积正是周期卷积取主值序列： ( )( )( )( )( ) CLL fnx ny nfn Rn )()(nRrLnf L r 所以使循环卷积等于线性卷积而不产生混淆的循环卷积等于线性卷积而不产生混淆的必 要条件是： LN+M-1 【例3.2.4】上例数据，将y(n)看成是某离散系统的单位脉 冲响应，x(n)是其输入，那么系统的零状态响应就是二 者的线性卷积。调用函数conv（x

31、,y），计算两个序列 线性卷积， 即可得到6+6-1=11个点的输出响应数据。 f(n)=-2，1，-5，30，10，41，77，67，55，52，42 在x(n)后面添加5个0，使得序列成为11个点，即 x(n)=-2，5，-1，3，4，7，0，0，0，0，0； n=010。然后DFT求出X(k)，k=010。 同样，y(n)后面添5个0，再经过DFT得到Y(k)； 最后求IDFTX(k)Y(k)而得到输出响应f(n)= x(n)*y(n)。这结果与直接卷积conv（x,y）一样。从 而实现了用DFT求取系统响应的目的。 课外练习: 设有两个序列，x（n）为N=4矩形序列，y（n） 为M=6

32、矩形序列，计算其线性卷积和8点循环卷 积，并比较结果是否相同？为什么？ 3.2 利用利用DFT做连续信号的频谱分析做连续信号的频谱分析 采样采样截短截短DFT )( a jX a() x t( )xn )( j eX ( )( ) ( ) NN x n xnR n )( j N eX )(kX N Nk j N eX /2 )( （1）混迭 对连续信号x(t)进行数字处理前，要进行采 样 n a nTttxt x )()()( 采样序列的频谱是连续信号频谱的周期延拓， 周期为fs，如采样率过低，不满足采样定理， fs2fh，则导致频谱混迭，使一个周期内的谱对 原信号谱产生失真，无法恢复原信号，进一步 的数字处理失去依据。 （2） 泄漏信号截短造成的频谱扩散现象 处理实际信号序列 x（n）时，一般总要将它 截断为一有限长序列，设长为N点，相当于乘以一 个矩形窗 w(n)=RN(n)。 矩形窗函数，其频谱有主瓣，也有许多副瓣，窗口 越大，主瓣越窄，当窗口趋于无穷大时，就是一个 冲击函数。 我们知道，时域的乘积对应频域的卷积，所 以，加窗后的频谱实际是原信号频谱与矩形窗函 原信号频谱与矩形窗函 数频谱的卷积数频谱的卷积，卷积的结果使频谱延伸到了主瓣 以外，且一直延伸到无穷。当窗口无穷大时，与冲 击函数的卷积才是其本身，这时无畸变，否则