function func_filt = butterworth(func, Fs, cutoff_f, varargin)
%BUTTERWORTH   Apply a frequency domain Butterworth filter.
%
% DESCRIPTION:
%       butterworth applies frequency domain Butterworth filter to the
%       input func using the user defined cutoff_f. A 5th order filter is
%       used by default.
%
% USAGE:
%       func_filt = butterworth(func, Fs, cutoff_f, varargin)
%
% 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.
%
%       'Order'     - filter order (default = 5)
%       'Plot'      - boolean controlling whether the amplitude spectrum
%                     is displayed before and after filtering (default =
%                     false)
%
% OUTPUTS:
%       func_filt   - filtered function
%
% ABOUT:
%       author      - Bradley E. Treeby
%       date        - 2nd December 2009
%       last update - 3rd December 2009
%
% 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 fft, smooth, spectrum

% 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/>. 

% define filter order
filter_order = 5;

% set optional input defaults
num_req_input_variables = 3;
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 'Order'
                filter_order = varargin{input_index + 1};
            case 'Plot'
                plot_filter = varargin{input_index + 1};    
            otherwise
                error('Unknown optional input');
        end
    end
end

% calculate the relative cutoff frequency
[f func_as] = spectrum(func, Fs, 'Window', true);
cutoff_f_rel = 0.5 * cutoff_f / max(f(:));

% create the filter
len = length(func);
f_filt = (-0.5:1/(len - 1):0.5);
LP_filt = 1./(1.0 + (f_filt./cutoff_f_rel).^(2*filter_order));

% apply the filter enforcing the correct phase offset
if mod(len, 2)
    func_filt = real(ifft(ifftshift(fftshift(fft(func)).*LP_filt)));
else
    func_filt = real(ifft(fftshift(fftshift(fft(func)).*LP_filt)));
end

% plot the amplitude spectrum and the filter if required
if plot_filter
    figure;
    [f_sc scale prefix] = scaleSI(max(f));
    subplot(2, 1, 1), plot(f*scale, func_as, 'k-');
    xlabel(['Frequency [' prefix 'Hz]']);
    ylabel('Signal Amplitude');
    hold on;
    xlim = get(gca, 'XLim');
    ylim = get(gca, 'YLim');
    plot(f*scale, fliplr(LP_filt(1:length(f)))*max(func_as(:)), 'b-');
    line([cutoff_f*scale, cutoff_f*scale], [0, ylim(2)], 'LineStyle','--', 'Color', 'b');
    title('Amplitude Spectrum Of Input Signal');
    legend('signal', 'filter', 'cutoff');
    
    % compute the amplitude spectrum of the filtered signal
    [f func_as] = spectrum(func_filt, Fs, 'Window', true);
    
    % plot
    subplot(2, 1, 2), plot(f*scale, func_as, 'k-');
    set(gca, 'XLim', xlim);
    xlabel(['Frequency [' prefix 'Hz]']);
    ylabel('Signal Amplitude');    
    title('Amplitude Spectrum Of Filtered Signal');
end