巴特沃斯滤波器MATLAB到C语言的嵌入式实现

1. 项目背景与核心价值

巴特沃斯滤波器作为经典的IIR滤波器类型，在信号处理领域有着广泛应用。其最大特点是通带内具有最大平坦的幅度特性，阻带内单调递减。这种特性使其特别适合对相位失真不敏感但要求通带平坦的应用场景，比如音频处理、生物信号采集、工业传感器信号调理等领域。

在实际工程中，我们经常面临这样的需求：算法工程师用MATLAB设计出完美的滤波器参数，但最终需要移植到嵌入式设备上用C语言实现。这个过程中涉及到系数计算、量化处理、定点化优化等一系列问题。很多工程师在这个转换过程中会遇到频率响应畸变、数值溢出、计算效率低下等典型问题。

这个项目要解决的正是这个工程实践中的痛点——如何把MATLAB设计的巴特沃斯滤波器完美地转化为可嵌入式部署的C代码。我将分享从理论设计到工程实现的完整链路，包括MATLAB的滤波器设计工具箱使用、系数提取方法、定点化技巧、以及C语言实现时的优化策略。

2. 巴特沃斯滤波器设计基础

2.1 滤波器参数确定

设计巴特沃斯滤波器首先需要明确几个关键参数：

• 通带截止频率(Wp)：通常以归一化频率表示(0-1对应0-Fs/2)
• 阻带截止频率(Ws)
• 通带最大衰减(Rp)：单位dB
• 阻带最小衰减(Rs)：单位dB
• 滤波器阶数(N)：可由MATLAB的buttord函数计算

在MATLAB中，典型的参数计算代码如下：

matlab复制Fs = 1000; % 采样率1kHz
Wp = 50/(Fs/2); % 通带50Hz
Ws = 100/(Fs/2); % 阻带100Hz
Rp = 3; % 通带波纹3dB 
Rs = 40; % 阻带衰减40dB
[N, Wn] = buttord(Wp, Ws, Rp, Rs); % 计算最小阶数和截止频率

2.2 传递函数形式

巴特沃斯滤波器的传递函数有两种形式：

1. 零极点增益形式(z,p,k)：

matlab复制[z,p,k] = butter(N, Wn);
[sos,g] = zp2sos(z,p,k); % 转换为二阶节形式

2. 直接系数形式(b,a)：

matlab复制[b,a] = butter(N, Wn);

对于高阶滤波器(N>2)，建议使用二阶节(SOS)形式，可以避免数值不稳定问题。这是因为高阶直接形式的滤波器在定点实现时容易产生系数敏感性问题，而二阶节串联形式具有更好的数值特性。

3. MATLAB设计实操详解

3.1 滤波器设计与分析

完整的MATLAB设计流程如下：

matlab复制% 1. 参数设定
Fs = 1000; Wp = 50/(Fs/2); Ws = 100/(Fs/2);
Rp = 3; Rs = 40;

% 2. 计算最小阶数
[N, Wn] = buttord(Wp, Ws, Rp, Rs);

% 3. 设计滤波器(使用SOS形式)
[sos, g] = butter(N, Wn, 'low');

% 4. 频率响应分析
freqz(sos, 512, Fs);
title(['N=' num2str(N) '阶巴特沃斯低通滤波器']);

设计完成后，关键是要分析：

• 通带波纹是否满足要求
• 阻带衰减是否达标
• 过渡带陡峭程度
• 相位非线性程度(巴特沃斯滤波器相位非线性较明显)

3.2 系数提取与格式化

将SOS系数转换为C语言可用的格式：

matlab复制% 提取二阶节系数
num_of_sections = size(sos,1);
coef = sos(:,1:3); % 分子系数
coef = [coef, sos(:,4:6)]; % 合并分母系数

% 生成C语言数组
fprintf('const float filter_coef[%d][6] = {\n', num_of_sections);
for i=1:num_of_sections
    fprintf('    {');
    fprintf('%.8e, %.8e, %.8e, %.8e, %.8e, %.8e', coef(i,:));
    if i~=num_of_sections
        fprintf('},\n');
    else
        fprintf('}\n');
    end
end
fprintf('};\n\n');
fprintf('const float filter_gain = %.8e;\n', g);

注意：MATLAB的sos矩阵排列顺序是[b0,b1,b2,1,a1,a2]，其中a0固定为1

4. C语言实现方案

4.1 直接型实现

对于二阶节形式的滤波器，每个section的实现如下：

