第3章〓离散傅里叶变换


第2章对序列的长度未加任何限制，既可是无限长序列又可是有限长序列。对于只能在0≤n≤N-1的有限范围内考虑的序列x(n)，特引入离散傅里叶变换(Discrete Fourier Transform，DFT)，它本身也是有限长度序列，而不是连续函数。离散傅里叶变换不仅在理论上有重要意义，而且有快速计算的方法——快速傅里叶变换(Fast Fourier Transform，FFT)。因而离散傅里叶变换在各种实现数字信号处理的方法中起着重要作用。

本章主要讨论傅里叶变换的4种形式； 离散傅里叶级数(DFS)推导及主要性质； 离散傅里叶变换的定义及性质； 频率采样； 用离散傅里叶变换对连续时间信号逼近的问题； 加权技术与窗函数。

从理论上说，离散傅里叶变换是傅里叶变换的一种可能形式。为了更好地理解这点，并不致发生混淆，本书先研究傅里叶变换的各种可能形式。




视频讲解


3.14种傅里叶变换◆

所谓傅里叶变换就是以时间为自变量的“信号”与以频率为自变量的“频谱”函数之间的某种变换关系。这种变换同样可以应用于其他有关物理或数学的各种问题中，并可以采用其他形式的变量。当自变量“时间”或“频率”取连续形式和离散形式的不同组合，就可以形成各种不同的傅里叶变换对，有些变换对是前文介绍过的，为前后一致，采用和前文相同的符号。讨论中所绘出的虚拟函数图只是为了清楚地说明各种特性，并不代表任何实际的变换对。

3.1.1非周期连续时间信号的傅里叶变换

非周期连续时间信号xa(t)的傅里叶变换Xa(jΩ)可以表示为


Xa(jΩ)=∫∞-∞xa(t)e-jΩtdt
(3.1.1)

逆变换为


xa(t)=12π∫∞-∞Xa(jΩ)ejΩtdΩ
(3.1.2)


式(3.1.1)及式(3.1.2)是大家所熟悉的非周期性连续时间信号及其频谱间的变换对，其非周期连续时间函数及其傅里叶变换的形式如图3.1.1所示。可以看到： 时域的连续函数造成频域的非周期谱，时域的非周期性造成频域的连续谱。结论是： 一个非周期连续时间函数对应一个非周期连续频率变换函数。




图3.1.1非周期连续时间函数及其傅里叶变换


3.1.2周期连续时间信号的傅里叶变换

周期为tp的周期性连续时间信号xa(t)的傅里叶变换是离散频率函数


X(mΩ)=1tp∫tp2-tp2xa(t)e-jmΩtdt
(3.1.3)

逆变换为


xa(t)=∑∞m=-∞X(mΩ)ejmΩt
(3.1.4)


式(3.1.3)及式(3.1.4)是被称为傅里叶级数的变换形式，式(3.1.3)所表示的积分是在xa(t)的一个周期内进行的。两相邻谱线分量之间的角频率增量与周期tp之间的关系可表示为


Ω=2π1tp=2πFF=1tp
(3.1.5)


式(3.1.3) 及式(3.1.4)两个函数的特性如图3.1.2所示。可以看到： 时域的连续函数造成频域的非周期谱； 频域函数的离散(采样)造成了时域函数的周期。结论是： 一个周期连续时间函数对应一个非周期离散频率变换函数。



图3.1.2周期连续时间函数及其傅里叶变换


3.1.3非周期离散时间信号的傅里叶变换

非周期离散时间信号x(n)的傅里叶变换X(ejω)是第2章讨论的序列及其频域表示的情况。其变换对为


X(ejω)=∑∞n=-∞x(n)e-jωn
(3.1.6)
x(n)=12π∫π-πX(ejω)ejωndω
(3.1.7)


式(3.1.6)及式(3.1.7)用数字域频率ω来表示变换对，并且式(3.1.7)是在X(ejω)的一个周期内求积分的。

采样频率fs与采样周期T的关系是


fs=1T


采样的角频率为Ωs=2πT,而采样数字域频率ωs=2π。

式(3.1.6)及式(3.1.7)两个函数的特性如图3.1.3所示。该变换对说明，时域的采样，对应频域函数的周期延拓(其周期在数字域频率恰为2π)。而时域函数的非周期对应频域函数的连续。



图3.1.3非周期离散时间信号及其傅里叶变换



3.1.4周期离散时间信号的傅里叶变换

按照时间变量和频率变量是连续还是离散的不同组合，可以推断，存在时间变量和频率变量都是离散的情况。时域采样会得到频域的周期性函数。由于在非周期连续傅里叶变换中时间t及频率f是对称的，所以在一对傅里叶变换式中将t与f对调之后，计算关系同样成立。因此在频域采样，将使时域信号得到周期延拓(关于频率采样理论下文还要介绍)。这样，第4种傅里叶变换对实际上是周期的离散时间信号与周期的离散频率信号间的变换对，如图3.1.4所示。这就是本书将要分析的离散傅里叶级数变换。



图3.1.4周期离散时间函数及其傅里叶变换



总结以上4种傅里叶变换对的形式，可得如下结论： 若一个函数在一个域内(时间域或频率域)是周期性的，则相应的在另一个域中的变换式必是采样的形式，即离散变量的函数； 反之，如果在一个域中的函数是采样的，则在另一个域中必是周期性函数。在一个域中函数的周期必是另一个域中两个采样点间增量的倒数。





视频讲解

3.2离散傅里叶级数◆

为了更好地理解离散傅里叶变换的概念，作为一种过渡，先简要地研究一下离散傅里叶级数。

3.2.1离散傅里叶级数变换的推导

当用数字计算机对信号进行频谱分析时，要求信号必须以离散值作为输入，而计算机输出所得的频谱值，自然也是离散的。因此在3.1节中介绍的傅里叶变换的可能形式中，只有第4种形式对于数字信号处理有实用价值。前3种形式中或者信号是时间的连续函数，或者频谱是频率的连续函数，或者信号及频谱二者都是变量的连续函数，因此都不适合数字计算机进行计算。要使前3种形式能用数字计算机进行计算，必须针对每种形式的具体情况，或者在时域与频域上同时采样，或者在时域上采样，或者在频域上采样。信号在时域上采样导致频率的周期函数，而在频域上采样导致时域的周期函数，最后都将使原时间函数和频率函数二者都成为周期离散的函数，即由于采样的结果，前3种形式最后都能变为第4种形式——离散傅里叶级数形式。现在以第3种傅里叶级数形式(见图3.1.3)为例来推导离散傅里叶级数变换。

为更清楚地表示式(3.1.6)所示的第3种形式傅里叶变换的周期性，下文在X(ejω)上加表示周期性的上标“~”，并重写如下


X～(ejω)=∑∞n=-∞x(n)e-jωn
(3.2.1)


设x～(n)的列长为N，则式(3.2.1)可写为


X～(ejω)=∑N-1n=0x(n)e-jωn


再对X～(ejω)采样，使其成为周期性离散频率函数，并导致时域序列x(n)周期化为x～(n)，如图3.1.4所示。图中时域采样间隔为T，在一个周期内采样点数为N。

由3.2节分析可知，在自变量为t及f的情况下,在一个域中对函数进行采样，两采样点间增量的倒数，必是另一个域中函数的周期。当序列的周期为NT时，对频谱采样的谱间距是1NT。以数字频率表示时，谱间距是


ωI=2πN
(3.2.2)

因此，以数字频率ω为变量的X(ejω)被离散化时，其变量ω成为


ω=kωI=2πNkk=0,1,…，N-1
(3.2.3)

所以离散周期序列x～(n)的傅里叶级数可写成


X～(k)= X～ejωω=2πNk=∑N-1n=0x～(n)e-j2πNkn，k=0,1,2，…，N-1
(3.2.4)

并将数字域频率简化为以k表示。上面两公式中，k为整数，而且由于X～(ejω)的周期是2π，所以k只可取0~N-1。这就是说，X～(k)有N个不同的值，X～(k)与x～(n)都是以N个采样值为一个周期的周期性函数。这一点在公式中看起来并不明显，事实上式(3.2.4)只能计算出N个独立的复数值。设k=r，其中r是任意整数，则由式(3.2.4)有


X～(r)=∑N-1n=0x～(n)e-j2πNrn

又设,k=r+N，注意到


e-j2πn(r+N)N=e-j2πnrNe-j2πn=e-j2πnrN

故有


X～(r+N)=∑N-1n=0x～(n)e-j2πN
nr= X～(r)

因此式(3.2.4)的X～(k)是以N个采样值为一个周期的周期性函数。

求离散傅里叶级数的逆变换，即从X～(k)求x～(n)。将正变换式(3.2.4)两边同乘以ej2πNkr，并对一个周期求和，即


∑N-1k=0X～(k)ej2πNkr=∑N-1k=0∑N-1k=0x～(n)e-j2πkn/N·ej2πNkr
=N∑N-1n=0x～(n)1N∑N-1k=0ej2πNkr(r-n)


根据正交定理


1N∑N-1k=0ej2πNk(r-n)=1，r=n
0，r≠n
(3.2.5)

则得


∑N-1k=0X～(k)ej2πNkr=N∑N-1n=0x～(n)n=r=Nx～(r)

以n替换r得


x～(n)=1N∑N-1k=0X～(k)ej2πNkn
(3.2.6)


与式(3.2.4)一样，式(3.2.6)所表达的也是一个以N为周期的周期序列


x～(n+mN)=1N∑N-1k=0X～(k)ej2πN(n+mN)k
=1N∑N-1k=0X～(k)ej2πNkn
= x～(n)


综上所述，离散周期序列的傅里叶级数变换对可表达为


X～(k)=DFS［x～(n)］=∑N-1n=0x～(n)e-j2πNkn
(3.2.7)
x～(n)=IDFS［X～(k)］=1N∑N-1k=0X～(k)ej2πNnk
(3.2.8)

其中，DFS［·］表示离散傅里叶级数正变换，IDFS［·］表示离散傅里叶级数逆变换。有时为了方便，令


WN=e-j2πN
(3.2.9)

并称之为WN因子，则式(3.2.7)及式(3.2.8)可表达为


X～(k)=DFS［x～(n)］=∑N-1n=0x～(n)WknN
(3.2.10)
x～(n)=IDFS［X～(k)］=1N∑N-1k=0X～(k)W-knN
(3.2.11)

3.2.2傅里叶级数的主要性质

离散傅里叶级数的某些性质对其成功地应用于信号处理问题是极其重要的。下面简要地介绍一下这些重要性质。

1. 线性特性

若有周期皆为N的两个离散周期序列x～1(n)和x～2(n)线性组合成一个新的周期序列x～3(n)


x～3(n)=ax～1(n)+bx～2(n)
(3.2.12)

则


X～3(k)=DFS［ax～1(n)+bx～2(n)］=aX～1(k)+bX～2(k)
(3.2.13)

式中，a、b为任意常数。线性特性可根据DFS的定义证明。

由于是线性组合，所以x～3(n)的周期长度不变，仍为N。X～3(k)也是周期为N的离散周期序列。

2. 序列移位

1) 时域移位

周期序列x～(n)左移m位后，得x～(n+m)，则


DFS［x～(n+m)］=W-mkNX～(k)
(3.2.14)


证明： 


DFS［x～(n+m)］=∑N-1n=0x～(n+m)WnkN换元i=n+m
=∑N-1+mi=mx～(i)WkiNW-mkN
=W-mkN∑N-1+mi=mx～(i)WkiN


由于x～(i)及WkiN都是以N为周期的周期函数，因此对i求和时，下限从m至上限N-1+m与下限从0至上限N-1是相同的。因此


∑N-1+mi=mx～(i)WkiN=∑N-1i=0x～(i)WkiN= X～(k)


所以


DFS［x～(n+m)］=W-mkNX～(k)


显然，大于周期的任何移位(也即m≥N)和短于周期的移位在时域上不能区分。

2) 频域移位

当将X～(k)左移l时，得X～(k+l)，则


IDFS［X～(k+l)］=WnlNx～(n)
(3.2.15)

可用与上文类似的方法证明式(3.2.15)。

3. 周期卷积特性

1) 时域卷积

若x～1(n)和x～2(n)是周期为N的两个周期序列，并分别具有离散傅里叶级数X～1(k)和X～2(k)，则傅里叶级数X～3(k)=X～1(k)·X～2(k)所对应的序列x～3(n)为


x～3(n)=∑N-1m=0x～1(m)x～2(n-m)=∑N-1m=0x～2(m)x～1(n-m)
(3.2.16)


证明： 


x～3(n)=IDFS［X～3(k)］=IDFS［X～1(k)·X～2(k)］
=1N∑N-1k=0X～1(k)X～2(k)W-nkN

将


X～1(k)=∑N-1m=0x～1(m)WmkN

代入，则


x～3(n)=1N∑N-1k=0∑N-1m=0x
～1(m)X～2(k)W-(n-m)kN
=∑N-1m=0x～1(m)1N∑N-1k=0X～2(k)W-(n-m)kN
=∑N-1m=0x～1(m)x～2(n-m)

只要将x～1(m)、x～2(n-m)简单换元，就可证明


x～3(n)=∑N-1m=0x～2(m)x～1(n-m)
(3.2.17)


式(3.2.17)是一个卷积公式，但不同于上文讨论过的线性卷积，二者的差别在于这里的卷积过程只限于一个周期以内，表现在式(3.2.17)中m取值取0~N-1，所以称为周期卷积。

图3.2.1给出了两个周期为N=7的序列x～1(n)和x～2(n)进行周期卷积的过程。



图3.2.1周期卷积


(1) 先在哑元坐标m上作x～1(m),x～2(m) (图中未画出)； 

(2) 将x～2(m)以m=0的垂直轴为轴作反卷得x～2(-m)，并进行移位，例如图3.2.1(d)、(e)、(f)所示的对应于n=0,1,2时的x～2(n-m)； 

(3) 在一个周期内将x～2(n-m)与x～1(m)对应点逐点相乘后求和，就得到相应于n=0,1,2时的x～3(n)，继续下去就得到了整个周期内的全部x～3(n)。

由于x～1(n)，x～2(n)都是周期为N的周期序列，因此x～3(n)也是周期为N的周期序列。

2) 频域卷积

对于周期序列的乘积，存在着频域的周期卷积。若


x～3(n)= x～1(n)x～2(n)

则


X～3(k)=DFS［x～3(n)］
=1N∑N-1l=0X～1(l)X～2(k-l)
=1N∑N-1l=0X～2(l)X～1(k-l)
(3.2.18)

根据DFS与IDFS变换的对称性不难证明式(3.2.18)。

