Lab 4 - IIR滤波器设计

数字信号处理实验
南京工业大学浦江学院 计算机与通信工程学院 | 通信工程专业

理论实验目的

🎵 生活场景：为什么需要学习模拟滤波器？

想象你正在开发一款音乐App的均衡器功能。用户滑动一个滑块，就能增强低音或高音。你在课堂上学到了脉冲响应不变法(IIM)和双线性变换法(BLT)可以设计IIR数字滤波器，但是……

  • 问题来了：这两种方法的输入都是一个"模拟滤波器 $H_a(s)$"，可是这个模拟滤波器从哪里来？
  • 课堂回顾：课堂讲过 $s = \frac{2}{T}\frac{z-1}{z+1}$ (BLT) 或 $h[n] = Th_a(nT)$ (IIM)，但 $H_a(s)$ 的设计公式太复杂了
  • 本实验目标：学会使用MATLAB快速设计模拟滤波器，完成IIR设计的"最后一块拼图"

在本实验中，我们将补充课堂未详细展开的模拟滤波器设计，并使用MATLAB完成完整的IIR滤波器设计流程。通过本实验，你将：

  • 理解为什么要"借用"模拟滤波器：直接设计数字IIR滤波器需要解非线性方程，而模拟滤波器已有百年成熟理论
  • 掌握三种经典模拟滤波器的特性：Butterworth（最平坦）、Chebyshev（等波纹）、Elliptic（最高效）
  • 学会使用MATLAB函数设计模拟原型：buttord、butter、cheby1、ellip等
  • 熟练应用双线性变换：将课堂学到的BLT公式用bilinear函数实现
  • 深入理解频率预翘曲：这是课堂公式背后的物理意义

复习知识衔接：从课堂到实验

2.1 你已经学会的知识

📚 课堂知识回顾

在理论课上，你已经学习了两种将模拟滤波器转换为数字滤波器的方法：

脉冲响应不变法 (IIM)

核心思想：让数字滤波器的脉冲响应等于模拟滤波器脉冲响应的采样值

$$h[n] = T \cdot h_a(nT)$$

优点：时域特性保持良好

缺点：有频率混叠问题

双线性变换法 (BLT)

核心思想：用代数变换将s平面映射到z平面

$$s = \frac{2}{T} \cdot \frac{z-1}{z+1}$$

优点：无频率混叠

缺点：有频率翘曲（需预翘曲补偿）

2.2 你需要补充的知识

❓ 缺失的一环：模拟滤波器 $H_a(s)$ 从哪里来？

课堂上，直接给出了模拟滤波器的传递函数，例如：

$$H_a(s) = \frac{1}{s^2 + 1.414s + 1} \quad \text{（2阶Butterworth低通）}$$

但实际工程中，我们需要根据设计指标（通带、阻带、衰减要求）来自动计算这个 $H_a(s)$。

问题：手工推导 Butterworth/Chebyshev/Elliptic 滤波器的传递函数需要：

  • 复杂的极点分布计算（涉及复数运算）
  • 切比雪夫多项式或椭圆函数的求解
  • 大量的代数运算，容易出错

解决方案：使用MATLAB的 butter、cheby1、ellip 函数自动完成！

2.3 本实验的完整流程

IIR滤波器设计完整流程 IIR滤波器设计完整流程（课堂 + 本实验） Step 1 设计指标 Fp, Fst, Rp, Rs ⭐ 本实验重点 Step 2: 模拟原型 • 预翘曲频率 • 计算阶次N • 设计Ha(s) 课堂已学 ✓ Step 3: BLT s → z 变换 bilinear() Step 4 数字滤波器 H(z) Step 5 验证与分析 freqz(), filter() 本实验核心任务： 学会使用 MATLAB 函数（buttord, butter, cheby1, ellip）自动设计模拟原型滤波器 这是课堂理论知识的实践补充，解决"模拟滤波器从哪里来"的问题

图1: IIR滤波器设计完整流程（橙色框为本实验重点）

理论经典模拟滤波器

3.1 为什么要"借用"模拟滤波器？

📻 历史背景：站在巨人的肩膀上

在数字信号处理诞生之前的100多年里，电气工程师们已经发展出了成熟的模拟滤波器理论：

  • 1930年：Stephen Butterworth 发表最大平坦滤波器理论
  • 1931年：Pafnuty Chebyshev 的等波纹逼近理论被应用于滤波器
  • 1958年：Wilhelm Cauer 完善了椭圆滤波器理论

