第3章
CHAPTER 3


离散傅里叶变换





计算机只能计算有限长离散序列,因此有限长序列在数字信号处理中就显得很重要,虽然可以用Z变换和序列傅里叶变换研究它,但是这两种变换无法直接利用计算机进行数值计算。针对序列“有限长”这一特点,可以导出一种更有效的变换——离散傅里叶变换(Discrete Fourier Transform,DFT)。作为有限长序列的一种表示方法,离散傅里叶变换除了在理论上相当重要之外,由于存在有效的快速算法——快速傅里叶变换(Fast Fourier Transform,FFT),因而在各种数字信号处理的算法中起着核心作用。本章主要讨论离散傅里叶级数和离散傅里叶变换定义性质,圆周卷积,利用DFT计算线性卷积和以及频域采样理论等内容。


视频讲解


3.1引言
傅里叶变换是以时间t,n为自变量的“信号”与以频率(Ω、f或ω )为自变量的“频谱”函数之间的一种变换关系。当自变量“时间”与“频率”为连续形式和离散形式的不同组合时,就形成了各种不同形式的傅里叶变换对,即“信号”与“频谱”的对应关系。
1. 连续时间、非周期信号的傅里叶变换
连续时间、非周期信号通过连续傅里叶变换(FT)得到非周期连续频谱密度函数。
正变换
Xa(jΩ)= ∫∞-∞xa(t)e-jΩtdt(31)

逆变换
xa(t)=12π∫∞-∞Xa(jΩ)ejΩtdΩ(32)



从图31的矩形脉冲及其频谱可以看出: 
(1) 时域连续函数造成频域是非周期的谱; 
(2) 时域的非周期造成频域是连续的谱。 


图31矩形脉冲及其频谱


2. 连续时间周期信号的傅里叶级数
周期为T的周期性连续时间函数xa(t)可展开成傅里叶级数Fn,是离散非周期性频谱。
正变换
Fn=1T∫T2-T2xa(t)e-jnΩtdt(33)

逆变换
xa(t)=∑∞n=-∞FnejnΩt(34)


由图32的周期矩形脉冲及其频谱可以看出: 
(1) 时域的连续函数造成频域是非周期的频谱函数; 
(2) 频域的离散频谱与时域的周期时间函数对应(频域采样,时域周期延拓)。


图32周期矩形脉冲及其频谱


3. 非周期离散信号的傅里叶变换
非周期离散的时间信号得到周期性连续的频率函数。
正变换
X(ejw)=∑∞n=-∞x(n)e-jωn(35)

逆变换
x(n)=12π∫π-πX(ejω)ejωndω(36)
由图33的非周期离散时间信号及其频谱可以看出
(1) 时域的离散造成频域的周期延拓; 
(2) 时域的非周期对应于频域的连续。


图33非周期离散时间信号及其频谱


4. 傅里叶变换的4种形式
通过上面分析,可以得出一般性规律: 一个域的离散对应另一个域的周期延拓,一个域的连续必定对应另一个域的非周期。上面的3种傅里叶变换对,都不适用在计算机上运算,因为至少在一个域(时域或频域)中函数是连续的。因为从数字计算角度,感兴趣的是时域及频域都是离散的情况,就是这里要谈到的离散傅里叶变换。
傅里叶变换的4种形式如下: 
(1) 连续时间、连续频率——傅里叶变换(FT),是非周期连续信号傅里叶变换,如图34(a)所示; 
(2) 连续时间、离散频率——傅里叶级数(FS); 是周期连续信号傅里叶级数,如图34(b)所示; 
(3) 离散时间、连续频率——序列傅里叶变换(DTFT); 是序列傅里叶变换,如图34(c)所示; 
(4) 离散时间、离散频率——离散傅里叶变换(DFT),是离散周期序列频谱,如图34(d)所示。


图34各种形式的傅里叶变换




视频讲解


3.2周期序列的离散傅里叶级数
3.2.1离散傅里叶级数定义

若离散时间序列x(n)为周期序列,则一定满足
x(n)=x(n+rN)
其中,N(正整数)为信号的周期,r为任意整数。为了与非周期序列区分,周期序列记作: x~(n)。 因为周期序列不是绝对可和,因此周期序列不能用傅里叶变换表示,但是周期序列可以用傅里叶级数(Discrete Fourier Series,DFS)表示,离散傅里叶级数定义为
x~(n)=1N∑N-1k=0X~(k)ej2πNnk=IDFS[X~(k)](37)

其中,X~(k)为周期序列傅里叶级数的系数,表示为
X~(k)=∑N-1n=0x~(n)e-j2πNnk=DFS[x~(n)](38)
为了书写方便,常令符号WN=e-j2πN,这样周期序列傅里叶变换对可以再次写为

正变换
X~(k)=DFS[x~(n)]=∑N-1n=0x~(n)WnkN

反变换
x~(n)=IDFS[X~(k)]=1N∑N-1k=0X~(k)W-nkN

MATLAB程序实现DFS过程如下: 

function [Xk]=dfs(xn,N)%计算(DFS)系数

%[Xk]=dfs(xn,N)

%Xk=在0<=n<=N-1之间的一个单周期信号

%N=xn的基本周期

n=[0:1:N-1];                %n的行向量

k=[0:1:N-1];                %k的行向量

WN=exp(-j*2*pi/N);         %Wn因子

nk=n'*k;                     %产生一个含nk值的N×N维矩阵

Xk=xn* WN.^nk;              %DFS 系数行向量

MATLAB程序实现IDFS过程如下: 

function [xn]=idfs(Xk,N)  %计算逆(DFS)系数

%[xn]=dfs(Xk,N)

% xn =周期信号在0<=n<=N-1之间的一个单周期信号

% Xk =在0<=k<=N-1间的DFS系数数组

%N=Xk的基本周期

n=[0:1:N-1];                 %n的行向量

k=[0:1:N-1];                 %k的行向量

WN=exp(-j*2*pi/N);          %Wn因子

nk=n'*k;                      %产生一个含nk值的N×N维矩阵

xn=Xk* (WN.^(-nk))/N;      %DFS 系数行向量

