实验五、双线性变换法设计IIR 数字滤波器一、实验目的：1、熟悉用双线性变换法设计IIR 数字滤波器的原理与方法。
2、掌握数字滤波器的计算机仿真方法。
3、熟悉Batterworth 滤波器设计方法及特点 二、实验原理（一）、IIR 数字滤波器的设计步骤：① 按照一定规则把给定的滤波器技术指标转换为模拟低通滤波器的技术指标； ② 根据模拟滤波器技术指标设计为响应的模拟低通滤波器；③ 跟据脉冲响应不变法和双线性不变法把模拟滤波器转换为数字滤波器；④ 如果要设计的滤波器是高通、带通或带阻滤波器，则首先把它们的技术指标转化为模拟低通滤波器的技术指标，设计为数字低通滤波器，最后通过频率转换的方法来得到所要的滤波器。
在MATLAB 中，经典法设计IIR 数字滤波器主要采用以下步骤：IIR 数字滤波器设计步骤（二）、用模拟滤波器设计数字滤波器的方法 1、冲激响应不变法：冲激响应不变法是从时域出发，要求数字滤波器的冲激响应h (n ) 对应于模拟滤波器h (t ) 的等间隔抽样。
优点：时域逼近良好；保持线性关系。
缺点：频域响应混叠。
只适用于限带低通滤波器和带通滤波器2、双线性变换法模拟滤波器原型 buttap,cheb1ap频率变换 模拟离散化 bilinear,impin varIIR 数字滤波器/Tπ/T π-3/Tπ3/Tπ-j ΩσjIm (z)Re(z)1S 平面Z 平面1S ~S ππ-将整个平面压缩变换到平面一个的带状区域优点：克服了频域混叠 缺点：高频时会引起畸变1）冲激响应不变法impinvar格式：[BZ,AZ]= impinvar （B,A,Fs ）功能：把具有[B,A]模拟滤波器传递函数模型转换为采样频率为Fs 的数字滤波器的传递函数模型[BZ,AZ]，Fs 默认值为1。
例：一个4阶的Butterworth 模拟低通滤波器的系统函数如下：12251)(234++++=s s s s s H a试用冲激响应不变法求出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+++=s s s s H ，试用双线性变换法求出数字Butterworth 数字低通滤波器的系统函数。
num=1;den=[1,sqrt(3),sqrt(2),1];[num1,den1]=bilinear(num,den,1)3） IIR 数字滤波器的频率变换实现 步骤：① 按一定的规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标② 根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N 和截止频率Wc ③ 利用最小阶数N 产生模拟低通原型④ 利用截止频率Wc 把模拟低通滤波器原型转换为模拟低通、高通、带通、带阻滤波器 ⑤ 利用冲激响应不变法或双线性变换法把模拟滤波器转换为数字滤波器 表一 IIR 滤波器阶次估计表二 模拟低通滤波器原型设计函数名 功能说明buttap Butterworth 模拟低通滤波器原型设计 cleb1ap Chebyshe v Ⅰ模拟低通滤波器原型设计 cheb2ap Chebyshe v Ⅱ模拟低通滤波器原型设计ellipap椭圆模拟低通滤波器原型设计表三 模拟滤波器变换函数函数名 功能说明 lp2bp 模拟低通转换为带通 lp2bs 模拟低通转换为带阻 lp2hp 模拟低通转换为高通 lp2lp改变模拟低通的截止频率（三）、数字滤波器的设计 1． 数字滤波器的设计参数滤波器的4个重要的通带、阻带参数为：p f ：通带截止频率（Hz ） s f ：阻带起始频率（Hz ）p R ：通带内波动（dB ），即通带内所允许的最大衰减；s R ：阻带内最小衰减设采样速率（即奈奎斯特速率）为N f ，将上述参数中的频率参数转化为归一化角频率参数：p ω：通带截止角频率（rad/s ） ，)2//(N p p f f =ω；s ω：阻带起始角频率（rad/s ） ，)2//(N s s f f =ω通过以上参数就可以进行离散滤波器的设计。
2、巴特沃斯滤波器设计1） 巴特沃斯滤波器阶数的选择：在已知设计参数p ω，s ω，p R ，s R 之后，可利用“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)的分母多项式系数(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 ，要求设计一个低通滤波器，p f =2100Hz ，s f =2500Hz ，p R =3dB ，s R =25dB 。
1、采样速率为10000Hz ，要求设计一个巴特沃斯带阻滤波器，p f =[1000Hz ，1500Hz]，s f =[1200Hz ，1300Hz]，p R =3dB ，s R =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([0 4000 -30 3]);grid on pha=angle(H)*180/pi; subplot(2,1,2)plot(F,pha);grid on %画出相频特性图 xlabel('Frequency(Hz)'); ylabel('phase');-30-20-10Frequency(Hz)M a g n i t u d e (d B )-200-100100200Frequency(Hz)p h a s e用双线性变换法wp=2100*2*pi; %利用ws=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++++=s s s s s H a ，编写MATLAB 程序实现从()s H a 设计3dB 截止频率为2π=c w ，设采样周期为T=1，的四阶低通巴特沃斯数字滤波器。
程序如下：步骤一：将设计内容题所给归一化巴特沃斯低通滤波器以3dB 截止频率为2π=c w 进行去归一化。
0000.169048,206568.132262.50000.16)(234++++=s s s s s H a步骤二：用双线性变化法将低通模拟滤波器)(s H a 变换为低通数字滤波器)(z H421210177.04860.010940.03759.05639.03759.00940.0)(-----++++++=zz z z z z H3、已知四阶归一化低通巴特沃斯模拟滤波器系统函数为12251)(234++++=s s s s s H a ，编写MATLAB 程序实现从()s H a 设计3dB 截止频率为4c w π=的四阶高通巴特沃斯数字滤波器。
设计程序如下：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]),grid on;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])Hz幅度归一化模拟低通滤波器Hz幅度去归一化模拟低通滤波器ω/π幅度(d B )数字低通滤波器四、实验内容1、采样速率为10000Hz ，要求设计一个巴特沃斯带阻滤波器，p f =[1000Hz ，1500Hz]，s f =[1200Hz ，1300Hz]，p R =3dB ，s R =30dB 。