
function [dataFiltered] = bandFilter( data, lowFreq, highFreq, Fs )

%BANDFILTER	Filter Medac data with FFT based bandpass.
%	[filteredData] = BANDFILTER( data, lowFreq, highFreq, Fs [10] )
%	returns the bandpass filtered 'data' in the frequencies
%	defined by 'lowFreq' and 'highFreq' (in Hz), and assuming the
%	sampling frequency defined by Fs (in KHz).


%--- Default parameters ---

if( nargin < 3), error( 'This function requieres at least 3 arguments.' ), end
if( nargin < 4 ), Fs = 10;, end


%--- Test compatibility of sampling frequency ---

testFreq( data, Fs )


%--- Filter parms ---

lengthTrace = size( data, 1 );
nCutLow = lowFreq * lengthTrace / (Fs * 1000);
nCutHigh = highFreq * lengthTrace / (Fs * 1000);


%--- Extract  first column from data ---

timestamps = data( :, 1 );


%--- Filter the data ---

F = fft( data );									% Start with fft

for k = 1:round( lengthTrace / 2 )			% Zero out selected components
   if( k < nCutLow | k > nCutHigh )
      F( k + 1, : ) = complex(0);					   % Set this component to zero...
      F( lengthTrace - k + 1, : ) = complex(0);		% ...for all channels
   end
end


%--- Restore timestamp column ---

dataFiltered = real( ifft(F) );
dataFiltered( :, 1 ) = timestamps;

