数字信号处理四：离散时间信号与系统的频域分析

数字信号处理是数字信号处理领域中专业人员必须掌握的一项基本技能。通过本内容的学习，可以掌握信号变换的核心概念及其应用，理解数字信号处理的基本原理和方法。重点内容包括：
信号变换概述

• 傅里叶变换（FT）、离散傅里叶变换（DFT）、Z变换（ZT）等变换在信号分析和处理中的重要性。
• 时域和频域之间的转换关系，以及不同变换在信号处理中的应用场景。
  DFT的计算
• 使用快速傅里叶变换（FFT）算法高效计算DFT。
• 通过补零等方法提高DFT计算效率，以及在Python numpy库中的实现方法。
  DFT与DTFT的关系
• DFT是DTFT在区间[0, 2π]上的N点等间隔采样。
• 通过具体例子验证了DFT与DTFT之间的关系，展示了DFT在不同点数下的频谱特性。

一、信号变换概述

信号是数字信号处理领域中最基础且最重要的核心概念。数字信号变换技术构成了数字信号处理领域中最基础且不可或缺的核心技术手段。因此，数字信号变换技术已成为数字信号处理领域中专业人员必须掌握的一项基础且不可或缺的技术技能。简而言之，数字信号变换技术是通过数学变换，将一个域中的信号映射到另一个域，以实现处理和操作的目的。常用的数字信号变换主要有：傅里叶变换（DTFT）、离散傅里叶变换（DFT）、Z变换（ZT）等。这些变换，都有着各自的理论和其应用背景。

二、傅里叶变换

傅里叶变换在信号分析与处理领域扮演着关键角色。有限长序列作为离散信号的一种，在数字信号处理中占据着至关重要的地位。离散傅里叶变换在理论上具有重要意义，同时，快速傅里叶变换提供了高效的计算方法。因此，在数字信号处理的各种运算方法中，傅里叶变换正发挥着核心作用。

2.1 傅里叶变换的几种形式

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

非周期连续时间信号x(t)的傅里叶变换

可以表示为:

逆变换为:

在这里，

是模拟角频率。可以看出，时域上连续的函数确实会导致频域呈现出非周期性的谱线，而时域上的非周期特性则会引致频域上出现连续的谱。由此可见，非周期连续时间函数必然对应于一个非周期连续的频域变换函数。

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

周期信号x(t)具有周期为T的连续性，其傅里叶变换在离散频域上具有明确的数学表达式，具体表示为：X(k) = \sum_{n=0}^{N-1} x(n) e^{-j2\pi kn/N}。

逆变换为

这就是经常称之为傅里叶级数的变换形式。在这里，

也是角频率的模拟量。通过观察可知，在时域中，连续的函数对应频率域中的非周期谱，而离散的频域函数则导致时域函数呈现周期性特征。
结论：连续的周期时间函数对应于一个非周期的离散频域变换函数。

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

非周期离散时间信号x(n)的傅里叶变换

可以表示为:

逆变换为:

在这里， w是数字频率，它和模拟角频率的关系为

可以看出，时域的取样在频域上表现为周期延拓，这直接源于时域函数的非周期特性。时域函数的非周期性特性直接导致频域中出现离散的谱线。非周期离散时间函数在频域上对应于一个连续的周期变换函数，这种对应关系直接反映了时域和频域之间的内在联系。

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

周期离散时间信号

的傅里叶变换——离散傅里叶级数，可以表示为:

逆变换为:

可以看出，时域的取样等同于频域的周期延拓，此外，时域函数的周期性导致频域的离散谱。总结：周期离散时间函数等同于单一周期离散频域变换函数。

2.2 离散傅里叶变换（DFT）

离散傅里叶级数属于周期序列，尽管如此，仍然存在计算上的困难。尽管离散傅里叶级数属于周期序列，但其仅含有N个独立的数值，因此可以通过对有限长序列进行周期延拓来揭示其特性。对于一个长度为N的有限长序列x(n)，即x(n)仅在n=0到n=(N-1)区间内具有非零值，而在其余点的值均为零，即

把序列x(n)以N为周期进行周期延拓得到周期序列 ，则有

所以，有限长序列x(n)的离散傅里叶变换（DFT）为

逆变换为

2.3 快速傅里叶变换（FFT）

在信号处理领域，离散傅里叶变换（DFT）具有关键地位，它是实现信号相关、滤波以及谱估计等核心操作的基础。当处理规模较大时，执行一个N点的DFT需要进行N²次复数乘法运算和N(N-1)次复数加法运算，其计算复杂度较高。1965年，J.W. Cooley和J.W. Tukey开创性地提出了快速傅里叶变换（FFT）算法，极大地简化了DFT的计算过程。

因子的周期性和对称性特征，构建了高效的快速算法框架，即为快速傅里叶变换(FFT)方法的理论基础。它并不是一种全新的独立变换方法，而是通过优化计算过程而显著降低计算复杂度的高效算法。这种高效的快速算法主要基于

周期性数据序列的频域特性，通过离散傅里叶变换的快速计算实现，从而在信号处理和数据分析领域得到了广泛应用。

