从传递函数到频域滤波：CModule模块化实现嵌入式IIR滤波器

1. 从传递函数到频域滤波的整体设计思路很多做嵌入式信号处理的朋友第一次接触“传递函数”这个词大概率是在控制理论或者信号与系统课上一堆 s 域、z 域的公式推来推去考完试就还给老师了。但真正到了项目里比如你要在单片机上做一个 IIR 滤波器、要分析一个电机控制环路的稳定性、或者要把一段传感器数据做频域降噪你会发现——传递函数这套东西其实是绕不开的。它不是一个学术符号而是一张“系统行为的地图”。我这次想聊的就是把这条链路完整走一遍从连续域的 Laplace 变换出发理解传递函数到底在描述什么再到离散域的 Z 变换搞清楚数字系统里它长什么样最后落到 CModule 这个模块化实现上把频域滤波真正跑起来。整条链路的核心目标只有一个让你能从一个数学表达式推导出一段可以烧进芯片、实时运行的滤波代码。为什么非要走这条链路而不是直接抄一个滤波器系数表因为抄来的系数你不知道它为什么长这样一旦采样率变了、截止频率要调、或者相位特性不满足要求你就抓瞎了。而理解了传递函数到频域滤波的映射关系你就能自己算系数、自己判断稳定性、自己决定用哪种结构实现。这是从“会用”到“能改”的分水岭。适合谁来参考这篇内容如果你写过 C 语言、用过至少一款 MCU 或 DSP、对傅里叶变换有一点点印象那就可以往下看。我会尽量把数学部分用工程语言翻译一遍不堆公式重点讲“这个式子在实际代码里对应哪一行”。CModule 在这里扮演的角色是把滤波器抽象成一个可配置、可复用的模块而不是每次都在主循环里手写差分方程。整个设计思路可以概括成三步走连续域建模 → 离散域映射 → 模块化实现。第一步用 Laplace 传递函数描述“我想要一个什么样的频率响应”第二步用 Z 变换和双线性变换等方法把这个连续传递函数变成离散的差分方程第三步用 CModule 把差分方程封装成 init、process、reset 这样的标准接口方便在工程里反复调用。下面我逐层拆开讲。2. Laplace 变换与传递函数连续域到底在描述什么2.1 传递函数的本质是一个“频率响应说明书”先抛开数学定义。你可以把传递函数 H(s) 想象成一台机器的“性格说明书”你给我一个频率的正弦输入我告诉你输出会变成多大幅度、偏移多少相位。s 这个复变量里实部对应衰减/增长虚部对应振荡频率。当我们在频域滤波的语境下讨论时通常关心的是 s jω 这条虚轴上的取值也就是系统对不同频率正弦信号的稳态响应。举个最朴素的例子一个 RC 低通滤波器它的传递函数是 H(s) 1 / (1 sRC)。这个式子告诉我们什么当频率很低ω 趋近 0分母趋近 1输出约等于输入增益为 1当频率很高分母的 sRC 项占主导增益随频率下降这就是“低通”的由来。截止频率出现在 ω 1/RC 处此时增益降到约 0.707也就是 -3dB 点。我之所以强调这个是因为很多人在做数字滤波时脑子里没有这个连续域的“原型”。你直接看一个离散的差分方程 y[n] ax[n] bx[n-1] - c*y[n-1]很难直观说出它的截止频率在哪。但如果你知道它是由某个 H(s) 映射来的一切就清晰了。传递函数是设计的起点差分方程只是实现的手段。2.2 极点零点决定了系统的“脾气”传递函数写成因式分解形式后分子多项式的根叫零点分母多项式的根叫极点。这两个东西决定了系统的全部行为。极点在左半平面系统稳定极点越靠近虚轴对应频率附近的响应越尖锐零点可以用来“压掉”某些频率。做频域滤波时我们其实就是在“摆放”极点和零点的位置让幅频响应长成我们想要的样子。这里有个实操中很容易踩的坑连续域稳定的系统离散化之后不一定稳定。因为映射方法比如双线性变换、冲激响应不变法会改变极点位置。我见过有人照着模拟电路算了个传递函数直接离散化烧进芯片结果输出发散查了半天才发现是离散化后极点跑到单位圆外了。所以从 Laplace 到 Z 变换这一步必须验证离散极点的模是否小于 1。2.3 为什么工程上偏爱二阶节Biquad而不是高阶直接实现理论上你可以设计一个 10 阶的传递函数然后直接写一个 10 阶差分方程。但实际工程里几乎没人这么干原因是数值精度。高阶多项式的系数对量化误差极其敏感定点数运算下很容易不稳定。标准做法是把高阶传递函数分解成多个二阶节Biquad级联每个二阶节单独实现这样每个环节的极点零点都容易控制数值特性也好得多。CModule 的设计正是基于这个思路一个 CModule 实例对应一个二阶节多个实例级联就构成高阶滤波器。这样模块的粒度刚刚好既不会太细碎也不会因为阶数太高而失控。后面讲实现时我会具体说这个结构怎么定义。3. Z 变换与离散映射从 s 域到 z 域的关键一步3.1 Z 变换在工程视角下就是“离散信号的 Laplace”连续信号做 Laplace 变换得到 s 域表达式离散序列做 Z 变换得到 z 域表达式两者在结构上高度相似。z 域里z^-1 代表一个采样周期的延迟这个物理意义非常重要。一个差分方程 y[n] b0x[n] b1x[n-1] - a1y[n-1]对应的 Z 域传递函数就是 H(z) (b0 b1z^-1) / (1 a1*z^-1)。你看系数直接对应上了这就是为什么差分方程和 Z 域传递函数可以互相翻译。理解这一点之后频域滤波的实现就变成了“求系数”的问题。你不需要在代码里做任何变换运算只需要把预先算好的 b0、b1、a1 填进 CModule然后每个采样点执行一次乘加运算即可。实时性就是这么保证的。3.2 双线性变换最常用的 s 到 z 映射方法从连续传递函数 H(s) 到离散传递函数 H(z)有几种映射方法工程上最常用的是双线性变换Bilinear Transform。它的映射关系是 s (2/T) * (1 - z^-1) / (1 z^-1)其中 T 是采样周期。这个变换的好处是把整个 s 左半平面映射到 z 平面的单位圆内稳定性得以保持而且没有频率混叠问题代价是频率轴有非线性压缩也就是所谓的“频率翘曲”。频率翘曲这件事必须提前补偿。双线性变换把连续频率 ω 映射到离散频率 Ω 的关系是 Ω 2arctan(ωT/2)/T当频率接近奈奎斯特频率时这种非线性越来越明显。所以如果你要设计一个截止频率为 1kHz 的滤波器采样率 8kHz不能直接把 ω 2π1000 代进去算而要先用预畸变公式把目标频率修正一下。我一般会在设计脚本里加一行预畸变计算避免截止频率偏移。3.3 从传递函数到差分方程的完整推导示例拿一个一阶低通举例。连续传递函数 H(s) ωc / (s ωc)其中 ωc 是截止角频率。用双线性变换代入 s (2/T)*(1-z^-1)/(1z^-1)经过整理可以得到H(z) (ωcT / (2 ωcT)) * (1 z^-1) / (1 ((ωcT - 2)/(ωcT 2))*z^-1)令 K ωcT / (2 ωcT)a1 (ωcT - 2)/(ωcT 2)那么差分方程就是y[n] Kx[n] Kx[n-1] - a1*y[n-1]你看最终落到代码里就是三个系数K、K、a1。CModule 要做的就是把这些系数存起来然后每个采样点算一次这个式子。这就是频域滤波在嵌入式里的真实面貌——没有玄学就是乘加。4. CModule 模块化实现让滤波器变成可复用积木4.1 为什么选择模块化而不是裸写差分方程我早期做项目时滤波器代码都是直接写在主循环里的几个全局变量存历史输入输出几行乘加完事。这样做在小项目里没问题但一旦系统里有三四个滤波器级联或者需要在不同采样率下切换代码就会变得一团糟。历史状态散落在各处复位逻辑不统一调试时根本分不清哪个变量属于哪个滤波器。CModule 的思路是把每个滤波器封装成一个结构体加一组函数。结构体里存系数和历史状态函数提供初始化、处理单点、批量处理、复位这几个标准接口。这样每个滤波器就是一个独立的对象可以创建多个实例互不干扰。级联时只需要把前一个的输出喂给后一个的输入逻辑非常清晰。4.2 二阶节结构体的字段设计与内存布局一个标准的 Biquad CModule 结构体通常包含这些字段系数 b0、b1、b2、a1、a2a0 归一化为 1状态变量 x1、x2、y1、y2分别代表前两个输入和前两个输出以及可选的增益字段。用 float 还是定点数取决于平台Cortex-M4F 以上带 FPU 的芯片直接用 float 就很方便没有 FPU 的话用 Q15 定点数更划算。内存布局上有个细节值得注意把系数和状态分开存放系数用 const 修饰放在 Flash 里状态放在 RAM 里。这样多个实例可以共享同一组系数比如立体声左右声道用同样的滤波参数只需要各自维护状态即可省 RAM。我在做音频处理时就用了这个技巧两个声道共一套系数RAM 占用直接减半。4.3 直接I型、直接II型与转置直接II型的取舍Biquad 有三种常见实现结构。直接I型最直观但需要 4 个状态变量直接II型把状态变量减到 2 个省内存但中间变量的数值范围可能较大转置直接II型在定点运算下溢出风险最小是工程上最推荐的结构。我一般默认用转置直接II型它的差分方程写法是y[n] b0x[n] s1 s1 b1x[n] - a1y[n] s2 s2 b2x[n] - a2*y[n]其中 s1、s2 是两个中间状态。这种结构下每个状态变量的动态范围都比较可控定点实现时不容易溢出。CModule 里我会把这个计算过程写成一个内联函数减少函数调用开销。5. 频域滤波的实操流程与参数计算5.1 滤波器设计的完整步骤清单从需求到可运行代码我通常按这个顺序走明确滤波需求低通、高通、带通还是带阻截止频率多少通带纹波和阻带衰减要求多少。选择原型巴特沃斯最平坦切比雪夫过渡带更陡但有纹波椭圆滤波器最陡但相位非线性严重。音频应用偏爱巴特沃斯因为相位失真小。确定阶数根据过渡带宽度和衰减要求查表或计算一般 2 到 4 阶够用。计算连续传递函数的极点和零点。用双线性变换离散化注意预畸变。分解成二阶节计算每个节的系数。把系数填入 CModule验证频响。这个流程里第 4 到第 6 步通常用脚本完成Python 的 scipy.signal 几行就能搞定。但理解每一步在做什么比会调库重要得多。5.2 采样率、截止频率与系数精度的关系采样率的选择直接影响滤波器的数值特性。采样率越高截止频率相对采样率越低双线性变换的频率翘曲越小但系数之间的差异也会变小定点运算时精度损失更明显。我一般建议截止频率不要低于采样率的 1/100否则系数会非常接近量化误差可能导致实际频响偏离设计值。系数精度方面如果用 float基本不用担心如果用 Q15 定点系数要缩放到 -32768 到 32767 范围内且要注意乘加后的移位处理。我踩过的坑是两个 Q15 数相乘得到 Q30累加几次后可能溢出 32 位必须及时右移。CModule 里我会在每次乘加后做饱和处理宁可轻微失真也不能让数据回绕。5.3 用 Python 脚本生成 CModule 系数表实际项目中我习惯写一个 Python 脚本输入滤波参数输出一个 C 头文件里面是系数数组。这样每次调参只需要改脚本重新生成不用手动改代码。脚本核心逻辑大概是import numpy as np from scipy import signal fs 8000.0 fc 1000.0 order 4 # 预畸变 wc 2 * np.pi * fc wc_pre 2 * fs * np.tan(wc / (2 * fs)) # 设计模拟原型 b_analog, a_analog signal.butter(order, wc_pre, analogTrue) # 双线性变换离散化 b, a signal.bilinear(b_analog, a_analog, fs) # 分解成二阶节 sos signal.tf2sos(b, a) for section in sos: print(fb0{section[0]:.8f}, b1{section[1]:.8f}, b2{section[2]:.8f}) print(fa0{section[3]:.8f}, a1{section[4]:.8f}, a2{section[5]:.8f})生成的系数直接填进 CModule 的初始化函数即可。注意 a0 通常归一化为 1如果脚本输出的 a0 不是 1要先把所有系数除以 a0。6. 常见问题与排查技巧实录6.1 输出发散或振荡先查极点和状态复位滤波器输出发散是最常见的问题原因通常有三类一是离散化后极点跑到单位圆外二是定点运算溢出三是状态变量没有正确初始化。排查顺序我建议先看极点用脚本算一下离散传递函数的极点模值大于 1 就是设计问题如果极点没问题再看是不是溢出把中间变量打印出来看有没有异常大的值最后检查复位逻辑确保每次启动或切换参数时历史状态清零。6.2 截止频率与设计值不符多半是预畸变没做很多人发现实际测出来的 -3dB 点和设计值差不少尤其是截止频率接近奈奎斯特频率时。这几乎可以肯定是双线性变换的频率翘曲导致的。解决办法就是在设计时做预畸变把目标频率用 tan 函数修正一下。我一般会在脚本里直接对比修正前后的频响确认误差在可接受范围内。6.3 级联滤波器增益异常注意每个节的增益分配多个二阶节级联时如果每个节都按单位增益设计整体增益可能在某些频率上叠加出很大的值。标准做法是把总增益分配到各个节或者单独用一个增益系数统一调整。我在做四阶低通时就遇到过通带增益偏高 6dB 的情况后来发现是两个节的 Q 值都设得比较高峰值叠加了。把 Q 值重新分配后就正常了。6.4 定点运算下的量化噪声系数和状态都要考虑定点实现时系数被量化会改变频响状态变量被量化会引入噪声。我的经验是系数用 Q15 一般够用但状态变量最好用 Q30 或更高精度避免累加误差。如果平台支持用 32 位定点做状态累加最后再饱和到 Q15 输出效果比全程 Q15 好很多。问题现象可能原因排查方法解决措施输出发散极点不稳定计算离散极点模值重新设计或换映射方法输出发散定点溢出打印中间变量加饱和处理或提高位宽截止频率偏移未做预畸变对比设计值与实测值加入预畸变计算通带增益异常级联增益叠加逐节测频响重新分配增益和 Q 值输出有噪声量化误差对比浮点与定点结果提高状态变量精度7. 我在实际项目中的几点体会做这类频域滤波项目最大的感受是数学推导和代码实现之间隔着一层“工程直觉”。课本教你的是理想情况下的公式但实际芯片有量化误差、有中断延迟、有内存限制。我现在的习惯是任何滤波器设计先在 Python 里跑通频响再用浮点 C 代码验证最后才转定点。每一步都对比结果确保偏差在预期内。另一个体会是关于 CModule 的接口设计。我早期喜欢把系数直接写死在处理函数里后来发现调试时改一个参数要重新编译整个工程效率太低。现在我会把系数做成运行时可配置的通过一个 init 函数传入这样调参时只需要改配置数组不用动核心逻辑。这个改动看似小但实际调试效率提升非常明显。最后分享一个小技巧验证滤波器是否正常工作最直接的方法不是看波形而是喂一个扫频信号进去把输出幅值记录下来画成曲线。这比看时域波形直观得多一眼就能看出通带、阻带和过渡带是否符合预期。我在 CModule 里专门留了一个测试模式可以旁路滤波直接输出方便对比。