从MATLAB到单片机：手把手教你用C语言移植巴特沃斯滤波器（附完整代码）
从MATLAB到单片机：手把手教你用C语言移植巴特沃斯滤波器（附完整代码）
在嵌入式系统开发中，数字信号处理是一个绕不开的话题。当你需要在资源受限的单片机上实现实时滤波功能时，直接从MATLAB仿真到C语言移植往往是最直接的路径。本文将聚焦巴特沃斯滤波器的实现，带你完整走通从系数计算到单片机部署的全流程。
1. 巴特沃斯滤波器基础与MATLAB设计
巴特沃斯滤波器的核心特性在于其通带内具有最大平坦的幅度响应。这种特性使其成为许多嵌入式应用的理想选择，尤其是对相位失真不敏感的场景。
在MATLAB中设计一个4阶低通巴特沃斯滤波器，采样频率1kHz，截止频率100Hz：
fs = 1000; % 采样频率(Hz) fc = 100; % 截止频率(Hz) order = 4; % 滤波器阶数 [b, a] = butter(order, fc/(fs/2), 'low');运行后会得到两组系数：
- b：分子系数（前馈部分）
- a：分母系数（反馈部分）
注意：MATLAB返回的a系数第一个元素总是1，实际实现时需要去除。
2. 系数转换与定点数优化
嵌入式系统往往面临浮点运算性能不足的问题，将系数转换为定点数是常见优化手段。
2.1 浮点系数导出
将MATLAB系数保存为C语言头文件：
fid = fopen('filter_coeffs.h','w'); fprintf(fid,'#define FILTER_ORDER %d\n', order); fprintf(fid,'const float B[] = {'); fprintf(fid,'%.8ff,', b(1:end-1)); fprintf(fid,'%.8ff};\n', b(end)); fprintf(fid,'const float A[] = {'); fprintf(fid,'%.8ff,', a(2:end-1)); fprintf(fid,'%.8ff};\n', a(end)); fclose(fid);2.2 定点数转换
对于STM32等Cortex-M系列单片机，Q格式定点数能显著提升性能：
#define Q 15 // Q15格式，16位有符号数 int16_t B_q15[FILTER_ORDER+1] = { (int16_t)(B[0] * (1<<Q)), (int16_t)(B[1] * (1<<Q)), // ...其余系数 }; int16_t A_q15[FILTER_ORDER] = { (int16_t)(A[0] * (1<<Q)), // ...其余系数 };提示：Q格式运算时需注意溢出问题，中间结果可能需要32位容器存储。
3. C语言实现方案对比
3.1 直接型(Direct Form I)实现
最直观的实现方式，但需要较多存储空间：
float direct_form_i(float input, float *x_buf, float *y_buf) { // 移位历史数据 for(int i = FILTER_ORDER; i > 0; i--) { x_buf[i] = x_buf[i-1]; y_buf[i] = y_buf[i-1]; } x_buf[0] = input; // 计算输出 float output = 0.0f; for(int i = 0; i <= FILTER_ORDER; i++) { output += B[i] * x_buf[i]; } for(int i = 0; i < FILTER_ORDER; i++) { output -= A[i] * y_buf[i+1]; } y_buf[0] = output; return output; }3.2 直接型II(Direct Form II)实现
更节省内存的实施方案，适合资源受限环境：
float direct_form_ii(float input, float *w) { // 计算中间变量 w[0] = input; for(int i = 1; i <= FILTER_ORDER; i++) { w[0] -= A[i-1] * w[i]; } // 计算输出 float output = 0.0f; for(int i = 0; i <= FILTER_ORDER; i++) { output += B[i] * w[i]; } // 移位状态变量 for(int i = FILTER_ORDER; i > 0; i--) { w[i] = w[i-1]; } return output; }内存占用对比表：
| 实现方式 | 状态变量数量 | 适合场景 |
|---|---|---|
| 直接型I | 2N+2 | 代码可读性优先 |
| 直接型II | N+1 | 内存受限系统 |
| 级联二阶节 | 2N | 高数值稳定性需求 |
4. STM32实战优化技巧
4.1 使用ARM DSP库加速
对于Cortex-M4/M7等带DSP指令集的芯片，利用CMSIS-DSP库可大幅提升性能：
#include "arm_math.h" arm_biquad_casd_df1_inst_f32 filter; float32_t stateBuffer[2*NUM_STAGES]; // 每个二阶节需要4个状态变量 void init_filter() { float32_t coeffs[5*NUM_STAGES] = { // b0, b1, b2, a1, a2 每二阶节5个系数 0.1234, 0.2468, 0.1234, -1.2345, 0.4567, // ...更多二阶节系数 }; arm_biquad_cascade_df1_init_f32(&filter, NUM_STAGES, coeffs, stateBuffer); } float process_sample(float input) { float output; arm_biquad_cascade_df1_f32(&filter, &input, &output, 1); return output; }4.2 实时性保障措施
- 定时器触发ADC采样：确保严格等间隔采样
- DMA双缓冲：实现无阻塞数据处理
- 优先级管理：
- ADC中断：最高优先级
- 滤波处理：中等优先级
- 结果输出：最低优先级
// 示例：STM32 HAL库配置 void HAL_ADC_ConvCpltCallback(ADC_HandleTypeDef* hadc) { float raw = (float)HAL_ADC_GetValue(hadc) * 3.3f / 4095.0f; float filtered = process_sample(raw); send_to_output(filtered); // 非阻塞方式发送 }5. 常见问题与调试技巧
5.1 数值不稳定问题
症状：输出逐渐饱和或出现NaN 解决方案：
- 采用级联二阶节结构
- 加入输出限幅
- 定期重置状态变量
#define MAX_OUTPUT 3.3f // 假设系统供电电压 float safe_filter(float input) { float output = process_sample(input); if(isnan(output)) { reset_filter_states(); return 0.0f; } return fmaxf(fminf(output, MAX_OUTPUT), -MAX_OUTPUT); }5.2 频率响应验证
使用信号发生器+示波器组合验证：
- 输入正弦扫频信号
- 测量输出幅度变化
- 绘制实际频率响应曲线
实用技巧：可以在代码中加入自测试模式，自动输出阶跃信号或扫频信号。
6. 完整代码示例
以下是一个针对STM32H7的完整实现：
// filter.h #pragma once #include <stdint.h> #define FILTER_ORDER 4 #define Q15_SHIFT 15 typedef struct { int16_t B[FILTER_ORDER+1]; // 分子系数 int16_t A[FILTER_ORDER]; // 分母系数 int32_t w[FILTER_ORDER+1]; // 状态变量(Q15) } IIR_Filter; void IIR_Init(IIR_Filter *filt); int16_t IIR_Process(IIR_Filter *filt, int16_t input);// filter.c #include "filter.h" void IIR_Init(IIR_Filter *filt) { // 示例系数(实际应从MATLAB导出) const float B_float[] = {0.00041655, 0.00166620, 0.00249930, 0.00166620, 0.00041655}; const float A_float[] = {-2.6861574, 2.4196556, -0.7301653}; for(int i=0; i<=FILTER_ORDER; i++) { filt->B[i] = (int16_t)(B_float[i] * (1<<Q15_SHIFT)); } for(int i=0; i<FILTER_ORDER; i++) { filt->A[i] = (int16_t)(A_float[i] * (1<<Q15_SHIFT)); } for(int i=0; i<=FILTER_ORDER; i++) { filt->w[i] = 0; } } int16_t IIR_Process(IIR_Filter *filt, int16_t input) { // 转换为Q15格式的中间变量 int32_t w0 = (int32_t)input << Q15_SHIFT; // 反馈部分计算 for(int i=1; i<=FILTER_ORDER; i++) { w0 -= (int32_t)filt->A[i-1] * filt->w[i] >> Q15_SHIFT; } // 前馈部分计算 int32_t output = 0; for(int i=0; i<=FILTER_ORDER; i++) { output += (int32_t)filt->B[i] * filt->w[i] >> Q15_SHIFT; } // 更新状态 for(int i=FILTER_ORDER; i>0; i--) { filt->w[i] = filt->w[i-1]; } filt->w[0] = w0; return (int16_t)(output >> Q15_SHIFT); }在STM32CubeIDE中集成时，建议将滤波处理放在定时器中断中，确保严格实时性。对于需要更高阶滤波器的场景，可以考虑将多个二阶节级联，每个二阶节使用上述结构独立实现。