第 12 章 DSP 工程应用 / DSP Applications

前九章建立了 DSP 的数学基础与算法体系。本章将视角转向工程现场——看看这些理论如何在音频、通信、生物医学、雷达与图像处理等领域落地。

10.1 音频信号处理

10.1.1 感知音频编码（Perceptual Audio Coding）

MP3、AAC、Opus 等格式的核心思想：利用人耳的掩蔽效应（Masking Effect），丢弃感知上不可闻的频率分量。

子带编码（Subband Coding）将音频划分为多个频带，对各子带独立量化：

$$ X[k] = \sum_{n=0}^{N-1} x[n] , w[n] , e^{-j2\pi kn/N} $$

其中 $w[n]$ 为窗函数（通常取 MDCT 的正弦窗）。每子带分配的比特数由掩蔽阈值决定：

$$ b_i = \max!\Big(0,;\Big\lfloor \log_2!\Big(1 + \frac{\sigma_i^2}{T_i^2}\Big)\Big\rceil\Big) $$

$\sigma_i^2$ 为第 $i$ 子带信号方差，$T_i$ 为该子带的掩蔽阈值。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# 伪代码：简化感知编码比特分配
def perceptual_bit_allocation(subband_powers, masking_thresholds, total_bits):
    smr = 10 * np.log10(subband_powers / masking_thresholds)  # Signal-to-Mask Ratio
    bits = np.zeros(len(smr))
    remaining = total_bits
    while remaining > 0:
        i = np.argmax(smr - bits)  # 分配给 SMR 增益最大的子带
        bits[i] += 1
        remaining -= 1
    return bits

10.1.2 数字均衡器（Digital Equalizer）

参数均衡器（Parametric EQ）基于二阶 IIR 滤波器（Biquad）：

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}} $$

系数 $b_0, b_1, b_2, a_1, a_2$ 由中心频率 $f_0$、增益 $G$（dB）和品质因数 $Q$ 计算（Robert Bristow-Johnson 公式）。图形均衡器则将多个 Biquad 级联。

10.1.3 降噪：谱减法与维纳滤波

谱减法（Spectral Subtraction） 假设噪声功率谱在无语音段可估计：

$$ |\hat{S}(k)|^2 = \max!\big(|Y(k)|^2 - \alpha \cdot |\hat{N}(k)|^2,; \beta \cdot |\hat{N}(k)|^2\big) $$

$\alpha$ 为过减因子，$\beta$ 为谱下限因子（抑制音乐噪声 / Musical Noise）。

维纳滤波（Wiener Filter） 在最小均方误差意义下最优：

$$ H_{\text{opt}}(k) = \frac{P_S(k)}{P_S(k) + P_N(k)} = \frac{\text{SNR}(k)}{1 + \text{SNR}(k)} $$

1
2
3
4

# 维纳滤波伪代码
def wiener_filter(noisy_stft, noise_psd, signal_psd):
    gain = signal_psd / (signal_psd + noise_psd + 1e-10)
    return noisy_stft * np.sqrt(gain)  # 幅度域
flowchart LR
    A[带噪语音] --> B[STFT 分析]
    B --> C[噪声估计
VAD]
    C --> D[计算增益函数]
    D --> E[谱修正]
    E --> F[ISTFT 合成]
    F --> G[增强语音]
    style A fill:#1565C0,color:#fff
    style D fill:#2E7D32,color:#fff
    style G fill:#C62828,color:#fff

10.2 通信系统中的 DSP

10.2.1 调制解调与脉冲成形

数字调制将比特映射到复数符号（如 QPSK、16-QAM），再经脉冲成形滤波器（Pulse Shaping Filter）限制带宽：

$$ s(t) = \sum_k a_k , p(t - kT_s) $$

$p(t)$ 通常取根升余弦（Root Raised Cosine, RRC）脉冲，滚降因子 $\alpha \in [0, 1]$。发射端与接收端各用一个 $\sqrt{RRC}$，级联后满足 Nyquist 无 ISI 准则。

10.2.2 信道均衡（Channel Equalization）

多径信道引入码间干扰（ISI）。线性均衡器设计为信道频率响应的逆：

$$ H_{\text{EQ}}(f) \approx \frac{1}{H_{\text{ch}}(f)} $$

实际中常用最小均方（LMS）自适应均衡器：

$$ \mathbf{w}(n+1) = \mathbf{w}(n) + \mu \cdot e^*(n) \cdot \mathbf{x}(n) $$

其中 $e(n) = d(n) - \mathbf{w}^H(n)\mathbf{x}(n)$ 为误差信号，$\mu$ 为步长。

判决反馈均衡器（DFE） 利用已判决符号消除后向 ISI，性能优于线性均衡器但存在差错传播风险。

10.2.3 OFDM 中的 FFT/IFFT

正交频分复用（OFDM）是 4G/5G/Wi-Fi 的核心多载波技术。关键洞察：IDFT 将频域数据符号映射为时域 OFDM 符号。

