function [y,fout] = blimit(x,fin,band)
% BLIMIT
%   Basebands a specified frequency band within data, shifting
%   it to zero to the width of the band.  Also removes the
%   delay caused by lowpass filtering.
%
% [y,fout] = blimit(x,fin,band)
% 
% x          input time series
% fin        sample rate of input
% band
%   .freqs   [flow,fhigh] of final bandwidth signal
%   .blp     (optional) filter coeffs for lowpass filter. If
%            omitted uses an order 11 Hamming window low-pass filter
%
% y          real output that has mapped [flow,fhigh] to [0,fhigh-flow]
% fout       sample rate of output
%
% USES: fir1, resample
%
% $Id: blimit.m,v 1.8 2004/08/16 16:29:38 lsf Exp $

% high and low edges of band
flow = min(band.freqs);
fhigh = max(band.freqs);

% midpoint frequency of band
fmid = mean(band.freqs);

% bandwidth
deltaf = abs(diff(band.freqs));

% resampling parameters
[p,q] = rat(2*deltaf/fin);
fout = fin*p/q;

% shift down by fmid
t = (0:(length(x)-1));
x = x.*exp(2*pi*i*t*fmid/fin);

% resample to set fs = 2*deltaf
x = [0 x];
x = resample(x,p,q); 
x = x(2:end);

% lowpass filter to remove frequencies outside band
% also remove shift casued by lowpass filter
if (~isfield(band,'blp'))
  band.blp = fir1(11,0.5);
end
fshift = floor(length(band.blp)/2); %shift caused by filter
x = [x zeros(1,fshift)];
x = filter(band.blp,1,x);
x = x(fshift+1:end);

% shift up by deltaf/2
t = (0:(length(x))-1);
x = x.*exp(-2*pi*i*t*(deltaf/2)/fout);
y = 2*real(x);                          % correct power in band

return