FDAtool设计IIR滤波器：二阶节系数导出与C语言定点实现

简介：这份PDF文档围绕MATLAB FDAtool在IIR数字滤波器设计中的参数生成与C语言代码导出展开，面向数字信号处理学习者、嵌入式开发人员及需要将滤波器移植到C环境的工程师。内容从角频率与采样频率的关系、通带与阻带截止频率等基本概念切入，对比FIR与IIR在阶数、线性相位和计算复杂度上的差异；并以8kHz采样率下去除50Hz电频干扰、保留80-3200Hz语音的带通滤波器为例，说明36阶、18个二阶节的实现结构。文档重点展示如何将FDAtool设计结果导出为C头文件，解释系数按二阶节排列、增益项设置以及转换为单段形式后的可读性，便于在嵌入式程序或其他软件环境中复用。包内共1个PDF文件，约1.08MB，目前已有1423人学习。适合希望理解IIR滤波器参数到C代码落地过程、查漏补缺的读者参考。

1. 手写 IIR 系数为什么容易翻车，FDAtool 能替你做什么

见过太多这类现场：C 代码本身没写错，二阶节的框图也是照着教科书抄的，跑起来截止频率偏了两三成，或者通带边缘鼓出一个包。根因通常不在实现层，而在设计层——做双线性变换时漏了频率预畸变，或者顺手把模拟原型的系数直接填进了数字滤波器。IIR 是有反馈的，极点离单位圆越近，系数上第四位小数的误差就能把频响推歪。

FDAtool 是 MATLAB 滤波器设计器的旧名字，新版入口叫 Filter Designer，命令行敲fdatool或filterDesigner都能起来（跟 MATLAB 下载安装教程里那些版本差异没关系，近十几年的版本都在）。它把通带阻带指标、结构选择、系数量化和 C 头文件导出串成一条链，最后交到你手上的是一张二阶节系数表加一个增益因子。适合两类人：要在 MCU、DSP 上跑实时滤波的固件工程师，和要做离线数据预处理、需要拿 MATLAB 画图核对频响的人。下面按设计、导出、C 实现、验证四步拆开，每步都给能直接跑的代码和参数。

2. 在 FDAtool 里设计 IIR：指标怎么填，导出为什么选二阶节

2.1 IIR 与 FIR 的取舍，先把阶数账算清楚

同样的过渡带宽度，IIR 的阶数通常只有 FIR 的十分之一到二十分之一。对 48 kHz 采样、要滤 50 Hz 工频那类场景，FIR 可能要几百阶，MCU 上每样本几百次乘加直接吃掉一半主频；换成椭圆 IIR 五六阶就够了。代价是相位非线性，如果下游要做波形相关或相位差测量，得再加全通均衡器。

2.2 FDAtool 主界面里真正要动的那几组参数

打开后先别急着点 Design Filter，从上到下把这几栏填对：

• Response Type：Lowpass / Highpass / Bandpass / Bandstop，先定这个。
• Design Method：选 IIR，再选 Butterworth（最平）、Chebyshev Type I（等波纹通带）、Elliptic（等波纹通带加阻带，阶数最低）。
• Filter Order：可以先 Specify order 手填，也可以选 Minimum order 让工具按指标算最小阶。
• Frequency Specifications：Units 一定要设成 Hz，然后填 Fs、Fpass、Fstop。踩坑最多的地方就是 Units 还是归一化，人却按 Hz 填了 100。
• Magnitude Specifications：Apass 单位 dB，Astop 单位 dB。椭圆滤波器两个都要填。

拿一组具体数字：Fs = 1000 Hz，Fpass = 100 Hz，Fstop = 150 Hz，Apass = 0.5 dB，Astop = 40 dB。这组指标下椭圆滤波器大概 5 阶，巴特沃斯要 8 阶以上。填完点 Design Filter，上方会画出幅频响应，看过渡带是不是从 100 到 150 掉到 -40 dB 以下。

提示：改完任何一栏都要重新点一次 Design Filter，界面上的曲线不会自动刷新，很多人以为改参数没生效。

2.3 用脚本复现同一组系数，才进得了版本库

GUI 点出来的东西没法 diff，也没法让同事复现。实际做法是用几行脚本把同一套指标算出来，把系数矩阵存成文件跟着代码一起提交。

