function y = hbp(hlow,hhigh,sp,N,x)
%HBP {Hlow,Hhigh}-Hour BandPass filter, Nth order filter
% HBP(Hlow,Hhigh,SP,N,X) filters the input vector X, at the sampling
% rate SP, using a Nth order bandpass digital Butterworth
% filter, passing information between the periods Hlow and Hhigh.
% Hlow must be strictly less than < Hhigh.
% All arguments are REQUIRED.  SP must be the sample rate 
% in hours. The MATLAB routine FILTFILT is used to perform 
% the actual filtering.
%
%   INPUT: hlow  - low cut-off (hours)
%          hhigh - high cut-off (hours)
%          sp    - sampling rate (hours)
%          N     - filter order (try 2 first!!)
%          x     - signal
%
%  OUTPUT: y  -  {h1-h2}-hour band passed version of input x
%
% CALL AS: y=hbp(10,40,sp,2,x); % for 10-40-hour band pass, 
%                               % 2nd order filter

% Brian Blanton, Fall 1998
%

if nargin==0 & nargout==0
   disp('Call as: y=hbp(Hlow,Hhigh,sample_rate,N,x);')  
   return
end

if hlow >= hhigh
   error('Hlow MUST BE less than Hhigh')
end

whigh=cutoff(sp,hlow)
wlow=cutoff(sp,hhigh)
%[b,a] = butter(N,[wlow,whigh],'bandpass');
[b,a] = butter(N,[wlow,whigh]);
y = filtfilt(b, a, x);