这些理论经过几十年的打磨，已经有了完整的设计公式和查表法。当我们需要设计数字滤波器时，最聪明的做法就是：

  1. 先用成熟的理论设计一个满足指标的模拟滤波器
  2. 再用IIM或BLT方法转换成数字滤波器

这就是"借用"策略——借用前人的智慧成果！

3.2 三种经典滤波器的特性对比

🎯 用生活场景理解三种滤波器

假设你在设计一个声音过滤系统，要阻挡高频噪声，保留人声：

🔵 Butterworth（巴特沃斯）- "追求自然的完美主义者"

  • 就像高档音响追求"原汁原味"，通带内绝对平坦，不会改变任何保留的频率成分
  • 代价：过渡带较宽，"分界"不够干脆，可能漏过一些不想要的频率
  • 适用：音频处理（人耳对幅度变化敏感）、医疗信号处理

🟢 Chebyshev I（切比雪夫I型）- "务实的效率追求者"

  • 允许通带有轻微波动（可控制在0.5dB以内），换来更陡峭的过渡
  • 就像接受"略有瑕疵但性价比更高"的产品
  • 适用：通信系统（数字调制对幅度波动不敏感）

🔴 Elliptic（椭圆）- "极致效率的工程师"

  • 通带和阻带都允许波动，换来最窄的过渡带和最低的阶数
  • 同样指标，阶数可能只有Butterworth的一半！
  • 适用：嵌入式/实时系统（计算资源有限，需要最低复杂度）

特性对比表

滤波器类型 通带特性 阻带特性 过渡带 阶数效率 典型应用
Butterworth 最大平坦 ✓ 单调下降 最宽 最低 音频、医疗
Chebyshev I 等波纹 单调下降 较窄 中等 通信、雷达
Elliptic 等波纹 等波纹 最窄 ✓ 最高 ✓ 嵌入式、实时

3.3 幅度响应数学表达式

Butterworth 滤波器

$$|H_a(j\Omega)|^2 = \frac{1}{1 + (\Omega/\Omega_c)^{2N}}$$

特点：$\Omega=0$ 处前 $2N-1$ 阶导数为零 → "最大平坦"

Chebyshev I 型滤波器

$$|H_a(j\Omega)|^2 = \frac{1}{1 + \epsilon^2 T_N^2(\Omega/\Omega_p)}$$

其中 $T_N(x)$ 是N阶切比雪夫多项式，$\epsilon$ 控制通带波纹幅度

Elliptic 滤波器

$$|H_a(j\Omega)|^2 = \frac{1}{1 + \epsilon^2 R_N^2(\Omega/\Omega_p)}$$

其中 $R_N(x)$ 是N阶椭圆有理函数（雅可比椭圆函数）

💡 为什么不用手算这些公式？

以Butterworth为例，要从幅度响应推导传递函数 $H_a(s)$，需要：

  1. 计算极点位置：$p_k = \Omega_c \cdot e^{j\frac{\pi}{2N}(2k+N-1)}$，$k=1,2,...,N$
  2. 选择左半平面的极点（保证稳定）
  3. 构造传递函数：$H_a(s) = \frac{\Omega_c^N}{\prod_{k}(s-p_k)}$

Chebyshev和Elliptic的计算更加复杂，涉及双曲函数和椭圆积分。

结论：在实验中，我们使用MATLAB的 butter、cheby1、ellip 函数自动完成这些计算！

三种滤波器幅度响应对比 三种IIR滤波器幅度响应对比（相同设计指标） 频率 Ω Ωp (通带) Ωs (阻带) 幅度 |H(jΩ)| (dB) 0 -Rs 通带区域 阻带区域 滤波器对比 Butterworth Chebyshev I Elliptic ← 过渡带：Elliptic最窄

图2: 三种经典IIR滤波器的幅度响应对比

3.4 MATLAB设计函数速查

步骤 Butterworth Chebyshev I Elliptic
计算阶次 [N,Wn] = buttord(Wp,Ws,Rp,Rs,'s') [N,Wn] = cheb1ord(Wp,Ws,Rp,Rs,'s') [N,Wn] = ellipord(Wp,Ws,Rp,Rs,'s')
设计模拟原型 [num,den] = butter(N,Wn,'s') [num,den] = cheby1(N,Rp,Wn,'s') [num,den] = ellip(N,Rp,Rs,Wn,'s')
双线性变换 [numd,dend] = bilinear(num,den,Fs)
频率响应 [H,w] = freqz(numd,dend,N)

