Introduction to Fourier Optics

Fourier optics offers a robust analytical approach to understanding light propagation through optical systems by employing Fourier analysis techniques on optical fields. This framework elegantly connects image formation and optical resolution to the transmission of spatial information via light waves. Our exploration begins with examining complex transmittance functions, which give us fundamental insights into how various samples shape optical wavefronts. From this foundation, we will progress to the essential principles of Fourier optics and the associated diffraction integrals.

Transmission

When light interacts with an optical component or object, its amplitude and phase can be modified. Following Saleh and Teich’s formalism, we can characterize this interaction using the complex transmission factor \(t(x,y)\), which is defined as the ratio of the output field amplitude to the input field amplitude at each point \((x,y)\) in a plane:

\[t(x,y) = \frac{U_{\text{out}}(x,y)}{U_{\text{in}}(x,y)}\]

This transmission factor is generally complex-valued, with its magnitude representing amplitude modulation and its phase representing phase modulation of the incident light.

For a thin lens, the primary effect is phase modulation. To derive the transmission function for a thin lens, we need to consider the optical path length through the lens at each point. Consider a planoconvex lens with one flat surface and one spherical surface of radius \(R\). The thickness of the lens varies with position according to:

\[d(x,y) = d_0 - \frac{(x^2+y^2)}{2R}\]

where \(d_0\) is the thickness at the center. As light passes through the lens, it experiences a phase delay proportional to the optical path length, which is the product of the refractive index \(n\) and the physical path length:

\[\phi(x,y) = k \cdot n \cdot d(x,y) - k \cdot d(x,y)_{\text{air}}\]

where \(k = 2\pi/\lambda\) is the wavenumber. Simplifying:

\[\phi(x,y) = k(n-1)d(x,y) = k(n-1)\left(d_0 - \frac{(x^2+y^2)}{2R}\right)\]

The first term represents a constant phase shift that we can ignore, and the second term gives us the position-dependent phase modulation. For a lens with focal length \(f\), the relationship between \(R\) and \(f\) is given by the lensmaker’s formula, which for a planoconvex lens simplifies to:

\[(n-1)\frac{1}{R} = \frac{1}{f}\]

Substituting this into our phase equation:

\[\phi(x,y) = -k(n-1)\frac{(x^2+y^2)}{2R} = -k\frac{(x^2+y^2)}{2f}\]

The complex transmission factor is then:

\[t(x,y) = \exp[i\phi(x,y)] = \exp\left[-i\frac{k}{2f}(x^2+y^2)\right]\]

This quadratic phase factor represents the position-dependent phase delay introduced by the lens, with greater delays at the thicker portions of the lens.

Code 
# Parameters
wavelength = 0.5  # arbitrary units
k = 2*np.pi/wavelength
f = 10  # focal length
x = np.linspace(-5, 5, 400)
z = np.linspace(0, 15, 400)
X, Z = np.meshgrid(x, z)

# Calculate the phase profile at the lens (z=0)
lens_phase = k*(x**2)/(2*f)

# Calculate the phase in the x-z plane after the lens
phase_xz = np.zeros_like(X)
for i, z_val in enumerate(z):
    if z_val == 0:
        phase_xz[i,:] = lens_phase
    else:
        phase_xz[i,:] = k*((X[i,:]**2)/(2*f) - z_val)

# Wrap phase to [-π,π] for visualization
phase_xz_wrapped = np.angle(np.exp(1j*phase_xz))

fig, ax = plt.subplots(1, 2, figsize=get_size(12,6))

ax[0].plot(x, lens_phase)
ax[0].set_xlabel('x')
ax[0].set_ylabel('Phase [rad]')
ax[0].set_title('Lens Phase Profile (z=0)')
ax[0].grid(True)

im = ax[1].imshow(phase_xz_wrapped, extent=[x.min(), x.max(), z.min(), z.max()],
           origin='lower', cmap='twilight', aspect='auto')
ax[1].set_xlabel('x')
ax[1].set_ylabel('z')
ax[1].set_title('Wavefront Phase in x-z Plane')
ax[1].axhline(y=f, color='r', linestyle='--', alpha=0.7)
ax[1].text(x.max()-0.5, f+0.3, 'f', color='red')

fig.set_layout_engine('constrained')
plt.show()
Figure 1— Phase modulation effect of a thin lens on an incident plane wave. (a) The quadratic phase profile introduced by the lens at z=0. (b) The wavefront shape in the x-z plane after passing through the lens (twilight colormap: phase from \(-\pi\) to \(+\pi\) rad), showing how the initially flat wavefront is transformed into a converging spherical wavefront.

This transmission function is crucial in Fourier optics as it allows us to mathematically model how a lens transforms an incident field. When placed in the path of a light wave, the lens modifies the wavefront according to this transmission factor. As will be shown in the next lecture, a lens performs an exact spatial Fourier transform of the input field at its back focal plane.

Generalization to Arbitrary Thickness Objects

For arbitrary thickness objects, we can extend our treatment beyond the thin-element approximation. When light propagates through a medium of varying thickness and refractive index, the transmission function becomes:

\[t(x,y) = A(x,y) e^{i\phi(x,y)}\]

where \(A(x,y)\) represents amplitude modulation (absorption or gain) and \(\phi(x,y)\) represents phase modulation. For a thick object, the phase shift is given by the path integral through the object:

\[\phi(x,y) = k \int_\text{path} [n(x,y,z) - n_0] dz\]

