计算机与信息工程学院验证性实验报告

专业：通信工程        年级/班级：20##级      20##—20##学年第一学期

一、实验目的

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

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

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

二、实验内容和要求

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

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

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

三、实验主要仪器设备和材料

计算机，MATLAB

四、实验方法

复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容，用双线性变换法设计数字滤波器系统函数。其中满足本实验要求的数字滤波器系统函数为：

                                              （4.1）

式中：

                            （4.2）

根据设计指标，调用MATLAB信号处理工具箱buttord和butter，也可以得到。

由公式（3.1）和（3.2）可见，滤波器由三个二阶滤波器、和级联而成，如图3-1所示。

编写滤波器仿真程序，计算对心电图采样序列x(n)的响应序列y(n)。

设y_k(n)为第k级二阶滤波器H_k(z)的输出序列，y_k-1(n)为输入序列，如图3-1所示。由（4.2）式可得到差分方程为：

（4.3）

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

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

本实验中用到的心电图信号采用序列x(n)

人体心电图信号在测量过程往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。在实验中，以x(n)作为输入序列，滤除其中的干扰成分。

五、实验程序及结果

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(1,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(1,2,2);

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];

k=1;

close all;

figure(1);

subplot(2,2,1);

n=0:55;     %更正

stem(n,x,'.');

axis([0,56,-100,50]);  %更正

hold on;

n=0:60;

m=zeros(61);

plot(n,m);

xlabel('n');

ylabel('x(n)');

title('心电图信号采集序列x(n)');

B=[0.09036,2*0.09036,0.09036];

A=[1.2686,-0.7051];

A1=[1.0106,-0.3583];

A2=[0.9044,-0.2155];

while(k<=3)

    y=filter(B,A,x);    %The function is to filte the singal x

    x=y;

    if k==2

        A=A1;

    end

    if k==3

        A=A2;

    end

    k=k+1;

end

subplot(2,2,3)

n=0:55;     %更正

stem(n,y,'.');

axis([0,56,-15,5]);

hold on;

n=0:60;

m=zeros(61);

plot(n,m);

xlabel('n');

ylabel('y(n)');

title('三级滤波后的心电图信号');

%求数字滤波器的幅频特性

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,w]=freqz(A,B1,100);

[H2,w]=freqz(A,B2,100);

[H3,w]=freqz(A,B3,100);

H4=H1.*(H2);

H=H4.*(H3);

mag=abs(H);

db=20*log10((mag+eps)/max(mag));

subplot(2,2,2)

plot(w/pi,db);

axis([0,0.5,-50,10]);

grid on;       %更正

title('滤波器的幅频响应曲线');

教师签名：

         年    月     日

 

第二篇：数字信号处理 实验 双线性变换法IIR数字滤波器设计

数字信号处理实验报告

实验名称：双线性变换法IIR数字滤波器设计

学号：               姓名：      

评语：   

成绩：  

实验目的

1、掌握用双线性变换法设计低通IIR数字滤波器的基本原理和算法。

2、掌握用双线性变换法设计高通和带通IIR数字滤波器的基本原理和算法。

3、进一步了解数字滤波器和模拟滤波器的频率响应特性。

实验原理与计算方法

1、双线性变换法设计IIR低通数字滤波器的基本原理和算法

双线性变换法设计数字滤波器，采用了二次映射的方法，就是先将整个s平面压缩到s₁平面的一个的横形条带范围内，然后再将这个条带映射到z平面上，就能建立s平面到z平面的一一对应关系。对于低通数字滤波器，映射关系为

　　            　                               (1)

其中T为抽样周期。

用双线性变换法设计低通IIR数字滤波器的基本步骤，首先根据设计要求确定相应的模拟滤波器的传递函数，再应用(1)式得数字滤波器的传递函数

　　　　　　　　　  　　　　　　　　　　　　　　      (2)

通常可以给定的参数为：低通数字滤波器通带边界频率、阻带边界频率和对应的通带衰减函数、阻带衰减函数。s₁平面中的模拟角频率与数字角频率的关系为线性关系，在计算模拟滤波器的阶数N、极点s_i和传递函数之前，应作预畸变处理

                      　　　 　　　　           (3)

