% omega = normal orientation
% F = frequency
% Ka = Bandwidth // normal
% Kb = Bandwidth \perp normal
% szx = width of the output in pixel
% szy = height of the output in pixel
function g=spGabor(omega, F, Ka, Kb, szx, szy)
	%omega = 10;   	%Omega = 0, orientation
	%F= 15;	 	%peak frequency cycle/degree

	% Gabor Half Magnitude Range
	% Freq: [F-F*Ka F+F*Ka];
	% Orientation bandwidth: [omega-tan-1(Kb) omega+tan-1(Kb);]
    C =sqrt(log(2)/pi);
    a = Ka*F/C;
	b = Kb*F/C;
	f = F*[cosd(omega); sind(omega)];

	% rotation, scale
	mscale=diag([a b]);
	mrot = [cosd(omega) sind(omega); -sind(omega) cosd(omega)];
	mtran =mscale*mrot; 

	%rasterize
	nDeg= 4; %FoV measured in degree, here we assume face is about 4 deg viewed in regular distance.
	%dx = 1/24;
	%dy = 1/24;
	dx = nDeg/szx;
	dy = nDeg/szy;

	xtick = [(-nDeg/2+dx):dx:nDeg/2];
	ytick = [(-nDeg/2+dx):dy:nDeg/2];

	frot = inv(mtran)'*f;
	c = exp(-pi * frot' * frot);		%zero mean constant
	for ix = 1:length(xtick);
		for iy =1:length(ytick)
			pos = [xtick(ix);ytick(iy)];
			h(iy, ix) = exp( j * 2 * pi * ( pos' * f )) - c;
			
			xrot = mtran*pos;
			w(iy, ix) = exp(-pi * xrot' * xrot);
		end
	end

	g =  h.*w;
	g = g/norm(g);
	return