where \(n(x,y,z)\) is the spatially varying refractive index within the object, \(n_0\) is the refractive index of the surrounding medium, and the integration is performed along the light path through the object.

This formulation accounts for complex three-dimensional objects where both the thickness and the refractive index may vary with position. For inhomogeneous media, we can express the transmission function as:

\[t(x,y) = \exp\left[ -\frac{1}{2}\alpha(x,y) + i k\int_0^{d(x,y)} [n(x,y,z)-n_0]\,dz \right]\]

where \(\alpha(x,y)\) is the absorption coefficient integrated along the path, and \(d(x,y)\) is the thickness at position \((x,y)\).

For many practical applications, this can be approximated by considering the effective phase and amplitude changes, leading to the more manageable form:

\[t(x,y) = \tau(x,y) e^{i k(n-n_0)d(x,y)}\]

where \(\tau(x,y)\) is the amplitude transmission coefficient accounting for reflection and absorption losses.

This mathematical framework will become crucially important later when we describe image formation from waves that have propagated through an object. The transmission function directly encodes how an object modifies both the amplitude and phase of the incident light field, which determines how the object appears in an imaging system. Different imaging modalities (such as bright-field, phase-contrast, or differential interference contrast microscopy) essentially measure different aspects of this complex transmission function, revealing different properties of the object being imaged.

Basic Definitions

The Fourier transform decomposes a function into its constituent frequency components. For a function \(f(x)\), its Fourier transform \(\hat{F}(\nu)\) is defined as:

\[\hat{F}(\nu) = \int_{-\infty}^{\infty} f(x)\, e^{-i2\pi\nu x}\, dx\]

The inverse Fourier transform reconstructs the original function:

\[f(x) = \int_{-\infty}^{\infty} \hat{F}(\nu)\, e^{i2\pi\nu x}\, d\nu\]

Here \(\nu\) denotes the spatial frequency in cycles per unit length. This “ordinary frequency” convention is standard in optics (Goodman, Saleh & Teich) and is used throughout this lecture. An alternative “angular” convention writes \(F(k) = \int f(x)e^{-ikx}dx\) with \(k = 2\pi\nu\) (units: rad/length); the two are related by \(\hat{F}(\nu) = F(2\pi\nu)\).

In optics, \(x\) typically represents a spatial coordinate and \(\nu\) the corresponding spatial frequency. When working with discrete data you will use the Discrete Fourier Transform (DFT), efficiently computed by the FFT algorithm:

Code 
import numpy as np
from scipy import fftpack

# Generate a simple signal
x = np.linspace(0, 1, 1000)  # spatial coordinate
f = np.sin(2*np.pi*5*x) + 0.5*np.sin(2*np.pi*10*x)  # 5 and 10 cycles/unit

# Compute the FFT
F = fftpack.fft(f)
freqs = fftpack.fftfreq(len(x), x[1]-x[0])  # cycles/length

# fftshift centres zero frequency
F_shifted = fftpack.fftshift(F)
freqs_shifted = fftpack.fftshift(freqs)

plt.figure(figsize=get_size(12, 5))
plt.subplot(121)
plt.plot(x, f)
plt.xlabel('position x')
plt.ylabel('amplitude')
plt.title('original signal')

plt.subplot(122)
plt.plot(freqs_shifted, np.abs(F_shifted))
plt.xlabel(r'spatial frequency $\nu$ (cycles/length)')
plt.ylabel(r'magnitude $|\hat{F}(\nu)|$')
plt.title('fourier spectrum')
plt.xlim(-20, 20)
plt.tight_layout()
plt.show()

Important Properties

  1. Linearity: \(\mathcal{F}\{af(x) + bg(x)\} = a\hat{F}(\nu) + b\hat{G}(\nu)\)

    The Fourier transform of a sum is the sum of the Fourier transforms, allowing us to analyze complex signals by breaking them into simpler components.

  2. Shift Theorem: \(\mathcal{F}\{f(x-a)\} = e^{-i2\pi\nu a}\hat{F}(\nu)\)

    A shift in the spatial domain corresponds to a linear phase change in the frequency domain.

  3. Convolution Theorem: \(\mathcal{F}\{f * g\} = \hat{F}(\nu) \cdot \hat{G}(\nu)\)

    Convolution in the spatial domain becomes pointwise multiplication in the frequency domain — and vice versa. This is the key theorem behind spatial filtering in the 4f system (covered in the next lecture).

  4. Parseval’s Theorem: \(\int |f(x)|^2 dx = \int |\hat{F}(\nu)|^2 d\nu\)

    The total energy (or power) is preserved under the Fourier transform.

Common Fourier Transform Pairs

Function Fourier Transform \(\hat{F}(\nu)\)
\(\delta(x)\) \(1\)
\(1\) \(\delta(\nu)\)
\(\text{rect}(x/a)\) \(a\,\text{sinc}(a\nu)\)
\(e^{-\pi x^2}\) \(e^{-\pi \nu^2}\)
\(\cos(2\pi a x)\) \(\tfrac{1}{2}[\delta(\nu-a) + \delta(\nu+a)]\)
\(e^{i2\pi a x}\) \(\delta(\nu - a)\)

Here \(\text{sinc}(u) = \sin(\pi u)/(\pi u)\). Understanding these transform pairs is essential for optical analysis: a rectangular aperture produces a sinc-function diffraction pattern, a Gaussian beam maintains its Gaussian profile under propagation, and a periodic grating (cosine) diffracts into exactly two discrete orders.

