信号处理基础：数字信号处理基础_（12）.谱估计技术

谱估计技术

引言

数字信号处理领域中的一个核心分支是谱估计技术。其主要目标是从有限观测数据中推断出信号的频谱特性。这种分析方法在多个领域有着广泛的应用实例：例如，在通信系统中用于优化信息传输，在雷达系统中用于目标检测，在声学研究中用于声音传播特性分析以及在生物医学工程领域用于生理信号特征提取等。借助频谱估计方法，则可深入探究信号在不同频率上的能量分布情况。

1. 傅里叶变换与频谱估计

在频谱估计中，傅里叶变换被视为一种基础工具，在时域与频域之间建立联系的能力使它成为分析信号的重要手段之一。它能够将信号从时域转换到频域的表现形式，在工程应用中具有广泛的应用价值。其中最常见的两种类型分别是离散傅里叶变换（DFT）和快速傅里叶变换（FFT）。

1.1 离散傅里叶变换（DFT）

作为处理离散时间信号并将其转换为频率域表示的关键数学工具之一，在工程学和计算机科学中广泛应用于信号分析和处理

X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} kn}

其中，x[n] 是离散时间信号，X[k] 是其频谱，N 是信号的长度，j 是虚数单位。

1.2 快速傅里叶变换（FFT）

快速傅里叶变换（FFT）可被视为DFT的一种优化算法。该算法通过降低计算规模，在实际应用中显著提升了频谱估计的可行性。FFT的时间复杂度为 O(N \log N) ，而DFT的时间复杂度则为 O(N^2) 。

1.3 Python 代码示例

下面是一个使用 Python 和 NumPy 库进行 FFT 的示例：

    import numpy as np
    import matplotlib.pyplot as plt
    
    # 生成一个示例信号
    fs = 1000  # 采样频率
    t = np.arange(0, 1, 1/fs)  # 时间向量
    f1 = 50  # 信号频率1
    f2 = 120  # 信号频率2
    x = 0.7 * np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)  # 合成信号
    
    # 计算 FFT
    X = np.fft.fft(x)
    X magn = np.abs(X)  # 频谱的幅值
    f = np.fft.fftfreq(t.shape[-1], d=1/fs)  # 频率向量
    
    # 绘制频谱图
    plt.figure(figsize=(10, 6))
    plt.plot(f, X_magn)
    plt.title('频谱图')
    plt.xlabel('频率 (Hz)')
    plt.ylabel('幅值')
    plt.xlim(0, fs/2)  # 只绘制正频率部分
    plt.grid(True)
    plt.show()

1.4 代码解释

生成示例信号 ：

fs 被设定为采样率，在本研究中被设定为 f_s = 1\,\mathrm{kHz}。
t 被定义为空间向量 \mathbf{t} 或者时间序列 \mathbf{T} ，其取值范围是从 t = 0 到 t = T （此处 T = 1\,\mathrm{s}），其中取样间隔被设定为 t_\mathrm{step} = \frac{1}{f_s}。
f_1$ 和 $f_2$ 分别表示两个待分析信号的基本频率参数值（此处分别对应于 $f_1 =50\,\mathrm{Hz}$ 和 $f_2=120\,\mathrm{Hz}$）。 x 被定义为空间向量 \mathbf{x}$ ，它是由两个不同频率的正弦波叠加而成。

计算 FFT ：

• 使用快速傅里叶变换算法对信号x进行处理。

• 通过numpy库中的abs函数获取频谱中各频率分量的幅度。

• 计算相应的频率值，并确定其对应的采样时间间隔。

绘制频谱图 ：

借助 matplotlib 绘制频谱图形会更加直观有效地呈现信号特性

2. 周期图法

周期图法是一种以傅里叶变换为基础的非参数型谱估计方法，在分析时间序列数据时被广泛应用。它通过计算信号的时间自相关序列并将其转换为频域信息来推断信号的能量分布特性。