【例31】设x~(n)为周期脉冲串,求x~(n)=∑∞r=-∞δ(n+rN)。
【解】因为对于0≤n≤N-1,x~(n)=δ(n),x~(n)的DFS系数为
X~(k)=∑N-1n=0x~(n)WnkN=∑N-1n=0δ(n)WnkN=1


在这种情况下,对于所有k值,X~(k)均相同,于是有
x~(n)=∑∞r=-∞δ(n+rN)=1N∑N-1k=0W-nkN=1N∑N-1k=0ej2πNnk

当n为N整数倍时结果为1,这正好是周期性脉冲串。
说明: 
1N∑N-1n=0ej2πNkn=1N1-ej2πNkN1-ej2πNk=1,k=mN,m为整数

0,其他(39)

式(39)一般称为复正弦序列正交特性。
【例32】已知周期序列x~(n)如图35所示,其周期N=10,求解其傅里叶级数系数X~(k)。 


图35周期脉冲序列(N=10)


【解】
X~(k)=∑10-1n=0x~(n)Wnk10=∑4n=0e-j2π10nk

这一有限求和有闭合形式
X~(k)=1-W5k101-Wk10=e-j4πk10sin(5πk/10)sin(5πk/10)

周期脉冲序列DFS如图36所示,
MATLAB程序实现DFS过程如下: 

>>xn=[1,1,1,1,1,0,0,0,0,0];N=10;

>>Xk=dfs(xn,N)

Xk = Columns 1 through 5

5.0000 + 0.0000i1.0000 - 3.0777i-0.0000 + 0.0000i1.0000 - 0.7265i
-0.0000 + 0.0000i

Columns 6 through 10

1.0000 - 0.0000i  -0.0000 - 0.0000i1.0000 + 0.7265i-0.0000 - 0.0000i
1.0000 + 3.0777i



图36周期脉冲序列DFS


下面给出L=5,N=20周期方波的离散傅里叶级数的MATLAB程序: 

L=5;N=20;k=[-N/2:N/2];     %方波参数

xn=[ones(1,L),zeros(1,N-L)];             %方波x(n)

Xk=dfs(xn,N);                               %DFS

magXk=abs([Xk(N/2+1:N) Xk(1:N/2+1)]);  %DFS幅度

subplot(2,2,1);stem(k,magXk);

axis([-N/2,N/2,-0.5,5.5]);

xlabel('k');ylabel('Xtilde(k)');

title('方波的DFS:L=5,N=20')

上述程序产生的图及其他情形如图37所示,注意到X~(k)是周期信号,图中只画出了从-N/2到N/2的部分。从图37可以看出,方波的DFS系数包络看起来像sinc函数,k=0时的幅度为L,同时函数的零点位于N/L(占空比的倒数)的整数倍处。


图37不同L和N周期方波的离散傅里叶级数幅度



正变换定义公式(38)中的周期序列X~(k)可看成是对x~(n)的第一个周期x(n)作Z变换,然后将Z变换在Z平面单位圆上按等间隔角2π/N采样而


图38Z平面单位圆上等间隔采样


得到,如图38所示。令
x(n)= x~(n)·RN(n)=x~(n),0≤n≤N-1

0,其他(310)

通常称x(n)为x~(n)的主值区序列,则x(n)的Z变换为
X(z)=∑∞n=-∞x(n)z-n=∑N-1n=0x~(n)z-n(311)
X~(k)=X(z)z=W-kN=ej 2πNk(312)


由于单位圆上的Z变换为序列傅里叶变换,周期序列X~(k)也可以解释为x~(n)的一个周期x(n)的傅里叶变换的等间隔采样。因为
X(ejω)=∑N-1n=0x(n)e-jωn=∑N-1n=0x~(n)e-jωn(313)

比较式(38)和式(313),可以看出
X~(k)= X(ejω)ω=2πk/N=∑N-1n=0x~(n)e-j2πNkn,k=0,1,2,…,N-1(314)

这相当于以2π/N的频率间隔对傅里叶变换进行采样。也就是说,非周期离散时间信号经序列傅里叶变换(DTFT),得周期连续谱函数,再经采样得周期离散频谱函数(DFS),过程如图39所示。


图39序列傅里叶变换与离散傅里叶级数关系


由于频域N点取样使得频域离散从而形成时域序列的周期化。采样频点间隔为
ω0=2πN(315)

数字频率为
ω=kω0=2πNk,k=0,1,2,…,N-1(316)

【例33】傅里叶级数系数X~(k)和周期信号x~(n)的一个周期的傅里叶变换之间的关系。
【解】图310为周期序列x~(n),它的一个周期为
x(n)=1,0≤n≤4
0,其他

x~(n)的一个周期的傅里叶变换为
X(ejω)=∑4n=0e-jωn=1-e-j5ω1-e-jω=e-j2ωsin(5ω/2)sin(ω/2)



图310周期序列


根据|x(ejω)|绘制一个周期的DTFT幅度谱如图311所示。
根据上面分析,傅里叶级数为
X~(k)=X(ejω)ω=2πk/10=e-j4πk10sin(5πk/10)sin(πk/10)




图311x~(n)一个周期的DTFT幅度谱


根据X~(k)绘制周期序列傅里叶级数DFS如图312所示,可以看出X~(k)为|X(ejω)|在[0,2π]的N点等间隔采样。



图312x~(n)一个周期的DTFT幅度谱和DFS


3.2.2离散傅里叶级数的性质
由于可以用采样Z变换解释DFS,因此它的许多性质与Z变换性质非常相似。但是,由于x~(n)和X~(k)两者都具有周期性,这就使它与Z变换性质还有一些重要的差别。此外,DFS在时域和频域之间具有严格的对偶关系,这是序列的Z变换表示所不具有的。设x~1(n)和x~2(n)皆是周期为N的周期序列,各自的DFS分别为X~1(k)、X~2(k)。
1. 线性
DFS[ax~1(n)+bx~2(n)]=aX~1(k)+bX~2(k)(317)

2. 序列的移位
DFS[x~(n+m)]=W-mkNX~(k)=ej2πNmkX~(k)(318)
DFS[WnlNx~(n)]= X~(k+l)(319)