c复制float butter_filter(float input, const float coef[6], float* w) {
    float output;
    
    // 计算输出
    output = coef[0]*input + coef[1]*w[0] + coef[2]*w[1];
    
    // 更新状态变量
    w[1] = w[0];
    w[0] = input - coef[4]*w[0] - coef[5]*w[1];
    
    return output;
}

// 多级串联调用
float process_sample(float x) {
    static float w[SECTIONS][2] = {0};
    float y = x;
    
    for(int i=0; i<SECTIONS; i++) {
        y = butter_filter(y, filter_coef[i], w[i]);
    }
    
    return y * filter_gain;
}

4.2 定点优化实现

对于资源受限的嵌入式系统，可以采用定点数实现：

c复制typedef int32_t fixed_point;
#define FRAC_BITS 28  // Q4.28格式

fixed_point float_to_fixed(float x) {
    return (fixed_point)(x * (1<<FRAC_BITS));
}

fixed_point fixed_mul(fixed_point a, fixed_point b) {
    return (fixed_point)(((int64_t)a * b) >> FRAC_BITS);
}

fixed_point butter_filter_fixed(fixed_point input, const fixed_point coef[6], fixed_point w[2]) {
    fixed_point output;
    
    // 64位中间结果防止溢出
    int64_t temp = (int64_t)coef[0]*input + (int64_t)coef[1]*w[0] + (int64_t)coef[2]*w[1];
    output = (fixed_point)(temp >> FRAC_BITS);
    
    // 更新状态变量
    w[1] = w[0];
    temp = (int64_t)input - ((int64_t)coef[4]*w[0] + (int64_t)coef[5]*w[1]);
    w[0] = (fixed_point)(temp >> FRAC_BITS);
    
    return output;
}

重要提示：定点数实现时要特别注意：

1. 乘法后的位数扩展(FRAC_BITS*2)
2. 中间结果使用更大位宽存储(如int64_t)
3. 系数的动态范围，必要时进行缩放

5. 工程实践中的关键问题

5.1 稳定性问题

巴特沃斯滤波器虽然是理论上无条件稳定的，但在数字实现时可能因以下原因失稳：

1. 系数量化误差：高Q值极点对系数误差敏感
2. 运算舍入误差：递归部分的误差累积

解决方案：

• 采用二阶节串联结构
• 增加保护位(如32位系数，64位累加器)
• 定期重置状态变量(如检测到溢出时)

5.2 频率响应偏差

定点实现时常见的频率响应偏差来源：

1. 系数量化：特别是极点附近系数
2. 运算截断：递归部分误差影响更大

调试方法：

1. 在MATLAB中模拟量化效果：

matlab复制q_coef = round(coef * 2^16) / 2^16; % 模拟16位量化
freqz(q_coef, 512, Fs);

2. 在C代码中输出中间状态，与MATLAB仿真对比

5.3 计算效率优化

针对不同处理器的优化策略：

ARM Cortex-M系列：

• 使用CMSIS-DSP库的biquad函数
• 启用硬件FPU(如有)
• 使用SIMD指令并行处理多个样本

DSP处理器：

• 利用MAC(乘累加)指令
• 使用循环缓冲减少内存访问
• 展开内层循环

典型优化代码示例(Cortex-M)：

c复制#include "arm_math.h"

arm_biquad_casd_df1_inst_f32 S;
float32_t pState[2*SECTIONS];
float32_t pCoeffs[5*SECTIONS]; // {b0,b1,b2,a1,a2}

void init_filter() {
    // 转换系数格式
    for(int i=0; i<SECTIONS; i++) {
        pCoeffs[5*i] = filter_coef[i][0];
        pCoeffs[5*i+1] = filter_coef[i][1];
        pCoeffs[5*i+2] = filter_coef[i][2];
        pCoeffs[5*i+3] = filter_coef[i][4];
        pCoeffs[5*i+4] = filter_coef[i][5];
    }
    
    arm_biquad_cascade_df1_init_f32(&S, SECTIONS, pCoeffs, pState);
}

float process_sample_optimized(float x) {
    float y;
    arm_biquad_cascade_df1_f32(&S, &x, &y, 1);
    return y * filter_gain;
}

6. 验证与测试方法

6.1 频域验证

1. 白噪声测试法：

matlab复制% MATLAB生成测试信号
x = randn(1,10000); % 白噪声
y = filter(b,a,x); % 理想输出
% 将x保存为文件供C程序读取

% C程序处理x后保存输出y_c
% MATLAB比较
[y_c, Fs] = audioread('output_c.bin');
figure;
pwelch(y, [], [], [], Fs); hold on;
pwelch(y_c, [], [], [], Fs);
legend('MATLAB','C实现');

2. 正弦扫描法：

