function func_filt = lowPassFilter(func, Fs, cutoff_f, varargin)
%LOWPASSFILTER   Filter input with low pass filter.
%
% DESCRIPTION:
%       lowPassFilter filters an input time domain signal using a filter
%       with a specified cut-off frequency, stop-band attenutation and
%       transition band width. It uses the Kaiser Windowing method to
%       design the FIR filter, which can be implemented either a zero or
%       linear phase filter.
%       
% USAGE:
%       filtered_signal = LowPassFilter(func, Fs, cutoff_f)
%       filtered_signal = LowPassFilter(func, Fs, cutoff_f, ...)
% 
% INPUTS:
%       func        - data to filter
%       Fs          - sampling frequency [Hz]
%       cutoff_f    - filter cutoff frequency [Hz]
%
% OPTIONAL INPUTS:
%       Optional 'string', value pairs that may be used to modify the
%       default computational settings.
%
%       'Plot'            - boolean controlling whether the amplitude spectrum
%                           is displayed before and after filtering (default =
%                           false)
%       'StopBandAtten'   - attenuation in decibels in the stop band
%                           (default = 60)
%       'TransitionWidth' - size of the transition based on the temporal
%                           sampling frequency (default = 0.1)
%       'ZeroPhase'       - boolean controlling whether a zero phase filter
%                           is used (default = false)
%     
% OUTPUTS:
%       filtered_signal - the filtered time signal
%
% ABOUT:
%       author          - B.T. Cox
%       date            - 31st December 2009
%       last update     - 20th January 2010
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009, 2010 Bradley Treeby and Ben Cox
%
% See also butterworth, filterTimeSeries

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
% 
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details. 
% 
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>. 

% set optional input defaults
num_req_input_variables = 3;
zero_phase = false;
transition_width = 0.1; % as proportion of sampling frequency
stop_band_atten = 60;   % [dB]
plot_filter = false;

% replace with user defined values if provided
if nargin < num_req_input_variables
    error('Incorrect number of inputs');
elseif ~isempty(varargin)
    for input_index = 1:2:length(varargin)
        switch varargin{input_index}
            case 'Plot'
                plot_filter = varargin{input_index + 1};
            case 'StopBandAtten'
                stop_band_atten = varargin{input_index + 1};
            case 'TransitionWidth'
                transition_width = varargin{input_index + 1};
            case 'ZeroPhase'
                zero_phase = varargin{input_index + 1};
            otherwise
                error('Unknown optional input');
        end
    end
end

% correct the stopband attenuation if a zero phase filter is being used
if zero_phase
    stop_band_atten = stop_band_atten/2;
end

% calculate the amplitude spectrum of the input if required for plotting
if plot_filter
    [f func_as] = spectrum(func, Fs);
end

% decide the filter order
N = ceil((stop_band_atten - 7.95) / (2.285*(transition_width*pi)));

% construct impulse response of ideal bandpass filter h(n), a sinc function
fc = cutoff_f/Fs; %normalised cut-off
n = -N/2:N/2-1;
h = 2*fc*sinc(2*pi*fc*n);

% compute Kaiser window parameter beta
if stop_band_atten > 50
    beta = 0.1102*(stop_band_atten - 8.7);
elseif stop_band_atten >= 21
    beta = 0.5842*(stop_band_atten - 21)^0.4 + 0.07886*(stop_band_atten - 21);
else
    beta = 0;
end

% construct the Kaiser smoothing window w(n)
m = 0:N-1;
w = real(besseli(0,pi*beta*sqrt(1-(2*m/N-1).^2)))/real(besseli(0,pi*beta));

% window the ideal impulse response with Kaiser window to obtain the FIR
% filter coefficients hw(n) 
hw = w.*h;

% add some zeros to allow the reverse (zero phase) filtering room to work
L = length(func);    % length of original input signal
func = [zeros(1,N) func];

% apply the filter
func_filt = filter(hw,1,func);
if zero_phase
    func_filt = fliplr(filter(hw,1,func_filt(L+N:-1:1)));
end

% remove the part of the signal corresponding to the added zeros
func_filt = func_filt(N+1:L+N);

% plot the amplitude spectrum and the filter if required
if plot_filter
    
    % plot the input signal
    figure;
    [f_sc scale prefix] = scaleSI(max(f));  %#ok<ASGLU>
    plot(f*scale, func_as, 'k-');
    xlabel(['Frequency [' prefix 'Hz]']);
    ylabel('Signal Amplitude');
    hold on;
    xlim = get(gca, 'XLim');
    ylim = get(gca, 'YLim');
    
    % compute the amplitude spectrum of the filtered signal
    [f func_as] = spectrum(func_filt, Fs);
    
    % plot the filtered signal
    plot(f*scale, func_as, 'b-');
    set(gca, 'XLim', xlim);
    xlabel(['Frequency [' prefix 'Hz]']);
    ylabel('Signal Amplitude');    

    % plot the filter cutoff frequency
    line([cutoff_f*scale, cutoff_f*scale], [0, ylim(2)], 'LineStyle','--', 'Color', 'k');
    legend('input signal', 'filtered signal', 'filter cutoff');
end