数字信号处理--用双线性变换法设计IIR 数字滤波器

1、 数字信号处理实验内容用双线性变换法设计IIR 数字滤波器一、实验目的：1、熟悉用双线性变换法设计IIR 数字滤波器的原理与方法。2、掌握数字滤波器的计算机仿真方法。3、熟悉Batterworth滤波器设计方法及特点二、实验原理：（一）、IIR数字滤波器的设计步骤： 按照一定规则把给定的滤波器技术指标转换为模拟低通滤波器的技术指标； 根据模拟滤波器技术指标设计为响应的模拟低通滤波器； 跟据脉冲响应不变法和双线性不变法把模拟滤波器转换为数字滤波器； 如果要设计的滤波器是高通、带通或带阻滤波器，则首先把它们的技术指标转化为模拟低通滤波器的技术指标，设计为数字低通滤波器，最后通过频率转换的方法来得到

2、所要的滤波器。在MATLAB中，经典法设计数字滤波器主要采用以下步骤：模拟滤波器原型buttap,cheb1ap频率变换模拟离散化bilinear,impinvarIIR数字滤波器IIR数字滤波器设计步骤（二）、用模拟滤波器设计数字滤波器的方法1、冲激响应不变法：冲激响应不变法是从时域出发，要求数字滤波器的冲激响应h(n) 对应于模拟滤波器h(t) 的等间隔抽样。优点：时域逼近良好；保持线性关系。缺点：频域响应混叠。只适用于限带低通滤波器和带通滤波器2、双线性变换法优点：克服了频域混叠缺点：高频时会引起畸变1）冲激响应不变法impinvar格式：BZ,AZ= impinvar（B,A,Fs）功

3、能：把具有B,A模拟滤波器传递函数模型转换为采样频率为Fs的数字滤波器的传递函数模型BZ,AZ，Fs默认值为1。例：一个4阶的Butterworth模拟低通滤波器的系统函数如下：试用冲激响应不变法求出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,de

4、n,Fs)功能：把模拟滤波器的传递函数模型转换为数字滤波器的传递函数模型。例：一个三阶的模拟Butterworth模拟低通滤波器的系统函数如下：，试用双线性变换法求出数字Butterworth数字低通滤波器的系统函数。num=1;den=1,sqrt(3),sqrt(2),1;num1,den1=bilinear(num,den,1)3） IIR数字滤波器的频率变换实现步骤： 按一定的规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标 根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N和截止频率Wc 利用最小阶数N产生模拟低通原型 利用截止频率Wc把模拟低通滤波器原型转换为模

5、拟低通、高通、带通、带阻滤波器 利用冲激响应不变法或双线性变换法把模拟滤波器转换为数字滤波器表一 IIR滤波器阶次估计函数名功能说明buttord计算Butterworth滤波器的阶次及截止频率cheb1ord计算Chebyshev滤波器的阶次cheb2ord计算Chebyshev滤波器的阶次ellipord计算椭圆滤波器的最小阶次表二 模拟低通滤波器原型设计函数名功能说明buttapButterworth模拟低通滤波器原型设计cleb1apChebyshev模拟低通滤波器原型设计cheb2apChebyshev模拟低通滤波器原型设计ellipap椭圆模拟低通滤波器原型设计表三 模拟滤波器变换

6、函数函数名功能说明lp2bp模拟低通转换为带通lp2bs模拟低通转换为带阻lp2hp模拟低通转换为高通lp2lp改变模拟低通的截止频率（三）、数字滤波器的设计1 数字滤波器的设计参数滤波器的4个重要的通带、阻带参数为：通带截止频率（Hz） ：阻带起始频率（Hz）：通带内波动（dB），即通带内所允许的最大衰减；：阻带内最小衰减设采样速率（即奈奎斯特速率）为，将上述参数中的频率参数转化为归一化角频率参数：通带截止角频率（rad/s） ，； ：阻带起始角频率（rad/s） ，通过以上参数就可以进行离散滤波器的设计。2、巴特沃斯滤波器设计1） 巴特沃斯滤波器阶数的选择：在已知设计参数，之后，可利用“b

