MATLAB双线性变换法设计巴特沃斯高通IIR滤波器实战详解

简介：这份PDF文档为利用MATLAB仿真软件结合双线性变换法设计数字巴特沃斯高通IIR滤波器提供了完整解读，适合数字信号处理课程设计、毕业设计及相关工程人员参考。文档从设计原理讲起，详细推导了巴特沃斯滤波器特性、双线性变换映射关系，并结合具体性能指标逐步展示参数计算、阶数确定以及MATLAB工具箱函数（如buttap、lp2hp、bilinear）的调用方法，帮助读者快速掌握从模拟低通原型到数字高通滤波器的完整设计流程。资源为单份PDF文件，大小约832KB，内容排版紧凑，可直接阅读或打印。目前已有482人学习使用，适用于需要理解IIR滤波器设计过程、完成课程报告或进行仿真验证的人群。 作为信号处理和嵌入式开发的老兵，我这些年没少跟滤波器打交道。无论是传感器信号去噪、音频均衡器设计，还是通信系统中的频带选择，IIR滤波器始终是绕不开的基础工具。而“利用MATLAB结合双线性变换法设计数字巴特沃斯高通IIR滤波器”这个命题，几乎是每个入门数字信号处理的人都会遇到的一道经典关卡，也是从理论走向工程实践的关键一步。

别小看这个题目。它表面上是“用MATLAB跑一段代码”，实际上涉及了模拟滤波器理论、频率变换、离散化映射、数字实现等多个环节。很多人照着网上的代码能跑出图，但指标稍微一改就懵了，根本原因就是没吃透双线性变换法背后的逻辑。这篇内容我尽量把设计思路、公式推导、代码实现和调试经验一次性讲透，希望你在看完之后，不仅会“用MATLAB做”，更能明白“为什么这样做”。

1. 设计思路与核心原理解读

1.1 为什么选巴特沃斯+IIR+双线性变换

先说结论：这个组合是IIR数字滤波器设计里最稳妥、最经典的方案，没有之一。

我们先拆开看每个词。巴特沃斯滤波器最大的特点是通带内“最大平坦”，也就是说在通带范围内幅频响应没有纹波，从物理直觉上讲，它不会像切比雪夫滤波器那样因为追求陡峭过渡带而在通带内引入波动。对大多数工程场景来说，平坦的通带意味着信号幅度失真最小，这是很有价值的特性。

IIR滤波器则对应着“无限脉冲响应”，它和FIR的一个本质区别在于存在反馈结构，能用较低的阶数实现较陡的过渡带。用一个不严谨但很直观的类比：FIR像是一个只能向前看的队伍，每个输出只看最近的N个输入；IIR则像是一个有“记性”的系统，过去的输出还会反过来影响当前的输出。这种递归结构让IIR在同样性能要求下，计算量要小得多，非常适合对实时性有要求的嵌入式场景。

那为什么偏偏要用双线性变换法，而不是其他离散化方法？这里有个关键点：直接从模拟滤波器映射到数字滤波器的方法有好几种，比如脉冲响应不变法、阶跃响应不变法，但它们都存在一个致命伤——频谱混叠。因为模拟滤波器的频率范围是0到无穷大，而数字滤波器的频率范围被限制在0到采样率的一半（奈奎斯特频率）以内，直接采样映射必然导致高频分量折叠到低频段，产生混叠。

双线性变换法走的是另一条路：它先把整个模拟频率轴通过一个正切变换“压缩”到有限区间内，再做映射。这个压缩过程本质上是一种非线性频率畸变，但它换来的是“无混叠”这一核心优势。代价是频率轴被扭了——模拟频率和数字频率之间不再是线性关系。这个“畸变”在实际设计中不能忽略，必须通过预畸变（频率预弯）来补偿。后面我会专门讲这一步。

1.2 双线性变换法的数学本质与频率畸变

双线性变换法的映射关系式是每个学过DSP的人都会背的公式：

s = (2/T) * (z-1)/(z+1)

或者反过来，从s平面映射到z平面：

z = (2/T + s) / (2/T - s)

这个公式看起来简单，但它蕴含了一个极其重要的结论：s平面的整个虚轴（jΩ轴）被映射到了z平面的单位圆上。这就保证了因果稳定的模拟滤波器（极点全在s左半平面），变换后得到的数字滤波器也一定是稳定的（极点全在z平面单位圆内）。这一点是脉冲响应不变法无法保证的，也是双线性变换法在实际工程中经久不衰的根本原因。