Code 
# Demonstrate the rect-sinc transform pair
def rect(x, width=1.0):
    return np.where(np.abs(x) <= width/2, 1.0, 0.0)

def sinc_norm(nu):
    result = np.ones_like(nu, dtype=float)
    nz = nu != 0
    result[nz] = np.sin(np.pi * nu[nz]) / (np.pi * nu[nz])
    return result

x_fine = np.linspace(-5, 5, 1000)
rect_func = rect(x_fine, width=2.0)

rect_ft = fftpack.fftshift(fftpack.fft(rect_func))
nu_vals = fftpack.fftshift(fftpack.fftfreq(len(x_fine), x_fine[1]-x_fine[0]))

theoretical_sinc = 2.0 * sinc_norm(2.0 * nu_vals)

plt.figure(figsize=get_size(12, 5))
plt.subplot(121)
plt.plot(x_fine, rect_func)
plt.xlabel('position x')
plt.title('rectangle function')

plt.subplot(122)
plt.plot(nu_vals, np.abs(rect_ft)/np.max(np.abs(rect_ft)), 'b-', label='FFT result')
plt.plot(nu_vals, np.abs(theoretical_sinc), 'r--', label=r'analytical $|\text{sinc}|$')
plt.xlabel(r'spatial frequency $\nu$')
plt.xlim(-2, 2)
#plt.legend()
plt.tight_layout()
plt.show()

Spatial Frequencies and Angular Spectrum

Building on our analysis of wave propagation through various optical elements, we now explore a fundamental concept in Fourier optics that connects spatial patterns to wave propagation directions. Just as a lens converts a plane wave into a converging spherical wave and a prism tilts the wavefront to change the propagation direction, complex spatial patterns decompose into multiple propagation directions — a relationship that will become essential when we discuss optical imaging systems, diffraction limits, and the resolution capabilities of microscopes and telescopes.

The Concept of Spatial Frequencies

Just as a temporal signal can be decomposed into frequency components through Fourier analysis, a spatial pattern or object can be represented as a superposition of spatial harmonic functions with different spatial frequencies. The spatial frequency \(\nu\) represents how rapidly the intensity or phase of an optical field changes with distance.

For a two-dimensional complex spatial harmonic function:

\[f(x,y) = A e^{i 2\pi(\nu_x x + \nu_y y)}\]

where \(\nu_x\) and \(\nu_y\) are the spatial frequencies in the \(x\) and \(y\) directions (in cycles per unit length) and \(A\) is the complex amplitude. Higher spatial frequencies correspond to finer details in an object, while lower spatial frequencies represent coarser features.

The function \(f(x,y)\) can be expressed as the Fourier transform of its spatial frequency components:

