ARTICLE DETAIL

资讯详情

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

一阶IIR滤波器差分方程推导与嵌入式实现:低通、高通、带通完整指南

一阶IIR滤波器差分方程推导与嵌入式实现:低通、高通、带通完整指南 一阶IIR滤波器在嵌入式信号处理里出现的频率远比很多人想象中高。做传感器数据采集的、搞电机控制的、写音频固件的几乎都会在某个时刻需要从ADC原始数据里把高频噪声压下去或者把直流偏置剥出来。这时候如果直接上FFT或者高阶FIR要么算力吃紧要么引入的群延迟让控制环路直接震荡。一阶IIR的好处就在这——一个乘法、两个加法、一个历史状态变量就能跑出相当可用的滤波效果在8位单片机上都能轻松实时运行。但问题也恰恰出在“简单”上。很多人从网上抄一段y a*x (1-a)*y_prev就用了低通勉强能跑一旦要改成高通或者带通就懵了不知道差分方程该怎么变形系数怎么定截止频率和采样率之间是什么关系。更麻烦的是一阶IIR的差分方程形式在不同资料里写法差异很大有的用b0/b1/a1系数有的用alpha有的直接给RC离散化公式初学者很容易被绕进去。这篇内容就是把这三种滤波器的差分方程从推导到落地完整拆一遍。我会从最基础的RC电路离散化讲起把低通、高通、带通的系数计算、代码实现、参数调优全部展开每个公式都给出推导路径和实际验证方法。适合正在做嵌入式信号处理、传感器数据滤波、音频预处理的开发者也适合想真正搞懂IIR滤波器底层逻辑的学生和爱好者。读完你应该能自己推导任意截止频率下的一阶系数并且知道在什么场景下该选哪种配置。1. 从RC电路到差分方程一阶IIR的数学根基1.1 为什么一阶IIR本质上是RC电路的数字化身理解一阶IIR最直观的路径不是从Z变换开始而是从模拟RC电路出发。一个最简单的RC低通电路输入电压Vi通过电阻R给电容C充电输出电压Vo取自电容两端。这个电路的微分方程是C * dVo/dt (Vi - Vo) / R整理一下dVo/dt (Vi - Vo) / (R*C)令时间常数tau R*C则dVo/dt (Vi - Vo) / tau这个式子的物理含义很清晰输出电压的变化率正比于输入和输出之间的差值。差值越大电容充电越快差值越小输出越接近输入变化越慢。这就是低通滤波的本质——输出“跟不上”输入的快速变化所以高频成分被抑制。现在把这个微分方程离散化。用前向差分近似导数(Vo[n] - Vo[n-1]) / Ts (Vi[n] - Vo[n-1]) / tau其中Ts是采样周期。整理得到Vo[n] Vo[n-1] (Ts/tau) * (Vi[n] - Vo[n-1])再令alpha Ts / (tau Ts)可以写成更常见的形式Vo[n] alpha * Vi[n] (1 - alpha) * Vo[n-1]这就是一阶IIR低通滤波器最经典的差分方程。alpha是滤波系数取值范围在0到1之间。alpha越大输出越跟随输入截止频率越高alpha越小滤波越重截止频率越低。注意这里用的是前向差分实际上后向差分和双线性变换会得到略有差异的系数表达式但在Ts tau的条件下三者结果几乎一致。工程上直接用这个公式就够了。1.2 截止频率与alpha的精确换算关系上面推导出的alpha Ts / (tau Ts)是一个近似式。更精确的关系需要从模拟滤波器的截止频率反推。RC低通的-3dB截止频率是fc 1 / (2 * pi * R * C) 1 / (2 * pi * tau)所以tau 1 / (2 * pi * fc)。代入alpha的表达式alpha Ts / (1/(2*pi*fc) Ts)分子分母同乘2*pi*fcalpha (2 * pi * fc * Ts) / (1 2 * pi * fc * Ts)这个公式是工程上最常用的。我一般记成alpha 2 * pi * fc / (fs 2 * pi * fc)其中fs 1/Ts是采样频率。这个式子在fc远小于fs时非常准确当fc接近fs/2时会有偏差但一阶IIR本来就不适合做接近奈奎斯特频率的滤波所以实际使用中问题不大。举个实际数字假设采样率fs 1000Hz想要截止频率fc 10Hz那么alpha 2 * 3.14159 * 10 / (1000 2 * 3.14159 * 10) 62.83 / 1062.83 0.0591也就是说每次新采样值只贡献约6%的权重剩下94%来自历史输出。这个滤波效果相当重输出会非常平滑但响应也慢——阶跃响应大约需要1/(2*pi*10) ≈ 16ms的时间常数实际稳定需要5倍时间常数约80ms。如果fc 100Hz则alpha 628.3 / 1628.3 0.386新采样贡献38.6%响应明显快很多。这就是截止频率和滤波强度之间的权衡。1.3 差分方程三种等价写法的对照在实际代码和资料中一阶IIR低通有三种常见写法它们数学上等价但代码风格不同写法差分方程系数含义适用场景Alpha形式y[n] α·x[n] (1-α)·y[n-1]α 2πfc/(fs2πfc)快速手写嵌入式常用标准系数形式y[n] b0·x[n] b1·x[n-1] - a1·y[n-1]b0b1α/2, a1-(1-α)需要通用滤波器框架时增量形式y[n] y[n-1] α·(x[n] - y[n-1])同上定点数实现避免乘法溢出增量形式在定点MCU上特别有用因为x[n] - y[n-1]通常是个小数值乘以alpha后不会溢出而alpha * x[n]在输入满量程时可能接近累加器上限。我在STM32的ADC滤波里基本都用增量形式。标准系数形式看起来多了一个b1项实际上在一阶低通里b0 b1 alpha/2是一种常见的归一化写法目的是让直流增益精确为1。但直接用Alpha形式时直流增益天然就是1因为alpha (1-alpha) 1所以不需要额外处理。2. 一阶低通差分方程的代码落地与参数整定2.1 浮点实现与定点实现的取舍先看浮点版本这是最直观的typedef struct { float alpha; float prev_output; } LowPassIIR; float lowpass_update(LowPassIIR *f, float input) { f-prev_output f-alpha * input (1.0f - f-alpha) * f-prev_output; return f-prev_output; } void lowpass_init(LowPassIIR *f, float fc, float fs) { f-alpha 2.0f * 3.14159265f * fc / (fs 2.0f * 3.14159265f * fc); f-prev_output 0.0f; }这段代码在Cortex-M4带FPU的芯片上跑单次更新大约十几个时钟周期1kHz采样率下CPU占用可以忽略不计。但如果你的MCU没有硬件浮点比如Cortex-M0或者8位AVR浮点运算会非常慢这时候就得用定点。定点实现的关键是把alpha转成Q格式整数。假设用Q15格式1位符号15位小数alpha的范围是0到1对应整数0到32767。更新公式改成typedef struct { int32_t alpha_q15; int32_t prev_output_q15; } LowPassIIR_Fixed; int32_t lowpass_update_fixed(LowPassIIR_Fixed *f, int32_t input_q15) { int32_t diff input_q15 - f-prev_output_q15; int32_t delta (f-alpha_q15 * diff) 15; f-prev_output_q15 delta; return f-prev_output_q15; }这里用的是增量形式diff是两个Q15数相减范围在-65536到65535之间需要32位中间变量。alpha_q15 * diff最大约32767 * 65536 ≈ 2^31刚好在int32范围内不会溢出。右移15位后得到增量累加到历史输出上。实操心得定点实现时prev_output_q15建议用int32_t而不是int16_t因为累加过程中可能暂时超出Q15范围用int16会截断出错。输出给后续处理时再饱和到int16。2.2 截止频率选取的工程经验法则截止频率选多少这是实际项目里最常纠结的问题。我的经验是分三步走第一步确定你要保留的信号带宽。比如测温度信号变化本身就在0.1Hz以下那截止频率设1Hz都嫌高。测振动关注的是几十到几百Hz截止频率就得相应提高。第二步截止频率至少比采样率低一个数量级。fc fs/10是一阶IIR的舒适区此时alpha 0.386滤波效果明显且相位失真可控。如果fc fs/5alpha会超过0.6滤波效果很弱而且离散化误差开始显现。第三步考虑相位延迟。一阶IIR在截止频率处的相位延迟约45度对应的时间延迟是1/(8*fc)秒。对于闭环控制系统这个延迟会直接影响相位裕度。如果控制环路带宽是100Hz你加一个10Hz低通在100Hz处的额外相位延迟虽然不大但加上去可能就让系统临界震荡了。我一般会做一个简单的表格来辅助决策应用场景典型信号带宽建议fc建议fsalpha值温度采集0.1-1Hz1-5Hz10-50Hz0.1-0.4称重传感器1-10Hz5-20Hz50-200Hz0.1-0.4电机电流环100-500Hz500-2000Hz5-20kHz0.1-0.4音频预处理20Hz-20kHz视需求44.1-48kHz0.001-0.1振动监测10-1000Hz1-5kHz10-50kHz0.1-0.5注意音频预处理那一行alpha特别小是因为fc/fs比值很小。比如fc100Hz, fs48kHzalpha 628/48628 ≈ 0.013非常小滤波很重。2.3 阶跃响应验证怎么确认滤波器真的在工作写完代码别急着上系统先做个阶跃响应测试。给滤波器输入一个从0跳到1的阶跃信号观察输出。一阶IIR低通的阶跃响应是指数上升y[n] 1 - (1-alpha)^n理论上输出达到63.2%需要的时间是tau 1/(2*pi*fc)达到95%需要3*tau达到99%需要5*tau。比如fc10Hztau15.9ms。如果fs1000Hz采样周期1ms那么达到63.2%需要约16个采样点达到95%需要约48个点达到99%需要约80个点。你可以用Python快速验证import numpy as np import matplotlib.pyplot as plt fs 1000 fc 10 alpha 2*np.pi*fc / (fs 2*np.pi*fc) n 200 x np.ones(n) y np.zeros(n) for i in range(1, n): y[i] alpha * x[i] (1-alpha) * y[i-1] plt.plot(y) plt.axhline(0.632, colorr, linestyle--) plt.axhline(0.95, colorg, linestyle--) plt.show()如果实测阶跃响应和理论曲线对不上大概率是alpha算错了或者代码里把1-alpha写成了alpha。这种低级错误我见过不止一次调试时先检查这个。3. 一阶高通差分方程的推导与实现细节3.1 从低通到高通用输入减输出一阶高通和低通的关系非常优雅高通输出 输入 - 低通输出。这个关系不是近似而是精确成立的。为什么因为低通滤波器提取的是信号的低频成分那么原始信号减去低频成分剩下的自然就是高频成分。用公式表示y_hp[n] x[n] - y_lp[n]其中y_lp[n]是低通输出。把低通的差分方程代入y_hp[n] x[n] - (alpha * x[n] (1-alpha) * y_lp[n-1])而y_lp[n-1] x[n-1] - y_hp[n-1]代入整理y_hp[n] x[n] - alpha * x[n] - (1-alpha) * (x[n-1] - y_hp[n-1]) (1-alpha) * x[n] - (1-alpha) * x[n-1] (1-alpha) * y_hp[n-1] (1-alpha) * (x[n] - x[n-1] y_hp[n-1])这就是一阶高通的标准差分方程。令beta 1 - alpha则y_hp[n] beta * (x[n] - x[n-1] y_hp[n-1])beta越接近1即alpha越小高通截止频率越低beta越小截止频率越高。截止频率和beta的关系与低通一致只是beta 1 - alpha。3.2 高通实现中的直流漂移问题高通滤波器有一个低通没有的麻烦直流漂移。因为高通本质上是抑制直流分量的如果输入信号里有一个缓慢变化的直流偏置高通输出会围绕零点上下漂移但不会稳定在零。更具体地说如果输入突然加了一个直流偏置高通输出会先跳变然后指数衰减到零。衰减时间常数同样是tau 1/(2*pi*fc)。在衰减过程中输出会有一个“尾巴”这个尾巴可能被后续处理误认为是有效信号。我在做心电信号预处理时踩过这个坑。心电信号本身有基线漂移用高通滤掉基线后输出确实围绕零了但每次基线突变比如患者动了一下高通输出就会有一个大幅度的瞬态持续好几百毫秒。后来解决办法是在高通后面加一个软限幅或者用更复杂的基线恢复算法。代码实现上高通和低通结构几乎一样typedef struct { float beta; float prev_input; float prev_output; } HighPassIIR; float highpass_update(HighPassIIR *f, float input) { float output f-beta * (input - f-prev_input f-prev_output); f-prev_input input; f-prev_output output; return output; } void highpass_init(HighPassIIR *f, float fc, float fs) { float alpha 2.0f * 3.14159265f * fc / (fs 2.0f * 3.14159265f * fc); f-beta 1.0f - alpha; f-prev_input 0.0f; f-prev_output 0.0f; }注意高通需要保存prev_input和prev_output两个状态变量比低通多一个。这是差分方程里x[n-1]项带来的。3.3 高通滤波器的阶跃响应与低频抑制验证高通的阶跃响应和低通互补输入阶跃输出先跳到一个峰值然后指数衰减到零。峰值大小等于beta * 阶跃幅度衰减时间常数同样是tau。验证高通是否正常可以输入一个直流信号输出应该趋近于零。再输入一个方波输出应该是尖峰形状方波跳变时出尖峰平顶时衰减到零。实际测试时我一般用两个频率的正弦波叠加一个低于截止频率一个高于截止频率。高通输出应该几乎只剩高频那个。如果低频成分还很明显说明beta太大了截止频率设低了。注意一阶高通的阻带衰减只有-20dB/十倍频程也就是说截止频率以下一个十倍频程处衰减只有20dB。如果低频干扰比信号大40dB一阶高通是不够的需要二阶或更高阶。这是很多人容易忽略的地方——一阶滤波器的滚降很缓。4. 一阶带通差分方程的构建与中心频率调整4.1 低通串联高通最直接的带通方案一阶带通最自然的实现方式是低通和高通串联。信号先过高通滤掉低频再过低通滤掉高频。剩下的就是中间频带。串联的顺序理论上可以互换但实际有讲究。如果先低通后高通低通会先把高频噪声压掉再进高通时数值范围更小定点实现时精度更好。如果先高通后低通高通输出可能有大瞬态低通会把这个瞬态平滑掉但瞬态本身已经消耗了动态范围。我一般先高通后低通因为高通在前可以尽早去掉直流偏置避免直流分量在低通里积累。但如果你用的是浮点顺序影响不大。串联后的差分方程是两个方程的组合y1[n] beta * (x[n] - x[n-1] y1[n-1]) // 高通 y[n] alpha * y1[n] (1-alpha) * y[n-1] // 低通代码上就是两个滤波器结构体串联typedef struct { HighPassIIR hp; LowPassIIR lp; } BandPassIIR; float bandpass_update(BandPassIIR *f, float input) { float hp_out highpass_update(f-hp, input); return lowpass_update(f-lp, hp_out); } void bandpass_init(BandPassIIR *f, float fc_low, float fc_high, float fs) { highpass_init(f-hp, fc_low, fs); lowpass_init(f-lp, fc_high, fs); }这里fc_low是高通截止频率fc_high是低通截止频率。带通的有效通带在fc_low和fc_high之间。中心频率近似为sqrt(fc_low * fc_high)带宽为fc_high - fc_low。4.2 中心频率与带宽的独立调节技巧串联方案的缺点是中心频率和带宽不能独立调节。改变fc_low或fc_high会同时影响中心频率和带宽。但在很多应用里我们只关心中心频率带宽可以宽一些。比如做音频的1kHz提示音检测中心频率1kHz带宽可以设500Hz到2kHz这样fc_low500Hz, fc_high2kHz。如果中心频率要调到2kHz可以保持带宽比例设fc_low1kHz, fc_high4kHz。如果确实需要独立调节可以用状态变量滤波器结构但那已经不是一阶IIR了。一阶带通的Q值中心频率/带宽最大也就1左右做不了窄带。窄带需要二阶或更高阶。我整理了一个快速参考表中心频率带宽需求fc_lowfc_high适用场景1kHz宽(1oct)500Hz2kHz音频提示音检测10Hz窄(0.5oct)7Hz14Hz呼吸信号提取50Hz窄(0.2oct)45Hz55Hz工频干扰提取100Hz宽(2oct)25Hz400Hz振动包络分析注意50Hz那一行fc_low45, fc_high55带宽只有10Hz中心频率50HzQ5。一阶串联做不到这么窄实际需要二阶带通。一阶带通的Q值上限大约在1.5左右再窄就需要更高阶。4.3 带通滤波器的实测频响与常见偏差一阶带通的频响是低通和高通频响的乘积。在中心频率处如果fc_low和fc_high相距足够远增益接近1。但如果两者靠近中心频率处增益会小于1因为低通和高通在中心频率处都有衰减。具体来说在中心频率f0 sqrt(fc_low * fc_high)处高通增益为|H_hp(f0)| (f0/fc_low) / sqrt(1 (f0/fc_low)^2)低通增益为|H_lp(f0)| 1 / sqrt(1 (f0/fc_high)^2)总增益是两者乘积。如果fc_high/fc_low 4两个倍频程中心频率处总增益约0.8。如果比值是2一个倍频程总增益约0.6。如果比值是1.5总增益只有约0.45。这意味着如果你需要带通在中心频率处增益为1得额外乘一个补偿系数。补偿系数就是上面总增益的倒数。我在实际项目里一般会在带通后面加一个可调增益根据实测频响标定。实测频响的方法用信号发生器产生不同频率的正弦波记录带通输出幅度画成曲线。或者用白噪声激励做FFT分析。后者更快但需要能跑FFT的平台。实操心得一阶带通的过渡带很缓从通带到阻带需要大约两个倍频程才能衰减20dB。如果干扰频率离通带很近一阶带通基本没用直接上二阶或FIR。选滤波器阶数之前先算一下干扰频率和通带的距离距离小于一个倍频程就别考虑一阶了。5. 系数计算中的数值精度与实时调参5.1 浮点系数的精度陷阱alpha的计算公式里有个除法alpha 2*pi*fc / (fs 2*pi*fc)。当fc远小于fs时alpha是个很小的数。比如fc0.1Hz, fs1000Hzalpha ≈ 0.000628。用float24位有效尾数存储精度约1e-7alpha的相对误差约0.02%可以接受。但如果fc0.01Hzalpha ≈ 6.28e-5float的相对误差还是0.02%左右但绝对误差已经到1e-8量级在累加过程中可能被淹没。更麻烦的是1-alpha。当alpha很小时1-alpha在float里几乎等于1alpha的信息可能丢失。比如alpha1e-81-alpha在float里就是1.0因为float的epsilon约1.2e-7。这时候低通滤波器实际上变成了y[n] y[n-1]输出永远不变。解决办法是用double或者把差分方程改写成增量形式y[n] y[n-1] alpha*(x[n] - y[n-1])。增量形式里alpha直接乘以差值不会出现1-alpha的精度损失。我在所有嵌入式代码里都用增量形式就是为了避开这个坑。5.2 运行时动态调整截止频率的正确姿势有些应用需要在运行时改变截止频率比如自适应滤波。直接改alpha行不行行但要注意状态变量的连续性。如果只是改alphaprev_output保持不变滤波器会平滑过渡到新的截止频率不会有突变。这是正确的做法。但如果重新初始化滤波器把prev_output清零输出会突然跳到零产生一个瞬态。这个瞬态可能被后续处理误判。所以动态调参时只改系数不清状态。代码上void lowpass_set_fc(LowPassIIR *f, float fc, float fs) { f-alpha 2.0f * 3.14159265f * fc / (fs 2.0f * 3.14159265f * fc); // 不清零 prev_output }如果alpha变化很大比如从0.01跳到0.5输出会有一个快速跟踪过程但不会震荡。一阶IIR是无条件稳定的只要alpha在0到1之间怎么调都不会发散。注意alpha必须严格在0到1之间。alpha0意味着输出永远不变alpha1意味着输出等于输入没有滤波。实际使用中alpha建议限制在0.001到0.999之间避免极端值导致的数值问题。5.3 多通道滤波时的状态管理实际项目里经常要同时滤多路信号比如三轴加速度计。每路信号需要独立的滤波器状态不能共用。常见错误是定义一个全局滤波器结构体然后三路信号轮流调用。这样第二路会继承第一路的历史状态输出完全错误。正确做法是每路一个结构体实例LowPassIIR accel_x_filter; LowPassIIR accel_y_filter; LowPassIIR accel_z_filter; lowpass_init(accel_x_filter, 10.0f, 1000.0f); lowpass_init(accel_y_filter, 10.0f, 1000.0f); lowpass_init(accel_z_filter, 10.0f, 1000.0f);如果通道数很多可以用数组#define NUM_CHANNELS 8 LowPassIIR filters[NUM_CHANNELS]; for (int i 0; i NUM_CHANNELS; i) { lowpass_init(filters[i], 10.0f, 1000.0f); }内存占用上每个低通滤波器两个floatalpha和prev_output8通道也就64字节可以忽略。高通多一个prev_input带通是两个滤波器串联内存翻倍。在RAM紧张的MCU上如果通道数很多可以考虑用定点压缩存储。6. 从理论到实战三个典型场景的完整配置6.1 场景一称重传感器的高频噪声抑制称重传感器输出的是缓慢变化的重量信号但机械振动和电气噪声会叠加高频干扰。典型信号带宽不到10Hz噪声可能到几百Hz。配置fs 100Hz每10ms采样一次fc 5Hz。计算alphaalpha 2*pi*5 / (100 2*pi*5) 31.4 / 131.4 0.239用增量形式实现float weight_filter(float input) { static float prev 0.0f; float alpha 0.239f; prev prev alpha * (input - prev); return prev; }阶跃响应时间常数tau 1/(2*pi*5) ≈ 32ms稳定时间约160ms。对于称重应用这个响应速度足够因为重量本身变化就慢。实测效果原始信号峰峰值噪声约50个ADC码滤波后降到5个码以内相当于10倍改善。代价是响应延迟160ms但称重不需要快速响应。实操心得称重传感器滤波后建议再做一次“去皮”处理把空载时的滤波输出存为零点后续读数减去零点。因为一阶IIR的直流增益是1但初始状态为零上电后需要一段时间才能跟踪到真实零点。去皮可以消除这个启动瞬态。6.2 场景二心电信号基线漂移去除心电信号带宽约0.05Hz到100Hz但基线漂移通常在0.5Hz以下。需要高通滤掉基线同时保留0.5Hz以上的心电成分。配置fs 500Hz高通fc 0.5Hz。计算betaalpha 2*pi*0.5 / (500 2*pi*0.5) 3.14 / 503.14 0.00624 beta 1 - alpha 0.99376高通实现float ecg_highpass(float input) { static float prev_input 0.0f; static float prev_output 0.0f; float beta 0.99376f; float output beta * (input - prev_input prev_output); prev_input input; prev_output output; return output; }这里beta非常接近1高通截止频率很低。阶跃响应的衰减时间常数tau 1/(2*pi*0.5) ≈ 318ms意味着基线突变后需要约1.6秒才能衰减到5%以下。对于心电监测这个恢复时间可以接受因为基线突变不频繁。实测效果基线漂移从±200个ADC码降到±10个码以内心电波形基本不变形。但每次患者移动导致的基线突变输出会有约1秒的瞬态需要后续算法标记为无效段。6.3 场景三音频1kHz提示音检测的带通预处理要从麦克风信号里检测1kHz提示音先用带通滤掉其他频率再做幅度检测。配置fs 16kHz带通fc_low 500Hzfc_high 2kHz。计算alpha_lp 2*pi*2000 / (16000 2*pi*2000) 12566 / 28566 0.440 alpha_hp 2*pi*500 / (16000 2*pi*500) 3141.6 / 19141.6 0.164 beta_hp 1 - 0.164 0.836带通实现typedef struct { float hp_beta; float hp_prev_in; float hp_prev_out; float lp_alpha; float lp_prev_out; } AudioBandPass; float audio_bandpass(AudioBandPass *f, float input) { // 高通 float hp_out f-hp_beta * (input - f-hp_prev_in f-hp_prev_out); f-hp_prev_in input; f-hp_prev_out hp_out; // 低通 f-lp_prev_out f-lp_prev_out f-lp_alpha * (hp_out - f-lp_prev_out); return f-lp_prev_out; }中心频率f0 sqrt(500*2000) 1000Hz正好是目标频率。中心频率处增益约0.8因为fc_high/fc_low 4后续检测阈值需要相应调整。实测效果1kHz提示音通过500Hz以下和2kHz以上的干扰被明显抑制。但如果是800Hz的干扰衰减只有约6dB可能误检。所以带通后面还需要做频率确认比如过零率检测确保确实是1kHz附近。注意音频应用里一阶带通的滚降太缓通常需要二阶或更高阶。这里用一阶只是为了演示配置方法实际产品里建议用二阶巴特沃斯带通滚降-40dB/十倍频程选择性好得多。7. 一阶IIR的边界与替代方案选择7.1 什么时候一阶IIR不够用一阶IIR的滚降是-20dB/十倍频程这是它最大的局限。如果干扰频率离通带只有一个倍频程一阶滤波只能提供约6dB衰减基本没用。这种情况下必须上二阶。判断标准很简单算一下干扰频率和截止频率的比值。如果比值小于4两个倍频程一阶滤波效果有限比值在4到10之间一阶可用但效果一般比值大于10一阶效果很好。另一个局限是过渡带形状。一阶IIR的过渡带是对数曲线没有纹波但也不够陡。如果应用要求通带平坦、阻带陡峭一阶做不到需要巴特沃斯或切比雪夫。还有群延迟。一阶IIR的群延迟在截止频率附近变化很大对于需要保持波形形状的应用比如心电、音频可能造成波形失真。这种情况下用线性相位FIR更好代价是更多计算量。7.2 一阶IIR与移动平均、中值滤波的对比移动平均滤波是另一种常见的低通方法实现简单但频响有旁瓣阻带衰减只有-13dB左右而且旁瓣会泄漏高频噪声。一阶IIR的频响是单调下降的没有旁瓣阻带衰减虽然也只有-20dB/十倍频程但至少是单调的。中值滤波对脉冲噪声特别有效但计算量大需要排序而且对高斯噪声效果不如IIR。实际项目里经常是IIR加中值组合先中值去掉脉冲再IIR平滑。我整理了一个对比表滤波方法计算量脉冲噪声抑制高斯噪声抑制相位失真适用场景一阶IIR极低差好中等实时嵌入式移动平均低差中等线性简单平滑中值滤波中高极好差非线性脉冲噪声二阶IIR低差很好较大需要陡滚降FIR高中等好可线性波形保真一阶IIR在计算量和滤波效果之间取得了很好的平衡这是它在嵌入式领域长盛不衰的原因。7.3 从一阶到二阶的平滑升级路径如果一阶不够用升级到二阶是最自然的路径。二阶IIR的差分方程是y[n] b0*x[n] b1*x[n-1] b2*x[n-2] - a1*y[n-1] - a2*y[n-2]系数计算比一阶复杂但有很多现成工具。我一般用Python的scipyfrom scipy import signal b, a signal.butter(2, fc/(fs/2), low) print(b, a)得到系数后直接填到代码里。二阶的滚降是-40dB/十倍频程比一阶好一倍。代价是多两个状态变量和几个乘法在Cortex-M4上仍然很轻松。从一阶升级到二阶时注意截止频率的定义。一阶的fc是-3dB点二阶巴特沃斯的fc也是-3dB点但二阶的过渡带更陡。如果原来一阶用fc10Hz升级到二阶时fc可以设得更高一些比如15Hz因为二阶在10Hz处的衰减已经比一阶大了。实操心得升级到二阶后阶跃响应会出现轻微过冲巴特沃斯约4%这是正常的。如果过冲不可接受可以用贝塞尔响应相位更线性但滚降更缓。选哪种响应取决于应用对相位失真和滚降的要求。8. 调试工具与验证方法8.1 用Python快速验证系数和频响写代码之前先用Python验证系数和频响可以省去大量调试时间。下面这段代码画出一阶低通、高通、带通的频响import numpy as np import matplotlib.pyplot as plt fs 1000 fc_lp 100 fc_hp 10 alpha_lp 2*np.pi*fc_lp / (fs 2*np.pi*fc_lp) alpha_hp 2*np.pi*fc_hp / (fs 2*np.pi*fc_hp) beta_hp 1 - alpha_hp f np.logspace(0, np.log10(fs/2), 500) w 2*np.pi*f/fs # 低通频响 H_lp alpha_lp / (1 - (1-alpha_lp)*np.exp(-1j*w)) # 高通频响 H_hp beta_hp * (1 - np.exp(-1j*w)) / (1 - beta_hp*np.exp(-1j*w)) # 带通频响 H_bp H_lp * H_hp plt.semilogx(f, 20*np.log10(np.abs(H_lp)), labelLowpass) plt.semilogx(f, 20*np.log10(np.abs(H_hp)), labelHighpass) plt.semilogx(f, 20*np.log10(np.abs(H_bp)), labelBandpass) plt.axhline(-3, colork, linestyle--) plt.legend() plt.grid() plt.show()这段代码可以直接跑看到三条曲线。低通在100Hz处-3dB高通在10Hz处-3dB带通在10Hz到100Hz之间平坦。如果曲线形状不对说明系数算错了。8.2 嵌入式端的实时调试技巧在MCU上调试滤波器最直接的方法是把输入输出通过串口或DAC输出。如果MCU有DAC把滤波前后的信号分别输出到两个DAC通道用示波器看非常直观。如果没有DAC可以用PWM加RC滤波模拟DAC或者通过串口把数据传到电脑上画图。串口传输时注意采样率不要太高115200波特率下每秒传1000个点已经接近极限。我常用的方法是在RAM里开一个环形缓冲区存最近N个输入输出点然后通过调试器读取。这样不影响实时性又能看到波形。N取256或512就够看一个阶跃响应了。还有一个技巧用GPIO翻转来测量滤波器执行时间。在滤波函数入口拉高GPIO出口拉低示波器看脉冲宽度。一阶IIR在72MHz的Cortex-M3上大约1-2微秒在带FPU的M4上不到1微秒。8.3 常见异常现象与排查对照现象可能原因排查方法输出不变alpha0或prev_output未更新检查alpha计算打印中间值输出等于输入alpha1或beta0检查截止频率是否大于fs/2输出震荡alpha超出0-1范围检查fc和fs的数值输出缓慢漂移高通直流漂移正常现象加后续处理多通道串扰共用滤波器状态每通道独立结构体定点输出溢出中间变量位宽不够用int32中间变量阶跃响应不对系数公式用错用Python验证频响这张表是我多年调试经验的总结遇到问题先对照排查大部分情况能快速定位。9. 个人实操体会与几个容易忽略的细节一阶IIR看起来简单但真正用好需要理解它的边界。我最大的体会是不要试图用一阶IIR做它做不到的事。需要陡滚降就上二阶需要线性相位就上FIR需要窄带就上二阶带通。一阶IIR的定位是“轻量级实时平滑”在这个定位内它非常出色超出这个定位就是给自己找麻烦。另一个体会是关于初始状态。很多人在初始化时把prev_output设为零如果输入信号本身有直流偏置滤波器需要很长时间才能跟踪到真实值。更好的做法是用第一个采样值初始化prev_output这样启动瞬态最小。代码上就是void lowpass_init_with_first_sample(LowPassIIR *f, float first_sample) { f-prev_output first_sample; }还有一个细节是采样率变化。如果系统支持多种采样率alpha必须随采样率重新计算。我见过一个项目采样率从1kHz切到10kHz但alpha没更新结果截止频率实际变成了原来的10倍滤波完全失效。这种问题很隐蔽因为代码不报错只是效果不对。最后分享一个快速估算技巧一阶IIR的建立时间到99%约等于5/(2*pi*fc) ≈ 0.8/fc秒。比如fc10Hz建立时间约80ms。这个估算在选截止频率时非常有用可以快速判断响应速度是否满足要求。一阶IIR的差分方程就这些内容低通、高通、带通的配置方法都已经展开。实际使用时先用Python验证频响再移植到嵌入式最后用阶跃响应确认。这套流程走下来基本不会出问题。
返回列表