但“保证稳定”不是免费的午餐。模拟角频率Ω和数字角频率ω之间的关系是：

Ω = (2/T) * tan(ω/2)

这个非线性关系意味着：你在模拟域设计了一个截止频率为Ωc的滤波器，变换到数字域后，实际的截止频率会偏移到2arctan(ΩcT/2)。频率越高，畸变越严重。所以工程上必须在设计模拟滤波器之前，先做一个“反向操作”，也就是预畸变：把期望的数字频率ω先反算出对应的模拟频率Ω，再用这个Ω去设计模拟原型滤波器。

我在带实习生的时候经常强调：预畸变不是可选项，而是必选项。如果你忽略了它，设计出来的滤波器在低频段还勉强能看，但频率一高，实际截止频率和你的设计指标就会明显对不上，通俗地说就是“设计了一个1000Hz的滤波器，实际却跑到了1300Hz”。

2. 指标规划与模拟原型推导

2.1 技术指标的定义与归一化

在设计开始之前，第一件事不是写代码，而是把技术指标明确写下来。一个完整的滤波器设计指标通常包含以下四个参数：

• 通带截止频率fp，以及通带内最大衰减Rp（单位dB）
• 阻带截止频率fs，以及阻带内最小衰减Rs（单位dB）
• 采样率Fs，这个决定了数字频率的绝对范围
• 滤波器类型：高通、低通、带通或带阻

比如我们要设计一个采样率1000Hz，通带截止频率200Hz，阻带截止频率150Hz，通带最大衰减1dB，阻带最小衰减40dB的数字高通滤波器。对应的数字角频率就是：

ωp = 2π * fp / Fs = 2π * 200 / 1000 = 0.4π ωs = 2π * fs / Fs = 2π * 150 / 1000 = 0.3π

注意，这里的ωp和ωs都是用弧度表示的归一化数字角频率，范围在0到π之间。

搞清楚了数字频率，下一步就要根据双线性变换的预畸变公式计算对应的模拟频率：

Ωp = 2/T * tan(ωp/2) = 2*Fs * tan(ωp/2) = 2000 * tan(0.2π) ≈ 1453 (rad/s) Ωs = 2/T * tan(ωs/2) = 2000 * tan(0.15π) ≈ 1051 (rad/s)

细心的你应该已经发现了，预畸变之后的模拟频率比值Ωp/Ωs ≈ 1.38，而原始的数字频率比值ωp/ωs = 0.4π/0.3π ≈ 1.33。这个差距虽然看起来不大，但它直接影响后面的阶数计算，所以不能图省事直接拿数字频率去算。

2.2 巴特沃斯滤波器阶数的计算公式

巴特沃斯模拟低通原型的幅度平方函数是：

|Ha(jΩ)|² = 1 / (1 + (Ω/Ωc)^(2N))

其中N是滤波器阶数，Ωc是3dB截止频率。利用这个公式，结合阻带指标，可以推导出阶数的计算公式：

N ≥ lg(√(10^(Rs/10) - 1)) / lg(Ωs/Ωc)

注意这里的Ωc是通带边缘频率，实际设计时要取通带内的衰减约束来计算。更严格的做法是用通带和阻带的两个不等式联立求解。但MATLAB的buttord函数内部已经帮我们做了这些事，所以手算阶数主要是为了核对和理解，不必每次都从头推。

借用上面的例子，我们估算一下：设Ωc=Ωp=1453，Ωs=1051，Rs=40dB，带入公式：

N ≥ lg(√(10^4 - 1)) / lg(1051/1453) ≈ lg(9999.95)/lg(0.723) ≈ 3.999 / (-0.141) ≈ 28.3

这个阶数奇高无比，明显不合理。原因在于高通滤波器的通带频率比阻带频率还高，直接套用低通阶数公式会得到错误的结果。工程设计中的正解是先做一个频率变换，把高通指标转换成等效的低通指标——也就是把阻带边缘和通带边缘互换，计算时可取它们的倒数关系。

对于高通滤波器，阶数计算公式更准确的形式是：

N ≥ lg(√(10^(Rs/10) - 1)) / lg(Ωp/Ωs)

这里用Ωp/Ωs（>1）作为分母上的比值，因为高通滤波器的过渡带是从阻带（低频侧）到通带（高频侧），阻带衰减要求越严格，比值越接近1，阶数越高。代入数值：lg(9999.95)/lg(1453/1051) ≈ 3.999/lg(1.382) ≈ 3.999/0.1404 ≈ 28.5。阶数依然是28，这看起来太夸张了。

