afdd1.m

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

function [nD1,dD1,numH,denH,n,a,e,fp] = afdd1(speck,ninc,filtype)

% afdd1.m  AFD Design alternative D1
% 3:10  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:
%% afdnell.m, afdhp.m, afdsimsp.m

if filtype == 'l'
 Fp = speck(1);
 Fs = speck(2);
 Kp = speck(3);
 Ks = speck(4);
 n  = afdnell(Fp,Fs,Kp,Ks);
 n  = n + ninc;
 a  = Fs/Fp;
 e  = Kp;
 fp = Fp;
 [numH,denH] = afdhp(n,a,e);
 nn = length(numH)-1:-1:0;
 nD1 = (ones(size(nn))/(2*pi*fp)).^nn .*numH;
 nd = length(denH)-1:-1:0;
 dD1 = (ones(size(nd))/(2*pi*fp)).^nd .*denH;
end

if filtype == 'h'
 Fp = speck(1);
 Fs = speck(2);
 Kp = speck(3);
 Ks = speck(4);
 n  = afdnell(Fp,Fs,Kp,Ks);
 n  = n + ninc;
 a  = Fs/Fp;
 e  = Kp;
 fp = Fs;
 [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;
 nD1 = (ones(size(nn))/(2*pi*fp)).^nn .*nD1h;
 nd = length(dD1h)-1:-1:0;
 dD1 = (ones(size(nd))/(2*pi*fp)).^nd .*dD1h;
end

if filtype == 'b'
 speck = afdsimsp(speck);
 Fs1 = speck(1);
 Fp1 = speck(2);
 Fp2 = speck(3);
 Fs2 = speck(4);
 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]);
 Kp = Kpb;
 n = afdnell(Fp,Fs,Kp,Ks);
 n  = n + ninc;
 a  = Fs/Fp;
 e  = Kp;
 fp = Fs;
 [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;
 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;
 nD1 = nBP;
 dD1 = dBP;
end

if filtype == 'r'
 speck = afdsimsp(speck);
 Fp1 = speck(1);
 Fs1 = speck(2);
 Fs2 = speck(3);
 Fp2 = speck(4);
 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]);
 Ks = Ksr;
 n = afdnell(Fp,Fs,Kp,Ks);
 n  = n + ninc;
 a  = Fs/Fp;
 e  = Kp;
 fp = Fp;
 [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;
 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;
 nD1 = nBP;
 dD1 = dBP;
end