function f_sm = smooth(f, kgrid, restore_max)
%SMOOTH     Smooth a matrix.
%
% DESCRIPTION:
%       Smooth filters an input matrix using an n-dimensional
%       frequency-domain Hanning window.
%
% USAGE:
%       f_sm = smooth(f, kgrid)
%       f_sm = smooth(f, kgrid, restore_max)
%
% INPUTS:
%       f           - spatial distribution to smooth
%       kgrid       - grid structure returned by makeGrid.m
%
% 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 - 17th July 2009
%       
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
%
% See also fft, ifft, fft2, ifft2, fftn, ifftn, makeGrid 

% check optional inputs
if nargin < 3
    restore_max = false;
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