Seminar 3 — Fourier Optics Workshop

Introduction

In this seminar, we dive into Fourier optics, a framework connecting spatial frequencies to diffraction patterns. You will analyze how apertures shape light through the Fourier transform, design spatial filters, and understand diffraction-limited imaging. These concepts are central to microscopy, imaging systems, and holography.


Pen & Paper Problems

Problem 1: Rectangular Aperture Diffraction

A rectangular aperture of width \(a\) is illuminated by a plane wave of wavelength \(\lambda\). The diffraction pattern is observed on a screen at distance \(L\) from the aperture.

Tasks: 1. The Fourier transform of a rectangular aperture (rect function) is a sinc function: \(\mathcal{F}[\text{rect}(x/a)] = a \, \text{sinc}(u)\) where \(u = \pi a \sin\theta / \lambda\). 2. At what angles do the first zeros of the diffraction pattern occur? (Hint: \(\text{sinc}(x) = 0\) for \(x = \pm n\pi\), \(n \neq 0\).) 3. For \(a = 100\) μm and \(\lambda = 633\) nm (He-Ne laser), calculate the angular positions of the first three minima. 4. At distance \(L = 1\) m, what are the linear positions of these minima on the screen? 5. Sketch the diffraction pattern (intensity vs. angle), noting the central bright spot (Airy rectangle) and side lobes.


Problem 2: Circular Aperture & Airy Pattern

A circular aperture of diameter \(D\) produces the famous Airy pattern in the Fraunhofer (far-field) limit.

Tasks: 1. The intensity of the Airy pattern is: \[I(\theta) = I_0 \left[ \frac{2 J_1(x)}{x} \right]^2, \quad x = \frac{\pi D \sin\theta}{\lambda}\] where \(J_1\) is the first-order Bessel function.

  1. The first minimum occurs at \(x \approx 3.83\). Calculate the angular radius of the Airy disk (first dark ring) for:

    • Optical microscope objective with NA = 0.65 and \(\lambda = 500\) nm
    • Expression: radius \(\approx 1.22 \lambda / D\) or \(\approx 0.61 \lambda / \text{NA}\) for a lens.
  2. For a microscope with NA = 0.65 and λ = 500 nm, what is the Airy disk radius in μm?

  3. Compare with the Rayleigh criterion: two point sources are resolved if separated by more than one Airy disk radius. What is the minimum resolvable distance?


Problem 3: Double-Slit Interference & Diffraction

Two slits, each of width \(a\), separated by distance \(d\), are illuminated by a plane wave. The pattern combines single-slit diffraction (envelope) and two-slit interference (fringes).

Tasks: 1. The combined intensity is: \[I(\theta) = I_0 \left[ \frac{\sin(\beta)}{\beta} \right]^2 \cos^2(\delta)\] where \(\beta = \pi a \sin\theta / \lambda\) (single slit) and \(\delta = \pi d \sin\theta / \lambda\) (slit separation).

  1. For \(a = 100\) μm, \(d = 0.5\) mm, \(\lambda = 633\) nm, identify:

    • Spacing between interference fringes (given by the cosine term)
    • Period of the diffraction envelope (given by the sinc term)
  2. Sketch the combined pattern, showing fringes and envelope.

  3. Spatial frequency: The spacing between fringes corresponds to what spatial frequency (cycles per mm)?


Problem 4: Spatial Filtering in 4f System (Challenge)

A 4f system (two lenses of focal length \(f\), separated by distance \(2f\)) acts as a Fourier transform processor. A spatial filter can be placed at the Fourier plane to selectively remove or enhance features.

Tasks: 1. Explain how a 4f system creates a Fourier plane at the midpoint between the two lenses. 2. A low-pass filter (opaque except near the center) removes fine details. Describe what you would see on the output screen if a small aperture is placed at the Fourier plane. 3. A high-pass filter blocks low spatial frequencies (large features). What happens to a grating pattern or text? 4. Design a notch filter to remove a periodic grating (e.g., unwanted stripes in an image). Where would you place the filter?


Problem 5: SLM as Programmable Spatial Filter

A spatial light modulator (SLM) can imprint an arbitrary phase or amplitude pattern on a wavefront.

