实验四IIR数字滤波器设计及软件实现实验报告

一、实验目的 

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

（2）学会调用MATLAB信号处理工具箱中滤波器设计函数（或滤波器设计分析工具fdatool）设计各种IIR数字滤波器，学会根据滤波需求确定滤波器指标参数。

（3）掌握IIR数字滤波器的MATLAB实现方法。

（3）通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念。

二、实验原理

设计IIR数字滤波器一般采用间接法（脉冲响应不变法和双线性变换法），应用最广泛的是双线性变换法。基本设计过程是：①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标； ②设计过渡模拟滤波器；③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。第六章介绍的滤波器设计函数butter、cheby1 、cheby2 和ellip可以分别被调用来直接设计巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。本实验要求读者调用如上函数直接设计IIR数字滤波器。

本实验的数字滤波器的MATLAB实现是指调用MATLAB信号处理工具箱函数filter对给定的输入信号x(n)进行滤波，得到滤波后的输出信号y(n）。

三、实验内容及步骤

（1）调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，该函数还会自动绘图显示st的时域波形和幅频特性曲线，如图10.4.1所示。由图可见，三路信号时域混叠无法在时域分离。但频域是分离的，所以可以通过滤波的方法在频域分离，这就是本实验的目的。

图10.4.1  三路调幅信号st的时域波形和幅频特性曲线

（2）要求将st中三路调幅信号分离，通过观察st的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率。要求滤波器的通带最大衰减为0.1dB,阻带最小衰减为60dB。

提示：抑制载波单频调幅信号的数学表示式为

其中，称为载波，f_c为载波频率，称为单频调制信号，f₀为调制正弦波信号频率，且满足。由上式可见，所谓抑制载波单频调幅信号，就是2个正弦信号相乘，它有2个频率成分：和频和差频，这2个频率成分关于载波频率f_c对称。所以，1路抑制载波单频调幅信号的频谱图是关于载波频率f_c对称的2根谱线，其中没有载频成分，故取名为抑制载波单频调幅信号。容易看出，图10.4.1中三路调幅信号的载波频率分别为250Hz、500Hz、1000Hz。如果调制信号m(t)具有带限连续频谱，无直流成分，则就是一般的抑制载波调幅信号。其频谱图是关于载波频率f_c对称的2个边带（上下边带），在专业课通信原理中称为双边带抑制载波 (DSB-SC) 调幅信号,简称双边带 (DSB) 信号。如果调制信号m(t)有直流成分，则就是一般的双边带调幅信号。其频谱图是关于载波频率f_c对称的2个边带（上下边带），并包含载频成分。

（3）编程序调用MATLAB滤波器设计函数ellipord和ellip分别设计这三个椭圆滤波器，并绘图显示其幅频响应特性曲线。

（4）调用滤波器实现函数filter，用三个滤波器分别对信号产生函数mstg产生的信号st进行滤波，分离出st中的三路不同载波频率的调幅信号y₁(n)、y₂(n)和y₃(n)， 并绘图显示y1(n)、y2(n)和y3(n)的时域波形，观察分离效果。

四、信号产生函数mstg清单

function st=mstg

%产生信号序列向量st,并显示st的时域波形和频谱

%st=mstg 返回三路调幅信号相加形成的混合信号，长度N=1600

N=1600   %N为信号st的长度。

Fs=10000;T=1/Fs;Tp=N*T; %采样频率Fs=10kHz，Tp为采样时间

t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;

fc1=Fs/10;      %第1路调幅信号的载波频率fc1=1000Hz,

fm1=fc1/10;     %第1路调幅信号的调制信号频率fm1=100Hz

fc2=Fs/20;     %第2路调幅信号的载波频率fc2=500Hz

fm2=fc2/10;     %第2路调幅信号的调制信号频率fm2=50Hz

fc3=Fs/40;      %第3路调幅信号的载波频率fc3=250Hz,

fm3=fc3/10;     %第3路调幅信号的调制信号频率fm3=25Hz

xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %产生第1路调幅信号

xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %产生第2路调幅信号

xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %产生第3路调幅信号

st=xt1+xt2+xt3;         %三路调幅信号相加

fxt=fft(st,N);          %计算信号st的频谱

%====以下为绘图部分，绘制st的时域波形和幅频特性曲线=

subplot(3,1,1)

plot(t,st);grid;xlabel('t/s');ylabel('s(t)');

axis([0,Tp/8,min(st),max(st)]);title('(a) s(t)的波形')

subplot(3,1,2)

stem(f,abs(fxt)/max(abs(fxt)),'.');grid;title('(b) s(t)的频谱')

axis([0,Fs/5,0,1.2]);

xlabel('f/Hz');ylabel('幅度')

五、实验程序框图如图10.4.2所示

六、滤波器参数及实验程序清单

1、滤波器参数选取

观察图10.4.1可知，三路调幅信号的载波频率分别为250Hz、500Hz、1000Hz。带宽（也可以由信号产生函数mstg清单看出）分别为50Hz、100Hz、200Hz。所以，分离混合信号st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的指标参数选取如下：

对载波频率为250Hz的条幅信号，可以用低通滤波器分离，其指标为

带截止频率Hz，通带最大衰减dB；

阻带截止频率Hz，阻带最小衰减dB，

对载波频率为500Hz的条幅信号，可以用带通滤波器分离，其指标为

带截止频率Hz，Hz，通带最大衰减dB；

阻带截止频率Hz，Hz，Hz，阻带最小衰减dB，

对载波频率为1000Hz的条幅信号，可以用高通滤波器分离，其指标为

带截止频率Hz，通带最大衰减dB；

阻带截止频率Hz，阻带最小衰减dB，

说明：（1）为了使滤波器阶数尽可能低，每个滤波器的边界频率选择原则是尽量使滤波器过渡带宽尽可能宽。

（2）与信号产生函数mstg相同，采样频率Fs=10kHz。

（3）为了滤波器阶数最低，选用椭圆滤波器。

按照图10.4.2 所示的程序框图编写的实验程序为exp4.m。

2、实验程序清单

%实验4程序

% IIR数字滤波器设计及软件实现

clear all;close all

Fs=10000;T=1/Fs;   %采样频率

%调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st

st=mstg;

%低通滤波器设计与实现

fp=280;fs=450;

wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;   %DF指标（低通滤波器的通、阻带边界频）

[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip(N,rp,rs,wp);      %调用ellip计算椭圆带通DF系统函数系数向量B和A

y1t=filter(B,A,st);     %滤波器软件实现

% 低通滤波器设计与实现绘图部分

figure(2);subplot(3,1,1);

myplot(B,A);  %调用绘图函数myplot绘制损耗函数曲线

yt='y_1(t)';

subplot(3,1,2);tplot(y1t,T,yt); %调用绘图函数tplot绘制滤波器输出波形

%带通滤波器设计与实现

fpl=440;fpu=560;fsl=275;fsu=900;

wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60;

[N,wp]=ellipord(wp,ws,rp,rs);    %调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip(N,rp,rs,wp); %调用ellip计算椭圆带通DF系统函数系数向量B和A

y2t=filter(B,A,st);     %滤波器软件实现

% 带通滤波器设计与实现绘图部分（省略）

%高通滤波器设计与实现

fp=890;fs=600;

wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;   %DF指标（低通滤波器的通、阻带边界频）

[N,wp]=ellipord(wp,ws,rp,rs);    %调用ellipord计算椭圆DF阶数N和通带截止频率wp

[B,A]=ellip(N,rp,rs,wp,'high'); %调用ellip计算椭圆带通DF系统函数系数向量B和A

y3t=filter(B,A,st);     %滤波器软件实现

% 高低通滤波器设计与实现绘图部分（省略）

七、实验程序运行结果

实验4程序exp4.m运行结果如图104.2所示。由图可见，三个分离滤波器指标参数选取正确，算耗函数曲线达到所给指标。分离出的三路信号y1(n)，y2(n)和y3(n)的波形是抑制载波的单频调幅波。

(a) 低通滤波器损耗函数及其分离出的调幅信号y₁(t)

(b) 带通滤波器损耗函数及其分离出的调幅信号y₂(t)

(c)高通滤波器损耗函数及其分离出的调幅信号y₃(t)

 图104. 实验4程序exp4.m运行结果

八、思考题简答

（1）请阅读信号产生函数mstg，确定三路调幅信号的载波频率和调制信号频率。

（2）信号产生函数mstg中采样点数N=800，对st进行N点FFT可以得到6根理想谱线。如果取N=1000，可否得到6根理想谱线？为什么？N=2000呢？请改变函数mstg中采样点数N的值，观察频谱图验证您的判断是否正确。

（3）修改信号产生函数mstg，给每路调幅信号加入载波成分，产生调幅（AM）信号，重复本实验，观察AM信号与抑制载波调幅信号的时域波形及其频谱的差别。

答：分析发现，st的每个频率成分都是25Hz的整数倍。采样频率Fs=10kHz=25×400Hz，即在25Hz的正弦波的1个周期中采样400点。所以，当N为400的整数倍时一定为st的整数个周期。因此，采样点数N=800和N=2000时，对st进行N点FFT可以得到6根理想谱线。如果取N=1000，不是400的整数倍，不能得到6根理想谱线

第二篇：实验四-IIR数字滤波器的设计

实验四 IIR数字滤波器的设计

1、实验目的

（1）        掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的具体方法和原理，熟悉双线性变换法和脉冲响应不变法设计低通、带通IIR数字滤波器的计算机编程；

（2）        观察双线性变换法和脉冲响应不变法设计的数字滤波器的频域特性，了解双线性变换法和脉冲响应不变法的特点和区别；

（3）        熟悉Butterworth滤波器、Chebyshev滤波器和椭圆滤波器的频率特性。

2、实验原理与方法

IIR数字滤波器的设计方法可以概括为如图所示，本实验主要掌握IIR滤波器的第一种方法，即利用模拟滤波器设计IIR数字滤波器，这是IIR数字滤波器设计最常用的方法。利用模拟滤波器设计，需要将模拟域的H_a(s)转换为数字域H(z)，最常用的转换方法为脉冲响应不变法和双线性变换法。

（1）脉冲响应不变法

用数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应h_a(t)，让h(n)正好等于h_a(t)的采样值，即

其中T为采样间隔。如果以H_a(s)及H(z)分别表示h_a(t)的拉氏变换及h(n)的Z变换，则

在MATLAB中，可用函数impinvar实现从模拟滤波器到数字滤波器的脉冲响应不变映射。

（2）双线性变换法

S平面与z平面之间满足下列映射关系

   或 

S平面的虚轴单值地映射于z平面的单位圆上，s平面的左半平面完全映射到z平面的单位圆内。双线性变换不存在频率混叠问题。

在MATLAB中，可用函数bilinear实现从模拟滤波器到数字滤波器的双线性变换映射。

双线性变换是一种非线性变换，即，这种非线性引起的幅频特性畸变可通过预畸变得到校正。

（3）设计步骤

IIR数字滤波器的设计过程中，模拟滤波器的设计是关键。模拟滤波器的设计一般是采用分布设计的方式，这样设计原理非常清楚，具体步骤如前文所述。MATLAB信号处理工具箱也提供了模拟滤波器设计的完全工具函数：butter、cheby1、cheby2、ellip、besself。用户只需一次调用就可完成模拟滤波器的设计，这样虽简化了模拟滤波器的设计过程，但设计原理却被屏蔽了。

模拟滤波器设计完成之后，利用impinvar或bilinear函数将模拟滤波器映射为数字滤波器，即完成了所需数字滤波器的设计。

下图给出了实际低通、高通、带通和带阻滤波器的幅频特性和各截止频率的含义。另外，为了描述过渡带的形状，还引入了通带衰减和阻带衰减的概念。

 

 

图 实际滤波器的幅频特性和各截止频率的含义

通带衰减：    dB

阻带衰减：    dB

在MATLAB信号处理工具箱中，通常用R_p和R_s来表示α_p和α_s。

3、实验内容

（1）参照教材5.5节所述滤波器设计步骤，利用双线性变换法设计一个Chebyshev I型数字高通滤波器，观察通带损耗和阻带衰减是否满足要求。已知滤波器的指标为f_p=0.3kHz，α_p=1.2dB，f_s=0.2kHz，α_s=20dB，T=1ms。

（2）已知f_p=0.2kHz，α_p=1dB，f_s=0.3kHz，α_s=25dB，T=1ms，分别用脉冲响应不变法和双线性变换法设计一个Butterworth数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。

（3）设计一个数字带通滤波器，通带范围为0.25π~0.45π，通带内最大衰减3dB，0.15π以下和0.55π以上为阻带，阻带内最小衰减为15dB，试采用Butterworth或ellip（椭圆）模拟低通滤波器设计。

（4）利用双线性变换法设计一个带宽为0.08π的10阶椭圆带阻滤波器以滤除数字频率为0.44π的信号，选择合适的阻带衰减值，画出幅度响应。产生下面序列的201个样本

，  n=0，2，…，200

并将它通过这个带阻滤波器进行处理（filter函数），讨论所得到的结果。

4、实验报告

（1）简述实验目的和实验原理。

（2）按实验步骤附上所设计的滤波器传递函数H(z)及相应的幅频特性曲线，定性分析所得到的图形，判断设计是否满足要求。

（3）总结脉冲响应不变法和双线性变换法的特点及设计全过程。

（4）收获与建议。

5、实验源程序

%用双线性变换法设计一个Chebyshev型高通滤波器程序如下

Rp=1.2;Rs=20;T=0.001;fp=300;fs=200;

%求出待设计的数字滤波器的边界频率

wp=2*pi*fp*T;ws=2*pi*fs*T;

%预畸变

wp1=(2/T)*tan(wp/2);ws1=(2/T)*tan(ws/2);

%设计模拟滤波器

[n,wn]=cheb1ord(wp1,ws1,Rp,Rs,'s');

[b,a]=cheby1(n,Rp,wn,'high','s');

%双线性变换

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi,db);

axis([0,1,-30,2]);

%用双线性变换法设计一个Butterworth型数字低通滤波器程序如下

Rp=1;Rs=25;T=0.001;fp=300;fs=200;

%求出待设计的数字滤波器的边界频率

wp=2*pi*fp*T;ws=2*pi*fs*T;

%预畸变

wp1=(2/T)*tan(wp/2);ws1=(2/T)*tan(ws/2);

%设计模拟滤波器

[n,wn]=buttord(wp1,ws1,Rp,Rs,'s');

[b,a]=butter(n,wn,'low','s');

%双线性变换

[bz,az]=bilinear(b,a,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi,db);

axis([0,1,-30,2]);

%用脉冲响应不变法设计一个Butterworth数字低通滤波器的程序如下：

%模拟滤波器的技术要求

wp=400*pi;ws=600*pi;Rp=1;Rs=25;

%求模拟滤波器的系统函数

[n,wn]=buttord(wp,ws,Rp,Rs,'s')

[b,a]=butter(n,wn,'s')

%求模拟滤波器的频率响应，w取（0~1000pi）rad/s

[db,mag,pha,w]=freqs_m(b,a,500*2*pi);

%绘图，为了使模坐标显示频率f （单位Hz），将原变量w（模拟角频率，单位为rad/s）进行了处理

plot(w/(2*pi),db,'LineWidth',2,'Color','b');

axis([0,500,-20,1]);hold on

%脉冲响应不变法

fs=1000;[bz,az]=impinvar(b,a,fs);

%求数字滤波器的频率响应

[db,mag,pha,w]=freqz_m(bz,az);

%绘图，为了与模拟滤波器的频响在同一坐标中绘出，需要将数字频率w转换成模拟频率f，转换公式为f=w*fs/2*pi

plot(0.5*fs*w/pi,db,'LineWidth',2,'Color','r');

axis([0,599,-20,1]);hold off

%采用ellip（椭圆）模拟低通滤波器设计，其程序如下：

%确定所需类型数字滤波器的技术指标

Rp=3;Rs=15;T=0.001;

wp1=0.25*pi;wp2=0.45*pi;ws1=0.15*pi;ws2=0.55*pi;

%将所需类型数字滤波器的技术指标转换成模拟滤波器的技术指标

wp3=(2/T)*tan(wp1/2);wp4=(2/T)*tan(wp2/2);

ws3=(2/T)*tan(ws1/2);ws4=(2/T)*tan(ws2/2);

%将所需类型数字滤波器的技术指标转换成模拟滤波器的技术指标,设计模拟滤波器

wp=[wp3,wp4];ws=[ws3,ws4];

[n,wn]=ellipord(wp,ws,Rp,Rs,'s');[z,p,k]=ellipap(n,Rp,Rs);[b,a]=zp2tf(z,p,k);

%频率更换

w0=sqrt(wp3*wp4);Bw=wp4-wp3;

[b1,a1]=lp2bp(b,a,w0,Bw);

%双线性变换法

[bz,az]=bilinear(b1,a1,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi,db);

axis([0,1,-50,2]);

%采用Butterworth模拟低通滤波器设计，其程序如下：

%确定所需类型数字滤波器的技术指标

Rp=3;Rs=15;T=0.001;

wp1=0.25*pi;wp2=0.45*pi;ws1=0.15*pi;ws2=0.55*pi;

%将所需类型数字滤波器的技术指标转换成模拟滤波器的技术指标

wp3=(2/T)*tan(wp1/2);wp4=(2/T)*tan(wp2/2);

ws3=(2/T)*tan(ws1/2);ws4=(2/T)*tan(ws2/2);

%将所需类型数字滤波器的技术指标转换成模拟滤波器的技术指标,设计模拟滤波器

wp=[wp3,wp4];ws=[ws3,ws4];

[n,wn]=buttord(wp,ws,Rp,Rs,'s');[z,p,k]=buttap(n);[b,a]=zp2tf(z,p,k);

%频率更换

w0=sqrt(wp3*wp4);Bw=wp4-wp3;

[b1,a1]=lp2bp(b,a,w0,Bw);

%双线性变换法

[bz,az]=bilinear(b1,a1,1/T);

[db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi,db);

axis([0,1,-50,2]);