% 设计指标 Fs = 1000; Fpass = 100; Fstop = 150; Ap = 0.5; Ast = 40; % 归一化到 Nyquist 频率 Wp = Fpass / (Fs/2); Ws = Fstop / (Fs/2); % 最小阶数与截止频率 [N, Wn] = ellipord(Wp, Ws, Ap, Ast); % N 为最小阶数，Wn 为归一化截止频率 [z, p, k] = ellip(N, Ap, Ast, Wn); % 零极点增益形式 % 转成二阶节，'up' 表示把增益平均分摊到各节，'inf' 表示用无穷范数缩放 [sos, g] = zp2sos(z, p, k, 'up', 'inf'); disp(sos); % 每行一个二阶节：b0 b1 b2 1 a1 a2 disp(g); % 分摊后通常接近 1 % 直接画频响核对 fvtool(sos, 'Fs', Fs);

ellipord返回的最小阶数是保证同时满足 Apass 和 Astop 的最低值，实际可以往上取一阶留余量。zp2sos的'up'参数是关键：不写它的话，全部增益都压在第一个二阶节上，那节的中间状态值会大出好几个数量级，定点实现时极易溢出。fvtool用的是同一个sos矩阵，画的曲线和 FDAtool 界面里那根应该完全重合，不重合说明指标填错了。

2.4 为什么导出环节要选 SOS 而不是 b/a

butter、ellip也能直接给你b、a两个多项式系数，看着更简单：y = filter(b, a, x)。问题在于高阶 IIR 的极点对系数扰动极其敏感。一个 8 阶滤波器的a系数里，某一位的变化可能让某一对极点跨越单位圆，滤波器直接发散。拆成四个二阶节之后，每节只有两个极点，灵敏度下降一到两个数量级。

另一个现实原因是定点。二阶节的a1理论范围是 (-2, 2)，a2在 (-1, 1)，量程固定，容易选 Q 格式；合成后的高阶a系数动态范围能到几十甚至上百，统一 Q 格式根本放不下。所以 FDAtool 的导出目标、zp2sos、以及 C 侧的滤波器结构，全都围绕二阶节来组织。

3. 从 FDAtool 到 C 语言文件：系数导出、Q 格式与头文件组织

3.1 FDAtool 自带的 Targets → Generate C header

设计完成后，菜单 Targets → Generate C header，会弹窗让你确认系数类型（double / single / fixed-point）。生成的文件一般叫fdacoefs.h，里面是几个二维数组B[NSEC][3]、A[NSEC][3]，再加一个各节缩放因子数组。它的好处是零手工劳动，坏处是格式固定、命名不便、数组长度写死在宏里。字段名和数组组织方式在不同 MATLAB 版本间有调整，以你手上那版生成的为准，别照抄别人的旧文件。

如果只是自己临时验证，直接用它生成的头文件完全够用。要进产品代码，我更倾向自己写生成脚本，控制命名、加上静态断言和注释。

3.2 用 fprintf 自己生成 .h 和 .c

下面这个函数把sos矩阵直接写成 C 头文件，浮点版和定点版都能出。核心是fprintf格式化输出，MATLAB 的文件读写在这里当模板引擎用。

function write_iir_header(fname, sos, g, q) % fname: 输出文件名；sos: N×6 二阶节矩阵；g: 总增益；q: 小数位数，0 表示浮点 fid = fopen(fname, 'w'); assert(fid > 0, '无法打开文件 %s', fname); nsec = size(sos, 1); fprintf(fid, '/* 由 MATLAB 生成，请勿手工修改 */\n'); fprintf(fid, '#ifndef IIR_COEF_H\n#define IIR_COEF_H\n\n'); fprintf(fid, '#define IIR_NSEC %d\n', nsec); fprintf(fid, '#define IIR_Q %d\n\n', q); if q == 0 % 浮点版本 fprintf(fid, 'static const float iir_sos[%d][5] = {\n', nsec); for k = 1:nsec % 注意：a0 恒为 1，不存储，只存 a1 a2 fprintf(fid, ' { %+.9ef, %+.9ef, %+.9ef, %+.9ef, %+.9ef },\n', ... sos(k,1), sos(k,2), sos(k,3), sos(k,5), sos(k,6)); end fprintf(fid, '};\n\n'); fprintf(fid, 'static const float iir_gain = %.9ef;\n', g); else % 定点版本，系数按 Qq 定标并四舍五入 scale = 2^q; fprintf(fid, 'static const int32_t iir_sos[%d][5] = {\n', nsec); for k = 1:nsec fprintf(fid, ' { %6d, %6d, %6d, %6d, %6d },\n', ... round(sos(k,1)*scale), round(sos(k,2)*scale), round(sos(k,3)*scale), ... round(sos(k,5)*scale), round(sos(k,6)*scale)); end fprintf(fid, '};\n\n'); fprintf(fid, 'static const int32_t iir_gain = %d;\n', round(g*scale)); end fprintf(fid, '\n#endif\n'); fclose(fid); end

调用方式和生成结果的形状：

