PSFSHIFT Fix problem with phase direction in Figure 1
Coherent and Incoherent Transfer Functions and Point-Spread Functions for a tilt in the pupil plane
by Chuck DiMarzio Northeastern University December 2009
!! 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
fcenter=1/128; % Center frequency for bar chart in cycles per pixel tilt=2*pi/4*(1/fcenter); % Wavefront tilt 1/4 cycle per fcenter 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)); xmax=length(xax)/2/sqrt(length(fx));
Compute the transfer functions and point--spread functions
theshift=floor(1/fcenter/4); % We want to shift a quarter cycle. psf1=(x==theshift)&(y==0); mask1=fftshift(ifft2(fftshift(psf1))); % Go to pupil and generate the mask % for the displacement mask1=mask1.*(fsq<fmax^2); % Make the aperture psf1=fftshift(ifft2(ifftshift(mask1))); % Recompute the actual psf1 psf1a=abs(psf1).^2; mask1a=fftshift(ifft2(fftshift(psf1a))); % incoherent OTF
Show a slice of the phase through the Transfer Functions
a1=(abs(mask1)>1e-6*max(max(abs(mask1)))).*angle(mask1); a1a=(abs(mask1a)>1e-6*max(max(abs(mask1a)))).*(angle(mask1a)); fig1=figure;plot(fax,unwrap(a1(1+floor(length(fx)/2),:)),... fax,unwrap(a1a(1+floor(length(fx)/2),:))); grid on; % why does the phase go in the opposite direction? moose=axis;axis([fmax*6*[-1,1],moose([3,4])]); xlabel('f_y, Spatial Freq., pixel^{-1}'); ylabel('PTF(f_x,0)');

Show a slice of the amplitudes through the Transfer Functions
fig2=figure;plot(fax,abs(mask1(1+floor(length(fx)/2),:)/max(max(mask1))),... fax,abs(mask1a(:,1+floor(length(fx)/2)))/max(max(mask1a))); grid on; moose=axis;axis([fmax*4*[-1,1],moose([3,4])]); xlabel('f_x, Spatial Freq., pixel^{-1}'); ylabel('MTF(f_x,0)');

Show the PTF with imagesc
fig3=figure;imagesc(fax,fax,a1);axis image;colormap(flipud(gray)); % CPTF axis(fmax*4*[-1,1,-1,1]); xlabel('f_x, pixel^{-1}'); ylabel('f_y, pixel^{-1}'); colorbar; fig4=figure;imagesc(fax,fax,a1a);axis image;colormap(flipud(gray)); % IPTF axis(fmax*4*[-1,1,-1,1]); xlabel('f_x, pixel^{-1}'); ylabel('f_y, pixel^{-1}'); colorbar;


Plot slices of the point--spread functions
fig5=figure;plot(xax,real(psf1(1+floor(length(fax)/2),:)/max(max(psf1))),... xax,real(psf1a(1+floor(length(fax)/2),:))/max(max(psf1a))); moose=axis;axis([xmax*8*[-1,1],moose([3,4])]); grid on; xlabel('x, Position, Pixels'); ylabel('PSF(x,0)');

Now do the imaging: Coherent object
Show object on top and scaled image on bottom for comparison
cobject=cos(2*pi*fcenter*x); cobjectf=fftshift(ifft2(fftshift(cobject))); % Go to pupil cimagef=cobjectf.*mask1; cimage=fftshift(ifft2(ifftshift(cimagef))); fig6=figure;imagesc(xax,xax,[real(cobject(1:length(fax)/2,:));... real(cimage(length(fax)/2+2:end,:))/max(max(abs(cimage)))]); colormap(flipud(gray));colorbar; axis image; xlabel('x, Position, Pixels'); ylabel('x, Position, Pixels');

Incoherent object
iobject=abs(cobject).^2; iimage=abs(cimage).^2; fig7=figure;imagesc(xax,xax,[real(iobject(1:length(fax)/2,:));... real(iimage(length(fax)/2+2:end,:))/max(max(abs(iimage)))]); colormap(flipud(gray));colorbar; axis image; xlabel('x, Position, Pixels'); ylabel('x, Position, Pixels');