7、uttord”命令可求出所需要的滤波器的阶数和3dB截止频率，其格式为：n，Wn=buttordWp，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，W

8、n，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，要求设计一个低通滤波器，=2100Hz，=2500Hz，=3dB，=25dB。用直接设计法程序如下：fn=8000;%采样频率 fp=2100; %通带截止频率fs=2500; %

9、阻带起始频率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)&#

10、39;) title('低通滤波器')axis(0 4000 -30 3);grid onpha=angle(H)*180/pi;subplot(2,1,2)plot(F,pha);grid on %画出相频特性图xlabel('Frequency(Hz)'); ylabel('phase');用双线性变换法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

11、模拟滤波器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：模拟原型直接变换法设计数字滤波器：已知四阶归一化低通巴特沃斯模拟滤波器系统函数为，编写

12、MATLAB程序实现从设计3dB截止频率为，设采样周期为T=1，的四阶低通巴特沃斯数字滤波器。程序如下：步骤一：将设计内容题所给归一化巴特沃斯低通滤波器以3dB截止频率为进行去归一化。步骤二：用双线性变化法将低通模拟滤波器变换为低通数字滤波器设计程序如下：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

13、(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('去归一化模拟低通滤波器'

14、);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)四、实验内容：1、采样速率为10000Hz，要求设计一个巴特沃斯带阻滤波器，=1000Hz，1500Hz，=1200Hz，1300Hz，=3dB，=30dB。提示：b,a=butter(N,Wc,

15、9;stop')2、 采样速率为10000Hz，要求设计一个带通滤波器，=1000Hz，1500Hz，=600Hz，1900Hz，=3dB，=20dB。3、已知四阶归一化低通巴特沃斯模拟滤波器系统函数为，编写MATLAB程序实现从设计3dB截止频率为的四阶高通巴特沃斯数字滤波器。4、设计低通滤波器，把输入信号中的三个信号分离出来。要求：画出滤波前后信号的波形及频谱及低通滤波器的幅频响应。提示：b,a=butter(n,Wn);%计算H(z)分子、分母多项式系数y=filter(b,a,x);%对输入的信号进行滤波五、思考题：用双线性变换法设计数字滤波器过程中，变换公式中T 的取值对设计

16、结果有无影响？为什么？程序如下：五、实验结果：1、程序如下：fn=10000;%采样频率 fp=1000,1500; %通带截止频率fs=600,1900; %阻带起始频率Rp=3; %通带最大衰减Rs=20;%阻带最小衰减Wp=fp/(fn/2);%计算归一化角频率Ws=fs/(fn/2);n,Wn=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率b,a=butter(n,Wn,'stop');%计算H(z)分子、分母多项式系数H,F=freqz(b,a,1000,8000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)subplot(2,

17、1,1)plot(F,20*log10(abs(H) %画出幅频特性图xlabel('Frequency(Hz)'); ylabel('Magnitude(dB)') title('带阻滤波器')axis(0 4000 -30 3);grid onpha=angle(H)*180/pi;subplot(2,1,2)plot(F,pha);grid on %画出相频特性图xlabel('Frequency(Hz)'); ylabel('phase');2、程序如下：fn=10000;%采样频率 fp=1000,150

18、0; %通带截止频率fs=600,1900; %阻带起始频率Rp=3; %通带最大衰减Rs=20;%阻带最小衰减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)'); yla

19、bel('Magnitude(dB)') title('带通滤波器')axis(0 4000 -30 3);grid onpha=angle(H)*180/pi;subplot(2,1,2)plot(F,pha);grid on %画出相频特性图xlabel('Frequency(Hz)'); ylabel('phase');3、程序如下：clear; T=1; fs=1/T; N=4;wc=pi/4; omegach=2*tan(wc/2)/T;%模拟滤波器的截止频率M=1; N=1,sqrt(5),2,sqrt(2),1;h,w=freqs(M,N,512); %模拟滤波器的幅频响应subplot(2,1,1);plot(w,20*log10(abs(h); axis