或
IDFS[X~(k+l)]=WnlNx~(n)=e-j2πNnlx~(n)(320)

3. 周期卷积
如果
Y~(k)= X~1(k)X~2(k)(321)

则
y~(n)=IDFS[Y~(k)]=∑N-1m=0x~1(m)x~2(n-m)(322)
或
y~(n)=∑N-1m=0x~2(m)x~1(n-m)(323)

两个周期都为N的周期序列x~1(n)和x~2(n),其卷积的结果也是周期为N的周期序列,求和只在一个周期上进行,即m为0到N-1,所以称为周期卷积。注意: n的取值不在N的范围内,得到结果后进行周期延拓,如图313所示。


图313两个周期序列(N=6)的周期卷积过程




视频讲解


3.3有限长序列离散傅里叶变换
3.3.1离散傅里叶变换定义

由于长度为N的有限长序列可以看作周期为N的周期序列的一个周期,利用DFS计算周期序列的一个周期,就可以得到有限长序列的离散傅里叶变换。设x(n)是长度为N的有限长序列,可以把它看作是周期为N的周期序列x(n)的一个主周期,而将x~(n)看作x(n)以N为周期进行周期延拓得到,即
x(n)=x~(n),0≤n≤N-1

0,其他

= x~(n)RN(n)(324)

同理
X~(k)=X((k))N(325)
X(k)= X~(k)RN(k)(326)
通常,把x~(n)的第一个周期n=0到n=N-1定义为主值区间,故x(n)是x~(n)的主值序列。而称x~(n)为x(n)的周期延拓(图314)。


图314周期序列、主值区间、主值序列之间关系


DFS与IDFS变换前后都是周期的、无限长的,但这里一个周期信息与其他周期信息相同,因而可以得到有限长序列的离散傅里叶变换的定义:
X(k)=DFT[x(n)]=∑N-1n=0x(n)WnkN,0≤k≤N-1

x(n)=IDFT[X(k)]=1N∑N-1k=0X(k)W-nkN,0≤n≤N-1(327)
(328)

注意,所处理的有限长序列都是作为周期序列的一个周期表示的。换句话说,离散傅里叶变换隐含周期性。 
【例34】已知序列x(n)=δ(n),求它的N点DFT。 
【解】单位脉冲序列的DFT很容易由DFT的定义式(327)得到  
X(k)=∑N-1n=0δ(n)WnkN=W0N=1,k=0,1,2,…,N-1

δ(n)的X(k)如图315所示。这是一个很特殊的例子,它表明对序列δ(n)来说,不论对它进行多少点的DFT,所得结果都是一个离散矩形序列。


图315δ(n)及其DFT


【例35】已知x(n)=cos(nπ/6)是一个长度N=12的有限长序列,求它的N点DFT。
【解】由DFT的定义式(327)可得 
X(k)=∑11n=0cosnπ6Wnk12=∑11n=012ejnπ6+e-jnπ6e-j2π12nk

=12∑11n=0e-j2π12n(k-1)+∑11n=0e-j2π12n(k+1)

利用复正弦序列正交特性式(39),则
∑N-1n=0ej2πN(k-m)n=N,k=m
0,k≠m

再考虑k的取值区间,可得
X(k)=6,k=1,11

0,其他
有限长余弦序列及其DFT如图316所示。


图316有限长余弦序列及其DFT



【例36】已知X(k)=3,k=0

1,1≤k≤9,求k=0,1,2,…,9点IDFT。 
【解】X(k)可以表示为X(k)=1+2δ(k),0≤k≤9。 
由于一个单位脉冲序列的DFT为常数
x1(n)=δ(n)X1(k)=DFT[x1(n)]=1
同样,一个常数的DFT是一个单位脉冲序列
x2(n)=1X2(k)=DFT[x2(n)]=Nδ(k)
X(k)=∑10-1n=0x(n)e-j2π10nk=1-e-j2πk1-e-j2π10k=10,k=0
0,k=1,2,…,9

=10δ(k)

由于
10δk→1
2δk→15

所以
x(n)=15+δ(n)
3.3.2DFT与序列傅里叶变换、Z变换的关系
若x(n)是一个有限长序列,长度为N,Z变换X(z)=∑N-1n=0x(n)z-n,则
X(z)z=W-kN=∑N-1n=0x(n)WnkN=DFT[x(n)]=X(k)(329)

z=W-kN=ej2πNk表明W-kN是Z平面单位圆上幅角为ω=2πNk的点,即将Z平面单位圆N等分后的第k点,如图317所示。所以X(k)是对X(z)在Z平面单位圆上N点等间隔采样值。此外,由于序列的傅里叶变换X(ejω)即是单位圆上的Z变换,根据式(329),DFT与序列傅里叶变换的关系为
X(k)=X(ejω)ω=2πNk=X(ejkωN)(330)

其中,ωN=2πN。
DFT的物理意义: 说明X(k)可看作序列x(n)的傅里叶变换X(ejω)在区间[0, 2π]上的N点等间隔采样,其采样间隔为ωN=2π/N。显而易见,DFT的变换区间长度N不同,表示对X(ejω)在区间[0, 2π]上的采样间隔和采样点数不同,所以DFT的变换结果也不同。 



图317DFT物理意义说明图示



【例37】有限长序列x(n)为
x(n)=1,0≤n≤4
0,其他
求其N=5点离散傅里叶变换X(k)。 
【解】将有限长序列x(n),如图318(a)所示,以N=5为周期将x(n)延拓成周期序列x~(n),如图318(b),x~(n)的DFSX~(k)与x(n)的DFT X(k)相对应如图318(c)和图318(d)所示。因为图318(b)中的序列在区间0≤n≤N-1上为常数值,所以可以得出



图318例37图


X~(k)=∑N-1n=0e-j(2πk/N)n=1-e-j2πk1-e-j(2πk/N)=N,k=0, ±N, ±2N, …
0,其他

在k=0和k=N的整数倍处才有非零的DFS系数X~(k)值。就是X(ejω)在频率ωk=2πk/N处的样本序列。x(n)的DFT对应于取X~(k)的一个周期而得到的有限长序列X(k)。
X(k)=∑5-1n=0x(n)e-j2π5nk=1-e-j2πk1-e-j2π5k=5,k=0
0,k=1, 2, 3, 4