Tasks: 1. Explain how an SLM can replace fixed spatial filters in a 4f system. 2. A programmable phase mask \(\phi(x, y)\) on an SLM modulates the incident field to \(E'(x,y) = E(x,y) e^{i\phi(x,y)}\). 3. If the SLM imprints a blazed grating phase pattern, what is the effect on the output? (Hint: a grating diffracts light into different orders.) 4. How could an SLM implement a generalized wavefront shaping to focus light onto arbitrary 3D patterns or generate spatial modes?


Python Exercises

Setup


Exercise 1: Aperture Diffraction Patterns

Simulate diffraction from rectangular and circular apertures using the Fourier transform.

Code
# Create a simple aperture in the spatial domain
size = 1024
aperture_size_pixels = 100

# Rectangular aperture
rect_aperture = np.zeros((size, size))
rect_aperture[size//2 - aperture_size_pixels//2 : size//2 + aperture_size_pixels//2,
              size//2 - aperture_size_pixels//4 : size//2 + aperture_size_pixels//4] = 1

# Circular aperture
y, x = np.ogrid[:size, :size]
cy, cx = size // 2, size // 2
radius = aperture_size_pixels // 2
circ_aperture = ((x - cx)**2 + (y - cy)**2 <= radius**2).astype(float)

# Compute Fourier transforms (diffraction patterns)
diffr_rect = np.abs(fftshift(fft2(rect_aperture)))**2
diffr_circ = np.abs(fftshift(fft2(circ_aperture)))**2

# Normalize for display (log scale)
diffr_rect_log = np.log10(diffr_rect + 1)
diffr_circ_log = np.log10(diffr_circ + 1)

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

# Row 1: Apertures
axes[0, 0].imshow(rect_aperture, cmap='gray')
axes[0, 0].set_title('Rectangular Aperture')
axes[0, 0].set_aspect('equal')

axes[0, 1].imshow(circ_aperture, cmap='gray')
axes[0, 1].set_title('Circular Aperture')
axes[0, 1].set_aspect('equal')

# Row 2: Diffraction patterns (log scale)
im_rect = axes[1, 0].imshow(diffr_rect_log, cmap='hot', extent=[-10, 10, -10, 10])
axes[1, 0].set_title('Diffraction: Rectangle')
axes[1, 0].set_xlabel('Spatial frequency')
axes[1, 0].set_ylabel('Spatial frequency')

im_circ = axes[1, 1].imshow(diffr_circ_log, cmap='hot', extent=[-10, 10, -10, 10])
axes[1, 1].set_title('Diffraction: Circle (Airy)')
axes[1, 1].set_xlabel('Spatial frequency')
axes[1, 1].set_ylabel('Spatial frequency')

plt.tight_layout()
plt.savefig('img/aperture_diffraction.png', dpi=150, bbox_inches='tight')
plt.close()

print("Diffraction patterns generated.")
print(f"Central intensity (rect): {diffr_rect[size//2, size//2]:.2e}")
print(f"Central intensity (circ): {diffr_circ[size//2, size//2]:.2e}")
Diffraction patterns generated.
Central intensity (rect): 2.50e+07
Central intensity (circ): 6.15e+07

Exercise 2: Airy Disk & Resolution

Calculate and plot the Airy pattern for optical imaging.

Code
# Airy pattern: intensity of circular aperture in far field
def airy_intensity(x, x_scale):
    """
    Airy pattern intensity.
    x: normalized coordinate
    x_scale: normalization factor
    """
    x_norm = x * x_scale
    # Avoid division by zero at x=0
    J1 = j1(x_norm)
    I = np.where(x_norm == 0, 1.0, (2*J1/x_norm)**2)
    return I

# Parameters
NA = 0.65  # numerical aperture
wavelength = 500e-9  # 500 nm (green)
lambda_um = wavelength * 1e6

# Airy disk radius
airy_radius = 1.22 * wavelength / (2 * NA)  # For a lens

# Rayleigh criterion: two points separated by one Airy radius are just resolved
rayleigh_limit = airy_radius

# Spatial coordinate
r = np.linspace(0, 10 * airy_radius, 1000)
x_scale = np.pi / airy_radius

# Calculate Airy intensity
I_airy = airy_intensity(r, x_scale)

# Plot line profile
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=get_size(14, 5))

ax1.plot(r*1e6, I_airy, 'b-', linewidth=2)
ax1.axvline(airy_radius*1e6, color='red', linestyle='--', alpha=0.7, label=f'Airy radius = {airy_radius*1e6:.2f} µm')
ax1.fill_between(r*1e6, 0, I_airy, alpha=0.2)
ax1.set_xlabel('Radial distance (µm)')
ax1.set_ylabel('Intensity (normalized)')
ax1.grid(True, alpha=0.3)
ax1.legend(fontsize=8)
ax1.set_xlim([0, 10*airy_radius*1e6])

# 2D Airy pattern
x2d = np.linspace(-15*airy_radius, 15*airy_radius, 500)
y2d = np.linspace(-15*airy_radius, 15*airy_radius, 500)
X2d, Y2d = np.meshgrid(x2d, y2d)
R2d = np.sqrt(X2d**2 + Y2d**2)

I_2d = airy_intensity(R2d, x_scale)

im = ax2.imshow(I_2d, extent=[-15, 15, -15, 15], cmap='hot', origin='lower')
ax2.set_xlabel('x (µm)')
ax2.set_ylabel('y (µm)')
ax2.set_title('2D Airy Pattern')
ax2.set_aspect('equal')
plt.colorbar(im, ax=ax2, label='Intensity')

plt.tight_layout()
plt.savefig('img/airy_pattern.png', dpi=150, bbox_inches='tight')
plt.close()

print(f"NA = {NA}")
print(f"Wavelength = {lambda_um:.3f} µm")
print(f"Airy disk radius = {airy_radius*1e6:.3f} µm")
print(f"Rayleigh resolution limit = {rayleigh_limit*1e6:.3f} µm")
print(f"First minimum at x ≈ 3.83 (normalized coordinate)")
/var/folders/t4/_9qps8wj56jc60nwkr3nrcr00000gn/T/ipykernel_10081/2119980280.py:11: RuntimeWarning:

invalid value encountered in divide
NA = 0.65
Wavelength = 0.500 µm
Airy disk radius = 0.469 µm
Rayleigh resolution limit = 0.469 µm
First minimum at x ≈ 3.83 (normalized coordinate)

Exercise 3: Double-Slit Pattern

Simulate the combined effect of single-slit diffraction and two-slit interference.

Code
def double_slit_pattern(theta, slit_width, slit_separation, wavelength):
    """
    Combined intensity from two slits.
    """
    # Single-slit diffraction (envelope)
    beta = np.pi * slit_width * np.sin(theta) / wavelength
    sinc_part = np.sinc(beta / np.pi)**2

    # Two-slit interference (fringes)
    delta = np.pi * slit_separation * np.sin(theta) / wavelength
    interference_part = np.cos(delta)**2

    # Combined pattern
    I = sinc_part * interference_part

    return I, sinc_part, interference_part

# Parameters
slit_width = 100e-6  # 100 µm
slit_separation = 0.5e-3  # 0.5 mm
wavelength = 633e-9  # He-Ne laser

# Angle range (in radians and degrees)
theta_max = 0.05  # 50 mrad, about 3 degrees
theta = np.linspace(-theta_max, theta_max, 1000)

I_combined, I_envelope, I_interference = double_slit_pattern(theta, slit_width, slit_separation, wavelength)

# Convert to degrees for plotting
theta_deg = theta * 180 / np.pi

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=get_size(12, 8))