【例3.1】计算6点的数字序列的DFS和卷积。

import numpy as np



N = 6

W6_1 = np.exp(-1j * 2 * np.pi / N * 1)

print("W6_1:", W6_1)

W6_2 = np.exp(-1j * 2 * np.pi / N * 2)

print("W6_2:", W6_2)



x = np.array(［14, 12, 10, 8, 6, 10］)

print("x:", x)

X = np.fft.fft(x)

print("X:", X)



# Define DFS function

def DFS(x, N, k):

X = np.zeros(N, dtype=complex)

for n in range(N):

X［k］ += x［n］ * np.exp(-1j * 2 * np.pi * k * n / N)

return X



X2 = DFS(x, 6, 0)

X3 = DFS(x, 6, 1)

print("X2:", X2)

print("X3:", X3)



x1 = np.array(［1, 1, 1, 1, 0, 0, 0］)

x2 = np.array(［0, 0, 1, 1, 1, 0, 0］)

x3 = np.convolve(x1, x2)

print("x3:", x3)


程序输出结果为

W6_1: (0.5000000000000001-0.8660254037844386j) 

W6_2: (-0.4999999999999998-0.8660254037844387j) 

x: ［14 12 10 8 6 10］ 

X: ［60.+0.00000000e+00j 9.-5.19615242e+00j 3.+1.73205081e+00j 0.+4.44089210e-16j 3.-1.73205081e+00j 9.+5.19615242e+00j］ 

X2: ［60.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j］ 

X3: ［0.+0.j 9.-5.19615242j 0.+0.j 0.+0.j 0.+0.j 0.+0.j］ 

x3: ［0 0 1 2 3 3 2 1 0 0 0 0 0］






视频讲解

3.3离散傅里叶变换◆

3.3.1DFT只有N个独立的复值

离散傅里叶级数是周期序列，仍不便于计算机计算。离散傅里叶级数虽是周期序列却只有N个独立的复值，只要知道它一个周期的内容，其他的内容也就知道了。式(3.2.10) 表明只要把一个周期内的x～(n)乘以对应的WnkN，这里的n在［0，N-1］内取值，而k可在(-∞，∞)内取值，可得任意k对应的X～(k)。式(3.2.11)表明，仅用X～(k)的一个周期的值，即k只取［0，N-1］区间内的值，就可得任意n，即n在(-∞，+∞)区间内取值时的x～(n)。如果同时限制式(3.2.10)中的k和式(3.2.11)中的n都只在区间［0,N-1］内取值，就得到了一个周期的x(n)和一个周期的X(k)间的对应关系。


X(k)=DFT［x(n)］=∑N-1n=0x(n)WknN，0≤k≤N-1
(3.3.1)
x(n)=IDFT［X(k)］=1N∑N-1k=0X(k)W-knN,0≤n≤N-1
(3.3.2)


式(3.3.1)和式(3.3.2)即称为有限长序列的离散傅里叶变换对。式(3.3.1) 称为离散傅里叶变换，简称DFT。式(3.3.2) 称为离散傅里叶逆变换，简称为IDFT。

3.3.2DFT隐含周期性

若长度为N的有限长序列x(n)是由对连续时间函数x(t)采样得来的，则频域上已经意味着以Ωs=2πT(或ω=2π)为周期作周期延拓。再对频域作一次采样，则时间序列x(n)按周期N 延拓成为一个周期性时间序列x～(n)。因此，利用DFT对有限列长N的时间序列展开，相当于对此序列作周期性处理。离散傅里叶变换对式(3.3.1)及式(3.3.2)，表面上看为非周期性的，但一定要想到它们隐含周期性。

3.3.3DFT是连续傅里叶变换的近似

离散傅里叶变换可以看成连续函数在时域、频域采样构成的变换。只要取出x～(n)的一个周期，乘以相应的内插函数就可恢复原连续函数。由下文分析可知，对频域也可取出X～(k)的一个周期，乘以相应的内插函数，就可恢复原连续的频率函数。DFT变换对可唯一地确定X～(k)的一个周期X(k)及x～(n)的一个周期x(n)。x(n)及X(k)都是长度为N的序列，都有N个独立复值，因而具有的信息是等量的，它们乘以相应的内插函数后，复原的连续函数也就完全确定了，因而离散傅里叶变换可以看作连续傅里叶变换的近似。x(n)及X(k)都是有限长序列，便于用数字计算机计算，这样对连续函数的处理就可以用离散采样的处理代替。这就是要采用离散傅里叶变换的原因。

【例3.2】求6点的数字序列的DFT和IDFT。

import matplotlib.pyplot as plt

import numpy as np

import math



# Number of sample points

N = 6

# sample spacing

T = 1.0 / 6.0

x = np.linspace(0.0, N*T, N, endpoint=False)

ys = x



# Build signal: Create a function to calculate the composite matrix M

def synthesis_matrix(N):

ts = np.arange(N) / N

fs = np.arange(N)

args = np.outer(ts, fs)

M = np.exp(1j * 2 * math.pi * args)

return M



# Define DFT positive transformation

def dft(ys):

N = len(ys)

M = synthesis_matrix(N)

amps = M.conj().transpose().dot(ys) # Calculate the weighted sum of frequency elements

return amps



# Define inverse DFT transform

def idft(ys):

N = len(ys)

M = synthesis_matrix(N)

amps = M.dot(ys) / N

return amps



print(dft(ys))

print(idft(ys))


程序输出结果为

［ 2.5+0.00000000e+00j -0.5+8.66025404e-01j -0.5+2.88675135e-01j

-0.5-3.06161700e-16j -0.5-2.88675135e-01j -0.5-8.66025404e-01j］

［ 0.41666667+0.00000000e+00j -0.08333333-1.44337567e-01j

-0.08333333-4.81125224e-02j -0.08333333+5.10269500e-17j

-0.08333333+4.81125224e-02j -0.08333333+1.44337567e-01j］



本程序中主要定义了3个函数： synthesis_matrix是用于构建信号，给定一系列振幅和频率，合成一个输入矩阵，dft和idft是根据DFT正逆变换的表达式对输入信号进行求值。

3.4离散傅里叶变换的性质◆

在下文的讨论中，假定x1(n)与x2(n)都是列长为N的有限长序列，它们的离散傅里叶变换分别为


X1(k)=DFT［x1(n)］
X2(k)=DFT［x2(n)］




视频讲解


3.4.1线性特性

若两个有限长序列x1(n)和x2(n)的线性组合为


x3(n)=ax1(n)+bx2(n)
(3.4.1)

则x3(n)的离散傅里叶变换为X1(k)与X2(k)的线性组合


X3(k)=DFT［ax1(n)+bx2(n)］
=aDFT［x1(n)］+bDFT［x2(n)］(3.4.2)
=aX1(k)+bX2(k)

式中a、b为任意常数。

注意，如果x1(n)列长为N1，x2(n)列长为N2，则x3(n)的最大列长为N3=max［N1,N2］。因而，离散傅里叶变换X3(k)必须按N=N3计算。例如，若N1<N2，则X1(k)就是序列x1(n)增补N2-N1个零点后的离散傅里叶变换。

3.4.2离散傅里叶逆变换的另一个公式

式(3.3.2)给出的离散傅里叶逆变换形式与正变换不同之处在于WN因子用负指数，且有一个比例系数1/N。离散傅里叶逆变换还有另一种形式，即


x(n)=1N∑N-1k=0X*(k)WnkN*0≤n≤N-1
(3.4.3)

证明： 


1N∑N-1k=0X*(k)WnkN*=1N∑N-1k=0X(k)W-nkN=x(n)
(3.4.4)

式(3.4.3)的逆变换公式中括号内WN因子用正指数，与正变换公式中WN指数一致。式(3.4.3)与正变换公式不同之处是用X(k)的共轭X*(k)来作运算，并在作出括号内的变换运算后，再取共轭并乘以常数1/N，得出时间序列x(n)。这种逆变换形式在运算上与正变换一样，这样在电子计算机上，只要编一个程序就可以既用来计算离散傅里叶变换，又用来计算它的逆变换。

【例3.3】采用正变换求6点的数字序列的DFT和IDFT。

import matplotlib.pyplot as plt

import numpy as np

import math



# Build signal: Create a function to calculate the composite matrix M

def synthesis_matrix(N):

ts = np.arange(N) / N

fs = np.arange(N)

args = np.outer(ts, fs)

M = np.exp(1j * 2 * math.pi * args)

return M



# Define DFT positive transformation

def dft(x):

N = len(x)

M = synthesis_matrix(N)

amps = M.conj().transpose().dot(x) # Calculate the weighted sum of frequency elements

return amps



# Define inverse DFT transform

def idft(X):

N = len(X)

M = synthesis_matrix(N)

amps = M.dot(X) / N

return amps



N = 6

x = np.array(［1, 2, 3, 4, 5, 6］)

print("x:", x)



X = dft(x)

print("X:", X)



X1 = np.conj (X)

print("X1:", X1)



X2 = np.conj (dft(X1)) * (1 / N)

print("X2:", X2)


程序输出结果为

x: ［1 2 3 4 5 6］

X: ［21.+0.00000000e+00j -3.+5.19615242e+00j -3.+1.73205081e+00j

 -3.-2.20436424e-15j -3.-1.73205081e+00j -3.-5.19615242e+00j］

X1: ［21.-0.00000000e+00j -3.-5.19615242e+00j -3.-1.73205081e+00j

 -3.+2.20436424e-15j -3.+1.73205081e+00j -3.+5.19615242e+00j］

X2: ［1.-4.29286236e-15j 2.+8.88178420e-16j 3.-1.11022302e-15j

 4.+2.36847579e-15j 5.+3.91517295e-15j 6.-4.58892184e-15j］



3.4.3对称定理

若x(n)的离散傅里叶变换为X(k)，则当时间序列具有频谱序列的形状X(n)时，其对应的离散傅里叶变换对如下


1NX(n)x(-k)=x(N-k)
(3.4.5)

式(3.4.5)说明X(n)的对应频谱序列具有原来的时间序列x(n)在时间上倒置的形状。

证明： 

因为


x(-n)=x(N-n)=1N∑N-1k=0X(k)W-(N-n)kN=1N∑N-1k=0X(k)WnkN

所以


x(-k)=x(N-k)=∑N-1n=0X(n)NWknN=DFTX(n)N


关于x(-n)=x(N-n)，X(-k)=X(N-k)的说明见下文离散傅里叶变换的奇偶性及对称性。

3.4.4反转定理

若x(n)的离散傅里叶变换为X(k)，则x(-n)的离散傅里叶变换为X(-k)。这可直接由离散傅里叶变换的定义得到证明。

3.4.5序列的总和

列长为N的时间序列x(n)中各采样值的总和等于其离散傅里叶变换X(k)在k=0时的值，即


X(k)k=0=∑N-1n=0x(n)WnkNk=0=∑N-1n=0x(n)
(3.4.6)

3.4.6序列的始值 

若序列的离散傅里叶变换为X(k)，则对应的时间序列x(n)的始值x(0)为频谱序列各采样值X(k)的总和除以N，即


x(0)=1N∑N-1k=0X(k)
(3.4.7)

3.4.7延长序列的离散傅里叶变换

把序列x(n),0≤n≤N-1，填充零值，人为地加长到rN，得到


g(n),0≤n≤rN-1

式中r为正整数，而


g(n)=x(n),0≤n≤N-1
0,N≤n≤rN-1
(3.4.8)

g(n)的离散傅里叶变换为


G(k)=DFT［g(n)］=∑rN-1n=0g(n)e-j2πnkrN
=∑N-1n=0x(n)e-j2πnkrN
=Xkrk=0,1,…,rN-1(3.4.9)

这意味着，g(n)的频谱G(k)与x(n)的频谱X(k)是相对应的，只不过G(k)的频谱间隔比X(k)的频谱间隔降低k/r。也就是说，若把序列x(n)填充零值而人为地加长以后再进行离散傅里叶变换，就可使得到的频谱更加细致。

若增加的长度并非N的整数倍，例如g(n)的长度为L>N，则由列长为L的序列g(n)的离散傅里叶变换G(k)，可得到序列x(n)的L根谱线，它比由X(k)得出的N根谱线要多。

【例3.4】计算对16点序列补零到32点和64点后的DFT。

import numpy as np

import matplotlib.pyplot as plt



N1 = 16

n = np.arange(0, 16)

x = np.sin(2*np.pi*n/N1) + np.cos(4*np.pi*n/N1)



N2 = 32

n2 = np.arange(0, 32)

x2 = np.concatenate((x, np.zeros(16)))



N3 = 64

n3 = np.arange(0, 64)

x3 = np.concatenate((x, np.zeros(48)))



L1 = np.arange(0, 16)

dft_16 = np.fft.fft(x, 16)

L2 = np.arange(0, 32)

dft_32 = np.fft.fft(x2, 32)

L3 = np.arange(0, 64)

dft_64 = np.fft.fft(x3, 64)



nx = np.arange(0, 16)

K = 512

dw = 2*np.pi/K

k = np.arange(0, 512)

X = np.dot(x, np.exp(1j*dw*nx.reshape(-1, 1)*k))



plt.subplot(4, 1, 1)

plt.plot(k*dw/(2*np.pi), np.abs(X))

plt.ylim(［0, 12］)



plt.subplot(4, 1, 2)

plt.stem(L1/N1, np.abs(dft_16))

plt.ylim(［0, 12］)



plt.subplot(4, 1, 3)

plt.stem(L2/32, np.abs(dft_32))

plt.ylim(［0, 12］)



plt.subplot(4, 1, 4)

plt.stem(L3/64, np.abs(dft_64))

plt.ylim(［0, 12］)



plt.show()



延长序列DFT如图3.4.1所示。



图3.4.1延长序列DFT





视频讲解

3.4.8序列的圆周移位
1. 圆周移位

一列长为N的有限长序列x(n)，在区间［0,N-1］内取非零值，如图3.4.2(a)所示。如令其沿坐标n右移m后仍在［0,N-1］区间内取值，如图3.4.2(b)所示，则将发生信息损失。为避免这种情况，可将x(n)以N为周期作周期延拓得


x～(n)=x((n))N

再把x～(n)偏移m得


x～(n+m)=x((n+m))N
(3.4.10)

上述过程如图3.4.2(c)，(d)所示。然后再对x～(n+m)仍在［0,N-1］区间上取值，即取移位周期序列x～(n+m)的主值序列x((n+m))NRN(n)，如图3.4.2(e)所示。



图3.4.2圆周移位过程图


这样，有限长序列x(n)的圆周移位定义为


x1(n)=x((n+m))NRN(n)
(3.4.11)


