The 2D Ambiguity Function (AF) and its relation to 1D Optical Transfer Function (OTF)
The Ambiguity Function (AF) is an useful tool for optical system analysis. This post is a basic introduction to AF, and how it can be useful for analyzing incoherent optical systems. We will see that the AF simultaneously contains all the OTFs associated with an rectangularly separable incoherent optical system with varying degree of defocus [24]. Thus by inspecting the AF of an optical system, one can easily predict the performance of the system in the presence of defocus. It has been used in the design of extendeddepthoffield cubic phase mask system.
NOTE:
This post was created using an IPython notebook. The most recent version of the IPython notebook can be found here.
To understand the basic theory, we shall consider a onedimensional pupil function, which is defined as:
The *generalized pupil function* associated with is the complex function given by the expression [1]:
where is the aberration function. Then, the amplitude PSF of an aberrated optical system is the Fraunhofer diffraction pattern (Fourier transform with the frequency variable equal to ) of the generalized pupil function, and the intensity PSF is the squared magnitude of the amplitude PSF [1]. Note that is the distance between the diffraction pattern/screen and the aperture/pupil.
For the purpose of analyzing optical systems using the ambiguity function, it is convenient to separate the defocus term ( is the wavefront focus error coefficient [3]) from the rest of the aberration terms in the aberration function :
Since the amplitude PSF is the Fourier transform of the pupil function (see above), and the amplitude transfer function (ATF) is the Fourier transform of the amplitude PSF, the ATF is proportional to a scaled pupil function . Specifically, the ATF , is:
The (one dimensional) optical transfer function (OTF) is defined as the normalized autocorrelation of the ATF:
writing in place of , we have
The simplified equation is written as follows:
The ambiguity function, which is used for radar analysis, is defined as the Fourier transform of the product :
Comparing the simplified expression of the aberrated OTF, , and the ambiguity function, , we see that:
This implies that the ambiguity function associated with a base function contains the OTF along the line (which has a slope , and passes through the origin in the plane as shown in the following figure.
Note:
1. The “base function” itself does not contain the defocus term, although it may contain other aberrations (see equation (3)).
2. Equation (6) is the “normalized autocorrelation” of the generalized pupil function, which implies that the maximum value of the function is 1. In equation (7) there is no explicit normalization and so the maximum value can be greater or less than 1. We must be mindful of this when comparing the AF to the OTF and prefer to chose appropriate amplitude for the base function (as shown below).
In fact, the whole 2D AF contains the OTFs for various values of defocusing arranged in a polar fashion. Of special interest is the infocus OTF for which and hence corresponds to the line along the “” axis. i.e.
Example: Ambiguity function of a diffraction limited 1D pupil.
If the only aberration is that of defocus, in the expression for the generalized pupil function. Therefore the basefunction which we consider is
First, we will consider the evaluation of the following integral of the general form:
where, the function is defined as:
The following figure, which shows the regions of overlap between and (with ), is useful to understand the finite limites of the integral.
When :
where, the function used above is the normalized that is defined as
When :
Note that we really didn’t have to break the integral into two parts. Instead we could just evaluate the integral between the limits and .
We can combine the two cases and write the equation more compactly as:
In our example, for which the base function is , the parameters , (so that the maximum value is 1) and . We can now write:
Now, from equation (8) we have . Therefore we can write the expression for the OTF directly from the AF.
Note:
In [2], the authors have used the *unnormalized* function that is defined as . Hence, there is an “extra” within the function in their paper. We use the normalized definition as it is used in [1], and in Numpy.
The figure below (after the code block) shows the function for . It can be seen in figure that roots of the (normalized) function occurs at . Therefore the zero value loci for the AF associated with the onedimensional rectangular pupil are found by equating .
from __future__ import print_function, division import numpy as np import scipy as sp import matplotlib.cm as cm import matplotlib.pyplot as plt from IPython.display import Image %matplotlib inline
x = np.linspace(8, 8, 150) y = np.sinc(x) fig, ax = plt.subplots(1,1) ax.plot(x, y, label='$sinc(x)$') ax.set_ylim(0.3, 1.02) ax.set_xlim(8, 8) ax.set_xlabel('x') ax.set_title('sinc(x)', y=1.02) rootsx = range(8, 0) + range(1,9) ax.scatter(rootsx, np.zeros(len(rootsx)), c='r', zorder=20) ax.grid() plt.show()
The following figures plot the AF of the one dimensional rectangular pupil and variation of OTF with focus error . In each figure, the plot on the left shows the ambiguity function and several lines corresponding to different defocus error . The zero value locai, are denoted by the gray dashed lines. The intersection of zero value locai with the line(s) , which in the OTF plot represents the frequencies where the OTF is zero, is shown by the cross (‘x’) markers. The plots on the right show the corresponding variation of the OTF with different.
def plot_1dRectAF(w20LambdaBy2, N=15, umin=2, umax=2, ymin=8, ymax=8): """rudimentary function to show the AF Parameters  w20LambdaBy2 : list of real values specifies the amount defocus error W_{20} in terms of lambda/2. The slope of the OTF associated with W_{20} in AF plane is 2*W20/lambda. N : integer number of zero value loci to plot """ u = np.linspace(umin, umax, 200) y = np.linspace(ymin, ymax, 400) uu, yy = np.meshgrid(u, y) # grid # Numpy's normalized sinc function = sin(pi*x)/(pi*x) af = (1  np.abs(uu)/2)*np.sinc(yy*(2  np.abs(uu))) # plot fig = plt.figure(figsize=(12, 7)) ax1 = fig.add_axes([0.12, 0, 0.42, 1.0]) # [*left*, *bottom*, *width*,*height*] ax2 = fig.add_axes([0.6, 0.12, 0.38, 0.76]) im = ax1.imshow(af, cmap=cm.bwr, origin='lower', extent=[umin, umax, ymin, ymax], vmin=1.0, vmax=1.0, aspect=2./6) plt.colorbar(im, ax=ax1, shrink=0.77, aspect=35) # zero value loci for n in range(1, N+1): zvl = n/(2.0  abs(u[1:1])) ax1.plot(u[1:1], zvl, color='#888888', linestyle='dashed') ax1.plot(u[1:1], zvl, color='#888888', linestyle='dashed') # OTF line on AF plane for elem in w20LambdaBy2: otfY = elem*u # OTF line in AF with slope 2w_{20}/lambda ax1.plot(u, otfY) # intersections def get_intersections(b): # b is tan(phi) or 2w_{20}/lambda n = np.linspace(1, np.floor(b), np.floor(b)) u1 = 1 + np.sqrt(1  n/b) u2 = 1  np.sqrt(1  n/b) y1 = u1*b y2 = u2*b u = np.hstack((u1, u2)) y = np.hstack((y1, y2)) return u, y for elem in w20LambdaBy2: intersectionsU, intersectionsY = get_intersections(elem) ax1.scatter(intersectionsU, intersectionsY, marker='x', c='k', zorder=20) # OTF plots for elem in w20LambdaBy2: otf = (1  np.abs(u)/2)*np.sinc(elem*u*(2  np.abs(u))) ax2.plot(u, otf, label='$W_{20}' + '= {}lambda/2$'.format(elem)) # axis settings ax1.set_xlim(umin, umax) ax1.set_ylim(ymin, ymax) ax1.set_title('2D AF of 1D rect pupil P(x)', y=1.01) ax1.set_xlabel('u', fontsize=14) ax1.set_ylabel('y', fontsize=14) ax2.axhline(y=0, xmin=2, xmax=2, color='#888888', zorder=0, linestyle='dashed') ax2.grid(axis='x') ax2.legend(fontsize=12) ax2.set_xlim(2, 2); ax2.set_ylim(0.2, 1.005) ax2.set_title("Optical Transfer Function", y=1.02) ax2.set_xlabel("u (scaled saptial frequency)", fontsize=14) #fig.tight_layout() plt.show()
Plots of AF and OTF for which for , i.e. . The equation for the OTF was obtained directly from the ambiguity function by replacing with .
The points of intersection between lines and each zero loci curve (for ) can be found by equating , which gives
For example, the line intersects the zero loci curve at , the line intersects the curve at , and the curve at , and so on. As we can see, there are odd number of intersections of the lines with the zero loci curves when the defocus error is an integer multiple of . This corresponds to an odd number of zerovalue OTF points in the OTF plots when the defocus error is an integer multiple of .
Note that we only found the abscissa ( coordinates) for the intersection points above.
plot_1dRectAF(w20LambdaBy2=[0, 1, 2, 3])
Plots of AF and OTF for which for , i.e. (the particular values were arbitrarily chosen, and the line is plotted for comparison). In this case we find that there are even number of intersections for each line with the zeroloci curves in the AF plot. Correspondingly there are even number of zerovalue OTF points for each OTF with focus error for .
plot_1dRectAF(w20LambdaBy2=[0, 0.5, 0.99, 1.5, 3.6])
Also, note from the above plots that the first occurrence of the OTF intersecting/crossing the zero value happens as . This point is also the known as the traditional or Hopkins criterion for misfocus [4].
Example: Ambiguity function of cubic phase mask
Cubic phase masks (CPM) has been used to make hybrid optical systems largely invariant to defocus, thus extending the depth of field of such systems. For details on CPM systems, refer to [4].
The expression for the ambiguity function of a cubic phase mask is given below [4]:
We will numerically evaluation the above expression for the AF of cubic phase mask (CPM), recognizing that the expression is of the form of inverse Fourier Transform of where,
Here are the steps for rendering the 2D AF for CPM:
 Create a vector

Create a vector (i.e. between , which is the maximum region of integration)

For every “u” in create the sequence where, using the following rule:
a. if , then evaluate where
b. if , then evaluate where

Take inverse Fourier Transform of each sequence (with as the transform variable). i.e. IFFT
The expression for the OTF for the CPM is then given by:
We will use numerical integration to generate the plots of OTFs for the CPM with various amounts of defocus.
from scipy.integrate import quad import warnings warnings.simplefilter(action='error', category=np.ComplexWarning) # Turn on the warning to ensure that the numerical integration is "reliable"? warnings.simplefilter(action='always', category=sp.integrate.IntegrationWarning) ifft = np.fft.fft fftshift = np.fft.fftshift fftfreq = np.fft.fftfreq # cubic phase mask parameters alpha = 90 gamma = 3 umin, umax = 2, 2 ymin, ymax = 60, 60 w20LambdaBy2 = [0, 5, 15] # amounts of defocus in units of wavelength (by 2) uVec = np.linspace(umin, umax, 300) N = 512 # number of samples along "t" ... and for FFT L = 1 def gut(t, alpha, gamma, u): return 0.5*np.exp(1j*alpha*((t + u/2)**gamma  (t  u/2)**gamma)) guy = np.empty((N, len(uVec))) #roi = np.empty((N, len(uVec))) # for debugging & seeing the region of integration t = np.linspace(2*L, 2*L, N) dt = (4*L)/(N1) for i, u in enumerate(uVec): g = np.zeros_like(t, dtype='complex64') if 1 <= u/2.0 < 0: mask = (t > (u/2  1))*(t < (u/2 + 1)) #roi[:, i] = mask.astype('float32') g[mask] = gut(t[mask], alpha, gamma, u) guy[:, i] = np.abs(fftshift(ifft(g))) elif 0 <= u/2.0 <= 1: mask = (t > (u/2  1))*(t < (u/2 + 1)) #roi[:, i] = mask.astype('float32') g[mask] = gut(t[mask], alpha, gamma, u) guy[:, i] = np.abs(fftshift(ifft(g))) # Normalize to make maximum value = 1 guyMax = np.max(np.abs(guy.flat)) guy = guy/guyMax yindex = fftshift(fftfreq(N, dt)) ymin, ymax = yindex[0], yindex[1] fig = plt.figure(figsize=(12, 7)) ax1 = fig.add_axes([0.12, 0, 0.5, 1.0]) # [*left*, *bottom*, *width*,*height*] ax2 = fig.add_axes([0.66, 0.23, 0.32, 0.54]) im = ax1.imshow(guy**0.8, cmap=cm.YlGnBu_r, origin='lower', extent=[umin, umax, ymin, ymax], vmin=0.0, vmax=1.0,aspect=1./40) plt.colorbar(im, ax=ax1, shrink=0.55, aspect=35) # OTF line in AF for elem in w20LambdaBy2: otfY = elem*uVec # OTF line in AF with slope 2w_{20}/lambda ax1.plot(uVec, otfY, alpha = 0.6, linestyle='solid') ax1.set_xlim(umin, umax) ax1.set_ylim(ymin, ymax) ax1.set_xlabel('u', fontsize=14) ax1.set_ylabel('y', fontsize=14) ax1.set_title('2D AF of 1D cpm', y=1.01) # Magnitude plots of the OTF of the cpm def otf_cpm(t, alpha, gamma, u, w20LamBy2): return (0.5*np.exp(1j*alpha*((t + u/2)**gamma  (t  u/2)**gamma)) *np.exp(1j*2*np.pi*u*w20LamBy2*t)) def complex_quad(func, a, b, **kwargs): """Compute numerical integration of complex function between limits a and b Adapted from the following SO post: stackoverflow.com/questions/5965583/usescipyintegratequadtointegratecomplexnumbers """ def real_func(x, *args): if args: return sp.real(func(x, *args)) else: return sp.real(func(x)) def imag_func(x, *args): if args: return sp.imag(func(x, *args)) else: return sp.imag(func(x)) real_integral = quad(real_func, a, b, **kwargs) imag_integral = quad(imag_func, a, b, **kwargs) return (real_integral[0] + 1j*imag_integral[0], real_integral[1:], imag_integral[1:]) for elem in w20LambdaBy2: Huw = np.empty_like(uVec, dtype='complex64') for i, u in enumerate(uVec): if 1 <= u/2.0 < 0: Huw[i] = complex_quad(func=otf_cpm, a=u/2  1, b=u/2 + 1, args=(alpha, gamma, u, elem))[0] elif 0 <= u/2.0 <= 1: Huw[i] = complex_quad(func=otf_cpm, a=u/2  1, b= u/2 + 1, args=(alpha, gamma, u, elem))[0] HuwMax = np.max(np.abs(Huw)) ax2.plot(uVec, np.abs(Huw)/HuwMax, label='$W_{20}' + '= {}lambda/2$'.format(elem)) ax2.legend(fontsize=12) ax2.set_ylabel("Magnitude of OTF", fontsize=14) ax2.set_xlabel("Spatial frequency, u", fontsize=14) ax2.set_title('Optical Transfer Function of cpm', y=1.01) plt.show()
We can see from the above plots of the OTF, that the OTFs are insensitive to large amounts of defoucs. More importantly, there are no regions of zero values within the passband (though the bandwidth is not really constant for large amounts of defocus). This property is extremely important from the point of view of the reconstruction filter (such as an inverse filter or an Wiener filter).
References
 Introduction to Fourier Optics, Joseph Goodman
 The Ambiguity Function as a polar display of the OTF, K.H. Brenner, A.W. Lohmann and J. OjedaCastaneda, Optics Communication, 1982
 Misfocus tolerance seen by simple inspection of the ambiguity function, H. Bartelt, J. OjedaCastaneda, and Enrique Sicre, Applied Optics, 1984.
 Extended depth of field through wavefront coding, Edward R. Dowski, Jr., and W. Thomas Cathey, Applied Optics, 1995.
I received a comment from Health Zhou via email. He suggested that the line “rootsx =range(8,0)+range(1,9)” should be change to “rootsx =np.hstack((range(8,0),range(1,9)))”. I think, although it doesn’t change the code functionally, he has a fair point!.