$$ x[n] = \frac{1}{N}\sum_{k=0}^{N-1} X[k] , e^{j2\pi kn/N}, \quad n = 0, 1, \ldots, N-1 $$

发射端用 IFFT（$O(N\log N)$），接收端用 FFT 解调。循环前缀（Cyclic Prefix, CP）将线性卷积转化为循环卷积，使单抽头频域均衡成为可能：

$$ Y[k] = H[k] \cdot X[k] + W[k] \quad \Rightarrow \quad \hat{X}[k] = \frac{Y[k]}{H[k]} $$

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# OFDM 发射链路伪代码
def ofdm_transmit(bits, n_subcarriers, cp_len):
    symbols = qam_map(bits)                        # QAM 映射
    groups = symbols.reshape(-1, n_subcarriers)
    ofdm_symbols = []
    for group in groups:
        time_domain = np.fft.ifft(group, n=n_subcarriers)  # IFFT
        cp = time_domain[-cp_len:]                           # 循环前缀
        ofdm_symbols.append(np.concatenate([cp, time_domain]))
    return np.concatenate(ofdm_symbols)

10.3 生物医学信号处理

10.3.1 ECG QRS 复合波检测

心电图（ECG）中 QRS 复合波对应心室去极化，是心率分析的基础。经典的 Pan-Tompkins 算法流程：

flowchart TD
    A[原始 ECG] --> B[带通滤波
5-15 Hz]
    B --> C[微分]
    C --> D[平方]
    D --> E[滑动窗口积分
N=30 samples]
    E --> F[自适应阈值检测]
    F --> G{峰值 > 阈值?}
    G -->|是| H[标记 R 峰]
    G -->|否| I[更新阈值]
    H --> I
    style A fill:#1565C0,color:#fff
    style H fill:#2E7D32,color:#fff
    style F fill:#6A1B9A,color:#fff

核心步骤的数学表达：

  1. 带通滤波：级联低通（$f_c = 15$ Hz）和高通（$f_c = 5$ Hz）
  2. 微分：$y[n] = x[n] - x[n-1]$（突出 QRS 的高斜率特征）
  3. 平方：$z[n] = y^2[n]$（非线性放大）
  4. 积分：$w[n] = \frac{1}{N}\sum_{i=0}^{N-1} z[n-i]$

10.3.2 EEG 频带分析

脑电图（EEG）的诊断信息集中在特定频带：

频带频率范围生理意义
Delta（δ）0.5–4 Hz深睡眠
Theta（θ）4–8 Hz困倦、冥想
Alpha（α）8–13 Hz闭眼放松
Beta（β）13–30 Hz主动思考
Gamma（γ）30–100 Hz高级认知

提取方法：设计各频带的带通 FIR 滤波器，或对分段信号做 FFT/PSD 估计：

$$ P_{xx}(f) = \frac{1}{K}\sum_{i=1}^{K} |X_i(f)|^2 $$

其中 $K$ 为 Welch 平均周期图中的分段数。

10.4 雷达与声纳信号处理

10.4.1 脉冲压缩与匹配滤波

简单脉冲雷达的距离分辨力受限于脉宽 $\tau$。线性调频（LFM / Chirp）信号通过**脉冲压缩（Pulse Compression）**实现高分辨力与远探测距离的兼顾：

$$ s(t) = \text{rect}!\Big(\frac{t}{T}\Big) \cdot e^{j\pi \mu t^2} $$

$T$ 为脉宽，$\mu = B/T$ 为调频斜率，$B$ 为带宽。匹配滤波器输出：

$$ y(\tau) = \int s(t) \cdot s^*(t-\tau) , dt \approx T \cdot \text{sinc}(\pi B \tau) $$

压缩后的脉冲宽度约为 $1/B$，压缩比 $= T \cdot B$（时间-带宽积）。

10.4.2 多普勒处理（Doppler Processing）

目标径向速度 $v_r$ 引起的多普勒频移：

$$ f_d = \frac{2 v_r}{\lambda} $$

对 $M$ 个脉冲回波在慢时间维（Slow Time）做 FFT，即 MTI/脉冲多普勒处理，构成距离-多普勒（Range-Doppler）图：

$$ Z[r, d] = \sum_{m=0}^{M-1} z[r, m] , e^{-j2\pi d \cdot m / M} $$

1
2
3
4
5
6
7

# 距离-多普勒处理伪代码
def range_doppler_map(raw_pulses, n_range_bins, n_pulses):
    # 快时间维 FFT → 距离
    range_profile = np.fft.fft(raw_pulses, axis=1)
    # 慢时间维 FFT → 多普勒
    rd_map = np.fft.fftshift(np.fft.fft(range_profile, axis=0), axes=0)
    return np.abs(rd_map)

10.5 图像处理中的 DSP

10.5.1 二维线性滤波

图像可视为二维离散信号 $x[m, n]$。二维卷积：