x1(n)仍是一长度为N的有限长序列，它实际上x(n)是向左移位，但m移到n负轴上的序列又从右边N-1处循环回来。这可想象为将序列x(n)排列在一个N等分的圆周上，且n=0的点取为正水平方向。由于列长为N，故相邻二点的间隔为2π/N角度，然后将圆顺时针旋转m个2π/N的角度，再从正水平方向开始依序取各点的序列值，即得x(n)的圆周移位序列x1(n)，因此得名为“圆周移位”，也称“循环移位”。其过程如图3.4.3所示。



图3.4.3序列的圆周移位


2. 有限长序列圆周移位定理

1) 时间移位定理

设X～(k)和X～1(k)分别表示周期序列x～(n)和x～1(n)=x～(n+m)的离散傅里叶级数。据式(3.2.14)有


X～1(k)=W-kmNX～(k)

由式(3.3.8)可得序列x(n)圆周移位所得的序列x1(n)的DFT为


X1(k)=W-kmN(k)X～(k)RN(k)=W-kmNX(k)
(3.4.12)

式(3.4.12)表示时间移位定理。它说明： 序列在时间上的时移将引起各根谱线产生相应的相移。

2) 频率移位定理(也称调制定理)

对于频域，有限长序列X(k)也可以认为是分布在一个N等分的圆周上。X(k)的圆周移位为X((k+l))NRN(k)，利用x(n)与X(k)的对称特性可证明。


IDFT［X((k+l))NRN(k)］=WnlNx(n)
(3.4.13)

式(3.4.13)即称为频率移位定理，有时也称为调制定理。此定理说明： 若频谱序列在频率上作移位l，则时间序列调制了lΩ的频率。这可以由下式看出


WnlNx(n)=e-j2πNnlx(n)=x(n)e-j2πNTnlT=x(n)e-jlΩnT
(3.4.14)


【例3.5】已知序列u(n)=［0 1 2 1 0］，求： 

(1) 循环右移2位； 

(2) 循环左移2位； 

(3) 序列反转； 

(4) 序列反转循环右移1位。

import numpy as np

import matplotlib.pyplot as plt



u = np.array(［0, 1, 2, 1, 0］)



# 循环右移2位

u1 = np.roll(u, 2)

plt.subplot(5, 1, 1)

plt.stem(u)

plt.xlabel('序号')

plt.ylabel('u(n)')



plt.subplot(5, 1, 2)

plt.stem(u1)

plt.xlabel('序号')

plt.ylabel('u1')



# 循环左移2位

u2 = np.roll(u, -2)

plt.subplot(5, 1, 3)

plt.stem(u2)

plt.xlabel('序号')

plt.ylabel('u2')



# 序列反转

u3 = np.flip(u)

plt.subplot(5, 1, 4)

plt.stem(u3)

plt.xlabel('序号')

plt.ylabel('u3')



# 序列反转循环右移1位

u4 = np.roll(np.flip(u), 1)

plt.subplot(5, 1, 5)

plt.stem(u4)

plt.xlabel('序号')

plt.ylabel('u4')



# 调整子图之间的间距

plt.subplots_adjust(hspace=1)



# 手动替换负号为正确显示负号

plt.rcParams［'axes.unicode_minus'］ = False

plt.xticks(fontproperties='SimHei')

plt.yticks(fontproperties='SimHei')



plt.tight_layout()

plt.show()



序列圆周移位如图3.4.4所示。



图3.4.4序列圆周移位






视频讲解

3.4.9圆周卷积及其与有限长序列的线性卷积关系


圆周卷积也称循环卷积。圆周卷积定理： 两个序列离散傅里叶变换的乘积等于此两个序列的圆周卷积的离散傅里叶变换，即若


X3(k)=X1(k)X2(k)

则


x3(n)=IDFT［X3(k)］=∑N-1m=0x1(m)x2((n-m))NRN(n)
=∑N-1m=0x2(m)x1((n-m))NRN(n)
(3.4.15)

证明： 

由于X1(k),X2(k)都隐含周期性，所以其乘积X3(k)同样隐含周期性。为此要证明式(3.4.15)可以应用周期卷积的结果。

式(3.4.15)的卷积可以看作是周期序列x～1(n)与x～2(n)卷积后再取其主值序列，即将X3(k)周期延拓


X～3(k)= X～1(k)X～2(k)

则


x～3(n)=IDFS［X～3(k)］=∑N-1m=0x～1(m)x～2(n-m)
=∑N-1m=0x1((m))Nx2((n-m))N

因0≤m≤N-1,x1((m))N=x1(m)，因此


x3(n)= x～3(n)RN(n)=∑N-1m=0x1(m)x2((n-m))NRN(n)
(3.4.16)

式(3.4.16)经过换元可得


x3(n)=∑N-1m=0x2(m)x1((n-m))NRN(n)


式(3.4.16)的卷积过程如图3.4.5所示，与图3.2.1的周期卷积相比，两者的卷积过程是一样的，只是这里只取卷积结果的主值序列。由于卷积过程只在主值区间内进行0≤m≤N-1，因此x2((n-m))NRN(n)实际上就是x2(m)的圆周移位，所以上述卷积称为圆周卷积，习惯上常用表示圆周卷积，以区别于线性卷积。式(3.4.16)的圆周卷积可表示为


x3(n)=x1(n)x2(n)
(3.4.17)



图3.4.5圆周卷积



必须指出，两列长为N的有限长序列的圆周卷积的结果也是一个列长为N的有限长序列，而不是像线性卷积那样列长为2N-1。

由于x3(n)与X3(k)的对称特性，若


x3(n)=x1(n)x2(n)

同样可以证明


X3(k)=DFT［x3(n)］
=1N∑N-1l=0X1(l)X2((k-l))NRN(k)
=1N∑N-1l=0X2(l)X1((k-l))NRN(k)(3.4.18)

式(3.4.18)称为频率卷积定理。

1. 圆周卷积和线性卷积的关系

实际问题求解线性卷积： 例如，信号x(n)通过系统h(n)其输出为线性卷积y(n)=x(n)*h(n)。下文将会讲到，圆周卷积由于可以采用快速傅里叶变换技术，所以比线性卷积运算速度快。如果x(n)及h(n)都是有限长序列，能否用圆周卷积来取代线性卷积而不失真呢？为此必须理解圆周卷积与线性卷积的关系，这是正确运用圆周卷积的关键。

假定x1(n)是列长为N的有限长序列，x2(n)是列长为M的有限长序列，二者的线性卷积为


x(n)=x1(n)*x2(n)=∑∞m=-∞x1(m)x2(n-m)

x(n)也是有限长序列，其列长可以这样来决定。

从x1(m)看，序号m的区间为


0≤m≤N-1


从x2(n-m)看，序号n-m的区间为


0≤n-m≤M-1

将上两个不等式相加，得


0≤n≤N+M-2

在上述区间外不是x1(m)=0就是x2(n-m)=0，因而x(n)=0。所以x(n)是一个n取0~N+M-2的序列，列长为N+M-1(两序列长度之和减1)。例如，图3.4.6中x1(n)为N=3的矩形序列，x2(n)为M=5 的矩形序列，两者的线性卷积x(n)是具有N+M-1=7的采样值的序列。



图3.4.6圆周卷积与线性卷积


对x1(n)和x2(n)进行列长为L的圆周卷积XL(n)=x1(n)x2(n)的情况。把列长为N的序列x1(n)后面补上L-N个零值点，把列长为M的序列x2(n)后面补上L-M个零值点，均构成列长为L的有限长序列。为了进行圆周卷积，先将x1(n)及x2(n)进行周期延拓


x～1(n)=∑∞q=-∞x1(n+qL)
x～2(n)=∑∞k=-∞x2(n+qL)

则


x～2(-m)=∑∞k=-∞x2(-m+kL)
x～2(n-m)=∑∞k=-∞x2(n-m+kL)

x～1(n)及x～2(n)的列长为L的周期卷积为


x～L(n)=∑L-1m=0∑∞q=-∞x1(m+qL)∑∞k=-∞x2(n-m+kL)
=∑L-1m=0∑∞k=-∞x1(m)x2(n-m+kL)


交换求和次序，有


x～L(n)=∑∞k=-∞∑L-1m=0x1(m)x2(n+kL-m)
=∑∞k=-∞x(n+kL)
(3.4.19)

式(3.4.19)表明： x1(n)和x2(n)的周期卷积是x1(n)和x2(n)的线性卷积x(n)以L为周期的延拓。

由前文分析可知，线性卷积x(n)具有N+M-1个非零序列值。因此，如果周期卷积的周期L<N+M-1,则x(n)的周期延拓就必然有一部分非零序列值要交叠起来，发生混淆。只有当L≥N+M-1时，才不会发生交叠，这时x(n)的周期延拓x～1(n)中的每一个周期L内，前N+M-1个序列值是序列x(n)的值，而剩下的L-(N+M-1)个点上序列值则是补充的零值。

圆周卷积就是周期卷积取主值序列，即


xL(n)=x1(n)x2(n)= x～L(n)RL(n)

因此


xL(n)=∑∞k=-∞x(n+kL)RL(n)

所以要使圆周卷积与线性卷积相等，而不产生混淆的必要条件是


L≥N+M-1
(3.4.20)

符合式(3.4.20)的条件，则xL(n)=x(n)，即


x1(n)x2(n)=x1(n)*x2(n)


图3.4.7(d)~图3.4.7(f)反映了当L取值不同时，线性卷积与圆周卷积的关系。图3.4.7(d)中L=5,这时混淆现象严重，使得xL(n)与x(n)完全不一样，而图3.4.7(e)中L=6,这时有两点(n=0,n=6)发生混淆失真，


图3.4.7圆周位移及线性位移比较



在图3.4.7(f)中，满足条件L≥N+M-1=7时，圆周卷积与线性卷积的结果相同。只要具体比较一下在这种条件下，圆周卷积与线性卷积的实际过程就会透彻地理解两者为什么相等。

如图3.4.7所示，L满足L≥N+M-1=7的条件，列长为N=3的序列x(n)已补上L-N=4个零值点(图3.4.7(a))； 列长为M=5的序列x2(n)已补充上L-M=2个零值点(图3.4.7(b))，则在L长的区间上对两个序列作圆周卷积时，由于x2(-m)有L-M个零值点，所以它在［0，L-1］区间上作圆周移位时，在［0，L-1］区间上移位的结果和x2(-m)的线性移位结果是一样的，如图3.4.7(c)，(d)，(e)，(f)，(g)，(h)所示。另一方面，虽然在［N，L-1］区间上圆周移位结果和线性移位结果可能会不同，但由于x1(n)在此区间均为零值点，所以也不会有输出。因而在区间［0，L-1］上的圆周卷积结果和x1(n)，x2(n)的线性卷积结果实际是一样的。

综上所述可以清楚看出： 有限长序列x1(n)和x2(n)的周期卷积是x1(n)和x2(n)的线性卷积以周期延拓。而圆周卷积是周期卷积的去主值序列，当满足条件L≥N+M-1时，圆周卷积与线性卷积的结果相同。

【例3.6】(1) 求序列x1=［1, 2, 3］和序列x2=［2, 3, 1, 2］的线性卷积。

import numpy as np

import matplotlib.pyplot as plt



def conv_m(x, h):

nx = len(x)

nh = len(h)

ny = nx + nh - 1

y = np.zeros(ny)# 初始化

for n in np.arange(ny):

y［n］ = 0

for m in np.arange(nh):

k = n - m + 1

if k >=1 and k <=nx:

y［n］ = y［n］ + h［m］ * x［k - 1］# 卷积



plt.rcParams［'font.sans-serif'］ = ［'SimHei'］

plt.rcParams［'axes.unicode_minus'］ = False# 显示中文标签

plt.stem(y)

plt.title('线性卷积')

plt.xlabel('n')

plt.ylabel('y')

plt.xlim(0, 8)

plt.show()



conv_m(np.array(［1, 2, 3］), np.array(［2, 3, 1, 2］))




运行程序，输出结果如图3.4.8所示。



图3.4.8有限长序列的线性卷积


(2) 求序列x1=［1, 2, 3］和序列x2=［2, 3, 1, 2］在N=5, 6, 7点时的圆周卷积。

import numpy as np

import matplotlib.pyplot as plt



def conv_m(x, h):# 线性卷积

nx = len(x)

nh = len(h)

ny = nx + nh - 1

y = np.zeros(ny)# 初始化

for n in np.arange(ny):

y［n］ = 0

for m in np.arange(nh):

k = n - m + 1

if k >=1 and k <=nx:

y［n］ = y［n］ + h［m］ * x［k - 1］# 卷积

return y



def cir_con(x1, x2, N):# 圆周卷积

nx1 = np.arange(0, len(x1))

nx2 = np.arange(0, len(x2))

x_1 = np.append(x1, np.zeros(N - len(x1)))

h_1 = np.append(x2, np.zeros(N - len(x2)))

y1 = conv_m(x_1, h_1)# 调用线性卷积函数

z_1 = np.append(np.zeros(N), y1［0:N］)

z_2 = np.append(y1［N:2*N］, np.zeros(N))

z = z_1［0:2*N-1］ + z_2［0:2*N-1］ + y1［0:2*N-1］

y = z［0:N］

ny = np.arange(0, N)



plt.rcParams［'font.sans-serif'］ = ［'SimHei'］

plt.rcParams［'axes.unicode_minus'］ = False# 显示中文标签

plt.stem(y)

plt.title('圆周卷积')

plt.xlabel('n')

plt.ylabel('y')

plt.xlim(0, 8)

plt.show()



cir_con(np.array(［1, 2, 3］), np.array(［2, 3, 1, 2］), 5)

cir_con(np.array(［1, 2, 3］), np.array(［2, 3, 1, 2］), 6)

cir_con(np.array(［1, 2, 3］), np.array(［2, 3, 1, 2］), 7)



运行程序，结果如图3.4.9所示。



图3.4.9有限长序列的圆周卷积


2. 圆周卷积在信号处理中的应用——重叠相加法与重叠保留法


可用两种办法计算线性系统的输出序列。

(1) 按线性卷积来计算； 

(2) 通过圆周卷积，用快速卷积法计算。

设输入x(n)的列长为N，系统的单位采样响应的列长为M，则二者线性卷积后的输出序列y(n)的列长为L=N+M-1。由圆周卷积定理可知，两序列离散傅里叶变换的乘积等于它们圆周卷积的离散傅里叶变换。分别取x(n)和h(n)的L≥N+M-1点的离散傅里叶变换X(k)和H(k)，将两者相乘再求其逆变换即得两序列的圆周卷积，此时圆周卷积与线性卷积相等。这种方法称为快速卷积法，因为可以利用快速傅里叶变换迅速而有效地求出卷积，甚至当N+M-1为30的时候，快速卷积法也比直接卷积法更有效。快速卷积法是一种重要的信号处理工具。

