ch44雙線性變換法北京交通大學(xué)陳后金教授信號(hào)處理課件

1、近代數(shù)字信號(hào)處理(Advanced Digital Signal Processing),電子信息工程學(xué)院,信號(hào)與圖像處理研究室,IIR數(shù)字濾波器設(shè)計(jì)的基本思想 模擬低通濾波器設(shè)計(jì) 模擬域頻率變換 脈沖響應(yīng)不變法 雙線性變換法 IIR數(shù)字濾波器的基本結(jié)構(gòu) 利用MATLAB設(shè)計(jì)IIR DF,IIR數(shù)字濾波器的設(shè)計(jì),問(wèn)題的提出 雙線性變換法的基本原理 雙線性變換法設(shè)計(jì)DF的步驟,雙線性變換法,問(wèn)題的提出,如何將模擬濾波器轉(zhuǎn)變?yōu)閿?shù)字濾波器?,1. 脈沖響應(yīng)不變法,2. 雙線性變換法,問(wèn)題的提出,采用脈沖響應(yīng)不變法,DF的頻譜有混疊,解決方法,采用雙線性變換法,上節(jié)例題利用AF-BW filter及脈

2、沖響應(yīng)不變法設(shè)計(jì)一DF，滿足 Wp=0.2p, Ws=0.6p, Ap2dB, As15dB 。,As = 14.2dB,雙線性變換法的基本原理,基本思想：將非帶限的模擬濾波器映射為最高頻率為的帶限模擬濾波器,模擬頻率與數(shù)字頻率的關(guān)系為,s域到z域的映射關(guān)系,雙線性變換法的基本原理,雙線性變換,雙線性變換法的基本原理,穩(wěn)定性分析,令s=+j，則有,雙線性變換法的基本原理,穩(wěn)定性分析,1) s0,|z|1,S域左半平面映射到z域單位元內(nèi),2) s=0,|z|=1,3) s0,|z|1,S域虛軸映射到z域單位圓上,S域右半平面映射到z域單位圓外,因果、穩(wěn)定的AF系統(tǒng)映射為因果、穩(wěn)定的DF系統(tǒng),雙線

3、性變換法的基本原理,W和w 的關(guān)系,雙線性變換法的基本原理,雙線性變換法的優(yōu)缺點(diǎn),缺點(diǎn)：幅度響應(yīng)不是常數(shù)時(shí)會(huì)產(chǎn)生幅度失真,優(yōu)點(diǎn)：無(wú)混疊,雙線性變換法設(shè)計(jì)DF的步驟,1. 將數(shù)字濾波器的頻率指標(biāo)Wk轉(zhuǎn)換為 模擬濾波器的頻率指標(biāo)wk,2. 由模擬濾波器的指標(biāo)設(shè)計(jì)模擬濾波器的H(s)。,3. 利用雙線性變換法，將H(s)轉(zhuǎn)換H(z)。,雙線性變換法設(shè)計(jì)DF的步驟,Wp,Ws,wp,ws,H(s),H(z),設(shè)計(jì)模擬 濾波器,雙線性變換,numd,dend = bilinear(num,den,Fs),num,den：AF分子、分母多項(xiàng)式的系數(shù)向量,Fs：抽樣頻率,numd,dend：DF分子、分母多

4、項(xiàng)式的系數(shù)向量,利用MATLAB,例: 利用BW型模擬低通濾波器和雙線性變換法設(shè)計(jì)滿足指標(biāo)Wp=p/3，Ap=3dB，N=1的數(shù)字低通濾波器，并與脈沖響應(yīng)不變法設(shè)計(jì)的DF比較。,解：設(shè)雙線性變換中的參數(shù)為T,(1) 將DF的頻率指標(biāo)轉(zhuǎn)換為AF的頻率指標(biāo),(2) 設(shè)計(jì)3dB截頻為wp的一階BW型模擬低通濾波器，即,N=1, wc = wp,故,例: 利用BW型模擬低通濾波器和雙線性變換法設(shè)計(jì)滿足指標(biāo)Wp=p/3，Ap=3dB，N=1的數(shù)字低通濾波器，并與脈沖響應(yīng)不變法設(shè)計(jì)的DF比較。,解：設(shè)雙線性變換中的參數(shù)為T,(3) 用雙線性變換法將模擬濾波器轉(zhuǎn)換為數(shù)字濾波器,結(jié)論：參數(shù)T的取值和最終的設(shè)計(jì)

