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