当输入序列x(n)极长时，为了尽快得到输出，不允许等x(n)全部集齐后再进行卷积,否则会使输出相对于输入有较长的延时。再者，如果N+M-1太大，则需要大量存储单元。为此，可把x(n)分段，分别求出每段的卷积，合在一起即得总的输出。这种分组进行卷积法可利用上文介绍的离散傅里叶变换实现，并可细分为重叠相加法、重叠保留法两种。下面针对图3.4.10所描绘的信号x(n)和列长为M的单位采样响应h(n)来解释这两种方法。



图3.4.10有限列长的单位采样响应及待处理的信号


1) 重叠相加法

h(n)的列长为M，信号为x(n)。将x(n)分解为几段之和，每段长为L点，如图3.4.10所示。L的选择是相当复杂的，一个好的经验是使L与M的数量级相同。


以xk(n)表示第k段x(n)，于是输入信号x(n)可表示成


x(n)=∑∞k=0xk(n)
(3.4.21)


式中


xk(n)=x(n)，kL≤n≤(k+1)L-1
0,其他
(3.4.22)

x(n)与h(n)的线性卷积等于xk(n)和h(n)卷积之和，即


x(n)*h(n)=∑∞k=0xk(n)*h(n)
(3.4.23)

式(3.4.23)和式中的每一项［xk(n)*h(n)］的列长为L+M-1，因此可将h(n)及xk(n)两个序列补0，都加长到L+M-1点，以利用L+M-1的DFT，通过圆周卷积得到线性卷积xk(n)*h(n)。为便于用基2FFT运算，一般可取L+M-1=2v。由于每一段的起点和前后相邻各段起点相隔L个点，而每一段信号xk(n)和h(n)卷积后列长为(L+M-1)，所以式(3.4.23)在求和时，每段卷积的最后M-1个点必须和下一段的前M-1个点重叠，如图3.4.11所示。各输入段xk(n)描绘于图3.4.11(a)，相加各xk(n)波形即可重新组成输入信号x(n)的波形。卷积后各段输出如图3.4.11(b)所示，而总的输出序列y(n)=x(n)*h(n)是通过相加图3.4.11(b)中各段构成的。重叠部分也要相加，故称为重叠相加法。重叠是由于每段输入序列xk(n)与单位采样响应h(n)的线性卷积后的列长长于x(n)的分段长度造成的。



图3.4.11重叠相加法图形


2) 重叠保留法

重叠保留法延长分段序列的办法不是补零，而是保留原来的输入序列值，且保留在每段的前端，如图3.4.12(b)中的虚线部分所示。这时如利用DFT实现h(n)和xk(n)的圆周卷积，则其每段卷积结果的前M-1个点不等于线性卷积值，需舍去。

为了清楚地看出这点，研究x(n)中的任意一段列长为N的序列xi(n)与列长为M的h(n)的圆周卷积情况


y′i(n)=xi(n)*h(n)=∑N-1m=0xi(m)h((n-m))NRN(n)
(3.4.24)

由于h(n)的列长为M，当在0≤n≤M-2范围内进行圆周位移时，h((n-m))N将在xi(m)的尾部出现有非零值，如图3.4.12(c)所示的n=1的情况。所以在这一部分0≤n≤M-2的y′i(n)值中将混入xi(m)尾部与h((n-m))N的卷积值，从而使y′i(n)不同于线性卷积结果。但是当n从M-1开始，直到N-1点，则有h((n-m))N=h(n-m)，如图3.4.12(d)，(e)所示。因此从n=M-1点后，圆周卷积值完全与线性卷积值一样，y′i(n)是正确的卷积值，因而每一段卷积运算结果的前M-1个值需去掉，如图3.4.12(f)所示。



图3.4.12用保留信号代替补零后局部的混叠现象


为了不造成输出信号的遗漏，对x(n)分段时，需使相邻两段有M-1个点的重叠(对于第一段x(n)由于没有前一段保留信号，则在其前填充M-1个零值点)，为此将xk(n)定义为


xk(n)=x［n+k(N-(M-1))-M+1］，0≤n≤N-1
0,其他n
(3.4.25)

式中已规定将每段的时间原点放在该段的起始点，而不是x(n)的原点。

这种分段方法描绘于图3.4.13(a)。每段和h(n)的圆周卷积以y′k(n)表示，如图3.4.13(b)所示。图中已标明每输出波形段开始的0≤n≤M-2部分须舍去，把相邻各输出段留下的序列衔接起来，就构成了最终的没有遗漏而又正确的输出，即


y(n)=∑∞k=0yk(n-k(N-M+1))

这里n为总的输出序列y(n)的序号。而式中


yk(n)=y′k(n),M-1≤n≤N-1
0,其他
(3.4.26)

每段输出的时间原点放在yk(n)的起始点，而不是y(n)的原点。



图3.4.13重叠保留法示意图



重叠保留法的名称是因为每个相继的输入段均由N-M+1个新点和前一段保留下来的M+1个点所组成而得名的。




视频讲解

3.4.10圆周相关定理

圆周相关也称循环相关。

若X3(k)=X*1(k)X2(k)

则


x3(n)=IDFT［X3(k)］=∑N-1i=0x*1(l)x2((l+n))NRN(n)
(3.4.27)

证明： 

由于X1(k)和X2(k)都隐含周期性，所以其乘积X3(k)同样隐含周期性。将X3(k)周期延拓，则


X～3(k)= X～*1(k)X～2(k)

即


x～3(n)=IDFS［X～3(k)］=1N∑N-1k=0X～*1(k)X～2(k)W-nkN


将


X～*1(k)=∑N-1l=0x～1(l)WlkN*

代入，则


x～3(n)=1N∑N-1k=0∑N-1l=0x～*i(l)X～2(k)W-(n+l)kN
=∑N-1i=0x～*i(l)1N∑N-1k=0X～2(k)W-(n+l)kN
=∑N-1i=0x～*1(l)x～2(n+l)
=∑N-1l=0x*i((l))Nx2((n+l))N

因0≤l≤N-1，x*1((l))N=x*1(l)故有


x3(n)= x～3(n)RN(n)=∑N-1l=0x*l(l)x2((n+l))NRN(n)
(3.4.28)


如果仅考虑x1(l)是实数，其共轭还是其本身，则得


x3(n)=∑N-1i=0x1(l)x2((n+l))NRN(n)

可以得到两个N点实序列x1(n)与x2(n)的N点圆周相关，而x2((n+l))NRN(n)正好是x2(n)的圆周移位，变量是n，且无须折叠。若定义x3(n)=∑∞m=-∞x1(m)x2(n+m)为线性相关，则两列长分别为N1及N2的线性相关，应等于将这两序列补零到N=N1+N2-1列长后的圆周相关。通过圆周相关定理可知，实信号的相关性，可借助于离散傅里叶正逆变换求得。当一个较短的序列与一个很长序列相关时，可类似圆周卷积的做法，采用分段相关的办法。

3.4.11帕塞瓦尔定理

若X(k)=DFT［x(n)］，则


∑N-1n=0x2(n)=1N∑N-1k=0|X(k)|2
(3.4.29)

证明： 

在圆周相关定理中，取x2(l)=x1(l)，且n=0，则


∑N-1i=0x1(l)x2(l)=1N∑N-1k=0X*(k)X(k)W0N

所以


∑N-1n=0x2(n)=1N∑N-1k=0|X(k)|2

帕塞瓦尔定理建立了离散函数时域能量和频域能量之间的关系，又称为能量定理。

3.4.12离散傅里叶变换的奇偶性及对称性
1. 周期性共轭对称和共轭反对称分量

第2章曾讨论过任一序列都可分解成共轭对称与共轭反对称两个分量之和。共轭对称分量为


xe(n)=12［x(n)+x*(-n)］
(3.4.30)

共轭反对称分量为


xo(n)=12［x(n)-x*(-n)］
(3.4.31)


在讨论有限长序列离散傅里叶变换对称性质时，一般不采用第2章给出的共轭对称和共轭反对称分量的定义。因为列长为N的序列X(k)其离散傅里叶变换的列长也为N，而其xe(n)及xo(n)的列长却都为(2N-1)(见图3.4.14(g)~图3.4.14(j))，难以由它们找出DFT的对称性。因为DFT的列长为N。



图3.4.14周期与非周期序列的共轭对称及共轭反对称分量以及周期性共轭对称及共轭反对称分量


对于周期为N的周期序列x～(n)，它的共轭对称分量和共轭反对称分量仍然都是周期性的，并且周期仍为N。可将x(n)分解成两个列长为N的有限长序列，其一对应于x～(n)的共轭对称分量的一个周期，记作xep(n)，另一对应于x～(n)的共轭反对称分量的一个周期，记作xop(n)。

为此，将x(n)以N为周期延拓成周期序列


x～(n)=x((n))N

则周期序列x～(n)的共轭对称和共轭反对称分量分别为


x～e(n)=12［x～(n)+x～*(-n)］(3.4.32)
x～o(n)=12［x～(n)-x～*(-n)］
(3.4.33)

取0~N-1的一个周期值，则得


xep(n)=12［x((n))N+x*((-n))N］RN(n)
=12［x(n)+x*(N-n)］
(3.4.34)
xop(n)=12［x((n))N-x*((-n))N］RN(n)
=12［x(n)-x*(N-n)］
(3.4.35)

又由式(3.4.32)及式(3.4.33)可得


x～(n)= x～e(n)+x～o(n)

因此


x(n)= x～(n)RN(n)
=［x～e(n)+x～o(n)］RN(n)
=xep(n)+xop(n)
(3.4.36)


序列xep(n)和xop(n)分别称为列长为N的序列x(n)的周期性共轭对称分量和周期性共轭反对称分量。当xep(n)和xop(n)为实序列时，分别称作周期性偶分量和周期性奇分量。这样的术语易使人误解，因为序列xep(n)和xop(n)并不是周期序列，它们只是分别表示周期序列x～e(n)和x～o(n)的一个周期，其间关系如图3.4.14(c)~图3.4.14(f)所示。

显然，xep(n)和xop(n)与式(2.3.3)定义的xe(n)和式(2.3.4)定义的xo(n)并不等价。然而不难证明，xep(n)和xop(n)与列长为(2N-1)的xe(n)和xo(n)有下述关系


xep(n)=［xe(n)+xe(n-N)］RN(n)
(3.4.37)
xop(n)=［xo(n)+xo(n-N)］RN(n)
(3.4.38)


以上讨论也完全适用于序列X(k)，并得出类似的表达式。

最后还须指出： 因为本书把有限长序列视为周期为N的时间序列中的一个周期，因此有


x(-n)=x(N-n)
(3.4.39)

同样地


X(-k)=X(N-k)
(3.4.40)

2. 奇偶序列的DFT

1) 奇序列的DFT

当序列是奇对称时，即


x(n)=-x(-n)=-x(N-n)
(3.4.41)

则其离散傅里叶变换也是奇对称的，即


X(k)=-X(-k)=-X(N-k)
(3.4.42)


证明： 


X(k)=∑N-1n=0x(n)WknN=∑N-1n=0［-x(-n)］W(-k)(-n)N=-X(-k)
∑N-1n=0［-x(N-n)］W(N-k)(N-n)N=-X(N-k)

式(3.4.42）证明中，使用了下述关系


W(N-k)(N-n)N=WN2N·W-kNN·W-nNN·WknN=WknN


2) 偶序列的DFT

当序列是偶对称时，即


x(n)=x(-n)=x(N-n)
(3.4.43)

则其离散傅里叶变换也是偶对称的，即


X(k)=X(-k)=X(N-k)
(3.4.44)

可用与奇序列的DFT相似的方法证明。

3. 共轭复序列的DFT

若x*(n)为x(n)的共轭复序列，则


DFT［x*(n)］=X*(N-k)
(3.4.45)

证明： 


DFT［x*(n)］=∑N-1n=0x*(n)WnkN=∑N-1n=0x(n)W-nkN*，0≤k≤N-1

由于


WnNN=e-j2πNnN=e-j2πn=1


故


DFT［x*(n)］=∑N-1n=0x(n)W(N-k)nN*
=X*((N-k))N
=X*(N-k)


请注意，k=0时不能使用等式X*((N-k))N=X*(N-k)，因为这时X*((N-0))N=X*(N-0)=X*(N)，而X(k)只有0≤k≤N-1范围内的N个值，所以已超出取值区间。式(3.4.45)的严格公式应该是


DFT［x*(n)］=X*((N-k))N
(3.4.46)

这样，当k=0时有X*((N-0))N=X*(0)的正确结果。因为X(k)可认为是分布在N等分的圆周上，它的末点即它的始点，也即X(N)=X(0)。因此仍采用式(3.4.45)的习惯形式。在下文有关对称性的讨论中，凡是X(N)都认为是X((N))N=X(0)。

4. 复序列的DFT

若有限长序列x(n)是一个复序列，设xr(n)及jxi(n)分别表示x(n)的实部与虚部，即



x(n)=xr(n)+jxi(n)
(3.4.47)

而


xr(n)=12［x(n)+x*(n)］
jxi(n)=12［x(n)-x*(n)］
(3.4.48)

以Xep(k)及Xop(k)分别表示实部及虚部序列的DFT，则


Xep(k)=DFT［xr(n)］
=12DFT［x(n)+x*(n)］
=12［X(k)+X*(N-k)］(3.4.49)
Xop(k)=DFT［jxi(n)］
=12DFT［x(n)-x*(n)］
=12［X(k)-X*(N-k)］(3.4.50)


根据线性特性


X(k)=Xep(k)+Xop(k)
(3.4.51)

注意，这里Xep(k)与Xop(k)均为复数，所以式(3.4.51)右端总的实部和总的虚部还是不能直接表示出来。

现在分析一下Xep(k)与Xop(k)的一些对称特性。由式(3.4.49)得


X*ep(N-k)=12［X(N-k)+X*(N-N+k)］*
=12［X*(N-k)+X(k)］
(3.4.52)

与式(3.4.49)比较，可知


Xep(k)=X*ep(N-k)
(3.4.53)

因此Xep(k)称为X(k)的周期性共轭对称分量，其含义在上文已经说明。由此可认为Xep(k)是分布在N等分圆周上，则以k=0为原点，Xep(k)在左半圆上的序列与右半圆上的序列是共轭对称的。图3.4.15直观地表示了Xep(k)为实数时的示意图。若Xep(k)为复数，则其共轭对称的含义是模相等、幅角相反，即