📌 函数参数说明

  • Wp, Ws: 预翘曲后的通带/阻带边界频率 (rad/s)
  • Rp: 通带最大衰减 (dB)，如 0.5 表示通带内最多衰减 0.5dB
  • Rs: 阻带最小衰减 (dB)，如 45 表示阻带至少衰减 45dB
  • 's': 表示设计模拟滤波器（s域），这是关键参数！
  • Fs: 采样频率 (Hz)

理论频率预翘曲：BLT的关键补充

4.1 问题的由来

在课堂上学习BLT时，你可能注意到了一个映射关系：

$$\omega = 2\arctan\left(\frac{\Omega T}{2}\right)$$

这意味着：模拟频率 $\Omega$ 和数字频率 $\omega$ 不是线性关系！

🗺️ 用地图投影理解频率翘曲

想象你要把一个球面（地球）展开成一个平面地图：

  • 靠近赤道的区域（低频）：变形较小，比较准确
  • 靠近极点的区域（高频）：被严重拉伸或压缩
  • 极点本身（$\omega = \pi$）：整个无限的高频区域被压缩到一个点！

BLT的频率映射就像这种投影：低频区域近似线性，高频区域严重"翘曲"。

4.2 预翘曲的解决方案

既然BLT会导致频率失真，我们就在设计模拟滤波器之前，先"反向变形"一下频率：

预翘曲公式：

$$\Omega_{prewarped} = \frac{2}{T}\tan\left(\frac{\omega}{2}\right) = 2F_s \cdot \tan\left(\frac{\pi F}{F_s}\right)$$

其中：$F$ 是目标数字频率 (Hz)，$F_s$ 是采样频率 (Hz)

✅ MATLAB 中的预翘曲实现

% 设计指标（数字域） Fs = 80000; % 采样频率 (Hz) Fp = 4000; % 通带边界 (Hz) Fst = 20000; % 阻带边界 (Hz) % Step 1: 转换为数字角频率 wp = 2*pi*Fp/Fs; % 通带数字角频率 (rad/sample) ws = 2*pi*Fst/Fs; % 阻带数字角频率 (rad/sample) % Step 2: 预翘曲 → 模拟角频率 Omega_p = (2*Fs) * tan(wp/2); % 通带模拟角频率 (rad/s) Omega_s = (2*Fs) * tan(ws/2); % 阻带模拟角频率 (rad/s) % 现在可以用 Omega_p 和 Omega_s 设计模拟滤波器了！
频率翘曲与预翘曲 双线性变换的频率翘曲与预翘曲补偿 ❌ 没有预翘曲 设计频率 实际位置↓ Ω (模拟) ω (数字) ✅ 使用预翘曲 预翘曲后 目标频率 ✓ Ω (模拟) ω (数字) 预翘曲原理： 先把目标频率"反向变形"，经过BLT的"正向变形"后，恰好回到正确位置

图3: 频率翘曲问题及预翘曲解决方案

实践交互演示：IIR滤波器频率响应

🎛️ 交互式滤波器设计演示

调整下面的参数，观察不同IIR滤波器的频率响应特性：

幅度响应

相位响应

实践示例讲解

完整示例：Butterworth 数字低通滤波器设计

设计要求：

  • 采样频率：$F_s = 80$ kHz
  • 通带边界：$F_p = 4$ kHz，最大衰减 $R_p = 0.5$ dB
  • 阻带边界：$F_{st} = 20$ kHz，最小衰减 $R_s = 45$ dB