matlab复制freqs = logspace(0, log10(Fs/2), 100);
gain = zeros(size(freqs));
for i=1:length(freqs)
    x = sin(2*pi*freqs(i)*(0:999)/Fs);
    y = process_samples(x); % C实现的处理函数
    gain(i) = rms(y(900:end))/rms(x(900:end));
end
semilogx(freqs, 20*log10(gain)); grid on;

6.2 时域验证

脉冲响应测试：

matlab复制x = [1, zeros(1,999)]; % 单位脉冲
y_matlab = filter(b,a,x);
y_c = process_samples(x); % C实现

figure;
stem(0:99, y_matlab(1:100)); hold on;
stem(0:99, y_c(1:100), 'r');
legend('MATLAB','C实现');

阶跃响应测试：

matlab复制x = ones(1,1000); % 单位阶跃
y_matlab = filter(b,a,x);
y_c = process_samples(x);

figure;
plot(y_matlab); hold on;
plot(y_c, 'r');
legend('MATLAB','C实现');

7. 高级话题与扩展

7.1 自动代码生成

对于频繁修改的设计，可以使用MATLAB Coder自动生成C代码：

matlab复制% 定义入口函数
function y = process_sample(x)
    persistent filter
    if isempty(filter)
        [sos,g] = butter(6, 0.2);
        filter = dfilt.df2sos(sos, g);
    end
    y = filter(x);
end

% 配置并生成代码
cfg = coder.config('lib');
cfg.GenerateReport = true;
codegen -config cfg process_sample -args {0}

生成的代码包含：

• 滤波器初始化函数
• 步进处理函数
• 状态存储结构体

7.2 动态参数调整

对于需要在线调整参数的场景：

1. 系数重计算优化：预先计算多组系数，运行时切换
2. 参数插值：在两组有效系数间平滑过渡
3. 使用Farrow结构实现连续可调

示例代码框架：

c复制typedef struct {
    float coef[6];
    float w[2];
} BiquadSection;

typedef struct {
    BiquadSection sections[MAX_SECTIONS];
    float gain;
    int section_count;
} ButterworthFilter;

void update_cutoff(ButterworthFilter* f, float new_wn) {
    // 重新计算所有系数 -- 实际应用中可能需要离线计算
    float new_coef[MAX_SECTIONS][6];
    design_butterworth(f->section_count, new_wn, new_coef);
    
    // 平滑过渡
    for(int s=0; s<f->section_count; s++) {
        for(int i=0; i<6; i++) {
            f->sections[s].coef[i] += 0.1*(new_coef[s][i] - f->sections[s].coef[i]);
        }
    }
}

7.3 多速率处理

对于截止频率远低于采样率的场景，可结合抽取/插值：

1. 先抽取降低采样率
2. 在低采样率下滤波
3. 再插值恢复采样率

优点：

• 减少计算量
• 降低对滤波器阶数的要求
• 减少量化误差影响

实现示例：

c复制#define DECIMATION 4

float process_sample_decim(float x) {
    static float buffer[DECIMATION];
    static int ptr = 0;
    static int phase = 0;
    float y = 0;
    
    buffer[ptr] = x;
    ptr = (ptr+1) % DECIMATION;
    
    if(++phase == DECIMATION) {
        phase = 0;
        float x_decim = fir_decimator(buffer); // 抗混叠FIR
        y = butter_lowpass(x_decim); // 工作在Fs/DECIMATION
        y = fir_interpolator(y); // 插值滤波器
    }
    
    return y;
}

8. 实际项目经验分享

在最近的心电信号(ECG)处理项目中，我们遇到了这样的需求：需要在STM32F4系列MCU上实现50Hz工频陷波和0.5-40Hz带通滤波。经过实践验证，以下方案效果最佳：

1. 陷波器采用二阶IIR：

matlab复制wo = 50/(Fs/2); bw = wo/35;
[b,a] = iirnotch(wo, bw);

2. 带通采用4阶巴特沃斯(2个二阶节)：

matlab复制[b,a] = butter(2, [0.5 40]/(Fs/2), 'bandpass');
sos = tf2sos(b,a);

关键经验：

• 陷波器先于带通滤波器应用，避免非线性相位影响QRS波检测
• 使用Q31定点数实现，在224MHz的STM32F427上仅消耗0.8% CPU资源
• 每个样本处理时间小于300时钟周期
• 通过定期(每5秒)重置状态变量，避免长期运行的误差累积

调试中发现的问题：最初使用直接型实现时，带通滤波器在高温环境下偶尔会出现不稳定。改为二阶节形式后问题消失，即使环境温度达到85°C也能稳定工作。