write_iir_header('iir_coef.h', sos, g, 0); % 浮点版 write_iir_header('iir_coef_q13.h', sos, g, 13); % Q13 定点版

生成出来的数组每行 5 个数，顺序是b0 b1 b2 a1 a2。这里故意丢掉每行的a0——它恒等于 1，写进代码只是浪费一次判读。round用四舍五入而不是截断，截断会让所有系数系统性偏小，低频段增益会整体下移。

3.3 Q 格式怎么选：a1 决定了下限

定点实现里最容易错的就是给a1选了不够的整数位。二阶节稳定时a1 ∈ (-2, 2)，a2 ∈ (-1, 1)，b系数经过缩放通常也落在 (-2, 2)。Q15 只能表示 (-1, 1)，a1一旦超过 1 就直接溢出回绕。所以定点方案里通行的选择是 Q13 或 Q14，留 2 位符号加整数位。

验证选型是否够用的方法很直接，在 MATLAB 里跑一行：

q = 13; assert(max(abs(sos(:,5))) < 2^(3-1), 'a1 超出 Q13 范围，请降低 Q 值'); assert(max(abs(sos(:,1:3)), [], 'all') < 2^(3-1), 'b 系数超出 Q13 范围');

2^(3-1)里的 3 是 Q13 的整数位数（含符号位），写成这样是为了让你改 Q 值时只改一处。系数一旦超范围，定点版本的频响会完全不对，但代码不会报错，只会输出一堆噪声——这是最难查的一类问题。

4. C 侧实现：级联二阶节的结构与文件读写对拍

4.1 直接 II 型转置：每节只要两个状态变量

级联二阶节推荐用直接 II 型转置（TDF-II），每个二阶节只需两个状态变量，中间不需要存x[n-1]、x[n-2]、y[n-1]、y[n-2]四个延迟单元。浮点版本实现如下，可以直接编译。

#include <stdint.h> #include "iir_coef.h" typedef struct { float s1; float s2; } biquad_state_t; /* 处理一个样本，按顺序穿过所有二阶节 */ float iir_process(biquad_state_t *st, float x) { for (int k = 0; k < IIR_NSEC; ++k) { const float b0 = iir_sos[k][0]; const float b1 = iir_sos[k][1]; const float b2 = iir_sos[k][2]; const float a1 = iir_sos[k][3]; const float a2 = iir_sos[k][4]; /* TDF-II：先算输出，再更新两个状态 */ float y = b0 * x + st[k].s1; st[k].s1 = b1 * x - a1 * y + st[k].s2; st[k].s2 = b2 * x - a2 * y; x = y; /* 本节输出作为下一节输入 */ } return x * iir_gain; }

三个要点。第一，iir_gain放在循环外面只乘一次，而不是摊进每节的b0，这样b0的量级可控，也少三次乘法。第二，状态更新必须在算出y之后，顺序颠倒的话滤波器结构就变成了直接 I 型，效果不一样。第三，iir_sos[k][3]、[4]里存的是a1、a2本身，不是它们的相反数——MATLAB 的sos格式列的是1 a1 a2，推导差分方程时y的系数要移项，符号别抄错。

定点版本只在数据类型和移位位置上不同：把float换成int32_t，每次累加后右移IIR_Q位，乘法用 64 位中间变量防溢出。

/* 定点版单个二阶节，acc 用 int64_t 收中间结果 */ int32_t biquad_fixed(biquad_state_t *s, int32_t x) { int64_t acc = (int64_t)iir_sos[k][0] * x + ((int64_t)s->s1 << IIR_Q); /* 中间省略 a1/a2 的移位乘加，实际代码需要完整展开 */ return (int32_t)(acc >> IIR_Q); }

(int64_t)s->s1 << IIR_Q是把状态变量左移到与乘积相同的定标上再做加法，少了这一步结果会差 2 的 Q 次方倍。中间结果必须用int64_t，两个 Q13 的 16 位系数相乘就已经是 26 位，累加几次就顶到 32 位边界。

4.2 用 C 读写文件跑批量验证数据

算法写完不能只喂几个正弦。常规做法是在 MATLAB 里生成一段测试信号写进二进制文件，C 读进来跑完整流程，结果再写文件，最后回 MATLAB 比对。这套流程涉及的 C 语言文件读写操作代码其实就三组函数。