|Xep(k)|=|Xep(N-k)|
arg［Xep(k)］=-arg［Xep(N-k)］
(3.4.54)

或说实部相等、虚部相反。



图3.4.15周期性共轭对称分量Xep(k)


根据式(3.4.50)可以推得


Xop(k)=-X*op(N-k)
(3.4.55)

Xop(k)称为X(k)的周期性共轭反对称分量，它表示Xop(k)分布在圆周上时，以k=0为中心，Xop(k)在左半圆上的序列与Xop(k)在右半圆上的序列是共轭反对称的。图3.4.16是Xop(k)为实数时的示意图。当Xop(k)为复数时，则应按式(3.4.55)的意义去理解共轭反对称，即实部相反、虚部相等。


|Xep(k)|=|Xep(N-k)|
arg［Xep(k)］=-arg［Xep(N-k)］
(3.4.56)



图3.4.16周期性共轭反对称分量Xop(k)




通过上面的分析可知，对于一个复时间序列的DFT来说，序列的实部对应于X(k)的周期性共轭对称分量Xep(k)； 序列的虚部对应于X(k)的周期性共轭反对称分量Xop(k)。这种结果是有一定意义的。在此可以把两个实数序列x1(n)和x2(n)组合为单一的复数函数x(n)，如式(3.4.47)所示。当算出复数表示的X(k)后，即可利用式(3.4.49)和式(3.4.50)将X(k)分成两个独立的分量Xep(k)和Xop(k)，它们分别对应于x1(n)和x2(n)的离散傅里叶变换。因此利用这种方法，在一次计算中可得出两个独立信号的变换。

根据x1(n)与x2(n)的对称特性，同样可以找到X(k)的实部、虚部与x(n)的周期性共轭对称分量Xep(k)、周期性共轭反对称分量Xop(k)的关系。请读者自行证明： Xep(k)的DFT是X(k)的实部Re［X(k)］，Xop(k)的DFT是X(k)的虚部jIm［X(k)］。

5. 虚实序列的DFT


1) 虚序列的DFT

当序列是纯虚序列，即x(n)=jxi(n)时，则由式(3.4.50)知其X(k)只有周期性共轭反对称分量，即X(k)=xop(k)，由式(3.4.54)可见，虚序列的离散傅里叶变换X(k)的实部是奇对称的，虚部是偶对称的。

2) 实序列的DFT

若x(n)是纯实数序列，即x(n)=xr(n)时，由式(3.4.49)知其X(k)只有周期性共轭对称分量，即X(k)=Xep(k)，由式(3.4.54)知，其模是偶函数，而相角是奇函数，或者说，其实部是偶对称的，而虚部是奇对称的。

上述两种情况不论哪一种都只要知道一半数目的X(k)，利用对称特性就可以得到另一半数目的X(k)。在DFT中利用这个特点，可以提高运算效率。

联合考虑序列的奇偶特性及虚实特性，还可证明如下特性： 

(1) 当序列是实的偶对称序列，则其离散傅里叶变换是实的偶对称的。

(2) 当序列是实的奇对称的序列，则其离散傅里叶变换是虚的奇对称的。

(3) 当序列是虚的偶对称的序列，则其离散傅里叶变换是虚的偶对称的。

(4) 当序列是虚的奇对称序列时，则其离散傅里叶变换是实的奇对称的。

在结束有关奇偶性及对称性的讨论时，为便于比较，将DFT的奇偶虚实特性列于表3.4.1中。


表3.4.1DFT的奇偶虚实特性



x(n)X(k)x(n)X(k)


偶序列偶序列实偶实偶
奇序列奇序列实奇虚奇
实实部为偶，虚部为奇虚偶虚偶
虚实部为奇，虚部为偶虚奇实奇



根据x(n)与X(k)变换关系的对称性，只要交换x(n)和X(k)列的标题，表中各项仍然是正确的。

应该指出，为了更有效地运用计算程序和更好地理解所得结果，对上述特性应深入地理解。

3.4.13可将离散傅里叶变换看作一组滤波器

序列x(n)的离散傅里叶变换


X(k)=DFT［x(n)］=∑N-1n=0x(n)WnkN
=x(0)+x(1)WkN+x(2)W2kN+x(3)W3kN+…+
x(N-1)W(N-1)kN
(3.4.57)

X(k)为序列x(n)中k频率分量的大小，可看作如下的卷积形式


X(k)=∑N-1n=0x(n)hk(N-1-n)
(3.4.58)


令系统的单位采样响应为


hk(n)=W(N-1-n)kN，n=0,1,2,…，N-1
0，n<0和n≥N
(3.4.59)

输入序列为


0,0,0,x(0),x(1),x(2),…,x(N-1),0,0,0,…

此输入序列加入单位采样响应为hk(n)的滤波器后，在(N-1)时刻的输出为


y(N-1)=∑N-1n=0x(n)hk(N-1-n)
=x(0)hk(N-1)+x(1)hk(N-2)+
x(2)hk(N-3)+…+x(N-1)hk(0)
=x(0)W(N-1-N+1)kN+x(1)W(N-1-N+2)kN+
x(2)W(N-1-N+3)kN+…+x(N-1)W(N-1)kN
=x(0)+x(1)WkN+x(2)W2kN+x(3)W3kN+…+
x(N-1)W(N-1)kN
(3.4.60)


比较式(3.4.57)和式(3.4.60)可见，二者的结果是一样的。因此，求序列x(n)的离散傅里叶变换X(k)的过程相当于以序列x(n)为输入，加到单位采样响应为hk(n)的滤波器，其在(N-1)时刻的输出就是X(k)。由式(3.4.59)可见，(N-1)与k值有关，一定的k值对应一定的单位采样响应，取不同的k值可得不同的单位采样响应，分别对应不同滤波特性的数字滤波器。把k值取在0≤k≤N-1区间，就可得N个不同滤波特性的数字滤波器。对x(n)实施N个滤波器的滤波运算就可分别得到N个不同k值的X(k)。所以离散傅里叶变换可看作一组滤波器。

下面分析这组滤波器的频率特性。为此先对hk(n)作z变换可得


Hk(z)=∑N-1n=0hk(n)z-n
=∑N-1n=0W(N-1-n)kNz-n
=W(N-1)kN∑N-1n=0(W-kNz-1)n
=W(N-1)kN1-z-NW-NkN1-z-1W-kN
(3.4.61)

由于WN=e-j2πN，同时令z=ejω，代入式(3.4.61)，得出滤波器的振幅频率特性为


|Hk(ejω)|=1-(ejω)-Nej2πNkN1-(ejω)-1ej2πNk1=
ejN2ω-2πNk-e-jN2ω-2πNkej12ω-2πNk-e-j12ω-2πNk

根据公式sinx=ejx-e-jx2j可得


|Hk(ejω)|=sinN2ω-2πNksin12ω-2πNk
(3.4.62)

若令x=12ω-2πNk，则这种振幅频率特性具有sinNxsinx的形式。在k=0时，


H0=sinNω2sinω2

也是一种sinNxsinx的频率响应特性，当ω=0时，H0=N，可作出其振幅特性如图3.4.17所示，第一副瓣电平约为-13.3dB。



图3.4.17k=0时，H0的归一化振幅特性


若k≠0，而分别取1,2,3,…，N-1，则由式（3.4.62）可得出一组滤波器的振幅特性如图3.4.18所示。图中N=8。各Hk的响应形式和H0一样，最大幅值也为N，差别在于主瓣移到了ω=2πNk上。这些滤波器中心频率的间隔为ω=2πN。



图3.4.18N=8时，不同k值的Hk的振幅特性


因此，离散傅里叶变换相当于用中心频率为ω=2πNk(k=0,1,2,…，N-1)，频率响应形式为sinNxsinx的N个滤波器对输入序列进行滤波。因为ω=ΩT=Ωfs，在ω=2πN时，有Ω=2πNfs=2πNT。故当ω=k2πN时这些滤波器在(N-1)时刻的输出即是频率分量为kΩ的傅里叶系数X(k)。当输入序列x(n)中所含频率分量不是Ω的整倍数，而是在(k-1)Ω和kΩ之间时，输入序列经傅里叶变换后，在各个滤波器的输出端都可能有它的输出响应，而在第k个及第k-1个滤波器的输出响应中，将有一个或两个都具有较大值。所以利用离散傅里叶变换进行谱分析时，分辨能力决定于频谱谱线间隔Ω=2πNT=2πf1N，并称f1N为频率分辨率或频率分辨单元。

离散傅里叶变换相当于一组滤波器，这一性质在脉冲多普勒雷达信号处理中得到了广泛的应用。

3.4.14DFT与z变换

列长为N的有限长序列x(n)的z变换为X(z)，而其离散傅里叶变换为X(k)，X(z)和X(k)二者均表示了同一有限长序列x(n)的变换，它们之间的关系是： 对z变换采样可得 DFT，而DFT的综合就是z变换。



图3.4.19DFT与z变换

由前文的讨论可知，有限长序列z变换的收敛域是全部z平面，自然也包含了z平面的单位圆。如果在单位圆上等间隔取N个点，如图3.4.19(a)所示，则在此N个点处的z变换值为


X(zk)=∑N-1n=0x(n)z-nz=zk=ej2πNk=W-kN
=∑N-1n=0x(n)WnkN
=DFT［x(n)］
=X(k)(3.4.63)

zk=W-kN=ej2πNk是平面单位圆上幅角为 ω=2πNk的点，即z平面上单位圆N等分后的第k点。所以x(n)DFT的N个系数X(k)也即x(n)的z变换X(z)单位圆上等距离的采样值，如图3.4.19(b)所示。z变换在单位圆上的值就是X(ejω)。所以也可以说，X(k)是序列傅里叶变换X(ejω)在相应点ω=2πNk=kωN上的采样值，其采样间隔为ω=2πNk，即


X(k)=X(ejkωN)
ωN=2πN

由此可见，z变换的采样可得DFT。




视频讲解

3.5频域采样◆

由上文讨论可知，采用DFT后实现了频域采样。对于任意一个频率特性能否用频率采样的办法去逼近呢？为此，首先应弄清它的限制，再研究经过频率采样后会有什么误差？如何消除误差？采样后所获得的频率特性怎样？

3.5.1对X(z)采样时采样点数的限制

对任一绝对可和的非周期序列x(n)的z变换X(z)在单位圆上，即对X(ejω)进行等距采样，则由式(3.4.63)得


X(k)=X(z)z=W-kN=∑∞n=-∞x(n)WnkN
(3.5.1)

实现频域采样以后，信息有没有损失？能不能用序列的频谱X(ejω)的采样值X(k)恢复出原序列x(n)？设频率采样后所得的X(k)恢复出的有限长序列为x′(n)，则


x′(n)=IDFT［X(k)］

为了易于看清上述问题，本书先从周期序列x～′(n)开始研究


x～′(n)=IDFS［X～(k)］=1N∑N-1k=0X～(k)W-knN=1N∑N-1k=0X(k)W-knN

为探索频率采样后所恢复出的序列x′(n)和原序列x(n)之间的关系，将式(3.5.1)代入


x～′(n)==1N∑N-1k=0∑∞m=-∞x(m)WmkNW-knN=∑∞m=-∞x(m)1N∑N-1k=0W(m-n)kN

由于


1N∑N-1k=0W(m-n)kN=
1,m=n+rN,r为任意整数
0,其他m

所以


x～′(n)=∑∞r=-∞x(n+rN)
(3.5.2)


式(3.5.2)说明x～′(n)是原序列x(n)以N为周期的周期延拓序列。在第2章看到了时域的采样造成频域的周期延拓，现又证明了频域上的采样，同样也造成时域的周期延拓。这正是傅里叶变换中时域频域对称关系的反映。

如果x(n)是列长为M的有限长序列，则当频域采样点数N<M时，即频域的采样间隔不够密，x(n)的周期延拓就会出现某些序列交叠在一起，产生混叠现象。这样就不可能从x(n)中提取一个周期不失真地恢复出原序列x(n)来。因此对于列长为M的有限长序列x(n)，频率采样不失真的条件是N≥M,这时有


x′(n)= x～′(n)RN(n)=∑∞r=-∞x(n+rN)RN(n)=x(n)(3.5.3)


当x(n)为无限长序列时，无论N取什么值，x′(n)都将不可能完全消除混叠误差, 只能随着采样点N的增加而逐渐接近x(n)。

3.5.2X(z)的内插公式

由上文分析可知，列长为N的有限长序列x(n)，可从单位圆上X(z)的N个采样值X(k)，即x(n)的DFT恢复，因而这N个X(k)也应该能完全表达整个X(z)函数及频响X(ejω)，也可以说DFT的综合就是z变换。这进一步揭示了DFT与z变换的关系。


x(n)=1N∑N-1k=0X(k)W-nkN
X(z)=∑N-1n=0x(n)z-n

将x(n)的表达式代入得


X(z)=∑N-1n=01N∑N-1k=0X(k)W-nkNz-n
=1N∑N-1k=0X(k)∑N-1n=0(W-kNz-1)n
=1N∑N-1k=0X(k)1-W-kNNz-N1-W-kNz-1

因


W-kNN=ej2πNkN=1

即


X(z)=1-z-NN∑N-1k=0X(k)1-W-kNz-1
(3.5.4)

式(3.5.4)就是X(z)的内插公式。在已知X(k)时，可根据内插公式求得任意z点的X(z)值，因此X(z)的N个采样点的X(k)值，包含了z变换的全部信息。式(3.5.4)可表示为式(3.5.5)的形式内插函数

X(z)=∑N-1k=0X(k)k(z)
(3.5.5)

其中k(z)为

k(z)=1N1-z-N1-W-kNz-1
(3.5.6)

将z=ejω代入，可以得到序列x(n)的频响


X(ejω)=∑N-1k=0X(k)k(ejω)
(3.5.7)

k(ejω)=1N1-e-jNω1-e-jω-k2πN
=1NsinωN2sinω-k2πN2·e-jNω2-ω2+kπN
(3.5.8)

可将k(ejω)表示为更明了的形式


k(ejω)=ω-k2πN
(ω)=1NsinωN2sinω2e-jωN-12
(3.5.9)

因此


X(ejω)=∑N-1k=0X(k)ω-k2πN
(3.5.10)

内插函数(ω)的振幅与相位特性如图3.5.1所示。当其变量ω=0(本采样点)时(ω)=1，在其余采样点(ω=i2πN,i=1,2,…，N-1)上，(ω)=0。因而可知，有以下关系


ω-k2πN=
1，ω=k2πN
0，ω=i2πN,i≠k



图3.5.1内插函数振幅特性与相位特性


