当前位置：首页 > news >正文

音频进阶学习二十五——脉冲响应不变法实现低通滤波器

news 来源：原创 2025/6/27 16:55:53

文章目录

前言
一、脉冲响应不变法
- 1.定义
- 2.模拟系统冲激响应的周期采样
- 3.模拟系统和数字系统的频域响应关系
- - 1）S域和Z域的关系
  - 2）幅频响应的关系
- 4.通过有理函数设计滤波器
- 5.总结
二、低通滤波器的设计实例
- 1.给定数字滤波器指标
- 2.转换模拟滤波器指标
- 3.模拟滤波器的设计
- 4.模拟滤波器的传递函数
- 5.转换为数字系统的传递函数
三、使用matlab实现
- 1.通过传递函数实现
- 2.使用matlab工具
总结

前言

上一篇文章中，我们介绍了滤波器的设计流程，并且提到了从模拟滤波器传递函数到数字滤波器传递函数的两种转变方法：脉冲响应不变法、双线性转换法。

本篇文章将使用脉冲响应不变法设计出一个低通滤波器。在文章的开始，将会对于脉冲响应不变法原理、频域、幅频响应做一个深入解析。然后根据给定的指标，通过巴特沃斯滤波器得到数字滤波器的传递函数，并且在计算机中进行低通滤波的验证。

|版本声明：山河君，未经博主允许，禁止转载

一、脉冲响应不变法

1.定义

保持模拟系统的冲激响应不变，而离散化后的数字滤波器具有与原模拟滤波器相同的冲激响应。即 $h_a[t]$ 是模拟滤波器的单位冲激响应，而 $h[n]=h_a[nT]$ ，其中 $T$ 为采样周期。从而实现 $H_c(s)$ 到 $H (z)$ 的转换。

整体过程如下：
在这里插入图片描述

注意，梳理下文中所有用到的公式：

$H_c(s)=\int_{-\infty}^{\infty}h_a(t)e^{-st}dt$ ：连续时间（模拟）系统的传递函数
$H_c(s)=\sum_{k=1}^N\frac{A_k}{s-s_k}$ ：模拟系统的传递函数为以部分分式展开的形式
$H_a(j\Omega)=H_c(s)|_{s=j\Omega}=\int_{-\infty}^{\infty}h_a(t)e^{-j\Omega t}dt$ ：传递函数的傅里叶变换，即模拟系统的幅频响应
$H(e^{j\omega})=\sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}$ ：离散系统的幅频响应
$H(z)=\sum_{n=-\infty}^{\infty}x[n]z^{-n}$ ：离散时间（数字）系统的传递函数
$H_s(s)=\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}h_a(t)\delta(t-nTs)e^{-st}dt$ ：离散时间系统的拉普拉斯变换（下文中推导）

接下来，我们就具体分析其中的步骤

2.模拟系统冲激响应的周期采样

由上文可知，首先设计利用冲激串 $s (t)$ 进行周期采样，采样周期为 $T_s$ ：
$s(t)=\sum_{n=-\infty}^{\infty}\delta(t-nTs)$
再对于模拟滤波器冲激响应进行采样：
$h[n]=h_a(t)|_{t=nT_s}=h_a[nT_s]=h_s[nT_s],\quad h_s(t)=h_a(t)s(t)$
由上文可知：连续时间信号的系统的传递函数，以及连续时间系统频率响应分别为
$H_c(s)=\int_{-\infty}^{\infty}h_a(t)e^{-st}dt, \quad H_a(j\Omega)=H_c(s)|_{s=j\Omega}$
其中， $s$ 是复数，为 $\sigma+j\Omega$ ，由此对于模拟系统的幅频响应的采样为：
$H_s(s)=\int_{-\infty}^{\infty}h_a(t)\sum_{n=-\infty}^{\infty}\delta(t-nTs)e^{-st}dt$
对于上式同样有线性性，那么上式为：
$H_s(s)=\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}h_a(t)\delta(t-nTs)e^{-st}dt, \quad t=nT_s\\ =\sum_{n=-\infty}^{\infty}h_a(nT_s)e^{-snT_s}$
此时对比Z变换 $H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}$ ，即当 $z=e^{sT_s}$ 时，两者相等，此时 $H(z)=H_s(s)|_{s=\ln(z)/T_s}$ 。这就是脉冲响应不变法的原理推导。