$$ y[m, n] = \sum_{k}\sum_{l} h[k, l] \cdot x[m-k, n-l] $$

常用 $3\times3$ 核示例：

  • 锐化：$h = \begin{bmatrix} 0 & -1 & 0 \ -1 & 5 & -1 \ 0 & -1 & 0 \end{bmatrix}$
  • 高斯模糊：$h = \frac{1}{16}\begin{bmatrix} 1 & 2 & 1 \ 2 & 4 & 2 \ 1 & 2 & 1 \end{bmatrix}$
  • Sobel 边缘检测：$G_x = \begin{bmatrix} -1 & 0 & 1 \ -2 & 0 & 2 \ -1 & 0 & 1 \end{bmatrix}$

10.5.2 二维频域分析

二维 DFT：

$$ X[u, v] = \sum_{m=0}^{M-1}\sum_{n=0}^{N-1} x[m, n] , e^{-j2\pi(um/M + vn/N)} $$

频域滤波的工程优势：大核空域卷积（如 $31\times31$ 高斯）可通过 FFT 的 $O(N^2\log N)$ 实现，远快于空域的 $O(N^2 K^2)$。

频域滤波流程：

$$ Y[u,v] = X[u,v] \cdot H[u,v] \quad \Rightarrow \quad y[m,n] = \text{IDFT}{Y[u,v]} $$

10.6 工具链与开发实践

10.6.1 Python 生态

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

# SciPy 滤波器设计示例
from scipy.signal import butter, filtfilt, welch
import numpy as np

# 设计 4 阶 Butterworth 带通滤波器 (5-15 Hz, fs=360 Hz)
b, a = butter(4, [5, 15], btype='band', fs=360)
ecg_filtered = filtfilt(b, a, ecg_raw)

# Welch PSD 估计
f, Pxx = welch(eeg_segment, fs=256, nperseg=512)
alpha_power = np.trapz(Pxx[(f >= 8) & (f <= 13)], f[(f >= 8) & (f <= 13)])

常用库：numpy（数组运算）、scipy.signal（滤波/谱分析）、librosa（音频）、pywt（小波变换）。

10.6.2 MATLAB

1
2
3
4

% OFDM 仿真片段
nSubcarriers = 64;
dataSymbols = qammod(randi([0 15], nSubcarriers, 1), 16, 'UnitAveragePower', true);
txSignal = ifft(dataSymbols, nSubcarriers);

MATLAB 的 Signal Processing Toolbox、Communications Toolbox、Phased Array System Toolbox 提供从算法验证到系统仿真的完整工具链。

10.6.3 嵌入式 DSP：ARM CMSIS-DSP

从仿真到硬件部署是 DSP 工程的关键一步。ARM Cortex-M 系列通过 CMSIS-DSP 库提供优化的信号处理函数：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

#include "arm_math.h"

// 1024 点实数 FFT
arm_rfft_fast_instance_f32 fft_inst;
arm_rfft_fast_init_f32(&fft_inst, 1024);
arm_rfft_fast_f32(&fft_inst, input, output, 0);  // 正变换

// FIR 滤波
arm_fir_instance_f32 fir_inst;
arm_fir_init_f32(&fir_inst, numTaps, coeffs, state, blockSize);
arm_fir_f32(&fir_inst, input, output, blockSize);

CMSIS-DSP 利用 Cortex-M 的 SIMD 指令（如 Helium / DSP 扩展）实现关键运算的向量化，单周期 MAC 操作使实时滤波成为可能。

10.6.4 开发流程总览

flowchart LR
    A[算法设计
MATLAB/Python] --> B[定点仿真
量化与溢出分析]
    B --> C[C/C++ 实现
CMSIS-DSP]
    C --> D[交叉编译
目标平台]
    D --> E[实时验证
示波器/逻辑分析仪]
    E --> F{性能达标?}
    F -->|否| B
    F -->|是| G[产品部署]
    style A fill:#1565C0,color:#fff
    style G fill:#2E7D32,color:#fff
    style F fill:#E65100,color:#fff

本章小结

应用领域核心 DSP 技术关键数学工具
音频编码感知编码、MDCT心理声学模型、比特分配
音频降噪谱减法、维纳滤波STFT、功率谱估计
通信（OFDM）IFFT/FFT、均衡IDFT、自适应滤波（LMS）
ECG 检测带通滤波、阈值检测微分-平方-积分
EEG 分析频带滤波、PSD 估计Welch 周期图
雷达匹配滤波、多普勒 FFT二维 FFT、模糊函数
图像处理二维卷积、频域滤波二维 DFT

DSP 的工程价值在于：用统一的数学框架（卷积、变换、估计）解决跨域的信号处理问题。从算法设计到嵌入式部署，工具链的成熟使得这一过程日益标准化。掌握本章技术要点的关键不在于记忆每个公式，而在于理解"时域-频域对偶"、“时宽带宽积”、“自适应优化"等贯穿各领域的共性思想。