2.1 周期图法的原理

周期图法的基本步骤如下：

1. 计算信号的自相关函数。
2. 对自相关函数进行傅里叶变换，得到频谱估计。

自相关函数的定义为：

R_{xx}[\tau] = \frac{1}{N} \sum_{n=0}^{N-1-\tau} x[n] x[n+\tau]

其中，\tau 是时延，N 是信号的长度。

2.2 Python 代码示例

下面是一个使用 Python 和 SciPy 库进行周期图法谱估计的示例：

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import periodogram
    
    # 生成一个示例信号
    fs = 1000  # 采样频率
    t = np.arange(0, 1, 1/fs)  # 时间向量
    f1 = 50  # 信号频率1
    f2 = 120  # 信号频率2
    x = 0.7 * np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)  # 合成信号
    
    # 计算周期图
    f, Pxx = periodogram(x, fs)
    
    # 绘制周期图
    plt.figure(figsize=(10, 6))
    plt.plot(f, Pxx)
    plt.title('周期图')
    plt.xlabel('频率 (Hz)')
    plt.ylabel('功率谱密度 (dB/Hz)')
    plt.grid(True)
    plt.show()

2.3 代码解释

生成示例信号 ：

• fs 被设定为采样频率，并设置成一个值等于一千赫兹。