\[f(x,y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \hat{F}(\nu_x, \nu_y)\, e^{i 2\pi(\nu_x x + \nu_y y)}\, d\nu_x\, d\nu_y\]

where \(\hat{F}(\nu_x, \nu_y)\) is the spatial frequency spectrum of \(f(x,y)\).

A concrete illustration: the 2D Fourier transform of a real image reveals how its spatial frequency content is distributed. Low-\(\nu\) components cluster near the centre of the spectrum (slow variations, overall shape) while high-\(\nu\) components lie at the periphery (sharp edges, fine texture).

Code 
from PIL import Image
import matplotlib.pyplot as plt
import numpy as np

# Load image and convert to grayscale
img = plt.imread('img/frank.png')
if len(img.shape) > 2:
    img_gray = np.mean(img, axis=2)
else:
    img_gray = img

# Compute 2D FFT
img_fft = np.fft.fftshift(np.fft.fft2(img_gray))
freq_x = np.fft.fftshift(np.fft.fftfreq(img_gray.shape[1]))
freq_y = np.fft.fftshift(np.fft.fftfreq(img_gray.shape[0]))

fig, ax = plt.subplots(1, 3, figsize=get_size(12, 5))

ax[0].imshow(img_gray, cmap='gray')
ax[0].set_title('(a) Original Image')
ax[0].set_axis_off()

ax[1].imshow(np.log(np.abs(img_fft) + 1), cmap='viridis',
             extent=[freq_x.min(), freq_x.max(), freq_y.min(), freq_y.max()])
ax[1].set_title('(b) Magnitude')
ax[1].set_xlabel(r'$\nu_x$ (cycles/pixel)')
ax[1].set_ylabel(r'$\nu_y$ (cycles/pixel)')
ax[1].set_xticks([])
ax[1].set_yticks([])

ax[2].imshow(np.angle(img_fft), cmap='twilight',
             extent=[freq_x.min(), freq_x.max(), freq_y.min(), freq_y.max()])
ax[2].set_title('(c) Phase')
ax[2].set_xlabel(r'$\nu_x$ (cycles/pixel)')
ax[2].set_ylabel(r'$\nu_y$ (cycles/pixel)')
ax[2].set_xticks([])
ax[2].set_yticks([])

plt.tight_layout()
plt.show()
Figure 2— Spatial frequency analysis of an image. (a) Original grayscale image. (b) Magnitude of the 2D Fourier transform (log scale) — low frequencies cluster at the centre, high frequencies at the periphery. (c) Phase of the Fourier transform.

Correspondence to Plane Wave Angular Components

In the previous subsection we decomposed an arbitrary field into spatial harmonics \(e^{i2\pi(\nu_x x + \nu_y y)}\). Each such harmonic is, in fact, the cross-sectional trace of a plane wave. A monochromatic plane wave

\[U(x,y,z) = U_0\,e^{i(k_x x + k_y y + k_z z)}\]

evaluated at \(z=0\) gives exactly \(U_0\,e^{i2\pi(\nu_x x + \nu_y y)}\) provided we identify

\[k_x = 2\pi\nu_x, \qquad k_y = 2\pi\nu_y\]

The dispersion relation \(|\mathbf{k}|^2 = k_x^2+k_y^2+k_z^2 = (2\pi/\lambda)^2\) then fixes the axial component:

\[k_z = \frac{2\pi}{\lambda}\sqrt{1 - \lambda^2(\nu_x^2+\nu_y^2)}\]

and the corresponding propagation angles with respect to the \(z\)-axis are

\[\sin\theta_x = \lambda\nu_x, \qquad \sin\theta_y = \lambda\nu_y\]

So each spatial frequency pair \((\nu_x,\nu_y)\) in an object field corresponds uniquely to a plane wave traveling at angles \((\theta_x,\theta_y)\). Higher spatial frequencies (fine object features) diffract to larger angles; lower spatial frequencies (coarse features) travel nearly paraxially.

Figure 3— Principle of plane wave angular decomposition. (Image taken from Saleh/Teich “Principles of Photonics”)

Each spatial frequency component of an object therefore diffracts light into a specific direction. Behind the object, the plane-wave component \((\nu_x,\nu_y)\) propagates and acquires the phase factor \(e^{ik_z z}\):

\[U(x,y,z) = U_0\, e^{i(k_x x + k_y y)}\, e^{ik_z z}\]

where the propagation constant along \(z\) depends on the transverse spatial frequencies:

\[k_z = \sqrt{k^2 - k_x^2 - k_y^2} = \frac{2\pi}{\lambda}\sqrt{1 - \lambda^2(\nu_x^2 + \nu_y^2)}\]

This expression shows that spatial frequencies \(\nu_x^2+\nu_y^2 < 1/\lambda^2\) give real \(k_z\) (propagating waves), while \(\nu_x^2+\nu_y^2 > 1/\lambda^2\) gives imaginary \(k_z\) (evanescent waves that decay exponentially along \(z\)). The wavelength therefore sets an absolute upper bound on the spatial frequencies that can propagate freely — this is the origin of the diffraction limit.

Angular spectrum method. Because the plane waves \(e^{i2\pi(\nu_x x+\nu_y y)}\) form a complete basis, any field at the input plane can be decomposed as:

\[U(x,y,0) = \iint_{-\infty}^{\infty} \hat{U}(\nu_x,\nu_y)\, e^{i2\pi(\nu_x x+\nu_y y)}\,d\nu_x\,d\nu_y\]

where \(\hat{U}(\nu_x,\nu_y) = \mathcal{F}\{U(x,y,0)\}\) is the 2D Fourier transform. Since each plane-wave component propagates independently and acquires the phase \(e^{ik_z z}\), the field at any downstream plane \(z\) is:

\[U(x,y,z) = \iint_{-\infty}^{\infty} \hat{U}(\nu_x,\nu_y)\,H(\nu_x,\nu_y;z)\, e^{i2\pi(\nu_x x+\nu_y y)}\,d\nu_x\,d\nu_y\]

where

\[H(\nu_x,\nu_y;z) = e^{ik_z z} = \exp\!\left[i\frac{2\pi z}{\lambda}\sqrt{1 - \lambda^2(\nu_x^2+\nu_y^2)}\right]\]

is the transfer function of free space. Propagation is thus a linear filtering operation in the spatial-frequency domain: multiply the input spectrum by \(H\) and inverse-Fourier-transform. This is the exact scalar-wave result — no paraxial approximation has been made. Expanding \(k_z\) to lower order in \(\nu\) produces the paraxial approximations: the Fresnel transfer function is derived just below, and the Fraunhofer far-field limit follows in Section 0.0.8. In the next lecture the same transfer function \(H\) is the starting point for analysing how optical systems filter spatial frequencies.

Code 
wavelength_tf = 0.5
k_tf = 2 * np.pi / wavelength_tf
z_tf = 1.0
nu_max = 3.0 / wavelength_tf  # show beyond the evanescent boundary

nu = np.linspace(-nu_max, nu_max, 2000)

# Compute kz (complex-valued to capture evanescent regime)
arg = 1 - wavelength_tf**2 * nu**2
kz_tf = (2 * np.pi / wavelength_tf) * np.sqrt(arg.astype(complex))

# Transfer function along nu_x with nu_y = 0
H_tf = np.exp(1j * kz_tf * z_tf)

fig, axes = plt.subplots(1, 2, figsize=get_size(12, 5))

# Boundary between propagating and evanescent
nu_boundary = 1.0 / wavelength_tf

# (a) Magnitude — use log scale to reveal exponential decay
axes[0].plot(nu, np.abs(H_tf))
axes[0].axvline(x=nu_boundary, color='k', linestyle='--', linewidth=0.8, alpha=0.6, label=r'$|\nu|=1/\lambda$')
axes[0].axvline(x=-nu_boundary, color='k', linestyle='--', linewidth=0.8, alpha=0.6)
axes[0].set_xlabel(r'$\nu_x$ (cycles/length)')
axes[0].set_ylabel(r'$|H|$')
axes[0].set_yscale('log')
axes[0].set_ylim(bottom=1e-15, top=2)
axes[0].set_title(r'(a) Magnitude $|H(\nu_x, 0; z)|$')
axes[0].legend(fontsize=6, loc='upper right')

# (b) Phase — unwrap to remove 2π jumps
phase_tf = np.unwrap(np.angle(H_tf))
axes[1].plot(nu, phase_tf)
axes[1].axvline(x=nu_boundary, color='k', linestyle='--', linewidth=0.8, alpha=0.6, label=r'$|\nu|=1/\lambda$')
axes[1].axvline(x=-nu_boundary, color='k', linestyle='--', linewidth=0.8, alpha=0.6)
axes[1].set_xlabel(r'$\nu_x$ (cycles/length)')
axes[1].set_ylabel(r'$\arg(H)$ [rad]')
axes[1].set_title(r'(b) Phase $\arg\{H(\nu_x, 0; z)\}$')
axes[1].legend(fontsize=6, loc='upper right')

fig.set_layout_engine('constrained')
plt.show()
Figure 4— Magnitude and phase of the free-space transfer function \(H(\nu_x, 0; z)\) along the \(\nu_x\) axis for \(\lambda=0.5\) and \(z=1\). Inside \(|\nu_x| < 1/\lambda\) the function oscillates (propagating waves); outside it decays exponentially (evanescent waves).

Paraxial limit: the Fresnel transfer function. When the propagation angles are small, so that the transverse spatial frequencies satisfy \(\lambda^2(\nu_x^2+\nu_y^2)\ll 1\), the square root in the exponent of \(H\) can be expanded:

\[\sqrt{1-\lambda^2(\nu_x^2+\nu_y^2)} \approx 1 - \frac{1}{2}\lambda^2(\nu_x^2+\nu_y^2)\]

Inserting this into \(H\) and writing \(k = 2\pi/\lambda\) gives the Fresnel transfer function:

\[H_\mathrm{F}(\nu_x,\nu_y;z) = e^{ikz}\,\exp\!\left[-i\pi\lambda z(\nu_x^2+\nu_y^2)\right]\]

The exact transfer function is a pure phase that follows \(\sqrt{1-\lambda^2\nu^2}\); the Fresnel approximation replaces that phase with a parabola in \(\nu\). The constant term \(e^{ikz}\) is the on-axis phase, and the quadratic term is the leading angle-dependent correction. The approximation is accurate for paraxial spatial frequencies and fails only as \(\nu\) approaches the evanescent cutoff \(1/\lambda\), where the true phase curve bends away from the parabola.

The same propagation step has an equivalent form in real space. The inverse Fourier transform of \(H_\mathrm{F}\) is the Fresnel propagator \(h(x,y;z)\). Using the 2D Gaussian transform pair \(\mathcal{F}^{-1}\{e^{-i\pi\lambda z(\nu_x^2+\nu_y^2)}\} = (i\lambda z)^{-1}\exp[i\pi(x^2+y^2)/\lambda z]\) together with \(\pi/\lambda z = k/2z\):

\[h(x,y;z) = \mathcal{F}^{-1}\{H_\mathrm{F}\} = \frac{e^{ikz}}{i\lambda z}\,\exp\!\left[\frac{ik}{2z}(x^2+y^2)\right]\]

Propagation over a distance \(z\) can therefore be carried out two equivalent ways: multiply the angular spectrum by \(H_\mathrm{F}\), or convolve the field with \(h\). In Section 0.0.8 we reach this same propagator \(h(x,y;z)\) from the opposite starting point, the Huygens–Fresnel principle in real space, confirming that the angular-spectrum and Huygens–Fresnel descriptions are one theory.

We saw this principle in action when analyzing diffraction gratings, where we decomposed the grating’s periodic structure into angular components using the grating vector. For a grating with period \(d\), the spatial frequency is \(\nu_x = 1/d\), and the directions of diffracted orders are given by:

\[\sin\theta_m = m\lambda/d = m\lambda\nu_x\]

where \(m\) is the diffraction order. This shows how the grating’s spatial frequency determines the angles of diffracted light, which is a specific application of the more general Fourier relationship between spatial frequencies and propagation angles.

Code 
x = np.linspace(-5, 5, 400)
y = np.linspace(-5, 5, 400)
X, Y = np.meshgrid(x, y)

spatial_freq = 1.0  # cycles per unit length
wavelength = 0.5    # arbitrary units

object_pattern = np.cos(2 * np.pi * spatial_freq * X)

diffraction_angle = np.rad2deg(np.arcsin(wavelength * spatial_freq))

row_index = len(y) // 2
object_row = object_pattern[row_index, :]
ft_1d = np.fft.fftshift(np.fft.fft(object_row))
ft_magnitude_1d = np.abs(ft_1d)

fx = np.fft.fftshift(np.fft.fftfreq(len(x), x[1]-x[0]))
theta_x = np.rad2deg(np.arcsin(np.clip(wavelength * fx, -1, 1)))

fig, axs = plt.subplots(1, 3, figsize=get_size(12, 4))

im0 = axs[0].imshow(object_pattern, extent=[x.min(), x.max(), y.min(), y.max()],
                    cmap='gray', origin='lower')
axs[0].set_xlabel('x')
axs[0].set_ylabel('y')

axs[1].plot(fx, ft_magnitude_1d)
axs[1].set_xlabel(r'$\nu_x$ (cycles/length)')
axs[1].set_ylabel('magnitude')
axs[1].set_xlim(-3, 3)
axs[1].axvline(x=spatial_freq, color='r', linestyle='--', alpha=0.7)
axs[1].axvline(x=-spatial_freq, color='r', linestyle='--', alpha=0.7)
axs[1].text(spatial_freq, 0.7*np.max(ft_magnitude_1d), "+νₓ", color='red', ha='left')
axs[1].text(-spatial_freq, 0.7*np.max(ft_magnitude_1d), "-νₓ", color='red', ha='right')

valid_indices = np.abs(fx) <= 1/wavelength
axs[2].plot(theta_x[valid_indices], ft_magnitude_1d[valid_indices])
axs[2].set_xlabel(r'$\theta_x$ [°]')
axs[2].set_ylabel('magnitude')
axs[2].axvline(x=diffraction_angle, color='r', linestyle='--', alpha=0.7)
axs[2].axvline(x=-diffraction_angle, color='r', linestyle='--', alpha=0.7)
axs[2].text(diffraction_angle, 0.7*np.max(ft_magnitude_1d[valid_indices]), "+θ", color='red', ha='left')
axs[2].text(-diffraction_angle, 0.7*np.max(ft_magnitude_1d[valid_indices]), "-θ", color='red', ha='right')

plt.tight_layout()
plt.show()
Figure 5— Visualization of a 2D object containing a single spatial frequency in the x-direction. (a) The object pattern showing sinusoidal variation along x with frequency νₓ. (b) The Fourier transform magnitude, showing two symmetric peaks at ±νₓ. (c) The corresponding angular spectrum: the spatial frequency νₓ maps to diffraction angles ±θ according to sin(θ) = λνₓ.

Wave Propagation Through Objects

When a plane wave propagating along the \(z\)-axis encounters an object, its wavefronts are modified according to the object’s transmission function. The visualizations below illustrate how different optical elements transform the incident wavefront. The first three cases (free space, lens, prism) are accurate paraxial representations; the fourth panel shows propagation of an arbitrary phase object computed with the angular spectrum method (the exact scalar-theory treatment): the field at \(z=0\) is decomposed into plane waves via a Fourier transform, each component is propagated with its own phase factor \(e^{ik_z z}\), and the field is reconstructed by the inverse transform.

Code 
wavelength = 0.5  # arbitrary units
k = 2*np.pi/wavelength
x = np.linspace(-5, 5, 400)
z = np.linspace(0, 15, 400)
X, Z = np.meshgrid(x, z)

fig, ax = plt.subplots(2, 2, figsize=get_size(12,12))
ax = ax.flatten()

# (a) Plane wave in free space
phase_free_space = k*Z
ax[0].imshow(np.angle(np.exp(1j*phase_free_space)),
             extent=[x.min(), x.max(), z.min(), z.max()],
             origin='lower', cmap='twilight', aspect='auto')
ax[0].set_xlabel('x')
ax[0].set_ylabel('z')
ax[0].set_title('(a) Free Space Propagation')

# (b) Converging lens
f = 10
lens_phase = np.zeros_like(X)
for i, z_val in enumerate(z):
    if z_val < 0.1:
        lens_phase[i,:] = k*(x**2)/(2*f)
    else:
        lens_phase[i,:] = k*(X[i,:]**2/(2*f) - z_val)
ax[1].imshow(np.angle(np.exp(1j*lens_phase)),
             extent=[x.min(), x.max(), z.min(), z.max()],
             origin='lower', cmap='twilight', aspect='auto')
ax[1].set_xlabel('x')
ax[1].set_ylabel('z')
ax[1].set_title('(b) Converging Lens')
ax[1].axhline(y=f, color='r', linestyle='--', alpha=0.7)
ax[1].text(x.max()-0.5, f+0.3, 'f', color='red')

# (c) Prism
tilt_angle = 10
prism_phase = np.zeros_like(X)
for i, z_val in enumerate(z):
    if z_val < 0.1:
        prism_phase[i,:] = k*np.tan(np.deg2rad(tilt_angle))*x
    else:
        prism_phase[i,:] = k*(np.tan(np.deg2rad(tilt_angle))*x + z_val/np.cos(np.deg2rad(tilt_angle)))
ax[2].imshow(np.angle(np.exp(1j*prism_phase)),
             extent=[x.min(), x.max(), z.min(), z.max()],
             origin='lower', cmap='twilight', aspect='auto')
ax[2].set_xlabel('x')
ax[2].set_ylabel('z')
ax[2].set_title('(c) Prism')

# (d) Arbitrary phase object — angular spectrum propagation
def arbitrary_phase(xv):
    return 2*np.sin(xv) + 0.5*np.sin(3*xv)

U0 = np.exp(1j * k * arbitrary_phase(x))

# Precompute angular spectrum decomposition
dx = x[1] - x[0]
fx_as = np.fft.fftshift(np.fft.fftfreq(len(x), dx))
kx_as = 2*np.pi*fx_as
kz_as = np.where(k**2 - kx_as**2 >= 0, np.sqrt(np.maximum(k**2 - kx_as**2, 0)), 0.0)
U0_ft = np.fft.fftshift(np.fft.fft(U0))

arb_phase_map = np.zeros((len(z), len(x)))
for i, z_val in enumerate(z):
    propagator = np.exp(1j * kz_as * z_val)
    U_z = np.fft.ifft(np.fft.ifftshift(U0_ft * propagator))
    arb_phase_map[i, :] = np.angle(U_z)

ax[3].imshow(arb_phase_map,
             extent=[x.min(), x.max(), z.min(), z.max()],
             origin='lower', cmap='twilight', aspect='auto')
ax[3].set_xlabel('x')
ax[3].set_ylabel('z')
ax[3].set_title('(d) Arbitrary Phase Object\n(angular spectrum method)')

plt.tight_layout()
plt.show()
Figure 6— Wavefront propagation after transmission through different optical elements. (a) Free space: flat wavefronts maintained. (b) Converging lens: spherical wavefronts converging at the focal point. (c) Prism: tilted wavefronts, changed propagation direction. (d) Arbitrary phase object: phase map computed using the angular spectrum method, showing how rapidly-varying phase features (high spatial frequencies) diffract to large angles and disperse faster.

The top panels confirm that free space and paraxial elements (lens, prism) produce the expected wavefront shapes. In the bottom-right panel, the initial sinusoidal phase modulation \(\phi(x) = k(2\sin x + 0.5\sin 3x)\) decomposes into a discrete set of plane-wave components. Because each component propagates at a different angle (\(k_z = \sqrt{k^2-k_x^2}\)), the phase pattern evolves and eventually breaks up as higher spatial frequencies diffract more strongly. This is diffraction at work — the reason why the geometric-optics picture of rigid wavefront propagation fails at fine scales.

Diffraction Integrals

Having established that each spatial frequency component of an object propagates as a distinct plane wave, we now derive the scalar diffraction integral that connects the field \(U(x',y',0)\) at the object plane to the field \(U(x,y,z)\) at any downstream plane. You encountered the Huygens–Fresnel principle in your wave optics studies; here we make it quantitative and connect it directly to Fourier transforms.

The Huygens–Fresnel Integral

The Huygens–Fresnel principle states that every point on a wavefront acts as a secondary source of spherical waves. The field at any downstream point is the coherent superposition of all these secondary waves. In the scalar approximation and for monochromatic light of wavenumber \(k = 2\pi/\lambda\), this gives the diffraction integral:

\[U(x,y,z) = \frac{1}{i\lambda}\iint_{-\infty}^{\infty} U(x',y',0)\,\frac{e^{ikr}}{r}\,\cos\theta\,dx'\,dy'\]

where \(r = \sqrt{(x-x')^2+(y-y')^2+z^2}\) is the distance from source point \((x',y',0)\) to observation point \((x,y,z)\), and \(\cos\theta \approx z/r\) is the obliquity factor (close to 1 for paraxial propagation).

The Fresnel Approximation

For paraxial propagation — when transverse displacements are small compared to the propagation distance (\(|x-x'|, |y-y'| \ll z\)) — we expand \(r\) to second order:

\[r \approx z\left[1 + \frac{(x-x')^2+(y-y')^2}{2z^2}\right]\]

Keeping this quadratic correction in the phase \(e^{ikr}\) but approximating \(r\approx z\) in the amplitude \(1/r\) yields the Fresnel diffraction integral:

\[U(x,y,z) = \frac{e^{ikz}}{i\lambda z}\iint U(x',y',0)\,\exp\!\left[\frac{ik}{2z}\left((x-x')^2+(y-y')^2\right)\right]\!dx'\,dy'\]

This is a convolution of \(U(x',y',0)\) with the Fresnel propagator \(h(x,y;z) = \frac{e^{ikz}}{i\lambda z}\exp\!\left[\frac{ik}{2z}(x^2+y^2)\right]\). Note that Gaussian beam propagation from Lecture 5 is a direct consequence of the Fresnel integral applied to a Gaussian aperture field.

This real-space propagator is identical to the inverse Fourier transform of the Fresnel transfer function \(H_\mathrm{F}\) found earlier from the angular-spectrum expansion. The two derivations begin at opposite ends: one expands the exact transfer function \(H\) in the frequency domain, the other expands the spherical-wave kernel \(e^{ikr}/r\) in real space. They converge on the same \(h(x,y;z)\). This is the concrete confirmation that the angular-spectrum method and the Huygens–Fresnel principle are two representations of a single scalar diffraction theory, linked by a Fourier transform: multiplying the spectrum by \(H_\mathrm{F}\) and convolving the field with \(h\) are the same operation.

Fraunhofer Diffraction: The Far-Field Fourier Transform

In the far field — when \(z \gg \frac{k\,a^2}{2} = \frac{\pi a^2}{\lambda}\), where \(a\) is the transverse extent of the aperture — the quadratic phase \(e^{ik(x'^2+y'^2)/2z}\) accumulated across the object is negligible. Setting it to 1, the Fresnel integral factorises into:

\[U(x,y,z) = \frac{e^{ikz}}{i\lambda z}\,e^{\frac{ik}{2z}(x^2+y^2)}\iint U(x',y',0)\,e^{-i2\pi\!\left(\frac{x}{\lambda z}x'+\frac{y}{\lambda z}y'\right)}\!dx'\,dy'\]

Recognising the integral as a 2D Fourier transform evaluated at \(\nu_x = x/(\lambda z)\) and \(\nu_y = y/(\lambda z)\):

\[U(x,y,z) = \frac{e^{ikz}}{i\lambda z}\,e^{\frac{ik}{2z}(x^2+y^2)}\,\hat{U}\!\left(\nu_x = \frac{x}{\lambda z},\,\nu_y = \frac{y}{\lambda z}\right)\]

The intensity of the Fraunhofer diffraction pattern is therefore proportional to the power spectrum of the aperture field:

\[I(x,y,z) \propto \left|\hat{U}\!\left(\frac{x}{\lambda z},\frac{y}{\lambda z}\right)\right|^2\]

This is the central result of Fourier optics: the far-field diffraction pattern is the Fourier transform of the aperture. The spatial scale in the diffraction pattern is set by \(\lambda z\) — a wider aperture produces a narrower diffraction pattern, quantifying the wave-optics uncertainty principle.

Code 
from scipy.special import j1

N = 4096
x_ap = np.linspace(-8, 8, N)
dx_ap = x_ap[1] - x_ap[0]

# (a) Single slit, full width a
a_slit = 2.0
slit1 = np.where(np.abs(x_ap) <= a_slit/2, 1.0, 0.0)

# (b) Double slit: each slit has full width w, centres at ±d/2
w_slit = 0.8
d_slit = 2.0
slit2 = np.where(
    (np.abs(x_ap - d_slit/2) <= w_slit/2) | (np.abs(x_ap + d_slit/2) <= w_slit/2),
    1.0, 0.0)

# Fraunhofer intensity via FFT
def fraunhofer_1d(aperture, dx):
    ft = np.fft.fftshift(np.fft.fft(np.fft.ifftshift(aperture))) * dx
    nu = np.fft.fftshift(np.fft.fftfreq(len(aperture), dx))
    return nu, np.abs(ft)**2

nu1, I1 = fraunhofer_1d(slit1, dx_ap)
nu2, I2 = fraunhofer_1d(slit2, dx_ap)

# Analytical patterns
I1_anal = (a_slit * np.sinc(a_slit * nu1))**2
# double slit: FT = w*sinc(w*nu)*2*cos(pi*d*nu)
I2_anal = (2 * w_slit * np.sinc(w_slit * nu2) * np.cos(np.pi * d_slit * nu2))**2

# (c) 2D circular aperture, diameter D
M = 512
xy2 = np.linspace(-1, 1, M)
X2, Y2 = np.meshgrid(xy2, xy2)
R2 = np.sqrt(X2**2 + Y2**2)
D_circ = 0.3  # diameter
circ = np.where(R2 <= D_circ/2, 1.0, 0.0)

ft2 = np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(circ)))
I_circ = np.abs(ft2)**2

