第 11 章 功率谱估计 / Power Spectrum Estimation

9.1 功率谱密度的概念

随机信号（Random Signal）不具有确定性的傅里叶变换，但其功率在频域上的分布——功率谱密度（Power Spectral Density, PSD）——是描述其频域特性的核心工具。

9.1.1 定义

对于宽平稳（Wide-Sense Stationary, WSS）随机过程 $x[n]$，功率谱密度定义为：

$$P_{xx}(\omega) = \lim_{N\to\infty} \frac{1}{N} E\left[|X_N(\omega)|^2\right]$$

其中 $X_N(\omega)$ 是长度为 $N$ 的截断信号的 DTFT：

$$X_N(\omega) = \sum_{n=0}^{N-1} x[n] e^{-j\omega n}$$

根据 Wiener-Khinchin 定理，PSD 等于自相关函数（Autocorrelation Function）的 DTFT：

$$P_{xx}(\omega) = \sum_{k=-\infty}^{\infty} r_{xx}[k] e^{-j\omega k}$$

其中 $r_{xx}[k] = E[x[n]x[n+k]]$。

9.1.2 估计的基本问题

实际中只能获得有限长观测数据，因此需要从 $N$ 个样本中估计 $P_{xx}(\omega)$。评价谱估计质量的三个核心指标：

这三个指标之间存在折中关系（Trade-off）：降低方差往往牺牲分辨率，反之亦然。

9.2 非参数方法（Non-parametric Methods）

非参数方法不假设信号模型，直接从数据估计 PSD。核心思想是利用 FFT 高效计算。

9.2.1 周期图法（Periodogram）

最直接的 PSD 估计方法：

$$\hat{P}{per}(\omega) = \frac{1}{N} |X(\omega)|^2 = \frac{1}{N} \left| \sum{n=0}^{N-1} x[n] e^{-j\omega n} \right|^2$$

实际计算中，在 $N$ 个等间距频率点 $\omega_k = 2\pi k / N$ 上用 FFT 实现。

1
2
3
4
5
6