%% 完整示例：Butterworth 数字低通滤波器设计 % 这个示例展示了从设计指标到最终数字滤波器的完整流程 %% Step 1: 定义设计指标（数字域） Fs = 80000; % 采样频率 (Hz) Fp = 4000; % 通带边界频率 (Hz) Fst = 20000; % 阻带边界频率 (Hz) Rp = 0.5; % 通带最大衰减 (dB) Rs = 45; % 阻带最小衰减 (dB) %% Step 2: 频率预翘曲（课堂BLT公式的实践应用） % 将数字频率转换为模拟频率，补偿双线性变换的频率畸变 wp = 2*pi*Fp/Fs; % 通带数字角频率 (rad/sample) ws = 2*pi*Fst/Fs; % 阻带数字角频率 (rad/sample) Omega_p = (2*Fs) * tan(wp/2); % 预翘曲后的通带模拟角频率 (rad/s) Omega_s = (2*Fs) * tan(ws/2); % 预翘曲后的阻带模拟角频率 (rad/s) fprintf('预翘曲结果:\n'); fprintf(' Omega_p = %.2f rad/s (原始: %.2f)\n', Omega_p, 2*pi*Fp); fprintf(' Omega_s = %.2f rad/s (原始: %.2f)\n', Omega_s, 2*pi*Fst); %% Step 3: 计算滤波器阶次（本实验核心 - 模拟滤波器设计） [N, Wn] = buttord(Omega_p, Omega_s, Rp, Rs, 's'); fprintf('滤波器阶次: N = %d\n', N); fprintf('截止频率: Wn = %.2f rad/s\n', Wn); %% Step 4: 设计模拟原型滤波器 [num_a, den_a] = butter(N, Wn, 's'); % 's' 表示模拟滤波器 fprintf('模拟滤波器 Ha(s) 设计完成\n'); %% Step 5: 双线性变换（课堂已学的BLT） [num_d, den_d] = bilinear(num_a, den_a, Fs); fprintf('数字滤波器 H(z) 设计完成\n'); %% Step 6: 验证与分析 [H, w] = freqz(num_d, den_d, 2048); f = w * Fs / (2*pi); % 转换为Hz % 绘制幅度响应 figure('Position', [100, 100, 900, 400]); plot(f, 20*log10(abs(H)), 'b-', 'LineWidth', 2); hold on; % 标记设计指标 plot([Fp Fp], [-80 5], 'g--', 'LineWidth', 1.5); % 通带边界 plot([Fst Fst], [-80 5], 'r--', 'LineWidth', 1.5); % 阻带边界 plot([0 Fp], [-Rp -Rp], 'g:', 'LineWidth', 1.5); % 通带衰减限 plot([Fst Fs/2], [-Rs -Rs], 'r:', 'LineWidth', 1.5); % 阻带衰减限 grid on; xlabel('频率 (Hz)'); ylabel('幅度 (dB)'); title(sprintf('Butterworth 低通滤波器 (N=%d)', N)); legend('幅度响应', '通带边界', '阻带边界', 'Location', 'southwest'); xlim([0 Fs/2]); ylim([-80 5]);

🔍 代码解析：对应课堂知识

代码步骤 对应课堂知识 说明
Step 2: 预翘曲 BLT的频率映射 $\omega = 2\arctan(\Omega T/2)$ 反向使用映射公式，补偿畸变
Step 3-4: 模拟滤波器 本实验新学内容 使用MATLAB自动设计 $H_a(s)$
Step 5: bilinear BLT公式 $s = \frac{2}{T}\frac{z-1}{z+1}$ MATLAB自动完成变换

实践练习一：Butterworth 数字低通滤波器设计

任务

任务描述

设计一个Butterworth数字低通滤波器，满足以下要求：

  • 采样频率：$F_s = 80$ kHz
  • 通带边界：$F_p = 4$ kHz，最大衰减 $R_p = 0.5$ dB
  • 阻带边界：$F_{st} = 20$ kHz，最小衰减 $R_s = 45$ dB
  • 使用双线性变换方法

任务要求

  1. 按照示例的6个步骤完成设计
  2. 计算并输出滤波器阶次 $N$ 和截止频率 $W_n$
  3. 绘制幅度响应图，标记通带边界、阻带边界和设计指标线
  4. 验证设计结果是否满足指标要求

预期结果

滤波器阶次 N = 4，通带平坦（最大衰减 < 0.5dB），阻带衰减 > 45dB

实践练习二：Chebyshev I 型数字低通滤波器设计

任务

任务描述

使用与练习1相同的设计指标，设计一个Chebyshev I型数字低通滤波器：

  • 采样频率：$F_s = 80$ kHz
  • 通带边界：$F_p = 4$ kHz，最大衰减 $R_p = 0.5$ dB
  • 阻带边界：$F_{st} = 20$ kHz，最小衰减 $R_s = 45$ dB

任务要求

  1. 使用 cheb1ord 和 cheby1 函数
  2. 注意：cheby1 需要额外的 $R_p$ 参数
  3. 绘制幅度响应并与练习1对比
  4. 回答以下问题：
    • Q1: 阶数是否比Butterworth更低？低多少？
    • Q2: 通带是否有波纹？波纹幅度大约是多少dB？
    • Q3: 过渡带是否更陡峭？

实践练习三：Elliptic 数字低通滤波器设计

任务

任务描述

继续使用相同的设计指标，设计一个Elliptic（椭圆）数字低通滤波器。

任务要求

  1. 使用 ellipord 和 ellip 函数
  2. 注意：ellip 需要同时指定 $R_p$ 和 $R_s$ 参数
  3. 制作综合对比表格：

    📌 根据练习1-3的实验结果填写下表（使用相同的设计指标）

    特性 Butterworth Chebyshev I Elliptic
    滤波器阶数 N ? ? ?
    通带特性 ? ? ?
    阻带特性 ? ? ?
    过渡带宽度 ? ? ?

