%       BANDPASS              05-04-89         REN-CHIEH LIEN
%       
%         BANDX = BANDPASS(X,DT,LCTIME,HCTIME,NPOLE,plotind,FILTERTYPE)
%        
%           BANDX : BANDPASSED TIME SERIES
%               X : ORIGINAL TIME SERIES
%          LCTIME : LOWBOUND CUTOFF TIME     
%          HCTIME : HIGHBOUND CUTOFF TIME     
%           NPOLE : # OF POLES FOR THE NONLINEAR FILTER
%==========================================================================
%
   function [w,minff,maxff,ff,hl]= ...
              lowpass(x,delt,lowcutt,npole,plotind,filtype),
%
%--------------------------------------------------------------------------
%            define the sampling frequency and cutoff frequency band
%--------------------------------------------------------------------------
%
     mx = mean(x);
     x= x-mx; 
     fs = 1/delt;
     lowband = 2*delt/lowcutt;
%     higband = 2*delt/higcutt; 
%     lband = [lowband higband];
     lband = lowband;
     time = (0:(length(x)-1))*delt;
%--------------------------------------------------------------------------
%             design Butterworth filter 
%--------------------------------------------------------------------------
     if ( filtype == 'b' | filtype == 'B');
         [cl,ch] = butter(npole,lband,'low');
     end
     if ( filtype == 'c' | filtype == 'C');
         [cl,ch] = cheby1(npole,0.5,lband,'low');
     end
     nfreq = 5000;
     ff = 0.5*fs*(0:nfreq-1)/(nfreq);
     hl = freqz(cl,ch,nfreq);
     hl = hl.*conj(hl);
     minind = min(find(hl>=0.8));
     maxind = max(find(hl>=0.8));
     
     minff = ff(minind);
     maxff = ff(maxind);
%
%--------------------------------------------------------------------------
%             apply filter twice to avoid phase shift
%--------------------------------------------------------------------------
%
     i=1:length(x);
     y = filter(cl,ch,x);
     %figure;plot(i,x,i,y)
     z0 = y(length(y):-1:1);
     z = filter(cl,ch,z0);
     w = z(length(z):-1:1);
%     figure;plot(i,w)
     xband = [lowband/2/delt lowband/2/delt]; 
     if (plotind == 'y' | plotind =='Y'),
        figure
        semilogx(ff(2:nfreq),hl(2:nfreq),xband,[0 1])
%        semilogx(ff(2:nfreq),hl(2:nfreq))
%        plot(1./ff(2:nfreq)./12,hl(2:nfreq),1./xband./12,[0 1 1 0])
        xlabel('FREQUENCY   [cpm]')
        ylabel('TRANSFER   FUNCTION')
	axis([-Inf Inf 0 1]);
	drawnow
    figure;plot(time,x,time,w)

     end
%    title('BANDPASS  FILTER  FOR  DIURNAL  TIDE')
%    clg
     end
