MATLAB Answers


Modeling the point spread function of a simple lens with defocus aberration

Asked by eor
on 28 Jun 2013
Latest activity Commented on by Vanessa
on 8 Jun 2017

I am trying to model the PSF of a simple lens with a defocus aberration. I started by finding the geometric optics prediction of the system blur and convolving it with the impulse response of the exit pupil to account for diffraction. This gave me physically believable results. I would like to verify that with a purely wave-optical approach. I have attempted to do this following Goodman's formulation of an aberrated spherical wave front in chapter 6 of Introduction to Fourier Optics:

P(x,y) = P(x,y)*exp[j*k*W(x,y)]

where P(x,y) is the aperture and W(x,y) = (-1/2)*((1/za)-(1/zi))*(x^2 + y^2). zi is the distance between the lens and the image plane. za is the distance between the lens and the point to which the aberrated wave is converging. As the difference between zi and za gets larger, the PSF should get wider (more blur). But in the code that follows, I am seeing no change as a function of za.

% Camera system parameters
zi = .0981;                 % Distance to image plane (m)
Feff = .078;                % Effective focal length of the system (m)
w = .001;                   % Radius of the aperture (m)
% Array paramters
M = 1280;                   % number of pixels wide
N = 1040;                   % number of pixels long
dx = 6e-6;                  % sample interval (m)
Lx = dx*M;                  % physical width (m)
Ly = dx*N;                  % physical length (m)
% Create array
x = -Lx/2:dx:Lx/2-dx;
y = -Ly/2:dx:Ly/2-dx;
[xx,yy] = meshgrid(x,y);
% Illumination
lambda = 0.4e-6;            % wavelength (m)
k = 2*pi/lambda;            % wavenumber (1/m)
% Center of focus of aberrated wave (this is what is changing)
za = .099;
% wave front aberration
kW =  -0.5*((1/za) - (1/zi))*(xx.*xx + yy.*yy);
% complex pupil function
P = circ(sqrt(xx.*xx + yy.*yy)/(2*w)) .* exp(1i*k*kW);
% transfer function
h = fftshift(fft2(fftshift(P))) * dx^2;
im = abs(h).^2;

"Circ" creates a binary mask representing the pupil. I have tried to evaluate the function in frequency space (an OTF approach) with identical results. I have also tried increasing the sampling rate with no change. I am aware of solutions using Zernike or Siedel coefficents, but would like to keep this as simple as possible for illustrative purposes. Any suggestions on the code or my logic? I could be totally off-base with any or all of it.

Apologies for the lengthy question. Your time is much appreciated.

  1 Comment


I am also trying to estimate the PSF of a lens and I wanted to ask if you could provide me with the 'circ' code you used, as I am having trouble getting results.

Kind regards,


Log in to comment.


No products are associated with this question.

0 Answers

Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

MATLAB Academy

New to MATLAB?

Learn MATLAB today!