视频讲解


3.4离散傅里叶变换的性质
本节讨论DFT的一些性质,它们本质上与周期序列的DFS概念有关,而且是由有限长序列及其DFT表示式隐含的周期性得出的。以下讨论的序列都是N点有限长序列,用DFT[·]表示N点DFT,且设DFT[x1(n)]=X1(k)、DFT[x2(n)]=X2(k)。
1. 线性
若两个有限长序列x1(n)和x2(n)的线性组合x3(n)=ax1(n)+bx2(n),则有
DFT[ax1(n)+bx2(n)]=aX1(k)+bX2(k)(331)

其中,a,b为任意常数。
说明: (1) 若x1(n)和x2(n)的长度均为N,则x3(n)的长度为N; 
(2) 若x1(n)和x2(n)的长度不等,x1(n)的长度为N1,x2(n)的长度为N2,则x3(n)的长度为N=max[N1,N2],离散傅里叶变换的长度必须按N计算。
2. 圆周移位
一个长度为N的有限长序列x(n)的圆周移位定义为y(n)=x((n+m))NRN(n),将图319(a)所示x(n)以N为周期进行周期延拓,得到周期序列x~(n)=x((n))N,如图319(b)所示; 再将x~(n)加以移位x((n+m))N=x~(n+m)如图319(c)所示,然后,再对移位的周期序列x~(n+m)取主值区间(n为0到N-1)上的序列值,即x((n+m))NRN(n)。所以,一个有限长序列x(n)的圆周移位序列y(n)仍然是一个长度为N的有限长序列。


图319圆周移位过程


3. 时域圆周移位定理
设x(n)是长度为N的有限长序列,y(n)为x(n)圆周移位,即
y(n)=x((n+m))NRN(n)(332)
则圆周移位后的DFT为
Y(k)=DFT[y(n)]=DFT[x((n+m))NRN(n)]=W-mkNX(k)(333)

4. 频域圆周移位定理
若X(k)=DFT[x(n)],则
IDFT[X((k+l))NRN(k)]=WnlNx(n)=e-j2πNnlx(n)(334)
上式称为频率移位定理,也称为调制定理,此定理说明时域序列的调制等效于频域的圆周移位。
5. 圆周卷积
设x1(n)和x2(n)都是点数为N的有限长序列(0≤n≤N-1),且有
Y(k)=X1(k)X2(k)
则
y(n)=IDFT[Y(k)]=∑N-1m=0x1(m)x2((n-m))NRN(n)

=∑N-1m=0x2(m)x1((n-m))NRN(n)(335)
卷积过程: 圆周卷积流程如图320所示,过程如图321所示。先将x1(n)和x2(n)补零,使得长度均为N点,并将变量n变成m,x2(m)周期化x2((m))N,再反转x2((-m))N,取主值序列x2((-m))NRN(m),对x2(m)圆周右移n,形成x2((n-m))NRN(m),当n=0,1,2,…,N-1时,分别将x1(m)与x2((n-m))NRN(m)相乘,并在m为0到N-1区间内求和,便得到圆周卷积y(n)。 


图320圆周卷积流程





图321圆周卷积过程示意图


特别要注意: 两个长度小于或等于N的序列的N点圆周卷积长度仍为N,这与一般的线性卷积不同。为了区别线性卷积,用*表示线性卷积,用
表示N点圆周卷积。 
y(n)=x1(n)x2(n)=∑N-1m=0x1(m)x2((n-m))NRN(n)(336)

利用时域与频域的对称性,可以证明频域圆周卷积定理。
若y(n)=x1(n)x2(n),
x1(n),x2(n)皆为N点有限长序列,则
Y(k)=DFT[y(n)]=1N∑N-1l=0X1(l)X2((k-l))NRN(k)
=1N∑N-1l=0X2(l)X1((k-l))NRN(k)=1NX1(k)X2(k)(337)
即时域序列相乘,乘积的DFT等于各个DFT的圆周卷积再乘以1N。 
6. 线性卷积与圆周卷积关系
时域圆周卷积在频域上相当于两序列的DFT的乘积,而计算DFT可采用它的快速算法——快速傅里叶变换(FFT),因此圆周卷积与线性卷积相比,计算速度可以大大加快。若序列xi(n)的长度为Ni(i=1,2),当N1与N2大小相当时,它们的线性卷积可以用L点DFT快速实现,其中 L≥N1+N2-1,如图322所示。


图322用圆周卷积和计算出线性卷积和的过程


实际问题大多总是要求解线性卷积。例如,信号通过线性时不变系统,其输出就是输入信号与系统的单位脉冲响应的线性卷积,如果信号以及系统的单位脉冲响应都是有限长序列,那么是否能用圆周卷积运算来代替线性卷积运算而不失真呢?下面就来讨论这个问题。

设x1(n)是N1点的有限长序列(0≤n≤N1-1),x2(n)是N2点的有限长序列(0≤n≤N2-1)。
1) 线性卷积
y1(n)=x1(n)x2(n)=∑∞m=-∞x1(m)x2(n-m)=∑N1-1m=0x1(m)x2(n-m)(338)
y1(n)是N1+N2-1点有限长序列,即线性卷积的长度等于参与卷积的两序列的长度之和减1。
2) 圆周卷积
先假设进行L点的圆周卷积,再讨论L取何值时,圆周卷积才能代表线性卷积。设y(n)=x1(n)x2(n)是两序列的L点圆周卷积,L≥max[N1, N2],这就要将x1(n)与x2(n)都看成是L点的序列。在这L点的序列值中,x1(n)只有前N1个是非零值,后L-N1个均为补充的零值。同样,x2(n)只有前N2个是非零值,后L-N2个均为补充的零值。则
y(n)=x1(n)x2(n)=∑L-1m=0x1(m)x2((n-m))LRL(n)(339)