# Combined pattern
ax1.plot(theta_deg, I_combined, 'b-', linewidth=2, label='Combined (fringes × envelope)')
ax1.plot(theta_deg, I_envelope, 'r--', linewidth=1.5, label='Diffraction envelope')
ax1.fill_between(theta_deg, 0, I_combined, alpha=0.2, color='blue')
ax1.set_ylabel('Intensity (normalized)')
ax1.grid(True, alpha=0.3)
ax1.legend(fontsize=8)
ax1.set_xlim([-3, 3])

# Fringes zoom
ax2.plot(theta_deg, I_combined, 'b-', linewidth=2)
ax2.set_xlabel('Angle (°)')
ax2.set_ylabel('Intensity (normalized)')
ax2.grid(True, alpha=0.3)
ax2.set_xlim([-1, 1])

plt.tight_layout()
plt.savefig('img/double_slit_pattern.png', dpi=150, bbox_inches='tight')
plt.close()

# Calculate fringe spacing
fringe_spacing_rad = wavelength / slit_separation
fringe_spacing_deg = fringe_spacing_rad * 180 / np.pi

print(f"Slit width = {slit_width*1e6:.1f} µm")
print(f"Slit separation = {slit_separation*1e3:.2f} mm")
print(f"Wavelength = {wavelength*1e9:.0f} nm")
print(f"Fringe spacing = {fringe_spacing_rad*1e3:.3f} mrad = {fringe_spacing_deg:.4f}°")
print(f"Diffraction envelope (first minimum) at ≈ {wavelength/slit_width*180/np.pi:.3f}°")
Slit width = 100.0 µm
Slit separation = 0.50 mm
Wavelength = 633 nm
Fringe spacing = 1.266 mrad = 0.0725°
Diffraction envelope (first minimum) at ≈ 0.363°

