STM32实战：零公式实现巴特沃斯低通滤波的工程化方案

在嵌入式开发中，信号处理往往被视为需要深厚数学功底的领域。当项目需要处理肌电、心电或加速度计等传感器信号时，许多开发者面对差分方程和Z变换公式望而却步。但事实上，现代工具链已经让数字滤波的实现变得像组装乐高积木一样简单——你只需要知道如何获取正确的"积木块"（滤波器系数）和如何拼接它们（C语言实现）。

1. 滤波器设计的捷径：MATLAB工具链实战

传统数字信号处理教材会从拉普拉斯变换开始，逐步推导到差分方程，这种理论路径对工程实践者并不友好。实际上，利用MATLAB的Filter Designer工具，我们可以完全跳过数学推导，直接获得可用的滤波器系数。

1.1 可视化设计三步法

打开MATLAB，在命令行输入：

filterDesigner

你会看到如下关键配置区域：

1. 响应类型 ：选择Lowpass（低通）
2. 设计方法 ：选择IIR → Butterworth
3. 频率参数 ：
   • Fs（采样频率）：根据你的ADC采样率设置（如1000Hz）
   • Fc（截止频率）：根据信号特性设置（如20Hz）

点击"Design Filter"后，工具会自动计算满足指标的 最小阶数 。对于STM32F4系列，建议从4阶开始测试，资源受限时可降至2阶。

1.2 系数导出技巧

在Filter Designer界面：

1. 点击"File" → "Export..."
2. 选择输出为"Second-Order Sections"（SOS）
3. 勾选"Optimize for code generation"

这会生成更适合嵌入式实现的二阶节系数。例如一个4阶巴特沃斯滤波器（20Hz截止，1000Hz采样）的典型系数：

Section	b0	b1	b2	a0	a1	a2
1	0.0039	0.0078	0.0039	1.0000	-1.815	0.8314
2	1.0000	2.0000	1.0000	1.0000	-1.793	0.8170

提示：a0总是1，实际存储时可以省略，减少Flash占用

2. STM32上的C语言实现艺术

获得系数后，我们需要将其转化为高效的C代码。不同于MATLAB的向量运算，嵌入式实现需特别注意实时性和内存管理。

2.1 二阶节直接II型实现

这是最节省内存的实现方式，每个二阶节只需要2个历史状态变量：

typedef struct {
    float b[3];  // 分子系数b0,b1,b2
    float a[2];  // 分母系数a1,a2 (a0=1)
    float w[2];  // 状态变量
} BiquadSection;

float biquadProcess(BiquadSection* section, float input) {
    float wn = input - section->a[0]*section->w[0] - section->a[1]*section->w[1];
    float output = section->b[0]*wn + section->b[1]*section->w[0] + section->b[2]*section->w[1];
    
    section->w[1] = section->w[0];
    section->w[0] = wn;
    
    return output;
}

2.2 全滤波器链式调用

对于4阶滤波器（2个二阶节），创建滤波器实例并级联：

BiquadSection filter[2] = {
    {{0.0039, 0.0078, 0.0039}, {-1.815, 0.8314}, {0,0}}, // 第1节
    {{1.0000, 2.0000, 1.0000}, {-1.793, 0.8170}, {0,0}}  // 第2节
};

float butterworthFilter(float input) {
    float stage1 = biquadProcess(&filter[0], input);
    return biquadProcess(&filter[1], stage1);
}

2.3 定点数优化技巧

对于M0/M3等无FPU的MCU，可采用Q格式定点数优化。例如Q15格式（16位有符号数，小数点在第15位后）：

typedef int16_t q15_t;

q15_t q15_biquadProcess(q15_t b[3], q15_t a[2], q15_t w[2], q15_t input) {
    int32_t wn = (int32_t)input - 
                ((int32_t)a[0]*w[0] >> 15) - 
                ((int32_t)a[1]*w[1] >> 15);
    
    int32_t out = ((int32_t)b[0]*wn >> 15) +
                 ((int32_t)b[1]*w[0] >> 15) + 
                 ((int32_t)b[2]*w[1] >> 15);
    
    w[1] = w[0];
    w[0] = (q15_t)wn;
    
    return (q15_t)out;
}

