afd3d5.m

AFD - Advanced Filter Design using MATLABMiroslav D. Lutovac, Dejan V. Tosicversion 1.00 released 15

function [nD5,dD5,numH,denH,n,a,e,fp] = afd3d5(speck,ninc,filtype)

% afd3d5.m  AFD Design alternative 3D5
% 2:41  31/1/99
%
%   Authors: Dejan V. Tosic, Miroslav D. Lutovac, 1999.02.08
%                 [e-mail]
%                 [e-mail]
%
%   Copyright (c) 1999 by Tosic & Lutovac
%   $Revision: 1.0 $  $Date: 1999/02/08 03:07:42 $
%
%   References:
%        Miroslav D. Lutovac, Dejan V. Tosic, Brian L. Evans
%           Advanced Filter Design for Signal Processing
%                   Using MATLAB and Mathematica

%% calls: 
%% afda2k.m, afdnminq.m, afdxminq.m, afdl.m, afdhp.m, afdsimsp.m

if filtype == 'l'
 Fp = speck(1);
 Fs = speck(2);
 speck(3) = afda2k(10/3*log10(1+speck(3)^2));
 speck(4) = afda2k(10/3*log10(1+speck(4)^2));
 Kpp= speck(3);
 Ks = speck(4);
 Kp = 1/Ks;
 n = afdnminq(Fp, Fs, Kpp, Ks);
 n  = n + ninc;
 [a,fn] = afdxminq(n, Fp, Fs, Kpp, Ks);
 e  = 1/sqrt(afdl(n,a));
 fp = Fp/fn;
 [numH,denH] = afdhp(n,a,e);
 nn = length(numH)-1:-1:0;
 nD51 = (ones(size(nn))/(2*pi*fp)).^nn .*numH;
 nd = length(denH)-1:-1:0;
 dD51 = (ones(size(nd))/(2*pi*fp)).^nd .*denH;
 nD5 = conv(nD51,nD51);
 dD5 = conv(dD51,dD51);
 nD5 = conv(nD5,nD51);
 dD5 = conv(dD5,dD51);
end

if filtype == 'h'
 Fp = speck(1);
 Fs = speck(2);
 speck(3) = afda2k(10/3*log10(1+speck(3)^2));
 speck(4) = afda2k(10/3*log10(1+speck(4)^2));
 Kpp = speck(3);
 Ks = speck(4);
 Kp = 1/Ks;
 n = afdnminq(Fp, Fs, Kpp, Ks);
 n  = n + ninc;
 [a,fn] = afdxminq(n, Fp, Fs, Kpp, Ks);
 e  = 1/sqrt(afdl(n,a));
 fp = Fs*fn;
 [numH,denH] = afdhp(n,a,e);
 if max(size(numH)) < max(size(denH))
   numH = [0 numH];
 end
 nD1h = fliplr(numH);
 dD1h = fliplr(denH);
 nn  = length(nD1h)-1:-1:0;
 nD51 = (ones(size(nn))/(2*pi*fp)).^nn .*nD1h;
 nd  = length(dD1h)-1:-1:0;
 dD51 = (ones(size(nd))/(2*pi*fp)).^nd .*dD1h;
 nD5 = conv(nD51,nD51);
 dD5 = conv(dD51,dD51);
 nD5 = conv(nD5,nD51);
 dD5 = conv(dD5,dD51);
end

