function [cosGabor] = make2DGabor(N,k0,k1,sigma)
%
% This function returns a 2D sine and cosine Gabors with standard 
% deviation sigma and center frequency (k0,k1) cycles per N samples.   
% The Gabors are fftshifted so that their origin is in the middle of the
% 2D matrix.

%  NOTE: The Fourier transform of a Gaussian with standard deviation sigma
%  is a Gaussian with standard deviation (N / (2*pi*sigma)).

sigma_k = N / (2*pi * sigma);

%  NOTE: The bandwidth of a filter is sometimes defined by half height.  
%  However, when the filter has a Gaussian shape in the frequency domain  
%  (centered at k2=sqrt(k0*k0+k1*k1)), one often defines bandwidth as follows:

% octavebandwidth = log2( (k2 + sigma_k) / (k2 - sigma_k) );
% display(['octave bandwidth is ',num2str(octavebandwidth)])

%  NOTE:  we could have defined the function differently by specifying
%  its octave bandwidth and then computing k0 via:
%
%  k2 = N/(2*pi*sigma) * (power(2,octavebandwidth) + 1)/(power(2,octavebandwidth) -1);
%
%  ????????????????????????????????????????

cos2D = fftshift(mk2Dcosine(N,k0,k1));

for x = 1:N
    %  Define the Gabor to be centered at (N/2+1).   To filter an image
    %  with this Gabor, you will need to shift it.
    shiftx = ((1:N) - (N/2 + 1));
    Gaussian = 1/(sqrt(2*pi)*sigma) * exp(- shiftx.*shiftx/ (2 * sigma*sigma) );
    Gaussian2D = Gaussian' * Gaussian;
end
cosGabor = Gaussian2D .* cos2D;