可以证明
y~(n)=∑∞r=-∞y1(n+rL)(340)
y~(n)为x1(n)与x2(n)线性卷积的周期延拓,周期也为L,定义为周期卷积。
前面已经分析了y1(n)具有N1+N2-1个非零值。因此可以看到,如果周期卷积的周期L<N1+N2-1,那么y1(n)的周期延拓必然有一部分非零序列值会交叠起来,从而出现混叠现象。只有在L≥N1+N2-1时,才没有交叠现象,而圆周卷积正是周期卷积取主值序列。 
y(n)=x1(n)x2(n)= y~(n)RL(n)(341)
因此
y(n)=∑∞r=-∞y1(n+rL)RL(n)(342)
所以要使圆周卷积等于线性卷积而不产生混叠的必要条件为L≥N1+N2-1。线性卷积与圆周卷积关系如表31所示。


表31线性卷积与圆周卷积的关系



圆 周 卷 积
线 性 卷 积
1. 针对DFT引出的一种表示方法
1. 信号通过线性系统时,信号输出等于输入与系统单位脉冲响应的卷积
2. 两序列长度必须相等,若不等则按要求补零
2. 两序列长度可以相等,也可以不相等
3. 卷积结果长度与两信号长度相等皆为N
3. 卷积结果长度N=N1+N2-1

下面以具体示例说明,两个有限长矩形序列x1(n)与x2(n),它们长度分别为N1=4、N2=5,分别计算线性卷积和圆周卷积结果如图323所示。其中,图322(c)为线性卷积结果,图322(d)~(f)分别为长度为6、8、10的圆周卷积结果,根据L≥N1+N2-1关系,确实可以看到当L=8时线性卷积与圆周卷积结果相等。


图323线性卷积与圆周卷积示例


【例38】一个有限长序列为xn=δ(n)+2δ(n-5),
(1) 计算序列x(n)的10点离散傅里叶变换。
(2) 若序列y(n)的DFT为Yk=ej2k2π10Xk,式中,X(k)是x(n)的10点离散傅里叶变换,求序列y(n)。
(3) 若10点序列y(n)的10点离散傅里叶变换是Yk=XkWk,式中, X(k)是序列x(n)的10点DFT,W(k)是序列w(n)的10点DFT。
wn=1,0≤n≤6
0,其他
求序列y(n)。
【解】(1) 由式(327)可求得x(n)的10点DFT: 
X(k)=∑N-1n=0x(n)WnkN=∑10-1n=0[δ(n)+2δ(n-5)]Wnk10

=1+2W5k10=1+2e-j2π105k=1+2(-1)k

(2) 根据移位性质DFT[x((n+m))NRN(n)]=W-mkNX(k),X(k)乘以一个W-mkN形式的复指数相当于是x(n)圆周移位m点。本题中m=2, x(n)向左圆周移位了2点,就有
yn=xn+210R10n=2δ(n-3)+δ(n-8)

(3) X(k)乘以W(k)相当于x(n)与w(n)的圆周卷积。为了进行圆周卷积,可以先计算线性卷积再将结果周期延拓并取主值序列。x(n)与w(n)的线性卷积为
zn=xnwn=1,1,1,1,1,3,3,2,2,2,2,2

根据式(342),圆周卷积为
y(n)=∑∞r=-∞z(n+10r)R10(n)

在0≤n≤9求和中,仅有序列z(n)和z(n+10)有非零值,用表列出z(n)和z(n+10)的值,对n=0, 1, 2, …, 9求和,得到






n0123456789
1011

Z(n)
111113322222
z(n+10)2200000000
00

y(n)
3311133222
——

所以10点圆周卷积为: yn=3,3,1,1,1,3,3,2,2,2,由于6+7-1=12>
10所以线性卷积不等于圆周卷积。
【例39】已知x1(n)={1,1,1,1,0,2},x2(n)={0,1,2,1,0,3},n=0,1,…,5。求圆周卷积y1(n)=x1(n)⑥x2(n),y2(n)=x1(n)⑨x2(n),y3(n)=x1(n)x2(n)。
【解】首先计算线性卷积和yl(n)={0,1,3,4,4,6,6,7,5,0,6},n=0,1, …,10。yl(n)长度为11。
然后,将yl(n)以6为周期进行周期延拓得到y~(n)。由于周期延拓的周期小于yl(n)的长度,所以,在周期延拓时每个周期会有(N1+N2-1)-L个混叠点。当L=6,N1+N2-1=11,每个周期有5个混叠点,得到


y~(n)={…,6,8,8,4,10,6,…}
最后取y~(n)的主值序列即为圆周卷积和


y1(n)=x1(n)⑥x2(n)={6,8,8,4,10,6},n=0,1,…,5
同理,可以得到


y2(n)=x1(n)⑨x2(n)={0,7,3,4,4,6,6,7,5},n=0,1,…,8
当L≥N1+N2-1时,y(n)=yl(n),所以


y3(n)=x1(n)x2(n)=yl(n)={0,1,3,4,4,6,6,7,5,0,6},n=0,1,…,10

7.  DFT形式下的帕塞瓦尔定理
∑N-1n=0x(n)y(n)=1N∑N-1k=0X(k)Y(k)(343)

证明:
∑N-1n=0x(n)y(n)= ∑N-1n=0x(n)1N∑N-1k=0Y(k)W-knN

=1N∑N-1k=0Y(k)∑N-1n=0x(n)WknN

=1N∑N-1k=0X(k)Y(k)
如果令y(n)=x(n)
∑N-1n=0x(n)x*(n)=1N∑N-1k=0X(k)X*(k)
即
∑N-1n=0|x(n)|2=1N∑N-1k=0|X(k)|2

表明一个序列在时域计算的能量与在频域计算的能量是相等的。


视频讲解


3.5频域采样理论
考虑一个任意的绝对可和的非周期序列x(n),它的Z变换为
X(z)=∑∞n=-∞x(n)z-n
由于绝对可和,所以其傅里叶变换存在且连续,故Z变换收敛域包括单位圆。如果对X(z)在单位圆上进行N点等距采样:  
X(k)=X(z)z=W-kN=∑∞n=-∞x(n)WnkN,k=0,1,2,…,N-1(344)