Exercise 4: Spatial Filtering in Fourier Domain

Demonstrate low-pass, high-pass, and notch filtering on a test image.

Code
# Create a test image with multiple features
size = 256
test_image = np.zeros((size, size))

# Add a grating pattern (periodic structure)
x_freq = np.linspace(0, 2*np.pi*10, size)
grating = 0.5 * (1 + np.cos(x_freq))
test_image += grating[np.newaxis, :] * 0.5

# Add a Gaussian blob (smooth feature)
y, x = np.ogrid[:size, :size]
cy, cx = size // 2, size // 2
blob = np.exp(-((x - cx)**2 + (y - cy)**2) / (2*50**2))
test_image += blob * 0.5

# Add some noise
np.random.seed(42)
test_image += 0.1 * np.random.randn(size, size)

# Normalize
test_image = (test_image - test_image.min()) / (test_image.max() - test_image.min())

# Compute Fourier transform
F = fftshift(fft2(test_image))
F_mag = np.abs(F)
F_mag_log = np.log10(F_mag + 1)

# Create filters
cy, cx = size // 2, size // 2
yy, xx = np.ogrid[:size, :size]
r = np.sqrt((xx - cx)**2 + (yy - cy)**2)

# Low-pass filter (keep center)
low_pass = r < 20
F_lowpass = F * low_pass
img_lowpass = np.abs(ifft2(ifftshift(F_lowpass)))

# High-pass filter (remove center)
high_pass = r > 10
F_highpass = F * high_pass
img_highpass = np.abs(ifft2(ifftshift(F_highpass)))

# Notch filter (remove specific region)
notch = np.ones((size, size), dtype=bool)
notch_x_range = (xx > cx - 5) & (xx < cx + 5) & (yy > cy - 30) & (yy < cy + 30)
notch[notch_x_range] = False
F_notch = F * notch
img_notch = np.abs(ifft2(ifftshift(F_notch)))

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

# Original image and its FFT
axes[0, 0].imshow(test_image, cmap='gray')
axes[0, 0].set_title('Original Image')
axes[0, 0].set_aspect('equal')

axes[0, 1].imshow(F_mag_log, cmap='hot')
axes[0, 1].set_title('Fourier Spectrum (log)')
axes[0, 1].set_aspect('equal')

axes[0, 2].imshow(low_pass, cmap='gray')
axes[0, 2].set_title('Low-pass Filter Mask')
axes[0, 2].set_aspect('equal')

# Filtered results
axes[1, 0].imshow(img_lowpass, cmap='gray')
axes[1, 0].set_title('Low-pass Result')
axes[1, 0].set_aspect('equal')

axes[1, 1].imshow(img_highpass, cmap='gray')
axes[1, 1].set_title('High-pass Result')
axes[1, 1].set_aspect('equal')

axes[1, 2].imshow(img_notch, cmap='gray')
axes[1, 2].set_title('Notch Filter Result')
axes[1, 2].set_aspect('equal')

plt.tight_layout()
plt.savefig('img/spatial_filtering.png', dpi=150, bbox_inches='tight')
plt.close()

print("Spatial filtering demonstration complete.")
print(f"Original image size: {test_image.shape}")
print(f"Low-pass: removes high-frequency details (sharp edges, noise)")
print(f"High-pass: removes low-frequency components (background, large structures)")
print(f"Notch: removes specific spatial frequencies (e.g., periodic artifacts)")
Spatial filtering demonstration complete.
Original image size: (256, 256)
Low-pass: removes high-frequency details (sharp edges, noise)
High-pass: removes low-frequency components (background, large structures)
Notch: removes specific spatial frequencies (e.g., periodic artifacts)