# 周期图法伪代码
def periodogram(x, N=None):
    N = N or len(x)
    X = fft(x, N)              # N 点 FFT
    P = (1.0 / len(x)) * abs(X)**2
    return P[:N//2+1]          # 返回单边谱

特点：

• 优点：计算简单，一次 FFT 搞定
• 缺点：方差约为 $P_{xx}^2(\omega)$（与数据长度 $N$ 无关！），一致性差
• 渐近特性：$\hat{P}{per}(\omega)$ 是 $P{xx}(\omega)$ 的渐近无偏估计，但非一致估计

9.2.2 Bartlett 法（分段平均周期图）

将长度为 $N$ 的数据分成 $K$ 段（每段 $M = N/K$ 个样本），分别计算周期图后取平均：

$$\hat{P}{Bart}(\omega) = \frac{1}{K} \sum{i=1}^{K} \hat{P}_{per}^{(i)}(\omega)$$

其中第 $i$ 段的周期图为：

$$\hat{P}{per}^{(i)}(\omega) = \frac{1}{M} \left| \sum{n=0}^{M-1} x[n + (i-1)M] e^{-j\omega n} \right|^2$$

方差降低为 $1/K$，但频率分辨率从 $\Delta f = 2\pi/N$ 恶化为 $\Delta f = 2\pi/M$。

1
2
3
4
5
6
7
8

# Bartlett 法伪代码
def bartlett(x, K):
    M = len(x) // K
    P_sum = zeros(M // 2 + 1)
    for i in range(K):
        segment = x[i*M : (i+1)*M]
        P_sum += periodogram(segment, M)
    return P_sum / K

9.2.3 Welch 法（重叠分段 + 加窗）

工程中最常用的非参数谱估计方法。在 Bartlett 法基础上引入两个改进：

1. 相邻段重叠（Overlap）：通常重叠 50%
2. 加窗（Windowing）：每段乘以窗函数（如汉宁窗）减少频谱泄漏

$$\hat{P}{Welch}(\omega) = \frac{1}{K} \sum{i=1}^{K} \hat{P}_{win}^{(i)}(\omega)$$

其中加窗后的周期图为：

$$\hat{P}{win}^{(i)}(\omega) = \frac{1}{U \cdot M} \left| \sum{n=0}^{M-1} w[n] \cdot x[n + iD] \cdot e^{-j\omega n} \right|^2$$

$D$ 为段间偏移量（$D = M/2$ 时重叠 50%），$U = \frac{1}{M}\sum_{n=0}^{M-1} w^2[n]$ 为窗函数归一化因子。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

# Welch 法伪代码
def welch(x, M, overlap=0.5, window='hann'):
    w = get_window(window, M)        # 生成窗函数
    U = sum(w**2) / M                 # 归一化因子
    D = int(M * (1 - overlap))        # 段间偏移
    K = (len(x) - M) // D + 1         # 段数

    P_sum = zeros(M // 2 + 1)
    for i in range(K):
        start = i * D
        segment = x[start:start+M] * w
        P_sum += abs(fft(segment, M)[:M//2+1])**2
    return P_sum / (K * U * M)

flowchart TD
    A[原始信号 x n 长度 N] --> B[分段
每段 M 点 重叠 50%]
    B --> C[每段加窗
Hann / Hamming]
    C --> D[每段 FFT]
    D --> E[计算模平方]
    E --> F[多段平均]
    F --> G[归一化输出
Welch PSD 估计]

Welch 法参数选择经验：

9.2.4 Blackman-Tukey 法（加窗自相关法）

先估计自相关函数，再加窗平滑后做 FFT：

$$\hat{P}{BT}(\omega) = \sum{k=-(M-1)}^{M-1} w[k] \cdot \hat{r}_{xx}[k] \cdot e^{-j\omega k}$$

其中 $\hat{r}_{xx}[k]$ 为样本自相关估计，$w[k]$ 为滞后窗（Lag Window，如三角窗或汉明窗），$M \ll N$。

原理：自相关估计的尾部（大滞后）方差大、可靠性低，加窗截断抑制了这些不可靠部分，等效于在频域做平滑，降低方差。

9.3 参数方法（Parametric Methods）

参数方法假设信号服从某种模型，通过估计模型参数间接获得 PSD。优势在于：短数据下仍可获得高分辨率。

9.3.1 AR 模型（AutoRegressive Model）

$p$ 阶 AR 模型（也称线性预测模型）：

$$x[n] = -\sum_{k=1}^{p} a_k , x[n-k] + w[n]$$

其中 $w[n]$ 是方差为 $\sigma_w^2$ 的白噪声。对应的 PSD 为：

$$P_{AR}(\omega) = \frac{\sigma_w^2}{\left| 1 + \sum_{k=1}^{p} a_k , e^{-j\omega k} \right|^2}$$

可见，只需估计 $p+1$ 个参数 ${a_1, \ldots, a_p, \sigma_w^2}$ 即可得到全频段的 PSD。

9.3.2 Yule-Walker 方程

对 AR 模型两边乘以 $x[n-m]$ 并取期望，得到 Yule-Walker 方程组：

$$\begin{bmatrix} r[0] & r[1] & \cdots & r[p-1] \ r[1] & r[0] & \cdots & r[p-2] \ \vdots & \vdots & \ddots & \vdots \ r[p-1] & r[p-2] & \cdots & r[0] \end{bmatrix} \begin{bmatrix} a_1 \ a_2 \ \vdots \ a_p \end{bmatrix} = -\begin{bmatrix} r[1] \ r[2] \ \vdots \ r[p] \end{bmatrix}$$

噪声方差：

$$\sigma_w^2 = r[0] + \sum_{k=1}^{p} a_k , r[k]$$

系数矩阵是 Toeplitz 结构，可用 Levinson-Durbin 递推 在 $O(p^2)$ 时间内高效求解。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

# Levinson-Durbin 递推伪代码
def levinson_durbin(r, p):
    a = [1.0]  # a_0 = 1
    sigma2 = r[0]
    for k in range(1, p + 1):
        # 反射系数
        rho = -sum(a[j] * r[k-j] for j in range(k)) / sigma2
        # 更新系数
        a_new = a + [0]
        for j in range(1, k):
            a_new[j] = a[j] + rho * a[k - j]
        a_new[k] = rho
        a = a_new
        sigma2 *= (1 - rho**2)
    return a[1:], sigma2  # 返回 AR 系数和噪声方差

9.3.3 Burg 算法

Burg 算法直接从数据估计反射系数，不经过自相关估计，因而对短数据更稳健。

核心思想：最小化前向和后向预测误差的平均功率：

$$\epsilon^{(p)} = \frac{1}{2(N-p)} \sum_{n=p}^{N-1} \left( |e_f^{(p)}[n]|^2 + |e_b^{(p)}[n]|^2 \right)$$

第 $p$ 阶反射系数的闭合解：

$$k_p = \frac{-2 \sum_{n=p}^{N-1} e_f^{(p-1)}[n] \cdot e_b^{(p-1)}[n-1]}{\sum_{n=p}^{N-1} \left( |e_f^{(p-1)}[n]|^2 + |e_b^{(p-1)}[n-1]|^2 \right)}$$

预测误差的递推更新：

$$e_f^{(p)}[n] = e_f^{(p-1)}[n] + k_p \cdot e_b^{(p-1)}[n-1]$$ $$e_b^{(p)}[n] = e_b^{(p-1)}[n-1] + k_p^* \cdot e_f^{(p-1)}[n]$$

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

# Burg 算法伪代码
def burg(x, p):
    N = len(x)
    ef = x.copy()  # 前向误差，初始为数据本身
    eb = x.copy()  # 后向误差
    a = []          # AR 系数
    sigma2 = sum(x**2) / N  # 初始误差功率

    for k in range(p):
        # 反射系数
        num = -2 * sum(ef[n] * eb[n-1] for n in range(k+1, N))
        den = sum(ef[n]**2 + eb[n-1]**2 for n in range(k+1, N))
        kappa = num / den

        # 更新 AR 系数（Levinson 递推）
        a_new = a + [0]
        for j in range(k):
            a_new[j] = a[j] + kappa * a[k-1-j]
        a_new[k] = kappa
        a = a_new

        # 更新误差
        ef_new = ef.copy()
        for n in range(k+1, N):
            ef_new[n] = ef[n] + kappa * eb[n-1]
            eb[n] = eb[n-1] + kappa * ef[n]  # 用旧 ef
        ef = ef_new
        sigma2 *= (1 - kappa**2)

    return a, sigma2

flowchart LR
    A[观测数据 x n] --> B[Burg 算法
逐阶递推]
    B --> C[反射系数 k_1 k_p]
    B --> D[AR 系数 a_1 a_p]
    B --> E[噪声方差 sigma_w²]
    D --> F[代入 AR PSD 公式]
    E --> F
    F --> G[高分辨率 PSD 估计]

9.3.4 AR vs MA vs ARMA

工程实践：AR 模型因其参数估计简单（Yule-Walker / Burg 均有高效算法）且对尖锐谱峰建模能力突出，在实际中应用最广泛。ARMA 模型虽更灵活，但参数估计是非线性问题，计算复杂度高。

AR 模型阶数选择：常用准则——

• AIC（Akaike Information Criterion）：$\text{AIC}(p) = N \ln \hat{\sigma}_w^2 + 2p$
• MDL（Minimum Description Length）：$\text{MDL}(p) = N \ln \hat{\sigma}_w^2 + p \ln N$

选使准则最小的 $p$ 值。MDL 更倾向于选择较低阶数，一致性更好。

9.4 性能对比

关键结论：

1. 非参数方法（Welch 为主）：数据充足时首选，结果稳健、无模型失配风险
2. 参数方法（AR/Burg 为主）：短数据、需要高分辨率时使用，但存在模型阶数选择和模型失配问题
3. Welch 法是工程中最通用的非参数方法，几乎是频谱分析仪的标配算法

9.5 工程应用

9.5.1 频谱监测（Spectrum Monitoring）

在无线电频谱管理中，需要实时监测频段占用情况。Welch 法因其稳健性和 FFT 的高效性，是认知无线电（Cognitive Radio）中频谱感知（Spectrum Sensing）的基础算法。

典型参数：采样率数十 MHz，FFT 点数 1024~4096，重叠 50%，汉宁窗。能量检测门限基于噪声功率谱估计设定。

9.5.2 振动分析（Vibration Analysis）

旋转机械故障诊断（如轴承缺陷、齿轮磨损）通过振动信号的功率谱特征识别故障类型。

• 轴承外圈缺陷：在特征频率 $f_{BPFO}$ 处出现谱峰及谐波
• 齿轮啮合故障：啮合频率 $f_{mesh}$ 及其边带增强
• 不平衡 / 不对中：1× 和 2× 转频处出现明显峰值

AR 谱估计在短数据段（如启停机瞬态过程）中优势明显，可追踪频率随转速的变化（阶次跟踪）。

9.5.3 生物信号处理（Biomedical Signal Processing）

脑电图（EEG） 功率谱分析是临床和研究中不可或缺的工具：

心电图（ECG）心率变异性（HRV）分析中，PSD 的低频（LF, 0.04–0.15 Hz）与高频（HF, 0.15–0.4 Hz）功率比值 LF/HF 反映交感-副交感神经平衡状态。Welch 法是 HRV 频域分析的标准方法。

flowchart TD
    A[原始信号采集] --> B{信号特征?}
    B -->|平稳 长数据| C[Welch 法
稳健 通用]
    B -->|平稳 短数据| D[AR Burg 法
高分辨率]
    B -->|非平稳| E[时频分析
STFT / 小波]
    C --> F[频谱特征提取]
    D --> F
    F --> G{应用场景}
    G --> H[频谱监测
占用检测]
    G --> I[振动分析
故障诊断]
    G --> J[生物信号
EEG / HRV]

9.6 本章小结

功率谱估计是随机信号频域分析的核心技术。非参数方法（尤其是 Welch 法）以其稳健性和高效性成为工程首选；参数方法（尤其是 AR 模型的 Burg 算法）在短数据、高分辨率需求场景中发挥重要作用。实际工程中，应根据数据长度、信号特性、分辨率需求和计算资源综合选择合适的方法。

进一步阅读：

• Stoica, P. & Moses, R., Spectral Analysis of Signals, Pearson, 2005
• Oppenheim, A.V. & Schafer, R.W., Discrete-Time Signal Processing, 3rd ed., Pearson, 2010, Ch. 10
• Hayes, M.H., Statistical Digital Signal Processing and Modeling, Wiley, 1996