从物理模型到代码：一阶与二阶RC滤波器的离散化推导与实践

1. 从物理模型到代码：RC滤波器的工程实践意义

在嵌入式系统和实时信号处理中，RC滤波器是最基础也最实用的工具之一。我第一次接触RC滤波器是在一个工业传感器项目中，当时需要处理来自压力传感器的信号，原始数据中混杂着高频噪声，导致控制算法频繁误动作。用硬件RC电路虽然简单，但调整参数需要更换电阻电容，非常麻烦。后来改用软件实现后，只需修改代码中的截止频率参数，就能快速适配不同工况。

软件实现的RC滤波器本质上是通过算法模拟硬件电路的行为。一阶RC滤波器对应着一阶微分方程，二阶则对应二阶微分方程。这种从物理模型到代码的转换过程，正是嵌入式工程师的必备技能。实际项目中，我遇到过采样周期不稳定导致滤波器效果变差的情况，也踩过浮点运算精度不足的坑，这些经验都会在后续章节分享。

2. 一阶RC低通滤波器的实现全流程

2.1 连续域模型与微分方程

一阶RC低通滤波器的硬件电路非常简单：一个电阻串联一个电容。根据基尔霍夫电压定律，可以得到：

V_in = V_R + V_C = R*i(t) + 1/C ∫i(t)dt

对两边求导并整理后得到标准的微分方程形式：

dV_out/dt + 1/(RC)*V_out = 1/(RC)*V_in

这个方程描述的就是电容电压V_out随时间变化的规律。在模拟电路中，当输入信号频率远低于1/(2πRC)时，信号能顺利通过；高于这个截止频率的信号则会被衰减。

2.2 离散化推导过程

在数字系统中，我们需要将这个连续时间方程转化为离散形式。采用后向差分法是最直观的方式：

dV_out/dt ≈ (V_out[k] - V_out[k-1])/T_s

代入原方程并整理：

V_out[k] = (1 - T_s/(RC + T_s)) * V_out[k-1] + (T_s/(RC + T_s)) * V_in[k]

令 α = T_s/(RC + T_s)，这就是我们常说的滤波系数。实际使用时更习惯用截止频率fc表示：

α = 2πfcT_s / (1 + 2πfcT_s)

2.3 C语言实现技巧

在嵌入式环境中实现时，有几点需要特别注意：

// 推荐使用宏定义实现，避免函数调用开销 #define LPF_ALPHA(fc, ts) (2*3.1415926f*(fc)*(ts)/(1 + 2*3.1415926f*(fc)*(ts))) #define LPF_UPDATE(y, x, a) ((y) = (1-(a))*(y) + (a)*(x)) // 使用时： float filtered_value = 0.0f; float new_sample = ADC_Read(); float alpha = LPF_ALPHA(10.0f, 0.001f); // fc=10Hz, Ts=1ms LPF_UPDATE(filtered_value, new_sample, alpha);

实测发现，在STM32F4系列MCU上，这种实现方式单个采样点的处理时间不到1μs。对于采样率1kHz的系统，计算负载几乎可以忽略不计。

2.4 参数选择与验证

在电机控制项目中，我曾用这种滤波器处理电流采样信号。关键经验是：

• 截止频率应设为信号带宽的3-5倍
• 采样周期T_s必须稳定，最好用硬件定时器触发
• 对于噪声特别大的场景，可以级联两个一阶滤波器

验证时可以用信号发生器+ADC采集，或者像下面这样用代码生成测试信号：

// 生成50Hz信号+100Hz噪声 for(int i=0; i<1000; i++) { float t = i*0.001f; float signal = sin(2*3.1415926f*50*t); float noise = 0.3f*sin(2*3.1415926f*100*t); float raw = signal + noise; LPF_UPDATE(filtered, raw, alpha); }

将原始信号和滤波后信号通过SWO或串口输出，用Python matplotlib绘制，能直观看到高频噪声被有效抑制。

3. 一阶RC高通滤波器的实现

3.1 数学模型差异点

高通滤波器的微分方程形式与低通类似，但输入输出关系不同：

V_out = V_in - V_C = V_in - 1/RC ∫V_out dt

