1、8 IIR Filter Design Butterworth Filter DesignStep1 ：設計一 c=0.5 的三階 Butterworth 類比原型 Filter z,p,k=buttap(N) %Matlab 提供之函數用以設計 c=1 的 N 階 Butterworth 類比（Analog）原型（ Prototype）Filter function b,a = u_buttap(N,Omegac); % 根據 buttap(N)，進一步使 c 不限為 1% Unnormalized Butterworth Analog Lowpass Filter Prototype% b
2、,a = u_buttap(N,Omegac);% b = numerator polynomial coefficients of Ha(s)% a = denominator polynomial coefficients of Ha(s)% N = Order of the Butterworth Filter% Omegac = Cutoff frequency in radians/secz,p,k = buttap(N);p = p*Omegac;k = k*OmegacN;B = real(poly(z);b0 = k;b = k*B;a = real(poly(p); func
3、tion C,B,A = sdir2cas(b,a); %將 H(s)直接式轉換成串接形式% DIRECT-form to CASCADE-form conversion in s-plane% C,B,A = sdir2cas(b,a);% C = gain coefficient% B = K by 3 matrix of real coefficients containing bks% A = K by 3 matrix of real coefficients containing aks% b = numerator polynomial coefficients of DIREC
4、T form% a = denominator polynomial coefficients of DIRECT formNa = length(a)-1; Nb = length(b)-1;% compute gain coefficient Cb0 = b(1); b = b/b0;a0 = a(1); a = a/a0;C = b0/a0;% Denominator second-order sections:p = cplxpair(roots(a); K = floor(Na/2);if K*2 = Na % Computation when Na is evenA = zeros
5、(K,3);for n=1:2:NaArow = p(n:1:n+1,:);Arow = poly(Arow);A(fix(n+1)/2),:) = real(Arow);endelseif Na = 1 % Computation when Na = 1A = 0 real(poly(p);else % Computation when Na is odd and 1A=zeros(K+1,3);for n=1:2:2*KArow = p(n:1:n+1,:);Arow = poly(Arow);A(fix(n+1)/2),:) = real(Arow);endA(K+1,:) = 0 re
6、al(poly(p(Na);end% Numerator second-order section:z = cplxpair(roots(b); K = floor(Nb/2);if Nb = 0 % Computation when Nb = 0B = 0 0 poly(z);elseif K*2 = Nb % Computation when Nb is evenB=zeros(K,3);for n=1:2:NbBrow = z(n:1:n+1,:);Brow = poly(Brow);B(fix(n+1)/2),:) = real(Brow);endelseif Nb = 1 % Com
7、putation when Nb = 1B = 0 real(poly(z);else % Computation when Nb is odd and 1B = zeros(K+1,3);for n=1:2:2*KBrow = z(n:1:n+1,:);Brow = poly(Brow);B(fix(n+1)/2),:) = real(Brow);endB(K+1,:) = 0 real(poly(z(Nb);end Matlab Program:N=3; OmegaC=0.5; b,a=u_buttap(N,OmegaC); C,B,A=sdir2cas(b,a);結果：a = 1.000
8、0 1.0000 0.5000 0.1250b = 0.1250A = 1.0000 0.5000 0.25000 1.0000 0.5000B = 0 0 1C = 0.1250Step1 ：設計一滿足以下規格的 Butterworth 類比原型 Filter通帶截止：p=0.2 ; 通帶漣波：Rp=7dB阻帶截止：s=0.3; 阻帶漣波：As=16dB註： 其實經初步計算結果，設計一 c 約 0.5 的三階 Butterworth （同上例） ，即可完成 function b,a = afd_butt(Wp,Ws,Rp,As); % 比 u_buttap 函數更廣義的設計函數% Analo
9、g Lowpass Filter Design: Butterworth% b,a = afd_butt(Wp,Ws,Rp,As);% b = Numerator coefficients of Ha(s)% a = Denominator coefficients of Ha(s)% Wp = Passband edge frequncy in rad/sec; Wp 0% Ws = Stopband edge frequncy in rad/sec; Ws Wp 0 % Rp = Passband ripple in +dB; (Rp 0)% As = Stopband attenuati
10、on in +dB; (As 0)if Wp H(z)function b,a = imp_invr(c,d,T)% Impulse Invariance Transformation from Analog to Digital Filter% -% b,a = imp_invr(c,d,T)% b = Numerator polynomial in z(-1) of the digital filter% a = Denominator polynomial in z(-1) of the digital filter% c = Numerator polynomial in s of t
11、he analog filter% d = Denominator polynomial in s of the analog filter% T = Sampling (transformation) parameter%R,p,k = residue(c,d);p = exp(p*T);b,a = residuez(R,p,k);b = real(b); a = real(a);例： T=0.1651)(2ssHMatlab: subplot(1,1,1);c = 1,1; d = 1,5,6; T = 0.1; Fs = 1/T;b,a = imp_invr(c,d,T)%b = 1.0
12、000 -0.8966%a = 1.0000 -1.5595 0.6065% Impulse response of the analog filtert = 0:0.01:3; ha,ta = impulse(c,d,t);subplot(2,1,1); plot(t,ha,t,0.*t,k); axis(0,3,-0.1,1);hold on% Impulse response of the digital filtern = 0:1:3/T; hn = filter(b,a,impseq(0,0,3/T);stem(n*T,hn); xlabel(time in sec); title
13、(Impulse Responses);hold off% Magnitude Response of the digital filterdb,magd,pha,grd,wd = freqz_m(b,a);% magnitude response of the analog filterdb,mags,pha,ws = freqs_m(c,d,2*pi*Fs);subplot(2,1,2); plot(ws/(2*pi),mags*Fs,wd/(2*pi)*Fs,magd)xlabel(frequency in Hz); title(Magnitude Responses);ylabel(M
14、agnitude); text(1.4,.5,Analog filter); text(1.5,1.5,Digital filter)0 1 5 302 4 6 8 1 c 0 1 2 3 4 5 6 7 8 9 0123z s e r 例： T=0.13)(sHaMatlab: c = 3; d = 1 3; T = 0.1; Fs = 1/T;b,a = imp_invr(c,d,T)結果：b = 3a = 1.0000-0.7408Step2 : 使用 s - z 轉換例： 使用 Butterworth 原型設計滿足以下規格之 LP Digital Filterwp=0.2; Rp=1d
15、Bws=0.3; As=15dB註： 以上給的是 H(z)頻譜的規格，需先轉回 H(s)上的規格，再利用 Step1 中的做法，設計出 Butterworth Analog Filter H(s)，之後再將 H(s)轉成 H(z) Matlab Program:% Digital Filter Specifications:wp = 0.2*pi; % digital Passband freq in Hzws = 0.3*pi; % digital Stopband freq in HzRp = 1; % Passband ripple in dBAs = 15; % Stopband at
16、tenuation in dB% Analog Prototype Specifications: Inverse mapping for frequenciesT = 1; % Set T=1OmegaP = wp * T; % Prototype Passband freqOmegaS = ws * T; % Prototype Stopband freqep = sqrt(10(Rp/10)-1); % Passband Ripple parameterRipple = sqrt(1/(1+ep*ep); % Passband RippleAttn = 1/(10(As/20); % S
17、topband Attenuation% Analog Butterworth Prototype Filter Calculation:cs,ds = afd_butt(OmegaP,OmegaS,Rp,As); % Impulse Invariance transformation:b,a = imp_invr(cs,ds,T);C,B,A = dir2par(b,a)% Plottingfigure(1); subplot(1,1,1)db,mag,pha,grd,w = freqz_m(b,a);subplot(2,2,1); plot(w/pi,mag); title(Magnitu
18、de Response)xlabel(frequency in pi units); ylabel(|H|); axis(0,1,0,1.1)set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickmode,manual,YTick,0,Attn,Ripple,1); gridsubplot(2,2,3); plot(w/pi,db); title(Magnitude in dB);xlabel(frequency in pi units); ylabel(decibels); axis(0,1,-40,5);set(gca,XTick
19、Mode,manual,XTick,0,0.2,0.3,1);set(gca,YTickmode,manual,YTick,-15,0); gridsubplot(2,2,2); plot(w/pi,pha/pi); title(Phase Response)xlabel(frequency in pi units); ylabel(pi units); axis(0,1,-1,1);set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickmode,manual,YTick,-1,0,1); gridsubplot(2,2,4); pl
20、ot(w/pi,grd); title(Group Delay)xlabel(frequency in pi units); ylabel(Samples); axis(0,1,0,10)set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickmode,manual,YTick,0:2:10); grid執行結果：* Butterworth Filter Order = 6 C = B = 1.8557 -0.6304-2.1428 1.14540.2871 -0.4466A = 1.0000 -0.9973 0.25701.0000 -1.0691 0.36991.0000 -1.2972 0.69490 2 3 10 1 s | 0 2 3 15 0B s 01 İ 2 3 68 y s