• 时间向量`t'是由零开始到一秒，并采用周期长度等于一除以fs的时间间隔进行采样的。

• 这两个信号的频率分别被定义成五个百赫与一百二十百赫。

• 合成信号x则由这两个正弦波组成。

计算周期图 ：

 * `periodogram(x, fs)` 计算信号 x 的周期图，返回频率向量 `f` 和功率谱密度 `Pxx`。

绘制周期图 ：

通过 matplotlib 库绘制周期图，并对图像进行标注处理

3. Welch 法

Welch 方法是一种优化型周期图法，在分析过程中将信号均分为若干重叠的小块，并对每一块进行快速傅里叶变换；接着计算各小块频谱并取平均值；这种方法能够有效降低周期图估计中的方差。

3.1 Welch 法的原理

Welch 法的基本步骤如下：

1. 信号被划分为多个连续的时间区间。
2. 为各个时间片段计算其傅里叶转换。
3. 确定各时间段内的功率谱密度。
4. 统计所有时间片段内频率响应的数据并求其平均。

3.2 Python 代码示例

下面是一个使用 Python 和 SciPy 库进行 Welch 法谱估计的示例：

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import welch
    
    # 生成一个示例信号
    fs = 1000  # 采样频率
    t = np.arange(0, 1, 1/fs)  # 时间向量
    f1 = 50  # 信号频率1
    f2 = 120  # 信号频率2
    x = 0.7 * np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)  # 合成信号
    
    # 计算 Welch 法
    f, Pxx = welch(x, fs, nperseg=256)
    
    # 绘制 Welch 法谱估计
    plt.figure(figsize=(10, 6))
    plt.plot(f, Pxx)
    plt.title('Welch 法谱估计')
    plt.xlabel('频率 (Hz)')
    plt.ylabel('功率谱密度 (dB/Hz)')
    plt.grid(True)
    plt.show()

3.3 代码解释

生成示例信号 ：

• 采样频率被定义为 fs，并赋值为 1000 Hz。

• 时间向量 t 被构造出来，在时间范围从 0 秒到 1 秒之间均匀分布。

• 时间向量 t 的间隔被设定为了 dt = 1/fs 秒。

• 时间向量 t 的长度被确定为了 N 点。

• 基于时间向量 t 进行计算得到的结果存储在变量 y 中。

• 合成信号 x 被构造出来。

• 合成信号 x 被用来分析两个不同频率的正弦波特性。

计算 Welch 法 ：

welch(x, fs, nperseg=256) 采用 Welch 方法计算信号 x 的功率谱密度估计值 Pxx。其中 nperseg 设定为 256 表示每个子段的长度。

绘制 Welch 法谱估计 ：

• 通过调用 matplotlib 绘制 Welch 法谱估计图。
• 分别利用 plt 的 xlabel 和 ylabel 方法设置 x轴为‘频率（Hz）’、y轴为‘功率谱密度（dB/Hz）’的标签。

4. 最小二乘法

该方法属于一种基于参数化的频谱估计技术。其核心原理在于通过构建信号模型并对其进行拟合来推断信号频谱。在实际应用中，我们通常涉及AR（自回归）、MA（滑动平均）以及ARMA（自回归滑动平均）等典型模型。

4.1 AR 模型

AR 模型基于线性组合和加性白噪声生成信号。该模型中参数 p 代表所使用的历史数据点数量。

4.2 Python 代码示例

下面是一个使用 Python 和 SciPy 库进行 AR 模型谱估计的示例：

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy import signal
    
    # 生成一个示例信号
    fs = 1000  # 采样频率
    t = np.arange(0, 1, 1/fs)  # 时间向量
    f1 = 50  # 信号频率1
    f2 = 120  # 信号频率2
    x = 0.7 * np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)  # 合成信号
    
    # 计算 AR 模型
    ar_order = 16  # AR 模型的阶数
    ar_coeffs = signal.arburg(x, ar_order)
    w, h = signal.freqz(1, ar_coeffs[0], whole=True)
    
    # 计算功率谱密度
    Pxx = (1.0 / (np.abs(h) ** 2))
    
    # 绘制 AR 模型谱估计
    f = w / (2 * np.pi) * fs
    plt.figure(figsize=(10, 6))
    plt.plot(f, 10 * np.log10(Pxx))
    plt.title('AR 模型谱估计')
    plt.xlabel('频率 (Hz)')
    plt.ylabel('功率谱密度 (dB/Hz)')
    plt.grid(True)
    plt.show()

4.3 代码解释

生成示例信号 ：

• fs 被设定为采样频率值（sampling frequency），其具体数值定位于1千赫兹（kHz）。

• 时间变量 t 被定义在一个持续时间为一秒钟的时间区间内（second），其采样周期（sampling period）等于 T_s = \frac{1}{f_s}。

• 参数 $f_1$ 和 $f_2$ 分别代表两个待分析信号的基本频率分量（fundamental frequency components），它们被赋以具体的数值值：50Hz 和 120Hz。

• 数据序列 $x$ 被构造作为一个合成波形（composite waveform），它是通过将两个不同频率的正弦波叠加而成（superposition）而得到的结果。

计算 AR 模型 ：

• 使用 signal 库中的 arburg 函数通过信号 x 估计其 AR 模型参数
  其中参数 ar_order 表示模型的阶数 并设定为 16

• 调用 signal 库中的 freqz 函数基于估计得到的 AR 参数 计算其频率响应
  并返回相应的频率向量 w 和系统函数 h

计算功率谱密度 ：

 * `Pxx = (1.0 / (np.abs(h) ** 2))` 计算功率谱密度。

绘制 AR 模型谱估计 ：

• 借助 Matplotlib 绘制 AR 模型的谱估计图。
• 通过公式将频率向量转换为 f = \frac{w}{2\pi} \cdot fs。
• 设置图像的x轴与y轴标签分别为'频率 (Hz)'与'功率谱密度 (dB/Hz)'。

5. 高分辨率谱估计方法

高频率谱估计技术通过提升频谱分辨力来优化谱估计的效果。主要采用MUSIC（Multiple Signal Classification）算法及ESPRIPT（Estimation of Signal Parameters via Rotational Invariance Techniques）算法进行频谱分析。

5.1 MUSIC 算法

MUSIC算法利用子空间方法实现对信号频率的有效估计，在处理过程中主要通过对观测数据的空间分解来分离出不同的信号成分。具体实施步骤如下：

1. 计算信号的协方差矩阵。
2. 对协方差矩阵进行特征值分解。
3. 估计信号的频率。

5.2 Python 代码示例

下面是一个使用 Python 和 SciPy 库进行 MUSIC 算法谱估计的示例：

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy import linalg
    
    # 生成一个示例信号
    fs = 1000  # 采样频率
    t = np.arange(0, 1, 1/fs)  # 时间向量
    f1 = 50  # 信号频率1
    f2 = 120  # 信号频率2
    x = 0.7 * np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)  # 合成信号
    
    # 计算协方差矩阵
    Rxx = np.cov(x)
    
    # 特征值分解
    eigenvalues, eigenvectors = linalg.eig(Rxx)
    
    # 选择信号子空间和噪声子空间
    signal_eigenvectors = eigenvectors[:, :2]
    noise_eigenvectors = eigenvectors[:, 2:]
    
    # 定义频率搜索范围
    f_search = np.linspace(0, fs/2, 500)
    
    # 计算 MUSIC 谱
    P_music = np.zeros_like(f_search)
    for f in f_search:
    e = np.exp(-1j * 2 * np.pi * f * t)
    P_music[f_search == f] = 1 / np.sum(np.abs(np.dot(noise_eigenvectors.T, e)) ** 2)
    
    # 绘制 MUSIC 谱估计
    plt.figure(figsize=(10, 6))
    plt.plot(f_search, 10 * np.log10(P_music))
    plt.title('MUSIC 谱估计')
    plt.xlabel('频率 (Hz)')
    plt.ylabel('功率谱密度 (dB/Hz)')
    plt.grid(True)
    plt.show()

5.3 代码解释

生成示例信号 ：

• 采样频率被定义为 fs，并赋值为 1 \times 10^{3} 赫兹。

• 时间向量 t 被定义为从零开始一直到一秒，并采用步长设定为 1/fs 的方式生成。

• 这两个信号的频率分别设定为 50 赫兹和 120 赫兹。

• 合成信号 x 被定义为此两个正弦波的叠加结果。

计算协方差矩阵 ：

 * `np.cov(x)` 计算信号 x 的协方差矩阵 `Rxx`。

特征值分解 ：

linalg.eig(Rxx) 执行协方差矩阵的特征值分解，并生成对应的特征值 eigenvalues 和特征向量 eigenvectors。

选择信号子空间和噪声子空间 ：

signal_eigenvectors 属于信号子空间中的特征向量。
* noise_eigenvectors 属于噪声子空间中的特征向量。

定义频率搜索范围 ：

 * `f_search` 是频率搜索范围，从 0 到 fs/2，共 500 个点。

计算 MUSIC 谱 ：

 * `P_music` 是 MUSIC 谱的估计结果。
 * 对于每个频率点 `f`，计算对应的谱值。

绘制 MUSIC 谱估计 ：

• 通过 matplotlib 来绘制 MUSIC 谱估计图。
• 指定 plt.xlabel('频率 (Hz)') 和 plt.ylabel('功率谱密度 (dB/Hz)') 以分别给 x 轴和 y 轴命名。

6. 相关性与互谱估计

相关性和互谱估计是频谱估计中的关键概念，在研究两个信号之间相互作用方面发挥重要作用。它们广泛应用于通信、雷达以及生物医学等多个领域，并且特别适用于解决信号检测与分离问题。

6.1 互谱估计

主要通过计算两个信号的互相关函数，并结合傅里叶变换技术来推断两者间的频谱关联

R_{xy}[\tau] = \frac{1}{N} \sum_{n=0}^{N-1-\tau} x[n] y[n+\tau]

在处理过程中，默认情况下我们假设输入数据为两路信号x[n]和y[n]（n=0,1,…,N-1），其中\tau表示时间延迟参数，并且数据序列的长度为N。互谱估计则能够揭示两路信号在不同频率点上的相位关系及其能量分布情况

6.2 Python 代码示例

下面是一个使用 Python 和 SciPy 库进行互谱估计的示例：

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import csd
    
    # 生成两个示例信号
    fs = 1000  # 采样频率
    t = np.arange(0, 1, 1/fs)  # 时间向量
    f1 = 50  # 信号频率1
    f2 = 120  # 信号频率2
    x = 0.7 * np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)  # 合成信号 x
    y = 0.5 * np.sin(2 * np.pi * f1 * t + np.pi/4) + np.sin(2 * np.pi * f2 * t + np.pi/8)  # 合成信号 y
    
    # 计算互谱
    f, Pxy = csd(x, y, fs)
    
    # 绘制互谱图
    plt.figure(figsize=(10, 6))
    plt.plot(f, np.abs(Pxy))
    plt.title('互谱图')
    plt.xlabel('频率 (Hz)')
    plt.ylabel('互谱密度 (dB/Hz)')
    plt.xlim(0, fs/2)  # 只绘制正频率部分
    plt.grid(True)
    plt.show()

6.3 代码解释

生成示例信号 ：

fs 被设定为采样频率并被指定为1000赫兹。
时间向量`t'从零持续到一秒，并以1/\text{fs}的时间间隔进行采样。
参数f_1和f_2分别代表两个信号的频率值分别为50赫兹与120赫兹。
合成信号x由两个基本正弦波叠加而成。
合成信号y同样由两个基本正弦波叠加而成但其相位存在细微差别。