3.模拟系统和数字系统的频域响应关系

1）S域和Z域的关系

上文中说当 $z=e^{sT_s}$ 时， $H(z)=H_s(s)|_{s=\ln(z)/T_s}$ ，我们知道：

$s=\sigma+j\Omega$
$z=re^{j\omega}$

则将上式转为实部与虚部 $re^{j\omega}=e^{(\sigma+j\Omega)T_s}=e^{\sigma T_s}e^{j\Omega T_s}$ ，此时就有：

$r=e^{\sigma T_s}$ ，决定了模值
$\omega=\Omega T_s$ ，决定了角度

然而由于指数函数在虚部方向是周期性的 $e^{sT_s+2\pi kT_s}=e^{sT_s} * 1$ ，即对于不同频率分量 $k$ 实际上映射到了同一点，这也就是多重映射：

$\sigma=0$ 时， $∣ z ∣ = 1$ ，此时S域中虚轴就变成了Z域上的单位圆
$\sigma<0$ 时， $∣ z ∣ < 1$ ，此时S域就映射到了单位圆内部，此时为因果稳定系统
$\sigma>0$ 时， $∣ z ∣ > 1$ ，此时S域就映射到了单位圆外

2）幅频响应的关系

从连续时间傅里叶变换CTFT到DTFT，即描述连续时间系统频率响应和离散时间系统频率响应关系存在以下公式：
$H(e^{j\omega})=\frac{1}{T_s}\sum_{k=-\infty}^{\infty}H_a(j\frac{\omega}{T_s}-j\frac{2\pi k}{T_s})$

混叠
根据采样定理：采样过程实际上是在时域对信号进行周期性地重复。采样后的信号的频域会变成周期性的。也就是在频域上，信号会以采样频率 $f_s=1/T_s$ 为周期出现重复。此时会发生：
- 连续时间的频率响应基本上是非带限的信号，具有无限带宽，所以在S域映射到Z域（如上图）时，高频信号会折叠到低频信号，但是这种折叠在低频时可以忽略不计
- 如果冲激响应的的最高频率很高，但是采样频率低，这种重复可能会导致幅频响应 $H(e^{j\omega})$ 的混叠，如下图

在这里插入图片描述
所以冲激响应不变法只适用于低通滤波器。

增益
我们知道，在实际应用中，语音是先通过模数转换为数字信号处理，然后通过数模转为为模拟信号播出，在数模转换时现在给出增益为 $T_s$ 。在滤波器中，我们知道到输出的频谱是输入的频谱乘上冲激响应频谱：
$Y(e^{j\omega})=X(e^{j\omega})H(e^{j\omega})$
而根据上式我们知道 $X(e^{j\omega})H(e^{j\omega})$ 对于频谱缩小了 $1/T_s^2$ ，即使数模转换时的增益时 $T_s$ 依旧缩小了 $1/T_s$ ，所以在实际对于连续时间冲激响应的周期采样实际上是对于 $T_sH(s)$ 采样。

4.通过有理函数设计滤波器

由上文可知，对于模拟系统的传递函数为以部分分式展开的形式表示为
$H_c(s)=\sum_{k=1}^N\frac{A_k}{s-s_k}$
其中：

$s_k$ ：系统的极点
$A_k$ ：每个极点对应的系数
$N$ ：代表系统的极点个数。

对于它的反变换得到 $h (t)$ ：
$h(t)=\sum_{k=1}^NA_ke^{s_kt}u(t)$
由上文可知对其进行采样，需要扩大增益为 $T_s$ ，即：
$h[n]=T_sh(t)|_{t=nT_s} \quad =>\quad h[n]=T_s\sum_{k=1}^NA_ke^{s_knT_s}u[nT_s]$
在以前的文章中音频进阶学习十三——Z变换二（有理z变换、稳定性与反变换），我们知道 $a^nu[n]$ 的Z变换为 $\frac{1}{1-a^z{-1}}$ ，所以上式的Z变换为
$H(z)=\sum_{k=1}^N\frac{T_sA_k}{1-e^{s_kT_s}z^{-1}}$