fig, axes = plt.subplots(3, 2, figsize=get_size(10, 14))

# Row 0: single slit
axes[0, 0].fill_between(x_ap, slit1, alpha=0.7)
axes[0, 0].set_xlim(-5, 5)
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('transmission')
axes[0, 0].set_title('(a) Single slit aperture')

axes[0, 1].plot(nu1, I1/I1.max(), label='FFT')
axes[0, 1].plot(nu1, I1_anal/I1_anal.max(), 'r--', label=r'sinc$^2$')
axes[0, 1].set_xlim(-1.5, 1.5)
axes[0, 1].set_xlabel(r'$\nu_x$ (cycles/length)')
axes[0, 1].set_ylabel('intensity (norm.)')
axes[0, 1].set_title('(b) Fraunhofer pattern')
#axes[0, 1].legend()

# Row 1: double slit
axes[1, 0].fill_between(x_ap, slit2, alpha=0.7)
axes[1, 0].set_xlim(-5, 5)
axes[1, 0].set_xlabel('x')
axes[1, 0].set_ylabel('transmission')
axes[1, 0].set_title('(c) Double slit aperture')

axes[1, 1].plot(nu2, I2/I2.max(), label='FFT')
axes[1, 1].plot(nu2, I2_anal/I2_anal.max(), 'r--', label='analytical')
axes[1, 1].set_xlim(-1.5, 1.5)
axes[1, 1].set_xlabel(r'$\nu_x$ (cycles/length)')
axes[1, 1].set_ylabel('intensity (norm.)')
axes[1, 1].set_title('(d) Fraunhofer pattern')
#axes[1, 1].legend()