计算互谱 ：

 * `csd(x, y, fs)` 计算信号 x 和 y 的互谱密度，返回频率向量 `f` 和互谱密度 `Pxy`。

绘制互谱图 ：

通过 matplotlib 绘制互谱图的可视化图形，并设定其显示范围与坐标标注。具体操作包括：设定 x 坐标轴的取值范围为 [0, fs/2] 以限定频谱分析仅显示正频率部分；并在 x 轴和 y 轴上分别标注 '频率 (Hz)' 和 '互谱密度 (dB/Hz)' 来命名坐标轴内容。

7. 现代谱估计方法

现代谱估计技术在提升信号精确度与分辨率达到一定程度上已取得显著成效，在应对非平稳信号以及多信号混杂的情况时展现出明显优势。这类方法主要包含基于子空间模型以及基于最大似然估计的技术

7.1 基于子空间的方法

通过子空间方法将信号空间划分为信号子空间与噪声子空间从而计算信号的频谱。主要采用MUSIC算法与ESPRIT算法这两种方法在处理多输入信号以及实现高分辨率谱估计方面表现优异

7.2 基于最大似然的方法

采用最大似然法利用最大化似然函数用于估计信号参数。此类方法特别适用于处理复杂的信号模型，并且表现出色。然而其计算复杂度相对较高。其中最常见的两种方法分别是Capon法以及MEMP（Maximum Entropy Method）法。

