ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

一阶IIR滤波器实战:差分方程系数计算与嵌入式C语言实现

一阶IIR滤波器实战:差分方程系数计算与嵌入式C语言实现 1. 一阶IIR滤波器到底在做什么1.1 从一个生活场景说起你拿手机录一段语音回放的时候发现底噪很大嘶嘶的声音让人难受。你想把它弄干净但又不想花太多计算资源。这时候一阶IIR滤波器就是最顺手的那把刀。它的核心逻辑特别朴素当前输出 一部分当前输入 一部分上一次输出。用数学写出来就是y[n] b0 · x[n] b1 · x[n-1] - a1 · y[n-1]这就是一阶IIR滤波器差分方程的标准形式。x[n]是当前采样进来的原始数据y[n]是滤波后输出的数据y[n-1]是上一次的输出结果。系数b0、b1、a1决定了这个滤波器到底是低通、高通还是带通。我第一次接触这个公式的时候也觉得抽象后来想明白了一件事它本质上就是一个带记忆的加权平均器。普通平均是取最近N个数的算术平均而IIR滤波器用反馈的方式让“记忆”以指数衰减的形式保留下来不需要存N个历史数据一个y[n-1]就够了。这就是它计算量小、内存占用低的根本原因。1.2 为什么一阶IIR在嵌入式里这么常见做过单片机ADC采样的人都知道原始信号里混杂的噪声有多烦人。你采一个温度传感器理论上变化应该很缓慢但实际读出来的值可能上下跳动好几度。这时候用一阶IIR低通滤波几行代码就能把数据磨平。它的优势非常明确计算量极小每次采样只需要2次乘法和2次加法Cortex-M0都能轻松跑内存占用极低只需要保存一个历史输出值y[n-1]和一个历史输入值x[n-1]参数调节直观截止频率和系数之间有明确的数学关系改一个系数就能调整滤波强度实时性好没有缓冲区来一个采样点算一个延迟固定且极小但它的局限也很明显滚降斜率只有-20dB/十倍频程过渡带比较宽。如果你需要更陡的截止特性一阶不够用得上二阶甚至更高阶。不过在大量实际场景里一阶已经能把问题解决得七七八八了。1.3 低通、高通、带通在一阶条件下的关系这里要先说清楚一个事实严格意义上的一阶IIR滤波器只能实现低通和高通。带通需要至少二阶。但为什么标题里提到了带通配置因为在实际工程中我们经常用“一阶低通 一阶高通”级联的方式来近似实现带通或者用一阶差分方程配合特定系数组合来获得带通效果。这个点很关键很多初学者看到“一阶带通”会困惑。我的理解是如果你要的是真正的带通滤波器至少需要两个极点也就是二阶结构。但如果你只是想要一个“去掉直流分量同时抑制高频噪声”的效果一个一阶差分方程配合合适的系数就能做到近似。下面这张表把三种类型在一阶条件下的实现方式理清楚滤波类型实现方式典型应用低通单个一阶IIR传感器去噪、ADC平滑高通单个一阶IIR去除直流偏置、交流耦合带通低通级联高通等效二阶音频分频、振动信号提取2. 差分方程系数怎么算出来的2.1 从模拟滤波器到数字滤波器的映射一阶IIR滤波器的系数不是拍脑袋定的它来自模拟滤波器原型的离散化。最常用的方法是双线性变换和后向差分。我平时用后向差分比较多因为计算简单而且在低频段精度够用。模拟一阶低通滤波器的传递函数是H(s) ωc / (s ωc)其中ωc 2πfcfc是截止频率。用后向差分法令s (1 - z⁻¹) / TT是采样周期代入后整理得到y[n] y[n-1] α · (x[n] - y[n-1])这就是嵌入式里最常见的一阶低通滤波公式。α是滤波系数取值范围0到1。α越小滤波越强但响应越慢α越大响应越快但滤波效果越弱。α和截止频率的关系是α 2πfc · T / (1 2πfc · T)其中T 1/fsfs是采样率。这个公式我建议你记住因为实际调参的时候全靠它。2.2 低通滤波系数的计算实例假设你的ADC采样率是1kHz你想让截止频率在10Hz左右。代入公式T 1/1000 0.001秒2πfc·T 2 × 3.14159 × 10 × 0.001 0.0628α 0.0628 / (1 0.0628) 0.059所以α取0.06左右。实际写代码的时候为了避开浮点运算可以用定点数α 6/100或者用移位操作α ≈ 1/16 0.0625。这样在单片机上跑起来飞快。我实测过α0.06的时候一个阶跃信号大概需要30到40个采样周期才能稳定到最终值。如果你觉得响应太慢把α调到0.1到0.2之间滤波效果会弱一些但响应更快。这个取舍完全取决于你的应用场景。2.3 高通滤波系数的推导一阶高通滤波器的差分方程形式是y[n] α · (y[n-1] x[n] - x[n-1])注意这里的α含义和低通不一样。高通滤波器的α越接近1截止频率越低α越接近0截止频率越高。从模拟高通传递函数H(s) s / (s ωc)出发用同样的后向差分法可以得到α 1 / (1 2πfc · T)这个公式和低通的α刚好是互补关系。低通的α_hp 1 - α_lp。所以如果你已经算好了低通系数高通系数直接减一下就行。2.4 带通配置的系数组合策略前面说了一阶带通实际上是低通和高通的级联。具体做法是先对信号做一阶高通截止频率fL再对结果做一阶低通截止频率fH要求fL fH两者之间的频段就是通带两个滤波器级联后总的传递函数是两个一阶传递函数相乘滚降斜率变成-40dB/十倍频程但严格来说这已经是二阶系统了。不过在工程语境下大家还是习惯叫它“一阶带通配置”因为每个环节都是一阶的。级联的时候有个坑要注意两个滤波器会互相影响。如果你先高通后低通高通输出的直流分量已经被去掉了低通的输入范围会变小实际截止特性会和单独设计时略有偏差。我的经验是设计时把两个截止频率拉开一些至少保持2倍以上的间隔这样互相影响可以忽略。3. 代码实现与参数调优3.1 低通滤波的C语言实现先看最基础的低通实现。这段代码可以直接抄到你的单片机项目里typedef struct { float alpha; float prev_output; } LowPassFilter; float lowpass_update(LowPassFilter *f, float input) { f-prev_output f-prev_output f-alpha * (input - f-prev_output); return f-prev_output; }初始化的时候把prev_output设成第一次的输入值避免启动时的冲击。alpha根据前面的公式算出来填进去。如果你用的是没有FPU的单片机把float换成int32_t用定点数运算typedef struct { int32_t alpha_q15; // Q15格式范围0~32767 int32_t prev_output; } LowPassFilterFixed; int32_t lowpass_update_fixed(LowPassFilterFixed *f, int32_t input) { int32_t diff input - f-prev_output; int32_t delta (f-alpha_q15 * diff) 15; f-prev_output delta; return f-prev_output; }Q15格式就是把0到1的浮点数映射到0到32767的整数。alpha0.06对应alpha_q15 0.06 × 32768 ≈ 1966。这样一次滤波只需要一次乘法和一次移位Cortex-M0都能在几个时钟周期内完成。3.2 高通滤波的代码实现高通滤波器的代码稍微多一行因为要保存历史输入typedef struct { float alpha; float prev_output; float prev_input; } HighPassFilter; float highpass_update(HighPassFilter *f, float input) { f-prev_output f-alpha * (f-prev_output input - f-prev_input); f-prev_input input; return f-prev_output; }注意prev_input每次都要更新这个很容易漏掉。我刚开始写的时候忘了更新prev_input结果滤波器变成了一个奇怪的积分器输出直接饱和了。排查了半天才发现是这行的问题。3.3 带通滤波的级联实现带通就是把两个滤波器串起来typedef struct { HighPassFilter hp; LowPassFilter lp; } BandPassFilter; float bandpass_update(BandPassFilter *f, float input) { float hp_out highpass_update(f-hp, input); return lowpass_update(f-lp, hp_out); }初始化的时候高通和低通的prev_output都设成0prev_input也设成0。如果信号有直流偏置高通输出刚开始会有个瞬态几个采样周期后就稳定了。3.4 参数调优的实操经验调参这件事理论公式给你一个起点但最终值一定要在实际信号上试。我一般按这个流程走先用公式算出理论α值作为初始值采集一段实际信号跑一遍滤波看输出波形如果噪声还是大把α调小20%到30%如果响应太慢把α调大20%到30%反复两三轮找到平衡点有个经验数据可以参考对于温度传感器这种变化极慢的信号α取0.01到0.05就够了对于音频信号α取0.1到0.3比较合适对于振动信号α可以取到0.3到0.5。还有一个技巧动态调整α。信号变化剧烈的时候用大α快速跟踪信号平稳的时候用小α深度滤波。实现方式是检测输入和输出的差值差值大就增大α差值小就减小α。这个在电机控制里特别有用。4. 实际项目中踩过的坑和排查方法4.1 滤波器发散和振荡问题一阶IIR滤波器理论上不会发散因为极点都在单位圆内。但实际写代码的时候如果系数算错了或者定点数溢出输出可能会越来越大直到饱和。我遇到过一次用Q15定点数实现低通alpha_q15设成了35000超过了32767的范围结果移位后符号位出错输出直接飞到负的最大值。排查方法很简单把alpha_q15打印出来看看是不是在0到32767之间。还有一种情况是级联滤波时前一级的输出超出了后一级的输入范围。比如高通输出可能有负值但低通滤波器如果用的是无符号整数就会出问题。解决办法是全程用有符号数或者在高通输出后加一个偏置。4.2 截止频率和实际效果对不上理论计算的截止频率是-3dB点但实际信号经过滤波后你可能感觉截止频率偏了。原因通常有两个一是采样率不稳定。如果你的ADC采样率实际是900Hz但你以为的1000Hz算出来的α对应的截止频率就会偏。解决办法是用定时器触发ADC确保采样周期精确。二是信号本身不是正弦波。-3dB的定义是基于正弦稳态响应的如果你的信号是方波或者脉冲滤波后的形状变化不能用截止频率直接描述。这时候要看群延迟和阶跃响应。4.3 启动瞬态的处理滤波器刚启动的时候prev_output是0但实际信号可能有个非零的初始值。这会导致输出从0慢慢爬到真实值形成一个瞬态。对于慢速信号这个瞬态可能持续好几秒。我的处理方法是第一次调用滤波函数时直接把prev_output设成input。这样启动瞬间输出就等于输入没有爬升过程。代码里加一个初始化标志位就行if (!f-initialized) { f-prev_output input; f-prev_input input; f-initialized 1; return input; }4.4 常见问题速查表现象可能原因解决方法输出一直为0alpha为0或prev_output未初始化检查alpha值初始化prev_output输出饱和到最大值定点数溢出或alpha过大检查alpha范围增加饱和保护滤波效果不明显alpha太大减小alpha重新计算截止频率响应太慢alpha太小增大alpha或改用动态alpha输出有直流偏置高通滤波器未正确初始化初始化prev_input和prev_output级联后效果变差两级截止频率太接近拉开截止频率间隔至2倍以上4.5 一个容易被忽略的细节数据类型选择float用起来方便但在没有FPU的平台上速度慢。我一般这样选有FPU的Cortex-M4/M7直接用float代码简洁无FPU的Cortex-M0/M3用Q15定点数速度快10倍以上8位单片机用Q7或Q8格式注意溢出保护定点数的移位操作要注意符号扩展。右移的时候有符号数是算术右移无符号数是逻辑右移。C语言里对signed类型做右移是implementation-defined保险的做法是先转成unsigned再移位或者用除法代替。5. 从一阶到二阶的扩展思路5.1 什么时候需要二阶一阶滤波器的滚降是-20dB/十倍频程如果你需要更陡的截止特性比如在截止频率后快速衰减一阶就不够了。典型的场景是音频分频器需要精确的频率分割抗混叠滤波需要在奈奎斯特频率前快速衰减窄带信号提取需要抑制邻近频道的干扰这时候需要二阶IIR滤波器滚降变成-40dB/十倍频程。二阶的实现方式和一阶类似只是差分方程多了一项y[n] b0·x[n] b1·x[n-1] b2·x[n-2] - a1·y[n-1] - a2·y[n-2]系数计算更复杂一些但思路是一样的从模拟原型出发用双线性变换或后向差分离散化。5.2 二阶低通的设计要点二阶低通滤波器有两种常见拓扑Sallen-Key和多重反馈。数字实现的话直接用双二阶节biquad结构。每个biquad就是一个二阶IIR级联多个biquad可以实现更高阶的滤波器。设计二阶低通的时候除了截止频率还要确定品质因数Q。Q决定了截止频率附近的峰值特性Q 0.707巴特沃斯最平坦的通带响应Q 0.707截止频率附近有峰值Q 0.707过渡带更平缓我一般先用巴特沃斯响应Q0.707然后根据实际效果微调。5.3 一阶和二阶的取舍一阶够用的时候不要上二阶。二阶的计算量是一阶的两倍多而且参数更多调参更麻烦。我判断的标准是如果噪声频率远高于信号频率一阶足够如果噪声频率接近信号频率需要二阶如果对相位线性度有要求考虑FIR而不是IIR一阶IIR的最大优势就是简单。在资源受限的嵌入式系统里简单就是可靠。我做过一个项目用一阶IIR滤波处理心率信号效果完全满足需求代码只有十几行。后来有人建议换成二阶我试了一下效果提升有限但代码复杂度翻倍最后还是用回了一阶。5.4 带通采样定理的关联热词里提到了带通采样定理这里简单说一下它和带通滤波器的关系。带通采样定理说的是如果一个带通信号的带宽是B中心频率是f0那么采样率只要大于2B就能无失真采样不需要满足奈奎斯特的2f0。这个定理在射频和通信领域用得很多。但要注意带通采样后信号的频谱会搬移后续处理需要配合带通滤波器把需要的频段提取出来。一阶带通配置在这里可以作为预滤波把带外噪声先压一压减轻后续处理的压力。6. 工程实践中的几个关键决策6.1 采样率怎么定采样率的选择直接影响滤波器的设计。采样率越高同样的截止频率对应的α越小滤波效果越强但计算量也越大。我的经验是采样率取截止频率的10到20倍比如你要滤除10Hz以上的噪声采样率取100Hz到200Hz就够了。取太高浪费计算资源取太低会导致混叠。如果信号本身带宽很宽比如音频信号20Hz到20kHz那采样率至少40kHz。这时候一阶IIR的截止频率如果设在20kHzα会接近1滤波效果很弱。这种情况下应该用更高阶的滤波器。6.2 系数用浮点还是定点这个决策取决于你的硬件平台和精度要求。我整理了一个对比表维度浮点定点精度高取决于位宽速度无FPU慢快速度有FPU快快代码复杂度低中溢出风险低高可移植性好一般我的建议是有FPU就用浮点没有FPU就用Q15定点。Q15的精度对于大多数传感器信号足够了16位有效位数对应约90dB的动态范围。6.3 滤波器放在系统的哪个位置滤波器在信号链中的位置很重要。一般有两种选择方案AADC采样后立即滤波优点是噪声在进入后续处理之前就被抑制了后续算法不会被噪声干扰。缺点是滤波器的延迟会影响控制环路的响应速度。方案B在特定处理环节前滤波比如在做FFT之前滤波或者在阈值判断之前滤波。优点是针对性强不影响其他环节。缺点是需要保存原始数据内存占用大。我通常选方案A因为一阶IIR的延迟很小对大多数应用没有影响。只有在高速控制环路里才会考虑把滤波器放在反馈路径之外。6.4 实测验证方法滤波器设计完之后一定要实测验证。我的验证流程是阶跃响应测试给一个阶跃输入看输出的上升时间和过冲正弦扫频测试从低到高扫频看幅频特性是否符合预期实际信号测试用真实传感器数据跑一遍看滤波效果长时间稳定性测试跑几个小时看输出是否漂移或发散阶跃响应测试最简单也最有效。如果阶跃响应不对后面的测试都不用做了。我一般用信号发生器产生一个方波频率设成截止频率的十分之一观察输出波形。如果输出是平滑的指数上升说明滤波器工作正常。正弦扫频测试需要信号发生器或者DAC输出。如果没有设备可以用软件生成正弦波数据直接喂给滤波函数把输出打印出来画图。Python的matplotlib画个伯德图很方便。实际信号测试是最关键的。理论计算再漂亮实际信号上效果不好就是白搭。我习惯把滤波前后的数据都保存下来用Excel或者Python画对比图直观地看滤波效果。长时间稳定性测试容易被忽略但很重要。有些滤波器在短时间测试时正常跑几个小时后因为浮点累积误差或者定点溢出输出会慢慢漂移。我遇到过一次定点数滤波跑了半天后输出偏了5%原因是每次移位都丢失低位累积误差越来越大。解决办法是定期用原始值重置滤波器状态或者用更高位宽的定点数。6.5 一个完整的参数计算示例最后给一个完整的计算示例把前面的内容串起来。需求ADC采样率1kHz信号带宽0到5Hz需要滤除5Hz以上的噪声同时去除0.1Hz以下的直流漂移。低通设计截止频率fc 5HzT 0.001s2πfcT 0.0314α_lp 0.0314 / 1.0314 0.0305取α_lp 0.03高通设计截止频率fc 0.1Hz2πfcT 0.000628α_hp 1 / (1 0.000628) 0.9994取α_hp 0.999级联顺序先高通后低通。因为高通先去掉直流低通的输入范围更集中滤波效果更好。验证阶跃响应上升时间约1/(2π×5) 32ms对应32个采样周期。高通的时间常数约1/(2π×0.1) 1.6s对应1600个采样周期。整体响应时间由高通决定约1.6秒稳定。这个配置我实际用过处理称重传感器的信号效果很稳。称重信号本身变化慢1.6秒的稳定时间完全可以接受滤波后的数据波动从±5g降到了±0.5g。7. 写在最后的一些个人体会一阶IIR滤波器看起来简单但真正用好需要理解它的数学本质和工程约束。我见过太多人直接抄一段代码alpha随便填个0.1然后抱怨滤波效果不好。问题不在代码在于没有根据实际信号特性去计算参数。我的建议是每次用一阶IIR之前先花五分钟算一下α再花十分钟在实际信号上验证。这十五分钟的投资能帮你省下几个小时的调试时间。还有一点不要迷信理论公式。理论计算给你一个起点但实际信号往往有各种非理想因素。我调过的一个电机电流信号理论算出来α0.05但实际用0.02效果更好因为电流信号里的噪声频率比预期的高。所以理论指导实践实践修正理论这个循环走两遍你就能找到最优参数。最后分享一个小技巧如果你不确定α取多少先用0.1跑一遍看输出波形。如果噪声还是大减半如果响应太慢加倍。两三次就能收敛到合适的值。这个方法虽然土但特别管用。
返回列表