数字信号处理实验 matlab版 用双线性变换法设计IIR数字滤波器.docx
《数字信号处理实验 matlab版 用双线性变换法设计IIR数字滤波器.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验 matlab版 用双线性变换法设计IIR数字滤波器.docx（22页珍藏版）》请在冰点文库上搜索。
数字信号处理实验matlab版用双线性变换法设计IIR数字滤波器
实验21用双线性变换法设计IIR数字滤波器
（完美格式版，本人自己完成，所有语句正确，不排除极个别错误，特别适用于山大，勿用冰点等工具下载，否则下载之后的word格式会让很多部分格式错误，）
XXXX学号处XXXX
一、实验目的
1.加深对双线性变换法设计IIR数字滤波器基本方法的了解。
2.掌握用双线性变换法设计数字低通、高通、带通、带阻滤波器的方法。
3.了解MATLAB有关双线性变换法的子函数。
二、实验容
1.双线性变换法的基本知识。
2.用双线性变换法设计IIR数字低通滤波器。
3.用双线性变换法设计IIR数字高通滤波器。
4.用双线性变换法设计IIR数字带通滤波器。
5.用双线性变换法设计IIR数字带阻滤波器。
三、实验环境
MATLAB7.0
四、实验原理
1.实验涉及的MATLAB子函数
Bilinear
功能：
双线性变换——将s域（模拟域）映射到z域（数字域）的标准方法，将模拟滤波器变换成离散等效滤波器。
调用格式：
[numd,dend]=bilinear（num,den,Fs）;将模拟域传递函数变换为数字域传递函数,Fs为取样频率。
[numd,dend]=bilinear（num,den,Fs,Fp）;将模拟域传递函数变换为数字域传递函数,Fs为取样频率,Fp为通带截止频率。
[zd,pd,kd]=bilinear（z,p,k,Fs）;将模拟域零极点增益系数变换到数字域,Fs为取样频率。
[zd,pd,kd]=bilinear（z,p,k,Fs,Fp）;将模拟域零极点增益系数变换到数字域,Fs为取样频率,Fp为通带截止频率。
[Ad,Bd,Cd,Dd]=bilinear（A,B,C,D,Fs）;将模拟域状态变量系数变换到数字域,Fs为取样频率。
2.双线性变换法的基本知识
双线性变换法是将整个s平面映射到整个z平面，其映射关系为
或
双线性变换法克服了脉冲响应不变法从s平面到z平面的多值映射的缺点，消除了频谱混叠现象。
但其在变换过程中产生了非线性的畸变，在设计IIR数字滤波器的过程中需要进行一定的预修正。
用MATLAB双线性变换法进行IIR数字滤波器设计的步骤（参见图19-1）与脉冲响应不变法设计的步骤基本相同：
（1）输入给定的数字滤波器设计指标；
（2）根据公式W＝2[]Ttan[JB（（]w[]2[JB））]进行预修正，将数字滤波器指标转换成模拟滤波器设计指标；
（3）确定模拟滤波器的最小阶数和截止频率；
（4）计算模拟低通原型滤波器的系统传递函数；
（5）利用模拟域频率变换法，求解实际模拟滤波器的系统传递函数；
（6）用双线性变换法将模拟滤波器转换为数字滤波器。
3.用双线性变换法设计IIR数字低通滤波器
例21-1采用双线性变换法设计一个巴特沃斯数字低通滤波器，要求：
wp＝0.25p，Rp＝1dB；ws＝0.4p，As＝15dB，滤波器采样频率Fs＝100Hz。
解程序如下：
wp=0.25*pi;
ws=0.4*pi;
Rp=1;As=15;
ripple=10^（-Rp/20）;Attn=10^（-As/20）;
Fs=100;T=1/Fs;
Omgp=（2/T）*tan（wp/2）;
Omgs=（2/T）*tan（ws/2）;
[n,Omgc]=buttord（Omgp,Omgs,Rp,As,'s'）;
[z0,p0,k0]=buttap（n）;
ba=k0*real（poly（z0））;
aa=real（poly（p0））;
[ba1,aa1]=lp2lp（ba,aa,Omgc）;
%注意,以上4行求滤波器系数ba1、aa1的程序,可由下一条程序替代
%[ba1,aa1]=butter（n,Omgc,'s'）;
[bd,ad]=bilinear（ba1,aa1,Fs）
[sos,g]=tf2sos（bd,ad）
[H,w]=freqz（bd,ad）;
dbH=20*log10（（abs（H）+eps）/max（abs（H）））;
subplot（2,2,1）,plot（w/pi,abs（H））;
ylabel（'|H|'）;title（'幅度响应'）;axis（[0,1,0,1.1]）;
set（gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]）;
set（gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]）;grid
subplot（2,2,2）,plot（w/pi,angle（H）/pi）;
ylabel（'＼phi'）;title（'相位响应'）;axis（[0,1,-1,1]）;
set（gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]）;
set（gca,'YTickMode','manual','YTick',[-1,0,1]）;grid
subplot（2,2,3）,plot（w/pi,dbH）;title（'幅度响应（dB）'）;
ylabel（'dB'）;xlabel（'频率（＼pi）'）;axis（[0,1,-40,5]）;
set（gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]）;
set（gca,'YTickMode','manual','YTick',[-50,-15,-1,0]）;grid
subplot（2,2,4）,zplane（bd,ad）;
axis（[-1.1,1.1,-1.1,1.1]）;title（'零极图'）;
程序运行结果如下：
n＝5
Omgc＝103.2016
bd＝0.00720.03620.07250.07250.03620.0072
ad＝1.0000－1.94341.9680－1.07020.3166－0.0392
sos＝1.00001.003601.0000－0.31930
1.00002.00221.00221.0000－0.69840.2053
1.00001.99420.99421.0000－0.92570.5976
g＝0.0072
频率特性如图21-1所示。
图21-1用双线性变换法设计的巴特沃斯数字低通滤波器特性
由频率特性曲线可知，该设计结果在通阻带截止频率处能满足Rp≤1dB、As≥15dB的设计指标要求，系统的极点全部在单位圆，是一个稳定的系统。
由n＝5可知，设计的巴特沃斯数字低通滤波器是一个5阶的系统，原型Ha（s）在s＝－∞处有5个零点，映射到z＝－1处。
这个巴特沃斯数字低通滤波器的级联型传递函数应为
4.用双线性变换法设计IIR数字高通滤波器
例21-2采用双线性变换法设计一个椭圆数字高通滤波器，要求通带fp＝250Hz，Rp＝1dB；阻带fs＝150Hz，As＝20dB，滤波器采样频率Fs＝1000Hz。
解程序如下：
fs=150;fp=250;Fs=1000;T=1/Fs;
wp=fp/Fs*2*pi;
ws=fs/Fs*2*pi;
Rp=1;As=20;
ripple=10^（-Rp/20）;
Attn=10^（-As/20）;
Omgp=（2/T）*tan（wp/2）;
Omgs=（2/T）*tan（ws/2）;
[n,Omgc]=ellipord（Omgp,Omgs,Rp,As,'s'）
[z0,p0,k0]=ellipap（n,Rp,As）;
ba=k0*real（poly（z0））;
aa=real（poly（p0））;
[ba1,aa1]=lp2hp（ba,aa,Omgc）;
[bd,ad]=bilinear（ba1,aa1,Fs）
[H,w]=freqz（bd,ad）;
dbH=20*log10（（abs（H）+eps）/max（abs（H）））;
subplot（2,2,1）,plot（w/2/pi*Fs,abs（H）,'k'）;
ylabel（'|H|'）;title（'幅度响应'）;axis（[0,Fs/2,0,1.1]）;
set（gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]）;
set（gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]）;grid
subplot（2,2,2）,plot（w/2/pi*Fs,angle（H）/pi*180,'k'）;
ylabel（'\phi'）;title（'相位响应'）;axis（[0,Fs/2,-180,180]）;
set（gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]）;
set（gca,'YTickMode','manual','YTick',[-180,0,180]）;grid
subplot（2,2,3）,plot（w/2/pi*Fs,dbH）;
title（'幅度响应（dB）'）;axis（[0,Fs/2,-40,5]）;
ylabel（'dB'）;xlabel（'频率（\pi）'）;
set（gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]）;
set（gca,'YTickMode','manual','YTick',[-50,-20,-1,0]）;grid
subplot（2,2,4）,zplane（bd,ad）;
axis（[-1.1,1.1,-1.1,1.1]）;title（'零极图'）;
运行结果如下：
n=3
Omgc=2.0000e+003
bd=0.2545-0.43220.4322-0.2545
ad=1.00000.18900.71970.1574
频率特性如图21-2所示。
图21-2用双线性变换法设计椭圆高通数字滤波器
由频率特性曲线可知，该设计结果在通阻带截止频率处能满足Rp≤1dB、As≥20dB的设计指标要求。
由n＝3可知，设计的椭圆数字高通滤波器是一个3阶的系统，极点全部在Z平面的单位圆，是一个稳定的系统。
这个高通滤波器的传递函数应为
5.用双线性变换法设计IIR数字带通滤波器
例21-3采用双线性变换法设计一个切比雪夫Ⅰ型数字带通滤波器，要求：
通带wp1＝0.3p，wp2＝0.7p，Rp＝1dB；阻带ws1＝0.2p，ws2＝0.8p，As＝20dB，滤波器采样周期Ts＝0.001s。
解程序如下：
wp1=0.3*pi;wp2=0.7*pi;
ws1=0.2*pi;ws2=0.8*pi;
Rp=1;As=20;
T=0.001;Fs=1/T;
Omgp1=（2/T）*tan（wp1/2）;Omgp2=（2/T）*tan（wp2/2）;
Omgp=[Omgp1,Omgp2];
Omgs1=（2/T）*tan（ws1/2）;Omgs2=（2/T）*tan（ws2/2）;
Omgs=[Omgs1,Omgs2];
bw=Omgp2-Omgp1;w0=sqrt（Omgp1*Omgp2）;
[n,Omgn]=cheb1ord（Omgp,Omgs,Rp,As,'s'）
[z0,p0,k0]=cheb1ap（n,Rp）;
ba1=k0*real（poly（z0））;
aa1=real（poly（p0））;
[ba,aa]=lp2bp（ba1,aa1,w0,bw）;
[bd,ad]=bilinear（ba,aa,Fs）
[H,w]=freqz（bd,ad）;
dbH=20*log10（（abs（H）+eps）/max（abs（H）））;
subplot（2,2,1）,plot（w/2/pi*Fs,abs（H）,'k'）;
ylabel（'|H|'）;title（'幅度响应'）;axis（[0,Fs/2,0,1.1]）;
set（gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]）;
set（gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]）;grid
subplot（2,2,2）,plot（w/2/pi*Fs,angle（H）/pi*180,'k'）;
ylabel（'\phi'）;title（'相位响应'）;axis（[0,Fs/2,-180,180]）;
set（gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]）;
set（gca,'YTickMode','manual','YTick',[-180,0,180]）;grid
subplot（2,2,3）,plot（w/2/pi*Fs,dbH）;
title（'幅度响应（dB）'）;axis（[0,Fs/2,-40,5]）;
ylabel（'dB'）;xlabel（'频率（\pi）'）;
set（gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]）;
set（gca,'YTickMode','manual','YTick',[-50,-20,-1,0]）;grid
subplot（2,2,4）,zplane（bd,ad）;
axis（[-1.1,1.1,-1.1,1.1]）;title（'零极图'）;
程序运行结果如下：
n=3
Omgn=1.0e+003*1.01913.9252
bd=0.07360.0000-0.22080.00000.2208-0.0000-0.0736
ad=1.00000.00000.97610.00000.85680.00000.2919
频率特性及零极点图形如图21-3所示。
图21-3用双线性变换法设计切比雪夫Ⅰ型带通数字滤波器
由频率特性曲线可知，该设计结果在通阻带截止频率处能满足Rp≤1dB、As≥20dB的设计指标要求。
由n＝3可知，由3阶的模拟低通原型用双线性变换法设计出来的切比雪夫Ⅰ型数字带通滤波器是一个6阶的系统，极点全部在Z平面的单位圆，是一个稳定的系统。
这个滤波器的传递函数应为
注意：
在使用［z0，p0，k0］＝cheb1ap（n，Rp）设计模拟低通原型时，需要输入通带衰减Rp，即切比雪夫Ⅰ型模拟低通原型是以通带衰减Rp为主要设计指标的。
因此，由模拟低通原型变为数字带通（或带阻）滤波器时，使用lp2bp（或lp2bs）语句要求输入模拟通带带宽W0和中心频率BW，应采用通带截止频率来计算，即
bw=Omgp2-Omgp1;w0=sqrt（Omgp1*Omgp2）;%[ZK（]模拟滤波器通带带
宽和中心频率
如果将例21-3改为：
采用双线性变换法设计一个切比雪夫Ⅱ型数字带通滤波器,其它条件不变,则需要修改下面几句程序：
bw=Omgs2-Omgs1;w0=sqrt（Omgs1*Omgs2）;%[ZK（]模拟滤波器阻带带
宽和中心频率
[n,Omgn]=cheb2ord（Omgp,Omgs,Rp,As,'s'）%计算阶数n和截止频率
[z0,p0,k0]=cheb2ap（n,As）;%设计归一化的模拟原型滤波器
采用阻带截止频率来计算W0和BW，是因为切比雪夫Ⅱ型模拟低通原型是以阻带衰减As为主要设计指标的。
程序运行结果如下：
n=3
Omgn=1.0e+003*
1.66416.0234
bd=0.25370.0000-0.47330.00000.47330.0000-0.2537
ad=1.0000-0.00000.00080.00000.4206-0.0000-0.0343
频率特性及零极点图形如图21-4所示。
图21-4用双线性变换法设计切比雪夫Ⅱ型带通数字滤波器
由程序数据和曲线可知，该设计结果在通阻带截止频率处能满足Rp≤1dB、As≥20dB的设计指标要求。
由n＝3可知，由3阶的模拟低通原型用双线性变换法设计出来的切比雪夫Ⅱ型数字带通滤波器是一个6阶的系统，极点全部在z平面的单位圆，是一个稳定的系统。
这个滤波器的传递函数应为
6.用双线性变换法设计IIR数字带阻滤波器
例21-4采用双线性变换法设计一个切比雪夫Ⅰ型数字带阻滤波器，要求：
下通带wp1＝0.2p，上通带wp2＝0.8p，Rp＝1dB；阻带下限ws1＝0.3p，阻带上限ws2＝0.7p，As＝20dB，滤波器采样频率Fs＝1000Hz。
解由题目可知，本例只是将例21-3的条件改为相反，即将原带通滤波器通带的频率区域改为带阻滤波器阻带的频率区域，将原带通滤波器阻带的频率区域改为带阻滤波器通带的频率区域。
程序只需作5句修改：
ws1=0.3*pi;ws2=0.7*pi;
wp1=0.2*pi;wp2=0.8*pi;
[ba,aa]=lp2bs（ba1,aa1,w0,bw）;
程序运行结果如下：
n=3
Omgn=1.0e+003*
0.64986.1554
bd=0.0736-0.00000.22080.00000.2208-0.00000.0736
ad=1.00000.0000-0.9761-0.00000.85680.0000-0.2919
频率特性及零极点图形如图21-5所示。
图21-5用双线性变换法设计切比雪夫Ⅰ型带阻数字滤波器
由程序数据和曲线可知，该设计结果在通阻带截止频率处能满足Rp≤1dB、As≥20dB的设计指标要求。
由3阶的模拟低通原型用双线性变换法设计出来的切比雪夫Ⅰ型数字带阻滤波器是一个6阶的系统，极点全部在z平面的单位圆，是一个稳定的系统。
这个滤波器的传递函数应为
如果将例21-4改为：
采用双线性变换法设计一个切比雪夫Ⅱ型数字带阻滤波器，其它条件不变，则在上面程序的基础上与例21-3一样，修改下面几句程序：
bw=Omgs2-Omgs1;w0=sqrt（Omgs1*Omgs2）;%[ZK（]模拟滤波器阻带带宽
和中心频率
[n,Omgn]=cheb2ord（Omgp,Omgs,Rp,As,'s'）[KG*2]%计算阶数n和截止
频率
[z0,p0,k0]=cheb2ap（n,As）;%设计归一化的模拟原型滤波器
程序运行结果如下：
n＝3
Omgn＝10.029339.8833
bd＝0.25370.00000.4733－0.00000.47330.00000.2537
ad＝1.0000－0.0000－0.0008－0.00000.4206－0.00000.0343
频率特性及零极点图形如图21-6所示。
图21-6用双线性变换法设计切比雪夫Ⅱ型带阻数字滤波器
由程序数据和曲线可知，该设计结果在通阻带截止频率处能满足Rp≤1dB、As≥20dB的设计指标要求。
由3阶的模拟低通原型用双线性变换法设计出来的切比雪夫Ⅱ型数字带阻滤波器是一个6阶的系统，极点全部在z平面的单位圆，是一个稳定的系统。
这个滤波器的传递函数应为
五、实验过程
1.用双线性变换法设计切比雪夫Ⅱ型数字滤波器，列出传递函数并描绘模拟和数字滤波器的幅频响应曲线。
（1）设计一个数字低通，要求：
通带wp＝0.2p，Rp＝1dB；阻带ws＝0.35p，As＝15dB，滤波器采样频率Fs＝10Hz。
（2）设计一个数字高通，要求：
通带wp＝0.35p，Rp＝1dB；阻带ws＝0.2p，As＝15dB，滤波器采样频率Fs＝10Hz。
解
（1）MATLAB程序如下：
wp=0.2*pi;
ws=0.35*pi;
Rp=1;As=15;
ripple=10^（-Rp/20）;Attn=10^（-As/20）;
Fs=10;T=1/Fs;
Omgp=（2/T）*tan（wp/2）;
Omgs=（2/T）*tan（ws/2）;
[n,Omgc]=cheb2ord（Omgp,Omgs,Rp,As,'s'）
[z0,p0,k0]=cheb2ap（n,As）;
ba=k0*real（poly（z0））;
aa=real（poly（p0））;
[ba1,aa1]=lp2lp（ba,aa,Omgc）;
[bd,ad]=bilinear（ba1,aa1,Fs）
[sos,g]=tf2sos（bd,ad）
[H,w]=freqz（bd,ad）;
dbH=20*log10（（abs（H）+eps）/max（abs（H）））;
subplot（2,2,1）,plot（w/pi,abs（H））;
ylabel（'|H|'）;title（'幅度响应'）;axis（[0,1,0,1.1]）;
set（gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]）;
set（gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]）;grid
subplot（2,2,2）,plot（w/pi,angle（H）/pi）;
ylabel（'\phi'）;title（'相位响应'）;axis（[0,1,-1,1]）;
set（gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]）;
set（gca,'YTickMode','manual','YTick',[-1,0,1]）;grid
subplot（2,2,3）,plot（w/pi,dbH）;title（'幅度响应（dB）'）;
ylabel（'dB'）;xlabel（'频率'）;axis（[0,1,-40,5]）;
set（gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]）;
set（gca,'YTickMode','manual','YTick',[-50,-15,-1,0]）;grid
subplot（2,2,4）,zplane（bd,ad）;
axis（[-1.1,1.1,-1.1,1.1]）;title（'零极图'）;
运行结果如图1-1所示。
图1-1
（2）MATLAB程序如下：
Fs=1000;T=1/Fs;
wp=0.35*pi;
ws=0.2*pi;
fp=wp/（2*pi）*Fs;
fs=ws/（2*pi）*Fs;
Rp=1;As=15;
ripple=10^（-Rp/20）;
Attn=10^（-As/20）;
Omgp=（2/T）*tan（wp/2）;
Omgs=（2/T）*tan（ws/2）;
[n,Omgc]=cheb2ord（Omgp,Omgs,Rp,As,'s'）
[z0,p0,k0]=cheb2ap（n,As）;
ba=k0*real（poly（z0））;
aa=real（poly（p0））;
[ba1,aa1]=lp2hp（ba,aa,Omgc）;
[bd,ad]=bilinear（ba1,aa1,Fs）
[H,w]=freqz（bd,ad）;
dbH=20*log10（（abs（H）+eps）/m
升级会员 