utility_zpk.m

Talking about adaptive filter digital

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