注意：系数导入前需乘以32767并取整，动态范围需控制在±1.0之间

3. 性能与效果的工程权衡

在实际项目中，滤波器的阶数和实现方式需要根据硬件资源灵活调整。以下是基于STM32F407（168MHz）的实测数据：

阶数	浮点运算时间(μs)	Q15运算时间(μs)	RAM占用(字节)
2	8.2	12.7	16
4	16.5	25.3	32
6	24.8	38.0	48

3.1 采样率与截止频率的黄金比例

经验表明，当采样频率(Fs)与截止频率(Fc)的比值超过10:1时，滤波效果显著提升。但在资源受限时，可接受的最低比例为：

Fs/Fc ≥ 5 （保证基本滤波效果）
Fs/Fc ≥ 10 （获得理想衰减特性）

3.2 阶数选择策略

根据信号特性推荐：

• 生物电信号（心电/肌电） ：4-6阶，50Hz工频干扰需要至少40dB衰减
• 加速度计信号 ：2-4阶，机械振动通常带宽较宽
• 音频信号 ：4-8阶，需要更陡峭的过渡带

4. 调试与优化实战指南

4.1 频率响应测试方法

无需专业设备，用代码自生成测试信号：

void freqResponseTest(float freq, uint32_t length) {
    float sampleInterval = 1.0f / SAMPLE_RATE;
    float input[length], output[length];
    
    // 生成正弦扫频信号
    for(uint32_t i=0; i<length; i++) {
        float t = i * sampleInterval;
        input[i] = sinf(2 * M_PI * freq * t);
        output[i] = butterworthFilter(input[i]);
    }
    
    // 计算幅值比（简易FFT）
    float in_max = 0, out_max = 0;
    for(uint32_t i=length/2; i<length; i++) {
        in_max = fmaxf(in_max, fabsf(input[i]));
        out_max = fmaxf(out_max, fabsf(output[i]));
    }
    
    printf("Freq: %.1fHz, Attenuation: %.2fdB\n", 
           freq, 20*log10f(out_max/in_max));
}

4.2 常见问题排查

问题1：输出信号幅度异常

• 检查系数符号：a系数通常为负值
• 确认历史状态变量初始化：首次调用前应清零

问题2：高频分量未被有效滤除

• 验证采样率与截止频率比例是否≥5:1
• 检查MATLAB设计时的Fs设置是否与代码一致

问题3：滤波器不稳定（输出震荡）

• 降低阶数尝试（如6阶→4阶）
• 改用更稳定的滤波器类型（如切比雪夫I型）

4.3 动态参数调整技巧

对于需要运行时调整截止频率的场景，可采用系数预计算表：

typedef struct {
    float cutoffFreq;
    float b[3];
    float a[2];
} FilterCoeffSet;

const FilterCoeffSet coeffTable[] = {
    {10.0f, {0.0019, 0.0038, 0.0019}, {-1.911, 0.915}},
    {20.0f, {0.0039, 0.0078, 0.0039}, {-1.815, 0.831}},
    // 更多预计算系数...
};

void updateFilterCoeff(BiquadSection* section, float cutoffFreq) {
    // 查找最近的预设频率
    uint8_t index = 0;
    float minDiff = fabsf(coeffTable[0].cutoffFreq - cutoffFreq);
    for(uint8_t i=1; i<sizeof(coeffTable)/sizeof(FilterCoeffSet); i++) {
        float diff = fabsf(coeffTable[i].cutoffFreq - cutoffFreq);
        if(diff < minDiff) {
            minDiff = diff;
            index = i;
        }
    }
    
    // 更新系数
    memcpy(section->b, coeffTable[index].b, 3*sizeof(float));
    memcpy(section->a, coeffTable[index].a, 2*sizeof(float));
}