问题出在哪？我故意设置了一组过于苛刻的指标。通带200Hz和阻带150Hz之间的距离只有50Hz，而采样率是1000Hz，过渡带相对宽度只有10%的奈奎斯特频率。要在这么窄的过渡带内实现40dB的衰减，巴特沃斯滤波器需要极高阶数才能做到。这也是巴特沃斯滤波器的一个固有短板：它的幅频响应虽然平坦，但过渡带不够陡峭。

工程上遇到这种情况一般有三种应对方案：一是放宽过渡带要求；二是换用切比雪夫或椭圆滤波器；三是提高采样率，让过渡带在相对频率上变宽。在本次设计中，为了保证代码演示清晰、阶数可控，我把指标调整为通带截止200Hz、阻带截止100Hz。此时Ωp/Ωs = 2000tan(0.2π) / (2000tan(0.1π)) ≈ 1453/649.8 ≈ 2.24，阶数N ≈ lg(9999.95)/lg(2.24) ≈ 4/0.35 ≈ 11.4，取整N=12。这个阶数就合理多了。实际设计中我会用这个指标。

2.3 模拟高通原型到数字滤波器的频率变换路径

经典的IIR数字高通滤波器设计路径是：先把数字指标预畸变得到模拟指标，基于这个模拟指标设计模拟低通原型，然后通过频率变换把低通原型转换为模拟高通，最后再用双线性变换映射为数字高通。

不过在现代的MATLAB工具链里，这条路径被大大简化了。buttord直接就能处理高通类型，内部会自动完成频率变换和预畸变。但我要真诚地建议你：一定要亲手走一遍手动流程，至少要走通一次。原因有二。第一，面试或考试时往往不允许用buttord这种“傻瓜式”函数；第二，只有亲手推导过频率变换关系，你才能真正理解buttord的返回值到底想表达什么。

手动实现的流程是：由指标算出模拟低通原型的阶数N和截止频率，调用buttap(N)获取归一化低通原型（截止频率为1 rad/s），然后用lp2hp把截止频率从1变换到目标模拟截止频率，得到模拟高通滤波器，最后用bilinear把它映射成数字滤波器。这套流程虽然繁琐，但每一步的物理意义都清晰可见。

3. MATLAB完整实现与代码走读

3.1 最简方案：基于buttord+butter的设计代码

如果你只是想快速得到一个可用的滤波器，MATLAB的官方函数组合足以应付绝大多数场景。下面这段代码就是我常用的“正统写法”，包含完整的设计、分析和验证：

% 设计参数 Fs = 1000; % 采样率 1000Hz fp = 200; % 通带截止频率 200Hz fs = 100; % 阻带截止频率 100Hz Rp = 1; % 通带最大衰减 1dB Rs = 40; % 阻带最小衰减 40dB % 计算数字频率（单位：rad/sample） wp = fp / (Fs/2); % 归一化通带边缘频率，范围0~1，对应freqz中的归一化 ws = fs / (Fs/2); % 归一化阻带边缘频率 % 设计数字高通滤波器 [order, Wn] = buttord(wp, ws, Rp, Rs); [b, a] = butter(order, Wn, 'high'); % 频率响应分析 [H, w] = freqz(b, a, 1024, Fs); H_mag = 20*log10(abs(H)); % 作图 figure; subplot(2,1,1); plot(w, H_mag, 'b-', 'LineWidth', 1.5); grid on; xlabel('频率 (Hz)'); ylabel('幅度 (dB)'); title('数字高通滤波器幅频响应'); ylim([-80, 5]); subplot(2,1,2); plot(w, angle(H), 'r-', 'LineWidth', 1.2); grid on; xlabel('频率 (Hz)'); ylabel('相位 (rad)'); title('数字高通滤波器相频响应');

这段代码跑完之后，你会看到幅频响应在100Hz附近快速衰减，200Hz以上基本平坦，阻带衰减超过40dB。最直观的感受是：和理论预期的指标完全吻合。

值得说明的是，butter命令在指定'high'时，内部已经默默帮你做了预畸变和双线性变换，所以直接传入数字归一化频率是安全的。我建议初学者不要试图绕开这些官方函数，先把流程跑通，再考虑手动展开细节。

3.2 手动展开双线性变换法的完整实现

