function f_sm = smooth(kgrid, f, restore_max)
%SMOOTH     Smooth a matrix.
%
% DESCRIPTION:
%       smooth filters an input matrix using an n-dimensional
%       frequency-domain Hanning window.
%
% USAGE:
%       f_sm = smooth(kgrid, f)
%       f_sm = smooth(kgrid, f, restore_max)
%
% INPUTS:
%       kgrid       - grid structure returned by makeGrid.m
%       f           - spatial distribution to smooth
%
% OPTIONAL INPUTS:
%       restore_max - Boolean controlling whether the maximum value is
%                     restored after smoothing (default = false)
%
% OUTPUTS:
%       f_sm        - smoothed spatial distribution
%
% USES:
%       numDim.m
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 29th April 2009
%       last update - 4th 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, ifft, fft2, ifft2, fftn, ifftn, makeGrid 

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

% check optional inputs
if nargin < 3
    restore_max = false;
end

% check for release B.0.1 inputs of smooth(f, kgrid, restore_max)
if ~isstruct(kgrid)
    disp('WARNING smooth function inputs changed, please update usage');
    f_sm = smooth(f, kgrid, restore_max);
    return
end

% extract the number of dimensions
switch numDim(f)
    case 1

        % create the filter
        Hann_filt = (0.5 - 0.5*cos(2*pi*(0:kgrid.Nx-1)/(kgrid.Nx-1))).';
        
        % apply the filter
        f_sm = abs(ifft(fftshift(fftshift(fft(f)).*Hann_filt)));
        
    case 2
        
        % create the filter
        Hann_x = (0.5 - 0.5*cos(2*pi*(0:kgrid.Nx-1)/(kgrid.Nx-1)))'; 
        Hann_z = (0.5 - 0.5*cos(2*pi*(0:kgrid.Nz-1)/(kgrid.Nz-1)))'; 
        Hann_filt = (Hann_x*Hann_z').';
        
        % apply the filter
        f_sm = abs(ifft2(fftshift(fftshift(fft2(f)).*Hann_filt)));
        
    case 3
        
        % create the filter
        Hann_x = (0.5 - 0.5*cos(2*pi*(0:kgrid.Nx-1)/(kgrid.Nx-1)))'; 
        Hann_y = (0.5 - 0.5*cos(2*pi*(0:kgrid.Ny-1)/(kgrid.Ny-1)))'; 
        Hann_z = (0.5 - 0.5*cos(2*pi*(0:kgrid.Nz-1)/(kgrid.Nz-1)))'; 
        Hann_filt_2D = (Hann_z*Hann_y');
        Hann_filt = zeros(kgrid.Nz, kgrid.Nx, kgrid.Ny);
        for index = 1:kgrid.Nx
            Hann_filt(:, index, :) = Hann_filt_2D(:, :)*Hann_x(index);
        end        
        
        % apply the filter
        f_sm = abs(ifftn(fftshift(fftshift(fftn(f)).*Hann_filt)));
        
end

% restore magnitude if required
if restore_max
    f_sm = (max(f(:))/max(f_sm(:)))*f_sm;
end