问题在于,这样采样以后是否仍能不失真地恢复原序列x(n)。也就是说,频率采样后从X(k)的反变换中所获得的有限长序列,即xNn=IDFT[Xk],能不能代表原序列x(n)?为此,先来分析X(k)的周期延拓序列X~(k)的离散傅里叶级数的反变换,令其为x~N(n)。
x~N(n)=IDFS[X~(k)]=1N∑N-1k=0X~(k)W-nkN=1N∑N-1k=0X(k)W-nkN(345)

将式(327)代入式(345),可得
x~N(n)=1N∑N-1k=0∑∞m=-∞x(m)WmkNW-nkN=∑∞m=-∞x(m)1N∑N-1k=0W(m-n)kN(346)
由于
1N∑N-1k=0W(m-n)kN=1,m=n+rN,r为任意整数
0,其他(347)

所以
x~N(n)=∑∞r=-∞x(n+rN)(348)

这说明由X~(k)得到的周期序列x~N(n)是原非周期序列x(n)的周期延拓,其时域周期为频域采样点数N。在第1.2节中已经知道,时域采样造成频域的周期延拓,这里又看到一个对称的特性,即频域采样同样会造成时域的周期延拓。
(1) 如果x(n)是有限长序列,点数为M,则当频域采样不够密,即当N<M时,x(n)以N为周期进行延拓,就会造成混叠。这时,从x~N(n)就不能不失真地恢复出原信号x(n)。因此,频域采样不失真的条件是频域采样点数N大于或等于时域采样点数M(时域序列长度),即满足
N≥M
此时可得到
xN(n)= x~N(n)RN(n)=∑∞r=-∞x(n+rN)RN(n)=x(n),N≥M(349)
也就是说,点数为N(或小于N)的有限长序列,可以利用它的Z变换在单位圆上的N个等间隔点上的采样值精确地表示。
(2) 如果x(n)不是有限长序列(即无限长序列),时域周期延拓后,必然造成混叠现象,因而一定会产生误差; 当n越大时信号衰减得越快,或频域采样越密(即采样点数N越大),则误差越小,即xN(n)越接近x(n)。
【例310】已知一个序列x(n)为5点矩形序列,其序列傅里叶变换如图324所示,分析频域上的5点抽样。


图324矩形序列及其DTFT


【解】序列x(n)时域是有限长非周期的,所以频域是连续信号。现在频域上进行抽样处理,使其频域离散化。按N=5点进行频域抽样,由于频域抽样会造成时域延拓相加,时域延拓的周期个数等于频域的抽样点数N=5,由于N=M,所以时域延拓恰好无混叠现象,如图325所示。


图325频域抽样时域周期延拓(N=M)


按N=4时进行抽样,由于N=4,序列长度为M=5,N<M,时域延拓后产生混叠现象如图326所示。


图326频域抽样时域周期延拓(N<M)


既然N个频域采样X(k)能不失真地代表N点有限长序列x(n),那么这N个采样值X(k)也一定能够完全地表达整个X(z)及频率响应X(ejω)。讨论如下: 
X(z)=∑N-1n=0x(n)z-n

由于
x(n)=1N∑N-1k=0X(k)W-nkN(350)
X(z)=∑N-1n=01N∑N-1k=0X(k)W-nkNz-n=1N∑N-1k=0X(k)∑N-1n=0W-nkNz-n

=1N∑N-1k=0X(k)1-W-NkNz-N1-W-kNz-1(351)

由于W-NkN=1,因此
X(z)=1-z-NN∑N-1k=0X(k)1-W-kNz-1(352)

这就是用N个频率采样X(k)表示X(z)的内插公式。它可以表示为
X(z)=∑N-1k=0X(k)Φk(z)(353)
Φk(z)=1N1-z-N1-W-kNz-1(354)

式(354)称为内插函数。令其分子为零,得
z=ej2πNr,r=0,1,2,…,k,…,N-1(355)


图327内插函数的零极点分布

即内插函数在单位圆的N等分点上(也即采样点上)有N个零点。而分母为零,则有z=W-kN=ej2πNk的一个极点,它将与第k个零点相抵消。因而,插值函数Φk(z)只在本身采样点r=k处不为零,在其他(N-1)个采样点r上(r=0,1,2,…, N-1,r≠k)都是零点(有N-1个零点)。而它在z=0处还有(N-1)阶极点,如图327所示。 

现在讨论频率响应,即求单位圆上z=ejω的Z变换。
X(ejω)=∑N-1k=0X(k)Φk(ejω)(356)
将z=ejω代入式(354),Φk(ejω)可以表示成更方便的形式:  
Φk(ejω)=1N1-e-jωN1-e-jω-k2πN=1NsinωN2sinω-2πNk2e-jN-12ω+kπN

=1NsinNω2-πNksinω2-πNkejkπN(N-1)e-jN-12ω
这样
Φk(ejω)=Φω-k2πN(357)

其中,
Φ(ω)=1Nsin(ωN/2)sin(ω/2)e-jN-12ω(358)

式(358)称为频域响应的内插函数。而式(356)又可改写为
X(ejω)=∑N-1k=0X(k)Φω-2πNk(359)


Φω-k2πN满足以下关系 
Φω-k2πN=1,ω=k2πN=ωk
0,ω=i2πN=ωi,i≠k(360)

也就是说,函数Φω-k2πN在本采样点ωk=k2πN上,Φωk-k2πN=1,而在其他采样点ωi=i2πN,i≠k上,函数Φωi-k2πN=0。整个X(ejω)就是由N个Φω-k2πN函数分别乘上X(k)后求和的。所以很明显,在每个采样点上X(ejω)就精确地等于X(k)(因为其他点的插值函数在这一点上的值为零,没有影响)即
X(ejω)ω=2πNk=X(k),k=0,1,2,…,N-1(361)
就是说,各采样点之间的X(ejω)值由各采样点的加权插值函数X(k)Φω-2πNk在所求ω点上的值的叠加得到的,如图328所示。 


图328内插函数幅度特性与相位特性(N=5)



【例311】频域采样定理的验证,给定信号如下: 
xn=n+1,0≤n≤13
27-n,14≤n≤26
0,其他 
编写程序分别对频谱函数X(ejω)=FT[x(n)]在区间0,2π上等间隔采样32和16点,得到
X32(k)=X(ejω)ω=2π32k,k=0,1,2,…,31
X16(k)=X(ejω)ω=2π16k,k=0,1,2,…,15