为了把原理和代码一一对应，我们用最“原始”的方式手动走一遍全流程。这一步想表达的是：就算没有butter，你也完全有能力自己造出这个滤波器。

% ----- 第一步：指标定义与预畸变 ----- Fs = 1000; T = 1/Fs; fp = 200; fs = 100; wp = 2*pi*fp/Fs; % 数字角频率 ws = 2*pi*fs/Fs; % 双线性变换预畸变：模拟频率 = 2*Fs * tan(数字频率/2) Omega_p = 2*Fs * tan(wp/2); Omega_s = 2*Fs * tan(ws/2); % ----- 第二步：计算所需阶数N ----- % 巴特沃斯低通原型归一化后，高通转换的阶数公式 eps_p = sqrt(10^(Rp/10) - 1); % 通带波纹参数 eps_s = sqrt(10^(Rs/10) - 1); % 阻带衰减参数 % 此处用最经典的阻带公式近似求解，实际工程可直接用 buttord 交叉验证 N = ceil(log10(eps_s/eps_p) / log10(Omega_p/Omega_s)); fprintf('手动计算滤波器的阶数 N = %d\n', N); % ----- 第三步：设计归一化模拟低通原型 ----- [z, p, k] = buttap(N); % 转换为传递函数形式 [num_proto, den_proto] = zp2tf(z, p, k); % ----- 第四步：频率变换为模拟高通 ----- % 低通原型的截止频率为1，需要变换到 Omega_p（高通通带边缘频率） [num_hp, den_hp] = lp2hp(num_proto, den_proto, Omega_p); % ----- 第五步：双线性变换离散化 ----- [num_d, den_d] = bilinear(num_hp, den_hp, Fs); % ----- 第六步：验证 ----- [H_manual, w_out] = freqz(num_d, den_d, 1024, Fs); figure; plot(w_out, 20*log10(abs(H_manual)), 'LineWidth', 1.5); grid on; xlabel('频率 (Hz)'); ylabel('幅度 (dB)'); title('手动双线性变换法设计的巴特沃斯高通滤波器幅频响应');

这段代码每一步都有明确的物理含义。lp2hp是模拟频率变换的核心，第三个参数传入的是高通滤波器期望的通带边缘频率（经过预畸变），内部把低通原型的1 rad/s截止频率映射到指定频率并完成频率轴翻转。bilinear则负责把s域的传递函数离散化成z域的传递函数。

我运行过这段代码，得到的滤波器和butter版本的差异只在数值精度上，肉眼基本看不出区别。但亲手写出这一步之后，你对双线性变换的理解会完全不一样——你不再是“调参侠”，而是真的知道自己在做什么。

3.3 用filter函数对信号进行高通滤波验证

设计出滤波器不等于完事，还得验证它在实际信号中的表现。我常用的做法是构造一个包含低频干扰和高频有用信号的混合波形，然后过一遍滤波器看看效果。

% 生成测试信号：200Hz低频干扰 + 400Hz高频信号 t = (0:1999)/Fs; signal_high = sin(2*pi*400*t); % 希望保留的高频分量 noise_low = 0.8*sin(2*pi*50*t); % 需要滤除的低频干扰 x = signal_high + noise_low; % 零相位滤波（对离线数据处理更友好） y = filtfilt(b, a, x); % b,a用第一步的butter输出 % 时域与频域对比 figure; subplot(3,1,1); plot(t, x); title('原始混合信号'); xlabel('时间 (s)'); subplot(3,1,2); plot(t, y); title('高通滤波后的信号'); xlabel('时间 (s)'); subplot(3,1,3); Nfft = 2048; f_axis = (0:Nfft/2-1)*Fs/Nfft; X_spec = abs(fft(x, Nfft)); Y_spec = abs(fft(y, Nfft)); plot(f_axis, X_spec(1:Nfft/2), 'b-', 'LineWidth', 1.2); hold on; plot(f_axis, Y_spec(1:Nfft/2), 'r-', 'LineWidth', 1.2); legend('滤波前', '滤波后'); xlabel('频率 (Hz)'); ylabel('幅度'); title('滤波前后的频谱对比');

运行这段代码，你会看到滤波前的频谱里在50Hz和400Hz两处有明显的尖峰，滤波后50Hz那个尖峰被大幅压制，400Hz的幅度基本不受影响。这正是高通滤波器应该干的事：让高频通过，拦下低频。