# Row 2: circular aperture
axes[2, 0].imshow(circ, extent=[-1, 1, -1, 1], cmap='gray', origin='lower')
axes[2, 0].set_xlabel('x')
axes[2, 0].set_ylabel('y')
axes[2, 0].set_title('(e) Circular aperture')

vmax = np.percentile(I_circ, 99.5)
axes[2, 1].imshow((I_circ), cmap='inferno', origin='lower', vmax=vmax*3.5)
axes[2, 1].set_xlabel(r'$\nu_x$ (arb.)')
axes[2, 1].set_ylabel(r'$\nu_y$ (arb.)')
axes[2, 1].set_xlim(200,300)
axes[2, 1].set_ylim(200,300)
axes[2, 1].set_title('(f) Airy pattern')

plt.tight_layout()
plt.show()
Figure 7— Fraunhofer diffraction for three aperture geometries. Each row pairs the aperture (left) with its far-field intensity pattern (right). (a/b) Single slit of full width \(a\): sinc² pattern. (c/d) Double slit (width \(w\), separation \(d\)): sinc² envelope modulated by \(\cos^2(\pi d\nu_x)\). (e/f) Circular aperture of diameter \(D\): 2D Airy pattern.

The single-slit pattern confirms the sinc² form: the first zeros appear at \(\nu_x = \pm 1/a\), so a narrower slit produces a wider diffraction pattern. The double-slit pattern shows the familiar Young interference fringes (from the \(\cos^2\) factor) inside a sinc² envelope (from the finite slit width). The circular aperture yields the Airy pattern — the 2D circularly symmetric analogue of sinc² — whose first zero defines the Rayleigh criterion for resolution derived formally in the next lecture.