用双线性变换法设计IIR数字滤波器

一、实验目的

1、熟悉用双线性变换法设计IIR数字滤波器的原理与方法。

2、掌握数字滤波器的计算机仿真方法。

3、通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。

二、实验内容

1、 用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为：在通带内频率低于0.2Π时，最大衰减小于1dB；在阻带内[0.3Π，Π]频率区间上，最小衰减大于15dB。

2、 以T=1s为采样间隔，打印出数字滤波器在频率区间[0，Π/2]上的幅频响应特性曲线。

3、 用所设计的滤波器对实际心电图信号采样序列(在本实验后面给出)进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。

三、实验步骤

1、 复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容，用双线性变换法设计数字滤波器系统函数H(z）。

  教材中满足本实验要求的数字滤波器函数

  式中：

，k=1，2，3

A=0.09036，

B1=1.2686，  C1=-0.7051

B2=1.0106，  C2=-0.3583

B3=0.9044，  C3=-0.2155

根据设计指标，调用MATLAB信号处理工具函数buttord和butter，也可得到H(z）。

由滤波器的函数可见，滤波器H(z）

由三个二阶滤波器H1(z）、H2(z）、H3(z）级联组成。

2、 编写滤波器仿真程序，计算H(z）对心电图信号采样序列 的响应序列y(n)。

yk(n)为第K阶滤波器Hk(z)的输出序列，Yk-1(n)为输入输入序列，根据滤波器H(z)的组成可得差分方程。

当K=1时所以H(z）对x(n)的总响序列y(n)可以用顺序迭代算法得到。即依次对k=1，2，3求解差分方程，最后得到y3(n)=y(n)仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序，也可直接调用MATLAB filter函数实现仿真。

3、 在通用计算机上运行仿真滤波程序，完成实验内容（2）和（3）。

人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本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}

代码如下所示：

%%
clear all；
close all;
T=1;    %可任取，一般取1
Fs=1/T;
Wp=0.2*pi;
Ws=0.3*pi;
Rp=1;
As=15;
%给出的是数字滤波器指标，先转换成原型模拟滤波器指标
Wap=2/T*tan(Wp/2);
Was=2/T*tan(Ws/2);
[N,Wn]=buttord(Wap,Was,Rp,As,'s'); %阶次估计函数(符合指标的阶数和截止频率)
[z0,p0,k0]=buttap(N);                        %归一化巴特沃斯低通模拟原型函数(z0,p0,k0分别为零点、极点、增益)
[Bap,Aap]=zp2tf(z0,p0,k0);                   %将系统函数的零极点转化为系统函数一般形式的系数，将零点极点转化为模拟低通滤波器的传递函数表示
[bs,as]=lp2lp(Bap,Aap,Wn);               %频率转换函数，将模拟归一化原型低通滤波器转化为要求设计的去归一化的模拟低通滤波器
[bz,az]=bilinear(bs,as,Fs);                  %双线性变换法进行模拟数字滤波器变换
%打印出H（z）的分子、分母多项式系数
bz;
az;
%以0.02pi 为采样间隔，在频率区间[0,pi/ 2]上的幅频响应特性曲线
[H,W]=freqz(bz,az);   %计算频率响应freqz数字滤波器频率响应
figure(1); %绘制幅频响应曲线
plot(W,abs(H));
figure(2);
title('对数幅频特性曲线图');
plot(W,20*log10(abs(H)));

%信号滤波
xn=[-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]; 
figure(3);
subplot(1,1,1);
stem(xn);
y=filter(bz,az,xn);%进行滤波，x为对哪个信号进行滤波 y为滤波输出结果
figure(4);
subplot(1,1,1);
stem(y);%对y产生的滤波

%%
clear all；
close all;
k=1;
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,0];
figure(1)
subplot(2,1,1)
n=0:56;
stem(n,x,'o');%绘制离散序列数据
axis([0 56 -100 80]);
hold on;
n=0:60;
m=zeros(61);%创建一个由零值组成的 61×61矩阵。
ylabel('x(n)');
title('心电图信号采样序列x(n)');
B=[0.09036 2*0.09036 0.09036];
A=[1.2686 -0.7051];
plot(n,m);
xlabel('n');
A1=[1.0106 -0.3583];
A2=[0.9044 -0.2155];
while(k<=3)
 y=filter(B,A,x);
 x=y;
 if k==2
     A=A1;
 end
 if k==3
A=A2;
 end
 k=k+1;
end
subplot(2,1,2)
n=0:56;
stem(n,y,'*');
axis([0 56 -15 15]);
hold on;
n=0:60;
m=zeros(61);
plot(n,m);
xlabel('n');
ylabel('y(n)');
title('双线性变换法设计IIR滤波后的心电图信号');
%双线性变换法设计ButterWorth数字滤波器幅频特性
A=[0.09036,0.1872,0.09036];
B1=[1,-1.2686,0.7051];
B2=[1,-1.0106,0.3583];
B3=[1,-0.9044,0.2155];
[H1,~]=freqz(A,B1,100);%数字滤波器频率响应
[H2,~]=freqz(A,B2,100);
[H3,w]=freqz(A,B3,100);
H4=H1.*(H2);
H=H4.*(H3);
mag=abs(H);%计算H的绝对值和复数
db=20*log10((mag+eps)/max(mag));
figure(2);
subplot(2,1,1)%将当前图窗划分为 2×1 网格，并在 第一行 指定的位置画图。
plot(w/pi,db);
axis([0,1.2,-160,10]);%更改坐标轴范围，使X轴的范围为0到1.2，y轴的范围从-160到10.
xlabel('w/π');%x轴坐标为w/Π
ylabel('20lg|H(ejw|');%Y轴坐标轴为20lg|H(ejw|
title('双线性变换法设计IIR滤波器的幅频响应曲线');%标题
%脉冲响应不变法设计ButterWorth数字滤波器幅频特性
C1=[0.2871 -0.4466];
C2=[-2.1428 1.1454];
C3=[1.8558 -0.6304];
D1=[1 -0.1297 0.6949];
D2=[1 -1.0691 0.3699];
D3=[1 -0.9972 0.2570];
[E1,~]=freqz(C1,D1,100); %计算频率响应freqz数字滤波器频率响应
[E2,~]=freqz(C2,D2,100); %计算频率响应freqz数字滤波器频率响应
[E3,w]=freqz(C3,D3,100); %计算频率响应freqz数字滤波器频率响应
E4=E1+E2;
E=E4+E3;
mag=abs(E);%计算E的绝对值
db=20*log10((mag+eps)/max(mag));
subplot(2,1,2)%将当前图窗划分为 2×1 网格，并在 第二行 指定的位置画图。
plot(w/pi,db);%做x为w/pi，y为db的二维线图。
axis([0,1,-50,10]);%更改坐标轴范围，使X轴的范围为0到1，y轴的范围从-50到10.
xlabel('w/π');
ylabel('20lg|H(ejw|');