数字信号处理用脉冲响应不变法和双线性变换法设计巴特沃斯滤波器MATLAB实现——实例

符号含义

例题

直接计算步骤得到H(s)

脉冲响应不变法

设计思路

1. 借助模拟滤波器设计方法设计滤波器的系统函数H(s)
2. 连续系统的系统函数H(s)用脉冲响应不变法转换成离散系统的系统函数H(z)

下面是巴特沃斯滤波器的幅频特性

clear
close all
clc

fs=1000;%采样频率
fc=200;%通带截止频率
fr=300;%阻带截止频率
T=0.001; %采样周期


%%%%%%%脉冲响应不变法
wp1=2*pi*fc;%通带截止频率
wr1=2*pi*fr;%阻带截止频率

[N1,wn1]=buttord(wp1,wr1,1,25,'s'); %[滤波器的阶数,截止频率]=buttord(通带截止频率,阻带截止频率,通带波纹,阻带衰减,滤波器类型s-带阻)
[B1,A1]=butter(N1,wn1,'s'); %[传递函数系数B1是分子的 A1是分母的]=butter(滤波器的阶数,截止频率,滤波器类型s-带阻)
[num1,den1]=impinvar(B1,A1,fs); %利用脉冲响应不变法将模拟滤波器转数字滤波器[数字滤波器传递函数num1/den1]=impinvar(模拟的分子系数B1,模拟的分母系数A1,采样频率)
[h1,w]=freqz(num1,den1); %[频率响应,角频率]=freqz(数字滤波器传递函数num1/den1)

函数解读：

buttord [n，Wn] = buttord（Wp，Ws，Rp，Rs）返回数字巴特沃思滤波器的最低阶n，通带纹波不超过Rp dB，阻带衰减至少Rs dB。
Wp和Ws分别是滤波器的通带和阻带边沿频率，归一化为0到1，其中1对应于πrad / sample。
还返回相应的截止频率Wn的标量（或矢量）。
要设计巴特沃斯滤波器，请使用输出参数n和Wn作为butter的输入

butter butter（n，Wn）返回归一化截止频率Wn的n阶低通数字巴特沃斯滤波器的传递函数系数。

impinvar impinvar (B1,A1,fs)创建了一个分子系数和分母系数分别为B1和A1的数字滤波器，其脉冲响应等于系数为B1和A1的模拟滤波器的脉冲响应，按1/fs缩放，其中fs为采样率。

freqz [h,w] = freqz(d,n) 返回数字滤波器 d 的 n 点复频率响应。

双线性变换法

双线性变换法要对通带截止频率和阻带截止频率做一个预畸变处理

$\Omega =\frac{2}{T}\tan \frac{{\Omega }_{1}T}{2}=\frac{2}{T}\tan \frac{2\pi fT}{2}$

%%%%%%%双线性变换法
wp2=2*fs*tan(2*pi*fc/(2*fs)); %做预畸变处理后的通带截止频率
wr2=2*fs*tan(2*pi*fr/(2*fs)); %做预畸变处理后的阻带截止频率
[N2,wn2]=buttord(wp2,wr2,1,25,'s'); %[滤波器的阶数,截止频率]=buttord(通带截止频率,阻带截止频率,通带波纹,阻带衰减,滤波器类型s-带阻)
[B2,A2]=butter(N2,wn2,'s'); %[传递函数系数B2是分子的 A2是分母的]=butter(滤波器的阶数,截止频率,滤波器类型s-带阻)
[num2,den2]=bilinear(B2,A2,fs); %利用双线性的方法将模拟滤波器转数字滤波器[数字滤波器传递函数num2/den2]=impinvar(模拟的分子系数B2,模拟的分母系数A2,采样频率)
[h2,w]=freqz(num2,den2); %[频率响应,角频率]=freqz(数字滤波器传递函数num2/den2)
f=w/(2*pi)*fs; %将角频率转换为频率

画图

plot(f,20*log10(abs(h1)),'-.',f,20*log10(abs(h2)),'-'); %画脉冲响应不变法和双线性变换法的图
axis([0,500,-80,10]);% 设置当前坐标轴x轴和y轴的限制范围
grid; %切换改变主网格线的可见性
xlabel('频率/Hz');ylabel('幅度/dB');title('巴特沃思数字低通滤波器'); 
legend('脉冲响应不变法','双线性变换法');

完整代码

clear
close all
clc

fs=1000;%采样频率
fc=200;%通带截止频率
fr=300;%阻带截止频率
T=0.001; %采样周期


%%%%%%%脉冲响应不变法
wp1=2*pi*fc;%通带截止频率
wr1=2*pi*fr;%阻带截止频率

[N1,wn1]=buttord(wp1,wr1,1,25,'s'); %[滤波器的阶数,截止频率]=buttord(通带截止频率,阻带截止频率,通带波纹,阻带衰减,滤波器类型s-带阻)
[B1,A1]=butter(N1,wn1,'s'); %[传递函数系数B1是分子的 A1是分母的]=butter(滤波器的阶数,截止频率,滤波器类型s-带阻)
[num1,den1]=impinvar(B1,A1,fs); %利用脉冲响应不变法将模拟滤波器转数字滤波器[数字滤波器传递函数num1/den1]=impinvar(模拟的分子系数B1,模拟的分母系数A1,采样频率)
[h1,w]=freqz(num1,den1); %[频率响应,角频率]=freqz(数字滤波器传递函数num1/den1)

%%%%%%%双线性变换法
wp2=2*fs*tan(2*pi*fc/(2*fs)); %做预畸变处理后的通带截止频率
wr2=2*fs*tan(2*pi*fr/(2*fs)); %做预畸变处理后的阻带截止频率
[N2,wn2]=buttord(wp2,wr2,1,25,'s'); %[滤波器的阶数,截止频率]=buttord(通带截止频率,阻带截止频率,通带波纹,阻带衰减,滤波器类型s-带阻)
[B2,A2]=butter(N2,wn2,'s'); %[传递函数系数B2是分子的 A2是分母的]=butter(滤波器的阶数,截止频率,滤波器类型s-带阻)
[num2,den2]=bilinear(B2,A2,fs); %利用双线性的方法将模拟滤波器转数字滤波器[数字滤波器传递函数num2/den2]=impinvar(模拟的分子系数B2,模拟的分母系数A2,采样频率)
[h2,w]=freqz(num2,den2); %[频率响应,角频率]=freqz(数字滤波器传递函数num2/den2)
f=w/(2*pi)*fs; %将角频率转换为频率

plot(f,20*log10(abs(h1)),'-.',f,20*log10(abs(h2)),'-'); %画脉冲响应不变法和双线性变换法的图
axis([0,500,-80,10]);% 设置当前坐标轴x轴和y轴的限制范围
grid; %切换改变主网格线的可见性
xlabel('频率/Hz');ylabel('幅度/dB');title('巴特沃思数字低通滤波器'); 
legend('脉冲响应不变法','双线性变换法');