7.3 Python 代码示例

下面是一个使用 Python 和 SciPy 库进行 Capon 方法谱估计的示例：

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy import linalg
    
    # 生成一个示例信号
    fs = 1000  # 采样频率
    t = np.arange(0, 1, 1/fs)  # 时间向量
    f1 = 50  # 信号频率1
    f2 = 120  # 信号频率2
    x = 0.7 * np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)  # 合成信号
    
    # 计算协方差矩阵
    Rxx = np.cov(x)
    
    # 特征值分解
    eigenvalues, eigenvectors = linalg.eig(Rxx)
    
    # 选择噪声子空间
    noise_eigenvectors = eigenvectors[:, 2:]
    
    # 定义频率搜索范围
    f_search = np.linspace(0, fs/2, 500)
    
    # 计算 Capon 谱
    P_capon = np.zeros_like(f_search)
    for f in f_search:
    e = np.exp(-1j * 2 * np.pi * f * t)
    Rxx_inv = linalg.inv(Rxx)
    P_capon[f_search == f] = 1 / np.dot(e.T, np.dot(Rxx_inv, e)).real
    
    # 绘制 Capon 谱估计
    plt.figure(figsize=(10, 6))
    plt.plot(f_search, 10 * np.log10(P_capon))
    plt.title('Capon 方法谱估计')
    plt.xlabel('频率 (Hz)')
    plt.ylabel('功率谱密度 (dB/Hz)')
    plt.grid(True)
    plt.show()