5.总结

从模拟系统通过脉冲响应不变法设计数字系统的总过程分为：

通过模拟系统的有理传递函数得到脉冲响应
通过模拟系统的脉冲响应采样得到数字系统的脉冲响应
通过数字系统的脉冲响应得到数字系统的有理传递函数

最终表示为：
$H(z)=\sum_{k=1}^N\frac{T_sA_k}{1-e^{s_kT_s}z^{-1}}$
其中需要注意的是：

在低频段可以保持远小于采样频率，而高频段会出现混叠，所以脉冲响应不变法更适合于低通滤波器
由于使得数字系统脉冲响应很好地逼近原始系统的频率响应，所以时域逼近良好
相位线性变换，可以保持线性相位
为了满足数模转换幅频相同，所以要对于系统进行 $T_s$ 的增益

二、低通滤波器的设计实例

在上一篇文章音频进阶学习二十四——滤波器设计方法中，给定了数字滤波器设计的总过程如下图：
在这里插入图片描述
现在结合实例来一步步设计一个数字滤波器

1.给定数字滤波器指标

现在给出通带截止频率和阻带截止频率，设计一个低通的数字滤波器
$\begin{cases} 0.89125 \leq |H(e^{j\omega})| < 1, \quad |\omega| \leq 0.2\pi \\ |H(e^{j\omega})| < 0.17783, \quad 0.3\pi\leq|\omega| \leq \pi \end{cases}$
即如下图
在这里插入图片描述
也就是说，如果采样率是1k，那么根据归一化频率 $\omega=2\pi f/\Omega$ ：

通带截止频率： $\frac{2\pi f}{0.2\pi}=1000, \quad f=100Hz$
通带截止频率： $\frac{2\pi f}{0.3\pi}=1000, \quad f=150Hz$

2.转换模拟滤波器指标

根据脉冲响应不变法， $\Omega=\omega/T_s$ ，则模拟滤波器的指标为：
$\begin{cases} 0.89125 \leq |H(j\Omega)| < 1, \quad |\Omega| \leq \frac{0.2\pi}{T_s} \\ |H(j\Omega)| < 0.17783, \quad \frac{0.3\pi}{T_s} \leq|\Omega| \leq \frac{\pi}{T_s} \end{cases}$
如下图所示
在这里插入图片描述

3.模拟滤波器的设计

根据上一篇文章中给定的巴特沃斯滤波器
$|H_a(j\Omega)|^2=\frac{1}{1+(\frac{\Omega}{\Omega_c})^{2N}}$
将指标代入
$\begin{cases} |H(j\frac{0.2\pi}{T_s})| \geq 0.89125 \\ |H(j\frac{0.3\pi}{T_s})| \leq 0.17783\end{cases} => \begin{cases} 1+(\frac{0.2\pi}{T_s\Omega_c})^{2N} =(\frac{1}{0.89125 })^2 \\ \\ 1+(\frac{0.3\pi}{T_s\Omega_c})^{2N} =(\frac{1}{0.17783})^2 \end{cases}$
对于该式求解：
$\begin{cases} 1+(\frac{0.2\pi}{T_s\Omega_c})^{2N} =(\frac{1}{0.89125 })^2 \\ \\ 1+(\frac{0.3\pi}{T_s\Omega_c})^{2N} =(\frac{1}{0.17783})^2 \end{cases} => \begin{cases} (\frac{0.2\pi}{T_s\Omega_c})^{2N} =0.2589\\ \\ (\frac{0.3\pi}{T_s\Omega_c})^{2N} =30.622\end{cases} =>(\frac{0.2\pi}{0.3\pi})^{2N}=0.008456$
此时取对数 $(\frac{4}{9})^N=0.008456 , \quad N=5.8858$ ，即得到的指标为：
$\begin{cases} 1+(\frac{0.2\pi}{T_s\Omega_c})^{12} =(\frac{1}{0.89125 })^2 \\ \\ 1+(\frac{0.3\pi}{T_s\Omega_c})^{12} =(\frac{1}{0.17783})^2 \end{cases} =>\Omega_c=0.7032/T_s$

