实验三:用双线性变换法设计 IIR 数字滤波器
(设计性 4 学时)
一.实验目的:
(1)熟悉用双线性变换法设计 IIR 数字滤波器的原理与方法。
(2)掌握数字滤波器的计算机仿真方法。
(3)通过观察对实际心电图信号的滤波作用,获得数字滤波的感性知识。
二.实验内容及步骤:
(1) 用双线性变换法设计一个巴特沃斯低通 IIR 数字滤波器,设计指标参数为:在通带内频
率低于 0.2pi 时,最大衰减小于 1dB;在阻带内[0.3pi , pi] 频率区间上,,最小衰减大于 15dB;
(2) 以 0.02pi 为采样间隔,打印出数字滤波器在频率区间[ 0, 0.5pi]上的幅频响应特性曲线;
(3) 用所设计的滤波器对实际心电图信号采样序列(在本实验后面给出)进行仿真滤波处理,
并分别打印出滤波前后的心电图波形图,观察总结滤波作用与效果。
(4)采用不同阶数的 Butterworth 低通滤波器,比较滤波效果。
三.实验步骤:
(1)复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计 IIR 数字滤波器的内容,按
照教材例 6.4.2,用双线性变换法设计数字滤波器系统函数 H(z)。
方法一:教材例 6.4.2 种已求出满足本实验要求的数字滤波器系统函数:
方法二:根据设计指标,调用 MATLAB 信号处理工具箱函数 buttord 和 butter,也可得到 H
(z)。
(2)编写滤波器仿真程序,计算 H(z)对心电图信号采样序列 x(n)的相应序列 y(n)。
(3)在通过计算机上运行仿真滤波程序,并调用通用绘图子程序,完成实验内容(2)和(3)。
本实验要用的 MATLAB 绘图函数参阅教材。
四.,思考题:
用双线性变换法设计数字滤波器过程中,变换公式:
s=
2
T
z1
z1
1
1
中 T 的取值,对设计结果有无影响? 为什么?
五.实验报告要求
(1)简述实验目的及原理;
(2)由所打印的特性曲线及设计过程简述双线性变换法的特点;
(3)对比滤波前后的心电图信号波形,说明数字滤波器的滤波过程与滤波作用;
(4) 简要回答思考题.
六:心电图信号采样序列 x(n):
人体心电图信号在测量过程中往往受到工业高频干扰,所以必须经过低通滤波处理后,
才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样式本 x(n),其
中存在高频干扰,在实验中,以 x(n)作为输入序列,滤除其中的干扰成分。
{ x(n) } = { -4 , -2, 0, -4, -6, -4, -2, -4, -6, -6, -4, -4, -6, -6, -2, 6, 12, 8, 0, -16
-38, -60, -84, -90, -66, -32, -4, -2, -4, 8,12, 12 , 10, 6, 6, 6, 4, 0, 0, 0
0, 0, -2, -4, 0, 0, 0, -2, -2, 0,0 , -2,, -2, -2, -2, 0}
参考程序:
一;
T=1;Fs=1/T;
wpz=0.2;wsz=0.3;
wp=2*tan(wpz*pi/2);ws=2*tan(wsz*pi/2);rp=1;rs=15;
[N,wc]=buttord(wp,ws,rp,rs,'s');
[B,A]=butter(N,wc,'s');
fk=0:1/512:1;wk=2*pi*fk;
Hk=freqs(B,A,wk);
subplot(2,2,1);
plot(fk,20*log10(abs(Hk)));grid on;
xlabel('\omega/\pi');ylabel('幅度(dB)');
axis([0,1,-100,5]);title('(b)');
[N,wc]=buttord(wpz,wsz,rp,rs);
[Bz,Az]=butter(N,wc);
wk=0:pi/512:pi;
Hz=freqz(Bz,Az,wk);
subplot(2,2,4);
plot(wk/pi,20*log10(abs(Hz)));grid on;
xlabel('\omega/\pi');ylabel('幅度(dB)');
axis([0,1,-100,5]);title('(b)');
二;
x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,
6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];
subplot(2,2,1);
n=0:55;
stem(n,x,'.');
xlabel('n');
ylabel('x(n)');
title('x(n)的脉冲响应');
A=0.09036;
b1=[A,2*A,A];
a1=[1,-1.2686,0.7051];
h1=filter(b1,a1,x);
[H1,w]=freqz(b1,a1,100);
b2=[A,2*A,A];
a2=[1,-1.0106,0.3583];
h2=filter(b2,a2,h1);
[H2,w]=freqz(b2,a2,100);
b3=[A,2*A,A];
a3=[1,-0.9044,0.2155];
h3=filter(b3,a3,h2);
[H3,w]=freqz(b3,a3,100);
subplot(2,2,2);
stem(n,h3,'.');
xlabel('n');
ylabel('y(n)');
title('通过滤波器 H1(z),H2(z),H3(z)后的 y3(n)函数');
subplot(2,2,3);
H4=H1.*(H2);
H=H4.*(H3);
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
plot(w/pi,db);
xlabel('ω/π');
ylabel('20log[Ha3(ejw)]');
title('通过滤波器 H1(z),H2(z),H3(z)后的对数频率响应 20log[Ha3(ejw)]函数');
grid;
figure(2);
N=1024;n=0:N/2-1;
Xk=fft(x,N);AXk=abs(Xk(1:N/2));
f=(0:N/2-1)*Fs/N;
f=f/Fs;
subplot(211);plot(f,AXk);title('x(n)的频谱');
xlabel('f');ylabel('| X(k) |');axis([0,0.5,0,400]);
Yk=fft(y,N);AYk=abs(Yk(1:N/2));
subplot(212);plot(f,AYk);title('y(n)的频谱');
xlabel('f');ylabel('| Y(k) |');axis([0,0.5,0,400])
x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,
6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];
subplot(2,2,1);
n=0:55;
stem(n,x,'.');
xlabel('n');
ylabel('x(n)');
title('x(n)的脉冲响应');
A=0.09036;
b1=[A,2*A,A];
a1=[1,-1.2686,0.7051];
h1=filter(b1,a1,x);
[H1,w]=freqz(b1,a1,100);
b2=[A,2*A,A];
a2=[1,-1.0106,0.3583];
h2=filter(b2,a2,h1);
[H2,w]=freqz(b2,a2,100);
b3=[A,2*A,A];
a3=[1,-0.9044,0.2155];
h3=filter(b3,a3,h2);
[H3,w]=freqz(b3,a3,100);
subplot(2,2,2);
stem(n,h3,'.');
xlabel('n');
ylabel('y(n)');
title('通过滤波器 H1(z),H2(z),H3(z)后的 y3(n)函数');
subplot(2,2,3);
H4=H1.*(H2);
H=H4.*(H3);
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
plot(w/pi,db);
xlabel('ω/π');
ylabel('20log[Ha3(ejw)]');
title('通过滤波器 H1(z),H2(z),H3(z)后的对数频率响应 20log[Ha3(ejw)]函数');
grid;
%x(n)的心电脉冲函数
x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,
-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0]
subplot(2,2,1);
n=0:55;
stem(n,x,'.');
xlabel('n');
ylabel('x(n)');
title('x(n)的心电脉冲函数');
%通过滤波器 H1(z)后的 y1(n)函数
A=0.09036;
b1=[A,2*A,A];
a1=[1,-1.2686,0.7051];
h1=filter(b1,a1,x);
[H1,w]=freqz(b1,a1,100);
%通过滤波器 H1(z),H2(Z)后的 y2(n)函数
b2=[A,2*A,A];
a2=[1,-1.0106,0.3583];
h2=filter(b2,a2,h1);
[H2,w]=freqz(b2,a2,100);
%通过滤波器 H1(z),H2(Z),H3(Z)后的 y3(n)函数
b3=[A,2*A,A];
a3=[1,-0.9044,0.2155];
h3=filter(b3,a3,h2);
[H3,w]=freqz(b3,a3,100);
subplot(2,2,2);
stem(n,h3,'.');
xlabel('n');
ylabel('y(n)');
title('通过滤波器 H1(Z),H2(Z),H3(Z)后的 y3(n)函数');
subplot(2,2,3);
H4=H1.*(H2);
H=H4.*(H3);
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
plot(w/pi,db);
xlabel('w/pi');
ylabel('20log[Ha3(ejw)]');
title('通过滤波器 H1(z)、H2(z)、H3(z)后的对数频率响应 20log[Ha3(ejw)]函数');
grid;