7.4 代码解释

生成示例信号 ：

• fs 被设定为采样频率，并指定值为 1000 Hz。

  • t 被定义为空间范围从 0 到 1 秒的时间向量。
  • f1 和 f2 被指定了分别为两个信号的频率参数值：50 Hz 和 120 Hz。
  • x 被描述为空间范围内的合成信号变量名，并指出其由两个正弦波叠加而成。

计算协方差矩阵 ：

 * `np.cov(x)` 计算信号 x 的协方差矩阵 `Rxx`。

特征值分解 ：

linalg.eig(Rxx) 被用来对 covariance matrix 进行 eigenvalue decomposition, 输出 eigenvalues 和 eigenvectors.

选择噪声子空间 ：

 * `noise_eigenvectors` 是噪声子空间的特征向量。

定义频率搜索范围 ：

 * `f_search` 是频率搜索范围，从 0 到 fs/2，共 500 个点。

计算 Capon 谱 ：

 * `P_capon` 是 Capon 谱的估计结果。
 * 对于每个频率点 `f`，计算对应的谱值。

绘制 Capon 谱估计 ：

通过 matplotlib 绘制 Capon 谱估计图。通过使用 plt.xlabel('频率 (Hz)') 和 plt.ylabel('功率谱密度 (dB/Hz)') 设置 x 轴和 y 轴的标签

8. 谱估计的应用

谱估计技术在多个领域中都有广泛的应用，以下是一些典型的应用场景：

8.1 通信系统

在通信领域中，谱估计被用来分析信道属性、识别信号以及检测干扰。例如，在实际应用中，在估计了信道的频率响应后，则可以通过实现信道均衡以及多用户检测来提高系统的性能。

8.2 雷达系统

在雷达系统中，在处理接收信号时会用到一种叫做谱估计的技术来完成目标检测和参数估计的任务。利用接收到的信号频谱的信息来进行分析研究后能够推断出目标的距离、速度以及运动方向。

8.3 生物医学信号处理

在生物医学信号处理领域中，谱估计被用来分析心电图（ECG）、脑电图（EEG）等类型的信号；基于频谱分析方法能够从信号中提取出其特征信息，并以实现对这些生理过程的诊断和监控目的为基础进行应用研究

8.4 音频处理

在音频处理领域内使用谱估计来分析音频信号的频谱特征。具体来说，在进行音频压缩、减少噪音干扰以及音乐信号处理时都会用到这种方法。通过频谱分析技术，在识别并去除噪音的同时能够有效提升声音质量。

8.5 机械故障诊断

在机械故障诊断领域中，谱估计被用作分析振动信号的关键工具。它不仅有助于识别设备运行中的特征频率及其振幅分布情况，在定位潜在的问题方面也表现出显著的优势。频谱分析能够帮助工程师准确判断设备的工作状态并预测可能出现的问题。通过这种技术手段，在早期发现问题的基础上采取必要的维护措施以防止设备进一步损坏是可行的

9. 结论

作为数字信号处理中的核心工具之一，在众多工程应用中发挥着关键作用的各种手段能够有效解析与信号频谱相关的各种特性；涵盖从经典傅里叶变换到现代高分辨率谱估计等各类方法；在电子工程、通信技术等多个领域均展现了显著的应用价值；具体应用场景以及信号特征等因素共同决定了哪种谱估计方法更为合适

本文旨在帮助读者全面掌握谱估计技术并将其应用于实际问题中