dfirauwi.m

ADSP TOOLBOX: Version 2.0 and gui m-files

function [n,fc,f3,p1] = dfirauwi(ty,a,samp,fp,fs,wind,p1)
% FIRAUWIN Auto window-based design
%
%       [N,Fc,F0,P1] = firauwin(TY,A,SF,FP,FS,WIND)
%       INPUTS:
%	TY = 'lp', 'bp', 'bs', 'hp', A = attenuation array [Ap,As],
%       SF=sampling freq., FP=passband edge(s), FS=stopband edge(s) (Hz),
%	WIND = window name in quotes (DEFAULT: 'hamming')
%
%       OUTPUTS:
%       N=filter length,
%       Fc=LPP filter cutoff
%       F0=center frequency (for BP and BS filters only)
%	P1 = parameter for Kaiser and Dolph(Cheby) window

% Author: Ashok Ambardar <[e-mail]>
% Date: 2/2/96 11:32AM
% Copyright (c) 1997 by Prentice-Hall, Inc.


%p1=[];
odd=0;
ap=a(1);
as=a(2);
ass=as;
if ty==4,
   odd=1;
end
if wind==5;  %Cheby
   odd=1;
   p1=as-5;
end

%BETA (p1) FOR KAISER WINDOW
as0=as;
if ty==4,
   as0=as+2;
end
dels=10^(-.05*as);
delp=(10^(0.05*ap)-1)/(10^(0.05*ap)+1);

if delp<dels,
   as0=-20*log10(delp);
end,

if as0>=50,
   b=0.1102*(as0-8.7)/pi;
elseif as0>21 & as0<50,
   b=(0.5842*(as0-21)^(0.4) + 0.07886*(as0-21))/pi;
else,
   b=0;
end

flg=0;
if wind==4,
   p1=b;
end  % Kaiser

[wp,ws,w3]=dfirlpp(ty,2*pi*fp,2*pi*fs,samp);

n=2*fix(pi*2/(ws-wp)/1.3)-3;
n=max(3,n);
wx=(wp+ws)/2;
sc0=1.25;
wx0=wx/sc0;
wc=2;
as=ass;
a2=0;
a3=2*ap;

%%%%% Finding N
while a2<as | a3>ap
      if odd==0,
         n=n+1;
      else,
         n=n+2;
      end
n0=fix(n/2);
n1=1:n0;
      if rem(n,2)==1,
         h=sin(wx*n1)./n1/pi;
         h=[fliplr(h) wx/pi h];
      else,
         h=sin(wx*(n1-.5))./(n1-.5)/pi;
         h=[fliplr(h) h];
      end
hwin = winwind(wind,n,p1);
hw=h.*hwin;
w0=[wp ws 0:ws/100:ws];
fr=polyval(hw,exp(sqrt(-1)*w0));
a0=-20*log10(abs(fr));
a0=a0-min(a0);
ftr=w0(1:2)/2/pi;
w0(1:2)=[];
f0=w0/2/pi;
atr=a0(1:2);
a2=a0(2);
a3=0.922*a0(1);
    if n>50,
       if flg<1,
          a3=.98*a0(1);
       end,
    end

a0(1:2)=[];
l0=length(a0);
   if ws>2.9*wp,
      wx=wx0;
   else,
      i=find(a0>ap);
      if isempty(i),
         k=1;
      else,
         k=i(1);
      end
      if k==1,
         sc=1;
      else,
         sc=wp/w0(k);
      end,
      wx=wx*sc;
   end

if wx>ws,
   wx=wx0;
end

end

%%%%% Finding Fc
wxx=.5*(wp+ws);
if ty==4,
   add0=.02;
else,
   add0=-.02;
   wxx=1.2*wxx;
end
wx=wxx;
add=add0;

n0=fix(n/2);
n1=1:n0;
hwin = winwind(wind,n,p1);
num=int2str(n);acp=2*ap;ac0=as;ap0=acp;


acomp=0;flag=0;

while acomp<as
     if flag==0,
        if as-acomp<10,
           add=add/10;
           flag=1;
         end,
     end

wx=wx+add;

if acomp-as>10 | wx<0.8*wp,
%   disp('Cannot meet requirements'),
   wx=wxx;
   break,
end

%if wx<0.8*wp;
%  disp('Cannot meet requirements'),
%   wx=wxx;
%   break,
%end

ac0=acomp;
ap0=acp;
if rem(n,2)==1,
   h=sin(wx*n1)./n1/pi;
   h=[fliplr(h) wx/pi h];
else,
   h=sin(wx*(n1-.5))./(n1-.5)/pi;
   h=[fliplr(h) h];
end,
hw=h.*hwin;
n0=fix(n/2);
if ty==1,
   hwn=hw;
end
if rem(n,2)==1,
   nn=-n0:n0;
else,
   nn=[-n0+0.5:-0.5, 0.5:n0-0.5];
end
if ty==2,
   hwn=hw.*((-1).^nn);
end,
if ty==3,
   hwn=2*cos(nn*w3).*hw;
end
if ty==4,
   hwn=-2*cos(nn*w3).*hw;
   hwn(n0+1)=1+hwn(n0+1);
end
hwn=hwn.*(abs(hwn)>10*eps);
w1=2*pi*fp/samp;
w2=2*pi*fs/samp;
l=length([w1 w2]);
w0=[w1 w2 0:pi/100:pi];
fr=polyval(hwn,exp(sqrt(-1)*w0));
mag=abs(fr);
a0=20*log10(mag);
fr(1:l)=[];
ftr=w0(1:l)/2/pi;
w0(1:l)=[];
atr=a0(1:l);
a0(1:l)=[];
ma=mag(1:l);
mag(1:l)=[];
am=max(a0);
f0=w0/2/pi;
a0=a0-am;
atr=atr-am;
l0=length(a0);
if length(atr)<3,
   acomp=abs(min(atr));
   acp=abs(max(atr));
else,
   acomp=min(abs(atr(3:4)));
   acp=max(abs(atr(1:2)));
end

end
fc=wx/2/pi;f3=w3/2/pi;