简介这份文档资料面向数字信号处理课程的学习者与实验者聚焦用双线性变换法设计IIR数字滤波器这一核心实验。内容围绕MATLAB实现展开涵盖双线性变换的基本原理、s域到z域的映射关系及预修正方法并给出低通、高通、带通、带阻四类滤波器的完整设计流程。资源包为单个doc文档约172KB内含实验目的、实验内容、MATLAB子函数bilinear的多种调用格式说明以及巴特沃斯数字低通滤波器的完整设计实例包括参数计算、程序代码、运行结果与频率响应分析。读者可据此掌握从指标预修正、阶数确定到bilinear变换与freqz验证的全过程理解频谱混叠消除与非线性畸变校正的要点。目前已有350人学习适合需要完成实验报告或巩固IIR滤波器设计方法的本科生与自学者参考。1. 双线性变换法设计IIR滤波器为什么你算出的系数总是不对很多人在数字信号处理实验里第一次用 MATLAB 设计 IIR 滤波器时都会遇到一个很玄学的现象按照教材公式一步步算filter函数跑出来的频响却和预期对不上要么截止频率偏了要么增益莫名其妙多了个零点几的系数。问题往往不在公式本身而在双线性变换法里那个容易被忽略的预畸变步骤以及模拟原型到数字域映射时分子分母系数的排列顺序。这个实验的核心链路其实很清晰给定数字滤波器的性能指标通带截止频率、阻带截止频率、通带最大衰减、阻带最小衰减先通过预畸变公式把数字频率映射回模拟频率再设计一个满足指标的模拟低通原型巴特沃斯、切比雪夫I型或II型最后用双线性变换把模拟传递函数映射成数字传递函数。MATLAB 把这个链路封装成了butter、cheby1、cheby2、ellip等函数但实验要求你理解每一步在做什么而不是只会调函数。适合谁看正在做数字信号处理课程实验的本科生、需要快速验证滤波器设计的研究生以及工作中偶尔要写滤波器系数但不想每次都翻教材的工程师。下面从原理到代码把这条链路拆开讲清楚。2. 预畸变与模拟原型双线性变换法绕不开的两步2.1 为什么双线性变换必须做预畸变双线性变换的本质是用 $s \frac{2}{T} \cdot \frac{1 - z^{-1}}{1 z^{-1}}$ 把 s 平面的虚轴映射到 z 平面的单位圆上。这个映射是非线性的具体表现为模拟频率 $\Omega$ 和数字频率 $\omega$ 之间满足$$\Omega \frac{2}{T} \tan\left(\frac{\omega}{2}\right)$$也就是说数字域上均匀分布的频率映射到模拟域后会被“压缩”到无穷远。如果你直接拿数字截止频率 $\omega_p$ 当作模拟截止频率 $\Omega_p$ 去设计原型映射回来之后截止频率一定会偏低。预畸变就是提前把这个非线性补偿掉$$\Omega_p \frac{2}{T} \tan\left(\frac{\omega_p}{2}\right)$$实际写代码时如果采样周期 $T$ 取 1MATLAB 默认做法预畸变公式简化为 $\Omega_p 2 \tan(\omega_p / 2)$。注意这里的 $\omega_p$ 是弧度制不是 Hz。很多人翻车就翻在把 Hz 直接代进去或者忘了乘 $\pi/180$。提示MATLAB 的butter等函数内部已经自动做了预畸变你传入归一化频率Nyquist 频率为 1即可。但实验如果要求手写双线性变换就必须自己算这一步。2.2 巴特沃斯、切比雪夫I型、切比雪夫II型怎么选模拟原型的选择直接决定滤波器的阶数和通带/阻带特性。下面这张表是我在实验里反复对比后整理的选型依据原型类型通带特性阻带特性阶数趋势相位线性度适用场景巴特沃斯最平坦无纹波单调衰减最高较好对通带平坦度要求高阶数不敏感切比雪夫I型等纹波单调衰减中等一般需要陡峭过渡带允许通带纹波切比雪夫II型最平坦等纹波中等一般需要陡峭过渡带允许阻带纹波椭圆等纹波等纹波最低较差阶数极度受限对相位无要求选型逻辑很简单先看阶数能不能接受。如果巴特沃斯算出来 12 阶切比雪夫I型只要 6 阶而你的 DSP 或 FPGA 资源紧张那就果断换切比雪夫。如果通带里不允许有任何纹波比如音频处理那就只能巴特沃斯。椭圆滤波器阶数最低但通带和阻带都有纹波相位非线性也最严重一般只在阶数卡得很死的时候用。2.3 手写双线性变换的完整步骤假设要设计一个数字低通滤波器采样率 $f_s 1000$ Hz通带截止频率 $f_p 100$ Hz阻带截止频率 $f_{st} 200$ Hz通带最大衰减 $R_p 1$ dB阻带最小衰减 $R_s 40$ dB。用巴特沃斯原型手写步骤如下第一步把 Hz 转成弧度制归一化频率。$\omega_p 2\pi f_p / f_s 0.2\pi$$\omega_{st} 2\pi f_{st} / f_s 0.4\pi$。第二步预畸变得到模拟角频率。取 $T 1$$\Omega_p 2\tan(0.1\pi) \approx 0.6498$$\Omega_{st} 2\tan(0.2\pi) \approx 1.4531$。第三步计算巴特沃斯原型阶数$$N \geq \frac{\log_{10}\left(\frac{10^{R_s/10} - 1}{10^{R_p/10} - 1}\right)}{2\log_{10}\left(\frac{\Omega_{st}}{\Omega_p}\right)}$$代入数值算出来 $N \geq 5.2$取 $N 6$。第四步求模拟低通原型传递函数 $H_a(s)$截止频率取 $\Omega_c \Omega_p / (10^{R_p/10} - 1)^{1/(2N)}$。第五步代入 $s 2\frac{1 - z^{-1}}{1 z^{-1}}$整理成 $H(z) \frac{\sum b_k z^{-k}}{\sum a_k z^{-k}}$ 的形式。这五步里第三步和第四步最容易算错。阶数公式里的 $\log_{10}$ 别写成自然对数$\Omega_c$ 的指数别漏掉分母的 2。我一般会在 MATLAB 里单独写一个小脚本验证每一步的中间结果而不是一口气算完再看最后频响。3. 用 MATLAB 把设计链路跑通从系数到频响3.1 直接调用 butter 函数的最小可用代码先给一个能直接跑的最小示例采样率 1000 Hz通带 100 Hz阻带 200 Hz通带衰减 1 dB阻带衰减 40 dB% 参数定义 fs 1000; % 采样率 Hz fp 100; % 通带截止频率 Hz fst 200; % 阻带截止频率 Hz Rp 1; % 通带最大衰减 dB Rs 40; % 阻带最小衰减 dB % 归一化频率Nyquist 频率为 1 Wp fp / (fs/2); Wst fst / (fs/2); % 计算巴特沃斯原型的最小阶数和截止频率 [N, Wn] buttord(Wp, Wst, Rp, Rs); % 设计数字滤波器返回分子分母系数 [b, a] butter(N, Wn); % 查看系数 disp(分子系数 b:); disp(b); disp(分母系数 a:); disp(a); % 画频响 figure; freqz(b, a, 1024, fs); title(IIR 低通滤波器频响);这段代码里buttord负责算阶数和截止频率butter负责生成系数。注意Wn是归一化频率不是 Hz。freqz的第三个参数是频率点数第四个参数是采样率传了采样率之后横轴直接显示 Hz不传就是归一化频率。3.2 系数排列顺序与 filter 函数的对应关系MATLAB 的filter(b, a, x)实现的差分方程是$$y[n] \frac{1}{a_1}\left(\sum_{k1}^{N_b} b_k x[n-k1] - \sum_{k2}^{N_a} a_k y[n-k1]\right)$$也就是说b和a都是按 $z^0, z^{-1}, z^{-2}, \ldots$ 的顺序排列的。a(1)通常是 1如果不是 1filter会自动归一化。很多人从模拟传递函数手推数字系数时容易把 $z^{-1}$ 的系数顺序写反导致频响完全不对。一个简单的验证方法用tf(b, a, 1/fs)建立传递函数模型和freqz的结果对比如果一致说明系数顺序没问题。注意如果a(1)不等于 1先手动除以a(1)再做后续分析避免数值误差累积。3.3 用 freqz 验证频响并检查关键指标设计完滤波器不能只看系数必须验证频响是否满足指标。下面这段代码在频响图上标注通带和阻带边界并计算实际衰减% 计算频响 [H, f] freqz(b, a, 2048, fs); % 计算幅度 dB H_dB 20 * log10(abs(H)); % 找通带边界和阻带边界对应的索引 [~, idx_p] min(abs(f - fp)); [~, idx_st] min(abs(f - fst)); % 输出实际指标 fprintf(通带边界 %.1f Hz 处衰减: %.2f dB\n, f(idx_p), H_dB(idx_p)); fprintf(阻带边界 %.1f Hz 处衰减: %.2f dB\n, f(idx_st), H_dB(idx_st)); % 画图并标注 figure; plot(f, H_dB); hold on; yline(-Rp, r--, 通带最大衰减); yline(-Rs, r--, 阻带最小衰减); xline(fp, g--, 通带边界); xline(fst, g--, 阻带边界); xlabel(频率 (Hz)); ylabel(幅度 (dB)); title(IIR 滤波器频响验证); grid on;跑完之后看两个数通带边界处的衰减应该大于 $-R_p$阻带边界处的衰减应该小于 $-R_s$。如果阻带衰减不够说明阶数算少了或者Wn选偏了。我一般会把buttord算出来的N和Wn打印出来再手动微调Wn看频响怎么变这样能直观理解阶数和截止频率对过渡带的影响。3.4 不同原型函数的参数差异butter、cheby1、cheby2、ellip的调用方式略有不同% 巴特沃斯butter(N, Wn) [b1, a1] butter(N, Wn); % 切比雪夫I型cheby1(N, Rp, Wn)Rp 是通带纹波 dB [b2, a2] cheby1(N, Rp, Wn); % 切比雪夫II型cheby2(N, Rs, Wn)Rs 是阻带衰减 dB [b3, a3] cheby2(N, Rs, Wn); % 椭圆ellip(N, Rp, Rs, Wn)同时指定通带纹波和阻带衰减 [b4, a4] ellip(N, Rp, Rs, Wn);切比雪夫I型的Rp直接决定通带纹波大小Rp越小纹波越小但阶数越高。切比雪夫II型的Rs决定阻带纹波注意它指定的是阻带最小衰减不是纹波峰峰值。椭圆滤波器的Rp和Rs同时约束通带和阻带阶数最低但相位最差。选哪个取决于你的指标里对纹波和阶数的容忍度。4. 避坑双线性变换法设计 IIR 的五个血泪教训4.1 频率单位混用导致截止频率偏移现象设计一个截止频率 100 Hz 的滤波器频响图上截止点跑到了 160 Hz 附近。原因buttord和butter接受的归一化频率是以 Nyquist 频率为 1 的如果直接把 Hz 传进去相当于把截止频率放大了 $f_s/2$ 倍。另一种情况是手写预畸变时把 Hz 当成弧度制代入 $\tan$ 函数。解决统一用W f / (fs/2)做归一化手写公式时先把 Hz 转成弧度制 $\omega 2\pi f / f_s$再代入预畸变公式。在代码里加一行assert(Wp 1 Wst 1)做保护。4.2 阶数取整方向搞反现象buttord算出来 $N 5.2$手动取整成 5结果阻带衰减差 3 dB 不达标。原因阶数公式算出来的是下限必须向上取整。取 5 意味着阶数不够过渡带变缓阻带衰减自然不够。解决用ceil而不是round或fix。MATLAB 的buttord内部已经做了向上取整手写公式时务必注意。如果向上取整后阶数太高换切比雪夫或椭圆原型。4.3 系数顺序与差分方程不对应现象filter跑出来的信号完全失真频响和freqz画出来的对不上。原因手推传递函数时把分子分母的系数顺序写反了或者把 $z^{-1}$ 的系数放在了 $z^0$ 的位置。MATLAB 的filter要求b和a都按 $z^0, z^{-1}, z^{-2}$ 排列。解决用tf(b, a, 1/fs)建立模型和freqz对比。如果两者不一致说明系数顺序有问题。另一个办法是用filter跑一个单位脉冲看冲激响应和impz是否一致。4.4 预畸变时采样周期 T 取值不一致现象手写双线性变换得到的系数和butter生成的系数不一样频响也有偏差。原因双线性变换公式里 $s \frac{2}{T} \cdot \frac{1 - z^{-1}}{1 z^{-1}}$如果预畸变时用了 $T 1$代入变换时又用了 $T 1/f_s$两次的 $T$ 不一致会导致频率映射错位。解决全程统一 $T$ 的取值。MATLAB 默认 $T 1$预畸变和变换都用这个值。如果非要显式写 $T$就在代码开头定义T 1后面全部引用这个变量。4.5 忽略数值精度导致高阶滤波器不稳定现象设计一个 12 阶以上的巴特沃斯滤波器filter跑出来的结果发散或者有持续振荡。原因高阶 IIR 滤波器的极点非常靠近单位圆直接用b和a系数实现时系数的微小量化误差会导致极点跑到单位圆外系统不稳定。解决高阶滤波器改用二阶节级联形式MATLAB 里用zp2sos把零极点转成二阶节再用sosfilt滤波。或者降低阶数换切比雪夫/椭圆原型。如果必须用高阶考虑用tf2sos转换后再分析。5. 从系数到落地二阶节级联与定点化验证5.1 为什么高阶滤波器要用二阶节直接形式Direct Form的 IIR 滤波器在阶数超过 6 之后系数量化误差会显著影响极点位置。一个 12 阶滤波器的分母系数可能有 6 个以上不同的数量级用双精度浮点还能撑住一旦转到定点 DSP 或 FPGA 上系数截断误差足以让滤波器崩溃。二阶节级联SOS把高阶传递函数分解成多个二阶节的乘积每个二阶节独立处理量化误差不会跨节累积。MATLAB 里转换和滤波的代码如下% 设计高阶滤波器 N 12; Wn 0.3; [b, a] butter(N, Wn); % 转换为二阶节形式 [sos, g] tf2sos(b, a); % 查看二阶节系数 disp(SOS 矩阵:); disp(sos); disp(整体增益 g:); disp(g); % 用 sosfilt 滤波 x randn(1000, 1); y sosfilt(sos, x); % 对比直接形式和 SOS 形式的频响 [H1, f] freqz(b, a, 1024); [H2, f] freqz(sos, 1024); figure; plot(f, 20*log10(abs(H1)), b, f, 20*log10(abs(H2)), r--); legend(直接形式, SOS 形式); xlabel(归一化频率); ylabel(幅度 (dB)); title(直接形式与 SOS 形式频响对比);sos矩阵每行是一个二阶节的系数格式为 $[b_0, b_1, b_2, a_0, a_1, a_2]$。g是所有二阶节的总增益。sosfilt按节依次滤波数值稳定性比filter好很多。如果要在定点上实现每个二阶节的系数单独量化误差可控。5.2 定点化验证用 fi 对象模拟量化误差如果目标平台是定点 DSP 或 FPGA需要在 MATLAB 里先模拟系数量化后的频响。用 Fixed-Point Designer 的fi对象可以快速验证% 假设系数用 16 位定点1 位符号位15 位小数位 word_length 16; frac_length 15; % 量化 SOS 系数 sos_fi fi(sos, 1, word_length, frac_length); g_fi fi(g, 1, word_length, frac_length); % 转回 double 看量化后的频响 sos_q double(sos_fi); g_q double(g_fi); % 对比量化前后频响 [H_ideal, f] freqz(sos, 1024); [H_quant, f] freqz(sos_q, 1024); figure; plot(f, 20*log10(abs(H_ideal)), b, f, 20*log10(abs(H_quant)), r--); legend(量化前, 量化后); xlabel(归一化频率); ylabel(幅度 (dB)); title(系数量化对频响的影响);跑完之后重点看两个地方通带纹波有没有变大阻带衰减有没有变差。如果量化后阻带衰减掉了 10 dB 以上说明字长不够需要增加小数位或者换用二阶节缩放scaling来优化动态范围。5.3 一个我常用的验证习惯每次设计完滤波器我不会只看freqz的图而是会跑三个额外检查第一用impz看冲激响应有没有发散第二用stepz看阶跃响应有没有持续振荡第三用白噪声跑一遍filter或sosfilt看输出有没有异常尖峰。这三个检查加起来不到十行代码但能提前发现 90% 的稳定性问题。还有一个习惯把设计参数和生成的系数存成.mat文件文件名带上日期和关键参数比如iir_lp_fp100_fst200_Rp1_Rs40_20250101.mat。下次要改指标时直接改参数重跑不用从头推公式。这个习惯帮我省了很多重复劳动也避免了“上次那个能跑的系数到底是多少”的后悔药问题。希望帮到你。本文还有配套的精品资源点击获取
