PSF Demonstration of Gaussian transfer function reducing sidelobes

SYNTAX: psf;

by Chuck DiMarzio Northeastern University December 2009

Requres fftaxis, fftaxisshift

!! This file may be copied, used, or modified for educational and
!! research purposes provided that this header information is not
!! removed or altered, and provided that the book is cited in
!! publications, as DiMarzio, Charles A., Optics for Engineers,
!! CRC Press, Boca Raton, FL, 2011.
!! http://www.crcpress.com
!! Other distribution is prohibited without permission.

Contents

Setup

xax=fftaxisshift(0:1023);
[x,y]=meshgrid(xax,xax);
rsq=x.^2+y.^2;
fax=fftaxisshift(fftaxis(xax));
[fx,fy]=meshgrid(fax,fax);
fsq=fx.^2+fy.^2;
ftop=max(max(fx));fmax=ftop/sqrt(length(fx))*4;
fgau=fmax/2;  % 1/e frequency of Gaussian

xmax=length(xax)/2/sqrt(length(fx)); % For plotting

Compute the transfer functions

mask1=single(fsq<(fmax^2));  % Low Pass filter for coherent transfer
                             % function.
mask2=exp(-fsq/fgau^2).*mask1;
fig1=figure;imagesc(fax,fax,mask1);
axis image;colormap(flipud(gray));colorbar;
axis(fmax*2*[-1,1,-1,1]);
xlabel('f_x, Frequency, Pixels^{-1}');
ylabel('f_y, Frequency, Pixels^{-1}');
fig2=figure;imagesc(fax,fax,mask2);
axis image;colormap(flipud(gray));colorbar;
axis(fmax*2*[-1,1,-1,1]);
xlabel('f_x, Frequency, Pixels^{-1}');
ylabel('f_y, Frequency, Pixels^{-1}');

Generate the object in the field plane and its FT in the pupil plane

edge=single(fx>0);
fig3=figure;imagesc(xax,xax,edge);
axis image;colormap(flipud(gray));colorbar;
axis(xmax*4*[-1,1,-1,1]);
xlabel('x, Position, Pixels');
xlabel('y, Position, Pixels');
edgep=fftshift(ifft2(fftshift(edge)));  % FT of the edge function
edgep1=edgep.*mask1;  % Multiply by the pupil functions
edgep2=edgep.*mask2;

Compute the images

edge1=fftshift(fft2(fftshift(edgep1)));
edge2=fftshift(fft2(fftshift(edgep2)));

Show results

fig4=figure;imagesc(xax,xax,real(edge1));
axis image;colormap(flipud(gray));colorbar;
axis(xmax*4*[-1,1,-1,1]);
xlabel('x, Position, Pixels');
xlabel('y, Position, Pixels');

fig5=figure;imagesc(xax,xax,real(edge2));
axis image;colormap(flipud(gray));colorbar;
axis(xmax*4*[-1,1,-1,1]);
xlabel('x, Position, Pixels');
xlabel('y, Position, Pixels');

Plot a slice

fig6=figure;plot(xax,real(edge2(floor(length(fax)/2)+1,:)),...
                 xax,real(edge1(floor(length(fax)/2)+1,:)),...
                 xax,real(edge(floor(length(fax)/2)+1,:)));
grid on;
xlabel('x, Position, Pixels');
ylabel('Field');
moose=axis;axis([xmax*4*[-1,1],moose([3,4])]);

Show the point--spread functions

psf1=fftshift(fft2(fftshift(mask1)));
fig7=figure;imagesc(xax,xax,real(psf1));
axis image;colormap(flipud(gray));colorbar;
axis(xmax*4*[-1,1,-1,1]);
xlabel('x, Position, Pixels');
xlabel('y, Position, Pixels');

psf2=fftshift(fft2(fftshift(mask2)));
fig8=figure;imagesc(xax,xax,real(psf2));
axis image;colormap(flipud(gray));colorbar;
axis(xmax*4*[-1,1,-1,1]);
xlabel('x, Position, Pixels');
xlabel('y, Position, Pixels');