4. 数字信号处理基础：DFT、FFT、数字滤波器与窗函数

各位同学，咱们今天聊点硬核的。数字信号处理，说白了就是让单片机听懂心电信号在说什么。我当年刚入行时，觉得这玩意儿就是数学公式堆砌，直到在除颤仪项目里被50Hz工频干扰折磨得欲仙欲死，才明白——不懂DSP，你连噪声都干不掉。

4.1 离散傅里叶变换（DFT）——从时域到频域的“翻译官”

心电信号是随时间变化的电压值，这叫时域。但很多特征在时域里看不清楚，比如肌电干扰、基线漂移。这时候就需要DFT，把它翻译到频域。

DFT的公式长这样：

X[k] = Σ(n=0 to N-1) x[n] · e^(-j·2π·k·n/N)

别被吓到。说白了就是：把一段信号跟不同频率的正余弦波做内积，看看哪个频率成分多。

我在项目中遇到过一个坑：采样率Fs和点数N选不好，频率分辨率就崩了。频率分辨率Δf = Fs / N。比如你采样率250Hz，取250个点，Δf就是1Hz。想分辨0.5Hz的呼吸波？不行，得加长采样时间。

注意：DFT计算量是O(N²)。N=1024时，要算100多万次复数乘法。在嵌入式MCU上跑？等着看门狗复位吧。

4.2 快速傅里叶变换（FFT）——嵌入式救星

FFT不是新变换，是DFT的快速算法。它利用旋转因子的对称性和周期性，把计算量降到O(N·log₂N)。N=1024时，FFT只需要约1万次运算，比DFT快两个数量级。

最常用的是基2时间抽取FFT。要求N必须是2的整数次幂，比如256、512、1024。

// 我常用的FFT实现片段（C语言，定点数优化版）
void fft_fixed(int16_t *real, int16_t *imag, uint16_t n) {
    // 位反转排序
    for (uint16_t i = 0; i < n; i++) {
        uint16_t j = bit_reverse(i, log2(n));
        if (i < j) {
            swap(&real[i], &real[j]);
            swap(&imag[i], &imag[j]);
        }
    }
    // 蝶形运算
    for (uint16_t len = 2; len <= n; len <<= 1) {
        // ... 旋转因子计算与蝶形合并
    }
}

我个人习惯用定点数做FFT。浮点运算在Cortex-M4上虽然快，但功耗高。定点数用Q15格式，一个int16_t搞定实部和虚部，速度翻倍。

经验之谈：除颤仪的心电采样率通常是250Hz或500Hz。做FFT时取256点或512点，刚好覆盖2-40Hz的诊断频段。别取太多点，否则低频噪声全进来了。

4.3 数字滤波器设计——给信号“洗澡”

心电信号里噪声多：工频50Hz、肌电干扰、电极运动伪迹。滤波器就是给信号洗澡的。

4.3.1 FIR滤波器——稳定但费资源

FIR（有限冲激响应）滤波器，结构简单：

y[n] = Σ(k=0 to M-1) b[k] · x[n-k]

优点：线性相位，不会把心电波的形态搞歪。缺点：阶数高，M=100时，每个点要算100次乘加。

我曾经在STM32F4上实现一个50Hz陷波器，用FIR需要128阶，跑一次要0.5ms。后来换成IIR，8阶搞定，0.05ms完事。

4.3.2 IIR滤波器——高效但需小心

IIR（无限冲激响应）滤波器，带反馈：

y[n] = Σ(k=0 to P) b[k]·x[n-k] - Σ(k=1 to Q) a[k]·y[n-k]

优点：阶数低，计算快。缺点：可能不稳定，相位非线性。

在除颤仪里，IIR常用于高通滤波（去基线漂移）和陷波（去工频）。我一般用巴特沃斯或切比雪夫II型，通带平坦，阻带衰减够用。

滤波器类型 阶数 计算量 相位特性 适用场景
FIR 高（50-200） 大 线性 心电波形保真要求高
IIR 低（4-8） 小 非线性 实时性要求高、资源受限
避坑指南：IIR滤波器系数要用双精度算好，再量化到定点数。我见过有人直接用单精度浮点算系数，量化后滤波器直接振荡——输出全是大波浪，差点把除颤仪误判成室颤。

4.4 窗函数应用——别让信号“漏”了

做FFT时，我们只能取一段有限长度的信号。这相当于给信号加了个矩形窗，会在频域产生频谱泄漏——能量从主瓣漏到旁瓣，把真实频率搞模糊了。

窗函数就是用来抑制泄漏的。常见的有：

  • 汉宁窗：主瓣宽，旁瓣衰减快。我常用，适合心电频谱分析。
  • 海明窗：跟汉宁类似，但旁瓣更低一点。
  • 布莱克曼窗：旁瓣衰减最猛，但主瓣最宽，频率分辨率下降。
  • 凯泽窗：可调参数β，灵活控制主瓣宽度和旁瓣衰减。
// 汉宁窗生成代码
void hanning_window(float *w, uint16_t n) {
    for (uint16_t i = 0; i < n; i++) {
        w[i] = 0.5f * (1.0f - cosf(2.0f * M_PI * i / (n - 1)));
    }
}

你想想看：如果心电信号里有个0.5Hz的呼吸波，你用矩形窗做FFT，它可能泄漏到0.3Hz和0.7Hz，把ST段分析搞乱。加个汉宁窗，泄漏就压下去了。

我的习惯：做心电频谱分析时，先加汉宁窗再做FFT。做QRS波检测时，不加窗，直接用原始信号，因为窗函数会模糊时域特征。

4.5 嵌入式实现要点——从理论到代码

理论说完了，咱们聊聊怎么在MCU上跑起来。

  1. 定点数优化：用Q15或Q31格式，避免浮点运算。ARM Cortex-M4有SIMD指令，一次算两个16位乘加，速度翻倍。
  2. 查表法：旋转因子sin/cos提前算好存ROM，别实时算。256点FFT的旋转因子表才2KB，省下大量CPU时间。
  3. 乒乓缓冲：用双缓冲区，一个采数据，一个做FFT。DMA填满一个缓冲区后触发中断，CPU处理另一个。这样采样和处理并行，不丢数据。
  4. 滤波器系数量化：IIR系数用Q15格式，注意防止溢出。我一般把系数归一化到±0.999，留点余量。
// 乒乓缓冲示例（伪代码）
uint16_t buf_a[256], buf_b[256];
volatile uint8_t active_buf = 0;

void DMA_IRQHandler(void) {
    if (active_buf == 0) {
        process_fft(buf_a);  // 处理A缓冲区
        DMA_set_buffer(buf_b); // DMA填B缓冲区
    } else {
        process_fft(buf_b);
        DMA_set_buffer(buf_a);
    }
    active_buf ^= 1;
}

嗯，这里要注意：FFT处理时间不能超过采样周期。250Hz采样，周期4ms。你的FFT必须在4ms内跑完，否则下一包数据来了，你还在算上一包，系统就崩了。

警告：医疗设备里，FFT和滤波器不能阻塞主循环。一定要用RTOS任务或中断驱动。我见过有人把FFT放在主循环里轮询，结果心电波形显示卡成PPT，除颤同步延迟——这要是真救人，后果不堪设想。

好了，数字信号处理基础就聊到这儿。下一章咱们讲心电特征提取——怎么从一堆噪声里揪出QRS波。到时候我会分享一个我踩过的坑：自适应阈值怎么调才不会被T波误触发。

← 上一章 📖 返回目录 下一章 →