www.gusucode.com > 利用MATLAB GUI设计滤波器界面,可以设计IIR滤波器 > AFD/Utility_zpk.m
function strFilterObject=Utility_zpk(strFilterObject,p2,p3,p4) % Utility_zpk is a subfile of the AnalogFilter GUI collection % % James C. Squire, 2002 % Assistant Professor, Virginia Military Institute % ver 1.0 % Utility_zpk returns finds the zeros, poles, and gain of an analog filter % The filter, strFilterObject, is a global variable %global strFilterObject %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Fminsearch Callbacks % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if nargin==4 % called as part of fminsearch for the elliptic filters switch p4 case 'krat' m=strFilterObject; krat=p2; m = min(1,max(m,0)); if abs(m) > eps & abs(m)+eps < 1 k = ellipke([m,1-m]); r = k(1)./k(2) - krat; elseif abs(m) <= eps % m==0 r = -krat; else % m==1 => r == inf, but can't for non-ieee machines r = 1e20; end strFilterObject = abs(r); case 'vrat' u=strFilterObject; ineps=p2; mp=p3; [s,c] = ellipj(u,mp); strFilterObject = abs(ineps - s/c); otherwise error(['Unknown fminsearch callback in ' mfilename]) end return end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Main Program % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Create local variables n = strFilterObject.nOrder; Wc = strFilterObject.fFc * 2* pi; Gp = strFilterObject.fGp; Gs = strFilterObject.fGs; Gain = strFilterObject.fGain; % Get N-th order normalized prototype switch lower(strFilterObject.sType) case 'bessel' [z,p,k] = NormalBessel(n); case 'butterworth' [z,p,k] = NormalButterworth(n); case 'chebychev i' [z,p,k] = NormalChebychevI(n,Gp); case 'chebychev ii' [z,p,k] = NormalChebychevII(n,Gs); case 'elliptic' [z,p,k] = NormalElliptic(n,Gp,Gs); otherwise error(['Unknown Type in ' mfilename]) end % Scale frequency and change to highpass if necessary switch lower(strFilterObject.sPurpose) case 'lowpass' [z,p,k] = denormalizeLP(z,p,k,Wc); case 'highpass' [z,p,k] = denormalizeHP(z,p,k,Wc); otherwise error(['Unknown Purpose in ' mfilename]) end % sort, pair, and fix roundoff errors strFilterObject.vZeros = cplxpair(z); strFilterObject.vPoles = cplxpair(p); strFilterObject.fK = k*Gain; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Bessel % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [z,p,k] = NormalBessel(n) z = []; switch n case 1 p = -1; k = 1; case 2 p = [-1.10160133059216 + 0.63600982475703*i; -1.10160133059216 - 0.63600982475703*i]; k = 1.61803398874989; case 3 p = [-1.04740916100894 + 0.99926443628064*i; -1.04740916100894 - 0.99926443628064*i; -1.32267579991044]; k = 2.77179327460633; case 4 p = [ -1.37006783055144 + 0.41024971749375*i;-1.37006783055144 - 0.41024971749375*i;-0.99520876435027 + 1.25710573945467*i; -0.99520876435027 - 1.25710573945467*i]; k = 5.25819901024417; case 5 p = [ -0.95767654856268 + 1.47112432073039*i;-0.95767654856268 - 1.47112432073039*i;-1.38087732586044 + 0.71790958762677*i;-1.38087732586044 - 0.71790958762677*i;-1.50231627144748]; k = 11.21283668537052; case 6 p = [-1.57149040361604 + 0.32089637422262*i; -1.57149040361604 - 0.32089637422262*i;-1.38185809759657 + 0.97147189071158*i;-1.38185809759657 - 0.97147189071158*i; -0.93065652294686 + 1.66186326894259*i;-0.93065652294686 - 1.66186326894259*i]; k = 26.62976888914601; case 7 p = [-0.90986778062347 + 1.83645135303639*i;-0.90986778062347 - 1.83645135303639*i; -1.37890321679547 + 1.19156677780065*i; -1.37890321679547 - 1.19156677780065*i; -1.61203876622612 + 0.58924450693146*i; -1.61203876622612 - 0.58924450693146*i; -1.68436817927317]; k = 69.22126457866750; case 8 p = [-1.75740840040164 + 0.27286757510224*i;-1.75740840040164 - 0.27286757510224*i;-1.63693941812687 + 0.82279562513969*i;-1.63693941812687 - 0.82279562513969*i;-1.37384121763736 + 1.38835657587755*i;-1.37384121763736 - 1.38835657587755*i;-0.89286971884713 + 1.99832584364129*i;-0.89286971884713 - 1.99832584364129*i]; k = 1.940261933033208e+002; case 9 p = [-0.87839927616096 + 2.14980052431333*i;-0.87839927616096 - 2.14980052431333*i; -1.36758830979291 + 1.56773371223728*i;-1.36758830979291 - 1.56773371223728*i;-1.65239648457884 + 1.03138956698440*i;-1.65239648457884 - 1.03138956698440*i;-1.80717053496210 + 0.51238373057492*i;-1.80717053496210 - 0.51238373057492*i;-1.85660050122801]; k = 5.801750529768317e+002; case 10 p = [-1.92761969137237 + 0.24162347097177*i;-1.92761969137237 - 0.24162347097177*i;-1.84219624452473 + 0.72725759775741*i;-1.84219624452473 - 0.72725759775741*i;-1.66181024136221 + 1.22110021857916*i;-1.66181024136221 - 1.22110021857916*i;-1.36069227838455 + 1.73350574266150*i;-1.36069227838455 - 1.73350574266150*i;-0.86575690170838 + 2.29260483098249*i;-0.86575690170838 - 2.29260483098249*i]; k = 1.836246770717193e+003; case 11 p = [-0.85451258135236 + 2.42805946691507*i;-0.85451258135236 - 2.42805946691507*i;-1.35348667738810 + 1.88829684475974*i;-1.35348667738810 - 1.88829684475974*i;-1.66719364224535 + 1.39596290364648*i;-1.66719364224535 - 1.39596290364648*i;-1.86736123888712 + 0.92311558295178*i;-1.86736123888712 - 0.92311558295178*i;-1.98016064530384 + 0.45959874382676*i;-1.98016064530384 - 0.45959874382676*i;-2.01670147345004]; k = 6.114589477119062e+003; case 12 p = [ -2.08464450693174 + 0.21916153519047*i;-2.08464450693174 - 0.21916153519047*i;-2.01994593307477 + 0.65899650071675*i;-2.01994593307477 - 0.65899650071675*i;-1.88564961973203 + 1.10381488121541*i;-1.88564961973203 - 1.10381488121541*i;-1.66980358885276 + 1.55880270084110*i;-1.66980358885276 - 1.55880270084110*i;-1.34617468029050 + 2.03399850827847*i;-1.34617468029050 - 2.03399850827847*i;-0.84437887285039 + 2.55718896920289*i;-0.84437887285039 - 2.55718896920289*i]; k = 2.131985288096244e+004; case 13 p = [-0.83515201045011 + 2.68080279689105*i;-0.83515201045011 - 2.68080279689105*i;-1.33888032408279 + 2.17202294710701*i;-1.33888032408279 - 2.17202294710701*i;-1.67045856262298 + 1.71167828638810*i;-1.67045856262298 - 1.71167828638810*i;-1.89898611748036 + 1.27211943651370*i;-1.89898611748036 - 1.27211943651370*i;-2.05058087608408 + 0.84338310857716*i;-2.05058087608408 - 0.84338310857716*i;-2.13764829462372 + 0.42041630675713*i;-2.13764829462372 - 0.42041630675713*i;-2.16608270567884]; k = 7.752998977990125e+004; case 14 p = [-2.23093074227585 + 0.20200027007892*i;-2.23093074227585 - 0.20200027007892*i;-2.17970952065207 + 0.60702982702828*i;-2.17970952065207 - 0.60702982702828*i;-2.07445158463977 + 1.01536708959161*i;-2.07445158463977 - 1.01536708959161*i;-1.90866457514154 + 1.43007973190289*i;-1.90866457514154 - 1.43007973190289*i;-1.66971016070877 + 1.85613874359928*i;-1.66971016070877 - 1.85613874359928*i;-1.33167919736814 + 2.30345534461641*i;-1.33167919736814 - 2.30345534461641*i;-0.82668133243230 + 2.79955221217277*i;-0.82668133243230 - 2.79955221217277*i]; k = 2.930813527521330e+005; case 15 p = [ -0.81885183033332 + 2.91396992652773*i;-0.81885183033332 - 2.91396992652773*i;-2.28342656234184 + 0.38982893598505*i;-2.28342656234184 - 0.38982893598505*i;-2.21352748734394 + 0.78142997785105*i;-2.21352748734394 - 0.78142997785105*i;-2.09319972155537 + 1.17691617609367*i;-2.09319972155537 - 1.17691617609367*i;-1.91558434659226 + 1.57926035556851*i;-1.91558434659226 - 1.57926035556851*i;-1.66794080077851 + 1.99338080832681*i;-1.66794080077851 - 1.99338080832681*i;-1.32461676006183 + 2.42914977398455*i;-1.32461676006183 - 2.42914977398455*i;-2.30637005662894]; k = 1.148461735754945e+006; case 16 p = [ -2.36834668175327 + 0.18833295674897*i;-2.36834668175327 - 0.18833295674897*i;-2.32647906263646 + 0.56573669095214*i;-2.32647906263646 - 0.56573669095214*i;-2.24099322267511 + 0.94547404667305*i;-2.24099322267511 - 0.94547404667305*i;-2.10799083457213 + 1.32955250584788*i;-2.10799083457213 - 1.32955250584788*i;-1.92038809489843 + 1.72088333292392*i;-1.92038809489843 - 1.72088333292392*i;-1.66542192928512 + 2.12434959299958*i;-1.66542192928512 - 2.12434959299958*i;-1.31771943728305 + 2.54979207329592*i;-1.31771943728305 - 2.54979207329592*i;-0.81157346724783 + 3.02449807602193*i;-0.81157346724783 - 3.02449807602193*i]; k = 4.653711491998211e+006; otherwise error(['Order of bessel filter in ' mfilename ' cannot be greater than 16']) end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Butterworth % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [z,p,k] = NormalButterworth(n) z = []; p = exp(i*(pi*(1:2:n-1)/(2*n) + pi/2)); p = [p; conj(p)]; p = p(:); if rem(n,2)==1 % n is odd p = [p; -1]; end k = real(prod(-p)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Chebychev I % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [z,p,k] = NormalChebychevI(n,Gp) epsilon = sqrt(10^(.1*Gp)-1); mu = asinh(1/epsilon)/n; p = exp(j*(pi*(1:2:2*n-1)/(2*n) + pi/2)).'; p = sinh(mu)*real(p) + j*cosh(mu)*imag(p); z = []; k = real(prod(-p)); if ~rem(n,2) % n is even so patch k k = k/sqrt((1 + epsilon^2)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Chebychev II % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [z,p,k] = NormalChebychevII(n,Gs) delta = 1/sqrt(10^(.1*Gs)-1); mu = asinh(1/delta)/n; if (rem(n,2)) m = n - 1; z = j ./ cos([1:2:n-2 n+2:2:2*n-1]*pi/(2*n))'; else m = n; z = j ./ cos((1:2:2*n-1)*pi/(2*n))'; end p = exp(j*(pi*(1:2:2*n-1)/(2*n) + pi/2)).'; p = sinh(mu)*real(p) + j*cosh(mu)*imag(p); p = 1 ./ p; k = real(prod(-p)/prod(-z)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Elliptic % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [z,p,k] = NormalElliptic(n,Gp,Gs) if n == 1, z = []; p = -sqrt(1/(10^(Gp/10)-1)); k = -p; return else % Zeros epsilon = sqrt(10^(0.1*Gp)-1); k1 = epsilon/sqrt(10^(0.1*Gs)-1); k1p = sqrt(1-k1^2); if k1p == 1 error('Cannot design filter, Gp and Gs specifications are too strict.') end wp = 1; if abs(1-k1p^2) < eps krat = 0; else capk1 = ellipke([k1^2,k1p^2]); krat = n*capk1(1)/capk1(2); end fopt = optimset('maxfunevals',250,'maxiter',250,'display','none'); m = fminsearch('Utility_zpk',.5,fopt,krat,[],'krat'); if m<0 | m>1 m = fminbnd('Utility_zpk',0,1,fopt,krat,[],'krat'); end capk = ellipke(m); ws = wp/sqrt(m); m1 = 1 - m; j = (1-rem(n,2)):2:n-1; [ij,jj] = size(j); [s,c,d] = ellipj(j*capk/n,m*ones(ij,jj)); is = find(abs(s) > eps); z = 1 ./(sqrt(m)*s(is)); z = i*z(:); z = [z ; conj(z)]; % Poles r = fminsearch('Utility_zpk', ellipke(1-m),fopt,1/epsilon,k1p^2,'vrat'); v0 = capk*r/(n*capk1(1)); [sv,cv,dv] = ellipj(v0,1-m); p = -(c.*d*sv*cv + i*s*dv)./(1-(d*sv).^2); p = p(:); if rem(n,2) ip = find(abs(imag(p)) < eps*norm(p)); [pm,pn] = size(p); pp = 1:pm; pp(ip) = []; p = [p ; conj(p(pp))]; else p = [p; conj(p)]; end % Gain k = real(prod(-p)/prod(-z)); if (~rem(n,2)) k = k/sqrt((1 + epsilon^2)); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Lowpass Denormalize % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [z,p,k] = denormalizeLP(z,p,k,w) Np = length(p); Nz = length(z); k = k*w^(Np-Nz); p = p*w; z = z*w; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Highpass Denormalize % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[z,p,k] = denormalizeHP(z,p,k,w) Np = length(p); Nz = length(z); k = abs( k * prod(z) / prod(p) ); z = w./z; p = w./p; if Np > Nz % add extra zeros at zero z = [zeros(Np-Nz,1); z]; end