该算法通过不断将长序列的DFT分解为几个短序列的DFT，从而降低运算次数。最常用的是基2FFT算法。该算法分为两类，即时间抽取法（Decimation－In－Time，简称DIT-FFT）和频率抽取法（Decimation－In－Frequency，简称DIF-FFT）。在numpy库的fft模块中，提供了丰富的数学函数，用于计算数据的快速傅里叶变换。这些函数包括fft、ifft、fft2、ifft2、fftshift和ifftshift等。当处理的数据长度为2的幂次时，采用基2算法可以显著提高计算速度。因此，建议尽可能使处理的数据长度为2的幂次，或者通过补零使其成为2的幂次。

三、离散时间LTI系统的频率响应特性分析

对于因果稳定的离散时间系统，系统的频率响应定义为：

其中，

称为离散时间系统的幅频特性；

称为离散时间系统的相频特性；

是以2π为周期的周期函数。因此，只要分析

在分析范围内的情况时，可以得出系统的整体频率特性。

语句格式：
w, H = signal.dfreqresp(sys,whole=False)
其中返回值w包含

范围内的N个频率等分点；

返回值H则是离散时间系统频率响应

在指定频率区间内计算出N个频率点的值。此外，函数还可以通过调用w, H = signal.dfreqresp(sys,whole=True)来实现另一种调用方式。主要区别在于，角频率范围扩展为整个频率轴。

在函数调用时，sys参数代表系统。可以通过scipy库中的signal模块中的TransferFunction函数来实现。

该系统通过信号库中的TransferFunction函数实现。其中，b和a分别代表分子和分母多项式，这些多项式均按z的正幂次排列，最后一项为常数项，并且以系数向量的形式表示。dt为采样时间参数。

绘制系统

的频率响应曲线。
程序代码示例如下：

    from scipy import signal
    import matplotlib.pyplot as plt
    import numpy as np
    
    b=[1,-0.96,0.9028]
    a=[1,-1.56,0.8109]
    sys = signal.TransferFunction(b, a,dt=0.05)
    w, H = signal.dfreqresp(sys,whole=True)
    
    plt.figure(figsize=(15,15))
    plt.subplot(211)
    plt.plot(w, np.abs(H), "b")
    plt.xlabel('ω(rad/s)')
    plt.ylabel('Magnitude')
    plt.title('|H(e^jw)|')
    plt.grid('on')
    
    plt.subplot(212)
    plt.plot(w, np.angle(H), "r")
    plt.xlabel('ω(rad/s)')
    plt.ylabel('Phase')
    plt.title('φ(w)')
    plt.grid('on')
    plt.show()

程序运行结果如下:

。

已知

，绘制系统的幅频特性和相频特性曲线。程序代码示例如下：

    from scipy import signal
    import matplotlib.pyplot as plt
    import numpy as np
    
    b=[1,0,0,0,0,0,0,0,-1]
    a=[1,0,0,0,0,0,0,0,0]
    sys = signal.TransferFunction(b,a,dt=0.05)
    w, H = signal.dfreqresp(sys,whole=True)
    
    plt.figure(figsize=(12,12))
    plt.subplot(211)
    plt.plot(w, np.abs(H), "b")
    plt.xlabel('ω(rad/s)')
    plt.ylabel('Magnitude')
    plt.title('|H(e^jw)|')
    plt.grid('on')
    
    plt.subplot(212)
    plt.plot(w, np.angle(H), "r")
    plt.xlabel('ω(rad/s)')
    plt.ylabel('Phase')
    plt.title('φ(w)')
    plt.grid('on')
    plt.show()

程序运行结果如图下。

基于数字信号处理中学习过的函数，我们可以绘制此系统函数的零极点分布图，具体分布情况如下。

通过图形可以看到，系统包含8个零点，均匀地分布在单位圆上，具有1个8阶极点，位于原点。该系统被称作梳状滤波器。

，其傅里叶变换为

，绘制其幅度谱和相位谱。

程序代码示例如下：

    from scipy import signal
    import matplotlib.pyplot as plt
    import numpy as np
    
    b=[1,0,0,0,-1]
    a=[1,-1,0,0,0]
    sys = signal.TransferFunction(b,a,dt=0.05)
    w, H = signal.dfreqresp(sys,whole=True)
    plt.figure(figsize=(10,10))
    plt.subplot(211)
    plt.plot(w, np.abs(H), "b")
    plt.xlabel('ω(rad/s)')
    plt.ylabel('Magnitude')
    plt.title('|X(e^jw)|')
    plt.grid('on')
    
    plt.subplot(212)
    plt.plot(w, np.angle(H), "r")
    plt.xlabel('ω(rad/s)')
    plt.ylabel('Phase')
    plt.title('arg[X(e^jw)]')
    plt.grid('on')
    plt.show()

程序运行结果如下图所示。

按照同样的方法，我们可以得到 R5(n)、R6(n)的幅度谱和相位谱，如下图所示。

四、DFT的计算