即函数ω-k2πN在本采样点ω=k2πN上为1，而在其他采样点ω=i2πN,i≠k上为0。整个X(ejω)正是N个ω-k2πN函数乘以加权X(k)之和。明显可见，各个采样点处的X(ejω)就等于各X(k)的值，因为其余采样点的内插函数在这里都为零值。而采样点之间的X(ejω)值则由各采样值乘以相应的内插函数延伸叠加形成。

内插函数的另一重要特点是具有如图3.5.1所示的线性相移特性。

综上所述，对于有限长序列的X(z)及X(ejω)可有时域序列x(n)与频域序列X(k)表示的两套表达式


X(z)=∑N-1n=0x(n)z-n
X(z)=∑N-1k=0X(k)k(z)

及


X(ejω)=∑N-1n=0x(n)e-jωn
X(ejω)=∑N-1k=0X(k)kω-k2πN


时域采样定理说明，一个频带有限的信号，可对其进行时域采样而不丢失任何信息，因而,可用数字技术加工时域信号。频域采样理论说明，时间有限的信号(有限长序列)，也可对其频域采样而不丢失任何信息。DFT的理论使信号不仅在时域，而且在频域也离散化了，因而开辟了在频域采用数字技术处理的领域，而快速傅里叶变换(FFT)是计算DFT的快速有效算法。




视频讲解

3.6用DFT对连续时间信号逼近的问题◆

DFT的数学性质是确切的，但在许多信号处理的应用中，却很少作为一个最终目的被采用。对DFT感兴趣主要是因为它是连续傅里叶变换的一个近似。为了利用DFT对连续时间信号xa(t)进行傅里叶分析，需先对xa(t)进行采样，得到x(n)，再对x(n)进行DFT得X(k)。X(k)是x(n)的傅里叶变换X(ejω)在频率区间［0,2π］上的N点等间隔采样。这里x(n)和X(k)均为有限长序列。但是，由傅里叶变换理论可知，若信号的持续时间为有限长，则其频谱无限宽； 若信号的频谱为有限宽，则其持续时间无限长。严格地讲持续时间有限的带限信号是不存在的。为能满足DFT的变换条件，实际上对频谱很宽的信号，为防止时域采样后产生频谱混叠失真，可用前置滤波器滤除幅度较小的高频分量，使连续时间信号的带宽小于折叠频率。对于持续时间很长的信号，采样点数太多以致无法存贮和计算，只好截取为有限列长进行DFT。从工程实际角度看，滤除幅度很小的高频分量和截去幅度很小部分的时间信号是允许的。由上述可见，用DFT对连续时间信号进行傅里叶分析必然是近似的，近似的准确程度严格地说是被分析波形的一个函数。特别要指出，两个变换之间的差异是因DFT需要对连续时间信号采样和截断为有限列长而产生的。这里主要有两个问题，即两种变换间相对数值的确定； 以及在计算的变换与所需的变换之间造成误差的3种可能现象： 混叠现象、栅栏效应和频谱泄漏。下面进行分别讨论。

3.6.1计算的变换与所需变换间相对数值的确定

假设连续时间非周期性信号的变换定义为式(3.1.1)和式(3.1.2)，而周期信号的傅里叶级数的定义则为式(3.1.3)和式(3.1.4)。DFT的计算可以按式(3.1.1)所定义的非周期信号的傅里叶变换，或式(3.1.3)所定义的周期信号的傅里叶系数的近似来进行。

如果采用DFT的基本定义式(3.3.1)去计算一个非周期性信号的傅里叶变换，则频谱的正常幅度电平等于用DFT计算所得的频谱分量乘以T。

如果已利用真实积分变换定义式(3.1.1)，或者如上所述，利用DFT乘以T来计算频谱函数，则可用包括因子1/N在内的式(3.3.2)来近似计算式(3.1.2)所表示的时间函数。在此情况下，最后将各时间分量乘以NF=f，即得正常幅度电平。所以从时间到频率，再从频率到时间的整个过程总共乘了T·NF=1也即幅度电平未受影响。

傅里叶级数系数的表示式(3.1.3)实际上是求被积式的平均值。如果采用式(3.3.1)的DFT公式来近似确定式(3.1.3)所表示的周期性函数的傅里叶级数系数时，频谱的正常幅度电平等于DFT所求出的频谱分量乘以T，再除以tp(tp为周期性时间函数的有效周期)。由于(T/tp)=(1/N)，因此，如果利用DFT的定义式(3.3.1)和式(3.3.2)来计算正规的傅里叶级数，则(1/N)因子最好放在DFT的正变换式中，而不是逆变换式中。

在许多实际问题中，各变换分量的实际电平往往无关紧要。确定信号特性的是频谱中各分量之间的相对电平，相对电平关系保持不变就可以了。

3.6.2计算的变换与所需变换间的误差
1. 混叠现象

混叠现象在上文已详细讨论过。避免混叠现象的唯一办法是保证采样频率足够高。这就要求在确定采样频率之前，对频谱的性质要有所了解。下面着重讨论应用DFT时，为避免混叠所必须考虑的一些重要参数关系。

假设离散时间信号是从连续时间函数采样得出的，或者为了分析的目的选取相应的连续时间信号作为参考。并假设所处理的信号是基带信号，在采样前已利用低通模拟滤波器进行前置滤波，以避免高于折叠频率的分量出现。这样为避免混叠现象，要求采样频率


fs≥2fh
(3.6.1)


fs为信号的最高频率，而采样周期T必须满足


T≤12fs
(3.6.2)


设F表示频率分量间的增量，它就是前文提到的频率分辨率F=fsN，tp为最小记录长度，也就是前文提到的周期性函数的有效周期。tp应按照所需的频率分辨率进行选择


tp=1F
(3.6.3)

式(3.6.2)和式(3.6.3)说明在高频容量与频率分辨率间存在着矛盾。增加高频容量，T就必然减小，在采样点数N给定的情况下，记录长度必然缩短，从而降低了频率的分辨率。相反，要提高记录长度，分辨率就必须要增加,在采样点数N给定时，必然导致T的增加，因而减少了高频的容量。

在高频容量fh与频率分辨率F两个参数中，保持其中一个不变而增加另一个的唯一办法，就是增加在一记录长度内的点数N,如果fh和F都已给定，则N必须满足


N≥2fhF
(3.6.4)


这是未采用任何特殊数据处理(例如加窗处理)情况下，为实现基本的DFT算法所必须满足的最低条件。

2. 栅栏效应

栅栏效应是由于用DFT计算频谱只限制为基频的整数倍而不可能将频谱视为一个连续函数而产生的。就意义而言，栅栏效应表现为当用DFT计算整个频谱时，就好像通过一个“栅栏”来观看一个图景一样，只能在离散点的地方看到真实图景。如果不附加任何特殊处理，在两个离散的变换线之间若有一个特别大的频谱分量时，将无法检测出来。

减少栅栏效应的方法就是在原记录末端添加一些零值点来改变时间周期内的点数，并保持记录不变。从而在保持原有频谱连续形式不变的情况下，变更了谱线的位置。这样，原来看不到的频谱分量就能移动到可见的位置上。

必须注意，当在记录信号末端添加零点时，所用窗函数的宽度不能由于增加了零点而按较长的长度选择，而必须按照数据记录的实际长度来选择窗函数。关于窗函数的有关知识将在3.7节介绍。

3. 频谱泄漏

实际工作往往需要把信号的观察时间限制在一定的时间间隔之内。设有一个延伸到无限远处的离散时间信号x1(n)，其频谱为X1(ejω)(仅表示出周期频谱中周期的一部分)如图3.6.1所示。由于无法等待足够长的时间取用无限个数据，因此需要选择一段时间信号进行分析。



图3.6.1信号截断时产生的频谱泄漏现象


取用有限个数据，即将信号截断的过程，就等于将信号乘以窗函数。如果窗函数是一个矩形窗函数，数据项数突然被截断，而窗内各项数据并不改变，如图3.6.1中x2(n)所示。时域的截断，在频域相当于所研究的波形的频谱X1(ejω)与矩形窗函数频谱周期卷积过程。这一卷积造成的失真频谱见图3.6.1中的X2(ejω)所示。

频谱分量从其正常频谱扩展开来，称为“泄漏”。如图3.6.2(a)所示的x′1(t)的频率f1
是fs/N的整倍数，即


图3.6.2频率为fs/N的整倍数及为fs/N的非整倍数信号




f1=kfsN=kNT


这说明在长度NT内信号有k整个周期。这时由x′1(t)构成的以NT为周期的周期信号是连续的。当x′2(t)的频率不是fs/N的整倍数时，在NT的处理长度内，就不是恰好为信号周期的整数倍，由x′2(t)以NT为周期进行周期延拓所得到的周期信号就出现了不连续点，如图3.6.2(b)所示，造成了频谱分量从其正常频谱扩展开来，在DFT等效滤波器组的各个副瓣将有大小不同的数值输出，再一次看到了频谱泄漏现象。应该指出，泄漏是不能与混叠完全分开的，因为泄漏将导致频谱的扩展，从而使频谱的最高频率超过折叠频率，造成混叠。因为无法取用无限个数据，所以在进行离散傅里叶变换时，时域中的截断是必须的，因此泄漏效应是离散傅里叶变换所固有的，必须设法进行抑制。




视频讲解

3.7窗函数和加权◆

由上文分析可知，要抑制“泄漏”可以通过窗函数加权抑制DFT的等效滤波器的振幅特性的副瓣，或用窗函数加权使有限长度的输入信号周期延拓后在边界上尽量减少不连续程度的方法实现。在FIR数字滤波器设计中，为获得有限长单位采样响应，需要用窗函数截断无限长单位采样响应序列； 在功率谱估计中也要用到窗函加权问题。由此可见窗函数加权技术在数字信号处理中的重要地位。

在用窗函数加权时，有一个问题必须注意。一般文献上给出的窗函数w(n)都是偶对称的时间序列，即


w(n)n=-N2,…,-1,0,1,…，+N2


这里，N是一个偶数，而窗函数w(n)共有N+1个采样值。但实际中常需要单边表示的窗函数，例如在N为偶数点的离散傅里叶变换时，处理区间为n=0~N-1。因此，必须把偶对称表示的窗函数向右平移N2点，让左端点与n=0重合。如除去右端的一个采样值(一般为零值)，则窗函数w(n)便在n=0,1,2，…，N-1点上定义，构成了适用于离散傅里叶变换的单边窗函数序列。也可将被加权的序列向左平移N2点，作用是一样的。位移N2点只影响相位特性，并不影响振幅特性。本节主要研究加权的作用及窗函数的性能。

3.7.1加权

现对离散傅里叶变换的等效滤波器的单位采样响应进行加权，即


hk(n)=w(n)W(N-1-n)kN

式中,w(n)称为窗函数，用不同的窗函数对单位采样响应进行加权可以得到等效滤波器不同的振幅频率特性。

若w(n)=1,n=0,…，N-1，称为矩形窗。这时加权后的单位采样响应与等效滤波器原有的单位采样响应一样


hk(n)=W(N-1-n)kN，n=0,1,2,…,N-1
0，n<0和n≥N
(3.7.1)


其频率响应为


Hk(ejω)=sinN2ω-2πNksin12ω-2πNke-jN-12ω+2πNk

其中，k为零所对应的频率特性，由下式表示


H(ejω)=sinNω2sinω2e-jN-12ω


上文已指出，这种频率特性的第一副瓣峰值比ω=0处的主瓣峰值低13.3dB,约为主瓣幅度的1/4~1/5。采用合适的窗函数进行加权，就能使对应的频率特性副瓣压降下来。例如，采用余弦加权


w(n)=cosπnN，n=-N2,…,-1,0,1,…,+N2
(3.7.2)


为适用于离散傅里叶变换加权，需将偶对称的余弦窗移位N2后采用正弦窗函数的形式，即


w(n)=cosπnN-π2=sinπnN
(3.7.3)


这时对应的加权单位采样响应为


hk(n)=sinπnNW(N-1-n)kNn=0,1,2,…,N-1
(3.7.4)


其频率特性可通过z变换求得


Hk(ejω)=Hk(z)z=ejω
=∑N-1n=0sinπnNW(N-1-n)kNz-n
=W(N-1)k∑N-1n=0sinπnNW-knNz-n


暂不考虑求和号外的相移因子W(N-1)kN，可求出


Hk(ejω)=∑N-1n=0sinπnNej2πNknz-n
=12j∑N-1n=0ejπnN-e-jπnNej2πNkn·e-jωn
=12j∑N-1n=0e-jnω-k+122πN
-12j∑N-1n=0e-jnω-k-122πN
(3.7.5)


把式(3.7.5)的最后两求和式总和起来，可得


Hk(ejω)=12jsinN2ω-k-122πNsin12ω-k+122πNe-jN-12ω-k+122πN-
sinN2ω-k-122πNsin12ω-k-122πNe-jN-12ω-k-122πN
=12jsinN2ω-k+122πNsin12ω-k+122πN-
sinN2ω-k-122πNsin12ω-k-122πNe-j(N-1)πNe-jN-12ω-k+122πN


由于e-jN-1Nπ=-ejπN，故有


Hk(ejω)=12sinN2ω-k+122πNsin12ω-k+122πNe-jπ2N+12sinN2ω-k-122πNsin12ω-k-122πNejπ2N×
e-jN-12ω-k+122πN+π2-π2N
(3.7.6)


在式(3.7.6)的相位中还应加上文提到过的相位因子才是Hk(ejω)的完整表达式。

由式(3.7.6)可知，余弦加权后的频率特性Hk(ejω)是由两相邻的中心频率间隔为2π/N的矩形窗加权的频率特性矢量合成的。由于N>1，故括号内的两个振幅频率特性间的相位差π/N很小，因此矢量合成接近于代数相加，如图3.7.1所示。由图3.7.1可见，这时两个相邻的振幅频率特性的副瓣接近相减。因此通过选择合适的窗函数确实可以压低等效滤波器频率特性副瓣，达到抑制“频谱泄漏”的目的。由图3.7.1还可以看出，加权后还会引起频率特性的主瓣加宽，峰值响应降低的副作用。



图3.7.1对离散傅里叶变换等效滤波器单位采样响应余弦加权后的频率响应


也可以在离散傅里叶变换式中对输入序列x(n)直接进行加权。这时，合适窗函数的加权作用是使被加权序列在边缘(n=0和n=N-1附近)比矩形窗函数圆滑而减小了陡峭的边缘所引起的副瓣分量。

3.7.2常用的窗函数

下面讨论几种常用的窗函数及其主要性能指标。各窗函数的幅度响应以分贝形式表示，定义为


WdB(ω)=20log10|W(ejω)|W(ej0) 
(3.7.7)


W(ejω)为w(n)的傅里叶变换,W(ej0)为该变换的直流值。

