AI时代软件工程变革：从敏捷开发到智能协作
2026/9/17 22:17:01
简介：这份 PDF 资料面向数字信号处理、嵌入式开发与语音通信方向的工程师及学生，围绕 MATLAB FDAtool 设计 IIR 滤波器并把参数导出为 C 语言文件这一实际问题展开。内容从角频率与采样频率的换算关系切入，说明通带、阻带截止频率（边沿频率）的设定方法，并对比 FIR 与 IIR 在阶数、计算量和相位特性上的差异；随后以一个采样率 8kHz、通带 80–3200Hz、用于滤除 50Hz 工频干扰的带通滤波器为例，演示在 FDAtool 中设定指标、由工具自动计算并得到 36 阶（18 个二阶节）滤波器的完整过程，还涉及自定义阶数、增益项对精度与稳定性的作用，以及把系数转换为单段形式、导出可直接在 C 代码中复用的头文件等关键环节。资源为单个 PDF 文件，压缩包约 1.08MB，结构紧凑、便于随查随用。目前已有 1423 人学习，适合希望把 MATLAB 滤波器设计成果迁移到嵌入式或其他编程环境的读者参考。
1. FDAtool 生成的系数，为什么不能直接手敲进 C 文件
做嵌入式音频、电机电流环或者传感器抗混叠的人，多半踩过同一个坑：在 MATLAB 里画出来的频响漂漂亮亮，系数抄进 C 代码一跑，实测曲线要么通带塌下去一块，要么高频直接自激。问题通常不在双线性变换推错了，而在两件事上。第一，高阶 IIR 的传递函数直接用 b、a 多项式实现，极点对系数的敏感度随阶数指数级上升，浮点截断一下极点就跑到单位圆外；第二，FDAtool 导出的系数是按 a0 归一化过的，手抄时忘了把 a0 除掉，整个滤波器的增益和极点位置都会偏。
FDAtool（较新版本里入口叫 Filter Designer）真正的价值不是画曲线，而是它能把设计好的滤波器拆成级联二阶节，导出成 C 头文件，让 MATLAB 里的仿真和 DSP 上跑的代码尽量对齐。这条链路适合两类人：一类是要把 MATLAB 做离散时间系统的结果搬到单片机上的人，另一类是需要在定点 DSP 上控制阶数的工程实现者。下面按"定参、选结构、导文件、写 C、复核"的顺序把这条路走一遍。
2. 用 FDAtool 设计 IIR：原型怎么选、参数怎么填
2.1 命令行先定参，GUI 后复核
GUI 适合调参，不适合复现。同一组参数，今天点出来一个结果，明天换个版本可能默认值就变了，所以我一般先用命令行把阶数和系数算出来，再打开filterDesigner用fvtool复核曲线。这样可以保证参数落在版本库里，随时能重跑。五个必填量是采样率fs、通带边缘fpass、阻带边缘fstop、通带波纹rp、阻带衰减rs。其中fpass/(fs/2)必须严格落在 0 到 1 之间，写成 1 或者大于 1，FDAtool 会直接提示归一化频率非法。
fs = 48000; % 采样率 Hz fpass = 3000; % 通带边缘 Hz fstop = 4000; % 阻带边缘 Hz rp = 0.5; % 通带波纹 dB rs = 60; % 阻带衰减 dB Wp = fpass/(fs/2); % 归一化到 Nyquist 频率 Ws = fstop/(fs/2); [N, Wn] = ellipord(Wp, Ws, rp, rs); % 求满足指标的最小阶数 [b, a] = ellip(N, rp, rs, Wn); % 直接型传递函数系数 [sos, g] = tf2sos(b, a); % 拆成二阶节，g 为总增益 fvtool(b, a, 'Fs', fs); % 看幅频、相频、群延迟 disp(['所需最小阶数 N = ', num2str(N)]);这段代码的逻辑是：ellipord先根据四个边界指标反推最小阶数，避免手工试阶数；ellip按这个阶数生成分子分母多项式；tf2sos再把它拆成若干二阶节并抽出总增益g。参数上最需要留意的是fpass和fstop的间距，两者越靠近，ellipord返回的N越大，乘加次数和定点溢出风险一起上涨。一般过渡带宽度小于通带边缘的 20% 时，就该认真考虑降采样或者换结构了。
2.2 四种 IIR 原型的取舍
同样一组频带指标，换原型得到的阶数和相位特性差别很大。椭圆滤波器能把阶数压到最低，但代价是通带和阻带同时等波纹，相位非线性最严重；巴特沃斯通带最平坦，阶数却常常高出一截。选型时先看你的约束是算力、延迟还是相位线性度。
| 原型 | 通带特性 | 阻带特性 | 过渡带陡度 | 相位非线性 | 典型场景 |
|---|---|---|---|---|---|
| Butterworth | 最平坦 | 单调下降 | 最缓 | 中等 | 生物信号、抗混叠 |
| Chebyshev I | 等波纹 | 单调下降 | 较陡 | 较差 | 窄过渡带音频均衡 |
| Chebyshev II | 单调下降 | 等波纹 | 较陡 | 较好 | 阻带抑制要求高 |
| Elliptic | 等波纹 | 等波纹 | 最陡 | 最差 | 阶数受限的定点 DSP |
实战里的经验值是：MCU 上没有硬件浮点、又必须把阶数压到 6 阶以内，优先椭圆；对相位失真敏感、能接受十几阶的，用巴特沃斯。Chebyshev II 常被忽略，它在"阻带必须干净、通带允许一点起伏"的场合其实很划算。选型不要只看曲线好不好看，要看阶数带来的乘加次数能不能塞进你的采样周期。
2.3 为什么要把高阶传递函数拆成 SOS
IIR 的极点在 z 平面上越靠近单位圆，多项式系数的一点点扰动就越容易被放大。八阶直接型的 a 系数里，某一位的尾数误差，可能让一对共轭极点在量化后跑到单位圆外，滤波器直接从"能滤"变成"会炸"。级联二阶节把一个大多项式拆成若干个二阶子系统的乘积，每个子系统的极点只由它自己的两三个系数决定，敏感度被摊薄，定点化后稳定性明显好过直接型。
拆成 SOS 之后还有两个附带好处：一是每个二阶节可以单独做饱和和舍入策略，二是可以用流水线或者 SIMD 并行处理多个节。需要注意tf2sos的默认配对策略并不总是最优，它在节顺序上只保证数值上的合理性，不保证定点实现里最优。如果定点噪声偏大，可以把tf2sos的第二个输出参数换成'down'或'up'调整配对顺序，再实测噪声底。
3. 把 FDAtool 的系数导出成可编译的 C 语言文件
3.1 Export to Workspace 与 Generate C header 的差别
FDAtool 的 File 菜单里有几条不同的出口，选错了后面要么格式对不上，要么没法进版本管理。Export to Workspace 是把设计对象、系数、SOS、增益导出成 MATLAB 变量，适合脚本继续处理；Generate C header 直接吐一个单精度浮点数组，适合临时验证，但它和 GUI 当前设置强绑定，别人拿到头文件也还原不出你当初的设计条件。
| 导出项 | 输出内容 | 适用阶段 | 需要留意 |
|---|---|---|---|
| Export to Workspace | 对象、b/a、SOS、g | 脚本继续加工 | 必须选 SOS 形式，别选 b/a |
| Generate C header | float 数组头文件 | 快速粘贴验证 | 无法追溯设计参数 |
| Generate MATLAB code | 可重跑的 M 脚本 | 进版本库 | 不是 C，需二次转换 |
| 自己 fopen/fprintf 生成 | 自定义 .h/.c | 量产工程 | 自由度最高，要自己定格式 |
我几乎不直接用 Generate C header，因为它给出的数组是单节的 b/a 形式，不是一个完整的二阶节矩阵，落到代码里还要自己重组。更稳妥的做法是把sos和g拿到手，自己用文件读写生成需要的 .h 和 .c，格式完全按目标工程来。
3.2 用 fopen 和 fprintf 自动生成 .h 与 .c
MATLAB 写 C 文件本质上就是一段普通的 C 语言文件读写操作，只是这段代码写在 MATLAB 里。把下面这个函数放到你的设计脚本后面调用，每次改参数重跑，头文件和源文件一起刷新，不会出现"代码里的系数和最新曲线对不上"这种事。
function export_sos_c(sos, g, outdir) % 把 SOS 矩阵和总增益写成 C 可编译的 .h / .c % sos: N x 6，每行 [b0 b1 b2 a0 a1 a2] nsec = size(sos, 1); sos(:, 4) = 1; % a0 统一归一化为 1 hf = fopen(fullfile(outdir, 'iir_sos.h'), 'w'); fprintf(hf, '#ifndef IIR_SOS_H\n#define IIR_SOS_H\n\n'); fprintf(hf, '#define IIR_SOS_SECTIONS %d\n\n', nsec); fprintf(hf, 'extern const float iir_gain;\n'); fprintf(hf, 'extern const float iir_sos[IIR_SOS_SECTIONS][6];\n\n#endif\n'); fclose(hf); cf = fopen(fullfile(outdir, 'iir_sos.c'), 'w'); fprintf(cf, '#include "iir_sos.h"\n\n'); fprintf(cf, 'const float iir_gain = %.9ef;\n\n', g); fprintf(cf, 'const float iir_sos[IIR_SOS_SECTIONS][6] = {\n'); for k = 1:nsec if k < nsec, sep = ','; else, sep = ''; end fprintf(cf, ' { %.9ef, %.9ef, %.9ef, 1.0f, %.9ef, %.9ef }%s\n', ... sos(k,1), sos(k,2), sos(k,3), sos(k,5), sos(k,6), sep); end fprintf(cf, '};\n'); fclose(cf); end逻辑上分两步：先写头文件，把节数和外部声明固定下来；再写源文件，把总增益和每个二阶节的六个数按行展开。%.9e保证单精度浮点能完整还原，位数少了会导致 C 里的系数和 MATLAB 不一致，位数多了又白占空间。sos(:,4)=1这一步别省，tf2sos返回的 a0 本来就该是 1，但如果你从别处拿到系数，不归一化就会出现整体增益偏差。调用时直接export_sos_c(sos, g, pwd)即可。
3.3 系数排列、a0 归一化与定点化约定
生成文件之前要把格式约定死，否则 C 侧一不小心就把 b1 当成 a1 用了。我习惯的排列是每行[b0 b1 b2 a0 a1 a2]，a0 恒为 1，处理时直接跳过第 4 列。如果目标平台是定点 DSP，还要决定 Q 格式：Q15 表示范围是 [-1,1)，精度 2^-15，适合 int16；Q31 精度 2^-31，适合 int32；带 FPU 的 MCU 直接用 float32 最省事。
| 格式 | 表示范围 | 量化精度 | 适用位宽 |
|---|---|---|---|
| Q15 | [-1, 1) | 2^-15 | int16 |
| Q31 | [-1, 1) | 2^-31 | int32 |
| float32 | 约 ±3.4e38 | 24 bit 有效 | 带 FPU 的 MCU |
定点化之前先看系数绝对值，如果某个 b0 已经接近或超过 1，直接按 Q15 放会溢出，这时要么把这一节的增益挪到前面的iir_gain里，要么整条链用 Q31。量化后一定要把系数读回 MATLAB 重画频响，量化误差对高 Q 值节的极点影响最大。
4. C 侧级联二阶节的实现与验证
4.1 转置直接型 II 的 C 实现
二阶节有几种实现形式，转置直接型 II（DF2T）状态量少、数值特性好，是定点实现里最常见的选择。它只需要两个状态量w[0]、w[1]，而且状态量的动态范围比直接型小。
#include "iir_sos.h" typedef struct { float w[IIR_SOS_SECTIONS][2]; /* 每个二阶节两个状态量 */ } iir_state_t; /* 转置直接型 II 单节处理，c = [b0 b1 b2 a0 a1 a2] */ static inline float biquad_df2t(const float c[6], float *w, float x) { float y = c[0] * x + w[0]; w[0] = c[1] * x - c[4] * y + w[1]; w[1] = c[2] * x - c[5] * y; return y; } /* 整条级联链，增益在入口统一缩放一次 */ float iir_process(iir_state_t *s, float x) { int k; x *= iir_gain; for (k = 0; k < IIR_SOS_SECTIONS; ++k) { x = biquad_df2t(iir_sos[k], s->w[k], x); } return x; }biquad_df2t里w[0]承担的是上一节输出的历史分量，w[1]是更早一个样本的延迟项。c[4]、c[5]对应 a1、a2，因为 a0 恒为 1，不用再做除法。iir_gain放在循环外做一次乘法，比在每个节里各自乘一遍少很多乘加，也避免中途某一节增益过冲。状态量w必须每个通道各有一份，左右声道共用一份会串音。
4.2 用阶跃与扫频对齐 MATLAB 和 C 的输出
代码写完不代表参数接对了。我一般用两组信号交叉验证：一段阶跃，看瞬态响应和超调；一段线性扫频，看幅频有没有出现不该有的凹陷。阶跃信号最能暴露极点量化误差，扫频则能暴露系数排列错位。MATLAB 侧用filter(sos, ...)得到参考输出，C 侧把同样的输入喂进iir_process，两条曲线叠在一起看。
% 阶跃输入，对比 MATLAB 参考输出 n = 4000; x = [ones(200,1); zeros(n-200,1)]; y_ref = filter(sos, 1, x); % sos 已在 2.1 求出 % C 侧把 iir_process 的输出存成 out.csv 后读回 y_c = readmatrix('out.csv'); plot(1:n, y_ref, 'b', 1:n, y_c, 'r--'); legend('MATLAB', 'C'); xlabel('sample'); ylabel('amplitude'); grid on;误差允许一个很小的量化底噪，但不能有系统性偏移。如果 C 侧整体放大或缩小了一个常数，多半是iir_gain丢了；如果只在某些频点偏差大，往高 Q 值那一节去查量化精度。阶跃响应检查完，再用扫频看通带边缘有没有提前滚降，那通常意味着过渡带参数填窄了。
4.3 Q 格式、溢出与二阶节顺序
定点实现里最常见的故障是中间结果溢出，而它往往不会立刻让程序崩溃，只是让输出偶尔冒出刺耳的爆音。DF2T 的状态量w[0]动态范围比输出大，用 Q15 时最容易被忽略。保险做法是给每节加饱和判断，或者干脆把状态量用 Q31 存、系数用 Q15 存，乘加后再移位。
二阶节的顺序也影响噪声。一般把极点离单位圆最远的那一节放最前面，让信号先衰减，再经过高 Q 值节，能减少中间溢出的概率。每节的增益分配也有讲究：如果某一节 b0 特别大，把它挪到整链的入口增益里，避免状态量被撑满。改完顺序后记得重跑一次阶跃对比，顺序变了输出在数值上应该几乎一致，只有瞬态底噪略有不同。
5. 量化后频响复核与系数重排技巧
定点化或者降到单精度之后，最值得做的事不是继续调参，而是把 C 里实际使用的系数读回 MATLAB，重画一条频响，和理想曲线叠在一起。这一步用freqz做，几行代码就能看出量化到底伤了多少。
% 把 C 里实际使用的系数按 Q15 量化后回灌验证 sos_q = sos; sos_q(:,1:3) = round(sos(:,1:3) * 2^15) / 2^15; sos_q(:,5:6) = round(sos(:,5:6) * 2^15) / 2^15; [h1, f] = freqz(sos, 4096, fs); [h2, ~] = freqz(sos_q, 4096, fs); plot(f, 20*log10(abs(h1)), 'b', f, 20*log10(abs(h2)), 'r--'); legend('理想系数', 'Q15 量化后'); xlabel('frequency (Hz)'); ylabel('magnitude (dB)'); grid on;判断标准很直接：通带内的两条曲线偏差应该在 0.1 dB 以内，阻带内的底噪抬高不超过几个 dB 就算合格。如果量化后通带边缘明显塌陷，说明某一节的 Q 值太高，Q15 撑不住，这时候要么把该节换成 Q31，要么降低rp放宽通带指标。反过来，如果阻带底噪抬高很多但通带几乎没变，多半是零点位置的量化误差，可以考虑把相邻两节的零点重新配对。
另一个实用技巧是二阶节的系数重排。tf2sos默认按数值稳定性配对，但定点实现里还有一层"哪一节的增益该被抽出去"的问题。把每节的 b0 单独看一遍，找一个最接近 1 的作为增益基准，其余的按比例缩放，能明显降低一节的动态范围压力。重排之后必须重新跑 4.2 的阶跃对比，确认输出和重排前在数值上一致，只有底噪级别的差别才算对。频响复核这件事不要只做一次，每次改量化位宽、改节顺序、改编译器浮点选项，都值得再画一次那条红色虚线。
本文还有配套的精品资源，点击获取