06 IIR(3) _ 双线性变换法(New).ppt

数字信号处理,2/68,脉冲响应不变法优点：时域逼近。使数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，即时域逼近良好。线性频率关系：模拟频率和数字频率之间呈线性关系=T。缺点：混叠失真效应。因此，只适用于限带的模拟滤波器（例如衰减特性很好的低通或带通滤波器），而且高频衰减越快，混叠效应越小；而对于高通和带阻滤波器，由于它们在高频部分不衰减，因此会产生混叠现象。,回顾,3/68,例1,【例】拟设计一个数字低通滤波器，指标如下：,【解】代码如下：T=1;%采样周期Wp=0.2*pi/T;Rp=1;Ws=0.35*pi/T;Rs=10;N,Wc=buttord(Wp,Ws,Rp,Rs,s);B,A=butter(N,Wc,s);%求Ha(s)Bz,Az=impinvar(B,A,1/T);%求H(z),4/68,例1,下面的代码分别绘制模拟滤波器和数字滤波器的频响特性：w=-4*pi:4*pi;h=freqs(B,A,w);%求Ha(s)的频响subplot(2,1,1),plot(w,20*log10(abs(h);gridon;xlabel(模拟频率();ylabel(幅度|Ha(j)|);axis(-10,10,-90,10);Hk=freqz(Bz,Az,w);%求H(z)的频响subplot(2,1,2),plot(w,20*log10(abs(Hk);gridon;xlabel(数字频率();ylabel(幅度|H()|);axis(-10,10,-90,10);,5/68,例1,频谱混叠,6/68,采用非线性频率压缩方法：首先，将S平面压缩映射到S1平面的横带区域；其次，通过标准变换关系将此横带区域变换到Z平面。,解决的思路,7/68,解决的思路,双线性变换法的基本思路：从频率响应出发，直接使数字滤波器的频率响应，逼近模拟滤波器的频率响应，进而求得H(z)。,在频域进行逼近,8/68,S平面的虚轴j压缩到S1平面j1轴上的段：,式中,T是采样间隔。,1基本原理,一、S平面S1平面频率压缩,当，则,即：,双线性变换法,9/68,将此关系解析延拓到整个S平面和S1平面。令j=s，j1=s1，则有：,1基本原理,双线性变换法,10/68,S1平面映射到Z平面：,从而得到S平面和Z平面的单值映射关系为：,1基本原理,二、S1平面Z平面数字化,S平面与Z平面之间的单值映射都是线性函数，故称双线性变换。,双线性变换法,11/68,三、性能论证,1基本原理,设：,1.代入，有：,即S平面的虚轴映射到Z平面的单位圆。,双线性变换法,需要验证的结论：S平面的虚轴映射到Z平面的单位圆上；位于S左半平面的极点应映射到Z平面的单位圆内。,12/68,2.其次，将代入，得：,因此,1基本原理,当0时，|z|1：S平面的右半平面映射到Z平面的单位圆外。,双线性变换法,系统函数经映射后稳定性不变,13/68,四、特点S平面整个j轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系：,1基本原理,双线性变换法,14/68,频率之间的非线性变换关系：,2存在的问题,双线性变换法,原来的频率关系：线性关系,15/68,2存在的问题,导致一些缺陷。阅读P185一、线性相位的DF经变换非线性相位的DF,16/68,2存在的问题,双线性变换法,二、DF的幅频响应相对于AD的幅频响应有畸变,例如：一个模拟微分器，它的幅度与频率是直线关系，但通过双线性变换后，就不可能得到数字微分器。,17/68,三、要求模拟滤波器的幅频响应必须是分段常数型即某一频率段的幅频响应近似于某一常数。,2存在的问题,双线性变换法,分段常数型AF经变换后，仍为分段常数型DF。分段边缘的临界频率点从AF转换到DF时，对应关系产生畸变。预畸变。,18/68,2存在的问题,预畸变将模拟滤波器的临界频率事先加以畸变，然后通过双线性变换后正好映射到所需要的频率上。利用关系式：将所要设计的数字滤波器临界频率点，变换成对应的模拟域频率，利用此设计模拟滤波器，再通过双线性变换，即可得到所需的数字滤波器，其临界频率正是。如图所示。,双线性变换法,19/68,2存在的问题,双线性变换法,20/68,由于s到z之间的简单代数关系，可以直接得到数字滤波器的系统函数：,频率响应也可用直接代换的方法得到,3数字化方法,双线性变换法,21/68,具体实现的二种方法：先将Ha(s)分解成并联或级联形式，再分别采用双线性变换。,3数字化方法,双线性变换法,22/68,预先求出H(z)与Ha(s)系数之间的关系式，列成表格，利用表格进行设计。,3数字化方法,双线性变换法,假设:,得：,23/68,3数字化方法,双线性变换法,系数关系表,24/68,例2,设计巴特沃斯低通数字滤波器，其技术指标如下：【解】：考虑预畸变，得模拟滤波器的技术指标：这里，假设T=1s。,25/68,例2,（1）确定参数,取：,26/68,例2,（2）求出极点（左半平面）,（3）构造系统函数,27/68,例2,或者，由N=6，直接查表，得：,再去归一化，得：,28/68,例2,（4）将变换成数字滤波器：,29/68,MATLAB函数,bilinear函数：实现双线性变换,阅读P189,30/68,例1（续）,代码如下：T=1;%采样周期Fs=1/T;Wp=2*tan(0.2*pi/2)/T;Rp=1;Ws=2*tan(0.35*pi/2)/T;Rs=10;N,Wc=buttord(Wp,Ws,Rp,Rs,s);B,A=butter(N,Wc,s);%求Ha(s)Bz,Az=bilinear(B,A,Fs);%求H(z),再看前面的例1，用双线性变换法设计一个数字低通滤波器，指标如下：,31/68,例1（续）,w=-pi:pi;h=freqs(B,A,w);subplot(2,1,1),plot(w,20*log10(abs(h);gridon;xlabel(模拟频率();ylabel(幅度|Ha(j)|);axis(-4,4,-50,4);Hk=freqz(Bz,Az,w);subplot(2,1,2),plot(w,20*log10(abs(Hk);gridon;xlabel(数字频率();ylabel(幅度|H()|);axis(-4,4,-50,10);,32/68,例1（续）,没有频谱混叠,33/68,本单元小结,设计模拟滤波器转换为数字滤波器脉冲响应不变法双线性变换法,34/68,本单元小结,设计模拟滤波器已知：Ap、As、s、p要求：确定滤波器阶次N和截止频率c。,35/68,本单元小结,MatLab提供的函数：N,Wc=buttord(Wp,Ws,Rp,Rs,s)B,A=butter(N,Wc,ftype,s),36/68,转换为数字滤波器,本单元小结,与的关系=T,37/68,将Ha(s)直接转换为数字滤波器H(z)。,本单元小结,脉冲响应不变法,极点传递时域逼近,38/68,脉冲响应不变法的优点：时域逼近。使数字滤波器的单位脉冲响应完全模仿模拟滤波器的单位冲激响应，即时域逼近良好。线性频率关系。模拟频率和数字频率之间呈线性关系=T。脉冲响应不变法的缺点：混叠失真效应因此，只适用于限带的模拟滤波器（例如衰减特性很好的低通或带通滤波器），而且高频衰减越快，混叠效应越小；而对于高通和带阻滤波器，由于它们在高频部分不衰减，因此会产生混叠现象。,本单元小结,39/68,本单元小结,频谱混叠,40/68,本单元小结,双线性变换法从频率响应出发，直接使数字滤波器的频率响应，逼近模拟滤波器的频率响应，进而求得H(z)。,频域逼近,S平面与Z平面之间的单值映射都是线性函数，故称双线性变换。,41/68,频率之间的非线性变换关系：,本单元小结,原来的频率关系：线性关系,42/68,MATLAB函数,Bz,Az=impinvar(b,a,Fs);Bz,Az=bilinear(b,a,Fs);,43/68,例3,已知某信号如下：请设计一个巴特沃斯数字低通滤波器，指标如下：采样频率：Fs=256Hz通带截止频率：fp=60Hz阻带截止频率：fs=70Hz通带最大衰减：Ap=1dB阻带最小衰减：As=25dB绘制滤波后信号的波形与频谱。,44/68,例3原始信号,clearall;Adc=2;%直流分量A1=3;%信号1的幅度A2=1.5;%信号2的幅度F1=50;%信号1的频率(Hz)F2=75;%信号2的频率(Hz)Fs=256;%采样频率(Hz)P1=-30;%信号1的相位(度)P2=90;%信号2的相位(度)N=256;%采样点数t=0:1/Fs:N/Fs;%采样时刻S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);figure(1);subplot(2,1,1);plot(S);title(滤波前的信号);Y=fft(S,N);%FFT变换Ayy=(abs(Y);%取模Ayy=Ayy/(N/2);%换算成实际幅度:An=A/(N/2),45/68,例3原始信号,Ayy(1)=Ayy(1)/2;%换算第1个点模值：A0=A/NF=(1:N-1)*Fs/N;%换算成实际频率值：F=(n-1)*fs/Nsubplot(2,1,2);plot(F(1:N/2),Ayy(1:N/2);%显示换算后结果,仅显示左半部分title(滤波前信号的频谱);holdon%标注关键点的坐标;fori=1:N/2ifAyy(i)0.01%实数一般不可能精确为0plot(F(i),Ayy(i),-o,MarkerFaceColor,r);%画实心圆点str=X:,num2str(F(i),Y:,num2str(Ayy(i);text(F(i),Ayy(i),str);%显示坐标endendholdoff,46/68,例3脉冲响应不变法,Fs=256;%采样频率fp=60;%通带截止频率fs=70;%阻带截止频率Rp=1;Rs=25;Wp=(fp/Fs)*2*pi;%临界频率采用角频率表示Ws=(fs/Fs)*2*pi;%临界频率采用角频率表示OmegaP=Wp*Fs;OmegaS=Ws*Fs;n,Wc=buttord(OmegaP,OmegaS,Rp,Rs,s);b,a=butter(n,Wc,s);Bz,Az=impinvar(b,a,Fs);,47/68,例3绘制频响图,freqz函数格式：H,W=freqz(B,A,N)在半圆（0）上取N个频率等分点H,W=freqz(B,A,N,whole)在单位圆（02）上取N个频率等分点其中：B、A：分别为系统函数分子、分母多项式的系数向量；N：默认值是512；W：包含了范围内的N个频率等分点（频率轴上）；H：包含了离散系统频响在N个频率等分点上的值。,48/68,例3绘制频响图,freqz另二种格式：H,F=freqz(B,A,N,Fs)H,F=freqz(B,A,N,whole,Fs)其中：Fs：采样频率（Hz）F：返回的频率向量（Hz）,不同版本稍有差异，使用时注意区分。,49/68,例3绘制频响图,figure(2);Hk,W=freqz(Bz,Az,N);%求H(z)的频响subplot(2,1,1);plot(W/pi,20*log10(abs(Hk);%对pi归一化xlabel(角频率();ylabel(幅度);title(滤波器频响);Hk,F=freqz(Bz,Az,N,Fs);%求H(z)的频响subplot(2,1,2);plot(F,20*log10(abs(Hk);xlabel(频率(Hz);ylabel(幅度);,50/68,例3脉冲响应不变法频响,51/68,例3双线性变换法,Fs=256;%采样频率fp=60;%通带截止频率fs=70;%阻带截止频率Rp=1;Rs=25;Wp=(fp/Fs)*2*pi;%临界频率采用角频率表示Ws=(fs/Fs)*2*pi;%临界频率采用角频率表示OmegaP=2*Fs*tan(Wp/2);%频率预畸OmegaS=2*Fs*tan(Ws/2);n,Wc=buttord(OmegaP,OmegaS,Rp,Rs,s);b,a=butter(n,Wc,s);Bz,Az=bilinear(b,a,Fs);,52/68,例3双线性变换法频响,53/68,例3脉冲响应不变法滤波效果,y=filter(Bz,Az,S);%进行滤波figure(3);subplot(2,1,2);plot(y);%滤波后信号的时域波形title(滤波后的信号);Y=fft(y,N);Ayy=(abs(Y);Ayy=Ayy/(N/2);Ayy(1)=Ayy(1)/2;F=(1:N-1)*Fs/N;subplot(2,1,2);plot(F(1:N/2),Ayy(1:N/2);%滤波后信号的频谱title(滤波后信号的频谱);,54/68,例3双线性变换法滤波效果,55/68,练习1高通,目标：提取第2个余弦信号：拟设计的高通滤波器的指标如下：通带截止频率：fp=70Hz阻带截止频率：fs=60Hz通带最大衰减：Ap=1dB阻带最小衰减：As=25dB,56/68,练习1高通,fp=70;%通带截止频率fs=60;%阻带截止频率Rp=1;Rs=25;Wp=(fp/Fs)*2*pi;Ws=(fs/Fs)*2*pi;OmegaP=2*Fs*tan(Wp/2);OmegaS=2*Fs*tan(Ws/2);n,Wc=buttord(OmegaP,OmegaS,Rp,Rs,s);b,a=butter(n,Wc,high,s);Bz,Az=bilinear(b,a,Fs);%一般用双线性变换法,57/68,练习2带通,目标：提取第1个余弦信号：拟设计的带通滤波器的指标如下：通带截止频率：fp1=40Hz，fp2=60Hz阻带截止频率：fs1=30Hz，fs2=70Hz通带最大衰减：Ap=1dB阻带最小衰减：As=25dB,58/68,练习2带通,fp1=40;%通带截止频率fs1=30;%阻带截止频率fp2=60;%通带截止频率fs2=70;%阻带截止频率Rp=1;Rs=25;Wp1=(fp1/Fs)*2*pi;Ws1=(fs1/Fs)*2*pi;Wp2=(fp2/Fs)*2*pi;Ws2=(fs2/Fs)*2*pi;Wp=Wp1,Wp2;%向量Ws=Ws1,Ws2;%向量,59/68,练习2带通,OmegaP=2*Fs*tan(Wp/2);OmegaS=2*Fs*tan(Ws/2);n,Wc=buttord(OmegaP,OmegaS,Rp,Rs,s);b,a=butter(2*n,Wc,s);%对于二维向量，默认是带通Bz,Az=bilinear(b,a,Fs);频响图为：,60/68,练习2带通（滤波效果）,61/68,练习3带阻,目标：滤除第1个余弦信号：拟设计的带阻滤波器的指标如下：通带截止频率：fp1=30Hz，fp2=70Hz阻带截止频率：fs1=40Hz，fs2=60Hz通带最大衰减：Ap=1dB阻带最小衰减：As=25dB,62/68,练习3带阻,fp1=30;%通带截止频率fs1=40;%阻带截止频率fp2=70;%通带截止频率fs2=60;%阻带截止频率Rp=1;Rs=25;Wp1=(fp1/Fs)*2*pi;Ws1=(fs1/Fs)*2*pi;Wp2=(fp2/Fs)*2*pi;Ws2=(fs2/Fs)*2*pi;Wp=Wp1,Wp2;Ws=Ws1,Ws2;,63/68,练习3带阻,OmegaP=2*Fs*tan(Wp/2);OmegaS=2*Fs*tan(Ws/2);n,Wc=buttord(OmegaP,OmegaS,Rp,Rs,s);b,a=butter(2*n,Wc,stop,s);%