#include <stdio.h> #include <stdlib.h> /* 读入 float 二进制文件，长度写回 n */ float *read_floats(const char *path, size_t *n) { FILE *fp = fopen(path, "rb"); if (!fp) { perror("fopen read"); return NULL; } fseek(fp, 0, SEEK_END); long bytes = ftell(fp); /* 文件总字节数 */ rewind(fp); *n = (size_t)bytes / sizeof(float); float *buf = (float *)malloc(bytes); if (fread(buf, sizeof(float), *n, fp) != *n) { fprintf(stderr, "short read on %s\n", path); free(buf); fclose(fp); return NULL; } fclose(fp); return buf; } /* 以文本形式写出，方便直接和 MATLAB 的 load 对接 */ int write_text(const char *path, const float *buf, size_t n) { FILE *fp = fopen(path, "w"); if (!fp) { perror("fopen write"); return -1; } for (size_t i = 0; i < n; ++i) fprintf(fp, "%.9g\n", buf[i]); fclose(fp); return 0; }

fseek+ftell+rewind这套是为了先拿到文件长度再一次性malloc，比循环fgetc快得多，也不会因为缓冲区大小猜错而截断数据。写文件用文本格式%.9g，float的有效十进制位大约是 9 位，%.6f会丢精度，导致后续比对时误差分不清是算法问题还是格式问题。

MATLAB 侧生成输入和比对的配套代码：

% 生成测试信号：白噪声 + 两个单音，覆盖通带和阻带 rng(0); x = 0.5*randn(20000,1) + sin(2*pi*80*(0:19999)'/Fs); fid = fopen('input.bin', 'wb'); fwrite(fid, x, 'float32'); % 头文件里 float 对应 float32 fclose(fid); % C 跑完后读回结果比对 y_ref = filter(sos, x); % 注意 filter 用 sos 时的调用形式见 doc yc = load('output.txt'); err = y_ref - yc; fprintf('RMS 误差 = %g，峰值误差 = %g\n', rms(err), max(abs(err)));

fwrite的'float32'必须和 C 侧的float对上。有些平台上double是 8 字节，写'double'而 C 用float读，会得到一串完全错位的数据，误差看起来像算法崩溃，实际只是类型不匹配。

4.3 对拍时误差多大才算正常

浮点实现下，单精度和双精度之间的 RMS 误差通常在 -100 dB 量级，用上面 0.5 幅度的信号，绝对误差应该在 1e-6 以下。如果误差到了 1e-3 这个量级，八成是系数抄错了或者状态更新顺序不对。定点实现另算：Q13 的系数量化会带来 -70 dB 左右的底噪，这是物理上限，不是 bug。

表格最后一行的意思是：如果你按 Q15 导出a1得到一堆错误结果，那不是代码问题，是格式选错了。

5. 定点溢出、极限环与实时性：几个能定位问题的技巧

定点 IIR 最典型的两种病：一是溢出，二是不输入信号了输出还在小幅振荡。后者叫极限环，根源是状态变量s1、s2被舍入到量化格点上，反馈回路靠这个残差自己维持住。定位方法很朴素：把输入置零，跑两千个样本，看输出尾部的峰峰值。如果它稳定在某个小常数上不衰减，就是极限环。

/* 极限环自检：零输入下的残余振荡 */ biquad_state_t st[IIR_NSEC] = {0}; float peak = 0.0f; for (int i = 0; i < 2000; ++i) { float y = i < 10 ? 1.0f : 0.0f; /* 10 个样本的脉冲，之后全零 */ y = iir_process(st, y); if (i > 500 && fabsf(y) > peak) peak = fabsf(y); } printf("极限环残余峰值 = %.3e\n", peak);

i > 500是为了跳过脉冲本身的响应尾巴，只看后期残余。浮点实现这个值应该在 1e-9 以下，定点 Q13 大概在 1e-4 到 1e-5，如果到了 1e-2 就说明 Q 值取得太低或者滤波器阶数被拆错了。

溢出的排查思路相反，要盯中间量而不是输出。把每个二阶节的y都打出来看峰值，超过 1.0 的节就是要处理的。手段有三种：把总增益分摊给每个二阶节（就是前面zp2sos的'up'参数干的事）、给累加器留保护位、以及改用饱和运算替代回绕。三者的效果差异很大：分摊增益不增加任何计算量，是最该先试的。

实时性上，每个二阶节每样本是 5 次乘加，6 阶滤波器 3 个节就是 15 次乘加加若干移位。Cortex-M4 上跑 48 kHz 单通道，这个量级基本不占什么主频。真正吃时间的是分支和函数调用——把iir_process里的for循环展开成宏、或者直接把每个节的代码写死，比开-O2更有效。用 ARM 平台的话，CMSIS-DSP 里的arm_biquad_cascade_df1_f32系列可以直接用，它传入的就是sos格式系数数组和节数，接口和你手上的数据结构几乎一一对应，切换成本很低。要用之前先确认a1、a2的符号约定，CMSIS 用的是标准差分方程形式，和前面 C 代码里y = b0*x + s1那套写法一致，系数不需要取反。

本文还有配套的精品资源，点击获取