对应的微分方程为：

dV_out/dt + 1/(RC)*V_out = dV_in/dt

离散化时要注意对输入信号的差分处理，最终得到：

V_out[k] = (1 - α)*V_out[k-1] + α*(V_in[k] - V_in[k-1])

3.2 实现中的常见问题

在实际编码时，我发现高通滤波器对初始条件更敏感。如果直接用上述公式，上电时V_in[-1]未知会导致输出异常。我的解决方案是：

float last_input = 0.0f; float hpf_output = 0.0f; void HPF_Update(float new_sample, float alpha) { hpf_output = (1-alpha)*hpf_output + alpha*(new_sample - last_input); last_input = new_sample; }

另一个问题是直流偏移。我曾用高通滤波去除ECG信号中的基线漂移，但发现相位延迟导致波形失真。后来改用截止频率0.5Hz的二阶高通滤波器解决了这个问题。

4. 二阶RC滤波器的进阶实现

4.1 双RC节模型的数学描述

二阶滤波器由两个RC节构成，传递函数更为复杂：

H(s) = 1 / (R1C1R2C2s² + (R1C1 + R2C2 + R1C2)s + 1)

采用双线性变换法离散化时，需要先将s域映射到z域：

s = (2/T_s)*(1-z⁻¹)/(1+z⁻¹)

经过一系列推导（过程较复杂，建议用符号计算工具辅助），最终得到差分方程：

V_out[k] = b0*V_in[k] + b1*V_in[k-1] + b2*V_in[k-2] - a1*V_out[k-1] - a2*V_out[k-2]

4.2 定点数优化技巧

在资源受限的MCU上，浮点运算可能成为瓶颈。我将算法改造成Q15格式定点数实现：

typedef int16_t q15_t; #define Q15_MUL(a,b) ((q15_t)(((int32_t)(a)*(b)) >> 15)) void IIR_Biquad(q15_t *y, q15_t x, const q15_t *coef) { static q15_t x_delay[2] = {0}; static q15_t y_delay[2] = {0}; q15_t acc = Q15_MUL(coef[0], x); // b0*x[n] acc += Q15_MUL(coef[1], x_delay[0]); // b1*x[n-1] acc += Q15_MUL(coef[2], x_delay[1]); // b2*x[n-2] acc -= Q15_MUL(coef[3], y_delay[0]); // a1*y[n-1] acc -= Q15_MUL(coef[4], y_delay[1]); // a2*y[n-2] // 更新延迟线 y_delay[1] = y_delay[0]; y_delay[0] = acc; x_delay[1] = x_delay[0]; x_delay[0] = x; *y = acc; }

测试发现，在Cortex-M0内核上，这种实现比浮点版本快5倍以上，而精度损失在可接受范围内。

5. 工程实践中的经验分享

在多个实际项目中验证过这些滤波器后，我总结出几条黄金法则：

1. 采样率选择：应至少是目标信号最高频率的10倍，同时截止频率不能超过采样率的1/5（奈奎斯特频率的1/2.5）

2. 参数调试：先用Python/scipy.signal生成理想响应曲线，再移植到嵌入式端。推荐我的调试步骤：

   • 用MATLAB计算理论系数
   • 在PC端验证算法逻辑
   • 移植到目标硬件后微调

3. 异常处理：加入以下保护措施：

   // 检测NaN/Inf if(!isfinite(filter_output)) { filter_output = 0; reset_delay_lines(); } // 抗溢出 if(filter_output > MAX_OUTPUT) filter_output = MAX_OUTPUT; if(filter_output < MIN_OUTPUT) filter_output = MIN_OUTPUT;

4. 性能优化：对于多通道系统，可以将滤波器设计为结构体形式：

   typedef struct { float alpha; float last_out; float last_in; // 高通需要 } Filter_1stOrder; void Filter_Update(Filter_1stOrder *f, float input) { f->last_out = (1-f->alpha)*f->last_out + f->alpha*input; }

最近在做一个无线传感器网络项目时，我将这些滤波器算法封装成可配置的模块，通过无线链路动态调整参数，极大方便了现场调试。这种从数学模型到实际产品的转化过程，正是嵌入式开发的魅力所在。