　　模拟滤波器的阶数N、极点s_i和传递函数的计算方法与冲激响应不变法相同，可以采用Butterworth逼近或Chebyshev逼近。

2、双线性变换法设计IIR高通、带通、带阻数字滤波器的基本原理和算法

由于双线性变换法获得的数字滤波器频率响应特性中不会出现混叠现象，因此可以适用于高通、带通和带阻滤波器的设计。IIR数字滤波器的设计通常要借助于模拟低通滤波器的设计，由原型低通滤波器到其他形式（高通、带通、带阻）IIR数字滤波器的频带变换有模拟频带变换法和数字频带变换法。

（1）模拟频带变换法

首先将给定的对数字滤波器(DF)的技术要求转换为一个低通模拟滤波器(AF)的技术要求，根据这种要求用某种逼近设计出原型的低通模拟滤波器(LP AF)，计算出模拟滤波器的阶数N、极点s_i和传递函数，再按照双线性变换的变换关系，将模拟滤波器的传递函数转换为数字滤波器的传递函数。

表8-1中列出了将给定的对数字滤波器(DF)的技术要求直接转换为对一个低通模拟滤波器(AF)的技术要求的频率预畸变校正关系和转换公式。

表8-1  双线性变换和频率预校正的计算公式

例：数字高通滤波器的设计

首先将给定的数字高通滤波器的技术指标根据公式转换为模拟低通滤波器的技术指标，利用cheb1ord(Wp,Ws,ap,as,'s')函数求出chebyshev模拟低通滤波器的阶数N，再利用cheb1ap(N,ap)函数求出模拟低通滤波器系统函数的零极点，zp2tf(z,p,k)函数将零极点转换为系统函数系数；然后利用lp2hp由模拟低通滤波器的系统函数得到模拟带通滤波器的系统函数，bilinear函数则用于实现双线性变换法由模拟带通滤波器系统函数计算数字数字带通滤波器系统函数的系数。

（2）数字频带变换法

首先将给定的对数字滤波器(DF)的技术要求转换为一个低通模拟滤波器(AF)的技术要求，用双线性变换法将原型的低通模拟滤波器(LP AF)映射为低通数字滤波器，再将数字低通滤波器根据相应的变换公式经频带变换到各型数字滤波器。

例：数字高通滤波器的设计

函数[bhp,ahp]=zmapping(blp,alp,Nz,Dz)用来实现从数字低通滤波器得到数字高通滤波器的有理函数。

%数字滤波器技术指标

>>wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;

%对应的模拟滤波器技术指标

>>T=1;Fs=1/T;Wp=(2/T)*tan(wp/2); =(2/T)*tan(ws/2);

>>[cs,cd]=afd_chb1(Wp,Ws,Rp,As); %Chebyshev模拟滤波器

>> [blp,alp]=bilinear(cs,cd,Fs)     %双线性变换

>>wphp=0.6*pi;   %数字高通滤波器截止频率

%低通－高通频带变换

>>alpha=-(cos((wplp+wphp)/2))/(cos((wplp+wphp)/2))

>>Nz=-[alpha,1];Dz=[1,alpha];

>> [bhp,ahp]=zmapping(blp,alp,Nz,Dz) %数字高通滤波器的系统函数系数

（3）IIR数字滤波器的设计

可利用MATLAB提供的函数直接设计相应的数字滤波器。

函数buttord和cheb1ord用来根据给定的技术指标求出滤波器的阶数N和边界频率wn，butter和cheby1则根据阶数和边界频率设计相应的数字滤波器。输入的参数不同则所设计的滤波器类型不同。

[N,wn]=buttord(wp,ws,Rp,As);

[N,wn]=cheb1ord(wp,ws,Rp,As);

[b,a]=butter(N,wn);

[b,a]=cheby1(N,Rp,wp);

实验内容及结果

（1）Chebyshev IIR数字带通滤波器满足如下技术指标

低阻带边界频率，高阻带边界频率，阻带衰减函数

低通带边界频率，高通带边界频率，通带波动