再分别对X32k和X16k进行32点和16点IFFT,得到
x32(n)=IFFT[X32(k)]32,n=0,1,2,…,31
x16(n)=IFFT[X16(k)]16,n=0,1,2,…,15

分别画出X(ejω)、X32k和X16k的幅度谱,并绘图显示x(n)、x32n和x16n的波形,进行对比和分析,验证总结频域采样理论。
提示: 频域采样用以下方法编程实现。
(1) 直接调用MATLAB函数FFT计算X32(k)=FFT[x(n)]32得到X(ejω)在0,2π的32点频率域采样。
(2) 抽取X32(k)的偶数点即可得到X(ejω)在0,2π的16点频率域采样X16(k),即X16(k)=X32(2k),k=0,1,2,…,15。
(3) 当然也可以按照频域采样理论,先将信号x(n)以16为周期进行周期延拓,取其主值区(16点),再对其进行16点DFT(FFT),得到的就是X(ejω)在0,2π的16点频率域采样X16(k)。
MATLAB代码如下: 

M=27;N=32;N=0:M;

%产生M长三角波序列x(n)

xa=0:floor(M/2);  xb= ceil(M/2)-1:-1:0; xn=[xa,xb];

Xk=fft(xn,1024);  %1024点FFT[x(n)], 用于近似序列x(n)的TF

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)

subplot(3,2,2);stem(n,xn,'.');box on

title('(b)三角波序列x(n)');xlabel('n');ylabel('x(n)');

axis([0,32,0,20]);k=0:1023;wk=2*k/1024;

subplot(3,2,1);plot(wk,abs(Xk));title('(a)FT[x(n)]');

xlabel('\omega/\pi');ylabel('|X(e^j^\omega)|');axis([0,1,0,200])

k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),'.');box on

title('(c) 16点频域采样');xlabel('k');ylabel('|X_1_6(k)|');

axis([0,8,0,200]);n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,'.');box on

title('(d) 16点IDFT[X_1_6(k)]');xlabel('n');ylabel('x_1_6(n)');

axis([0,32,0,20]);k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),'.');box on

title('(e) 32点频域采样');xlabel('k');ylabel('|X_3_2(k)|');

axis([0,16,0,200]);n1=0:N-1;subplot(3,2,6);stem(n1,x32n,'.');box on

title('(f) 32点IDFT[X_3_2(k)]');xlabel('n');ylabel('x_3_2(n)');

axis([0,32,0,20])



function tstem(xn,yn)

%时域序列绘图函数

% xn:信号数据序列,yn:绘图信号的纵坐标名称(字符串)

n=0:length(xn)-1;

stem(n,xn,'.');box on

xlabel('n');ylabel(yn);

axis([0,n(end),min(xn),1.2*max(xn)])


结果如图329所示。图329(a)和(b)分别为X(ejω)和xn的波形; 图329(c)和(d)分别为X(ejω)的16点采样X16k和x16n=IDFTX16k波形; 图329(e)和(f)分别为X(ejω)的32点采样X32k和x32n=IDFTX32k波形图; 由于实序列DFT满足共轭对称性,因此频域图仅画出0,π上的幅频特性波形。本例中x(n)的长度M=26。从图中可以看出,当采样点数N=16<M时,x16n确实等于原三角序列xn以16为周期的周期延拓序列的主值序列。由于存在时域混叠失真,因而x16n≠xn; 当采样点数N=32>M时,无时域混叠失真,x32n=IDFTX32k=xn。


图329例311图



视频讲解


3.6MATLAB应用实例
【例312】x[k]=cos(2πrk/N),    N=16,    r=4,利用MATLAB计算16点序列x[k]的512点DFT,如图330所示。


图330x[k]的512点DFT


MATLAB代码如下: 

N = 16;k = 0:N-1;

L = 0:511;x = cos(2*pi*k*4./16);

X = fft(x);plot(k/16,abs(X),'o');

hold on;XE = fft(x,512);

plot(L/512,abs(XE))

xlabel('归一化频率');ylabel('幅度');


【例313】离散傅里叶变换X(k)=∑N-1n=0x(n)WnkN(矩阵相乘的方法)。
MATLAB代码如下: 

function [Xk]=dft(xn,N)

n=[0:1:N-1];

k=[0:1:N-1];

WN=exp(-j*2*pi/N);

nk=n'*k;

WNnk=WN.^(nk);

Xk=xn*WNnk;

【例314】逆离散傅里叶变换x(n)=1N∑N-1k=0X(k)W-nkN。
MATLAB代码如下: 

function [xn]=idft(Xk,N)

n=[0:1:N-1];

k=[0:1:N-1];

WN=exp(-j*2*pi/N);

nk=n'*k;

WNnk=WN.^(-nk);

【例315】信号的傅里叶分解与合成。
MATLAB代码如下: 

clear all

N = 256; dt = 0.05; % data numbers and sampling intervel,sampling frequence is 20Hz

n=0:N-1;t=n*dt; %序号序列和时间序列

x1=sin(2*pi*t);x2=0.5*sin(2*pi*5*t);

x=sin(2*pi*t)+0.5*sin(2*pi*5*t); %signals add

m=floor(N/2)+1; %down for integer

a=zeros(1,m);b=zeros(1,m);

for k=0:m-1

for ii=0:N-1

a(k+1)=a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N);%matlab's array index must be increase from 1

b(k+1)=b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N);

end

c(k+1)=sqrt(a(k+1).^2+b(k+1).^2);

end

if(mod(N,2)~=1)a(m)=a(m)/2;end

for ii=0:N-1

xx(ii+1)=a(1)/2;

for k=1:m-1;

xx(ii+1)=xx(ii+1)+a(k+1)*cos(2*pi*k*ii/N)+b(k+1)*sin(2*pi*k* ii/N);

end

end

subplot(2,2,1),plot(t, x1); title('正弦信号1'), xlabel('时间/s');

subplot(2,2,2),plot(t, x2); title('正弦信号2'), xlabel('时间/s');

