function [y,h]=gausslp(x,Sigma,nstd)
%GAUSSLP  Gaussian lowpass filter
%	y=gausslp(x,Sigma,nstd) applies a Gaussian lowpass filter with
%	passband width determined by Sigma. If Sigma is of length two, the
%	first value is applied vertically and the second horizontally.
%	[y,h]=gausslp(...) return also the 2D frequency responce of the
%	filter.
%
%	The Farrokhnia approach employs Sigma=1/(2*sqrt(2)*CtrFreq) if the
%	corresponding Gabor filter has a center frequency CtrFreq.

if nargin<3,
  nstd=2;
end

% Generate the filter unit pulse responses
limit=nstd*Sigma;
% Vertical filter
i=-limit(1):limit(1);
h_v=1/(sqrt(2*pi)*Sigma(1)) * exp(-1/2*(i.^2)/Sigma(1)^2)';
% Horizontal filter
if length(Sigma)>1
  i=-limit(2):limit(2);
  h_h=1/(sqrt(2*pi)*Sigma(2)) * exp(-1/2*(i.^2)/Sigma(2)^2);
else
  h_h=h_v';
end

% Normalize the filters
%% h_h=h_h/sum(abs(h_h));
%% h_v=h_v/sum(abs(h_v));

% Filter the image
y=filter2(h_h,filter2(h_v,x));

h=h_v * h_h;










