function h = compute_psf(N, sigma)

[x, y] = meshgrid(1:N);
h = exp( -( (x-1-N/2).^2 + (y-1-N/2).^2)/(2*sigma*sigma) );

% normalize to unit area (DC gain = 1)
area = sum(sum(h));
h = h / area;