function [xf,co,npol,npas,tipe,HABS2,F,EPB]=...
    bandpass(x,Fs,colo,cohi,npol,npas,tipe)
% [xf,HABS2,F,EPB]=BANDPASS(x,Fs,colo,cohi,npol,npas,tipe)
%
% Filters signal 'x' with filter 'tipe' and corner
% frequencies 'cohi' and 'cohi' in Hz with 'npol' the 
% number of poles and in 'npas' passes. Sampling rate is 'Fs' in Hz.
%
% INPUT:
%
% x         The signal
% Fs        Its sampling frequency
% colo      The lower corner frequency
% cohi      The higher corner frequency
% npol      The number of poles
% npas      The number of passes
% tipe      The filter name
%
% OUTPUT:
%
% xf      The filtered signal
%
% Compare in SAC bp butter co 0.05 5 n 2 p 1
%
% Returns the npas frequency response and the effective pass band for
% one or two passes (3 dB level)
%
% You'll see that plot(F,decibel(HABS2)) (this is what freqz plots)
% shows how' you concentrate between cohi and colo at the 
% 3 dB-level
%
% Last modified by fjsimons-at-alum.mit.edu, Nov 21th, 2004

defval('npol',2)
defval('npas',1)
defval('colo',0.05)
defval('cohi',0.50)
defval('Fs',110)
defval('tipe','butter')

%disp(sprintf('BANDPASS %3.3f-%3.3f Hz %i pass %i poles %s',...
%	     colo,cohi,npas,npol,tipe))
						
% Corner frequency is in Hertz, now it is as a fraction of
% half the sampling rate.
Wn=2*[colo cohi]/Fs;

if Wn(2)>=1
  Wn(2)=0.99;
  warning('Frequencies adjusted to keep within the Nyquist rate')
end

if diff(Wn)<0.01
  warning(sprintf('%s\n%s\n%s\n%s',...
  ['Situations that seem to require an exceptionally narrow band'],...
  ['filter can be handled more reliably by decimation, filtering'],...
  ['with a filter of more  moderate band width, and interpolation'],...
  ['to the original sampling rate. (SAC manual)']))
end

[B,A]=feval(tipe,npol,Wn);

[H,F]=freqz(B,A,512,Fs);
HABS2=abs(H).^2;

xf=filter(B,A,detrend(x(:)));

if npas==2
  xf=flipud(filter(B,A,detrend(flipud(xf(:)))));  
  HABS2=HABS2.^2;
end

warning off
EPB=bpmin(decibel(HABS2),F,3);
warning on

co=[colo cohi];