References
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
The Ambiguity Function (AF) is an useful tool for optical system analysis. As we will see, the AF simultaneously contains all the OTFs associated with an rectangularly separable incoherent optical system with varying degree of defocus [2-4]. 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 extended-depth-of-field cubic phase mask system.
To understand the basic theory, we shall consider a one-dimensional pupil function, which is defined as:
$$ (1) \hspace{40pt} P(x) = \left \{ \begin{array}{cc} 1 & \,\, \text{if} \, & |x| \leq 1, \\ 0 & \,\, \text{if} \, & |x| > 1, \end{array}\right .\ $$The generalized pupil function associated with $P(x)$ is the complex function $\mathcal{P}(x)$ given by the expression [1]:
$$ (2) \hspace{40pt} \mathcal{P}(x) = P(x)e^{jkW(x)} $$where $W(x)$ is the aberration function. Then, the amplitude PSF of an aberrated optical system is the Fraunhofer diffraction pattern (Fourier transform with the frequency variable $f_x$ equal to $x/\lambda z_i$) of the generalized pupil function, and the intensity PSF is the squared magnitude of the amplitude PSF [1]. Note that $z_i$ 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 ($W_{20}$ is the wavefront focus error coefficient [3]) from the rest of the aberration terms in the aberration function $W(x)$:
$$ (3) \hspace{40pt} \begin{array}{rl} \mathcal{P}(x) & = & P(x)e^{jkW(x)} \\ & = & P(x)e^{jk(\tilde{W}(x) + W_{20}x^2)} \\ & = & \mathcal{P}_o(x)e^{jkW_{20}x^2} \end{array} $$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 $\mathcal{P}$. Specifically, the ATF $H(f_x)$, is:
$$ (4) \hspace{40pt} H(f_x) = \mathcal{P}(\lambda z_i f_x) $$The (one dimensional) optical transfer function (OTF) is defined as the normalized autocorrelation of the ATF:
$$ (5) \hspace{40pt} \mathcal{H}(f_x) = \frac{\int\limits_{-\infty}^{\infty} H \left(p + \frac{f_x}{2}\right) H^* \left(p - \frac{f_x}{2} \right) dp}{\int\limits_{-\infty}^{\infty} |H(p)|^2 dp} $$$$ \hspace{50pt} \mathcal{H}(f_x) = \frac{\int\limits_{-\infty}^{\infty} \mathcal{P} \left(q + \frac{\lambda z_i f_x}{2}\right) \mathcal{P}^* \left(q - \frac{\lambda z_i f_x}{2} \right) dq} {\int\limits_{-\infty}^{\infty} |\mathcal{P}(q)|^2 dq} $$writing $u$ in place of $\lambda z_i f_x$, we have
$$ \hspace{50pt} \mathcal{H}(u) = \frac{\int\limits_{-\infty}^{\infty} \mathcal{P} \left(q + \frac{u}{2}\right) \mathcal{P}^* \left(q - \frac{u}{2} \right) dq} {\int\limits_{-\infty}^{\infty} |\mathcal{P}(q)|^2 dq} $$$$ \hspace{50pt} \mathcal{H}(u) = \frac{\int\limits_{-\infty}^{\infty} \mathcal{P}_o\left(q + \frac{u}{2}\right) e^{jkW_{20}(q + u/2)^2} \mathcal{P}_o^* \left(q - \frac{u}{2} \right) e^{-jkW_{20}(q - u/2)^2} dq} {\int\limits_{-\infty}^{\infty} |\mathcal{P}_o(q)e^{jkW_{20}q^2}|^2 dq} $$$$ \hspace{50pt} \mathcal{H}(u) = \frac{\int\limits_{-\infty}^{\infty} \mathcal{P}_o\left(q + \frac{u}{2}\right) \mathcal{P}_o^* \left(q - \frac{u}{2} \right) e^{jkW_{20} [(q + u/2)^2 - (q - u/2)^2 ]} dq} {\int\limits_{-\infty}^{\infty} |P(q)|^2 dq} $$The simplified equation is written as follows:
$$ (6) \hspace{40pt} \mathcal{H}(u; W_{20}) = 0.5 \int\limits_{-\infty}^{\infty} \mathcal{P}_o\left(q + \frac{u}{2}\right) \mathcal{P}_o^* \left(q - \frac{u}{2} \right) e^{j 2\pi \left(\frac{2 W_{20}}{\lambda} \right)uq} dq $$The ambiguity function, which is used for radar analysis, is defined as the Fourier transform of the product $f(y + u/2) f^*(y - u/2)$:
$$ (7) \hspace{40pt} A(u, y) = \int\limits_{-\infty}^{\infty}f\left(t + \frac{u}{2} \right) f^*\left(t - \frac{u}{2} \right)e^{j2\pi yt} dt $$Comparing the simplified expression of the aberrated OTF, $\mathcal{H}(u; W_{20})$, and the ambiguity function, $A(u, y)$, we see that:
$$ (8) \hspace{40pt} \mathcal{H}(u; W_{20}) = A \left(u, 2 u \frac{W_{20}}{\lambda} \right) $$This implies that the ambiguity function $A(u, y)$ associated with a base function $\mathcal{P}_o(u)$ contains the OTF $\mathcal{H}(u; W_{20})$ along the line $y = \frac{2 W_{20}}{\lambda} u$ (which has a slope $\tan (\phi) =\frac{2 W_{20}}{\lambda}$, and passes through the origin in the $y-u$ plane as shown in the following figure.
Notes
The "base function" itself does not contain the defocus term, although it may contain other aberrations (see equation (3)).
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 $f(u)$ (as shown below).
Image('images/relationAFandOTF.png', embed=True, width=300)
In fact, the whole 2-D AF contains the OTFs $\mathcal{H}(u; W_{20})$ for various values of defocusing $W_{20}$ arranged in a polar fashion. Of special interest is the in-focus OTF for which $W_{20}=0$ and hence corresponds to the line along the "$u$" axis. i.e. $\mathcal{H}(u; 0) = A(u,0)$
If the only aberration is that of defocus, $\tilde{W}(x)=0$ in the expression for the generalized pupil function. Therefore the base-function which we consider is
$$ (9) \hspace{40pt} \mathcal{P}_o(x) = P(x) = \left \{ \begin{array}{cc} 1 & \,\, \text{if} \, & |x| \leq 1, \\ 0 & \,\, \text{if} \, & |x| > 1, \end{array}\right .\ $$First, we will consider the evaluation of the following integral of the general form:
$$ (10) \hspace{40pt} I(a, y) = \int\limits_{-\infty}^{\infty}f(t + a) f^*(t - a)e^{j2\pi yt} dt $$where, the function $f(t)$ is defined as:
$$ \hspace{50pt} f(t) = \left \{ \begin{array}{cc} A & \,\, \text{if} \, & |t| \leq L, \\ 0 & \,\, \text{if} \, & |t| > L, \end{array}\right .\ $$The following figure, which shows the regions of overlap between $f(t-a)$ and $f(t+a)$ (with $A=1$), is useful to understand the finite limits of the integral.
Image('images/ambiguityFnKernel.png', embed=True)
When $0 < a < L$:
$$ \begin{array}{cl} I(a, y) & = & \int\limits_{a-L}^{-a+L} A^2 e^{j2\pi yt} dt \\ & = & A^2 \left| \frac{e^{j2\pi yt}}{j2\pi y} \right|_{a-L}^{-a+L} \\ & = & A^2 \frac{1}{\pi y} \left[ \frac{e^{j2\pi (L-a)y} - e^{-j2\pi (L-a)y} }{2j} \right] \\ & = & A^2 \frac{1}{\pi y} \sin{( 2\pi (L-a)y)} \\ & = & 2(L-a)A^2 \frac{\sin{( 2\pi (L-a)y)}} {2\pi (L-a) y} \\ & = & 2(L-a)A^2 \text{sinc}[2(L-a)y] \end{array} $$where, the $\text{sinc}$ function used above is the normalized $\text{sinc}$ that is defined as $\text{sinc}(x)=sin(\pi x)/(\pi x)$
When $-L < a \leq 0$:
$$ \begin{array}{cl} I(a, y) & = & \int\limits_{-a-L}^{a+L} A^2 e^{j2\pi yt} dt \\ & = & A^2 \left| \frac{e^{j2\pi yt}}{j2\pi y} \right|_{-a-L}^{a+L} \\ & = & A^2 \frac{1}{\pi y} \left[ \frac{e^{j2\pi (L+a)y} - e^{-j2\pi (L+a)y} }{2j} \right] \\ & = & A^2 \frac{1}{\pi y} \sin{( 2\pi (L+a)y)} \\ & = & 2(L+a)A^2 \frac{\sin{( 2\pi (L+a)y)}} {2\pi (L+a) y} \\ & = & 2(L+a)A^2 \text{sinc}[2(L+a)y] \end{array} $$Note that we really didn't have to break the integral into two parts. Instead we could just evaluate the integral between the limits $-(L - |a|)$ and $(L - |a|)$.
We can combine the two cases and write the equation more compactly as:
$$ (11) \hspace{40pt} \begin{array}{ll} I(a, y) & = &\int\limits_{-\infty}^{\infty}f(t + a) f^*(t - a)e^{j2\pi yt} dt \\ & = & 2(L - |a|)A^2 \text{sinc}[2(L - |a|)y], \end{array} \,\,\, \text{for } |a| \leq L $$In our example, for which the base function is $\mathcal{P}_o = P(u)$, the parameters $a=u/2$, $A=\frac{1}{\sqrt{2}}$ (so that the maximum value is 1) and $L=1$. We can now write:
$$ (12) \hspace{40pt} \begin{array}{ll} A(u, y) & = &\int\limits_{-\infty}^{\infty}\mathcal{P}_o(t + a) \mathcal{P}_o^*(t - a)e^{j2\pi yt} dt \\ & = & \left(1 - \frac{|u|}{2}\right)\text{sinc} \left[y(2 - |u|)\right] \end{array} $$Now, from equation (8) we have $\mathcal{H}(u; W_{20}) = A(u, 2 u W_{20}/\lambda)$. Therefore we can write the expression for the OTF directly from the AF.
$$ (13) \hspace{40pt} \mathcal{H}(u; W_{20}) = \left(1 - \frac{|u|}{2} \right) \text{sinc} \left[ \frac{2 u W_{20}}{\lambda} (2 - |u|) \right] $$Note:
In [2], the authors have used the un-normalized $\text{sinc}$ function that is defined as $\text{sinc}(x)=\sin(x)/x$. Hence, there is an "extra" $\pi$ within the $\text{sinc}$ function in their paper. We use the normalized definition $\text{sinc}$ as it is used in [1], and in Numpy.
The figure below shows the function $\text{sinc}(x)$ for $x\in[-8,8]$. It can be seen in figure that roots of the (normalized) $\text{sinc}$ function occurs at $x=n, n\neq 0, n \in \mathbb{Z}$. Therefore the zero value loci for the AF associated with the one-dimensional rectangular pupil are found by equating $y(2 - |u|)=n$.
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 $W_{20}$. In each figure, the plot on the left shows the ambiguity function and several lines corresponding to different defocus error $W_{20}$. The zero value locai, $y = n/(2 - |u|)$ are denoted by the gray dashed lines. The intersection of zero value locai with the line(s) $y = 2uW_{20}/\lambda$, 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('2-D AF of 1-D 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 $W_{20} = n\frac{\lambda}{2}$ for $n=0, 1, 2, 3, \dots$, i.e. $n \in \mathbb{Z}$. The equation for the OTF was obtained directly from the ambiguity function by replacing $y$ with $2uW_{20}/\lambda$.
The points of intersection between lines $y = 2uW_{20}/\lambda$ and each zero loci curve $y = n/(2 - |u|)$ (for $n > 0$) can be found by equating $n/(2 - |u|) = 2uW_{20}/\lambda$, which gives $u = 1 \pm \sqrt{1 -\frac{n \lambda}{2 W_{20}}}$
For example, the line $y=2u \frac{W_{20}}{\lambda}\Big|_{W_{20}=1\lambda/2}$ intersects the zero loci curve $y = \frac{n}{(2 - |u|)}\Big|_{n=1}$ at $u=1$, the line $y=2u \frac{W_{20}}{\lambda}\Big|_{W_{20}=2\lambda/2}$ intersects the curve $y = \frac{n}{(2 - |u|)}\Big|_{n=1}$ at $u=1 \pm \sqrt{1 - 1/2}$, and the curve $y = \frac{n}{(2 - |u|)}\Big|_{n=2}$ at $u=1$, 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 $W_{20}$ is an integer multiple of $\lambda/2$. This corresponds to an odd number of zero-value OTF points in the OTF plots when the defocus error $W_{20}$ is an integer multiple of $\lambda/2$.
Note that we only found the abscissa ($u$ coordinates) for the intersection points above.
plot_1dRectAF(w20LambdaBy2=[0, 1, 2, 3])
Plots of AF and OTF for which $W_{20} = n\frac{\lambda}{2}$ for $n=0.5, 0.99, 1.4, 3.6, \dots$, i.e. $n \in \mathbb{R}, n \notin \mathbb{Z}$ (the particular values were arbitrarily chosen, and the $W_{20}=0$ line is plotted for comparison). In this case we find that there are even number of intersections for each line with the zero-loci curves in the AF plot. Correspondingly there are even number of zero-value OTF points for each OTF with focus error $W_{20}=n\lambda/2$ for $n \in \mathbb{R}, n \notin \mathbb{Z}, n \neq 0$.
plot_1dRectAF(w20LambdaBy2=[0, 0.5, 0.99, 1.5, 3.6])