5、結(jié)果無(wú)關(guān)。,為簡(jiǎn)單起見(jiàn)，一般取T=2,例: 利用BW型模擬低通濾波器和雙線性變換法設(shè)計(jì)滿足指標(biāo)Wp=p/3，Ap=3dB，N=1的數(shù)字低通濾波器，并與脈沖響應(yīng)不變法設(shè)計(jì)的DF比較。,解：,雙線性變換法設(shè)計(jì)的DF的系統(tǒng)函數(shù)為,脈沖響應(yīng)不變法設(shè)計(jì)的DF的系統(tǒng)函數(shù)為,令z=ejW ，可分別獲得兩者的幅度響應(yīng)。,例: 利用BW型模擬低通濾波器和雙線性變換法設(shè)計(jì)滿足指標(biāo)Wp=p/3，Ap=3dB，N=1的數(shù)字低通濾波器，并與脈沖響應(yīng)不變法設(shè)計(jì)的DF比較。,脈沖響應(yīng)不變法,雙線性變換法,脈沖響應(yīng)不變法存在頻譜混疊，所設(shè)計(jì)的DF不滿足給定指標(biāo)。而雙線性變換法不存在頻譜混疊，所設(shè)計(jì)的DF滿足給定指標(biāo)。,3dB

6、,例: 利用BW型模擬低通濾波器和雙線性變換法設(shè)計(jì)滿足指標(biāo)Wp=p/3，Ap=3dB，N=1的數(shù)字低通濾波器，并與脈沖響應(yīng)不變法設(shè)計(jì)的DF比較。,H雙(z)和H脈(z)幅度響應(yīng)比較的MATLAB實(shí)現(xiàn),Wp=pi/3; b=1-exp(-Wp);b1=tan(Wp/2)*1 1; a=1 -exp(-Wp);a1=1+tan(Wp/2) tan(Wp/2)-1; w=linspace(0,pi,512); h=freqz(b,a,w);h1=freqz(b1,a1,w); plot(w/pi,(abs(h),w/pi,(abs(h1) ); xlabel(Normalized frequency

7、); ylabel(Amplitude); set(gca,ytick,0 0.7 1); set(gca,xtick,0 Wp/pi 1); grid;,例：利用AF-BW filter及雙線性變換法設(shè)計(jì)一DF，滿足 Wp=0.2p, Ws=0.6p, Ap2dB, As15dB,解：,(1) 將數(shù)字低通指標(biāo)轉(zhuǎn)換成模擬低通指標(biāo)，取T=2,Ap2db, As15db,(2) 設(shè)計(jì)模擬低通濾波器 （BW型）,=2,=0.585 1,例：利用AF-BW filter及雙線性變換法設(shè)計(jì)一DF，滿足 Wp=0.2p, Ws=0.6p, Ap2dB, As15dB,解：,(3) 用雙線性變換法將模擬低通

8、濾波器轉(zhuǎn)換成數(shù)字低通濾波器,例：利用AF-BW filter及雙線性變換法設(shè)計(jì)一DF，滿足 Wp=0.2p, Ws=0.6p, Ap2dB, As15dB,%Design DF BW low-pass filter using impulse invariance %DF BW LP specfication Wp=0.2*pi; Ws=0.6*pi; Ap=2; As=15; T=2;Fs=1/T; %Sampling frequency(Hz) %Analog Butterworth specfication wp=2*tan(Wp/2)/T;ws=2*tan(Ws/2)/T; %dete

9、rmine the order of AF filter and the 3-dB cutoff frequency N,wc=buttord(wp,ws,Ap,As,s) %determine the AF-BW filter numa,dena=butter(N,wc,s),例：利用AF-BW filter及雙線性變換法設(shè)計(jì)一DF，滿足 Wp=0.2p, Ws=0.6p, Ap2dB, As15dB,%determine the DF filter numd,dend=bilinear(numa,dena,Fs) %plot the frequency response w=linspac

10、e(0,pi,1024); h=freqz(numd,dend,w); plot(w/pi,20*log10(abs(h); axis(0 1 -50 0);grid; xlabel(Normalized frequency); ylabel(Gain,dB); %computer Ap As of the designed filter w=Wp Ws;h=freqz(numd,dend,w); fprintf(Ap= %.4fn,-20*log10( abs(h(1); fprintf(As= %.4fn,-20*log10( abs(h(2);,例：利用AF-BW filter及雙線性變

11、換法設(shè)計(jì)一DF，滿足 Wp=0.2p, Ws=0.6p, Ap2dB, As15dB,Ap= 0.3945 As= 15.0000,例：利用AF-BW filter及雙線性變換法設(shè)計(jì)一DF，滿足 Wp=0.2p, Ws=0.6p, Ap2dB, As15dB,將雙線性變換法與脈沖響應(yīng)不變法所設(shè)計(jì)DF的結(jié)果比較。,雙線性變換 Ap= 0.3945 As= 15.0000,脈沖響應(yīng)不變法,雙線性變換法,脈沖響應(yīng)不變 Ap= 1.1187 As= 12.3628,非低通IIR數(shù)字濾波器的設(shè)計(jì),Wp,Ws,wp,ws,H(z),w =W/T,脈沖響應(yīng)不變法,雙線性變換法,H(s),模擬頻 率變換,設(shè)計(jì)

12、原型 低通濾波器,復(fù)頻率變換,注意： 脈沖響應(yīng)不變法不能設(shè)計(jì)高通和帶阻數(shù)字濾波器,方法一,非低通IIR數(shù)字濾波器的設(shè)計(jì),方法二,Wp,Ws,H(z),w =W/T,脈沖響應(yīng)不變法,雙線性變換法,數(shù)字頻 率變換,設(shè)計(jì)原型 低通濾波器,z域 變換,例：試設(shè)計(jì)滿足下列指標(biāo)的BW型數(shù)字帶阻濾波器 Wp1=2.8113rad/s, Wp2=2.9880rad/s, Ap1dB , Ws1=2.9203rad/s, Ws2=2.9603rad/s, As 10dB 。,解:,脈沖響應(yīng)不變法不適合設(shè)計(jì)數(shù)字帶阻濾波器，因此采用雙線性變換法設(shè)計(jì)。,(1) 將數(shù)字帶阻濾波器指標(biāo)轉(zhuǎn)換成模擬帶阻濾波器指標(biāo),取T=2，

13、利用,得模擬帶阻指標(biāo)為,wp1=6rad, wp2=13rad, ws1=9rad, ws2=1rad1,Ap1dB, As 10dB,例：試設(shè)計(jì)滿足下列指標(biāo)的BW型數(shù)字帶阻濾波器 Wp1=2.8113rad/s, Wp2=2.9880rad/s, Ap1dB , Ws1=2.9203rad/s, Ws2=2.9603rad/s, As 10dB 。,解:,(2) 將模擬帶阻濾波器指標(biāo)轉(zhuǎn)換成模擬低通濾波器指標(biāo),Ap1dB, As 10dB,wp1=6rad, wp2=13rad, ws1=9rad, ws2=1rad1,Ap1dB, As 10dB,模擬帶阻指標(biāo),例：試設(shè)計(jì)滿足下列指標(biāo)的BW型

14、數(shù)字帶阻濾波器 Wp1=2.8113rad/s, Wp2=2.9880rad/s, Ap1dB , Ws1=2.9203rad/s, Ws2=2.9603rad/s, As 10dB 。,解:,(3) 設(shè)計(jì)原型BW型模擬低通濾波器,Ap1dB, As 10dB,原型模擬低通指標(biāo),原型模擬低通濾波器的系統(tǒng)函數(shù)為,例：試設(shè)計(jì)滿足下列指標(biāo)的BW型數(shù)字帶阻濾波器 Wp1=2.8113rad/s, Wp2=2.9880rad/s, Ap1dB , Ws1=2.9203rad/s, Ws2=2.9603rad/s, As 10dB,解:,(4) 由復(fù)頻率變換將原型模擬低通轉(zhuǎn)換為模擬帶阻濾波器,(5) 由雙

15、線性變換模擬帶阻濾波器轉(zhuǎn)換成數(shù)字帶阻濾波器,利用MATLAB實(shí)現(xiàn)IIR數(shù)字濾波器,確定數(shù)字濾波器的階數(shù)及3dB截頻Wc N, Wc = buttord(Wp, Ws, Ap, As) 其中Wp, Ws為歸一化角頻率。 例Wp=0.1p, 則Wp=0.1 若為帶通或帶阻濾波器，則 Wp=Wp1, Wp2; Ws=Ws1, Ws2,BW型數(shù)字濾波器,利用MATLAB實(shí)現(xiàn)IIR數(shù)字濾波器,BW型數(shù)字濾波器,確定DF系統(tǒng)函數(shù)分子、分母多項(xiàng)式 低通 num,den = butter(N,Wc) 高通 num,den = butter(N,Wc,high) 帶通 num,den = butter(N,Wc

16、) 其中Wc=W1, W2 帶阻 num,den = butter(N,Wc,stop) 其中Wc=W1, W2,利用MATLAB實(shí)現(xiàn)IIR數(shù)字濾波器,CB I型數(shù)字濾波器,確定數(shù)字濾波器的階數(shù)及參數(shù)Wc N, Wc = cheb1ord(Wp, Ws, Ap, As) DF系統(tǒng)函數(shù)分子、分母多項(xiàng)式的確定 低通 num,den = cheby1(N,Ap,Wc) 高通 num,den = cheby1(N, Ap, Wc,high) 帶通 num,den = cheby1(N, Ap, Wc) 帶阻 num,den = cheby1(N, Ap, Wc,stop) * 帶通、帶阻中W=W1,

17、W2,利用MATLAB實(shí)現(xiàn)IIR數(shù)字濾波器,CB II型數(shù)字濾波器,確定數(shù)字濾波器的階數(shù)及參數(shù)Wc N, Wc = cheb2ord(Wp, Ws, Ap, As) DF系統(tǒng)函數(shù)分子、分母多項(xiàng)式的確定 低通 num,den = cheby2(N,As,Wc) 高通 num,den = cheby2(N, As, Wc,high) 帶通 num,den = cheby2(N, As, Wc) 帶阻 num,den = cheby2(N, As, Wc,stop) * 帶通、帶阻中W=W1, W2,利用MATLAB實(shí)現(xiàn)IIR數(shù)字濾波器,橢圓(C)型數(shù)字濾波器,確定數(shù)字濾波器的階數(shù)及參數(shù)Wc N,

18、Wc = ellipord(Wp, Ws, Ap, As) DF系統(tǒng)函數(shù)分子、分母多項(xiàng)式的確定 低通 num,den = ellip(N,Ap,As,Wc) 高通 num,den = ellip(N,Ap,As,Wc,high) 帶通 num,den = ellip(N,Ap,As,Wc) 帶阻 num,den = ellip(N,Ap,As,Wc,stop) * 帶通、帶阻中W=W1, W2,MATLAB中IIR設(shè)計(jì)的主要函數(shù),例：利用MATLAB實(shí)現(xiàn)數(shù)字帶阻濾波器 Wp1=2.8113rad/s, Wp2=2.9880rad/s, Ap1dB , Ws1=2.9203rad/s, Ws2=

19、2.9603rad/s, As 10dB。,Wp=2.813,2.9880;Ws=2.9203,2.9603; Ap=1;As=10; N,Wc=buttord(Wp/pi,Ws/pi,Ap,As) numd,dend=butter(N,Wc,stop),運(yùn)行結(jié)果 N =2 numd =0.9522, 3. 7327, 5.5624, 3.7327, 0.9522 dend =1.0000, 3.8242 , 5.5601, 3.6412, 0.9067,利用IIR數(shù)字濾波器處理脈搏信號(hào), 中醫(yī)或西醫(yī)領(lǐng)域的多種疾病診斷分析（心血管疾病等） 疲勞狀態(tài)預(yù)測(cè)（司機(jī)疲勞駕駛等） 情緒變化分析（學(xué)生、孕

20、婦情緒監(jiān)測(cè)） 心理學(xué)分析（測(cè)謊、犯罪等心理研究）,采集脈搏信號(hào)，經(jīng)數(shù)字濾波、特征提取等處理與分析，可實(shí)現(xiàn)多種應(yīng)用：,清醒狀態(tài)脈搏信號(hào)及頻譜,疲勞狀態(tài)脈搏信號(hào)及頻譜,時(shí)域：波形密集 頻域：波動(dòng)較大,時(shí)域：波形稀疏 頻域：波動(dòng)平緩,疲勞檢測(cè),利用IIR數(shù)字濾波器處理脈搏信號(hào),人體脈搏,脈搏信號(hào)采集系統(tǒng),輸出信號(hào)往往含有噪聲,基線漂移、人體呼吸等低頻干擾,50Hz工頻干擾,其他高頻干擾,信號(hào) 譜分析,脈搏信號(hào)處理,設(shè)計(jì)濾波器 去噪,特征檢測(cè) 相關(guān)應(yīng)用,需要對(duì)采集的脈搏信號(hào)進(jìn)行去噪處理!,如何處理?,利用IIR數(shù)字濾波器處理脈搏信號(hào),信號(hào)譜分析實(shí)現(xiàn),脈搏信號(hào)頻譜,%信號(hào)采樣率參數(shù) fs=200; %A

21、/D設(shè)備 T=1/fs; %脈搏時(shí)域信號(hào) x=load(x.txt); M=length(x); t=0:M-1*T; figure(1);plot(t,x); xlabel(time/(s);ylabel(x(t); title(脈搏時(shí)域信號(hào));grid on; %FFT頻譜分析 N=M*2-1; X=fft(x,N); figure(2); f=-M+1:M-1*fs/N; plot(f,abs(fftshift(X); xlabel(frequence/(Hz);ylabel(amplitude); title(脈搏頻譜);grid on;,脈搏時(shí)域波形,利用IIR數(shù)字濾波器處理脈搏信號(hào)

22、,信號(hào)頻譜,脈搏信號(hào)頻譜分析：,結(jié)論： 脈搏信號(hào)頻率范圍： 1-20Hz 噪聲信號(hào)頻率范圍： 0-1Hz，20Hz以上,方案一： 帶通濾波器 方案二： 高通和低通濾波器級(jí)聯(lián),基線漂移等低頻干擾，頻率小于1Hz,工頻干擾，頻率為50Hz,其他高頻干擾20Hz以上,如何去噪?,時(shí)域信號(hào),正常人的脈搏信號(hào)在120Hz的頻率范圍內(nèi),利用IIR數(shù)字濾波器處理脈搏信號(hào),濾波器指標(biāo)確定：濾除0-1Hz和20Hz以上頻率范圍噪聲,過(guò)渡帶小，濾波器階數(shù)高， 實(shí)現(xiàn)成本高。,1.帶通濾波器指標(biāo),2.高通與低通濾波器級(jí)聯(lián)指標(biāo),實(shí)際工程中常采用高通濾波器與低通濾波器級(jí)聯(lián)方案。,高通濾波器,低通濾波器,利用IIR數(shù)字濾波

23、器處理脈搏信號(hào),數(shù)字濾波器實(shí)現(xiàn),%Butterworth高通濾波器濾波； Wp1= 0.9/(fs/2); Ws1= 0.1/(fs/2); n,Wn=buttord(Wp1,Ws1,0.5,40); b,a = butter(n,Wn,high); figure(3);freqz(b,a,512,200); title(Butterworth Highpass Filter) y=filtfilt(b,a,x);figure(4);plot(t,y); xlabel(time/(s);ylabel(x(t);title(Butterworth高通濾波后脈搏信號(hào));grid on; Y=fft

24、(y,N);figure(5);plot(f,abs(fftshift(Y); xlabel(frequence/(Hz);ylabel(amplitude);title(FFT(幅頻曲線);grid on; % Butterworth低通濾波器濾波； Wp = 30/(fs/2);Ws = 50/(fs/2); Rp = 3; Rs = 60; n,Wc = buttord(Wp,Ws,Rp,Rs); b,a = butter(n,Wc,low); figure(6);freqz(b,a,512,200); title(Butterworth Lowpass Filter) q=filter

25、(b,a,y);figure(7);plot(t,q); xlabel(time/(s);ylabel(x(t);title(Butterworth低通濾波脈搏信號(hào));grid on; Y=fft(q,N);figure(8);plot(f,abs(fftshift(Y); xlabel(frequence/(Hz);ylabel(amplitude);title(FFT(幅頻曲線);grid on;,利用IIR數(shù)字濾波器處理脈搏信號(hào),Butterwoth,Chebyshev-I,Elliptic,N=2,N=3,N=3,階數(shù)最小 相位近似線性 低頻噪聲濾除效果好,利用IIR數(shù)字濾波器處理脈搏信號(hào),N=3,Butterwoth,Chebyshev-I,Elliptic,高頻噪聲濾除效果好,N=2,N=3,利用IIR數(shù)字濾波器處理脈搏信號(hào),時(shí)域信號(hào)對(duì)比,信號(hào)頻譜對(duì)比,原始信號(hào),濾波后信號(hào),原始信號(hào)頻譜,濾波后信號(hào)頻譜,利用IIR數(shù)字濾波器處理脈搏信號(hào),信號(hào)的時(shí)域波形,加噪信號(hào)的時(shí)域波形,問(wèn)