if filtype == 'b'
 speck = afdsimsp(speck);
 Fs1 = speck(1);
 Fp1 = speck(2);
 Fp2 = speck(3);
 Fs2 = speck(4);
 speck(5) = afda2k(10/3*log10(1+speck(5)^2));
 speck(6) = afda2k(10/3*log10(1+speck(6)^2));
 speck(7) = afda2k(10/3*log10(1+speck(7)^2));
 Ks1 = speck(5);
 Kpb = speck(6);
 Ks2 = speck(7);
 B = Fp2 - Fp1;
 Fr = sqrt(Fp1*Fp2);
 Fp = Fs1;
 Fs = Fp*(Fs2-Fs1)/B;
 Ks = max([Ks1 Ks2]);
 n = afdnminq(Fp, Fs, Kpb, Ks);
 n  = n + ninc;
 [a,fn] = afdxminq(n, Fp, Fs, Kpb, Ks);

 e  = 1/sqrt(afdl(n,a));
 fp = Fp/fn;
 Kp = 1/Ks;
 fp = (Fp2-Fp1);
 [nD1l,dD1l] = afdhp(n,a,e);
 nn = length(nD1l)-1:-1:0;

 nD1 = (ones(size(nn))/(2*pi*fp)).^nn .*nD1l;
 nd = length(dD1l)-1:-1:0;
 dD1 = (ones(size(nd))/(2*pi*fp)).^nd .*dD1l;

 if max(size(nD1l)) < max(size(dD1l))
   nD1l = [0 nD1l];
 end
 nn = max(size(nD1l));
 lb = zeros(nn,2*nn-1);
 w2 = 4*pi^2*Fp1*Fp2;
 fp = (Fp2-Fp1)/fn;
 w1 = 2*pi*fp;
 ns = 0:length(nD1);

 byquad = 1;
 for ind = 1:nn
   lb(ind,nn-ind+1:nn-ind+max(size(byquad))) = byquad;
   byquad = conv(byquad,[1 0 w2]);
 end
 W1 = (ones(size(ns))*w1).^ns;
 W1 = fliplr(W1);
 nBP = 0;
 dBP = 0;
 for ind=1:nn
  lb0(ind,:) =  lb(ind,:) * W1(ind);
  nBP = nBP + lb0(ind,:) * nD1l(nn-ind+1);
  dBP = dBP + lb0(ind,:) * dD1l(nn-ind+1);
 end
 numH = dD1l;
 denH = dD1l;
 nD51 = nBP;
 dD51 = dBP;
 nD5 = conv(nD51,nD51);
 dD5 = conv(dD51,dD51);
 nD5 = conv(nD5,nD51);
 dD5 = conv(dD5,dD51);
end

if filtype == 'r'
 speck = afdsimsp(speck);
 Fp1 = speck(1);
 Fs1 = speck(2);
 Fs2 = speck(3);
 Fp2 = speck(4);
 speck(5) = afda2k(10/3*log10(1+speck(5)^2));
 speck(6) = afda2k(10/3*log10(1+speck(6)^2));
 speck(7) = afda2k(10/3*log10(1+speck(7)^2));
 Kp1 = speck(5);
 Ksr = speck(6);
 Kp2 = speck(7);
 B = Fs2 - Fs1;
 Fr = sqrt(Fs1*Fs2);
 Fp = Fp1;
 Fs = Fp1*(Fp2-Fp1)/B;
 Kp = min([Kp1 Kp2]);
 n = afdnminq(Fp, Fs, Kp, Ksr);
 n  = n + ninc;
 [a,fn] = afdxminq(n, Fp, Fs, Kp, Ksr);

 e  = 1/sqrt(afdl(n,a));
 fp = Fp;
  Kp = 1/Ks;
 [nD1l,dD1l] = afdhp(n,a,e);
 nn = length(nD1l)-1:-1:0;

 nD1 = (ones(size(nn))/(2*pi*fp)).^nn .*nD1l;
 nd = length(dD1l)-1:-1:0;
 dD1 = (ones(size(nd))/(2*pi*fp)).^nd .*dD1l;

 if max(size(nD1l)) < max(size(dD1l))
   nD1l = [0 nD1l];
 end
 nn = max(size(nD1l));
 lb = zeros(nn,2*nn-1);
 w2 = 4*pi^2*Fs1*Fs2;
 fp = (Fp2-Fp1)*fn;
 w1 = 2*pi*fp;
 ns = 0:length(nD1);

 byquad = 1;
 for ind = 1:nn
   lb(ind,nn-ind+1:nn-ind+max(size(byquad))) = byquad;
   byquad = conv(byquad,[1 0 w2]);
 end
 W1 = (ones(size(ns))*w1).^ns;
 W1 = fliplr(W1);
 for ind=1:nn
  lb0(ind,:) =  lb(ind,:) * W1(ind);
 end
 nD1l = fliplr(nD1l);
 dD1l = fliplr(dD1l);
 nBP = 0;
 dBP = 0;
 for ind=1:nn
  nBP = nBP + lb0(ind,:) * nD1l(nn-ind+1);
  dBP = dBP + lb0(ind,:) * dD1l(nn-ind+1);
 end
 numH = dD1l;
 denH = dD1l;
 nD51 = nBP;
 dD51 = dBP;
 nD5 = conv(nD51,nD51);
 dD5 = conv(dD51,dD51);
 nD5 = conv(nD5,nD51);
 dD5 = conv(dD5,dD51);
end