1. 矩形窗(rectangular window)


w(n)=1,n=-N2,…,-1,0,1，…,N2 
(3.7.8)

式中,N为偶数。对序列的截断，实际上就是加矩形窗，其对应的频谱为


W(ejω)=sin(N+1)2ωsinω2
(3.7.9)

上述窗函数序列长为N+1个点。例如取偶数点，即令


w(n)=1,n=-N2,…,-1,0,1，…，N2
(3.7.10)

对应的谱函数为


W(ejω)=WR(ejω)=∑N2-1n=-N2e-jωm=ejωN2∑N-1n=0e-jωn=sinNω2sinω2
ejω2
(3.7.11)

以后本书就取式(3.7.11)所表示的函数为矩形窗频谱。

如果用幅度函数W(ω)与相位函数e-jωα来表示窗函数的频谱W(ejω)，即


W(ejω)=W(ω)e-jωα 

则上述窗函数的频谱幅度函数W(ω)具有sinNxsinx的形式，其频谱的主瓣宽度，以两个零交点之间的间隔计算，为2×2πN,第一副瓣电平比主瓣峰值低13dB左右。图3.7.2给出了矩形窗函数及其对应的频谱函数。本节所述的所有窗函数及其频谱图，都是根据相应的窗函数及其频谱函数表达式在计算机上绘成的。

单边表示的矩形窗函数为


w(n)=1,n=0,1,2,…,N-1 
(3.7.12)


它所对应的频谱函数为


WR(ejω)=sinNω2sinω2e-jN-12ω 
(3.7.13)

可以看出与式(3.7.11)相比，有一个序列移位N2点所出现的相移因子e-jN2ω。

【例3.7】绘出矩形窗及其幅度响应。

import numpy as np

from scipy import signal, fft

import matplotlib.pyplot as plt

from matplotlib.ticker import MaxNLocator



# 计算矩形窗的wn值和幅度响应

N = 51# 矩形窗长度

wn = signal.windows.boxcar(N)# 矩形窗的wn值

N0 = 2048

N1 = int(N0 / 2)

Ha = np.abs(fft.fft(wn, N0)) + 1e-10

Ha = Ha / np.max(Ha)

Ar = 20 * np.log10(Ha)

freq = np.linspace(0, 1, N1)



# 绘制矩形窗

fig, ax = plt.subplots()

ax.stem(wn, basefmt="")

ax.set_title('矩形窗')

ax.set_xlabel('n')

ax.set_ylim(［0, 1.5］)

plt.rcParams［'font.sans-serif'］ = ［'SimHei'］# 用来正常显示中文标签

fig.savefig('./win_rec1.png', dpi=500)



# 绘制矩形窗的幅度响应

fig, ax = plt.subplots()

ax.plot(freq, Ar［:N1］)

ax.grid()

ax.set_title('矩形窗的幅度响应')

ax.set_xlabel('k')

ax.set_xlabel(r'$ ＼omega / ＼pi $')

ax.set_ylabel(r'$ 20log_{10}| H(＼omega) | $')

ax.set_xlim(［0, 1］)

ax.set_ylim(［-100, 1］)

ax.xaxis.set_major_locator(MaxNLocator(11))

ax.yaxis.set_major_locator(MaxNLocator(11))

plt.rcParams［'axes.unicode_minus'］ = False# 用来显示负号



矩形窗的ω(n)和幅度响应绘制如图3.7.2所示。



图3.7.2矩形窗及其频谱幅度宽度


2. 三角形窗(triangular window)

三角形窗又称巴特利特(Bartlett)窗。偶对称表达的三角形窗定义为


w(n)=1.0-|n|N2,n=-N2,…,-1,0,1，…，N2
(3.7.14)

单边表示的三角形窗定义为


w(n)=
nN2,n=0,1,2,…,N2
w(N-n),n=N2,N2+1,…,N-1
(3.7.15)

它所对应的频谱函数为


W(ejω)=N2sinN4ωsinω22e-jN2-1ω
(3.7.16)

其零交点之间的主瓣宽度是矩形窗的两倍，第一个副瓣电平比主瓣峰值低26dB左右。三角形窗是最简单的，频谱函数W(ejω)为非负的一种窗函数。三角形窗序列及其频谱如图3.7.3所示。

【例3.8】绘出三角形窗。

import numpy as np

from scipy import signal, fft

import matplotlib.pyplot as plt

from matplotlib.ticker import MaxNLocator



# 计算三角形窗的wn值和幅度响应

N = 51# 三角形窗长度

wn = signal.windows.bartlett(N)# 三角形窗的wn值

N0 = 2048

N1 = int(N0 / 2)

Ha = np.abs(fft.fft(wn, N0)) + 1e-10

Ha = Ha / np.max(Ha)

Ar = 20 * np.log10(Ha)

freq = np.linspace(0, 1, N1)



# 绘制三角形窗

fig, ax = plt.subplots()

ax.stem(wn, basefmt="")

ax.set_title('三角形窗')

ax.set_xlabel('n')

ax.set_ylim(［0, 1.5］)

plt.rcParams［'font.sans-serif'］ = ［'SimHei'］# 用来正常显示中文标签

fig.savefig('./win_tri1.png', dpi=500)



# 绘制三角形窗的幅度响应

fig, ax = plt.subplots()

ax.plot(freq, Ar［:N1］)

ax.grid()

ax.set_title('三角形窗的幅度响应')

ax.set_xlabel('k')

ax.set_xlabel(r'$ ＼omega / ＼pi $')

ax.set_ylabel(r'$ 20log_{10}| H(＼omega) | $')

ax.set_xlim(［0, 1］)

ax.set_ylim(［-100, 1］)

ax.xaxis.set_major_locator(MaxNLocator(11))

ax.yaxis.set_major_locator(MaxNLocator(11))

plt.rcParams［'axes.unicode_minus'］ = False# 用来显示负号



三角形(巴特利特)窗的w(n)和幅度响应绘制如图3.7.3所示。



图3.7.3三角形窗及其频谱幅度函数


3. 汉宁窗(Hanning window)

汉宁窗也称余弦平方窗或升余弦窗，其偶对称表示式为


w(n)=cos2nNπ
=121+cos2nNπ
=12+12cos2nNπ,n=-N2,…,-1,0,1,…，+N2
(3.7.17)

单边表示为


w(n)=sin2nNπ=121-cos2nNπn=0,1,2,…,N-1
(3.7.18)

由式(3.7.17)，据欧拉公式可得


w(n)=12+1212ej2nNπ+12e-j2nNπ
=12+14ej2nNπ+14e-j2nNπ
(3.7.19)

已知矩形窗的频谱函数是


w(n)=1WR(e-jω)=sinNω2sinω2ejω2
(3.7.20)

由变换的频移定理可得


w(n)=1·ej2πNnWRejω-2πN
w(n)=1·e-j2πNnWRejω+2πN
(3.7.21)

将频谱函数表达为


W(ejω)=W(ω)e-jωα


这样，利用上述关系可以得到用矩形窗的频谱幅度函数WR(ω)
来表示偶对称式汉宁窗的频谱函数


W(ω)≈12WR(ω)+14WRω-2πN+14WRω+2πN
(3.7.22)


该频谱特性如图3.7.4所示，是由3个互有频移的不同幅值的矩形窗频谱幅度函数相加而成，这将使副瓣大为抵消，能量更有效地集中在主瓣内，但却使主瓣加宽了一倍。



图3.7.4汉宁窗频谱


把N点偶对称表示的汉宁窗函数右移N/2点，可得单边表示的汉宁窗，其对应的频谱函数为


W(ejω)=12WR(ω)e-jN2ω+14WRω-2πNe-jN2ω-2πN+14WRω+2πNe-jN2ω+2πN
=12WR(ω)-14WRω-2πN-14WRω+2πNe-jN2ω
(3.7.23)


注意，与偶对称表示的汉宁加权不同，式(3.7.23)中两个频移的矩形窗幅度函数前的符号是负的。在离散傅里叶变换中，直接实现对输入序列x(n)的汉宁窗加权时，可不在时域进行x(n)与w(n)相乘，而在输出的频谱序列X(k)上进行线性组合来实现。已知汉宁加权后的离散傅里叶变换XW(k)为


XW(k)=∑N-1n=0w(n)x(n)e-j2πNkn
=∑N-1n=012-14ej2πNn-14e-j2πNnx(n)e-j2nNkn
=12∑N-1n=0x(n)e-j2πNkn-14∑N-1n=0x(n)e-j2nN(k-1)n-
14∑N-1n=0x(n)e-j2πN(k+1)n(3.7.24)

由式(3.7.24)可知，汉宁窗加权的离散傅里叶变换输出XW(k)是矩形窗加权的离散傅里叶变换X(k)的线性组合，即


XW(k)=12X(k)-14X(k-1)-14X(k+1)
=12X(k)-12［X(k-1)+X(k+1)］
(3.7.25)

这种输出X(k)的组合需要附加2N次复加及2N次右移(实现乘1/2)操作来实现。汉宁窗加权的这些特点，在快速离散傅里叶变换运算中特别受到注意。其实，汉宁窗是一族窗函数中的一个，这族窗函数的偶对称表示式为


w(n)=cosαnπNn=-N2，…,-1,0,1,…，N2
(3.7.26)

单边表示为


w(n)=sinαnπNn=0,1,2,…，N-1
(3.7.27)


一般取α=1,2,3,4。当α=1时，是余弦窗。α=2，是上文介绍的汉宁窗。当α越大时，窗函数cosα(x)序列越平滑，对应的频谱函数的副瓣电平就越下降，副瓣跌落越快速。但主瓣则变得越来越宽。图3.7.5给出了余弦窗及其频谱。图3.7.6给出了汉宁窗及其频谱。



图3.7.5余弦窗及其频谱幅度函数


【例3.9】绘出余弦窗及其频谱幅度函数。

import numpy as np

from scipy import signal, fft

import matplotlib.pyplot as plt

from matplotlib.ticker import MaxNLocator



# 计算余弦窗的wn值和幅度响应

N = 51# 余弦窗长度

wn = signal.windows.cosine(N)# 余弦窗的wn值

N0 = 2048

N1 = int(N0 / 2)

Ha = np.abs(fft.fft(wn, N0)) + 1e-10

Ha = Ha / np.max(Ha)

Ar = 20 * np.log10(Ha)

freq = np.linspace(0, 1, N1)



# 绘制余弦窗

fig, ax = plt.subplots()

ax.stem(wn, basefmt="")

ax.set_title('余弦窗')

ax.set_xlabel('n')

ax.set_ylim(［0, 1.5］)

plt.rcParams［'font.sans-serif'］ = ［'SimHei'］# 用来正常显示中文标签

fig.savefig('./win_cos1.png', dpi=500)



# 绘制余弦窗的幅度响应

fig, ax = plt.subplots()

ax.plot(freq, Ar［:N1］)

ax.grid()

ax.set_title('余弦窗的幅度响应')

ax.set_xlabel('k')

ax.set_xlabel(r'$ ＼omega / ＼pi $')

ax.set_ylabel(r'$ 20log_{10}| H(＼omega) | $')

ax.set_xlim(［0, 1］)

ax.set_ylim(［-120, 1］)

ax.xaxis.set_major_locator(MaxNLocator(11))

ax.yaxis.set_major_locator(MaxNLocator(13))

plt.rcParams［'axes.unicode_minus'］ = False# 用来显示负号



【例3.10】绘出汉宁窗及其幅度频谱函数。

import numpy as np

from scipy import signal, fft

import matplotlib.pyplot as plt

from matplotlib.ticker import MaxNLocator



# 计算汉宁窗的wn值和幅度响应

N = 51# 汉宁窗长度

wn = signal.windows.hann(N)# 汉宁窗的wn值

N0 = 2048

N1 = int(N0 / 2)

Ha = np.abs(fft.fft(wn, N0)) + 1e-10

Ha = Ha / np.max(Ha)

Ar = 20 * np.log10(Ha)

freq = np.linspace(0, 1, N1)



# 绘制汉宁窗

fig, ax = plt.subplots()

ax.stem(wn, basefmt="")

ax.set_title('汉宁窗')

ax.set_xlabel('n')

ax.set_ylim(［0, 1.5］)

plt.rcParams［'font.sans-serif'］ = ［'SimHei'］# 用来正常显示中文标签

fig.savefig('./win_han1.png', dpi=500)



# 绘制汉宁窗的幅度响应

fig, ax = plt.subplots()

ax.plot(freq, Ar［:N1］)

ax.grid()

ax.set_title('汉宁窗的幅度响应')

ax.set_xlabel('k')

ax.set_xlabel(r'$ ＼omega / ＼pi $')

ax.set_ylabel(r'$ 20log_{10}| H(＼omega) | $')

ax.set_xlim(［0, 1］)

ax.set_ylim(［-120, 1］)

ax.xaxis.set_major_locator(MaxNLocator(11))

ax.yaxis.set_major_locator(MaxNLocator(13))

plt.rcParams［'axes.unicode_minus'］ = False# 用来显示负号



汉宁窗的w(n)和幅度响应绘制如下： 



图3.7.6汉宁窗及其幅度频谱函数


4. 汉明窗(Hamming window)

汉明窗也称改进的升余弦窗。受汉宁窗启发，可设法调整相邻矩形窗频谱的大小以便更好地获取副瓣相消，为此可设加权窗函数为


w(n)=α+(1-α)cos2πNn 
(3.7.28)

对应的频谱函数为


W(ω)=αWR(ω)+12(1-α)WRω-2πN+WRω+2πN
(3.7.29)

改变可调整的比例系数α就可改变相邻3个矩形窗频谱幅度的大小，若选α=0.54，就是所谓 “汉明窗”可使副瓣电平有显著的改善。

汉明窗的偶对称表示为


w(n)=0.54+0.46cos2πNn，n=-N2,…，-1,0,1,…,N2
(3.7.30)


单边表示为


w(n)=0.54-0.46cos2πNn，n=0，1,2,…,N-1
(3.7.31)


所得到的频谱幅度函数为


W(ω)=0.54WR(ω)+0.23WRω-2πN+WRω+2πN
(3.7.32)

结果达到99.96%的能量集中在主瓣内，在与汉宁窗相等的主瓣宽度下，获得了更好的副瓣抑制。若选α=0.53856，则副瓣电平是-43dB。图3.7.7给出了汉明窗及其频谱，可以看到，在第一副瓣处出现了很深的凹陷。汉明加权后的离散傅里叶变换也可用X(k)来表达


Xw(k)=0.54X(k)-0.23［X(k-1)+X(k+1)］
(3.7.33)