对长度为M的有限长序列进行N点离散傅里叶变换（其中N≥M），可以直接调用numpy库中fft模块中的fft函数和ifft函数。如果N等于2的整数次幂，可以采用快速傅里叶变换算法来加速计算。

用于计算DFT正变换的函数是fft，其调用格式为np.fft.fft(x)，其中x代表时间域的序列。用于计算DFT逆变换的函数是ifft，其调用格式为np.fft.ifft(X)，其中X代表频率域的频谱值。

已知

，求 (n)的4点和8点DFT，并绘制幅度谱。
程序代码示例如下：

    import matplotlib.pyplot as plt
    import numpy as np
    from matplotlib.font_manager import FontProperties
    
    font_set = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=12)
    N=4     # DFT的长度
    x=[1,1,1,1]
    N1=len(x)        # 序列x(n)的原始长度
    xn=np.zeros(N)
    xn[0:N1]=x        # 通过补零，使得x(n)的长度等于N
    
    XK = np.fft.fft(xn)   # 计算FFT
    print('X(k)=',XK)    # 打印DFT的结果
    n=np.arange(N)
    k=np.arange(N)
    
    plt.figure(figsize=(12,12))
    plt.subplot(211)
    plt.stem(n,xn, "r")
    plt.xlabel('n')
    plt.ylabel('x(n)')
    plt.title('R4(n)')
    plt.grid('on')
    
    plt.subplot(212)
    plt.stem(k,np.abs(XK))
    plt.xlabel('k')
    plt.ylabel('|X(k)|')
    plt.title('4点DFT结果|X(k)|',fontproperties=font_set)
    plt.grid('on')
    plt.show()

程序运行结果如下图所示。

将程序中N的值改为8，可得8点DFT的结果，如下图所示。

如果要求x(n)的L点DFT，只需将程序中N的值改为相应的数值即可。
已知

，求x(n)的6点DFT结果X(k)。
程序代码示例如下：

    import numpy as np
    xn=[2,1,5,4,5,1]
    XK = np.fft.fft(xn)      #计算6点DFT
    print(XK)

程序运行结果为：[18.+0.j -6.+0.j 0.+0.j 6.+0.j 0.+0.j -6.+0.j]

根据计算结果可以看出，在这个例题中，因为 是实数序列且具有偶对称性，所以X(k)同样具有实数属性和偶对称性。已知条件为 ，求X(k)的6点逆离散傅里叶变换结果。程序代码示例如下：

    import numpy as np
    XK=[18,-6,0,6,0,-6]
    xn = np.fft.ifft(XK) #计算6点IDFT
    print(xn)

程序运行结果为：[2.+0.j 1.+0.j 5.+0.j 4.+0.j 5.+0.j 1.+0.j]

从结果可知，

在上一题的示例中，我们观察到频谱关于N/2点呈现对称性。为了将频域数据中的直流分量（即频率为0处的值）平移到频谱的中间位置，可以应用函数np.fft.fftshift，其调用格式为：np.fft.fftshift(XK)。在进行逆变换之前，需要将频谱再平移回去，这一步骤需要用到函数np.fft.ifftshift，其调用格式为：np.fft.ifftshift(XK)。已知条件为：

计算x(n)的6点离散傅里叶变换结果X(k)，并将其频域的直流成分调整至频谱的中心位置。

    import numpy as np
    xn=[2,1,5,4,5,1]
    XK = np.fft.fft(xn)  #计算6点DFT
    XK_shift=np.fft.fftshift(XK)
    print(XK)
    print(XK_shift)

程序运行结果为：
[18.+0.j -6.+0.j 0.+0.j 6.+0.j 0.+0.j -6.+0.j]
[ 6.+0.j 0.+0.j -6.+0.j 18.+0.j -6.+0.j 0.+0.j]

五、DFT与DTFT的关系

x(n)的N点离散傅里叶变换（DFT）等于其在区间[0, 2π]上进行N点等间隔采样得到的离散时间傅里叶变换（DTFT）。通过例题可以验证这一概念。

已知

，其傅里叶变换为

绘制其DTFT的幅度谱和N点DFT的幅度谱。（N分别取8、16、32、64和128）当N=8时，程序示例代码如下：

    from scipy import signal
    import matplotlib.pyplot as plt
    import numpy as np
    
    b=[1,0,0,0,-1]
    a=[1,-1,0,0,0]
    sys = signal.TransferFunction(b,a,dt=0.05)
    w, H = signal.dfreqresp(sys,whole=True)
    N=8               # DFT的长度
    x1=[1,1,1,1]
    N1=len(x1)
    x=np.zeros(N)
    x[0:N1]=x1
    n1=np.arange(N)
    n=(2*np.pi)/N
    k=n*n1
    XK = np.fft.fft(x)    #计算FFT
    plt.figure(figsize=(8,4))
    
    plt.plot(w, np.abs(H), "r")
    plt.xlabel('ω(rad/s), k')
    plt.ylabel('Magnitude')
    plt.title('|X(e^jw)|, |X(k)|')
    plt.grid('on')
    plt.stem(k,np.abs(XK))
    plt.show()

程序运行结果如下图所示。

当N分别取16，32，64，128时，结果如下图所示。