fft_upsample( in_m,upsample_factor )

function out_m = fft_upsample( in_m,upsample_factor )
%
% fft_upsample  -   upsamples a given input "in_m" by an INTEGER factor using the Fourier Domain
%
% Format:   out_m = fft_upsample( in_m,upsample_factor )
%
% Input:    in_m            -   a given distribution to be upsampled
%           upsample_factor -   upsampling factor (integer, rounding towards zero if not) 
%
% Output:   out_m           -   the given distribution upsampled by factor: floor(upsample_factor)
%


% 
% Theory:   the upsampling is done by zero-padding in the input domain BETWEEN the samples,
%           then at the fourier domain, taking the single spectrum (out of the repetition of spectrums)
%           i.e. low pass with Fcutoff=PI/upsample_factor, zeroing the rest of the spectrum
%           and then doing ifft to the distribution.
%           since we have a zero padding operation in time, we need to multiply by the fourier gain.
%
%              +-----------+     +-------+     +---------+     +--------+     +--------+
%   y[n,m] --> | up-sample | --> |  FFT  | --> |   LPF   | --> | * Gain | --> |  IFFT  | --> interpolated 
%              | factor M  |     +-------+     | Fc=PI/M |     +--------+     +--------+
%              +-----------+                   +---------+
%
%           this operation is the same as the following one (which has less operations):
%
%              +-------+     +--------+     +--------------+     +--------+
%   y[n,m] --> |  FFT  | --> | * Gain | --> | Zero Padding | --> |  IFFT  | --> interpolated 
%              +-------+     +--------+     +--------------+     +--------+
%
%           NOTE THAT, the zero-padding must be such that the D.C. ferquency remains the D.C. frequency
%           and that the zero padding is applied to both positive and negative frequencies. 
%           The zero padding actually condences the frequency -> which yields a longer series in the 
%           image domain, but without any additional information, thus the operation must be an interpolation.
%
%           We would like to interpolate the data BETWEEN the samples of the given distribution
%           therefore, the output array should be of size: (input_array_size-1)*up_sample_factor + 1
%


% do rounding
upsample_factor = floor( upsample_factor + 10*eps );

% calculate target sizes
initial_x_size  = size(in_m,2);
initial_y_size  = size(in_m,1);
target_x_size   = initial_x_size*upsample_factor;
target_y_size   = initial_y_size*upsample_factor;
fourier_gain    = (upsample_factor*upsample_factor);

% calculate user sizes
center_x    = floor(target_x_size/2);
center_y    = floor(target_y_size/2);
span_x      = [1:initial_x_size]-ceil(initial_x_size/2);
span_y      = [1:initial_y_size]-ceil(initial_y_size/2);

% transform to fourier domain
spectrum_m = fftshift( fft2( in_m )*fourier_gain );

% zero padding (=condensing the frequency domain)
padded_specturm_m = zeros( target_y_size,target_x_size );
padded_specturm_m( center_y+span_y+1,center_x+span_x+1 ) = spectrum_m; 

% inverse transform to reveal the interpolated series in space domain. 
out_m = ifft2( ifftshift( padded_specturm_m ) );

% we want values which are ONLY between the given points (interpolation and not extrapolation)
out_m = out_m(1:end-upsample_factor+1,1:end-upsample_factor+1);