双线性变换法设计IIR数字滤波器-matlab.doc
Four short words sum up what has lifted most successful individuals above the crowd: a little bit more.-author-date双线性变换法设计IIR数字滤波器-matlab双线性变换法设计IIR数字滤波器-matlab实验四 双线性变换法设计IIR数字滤波器一、 实验目的1. 掌握利用双线性变换法设计IIR数字滤波器的原理和具体方法。2. 深入理解利用双线性变换法设计IIR数字滤波器的优缺点及使用范围。二、 实验内容1. 利用巴特沃思模拟滤波器,通过双线性变换法设计巴特沃思数字滤波器,数字滤波器的技术指标为采样周期为T=1。程序代码T=2; %设置采样周期为2fs=1/T; %采样频率为周期倒数Wp=0.25*pi/T; Ws=0.35*pi/T; %设置归一化通带和阻带截止频率Ap=20*log10(1/0.9);As=20*log10(1/0.18); %设置通带最大和最小衰减 N,Wc=buttord(Wp,Ws,Ap,As,'s'); %调用butter函数确定巴特沃斯滤波器阶数B,A=butter(N,Wc,'s'); %调用butter函数设计巴特沃斯滤波器W=linspace(0,pi,400*pi); %指定一段频率值 hf=freqs(B,A,W); %计算模拟滤波器的幅频响应 subplot(2,1,1);plot(W/pi,abs(hf)/abs(hf(1); %绘出巴特沃斯模拟滤波器的幅频特性曲线 grid on;title('巴特沃斯模拟滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');D,C= bilinear (B,A,fs); %调用双线性变换法 Hz=freqz(D,C,W); %返回频率响应 subplot(2,1,2);plot(W/pi,abs(Hz)/abs(Hz(1); %绘出巴特沃斯数字低通滤波器的幅频特性曲线 grid on;title('巴特沃斯数字滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');进行试验输出图像：实验分析通过查workspace可以得出脉冲响应为：滤波器幅度响应：滤波器相位响应：2. 利用巴特沃思模拟滤波器，通过双线性变换法设计数字带阻滤波器，数字滤波器的技术指标为采样周期为T=1。程序代码T=1; %设置采样周期为1fs=1/T; %采样频率为周期倒数wp=0.25*pi,0.75*pi; ws=0.35*pi,0.65*pi; Wp=(2/T)*tan(wp/2); Ws=(2/T)*tan(ws/2); %设置归一化通带和阻带截止频率Ap=20*log10(1/0.9);As=20*log10(1/0.18); %设置通带最大和最小衰减 N,Wc=buttord(Wp,Ws,Ap,As,'s');%调用butter函数确定巴特沃斯滤波器阶数B,A=butter(N,Wc, 'stop','s'); %调用butter函数设计巴特沃斯滤波器W=linspace(0,2*pi,400*pi); %指定一段频率值 hf=freqs(B,A,W); %计算模拟滤波器的幅频响应 subplot(2,1,1);plot(W/pi,abs(hf); %绘出巴特沃斯模拟滤波器的幅频特性曲线 grid on;title('巴特沃斯模拟滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');D,C=bilinear(B,A,fs); %调用双线性变换法 Hz=freqz(D,C,W); %返回频率响应 subplot(2,1,2);plot(W/pi,abs(Hz); %绘出巴特沃斯数字带阻滤波器的幅频特性曲线grid on;title('巴特沃斯数字滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');进行试验输出图像：实验分析通过查workspace可以得出脉冲响应为：滤波器幅度响应：滤波器相位响应：3. 通过切比雪夫模拟滤波器，利用双线性变换法设计数字高通滤波器，数字滤波器的技术指标为采样周期为T=1，并观察T分别为0.01，0.1，0.25，0.5，0.8时数字滤波器的频率响应。程序代码T=1; %设置采样周期为1fs=1/T; %采样频率为周期倒数 Wp=(2/T)*tan(0.75*pi/2); Ws=(2/T)*tan(0.65*pi/2); %设置归一化通带和阻带截止频率Ap=20*log10(1/0.9);As=20*log10(1/0.18); %设置通带最大和最小衰减 N,Wc=cheb1ord(Wp,Ws,Ap,As,'s');%调用cheb1ord函数确定切比雪夫I型滤波器阶数B,A=cheby1(N,Ap,Wc, 'high','s'); %调用cheby1函数设计切比雪夫I型滤波器W=linspace(0,7*pi,400*pi); %指定一段频率值 hf=freqs(B,A,W); %计算模拟滤波器的幅频响应 subplot(2,1,1);plot(W/pi,abs(hf); %绘出切比雪夫I型模拟滤波器的幅频特性曲线 grid on;title('切比雪夫I型模拟滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');w=linspace(0,pi,400*pi); %指定一段频率值D,C=bilinear(B,A,fs); %调用双线性变换法 Hz=freqz(D,C,w); %返回频率响应 subplot(2,1,2);plot(w/pi,abs(Hz); %绘出切比雪夫I型数字带阻滤波器的幅频特性曲线grid on;title('切比雪夫I型数字滤波器');xlabel('Frequency/Hz');ylabel('Magnitude');进行试验输出图像如下：实验分析通过查workspace可以得出脉冲响应为：滤波器幅度响应：滤波器相位响应：当T分别为0.01，0.1，0.25，0.5，0.8时的数字滤波器的频率响应如下图像所示：其中，T=0.01“ ”T=0. 1“-”T=0.25“··”T=0.5“”T=0.8“ ”由图像可以看出，当采样周期从0.01s增长到0.8s时，数字滤波器的幅度响应没有任何变化。故对于用双线性变换法设计的切比雪夫高通数字滤波器，由于事先做了预畸处理，消除了频率映射时引起的非线性关系，所得数字滤波器的幅度响应较为稳定。对比用双线性变换法和脉冲响应不变法得到的数字高通滤波器幅度响应形状可知：双线性变换法不会出现由于高频部分超过折叠频率而混叠到低频部分去的现象。这意味着，模拟滤波器的全部频率响应被压缩于等效的数字频率范围（0，）之内。这是双线性变换法的主要优点，虽然这是以频率的严重非线性关系为代价的，但因为可以用预畸的方法进行补偿，所以仍可以得到比较理想的数字滤波器幅度响应各种设计方法的MATLAB实现 下面以一带通数字滤波器的设计为例,说明基于MATLAB的三种实现方法。3.1 基于模拟低通原型的MATLAB实现%通过模拟低通原型滤波器进行数字带通设计程序：fp=480 %模拟低通通带上限频率fs=520 %模拟低通阻带下限频率wp=2*pi*fp %模拟低通通带上限角频率ws=2*pi*fs %模拟低通阻带下限角频率rp=3 %通带波动rs=20 %阻带衰减%巴特沃斯模拟低通原型滤波器设计n,wn=buttord(wp,WS,rp ,rs,'s')z,p,k=buttap(n) %模拟低通原型零、极点系数b1,a1=zp2tf(z,p,k) %零、极点系数转换为传递函数%巴特沃斯模拟低通原型滤波器频率响应hl,w1=freqs(b1,a1)mag1=abs(h1)%模拟低通原型滤波器幅频特性曲线subplot(221);semilogx(w1,mag1)fw=40 %模拟带通滤波器带宽频率bw=2*pi*fw %模拟带通滤波器带宽角频率% 由模拟低通原型变换为模拟带通滤波器b2,a2=lp2bp(b1,a1,wn,bw) %模拟带通滤波函数系数%巴特沃斯模拟带通滤波器频率响应h2,w2= freqs(b2,a2)mag2 =abs(1l2)%模拟带通滤波器幅频特性曲线(db)subplot(222);plot(20*log10(mag2)% 冲击响应不变法进行离散化设计fo=2000 %采样频率bz,az=impinvar(b2,a2,2000) %数字带通滤波函数系数%巴特沃斯型数字带通滤波器频率响应hz,w= freqz(bz,az)magz=abs(hz)phz=unwrap(angle(hz)subplot(223);plot(magz) %数字带通滤波器幅频特性曲线subplot(224);flot(plot) %数字带通滤波器相频特性曲线 3.2 基于合适类型模拟滤波器的MATLAB实现%通过合适类型模拟滤波器进行数字带通设计程序fp= 480,520;fs=450,550 %模拟通带、阻带频率wp=480,520*pi*2 %模拟通带角频率ws=450,550*pi*2 %模拟阻带角频率rp=3;rs=20 %通带波动、阻带衰减% 巴特沃斯型模拟带通滤波器设计n,wn=buttord (wp,ws,rp,rs,'s')b,a=butter(n,wn,'s') %模拟带通滤波函数系数%巴特沃斯型模拟带通滤波器频率响应ha,w= freqs(b,a)ma=abs(ha);pha=unwrap(angle(ha)subplot(421);plot(w/(2*pi),ma) %模拟幅频曲线subplot(423);plot(w/(2 pi),pha) %模拟相频曲线% 冲击响应不变法进行离散化设计fo=5000 %采样频率bn,an=impinvar(b,a,5000) %数字带通滤波函数系数%巴特沃斯型数字带通滤波器频率响应hz,w=freqz(bn,an)mz=abs(hz);phz=unwrap(angle(hz)subplot(422);plot(w,mz) %数字滤波器幅频曲线subplot(424);plot(w,phz) %数字滤波器相频曲线hi=impz(bn,an) %数字滤波器冲击响应subplot(425),plot(hi) %冲击响应曲线n=0:300;t=n/foxl=2*square(2*pi*500*t) %500Hz方波信号subplot(426);plot(x1) %500Hz方波波形yi=conv(hi,x1) %时域卷积输出subplot(427);plot(yi) %卷积输出波形y1=filter(bn,an,x1) %数字滤波函数输出subplot(428);plot(y1) %数字滤波器输出波形 3.3 基于直接原型变换法的MATLAB实现%数字带通滤波器直接设计程序fp= 480,520;fs=450,550 %模拟通带、阻带频率rp=3;rs=20 %通带波动、阻带衰减fo=10000 %采样频率%频率指标变换wp=2*pi*fp/f0 %数字通带频率ws=2*pi*fs/fo %数字阻带频率%切比雪夫1型数字带通滤波器直接设计n,wn=cbeblord(wp/pi,ws/pi,rp,rs)b,a=chebyl(n,rp,wn) %数字带通滤波器系数%切比雪夫1型数字带通滤波器频率响应h,w= freqz(b,a,128,10000)mag=abs(h;pha=unwrqp(angle(h)subplot(321);plot(w,mag) %幅频曲线subplot(322);plot(w,pha) %相频曲线hi=impz(b,a) %冲击响应subplot(324);plot(hi) %响应曲线n=0:500;t=n/fcx1=2*square(2*pi*500*t) %500Hz方波信号subplot(323);plot(t,x1) %500Hz方波波形yi=conv(hi,x1) %时域卷积输出subplot(326);plot(yi) %卷积输出波形y1=filter(b,a,x1) %数字滤波函数输出subplot(325); stem(y1) %数字滤波器输出波形 4 结束语 基于MATLAB的信号处理工具箱为数字滤波器设计带来了全新的实现手段,设计快捷方便,仿真波形直观。上述三种设计方案均可实现设计指标,但以直接原型变换法最为简便。实际应用中,数字滤波器也可以对连续时间信号进行处理,但需要先对连续信号进行A/D变换,经数字滤波后,再经D/A转换得到所需要的连续信号(此处略)。参考文献1 倪养华,等数字信号处理与实现上海：上海交通大学出版社,19982 施阳等MATLAB语言工具箱西安：西北工业大学出版-