抽样频率,记录所得的模拟滤波器的阶数N，画出模拟滤波器和数字滤波器的频率响应的幅频和相频特性曲线。N =2

实验代码：

fs1=100;

fs2=600;

fp1=200;

fp2=400;

fsa=2000;

As=18;

Rp=2;

T=1./fsa;                      %对应的模拟滤波器技数值

w1=2.*pi.*(fp1./fsa);        %Chebyshev模拟滤波器

w2=2.*pi.*(fp2./fsa);

wp1=2*pi*fp1*T;

wp2=2*pi*fp2*T;

ws2=2.*pi.*(fs2./fsa);

cosw0=(sin(w1+w2))./(sin(w2)+sin(w1));

w0=acos(cosw0);

bw=wp2-wp1;

Wp=(cosw0-cos(w2))./sin(w2);

Ws=(cosw0-cos(ws2))./sin(ws2);

[N,omgn]=cheb1ord(Wp,Ws,Rp,As,'s'); %返回模拟低通滤波器阶数N和边界频率n

[z,p,k]=cheb1ap(N,Rp);                 %得系统函数零极点

[blp,alp]=zp2tf(z,p,k);                %由零极点得系数

[bhp,ahp]=lp2bp(blp,alp,w0,bw);      %模拟低通到模拟带通

[bdf,adf]=bilinear(bhp,ahp,1);       %双线性变换将模拟带通滤波器转换成数字带通滤波器

[BPA,wa]=freqs(bhp,alp,fsa);

[BPD,wd]=freqz(bdf,adf,fsa);

subplot(2,2,1);

plot(abs(BPA));

title('模拟带阻滤波器幅频特性');

subplot(2,2,2);

plot(angle(BPA));

title('模拟带阻滤波器相频特性');

subplot(2,2,3);

plot(abs(BPD));

title('数字带阻滤波器幅频特性');

subplot(2,2,4);

plot(angle(BPD));

title('数字带阻滤波器相频特性')

实验截图：

（2）Chebyshev IIR数字带阻滤波器满足如下技术指标

低阻带边界频率，高阻带边界频率，阻带衰减函数

低通带边界频率，高通带边界频率，通带波动

抽样频率,记录所得的模拟滤波器的阶数N，画出模拟滤波器和数字滤波器的频率响应的幅频和相频特性曲线。N =2

实验代码：

fs1=1000;

fs2=2000;

fp1=500;

fp2=3000;

fsa=10000;

As=18;

Rp=2;

T=1./fsa;             %对应的模拟滤波器计数值

w1=2.*pi.*(fp1./fsa);                %Chebyshev模拟滤波器

w2=2.*pi.*(fp2./fsa);

wp1=2*pi*fp1*T;

wp2=2*pi*fp2*T;

ws2=2.*pi.*(fs2./fsa);

cosw0=(sin(w1+w2))./(sin(w2)+sin(w1));

w0=acos(cosw0);

bw=wp2-wp1;

Wp=(cosw0-cos(w2))./sin(w2);

Ws=(cosw0-cos(ws2))./sin(ws2);

[N,omgn]=cheb1ord(Wp,Ws,Rp,As,'s');       %返回模拟低通滤波器阶数N和边界频率n [z,p,k]=cheb1ap(N,Rp);                       %得系统函数零极点

[blp,alp]=zp2tf(z,p,k);                      %由零极点得系数

[bhp,ahp]=lp2bs(blp,alp,w0,bw);             %模拟低通到模拟带通

[bdf,adf]=bilinear(bhp,ahp,1);       %双线性变换将模拟带通滤波器转换成数字带通滤波器

[BPA,wa]=freqs(bhp,alp,fsa);

[BPD,wd]=freqz(bdf,adf,fsa);

subplot(2,2,1);

plot(abs(BPA));

title('模拟带阻滤波器幅频特性');

subplot(2,2,2);

plot(angle(BPA));

title('模拟带阻滤波器相频特性');

subplot(2,2,3);

plot(abs([BPD]));

title('数字带阻滤波器幅频特性');

subplot(2,2,4);

plot(angle([BPD]));

title('数字带阻滤波器相频特性')

实验截图：