实践练习四：实际应用中的滤波器选择

任务 (思考题)

场景1：音频均衡器设计

你正在为一款音乐播放App设计均衡器。采样率44.1 kHz，需要一个低通滤波器消除20 kHz以上噪声。

  • Q1: 你会选择哪种滤波器？为什么？
  • Q2: 通带波纹对音质有什么影响？
  • Q3: 相位失真对音质有什么影响？如何减轻？

场景2：通信信道选择

需要从多个相邻信道中选出900-910 kHz的目标信道，抑制相邻信道干扰。

  • Q1: 这种场景更适合哪种滤波器？
  • Q2: 如果Butterworth需要12阶，Elliptic只需6阶，你如何选择？
  • Q3: 设计带通滤波器需要哪些MATLAB函数？

场景3：嵌入式实时处理

在手机芯片上实现实时语音滤波，计算资源有限。

  • Q1: 滤波器阶数与计算复杂度的关系是什么？
  • Q2: 阶数从12降到6意味着什么？（计算量、功耗、延迟）
  • Q3: 如何在实时性和滤波效果之间权衡？

实验报告要求

📝 报告内容

  1. 实验目的和原理：说明本实验如何补充课堂知识
  2. 设计过程：详细记录每个练习的设计步骤和MATLAB代码
  3. 结果与分析：三种滤波器的对比表格、频率响应图
  4. 思考题回答：对练习4的问题给出分析
  5. 结论：总结不同应用场景下如何选择滤波器

附录脉冲响应不变法（IIM）补充说明

A.1 为什么本实验主要使用BLT而非IIM？

在课堂上，你学习了两种将模拟滤波器转换为数字滤波器的方法。本实验选择双线性变换法(BLT)作为主要方法，原因如下：

比较项 脉冲响应不变法 (IIM) 双线性变换法 (BLT)
频率混叠 ❌ 有混叠问题（高频分量会折叠到低频） ✅ 无混叠（整个模拟频率轴映射到 $-\pi$ ~ $\pi$）
频率映射 ✅ 线性映射 $\omega = \Omega T$ ⚠️ 非线性映射（需要预翘曲补偿）
适用滤波器类型 ❌ 仅低通、带通（高通和带阻有问题） ✅ 所有类型（低通、高通、带通、带阻）
时域特性保持 ✅ 脉冲响应完全对应 ⚠️ 脉冲响应有差异
工程应用 较少使用 主流方法

A.2 IIM的核心原理回顾

基本思想：让数字滤波器的脉冲响应等于模拟滤波器脉冲响应的采样值

$$h[n] = T \cdot h_a(nT)$$

其中 $T = 1/F_s$ 是采样周期

频域映射：

$$H(e^{j\omega}) = \sum_{k=-\infty}^{\infty} H_a\left(j\frac{\omega - 2\pi k}{T}\right)$$

这个求和表明：模拟滤波器的频率响应会以 $2\pi/T$ 为周期产生混叠

A.3 IIM的频率混叠问题图示

脉冲响应不变法的频率混叠问题 原始模拟滤波器 Ha(jΩ) Ω IIM转换后 H(ejω) - 有混叠! ω 混叠! 混叠! IIM ⚠️ IIM的致命缺陷：频率混叠 当模拟滤波器的带宽超过 Fs/2 时，高频分量会"折叠"到低频区域，造成失真 因此 IIM 只适用于带限的低通/带通滤波器，且需要足够高的采样率

图A1: IIM的频率混叠问题示意图

A.4 MATLAB中的IIM实现

📝 MATLAB函数：impinvar

虽然本实验不要求使用IIM，但如果你想尝试，MATLAB提供了 impinvar 函数：

% 使用脉冲响应不变法（仅供参考，本实验不要求） [num_d_iim, den_d_iim] = impinvar(num_a, den_a, Fs); % 对比 BLT 方法 [num_d_blt, den_d_blt] = bilinear(num_a, den_a, Fs);

A.5 何时考虑使用IIM？

🤔 IIM的适用场景

虽然BLT是主流方法，但在某些特殊情况下IIM可能有优势：

  • 需要精确的时域响应：例如模拟某个特定的模拟系统
  • 低通滤波器 + 高采样率：混叠可忽略不计
  • 线性频率映射更重要：不希望频率翘曲影响设计

但在大多数工程应用中，BLT因其无混叠和通用性而更受欢迎。