%subplot(4,1,1),plot((0:N-1)*dt,xx);title('Composed sitgnal');

subplot(2,2,3),plot(t, x); title('合成信号'), xlabel('时间/s');

subplot(2,2,4),plot((0:m-1)/(N*dt),c),title('傅里叶变换'),xlabel('频率/Hz'), ylabel('幅度');

信号傅里叶合成与分解如图331所示。


图331例315图


此处的1Hz和5Hz的振幅与原来信号振幅不完全一致,是由于数据采样点较少导致的,即N较小。
【例316】补零序列的离散傅里叶变换。
MATLAB代码如下:

n=0:4; x=[ones(1,5)];   %产生矩形序列

k=0:999;w=(pi/500)*k;

X=x*(exp(-j*pi/500)).^(n'*k);  %计算离散时间傅里叶变换

Xe=abs(X);    %取模

subplot(3,2,1);stem(n,x);ylabel('x(n)'); %画出矩形序列

subplot(3,2,2);plot(w/pi,Xe);ylabel('|X(ejw)|'); %画出离散时间傅里叶变换

N=10;x=[ones(1,5),zeros(1,N-5)]; %将原序列补零为10长序列

n=0:1:N-1; X=dft(x,N);  %进行DFT

magX=abs(X); k=(0:length(magX)'-1)*N/length(magX);

subplot(3,2,3);stem(n,x);ylabel('x(n)');   %画出补零序列

subplot(3,2,4);stem(k,magX);  %画出DFT结果

axis([0,10,0,5]);ylabel('|X(k)|');

N=20;x=[ones(1,5),zeros(1,N-5)]; %将原序列补零为20长序列

n=0:1:N-1; X=dft(x,N);  %进行DFT

magX=abs(X); k=(0:length(magX)'-1)*N/length(magX);

subplot(3,2,5);stem(n,x);ylabel('x(n)');   %画出补零序列

subplot(3,2,6);stem(k,magX);  %画出DFT结果

axis([0,20,0,5]);ylabel('|X(k)|');

补零序列的离散傅里叶变换如图332所示。


图332补零序列的离散傅里叶变换



序列末端补零后,尽管信号的频谱不会变化,但对序列做补零后再做L点DFT,计算出的频谱实际上是原信号频谱在[0,2π)区间上L个等间隔采样,从而增加了对真实频谱采样的点数,并改变了采样点的位置,将会显示出原信号频谱更多的细节。故而数据后面补零可以克服栅栏效应。在第4章快速傅里叶变换中将详细介绍。
【本章习题】
3.1填空题
(1) 已知一个长度为N的序列xn,它的离散时间傅里叶变换为Xejω,它的N点离散傅里叶变换Xk是关于Xejω的点等间隔。
(2) DFT与DFS有密切关系,因为有限长序列可以看成周期序列的,而周期序列可以看成有限长序列的。
(3) 对长度为N的序列xn圆周移位m位得到的序列用xmn表示,其数学表达式为xmn= 。
(4) 设序列x(n)的N点DFT为X(k),则x((n+m))NRN(n)的N点DFT为。
(5) 某序列的DFT表达式为X(k)=∑N-1n=0x(n)WknN,由此可以看出,该序列时域的长度为,变换后数字频域上相邻两个频率样点的间隔是。
(6) 圆周卷积可被看作是周期卷积的; 圆周卷积的计算是在区间中进行的,而线性卷积不受此限制。
(7) 有限长序列Xz与Xk的关系,Xk与X(ejω)的关系。 
3.2选择题
(1) 序列x1n的长度为4,序列x2n的长度为3,则它们线性卷积的长度是(),5点圆周卷积的长度是()。

A. 5,5B.  6,5C.  6,6D.  7,5
(2) 下面描述中最适合离散傅里叶变换DFT的是()。
A. 时域为离散序列,频域也为离散序列
B. 时域为离散有限长序列,频域也为离散有限长序列
C.  时域为离散无限长序列,频域为连续周期信号
D. 时域为离散周期序列,频域也为离散周期序列
(3) 若序列的长度为M,要能够由频域抽样信号X(k)恢复原序列,而不发生时域混叠现象,则频域抽样点数N需满足的条件是()。
A. N≥MB. N≤MC. N≤2MD. N≥2M
3.3x(n)和h(n)是如下给定的有限序列x(n)={5,2,4,-1,2},h(n)={-3,2,-1}。
(1) 计算x(n)和h(n)的线性卷积y(n)= x(n)*h(n); 
(2) 计算x(n)和h(n)的6点循环卷积y1(n)= x(n)⑥h(n); 
(3) 计算x(n)和h(n)的8点循环卷积y2(n)= x(n)⑧h(n); 
比较以上结果,有何结论?
3.4证明W(N-n)kN=W-nkN=(WnkN)*。
3.5对有限长序列x(n)=1,0,1,1,0,1的Z变换X(z)在单位圆上进行5等分取样,得到取样值X(k),即X(k)=X(z)z=W-k5,k=0,1,2,3,4求X(k)的逆傅里叶变换x1(n)。
3.6试用定义计算周期为5,且一个周期内x(n)={2-,1,3,0,4}的序列x~(n)的DFS。
3.7设
x(n)=1,n=0,1

0,其他
将x(n)以4为周期进行周期延拓,形成周期序列x~(n),画出x(n)和x~(n)的波形,求出x~(n)的离散傅里叶级数X~(k)和傅里叶变换。
3.8已知序列x(n)=δ(n)+2δ(n-2)+δ(n-3),若y(n)是x(n)与其本身的4点循环卷
积,求y(n)及其4点DFT Y(k)。
3.9序列x(n)为
x(n)=2δ(n)+δ(n-1)+δ(n-3)
计算x(n)的5点DFT,然后对得到的序列求平方: 
Y(k)=X2(k)
求Y(k)的5点DFT反变换y(n)。
3.10设序列x(n)={1,2,3,2,1,0},v(n)={3,2,1,0,1,2}。
(1) 求x(n)的傅里叶变换X(ejω);
(2) 求v(k)=DFT[v(n)]6;
(3) 请解释v(k)与X(ejω)之间的关系。