Exercise 5: 4f Spatial Filter System

Simulate a 4f imaging system with an adjustable spatial filter.

Code
# Create input image with text-like pattern
size = 256
input_image = np.zeros((size, size))

# Add some lines (simulate text or features)
input_image[50:100, 40:50] = 1
input_image[80:90, 60:150] = 1
input_image[120:130, 40:200] = 1

# Add a smooth background
y, x = np.ogrid[:size, :size]
cy, cx = size // 2, size // 2
background = 0.3 * np.exp(-((x - cx)**2 + (y - cy)**2) / (2*100**2))
input_image += background

# Normalize
input_image = (input_image - input_image.min()) / (input_image.max() - input_image.min() + 1e-6)

# Simulate 4f system: FFT -> filter -> inverse FFT
F_input = fftshift(fft2(input_image))
F_mag = np.abs(F_input)
F_mag_log = np.log10(F_mag + 1)

# Create aperture in Fourier plane (spatial filter)
cy, cx = size // 2, size // 2
yy, xx = np.ogrid[:size, :size]
r = np.sqrt((xx - cx)**2 + (yy - cy)**2)

# Vary filter radius
filter_radii = [10, 30, 50]
outputs = []

fig, axes = plt.subplots(2, 4, figsize=get_size(15, 8))

# Show input and spectrum
axes[0, 0].imshow(input_image, cmap='gray')
axes[0, 0].set_title('Input Image')
axes[0, 0].set_aspect('equal')

axes[1, 0].imshow(F_mag_log, cmap='hot', vmin=0, vmax=10)
axes[1, 0].set_title('Fourier Spectrum')
axes[1, 0].set_aspect('equal')

# Apply different filters
for idx, radius in enumerate(filter_radii):
    # Circular aperture at Fourier plane
    aperture = r < radius

    # Filter the spectrum
    F_filtered = F_input * aperture

    # Inverse FFT
    output = np.abs(ifft2(ifftshift(F_filtered)))

    # Show filter mask
    axes[0, idx+1].imshow(aperture, cmap='gray')
    axes[0, idx+1].set_title(f'Filter (radius = {radius})')
    axes[0, idx+1].set_aspect('equal')

    # Show output
    axes[1, idx+1].imshow(output, cmap='gray')
    axes[1, idx+1].set_title(f'Output (r = {radius})')
    axes[1, idx+1].set_aspect('equal')

plt.tight_layout()
plt.savefig('img/fourier_filter_system.png', dpi=150, bbox_inches='tight')
plt.close()

print("4f system simulation complete.")
print("Smaller filter radius → more low-pass filtering (blur)")
print("Larger filter radius → more high-frequency content preserved (sharp)")
4f system simulation complete.
Smaller filter radius → more low-pass filtering (blur)
Larger filter radius → more high-frequency content preserved (sharp)

Solutions Guide

Problem 1: The first zeros of the sinc function occur at \(u = \pm \pi, \pm 2\pi, ...\) For a 100 μm aperture and 633 nm, the first minimum is at approximately 3.6°.

Problem 2: The Airy disk radius is \(r \approx 1.22 \lambda / D\) or equivalently \(r \approx 0.61 \lambda / \text{NA}\). For NA = 0.65 and λ = 500 nm, this gives approximately 0.47 μm.

Problem 3: The fringe spacing is determined by slit separation, while the envelope width is determined by slit width. With values given, fringes are very closely spaced compared to the slow diffraction envelope.

Problem 4: A 4f system maps spatial frequencies to spatial positions in the Fourier plane. Low-pass filtering (small aperture) removes high spatial frequencies and blurs the image. High-pass filtering (large opaque region near axis) removes low frequencies and enhances edges.

Problem 5: SLMs enable programmable spatial filtering, allowing dynamic wavefront shaping, beam steering, and even calculation of arbitrary Fourier transforms in real time.


Summary

Fourier optics unifies diffraction theory with spatial frequency filtering, providing powerful tools for image processing and optical design. Understanding these concepts is essential for advanced topics in microscopy, holography, and optical communication.