FFT Calculation of Point Spread Function

Airy pattern (ideal point spread function)

FFT	exact

N = 128;
k = 8;
pupil = make_pupil(N,2*k);    % construct the pupil function
rmax = N/(4*k);
w = zeros(N,N);               % perfect wavefront
z = psf(pupil,w);             % generate the spread function
ncen = 1+N/2;
zp = z(ncen:ncen+30,ncen);    % take the central portion
xp = (0:30)/(2*k);            % generate x-coordinates

% calculate exact function
xe = linspace(0,30/(2*k),500);
ze = somb(2*xe);
ze = ze.*ze;

% plot the results
plot(xe,ze,xp,zp,'o');
xlabel('radius/\epsilon_o');
ylabel('irradiance');

% generate and save the images.
% Use logrithmic transformation for output image
imshow(log_image(z,5));
imwrite(255*log_image(z,5),gray(255),'airy.bmp');
imwrite(255*pupil,gray(255),'pupil.bmp');

Coma pattern (5 waves)

N = 256;
k = 2;
[pupil,x, rsq] = make_pupil(N,2*k);
rmax = N/(4*k);
w = (rsq-1.5) .* x';
z = psf(pupil,5.0*w);
out = 10*z;
idx = find(out>1.0);
out(idx)=1.0;
imshow(out);
imwrite(255*out,gray(255),'coma.bmp');
imwrite(255*pupil,gray(255),'pupil2.bmp');

Pupil construction

function [p, x, rsq] = make_pupil(N,d)
%
% MAKE_PUPIL [p, x, rsq] = make_pupil(N,d)
%
%	r is the FFT cell size (in normalized pupil coordinates)
%  N is the number of points in the cell (should be power of 2)

r = d/2;
x = linspace(-r,r,N+1);
x = ones(N,1)*x(1:N);
xsq = x.*x;
rsq = xsq+xsq';
p = rsq<=1.0;

PSF calculation

The final image is scaled to 1

function z = psf(pupil,w)


z = pupil.*exp(2*pi*j*w);
z = fftshift(z);
z = fft2(z);
z = fftshift(z);
z = z.*conj(z);
zmax = max(max(z));
z = z/zmax;

Logarithmic Image Transformation

function  out = log_image(in, decades)
%
%  out = log_image(in, decades)
%

out = 1+log10(in)/decades;
out = out.*(out>0.0);

Somb function

function z=somb(x)
%SOMB 2*j1(pi*x)/(pi*x) function.
%   SOMB(X) returns a matrix whose elements are the somb of the elements 
%   of X, i.e.
%        y = 2*j1(pi*x)/(pi*x)    if x ~= 0
%          = 1                    if x == 0
%   where x is an element of the input matrix and y is the resultant
%   output element.  

%   Author(s): J. Loomis, 6-29-1999

z=ones(size(x));
i=find(x);
z(i)=2.0*besselj(1,pi*x(i))./(pi*x(i));

Maintained by John Loomis, last updated 29 June 1999