实验三IIR数字滤波器设计

数字信号处理实验报告1实验三IIR数字滤波器设计一、实验目的1、掌握利用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及具体方法。2、掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及使用范围。二、实验设备与环境计算机、MATLAB软件环境三、实验基础理论1、数字滤波器和模拟滤波器的一些指标𝛿1为通带响应中的容限，𝛿2为阻带的容限，𝑅𝑝为通带波动，𝐴𝑠为阻带衰减，且𝑅𝑝=−20𝑙𝑜𝑔101−𝛿11+𝛿1𝐴𝑠=−20𝑙𝑜𝑔10𝛿21+𝛿1ε为通带波动系数，Ω𝑝是单位为rad/s的带通截止频率，A为以dB为单位的阻带衰减参数，Ω𝑠𝑡是单位为rad/s的阻带截止频率。参数ε和A与滤波器的带通波动𝑅𝑝及阻带衰减𝐴𝑠之间有如下关系Rp=−10log1011+ε2ε=√10Rp/10−1𝐴𝑠=−10𝑙𝑜𝑔101𝐴2𝐴=10𝐴𝑠/102、模拟原型滤波器（1）巴特沃斯滤波器设计方法一：在一般设计中，先把Ω𝑐选择为1rad/s，使频率得到归一化，归一化后的巴特沃斯滤波器系统函数的极点分布及分母多项式系数等都有现成的表格可查。然后将归一化巴特沃斯滤波器系统函数中的变量s用s/Ω𝑐替换以后，即可得到任意Ω𝑐的非归一化的巴特沃斯滤波器。由于巴特沃斯滤波器系统函数不存在零点，则将归一化崩塌沃斯滤波器系统函数的极点倍乘Ω𝑐，增益倍乘以Ω𝑐，即可实现上述由归一化巴特沃斯滤波器到非归一化巴特沃斯滤波器的变换。方法二：在MATLAB中还可以直接利用butter函数设计巴特沃斯模拟原型滤波器，butter函数既可以用于设计模拟巴特沃斯滤波器也可用于设计数字巴特沃斯滤波器，用于设计模拟巴特沃斯滤波器时，其调用格式为[b,a]=butter(N,Wc,’s’)，N表示滤波器的阶数，Wc表示截止频率，‘s’表示设计模拟滤波器。设计方程N=⌈𝑙𝑜𝑔10[(10𝑅𝑝10−1)/(10𝐴𝑠10−1)]2𝑙𝑜𝑔10(Ω𝑝/Ω𝑠𝑡)⌉由于N比实际需要的要大，所以在Ω𝑝或Ω𝑠𝑡处将满足或超过指标，一般选择下式结果作为Ω𝑐的值，Ω𝑐=Ω𝑝√10𝑅𝑝10−12𝑁数字信号处理实验报告2（2）切比雪夫低通滤波器切比雪夫I型：方法一：MATLAB提供了函数[z,p,k]=cheb1ap(N,𝑅𝑝)来设计一个阶数为N,通带波动为𝑅𝑝的归一化（截止频率为1）切比雪夫I型模拟原型滤波器，它返回数组z和零极点p以及增益k。切比雪夫I型滤波器系统函数不存在零点，将归一化切比雪夫I型滤波器系统函数的极点倍乘以Ω𝑐，增益倍乘以s=0时非归一化滤波器系统函数分母多项式和归一化滤波器系统函数多项式的壁纸，即可实现由归一化切比雪夫I型滤波器到非归一化切比雪夫I型滤波器的变换。方法二：在MATLAB中还可以直接利用cheby1函数设计切比雪夫I型模拟原型滤波器。cheby1函数既可用于设计模拟切比雪夫I型滤波器也可用于设计数字切比雪夫I型滤波器，用于设计模拟切比雪夫I型滤波器。设计方程：ε=√10𝑅𝑝/10−1𝐴=10𝐴𝑠/10Ω𝑐=Ω𝑝N=⌈arccosh(√𝐴2−1/ε)arccosh(Ω𝑠𝑡/Ω𝑝)⌉切比雪夫II型：方法一：MATLAB提供了函数[z,p,k]=cheb2ap(N,𝐴𝑠)来设计一个阶数为N，阻带衰减为𝐴𝑠的归一化（截止频率为1）切比雪夫II型模拟原型滤波器，它返回数组z和零极点p以及增益k。切比雪夫II型滤波器系统函数存在零点，可以将归一化切比雪夫II型滤波器系统函数的极点和零点倍乘Ω𝑐，增益倍乘s=0时非归一化的有理函数和归一化的有理函数的比值，即可实现由归一化切比雪夫II型滤波器到非归一化切比雪夫II型滤波器的变换。方法二：在MATLAB中还可以直接利用cheby2函数设计切比雪夫II型模拟原型滤波器。cheby2函数既可用于设计模拟切比雪夫II型滤波器也可用于设计数字切比雪夫II型滤波器，用于设计模拟切比雪夫II型滤波器时。（3）椭圆滤波器方法一：MATLAB提供函数[z,p,k]=ellipap(N,𝑅𝑝,𝐴𝑠)来设计一个阶数为N，通带波动为𝑅𝑝，阻带衰减为𝐴𝑠的归一化（截止频率为）椭圆模拟原型滤波器，它返回数组z和零极点p以及增益k，将归一化椭圆滤波器系统函数的极点和零点倍乘Ω𝑐，增益倍乘s=0时非归一化的有理函数和归一化的有理函数的比值，即可实现由归一化切比雪夫II型滤波器到非归一化椭圆滤波器的变换。方法二：在MATLAB中还可以直接利用ellip函数设计椭圆模拟原型滤波器。ellip函数既可用于设计模拟椭圆滤波器也可用于设计数字椭圆滤波器，用于设计模拟椭圆滤波器时。设计公式：Ω𝑐=Ω𝑝N=⌈𝐾(𝑘)𝐾(√1−𝑘12)𝐾(𝑘1)𝐾(√1−𝑘2)⌉其中，k=Ω𝑝/Ω𝑠𝑡,𝑘1=𝜀√𝐴2−1,K(x)=∫𝑑𝜃√1−𝑥2𝑠𝑖𝑛2𝜃𝜋/203、模拟滤波器到数字滤波器的变换（1）脉冲响应不变法：数字信号处理实验报告3基本原理：从时域响应出发，使数字滤波器的单位脉冲响应h(n)模仿模拟滤波器的单位冲激响应ℎ𝑎(t),h(n)等于ℎ𝑎(t)的取样值。变换方法：zaaHshthnThnHz拉普拉斯逆变换时域采样变换将数字滤波器性能指标转换成响应的模拟滤波器频率指标Ω𝑝=𝜔𝑝𝑇Ω𝑠𝑡=𝜔𝑠𝑡𝑇将H(s)转换为111kNksTkTAHzez方法一：利用residue函数和residuez函数实现脉冲响应不变法[r,p,k]=residue(b,a),[b,a]=residue(r,p,k)实现H(s)的分解[r,p,k]=residuez(b,a),[b,a]=residuez(r,p,k)实现H(z)的分解方法二：MATLAB提供了impinvar函数采用脉冲响应不变法实现模拟滤波器到数字滤波器的变换，其只用方法如下：[bz,az]=impinvar(b,a,fs)采用脉冲响应不变法将模拟滤波器系统函数的系数向量b和a变换成为数字滤波器系统函数的系数向量bz和az，fs为采样率。[bz,az]=impinvar(b,a)采样频率默认为1的情况下，采用脉冲响应不变法将模拟滤波器变换成为数字滤波器。（2）双线性变换法：基本原理：从频率响应出发，直接使数字滤波器的频率响应逼近模拟滤波器的频率响应，进而求得数字滤波器的系统函数H(z)。变换方法：根据频率的非线性关系，确定预畸的模拟滤波器频率指标Ω𝑝=2𝑇tan(𝜔𝑝2)Ω𝑠𝑡=2𝑇tan(𝜔𝑠𝑡2)将s替换为s=1−𝑧−11+𝑧−1即可得H(z).四、实验内容设采样频率为fs=10kh,设计数字低通滤波器，满足下列指标：通带截止频率fp=1khz,通带波动Rp=1db阻带截止频率fst=1.5khz,阻带波动As=15db要求设计巴特沃斯型、切比雪夫Ⅰ型、切比雪夫Ⅱ型、椭圆滤波器，并分别利用脉冲响应不变变换法和双线性变换法设计。（1）巴特沃斯型程序代码：wp=1000*2*pi;ws=1500*2*pi;Rp=1;As=15;N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(wp/ws)))数字信号处理实验报告4OmegaC=wp/((10^(Rp/10)-1)^(1/(2*N)))[b,a]=butter(N,OmegaC,'s');fs=10000;[bz,az]=impinvar(b,a,fs);[H,w]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');程序运行结果：N=6OmegaC=7.0321e+03MATLAB图形：00.5100.51()|H(ej)|00.51-80-60-40-200()|H(ej)|(dB)00.51-1-0.500.51()PhaseofH(ej)()数字信号处理实验报告5（2）切比雪夫Ⅰ型程序代码：Rp=1;As=15;wp=1000*2*pi;ws=1500*2*pi;E=((10^(Rp/10))-1)^0.5;A=10^(As/20);N=ceil((acosh(((A^2-1)^0.5)/E))/(acosh(ws/wp)))[b,a]=cheby1(N,Rp,wp,'s');fs=10000;[bz,az]=impinvar(b,a,fs);[H,w]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');运行结果：N=4MATLAB图形：数字信号处理实验报告6（3）切比雪夫Ⅱ型程序代码：Rp=1;As=15;wp=1000*2*pi;ws=1500*2*pi;E=((10^(Rp/10))-1)^0.5;A=10^(As/20);N=ceil((acosh(((A^2-1)^0.5)/E))/(acosh(ws/wp)))fs=10000;[b,a]=cheby2(N,As,ws,'s');[bz,az]=impinvar(b,a,fs);[H,w]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|');subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;xlabel('\omega(\pi)');ylabel('|H(e^j^\omega)|(dB)');subplot(223);plot(w/pi,angle(H)/pi);gridon;00.5100.511.5()|H(ej)|00.51-80-60-40-200()|H(ej)|(dB)00.51-1-0.500.51()PhaseofH(ej)()数字信号处理实验报告7xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');运行结果：N=4MATLAB图形：（4）椭圆型程序代码：Rp=1;As=15;wp=1000*2*pi;ws=1500*2*pi;E=((10^(Rp/10))-1)^0.5;A=10^(As/20);k=wp/ws;k1=E/((A^2-1)^0.5);N=ceil(ellipke(k)*ellipke((1-k1^2)^0.5)/(ellipke