实验五：双线性变换法的设计IIR-数字滤波器

1实验五、双线性变换法设计IIR数字滤波器一、实验目的：1、熟悉用双线性变换法设计IIR数字滤波器的原理与方法。2、掌握数字滤波器的计算机仿真方法。3、熟悉Batterworth滤波器设计方法及特点二、实验原理（一）、IIR数字滤波器的设计步骤：①按照一定规则把给定的滤波器技术指标转换为模拟低通滤波器的技术指标；②根据模拟滤波器技术指标设计为响应的模拟低通滤波器；③跟据脉冲响应不变法和双线性不变法把模拟滤波器转换为数字滤波器；④如果要设计的滤波器是高通、带通或带阻滤波器，则首先把它们的技术指标转化为模拟低通滤波器的技术指标，设计为数字低通滤波器，最后通过频率转换的方法来得到所要的滤波器。在MATLAB中，经典法设计IIR数字滤波器主要采用以下步骤：IIR数字滤波器设计步骤（二）、用模拟滤波器设计数字滤波器的方法1、冲激响应不变法：冲激响应不变法是从时域出发，要求数字滤波器的冲激响应h(n)对应于模拟滤波器h(t)的等间隔抽样。优点：时域逼近良好；保持线性关系。缺点：频域响应混叠。只适用于限带低通滤波器和带通滤波器2、双线性变换法模拟滤波器原型buttap,cheb1ap频率变换模拟离散化bilinear,impinvarIIR数字滤波器/T/T3/T3/TjjIm(z)Re(z)1S平面Z平面1S~STT将整个平面压缩变换到平面一个的带状区域2优点：克服了频域混叠缺点：高频时会引起畸变1）冲激响应不变法impinvar格式：[BZ,AZ]=impinvar（B,A,Fs）功能：把具有[B,A]模拟滤波器传递函数模型转换为采样频率为Fs的数字滤波器的传递函数模型[BZ,AZ]，Fs默认值为1。例：一个4阶的Butterworth模拟低通滤波器的系统函数如下：12251)(234sssssHa试用冲激响应不变法求出Butterworth模拟低通数字滤波器的系统函数。num=1;den=[1,sqrt(5),2,sqrt(2),1];[num1,den1]=impinvar(num,den)2）双线性变换法bilinear格式一：[Zd,Pd,Kd]=bilinear(Z,P,K,Fs)功能：把模拟滤波器的零极点模型转换成数字滤波器的零极点模型，Fs是采样频率格式二：[numd,dend]=bilinear(num,den,Fs)功能：把模拟滤波器的传递函数模型转换为数字滤波器的传递函数模型。例：一个三阶的模拟Butterworth模拟低通滤波器的系统函数如下：1231)(23ssssH，试用双线性变换法求出数字Butterworth数字低通滤波器的系统函数。num=1;den=[1,sqrt(3),sqrt(2),1];[num1,den1]=bilinear(num,den,1)3）IIR数字滤波器的频率变换实现步骤：①按一定的规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标②根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N和截止频率Wc③利用最小阶数N产生模拟低通原型④利用截止频率Wc把模拟低通滤波器原型转换为模拟低通、高通、带通、带阻滤波器⑤利用冲激响应不变法或双线性变换法把模拟滤波器转换为数字滤波器表一IIR滤波器阶次估计函数名功能说明buttord计算Butterworth滤波器的阶次及截止频率cheb1ord计算ChebyshevⅠ滤波器的阶次cheb2ord计算ChebyshevⅡ滤波器的阶次ellipord计算椭圆滤波器的最小阶次/T/TjjIm(z)Re(z)1S平面Z平面jS平面1113表二模拟低通滤波器原型设计函数名功能说明buttapButterworth模拟低通滤波器原型设计cleb1apChebyshevⅠ模拟低通滤波器原型设计cheb2apChebyshevⅡ模拟低通滤波器原型设计ellipap椭圆模拟低通滤波器原型设计表三模拟滤波器变换函数函数名功能说明lp2bp模拟低通转换为带通lp2bs模拟低通转换为带阻lp2hp模拟低通转换为高通lp2lp改变模拟低通的截止频率（三）、数字滤波器的设计1．数字滤波器的设计参数滤波器的4个重要的通带、阻带参数为：pf：通带截止频率（Hz）sf：阻带起始频率（Hz）pR：通带内波动（dB），即通带内所允许的最大衰减；sR：阻带内最小衰减设采样速率（即奈奎斯特速率）为Nf，将上述参数中的频率参数转化为归一化角频率参数：p：通带截止角频率（rad/s），)2//(Nppff；s：阻带起始角频率（rad/s），)2//(Nssff通过以上参数就可以进行离散滤波器的设计。2、巴特沃斯滤波器设计1）巴特沃斯滤波器阶数的选择：在已知设计参数p，s，pR，sR之后，可利用“buttord”命令可求出所需要的滤波器的阶数和3dB截止频率，其格式为：[n，Wn]=buttord[Wp，Ws，Rp，Rs]，其中Wp，Ws，Rp，Rs分别为通带截止频率、阻带起始频率、通带内波动、阻带内最小衰减。返回值n为滤波器的最低阶数，Wn为3dB截止频率。2）巴特沃斯滤波器系数计算：由巴特沃斯滤波器的阶数n以及3dB截止频率Wn可以计算出对应传递函数H(z)的分子分母系数，MATLAB提供的命令如下：(a)巴特沃斯低通滤波器系数计算：[b，a]=butter(n,Wn)，其中b为H(z)的分子多项式系数，a为H(z)的分母多项式系数4(b)巴特沃斯高通滤波器系数计算：[b，a]=butter(n，Wn，’High’)(c)巴特沃斯带通滤波器系数计算：[b，a]=butter(n，[W1，W2])，其中[W1，W2]为截止频率，是2元向量，需要注意的是该函数返回的是2*n阶滤波器系数。(d)巴特沃斯带阻滤波器系数计算：[b，a]=butter(ceil(n/2)，[W1，W2]，’stop’)，其中[W1，W2]为截止频率，是2元向量，需要注意的是该函数返回的也是2*n阶滤波器系数。三、巴特沃斯滤波器设计实例：例题1：采样速率为8000Hz，要求设计一个低通滤波器，pf=2100Hz，sf=2500Hz，pR=3dB，sR=25dB。1、采样速率为10000Hz，要求设计一个巴特沃斯带阻滤波器，pf=[1000Hz，1500Hz]，sf=[1200Hz，1300Hz]，pR=3dB，sR=30dB。用直接设计法程序如下：fn=8000;%采样频率fp=2100;%通带截止频率fs=2500;%阻带起始频率Rp=3;%通带最大衰减Rs=25;%阻带最小衰减Wp=fp/(fn/2);%计算归一化角频率Ws=fs/(fn/2);[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率[b,a]=butter(n,Wn);%计算H(z)分子、分母多项式系数[H,F]=freqz(b,a,1000,8000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)subplot(2,1,1)plot(F,20*log10(abs(H)))%画出幅频特性图xlabel('Frequency(Hz)');ylabel('Magnitude(dB)')title('低通滤波器')axis([04000-303]);gridonpha=angle(H)*180/pi;subplot(2,1,2)plot(F,pha);gridon%画出相频特性图xlabel('Frequency(Hz)');ylabel('phase');05001000150020002500300035004000-30-20-100Frequency(Hz)Magnitude(dB)低通滤波器05001000150020002500300035004000-200-1000100200Frequency(Hz)phase用双线性变换法wp=2100*2*pi;%利用5ws=2500*2*pi;Rp=3;Rs=25;Fs=8000;Ts=1/Fs;%选择滤波器的最小阶数[N,Wn]=buttord(wp,ws,Rp,Rs,'s');%创建butterworth模拟滤波器[Z,P,K]=buttap(N);%把滤波器零极点模型转化为传递函数模型[Bap,Aap]=zp2tf(Z,P,K);%把模拟滤波器原型转换成截至频率为Wn的低通滤波器[b,a]=lp2lp(Bap,Aap,Wn);%用双线性变换法实现模拟滤波器到数字滤波器的转换[bz,az]=bilinear(b,a,Fs);%绘制频率响应曲线[H,W]=freqz(bz,az);plot(W*Fs/(2*pi),abs(H));gridxlabel('频率/Hz')ylabel('幅度')例题2：模拟原型直接变换法设计数字滤波器：已知四阶归一化低通巴特沃斯模拟滤波器系统函数为16131.24142.36131.21234sssssHa，编写MATLAB程序实现从sHa设计3dB截止频率为2cw，设采样周期为T=1，的四阶低通巴特沃斯数字滤波器。程序如下：步骤一：将设计内容题所给归一化巴特沃斯低通滤波器以3dB截止频率为2cw进行去归一化。0000.169048,206568.132262.50000.16)(234sssssHa步骤二：用双线性变化法将低通模拟滤波器)(sHa变换为低通数字滤波器)(zH6421210177.04860.010940.03759.05639.03759.00940.0)(zzzzzzH3、已知四阶归一化低通巴特沃斯模拟滤波器系统函数为12251)(234sssssHa，编写MATLAB程序实现从sHa设计3dB截止频率为4cw的四阶高通巴特沃斯数字滤波器。设计程序如下：clear;T=1;fs=1/T;N=4;wc=pi/2;omegach=2*tan(wc/2)/T;%模拟滤波器的截止频率M=1;N=[1,2.6131,3.4142,2.6131,1];[h,w]=freqs(M,N,512);%模拟滤波器的幅频响应subplot(2,1,1);plot(w,20*log10(abs(h)));axis([0,10,-90,0]),gridon;xlabel('Hz');ylabel('幅度');title('归一化模拟低通滤波器');[Ms,Ns]=lp2lp(M,N,omegach);%对低通滤波器进行频率变换[hs,ws]=freqs(Ms,Ns,512);%模拟滤波器的幅频响应subplot(2,1,2);plot(ws,20*log10(abs(hs)));grid;axis([0,10,-90,0]);xlabel('Hz');ylabel('幅度');title('去归一化模拟低通滤波器');[Mz,Nz]=bilinear(Ms,Ns,1/T);%对模拟滤波器双线性变换[h1,w1]=freqz(Mz,Nz);%数字滤波器的幅频响应figureplot(w1/pi,20*log10(abs(h1)));grid;xlabel('ω/π');ylabel('幅度(dB)');title('数字低通滤波器');axis([0,1,-160,0])7012345678910-500Hz幅度归一化模拟低通滤波器012345678910-500Hz幅度去归一化模拟低通滤波器00.10.20.30.40.50.60.70.80.91-150-100-500ω/π幅度(dB)数字低通滤波器四、实验内容1、采样速率为10000Hz，要求设计一个巴特沃斯带阻滤波器，pf=[1000Hz，1500Hz]，sf=[1200Hz，1300Hz]，pR=3dB，sR=30dB。程序：fn=10000;%采样频率fp=[1000,1500];%通带截止频率fs=[1200,1300];%阻带起始频率Rp=3;%通带最大衰减Rs=30;%阻带最小衰减Wp