4.模拟滤波器的传递函数

在上文中说过，对于模拟系统的频率响应和传递函数关系为： $H_a(j\Omega)=H_c(s)|_{s=j\Omega}$ ,
由音频进阶学习十八——幅频响应相同系统、全通系统、最小相位系统中对于幅频响应相同系统的特点，对于巴特沃斯滤波器存在 $H(s)H^*(-s^*)=\frac{1}{1+(\frac{s}{j\Omega_c})^{2N}}$ ，由就可以得到对于给定指定模拟系统的极点
$(\frac{s}{j\Omega_c})^{12}=-1 \quad =>\quad \frac{s}{j\Omega_C}=e^{j\frac{2\pi k}{12}}=>s=j\Omega_ce^{j\frac{2\pi k}{12}}$
此时可知 $k = 0, 1, 2, ..., 11$ ，也就是说共有12个极点，且均匀地分布在一个以原点为中心，半径为 $\Omega_c$ 的圆上。所以每个极点的角度为 $2\pi/12=\pi/6$ .根据上文可知，S域的左半平面映射到Z平面保证在单位圆内，由此选取所左边平面的一组6个极点为：
$\begin{cases} [-0.182\pm j(0.679)]/T_s \\ [-0.497\pm j(0.497)]/T_s \\ [-0.679\pm j(0.182)]/T_s \end{cases}$
如同下图
在这里插入图片描述
将极点代入到模拟系统的传递函数为以部分分式展开为连加形式：
$H_c(s)=\sum_{k=1}^N\frac{A_k}{s-s_k}=\\ \frac{0.7032/T_s \times 0.7032/T_s}{(s-(-0.182 - j(0.679)/T_s)(-0.182 + j(0.679)/T_s)} \\ \times \frac{0.7032/T_s \times 0.7032/T_s}{(s-(-0.497- j(0.497)/T_s)(-0.497+ j(0.497)/T_s)} \\ \times \frac{0.7032/T_s \times 0.7032/T_s}{(s-(-0.679- j(0.182 )/T_s)(-0.679+ j(0.182 )/T_s)}$

5.转换为数字系统的传递函数

根据 $H(z)=\sum_{k=1}^N\frac{T_sA_k}{1-e^{s_kT_s}z^{-1}}$ ，将 $A_k,s_k$ 代入，得到数字滤波器的传递函数为
$H(z)=\sum_{k=1}^N\frac{A_k}{s-s_k}=\\ \frac{0.1435-0.2486i}{1-(0.6485- 0.5237i)z^{-1}}+ \frac{0.1435+0.2486i}{1-(0.6485+0.5237i)z^{-1}}\\ \frac{-1.0714-0.0001i}{1-(0.5346+0.2901i)z^{-1}}+ \frac{-1.0714+0.0001i}{1-(0.5346-0.2901i)z^{-1}} \\ \frac{0.9278-1.6093i}{1-(0.4986+0.0916i)z^{-1}}+ \frac{0.9278+1.6093i}{1-(0.4986-0.0916i)z^{-1}}$

三、使用matlab实现

1.通过传递函数实现

对于上文最终得到的传递函数，使用matlab画出频率响应，并给出信号为 50Hz 和 200Hz 叠加信号，进行滤波器过滤

% 合并共轭对为二阶分式（实数系数）

% 第一对：项1和项2
A1 = 0.1435 - 0.2486i;  % 分子系数
p1 = 0.6485 - 0.5237i;  % 极点

% 计算二阶分式的分子和分母
num_biquad1 = [2*real(A1), -2*real(A1*conj(p1))];  % 分子多项式
den_biquad1 = [1, -2*real(p1), abs(p1)^2];         % 分母多项式

% 第二对：项3和项4
A2 = -1.0714 - 0.0001i;
p2 = 0.5346 + 0.2901i;

num_biquad2 = [2*real(A2), -2*real(A2*conj(p2))];
den_biquad2 = [1, -2*real(p2), abs(p2)^2];

% 第三对：项5和项6
A3 = 0.9278 - 1.6093i;
p3 = 0.4986 + 0.0916i;

num_biquad3 = [2*real(A3), -2*real(A3*conj(p3))];
den_biquad3 = [1, -2*real(p3), abs(p3)^2];

% 创建传递函数并相加（实数系数）
sys_total = tf(num_biquad1, den_biquad1, 1, 'Variable', 'z^-1') + ...
            tf(num_biquad2, den_biquad2, 1, 'Variable', 'z^-1') + ...
            tf(num_biquad3, den_biquad3, 1, 'Variable', 'z^-1');

% 验证系数为实数
disp('分子系数：');
disp(sys_total.Numerator{1})
disp('分母系数：');
disp(sys_total.Denominator{1})

% 计算频率响应（采样率16kHz，实际设置Ts=1/16000）
Fs = 1000;
N_freq = 1024;
[H_total, w] = freqz(sys_total.Numerator{1}, sys_total.Denominator{1}, N_freq, Fs);

% 归一化频率（0到π）
f_normalized = w / (Fs/2) * pi; 

% 绘制归一化幅频响应
figure;
plot(f_normalized/pi, 20*log10(abs(H_total))); % 归一化为 ×π rad/sample
xlabel('Normalized Frequency (×π rad/sample)');
ylabel('Magnitude (dB)');
title('Normalized Frequency Magnitude Response');
grid on;


% 进行滤波测试
t = 0:1/Fs:1; % 1秒信号
x = sin(2*pi*50*t) + sin(2*pi*200*t); % 50Hz 和 200Hz 叠加信号
y = filter(b, a, x); % 通过滤波器

% 画出滤波前后的信号
figure;
subplot(2,1,1);
plot(t, x);
title('输入信号 (50Hz + 200Hz)');
xlabel('时间 (s)');
ylabel('幅值');
grid on;

subplot(2,1,2);
plot(t, y);
title('滤波后信号');
xlabel('时间 (s)');
ylabel('幅值');
grid on;

在这里插入图片描述

2.使用matlab工具

matlab自带的工具转换，只需要输入指标

clc; clear; close all;

% 设计巴特沃斯滤波器
N = 3; % 阶数 (可能需要调整)
Wp = 0.2 * pi;  % 归一化通带频率
Ws = 0.3 * pi;  % 归一化阻带频率
Rp = -20 * log10(0.89125); % 通带最大衰减 (dB)
Rs = -20 * log10(0.17783); % 阻带最小衰减 (dB)

% 设计模拟低通巴特沃斯滤波器
[~, Wn] = buttord(Wp/pi, Ws/pi, Rp, Rs);
[b, a] = butter(N, Wn);

Fs=1000;
% 计算数字滤波器频率响应
[H, w] = freqz(b, a, 1024, Fs);
f_normalized = w / (Fs/2) * pi; 

% 绘制幅频响应
figure;
subplot(2,1,1);
plot(f_normalized/pi, abs(H)); grid on;
xlabel('归一化频率 (\times\pi rad/sample)');
ylabel('幅值');
title('滤波器幅频响应');


% 进行滤波测试
t = 0:1/Fs:1; % 1秒信号
x = sin(2*pi*50*t) + sin(2*pi*200*t); % 50Hz 和 200Hz 叠加信号
y = filter(b, a, x); % 通过滤波器

% 画出滤波前后的信号
figure;
subplot(2,1,1);
plot(t, x);
title('输入信号 (50Hz + 200Hz)');
xlabel('时间 (s)');
ylabel('幅值');
grid on;

subplot(2,1,2);
plot(t, y);
title('滤波后信号');
xlabel('时间 (s)');
ylabel('幅值');
grid on;