简介这份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 五六阶就够了。代价是相位非线性如果下游要做波形相关或相位差测量得再加全通均衡器。维度IIR巴特沃斯/切比雪夫/椭圆FIR同等过渡带所需阶数低椭圆 58 阶常见高常 60 阶以上相位非线性需均衡可严格线性稳定性依赖极点位置量化敏感无反馈恒稳每样本乘加少多定点风险高需 Q 格式设计低2.2 FDAtool 主界面里真正要动的那几组参数打开后先别急着点 Design Filter从上到下把这几栏填对Response TypeLowpass / Highpass / Bandpass / Bandstop先定这个。Design Method选 IIR再选 Butterworth最平、Chebyshev Type I等波纹通带、Elliptic等波纹通带加阻带阶数最低。Filter Order可以先 Specify order 手填也可以选 Minimum order 让工具按指标算最小阶。Frequency SpecificationsUnits 一定要设成 Hz然后填 Fs、Fpass、Fstop。踩坑最多的地方就是 Units 还是归一化人却按 Hz 填了 100。Magnitude SpecificationsApass 单位 dBAstop 单位 dB。椭圆滤波器两个都要填。拿一组具体数字Fs 1000 HzFpass 100 HzFstop 150 HzApass 0.5 dBAstop 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/abutter、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 位符号加整数位。格式总位数整数位含符号可表示范围最小分辨率适合场景Q15161[-1, 1)3.05e-5只有 b 系数的 FIRQ14162[-2, 2)6.10e-5极点在原点附近的 IIRQ13163[-4, 4)1.22e-4一般 IIR 二阶节推荐Q12164[-8, 8)2.44e-4极点很靠近单位圆、a1 接近 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; }fseekftellrewind这套是为了先拿到文件长度再一次性malloc比循环fgetc快得多也不会因为缓冲区大小猜错而截断数据。写文件用文本格式%.9gfloat的有效十进制位大约是 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。实现方式系数存储预期 RMS 误差相对输入幅度 0.5常见异常原因双精度对单精度double / float 1e-7增益因子漏乘定点 Q13int16约 1e-4系数溢出回绕、移位方向错定点 Q13 且加饱和int16约 1e-4中间累加未用 64 位定点 Q15int16通常发散a1 超出 (-1,1)表格最后一行的意思是如果你按 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那套写法一致系数不需要取反。本文还有配套的精品资源点击获取
简介这份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 五六阶就够了。代价是相位非线性如果下游要做波形相关或相位差测量得再加全通均衡器。维度IIR巴特沃斯/切比雪夫/椭圆FIR同等过渡带所需阶数低椭圆 58 阶常见高常 60 阶以上相位非线性需均衡可严格线性稳定性依赖极点位置量化敏感无反馈恒稳每样本乘加少多定点风险高需 Q 格式设计低2.2 FDAtool 主界面里真正要动的那几组参数打开后先别急着点 Design Filter从上到下把这几栏填对Response TypeLowpass / Highpass / Bandpass / Bandstop先定这个。Design Method选 IIR再选 Butterworth最平、Chebyshev Type I等波纹通带、Elliptic等波纹通带加阻带阶数最低。Filter Order可以先 Specify order 手填也可以选 Minimum order 让工具按指标算最小阶。Frequency SpecificationsUnits 一定要设成 Hz然后填 Fs、Fpass、Fstop。踩坑最多的地方就是 Units 还是归一化人却按 Hz 填了 100。Magnitude SpecificationsApass 单位 dBAstop 单位 dB。椭圆滤波器两个都要填。拿一组具体数字Fs 1000 HzFpass 100 HzFstop 150 HzApass 0.5 dBAstop 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/abutter、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 位符号加整数位。格式总位数整数位含符号可表示范围最小分辨率适合场景Q15161[-1, 1)3.05e-5只有 b 系数的 FIRQ14162[-2, 2)6.10e-5极点在原点附近的 IIRQ13163[-4, 4)1.22e-4一般 IIR 二阶节推荐Q12164[-8, 8)2.44e-4极点很靠近单位圆、a1 接近 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; }fseekftellrewind这套是为了先拿到文件长度再一次性malloc比循环fgetc快得多也不会因为缓冲区大小猜错而截断数据。写文件用文本格式%.9gfloat的有效十进制位大约是 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。实现方式系数存储预期 RMS 误差相对输入幅度 0.5常见异常原因双精度对单精度double / float 1e-7增益因子漏乘定点 Q13int16约 1e-4系数溢出回绕、移位方向错定点 Q13 且加饱和int16约 1e-4中间累加未用 64 位定点 Q15int16通常发散a1 超出 (-1,1)表格最后一行的意思是如果你按 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那套写法一致系数不需要取反。本文还有配套的精品资源点击获取