它不能像汉宁窗那样用简单的右移来代替与系数作乘法，而需做实数乘复数的乘法。

【例3.11】绘出汉明窗及其幅度频谱函数。

import numpy as np

from scipy import signal, fft

import matplotlib.pyplot as plt

from matplotlib.ticker import MaxNLocator



# 计算汉明窗的wn值和幅度响应

N = 51# 汉明窗长度

wn = signal.windows.hamming(N)# 汉明窗的wn值

N0 = 2048

N1 = int(N0 / 2)

Ha = np.abs(fft.fft(wn, N0)) + 1e-10

Ha = Ha / np.max(Ha)

Ar = 20 * np.log10(Ha)

freq = np.linspace(0, 1, N1)



# 绘制汉明窗

fig, ax = plt.subplots()

ax.stem(wn, basefmt="")

ax.set_title('汉明窗')

ax.set_xlabel('n')

ax.set_ylim(［0, 1.5］)

plt.rcParams［'font.sans-serif'］ = ［'SimHei'］# 用来正常显示中文标签

fig.savefig('./win_ham1.png', dpi=500)



# 绘制汉明窗的幅度响应

fig, ax = plt.subplots()

ax.plot(freq, Ar［:N1］)

ax.grid()

ax.set_title('汉明窗的幅度响应')

ax.set_xlabel('k')

ax.set_xlabel(r'$ ＼omega / ＼pi $')

ax.set_ylabel(r'$ 20log_{10}| H(＼omega) | $')

ax.set_xlim(［0, 1］)

ax.set_ylim(［-100, 1］)

ax.xaxis.set_major_locator(MaxNLocator(11))

ax.yaxis.set_major_locator(MaxNLocator(11))

plt.rcParams［'axes.unicode_minus'］ = False# 用来显示负号



汉明窗的w(n)和幅度响应绘制如图3.7.7所示。



图3.7.7汉明窗及其幅度频谱函数


5. 布拉克曼窗(Blackman window)

布拉克曼窗也称二阶升余弦窗。汉宁、汉明加权都由3个中心频率不同的矩形窗频谱线性组合而成。布拉克曼利用更多的矩形窗频谱线性组合构成布拉克曼窗，其偶对称表示为


w(n)=∑K-1m=0amcos2πNmn，n=-N2,…,-1,0,1，…，N2
(3.7.34)

单边表示为


w(n)=∑K-1m=0(-1)mamcos2πNmn，n=0,1,2，…，N-1
(3.7.35)

单边表示的布拉克曼窗幅度频谱函数为


W(ω)=∑K-1m=0(-1)mam2WRω-2πNm+WRω+2πNm
(3.7.36)

这个窗函数中系数的选择应满足以下约束条件


∑K-1m=0am=1.0
(3.7.37)

因此，汉宁、汉明窗是a0和a1不为零，而其他系数都为零的布拉克曼窗。假如布拉克曼窗有K个非零的系数am，则其振幅频谱将由(2K-1)个中心频率不同的矩形窗频谱线性组合而成。显然，要使窗函数频谱的主瓣宽度窄，则K值不能选得很大。布拉克曼找到K=3时在ω=3.5(2π/N)及ω=4.5(2π/N)处出现零点的窗序列系数am(m=0,1,2)，其准确值及近似值如下式所示


a0=793810608≈0.42
a1=924018608≈0.5
a2=143018608≈0.08
(3.7.38)

按照式(3.7.38)的系数近似得出的窗函数，称为布拉克曼窗，其偶对称表示为


w(n)=0.42+0.50cos2πNn+0.08cos2πN2n
n=-N2,…,-1,0,1，…，N2
(3.7.39)

单边表示为


w(n)=0.42-0.50cos2πNn+0.08cos2πN2n
n=0,1,2,…，N-1
(3.7.40)

单边表示的频谱幅度函数为


W(ω)=0.42WR(ω)+0.25WRω-2πN+WRω+2πN+
0.04WRω-4πN+WRω+4πN
(3.7.41)

这样可以得到更低的副瓣，但主瓣宽度却进一步加宽到矩形窗的3倍。

【例3.12】绘出布拉克曼窗及其幅度频谱函数。

import numpy as np

from scipy import signal, fft

import matplotlib.pyplot as plt

from matplotlib.ticker import MaxNLocator



# 计算布拉克曼窗的wn值和幅度响应

N = 51# 布拉克曼窗长度

wn = signal.windows.blackman(N)# 布拉克曼窗的wn值

N0 = 2048

N1 = int(N0 / 2)

Ha = np.abs(fft.fft(wn, N0)) + 1e-10

Ha = Ha / np.max(Ha)

Ar = 20 * np.log10(Ha)

freq = np.linspace(0, 1, N1)



# 绘制布拉克曼窗

fig, ax = plt.subplots()

ax.stem(wn, basefmt="")

ax.set_title('布拉克曼窗')

ax.set_xlabel('n')

ax.set_ylim(［0, 1.5］)

plt.rcParams［'font.sans-serif'］ = ［'SimHei'］# 用来正常显示中文标签

fig.savefig('./win_blac1.png', dpi=500)



# 绘制布拉克曼窗的幅度响应

fig, ax = plt.subplots()

ax.plot(freq, Ar［:N1］)

ax.grid()

ax.set_title('布拉克曼窗的幅度响应')

ax.set_xlabel('k')

ax.set_xlabel(r'$ ＼omega / ＼pi $')

ax.set_ylabel(r'$ 20log_{10}| H(＼omega) | $')

ax.set_xlim(［0, 1］)

ax.set_ylim(［-150, 1］)

ax.xaxis.set_major_locator(MaxNLocator(11))

ax.yaxis.set_major_locator(MaxNLocator(16))

plt.rcParams［'axes.unicode_minus'］ = False# 用来显示负号



布拉克曼窗的w(n)和幅度响应绘制如图3.7.8所示。



图3.7.8布拉克曼窗及其幅度频谱函数


6. 最优化窗 

最优化窗是指在某种优化准则下得出的窗函数。

1) 高斯窗(Gaussian window)

已经知道，任何一种信号时宽T、频宽W之积TW满足下式


TW≥14π
(3.7.42)


在最小时宽频宽积TW=14π的优化准则下得出的是高斯波形，因此，可选择高斯波形作窗函数。高斯波形是一直延伸到无穷大的，但是若在波形3倍均方值的地方进行截断, 则误差很小。高斯窗如下式所示


w(n)=e-12anN/220≤|n|≤N-12
(3.7.43)


有限列长的高斯形窗相当于高斯形窗与矩形窗相乘。因此，有限列长的高斯形窗的频谱函数应是高斯窗频谱与矩形窗频谱的卷积，即


W(ω)=122πae-12ωa2*WR(ω)
≈N22πae-12ωa2，当a>2.5,ω很小时
(3.7.44)


a是标准差值的倒数，是频谱宽度的一种度量。图3.7.9为a=3时的高斯窗及其频谱。

【例3.13】绘出a=3时的高斯窗及其频谱幅度函数。

import numpy as np

from scipy import signal, fft

import matplotlib.pyplot as plt

from matplotlib.ticker import MaxNLocator



# 计算高斯窗的wn值和幅度响应

N = 51# 高斯窗长度

wn = signal.windows.gaussian(N, (N-1)/(2*3))# 高斯窗的wn值

N0 = 2048

N1 = int(N0 / 2)

Ha = np.abs(fft.fft(wn, N0)) + 1e-10

Ha = Ha / np.max(Ha)

Ar = 20 * np.log10(Ha)

freq = np.linspace(0, 1, N1)



# 绘制高斯窗

fig, ax = plt.subplots()

ax.stem(wn, basefmt="")

ax.set_title('高斯窗(a = 3)')

ax.set_xlabel('n')

ax.set_ylim(［0, 1.5］)

plt.rcParams［'font.sans-serif'］ = ［'SimHei'］# 用来正常显示中文标签

fig.savefig('./win_gauss1.png', dpi=500)



# 绘制高斯窗的幅度响应

fig, ax = plt.subplots()

ax.plot(freq, Ar［:N1］)

ax.grid()

ax.set_title('高斯窗的幅度响应(a = 3)')

ax.set_xlabel('k')

ax.set_xlabel(r'$ ＼omega / ＼pi $')

ax.set_ylabel(r'$ 20log_{10}| H(＼omega) | $')

ax.set_xlim(［0, 1］)

ax.set_ylim(［-100, 1］)

ax.xaxis.set_major_locator(MaxNLocator(11))

ax.yaxis.set_major_locator(MaxNLocator(11))

plt.rcParams［'axes.unicode_minus'］ = False# 用来显示负号



当N=51,α=3,σ=N-12α=506时，高斯窗的w(n)和幅度响应绘制如图3.7.9所示。



图3.7.9a=3时的高斯窗及其频谱幅度函数


2) 凯塞窗(Kaiser window)

凯塞窗的最优化准则是： 对于有限的信号能量，要求确定一个有限时宽T的信号波形，它使得频宽W内的能量最大。凯塞窗的定义如下


w(n)=I0πa1.0-nN/22I0(πa)，0≤|n|≤N2
(3.7.45)

式中，I0(x)=1+∑∞k=11k!x2k2是零阶第一类修正贝塞尔函数。凯塞窗的频谱幅度函数可近似为


W(ω)≈NI0(aπ)sinha2π2-Nω22a2π2-Nω22
(3.7.46)

上面公式中的πa是时宽、频宽积的一半。

【例3.14】绘出凯塞窗及其幅度响应。

import numpy as np

from scipy import signal, fft

import matplotlib.pyplot as plt

from matplotlib.ticker import MaxNLocator



# 计算凯塞窗的wn值和幅度响应

N = 51# 凯塞窗长度

wn = signal.windows.kaiser(N, 14)# 凯塞窗的wn值

N0 = 2048

N1 = int(N0 / 2)

Ha = np.abs(fft.fft(wn, N0)) + 1e-10

Ha = Ha / np.max(Ha)

Ar = 20 * np.log10(Ha)

freq = np.linspace(0, 1, N1)



# 绘制凯塞窗

fig, ax = plt.subplots()

ax.stem(wn, basefmt="")

ax.set_title('凯塞窗(πa = 14)')

ax.set_xlabel('n')

ax.set_ylim(［0, 1.5］)

plt.rcParams［'font.sans-serif'］ = ［'SimHei'］# 用来正常显示中文标签

fig.savefig('./win_kaiser1.png', dpi=500)



# 绘制凯塞窗的幅度响应

fig, ax = plt.subplots()

ax.plot(freq, Ar［:N1］)

ax.grid()

ax.set_title('凯塞窗的幅度响应(πa = 14)')

ax.set_xlabel('k')

ax.set_xlabel(r'$ ＼omega / ＼pi $')

ax.set_ylabel(r'$ 20log_{10}| H(＼omega) | $')

ax.set_xlim(［0, 1］)

ax.set_ylim(［-150, 1］)

ax.xaxis.set_major_locator(MaxNLocator(11))

ax.yaxis.set_major_locator(MaxNLocator(16))

plt.rcParams［'axes.unicode_minus'］ = False# 用来显示负号



当πa=14时，凯塞窗的w(n)和幅度响应绘制如图3.7.10所示。



图3.7.10凯塞窗及其幅度响应(πa=14)




图3.7.10(续)


3.8习题◆

1. 图3.8.1的序列x～1(n)是具有周期为4的周期性序列。试确定傅里叶级数的系数X～1(k)。

2. 在图3.8.2中表示了两个周期都为6的周期性序列，确定这两个序列的周期卷积的结果x～3(n),并画出草图。




图3.8.1题1图




图3.8.2题2图




3. 在图3.8.3中表示了几个周期性序列x(n),这些序列可以按离散傅里叶级数表示为


x～(n)=1N∑N-1k=0X～(k)ej2πNkn


(1) 对于哪些序列可以通过选择时间起点使得所有X～(k)是实数？

(2) 对于哪些序列可以通过选择时间起点使得除X～(0)外的所有X～(k)是虚数？

(3) 对于哪些序列能够做到X～(k)=0,k=±2,±4,…？



图3.8.3题3图


4. 试证明下面列出的周期序列离散傅里叶级数的对称特性。在证明中，可以利用离散傅里叶级数的定义及任何已证明的性质。例如在证明性质(3)时可以利用性质(1)和(2)。

序列离散傅里叶级数

(1) x～*(n)X～*(-k)

(2) x～*(-n)X～*(k)

(3) Re［x～(n)］X～e(k)

(4) jIm［x～(n)］X～o(k)

根据以上证明的性质，证明对于实数周期序列x～(n),离散傅里叶级数的下列对称特性成立： 

(1) Re［X～(k)］=Re［X～(-k)］

(2) jIm［X～(k)］=-Im［X～(k)］

(3) |X～(k)|=|X～(-k)|

(4) arg［X～(k)］=-arg［X～(-k)］

5. 求下列序列的DFT： 

(1) 11-1-1

(2) 1j-1-j

(3) x(n)=cn,0≤n≤N-1

(4) x(n)=sin2πnN,0≤n≤N-1

6. 在图3.8.4中表示了有限长序列x(n)，画出序列x1(n)和x2(n)的草图(注意： x1(n)和x2(n)是x(n)圆周移位的两个点)。


x1(n)=x((n-2))4R4(n)
x2(n)=x((-n))4R4(n)



图3.8.4题6图


7. 在图3.8.5中表示了两个有限长序列，试画出它们的6点圆周卷积。



图3.8.5题7图


8. 有限长序列的DFT对应于序列在单位圆上的z变换的取样。例如一个10点序列x(n)的DFT对应于图3.8.6表示的10个等间隔点上X(z)的采样。希望找出在图3.8.7所示的围线上X(z)的等间隔采样，即X(z)z=0.5ej2πk10+π10。证明如何修改x(n)以获得一个序列x1(n)致使x1(n)的DFT对应于所希望的X(z)的采样。




图3.8.6题8图1




图3.8.7题8图2




9. 列长为8的一个有限长序列具有8点离散傅里叶变换X(k)，如图3.8.8所示。列长为16点的一个新的序列y(n)定义为


y(n)=
xn2，n为偶数
0，n为奇数




图3.8.8题9图1


试从图3.8.9的几个图中选出相当于y(n)的16点离散傅里叶变换序列图。



图3.8.9题9图2



10. 图3.8.10表示一个4点序列x(n)。

(1) 试绘出x(n)同x(n)线性卷积略图； 

(2) 试绘出x(n)同x(n)4点圆周卷积略图； 

(3) 试绘出x(n)同x(n)10点圆周卷积略图； 

(4) 若x(n)同x(n)的某个N点圆周卷积与其线性卷积相同，试问此时N点的最小值是多少？



图3.8.10题10图