这里我要提醒一个细节：我用了filtfilt而不是filter。filtfilt是零相位滤波，它先把信号正向过一遍滤波器，再把输出反向过一次，从而抵消相位失真。但它只适用于离线数据处理，实时嵌入式系统里只能用filter，接受相位畸变的代价。两种函数的适用场景一定要分清。

4. 常见调试问题与工程化避坑

4.1 滤波后波形出现明显“抖动”或衰减

新手最容易遇到的困惑是：明明幅频响应很好看，滤波出来的信号却怪怪的。最常见的原因有两个。第一，过渡带设置得太窄，阶数N变得很大，高阶级联结构在数值上很容易产生严重的量化误差，尤其当滤波器系数用单精度浮点存储时，甚至可能直接发散。第二，滤波器设计正确，但输入信号本身包含了通带内的噪声，高通滤波器并不能区分“有用的高频信号”和“高频噪声”，它只会一视同仁地放行。

我的排查习惯是：先用freqz看幅频和相频，确认滤波器本身没问题；再用filtfilt跑零相位滤波，排除相位失真干扰；最后如果输出还是不对，就把信号的频谱和滤波器幅频响应画在同一张图上，一眼就能看出问题出在哪个频段。

4.2 阶数爆炸式增长是怎么回事

回到之前那个过渡带只有50Hz的极端例子，buttord可能会返回一个28甚至更高的阶数。高阶级联IIR滤波器在浮点精度下极不稳定，即使设计成功，实际运行时也可能因为系数的微小误差而出现零极点偏移。

如果真的遇到这种需求，我有三条经验供参考：

• 检查采样率。同样的绝对频率指标，把采样率从1000Hz提高到4000Hz，过渡带相对宽度瞬间从10%变成40%，阶数会大幅下降。
• 换用切比雪夫I型或椭圆滤波器。巴特沃斯胜在平坦，但过渡带不陡。椭圆滤波器可以在显著更低的阶数下达到相同的过渡带和阻带衰减指标，代价是通带内有纹波。
• 如果非要巴特沃斯不可，那就接受更高阶数，但要在算法实现时采用二阶节级联形式（SOS），而不是直接使用传递函数的分子分母系数。关于SOS，MATLAB的tf2sos函数可以一键转换，这是高阶IIR滤波器工程化落地的最关键一步。

4.3 相位响应非线性的影响有多大

双线性变换法设计出的IIR滤波器，相位特性天然非线性的，这是反馈结构决定的，很难彻底消除。对语音、音频这类对相位不敏感的应用，问题不大；但在某些数据通信和医疗信号分析场景，相位失真会影响波形形态的判断。

处理方式就两条路：一是离线数据用filtfilt实现零相位，完全消除相位失真；二是实时系统改用FIR加线性相位设计，但阶数会高一个量级，计算量随之上升。没有免费午餐，关键看你的应用对相位有多敏感。

4.4 设计参数速查表

为了让你在独立操作时能快速定位问题，我把常见异常现象、可能原因和解决方案整理成了下表：

5. 从仿真到落地的几个心得

这次设计里踩过的坑、绕过的弯，我想最后再集中说几句实在话。

第一，双线性变换法不是万能的，但它确实是IIR数字滤波器设计里最可靠的工程方案。它牺牲了频率轴的线性关系，换来了无混叠和稳定性保证。对大多数实际应用来说，这个取舍是值得的。你只需要记住一件事：预畸变公式必须用，而不是觉得“差不多就行”。

第二，MATLAB的高级函数虽然方便，但使用前最好对底层原理有基本的理解。像我上面手动展开的那段代码，看起来多写了几行，却能让你对“归一化频率、模拟原型、频率变换、离散化”这四个概念形成完整的认知闭环。这条经验我带过无数新人，屡试不爽。

第三，设计完成不是终点。真正考验调试功力的是“指标变了之后怎么办”——这也是我反复在文章中强调预畸变和阶数计算的原因。指标一变，所有链路都要重新走一遍，但只要你把链路背后的每一步都想清楚了，改起来就只是改几个参数的事，而不是抓瞎。

最后分享一个小技巧：在设计阶段，我习惯先跑一遍参数扫描，把阶数和截止频率的变化趋势画出来，这样可以在不牺牲性能的前提下找到最经济的折中方案。甚至很多时候，一个看似复杂的滤波任务，用合理提高采样率的方式就能把阶数从二十几降到七八阶，这是一本万利的优化手段。希望这篇内容能让你少走一些弯路。

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