EE5621 Homework Explanations

This document provides detailed explanations for homework problems from EE5621: Physical Optics, based on Goodman’s Introduction to Fourier Optics.

---

Homework 1: Fourier Analysis Foundations

What is the Fourier Transform?

The Fourier transform is a map that takes a function in one domain and produces a function in another domain:

$$ \mathcal{F}: f(x) \mapsto F(\nu) $$

Concretely, the map is defined by this integral:

$$ F(\nu) = \int_{-\infty}^{\infty} f(x) \, e^{-j2\pi \nu x} \, dx $$

	Input	Output
Forward transform $\mathcal{F}$	$f(x)$ — function of position/time	$F(\nu)$ — function of frequency
Inverse transform $\mathcal{F}^{-1}$	$F(\nu)$ — function of frequency	$f(x)$ — function of position/time

The two domains describe the same underlying object — just from different perspectives:

Domain	What it shows	Axes
Spatial/Time	Value at each location/moment	Position $x$ or time $t$
Frequency	Amount of each sinusoidal component	Frequency $\nu$ or $f$

Why do we care?

1. Simplifies convolution: Convolution in the spatial domain becomes multiplication in the frequency domain. Filtering an image is just multiplying spectra.

2. Reveals hidden structure: Periodic patterns that are hard to see spatially become obvious peaks in the frequency domain.

3. Connects to physics: Diffraction, interference, and imaging are naturally described in terms of spatial frequencies. A lens performs an optical Fourier transform.

4. Enables efficient computation: The Fast Fourier Transform (FFT) makes frequency analysis computationally practical.

---

Fourier Transforms in Optical Imaging

In imaging, the Fourier transform maps from the spatial domain (brightness at each location) to the spatial frequency domain (how much of each grating pattern is present):

$$ \mathcal{F}: f(x, y) \mapsto F(\nu_x, \nu_y) $$

What Does the Frequency Domain Mean for an Image?

Each point $(\nu_x, \nu_y)$ in the frequency domain represents a specific 2D sinusoidal grating — stripes at a particular spacing and orientation. The value $F(\nu_x, \nu_y)$ tells you the amplitude and phase of that pattern in the image.

Feature in image	Location in frequency domain
Uniform regions, slow gradients	Near origin (low $\nu_x$, $\nu_y$)
Edges, fine textures, sharp details	Far from origin (high $\nu_x$, $\nu_y$)
Vertical stripes/edges	Along the $\nu_x$ axis
Horizontal stripes/edges	Along the $\nu_y$ axis
Periodic textures (fabric, bricks)	Discrete peaks at the texture’s frequency

Why Lenses Perform Fourier Transforms

This isn’t just mathematical convenience — light physically computes Fourier transforms through diffraction.

When coherent light passes through an object and then through a lens:

1. The object modulates the light’s amplitude/phase at each point $(x, y)$
2. Diffraction causes each point to spread into waves traveling in all directions
3. The lens focuses waves traveling in the same direction to the same point
4. At the back focal plane, the light distribution is the Fourier transform of the object

This is why Fourier optics is so powerful: the mathematics directly describes what light does. Diffraction is Fourier transformation.

Filtering in the Frequency Domain

Since the frequency domain separates image features by their spatial scale, we can selectively modify them:

Operation	Frequency domain action	Effect
Blur	Multiply by low-pass filter (attenuate high $\nu$)	Removes fine detail, smooths edges
Sharpen	Boost high frequencies	Enhances edges, amplifies noise
Remove periodic noise	Zero out specific frequency peaks	Eliminates repeating patterns
Edge detection	High-pass filter (attenuate low $\nu$)	Keeps only rapid variations

The Complete Imaging Pipeline

A coherent optical imaging system performs:

$$ \text{Object } f(x,y) \xrightarrow{\text{Lens 1}} \text{Spectrum } F(\nu_x, \nu_y) \xrightarrow{\times H(\nu_x,\nu_y)} \text{Filtered spectrum} \xrightarrow{\text{Lens 2}} \text{Image } g(x,y) $$

where $H(\nu_x, \nu_y)$ is the transfer function of any filtering or apertures in the system.

This is why understanding Fourier transforms is essential for optics: every lens, aperture, and optical element can be understood as an operation in the frequency domain.

---

Problem 1: Properties of Delta Functions (Goodman 2.1)

The Dirac delta function $\delta(x)$ is not a function in the traditional sense — it’s a distribution or generalized function. It’s defined by its behavior under integration:

$$ \int_{-\infty}^{\infty} f(x) \delta(x - x_0) \, dx = f(x_0) $$

This is called the sifting property: the delta “picks out” the value of $f$ at the location of the spike.

---

Part (a): Scaling the 2D Delta Function

$$ \delta(ax, by) = \frac{1}{|ab|} \delta(x, y) $$

What this means: When you “squeeze” the coordinates by factors $a$ and $b$, the delta function must compensate to maintain unit volume under the integral.

The derivation uses two facts:

1. 1D scaling property: $\delta(cx) = \frac{1}{|c|} \delta(x)$

2. Separability: The 2D delta is the product of two 1D deltas:

$$ \delta(x, y) = \delta(x) \cdot \delta(y) $$

Combining these:

$$ \delta(ax, by) = \delta(ax) \cdot \delta(by) = \frac{1}{|a|}\delta(x) \cdot \frac{1}{|b|}\delta(y) = \frac{1}{|ab|}\delta(x,y) $$

---

Part (b): The Dirac Comb (Shah Function)

$$ \text{comb}(ax) \cdot \text{comb}(by) = \frac{1}{|ab|} \sum_{n=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} \delta\left(x - \frac{n}{a}, y - \frac{m}{b}\right) $$

What is the Dirac comb?

The Dirac comb is an infinite train of equally-spaced delta functions:

$$ \text{comb}(x) = \sum_{n=-\infty}^{\infty} \delta(x - n) $$

It’s also called the Shah function (denoted Ш) because it looks like a row of spikes at every integer.

The scaling effect:

When you write $\text{comb}(ax)$, you’re changing the spacing between spikes:

• If $a > 1$: spikes are closer together (spacing = $1/a$)
• If $a < 1$: spikes are farther apart (spacing = $1/a$)

The $1/|ab|$ factor in front compensates for the changed density of spikes to preserve normalization.

---

Problem 2: Fourier Transform Relations (Goodman 2.2)

These problems establish fundamental transform pairs that appear repeatedly in optics and signal processing.

---

Part (a): Rectangle to Sinc

$$ \mathcal{F}\{\text{rect}(x) \cdot \text{rect}(y)\} = \text{sinc}(f_X) \cdot \text{sinc}(f_Y) $$

The rectangle function:

$$ \text{rect}(x) = \begin{cases} 1 & |x| < 1/2 \\ 0 & |x| > 1/2 \end{cases} $$

It’s a “top hat” or “boxcar” function — constant over a unit interval, zero elsewhere.

The sinc function:

$$ \text{sinc}(x) = \frac{\sin(\pi x)}{\pi x} $$

This oscillating function with decaying envelope is ubiquitous in signal processing.

Why does rect transform to sinc?

Computing the 1D Fourier transform directly:

$$ \mathcal{F}\{\text{rect}(x)\} = \int_{-1/2}^{1/2} e^{-j2\pi f_X x} \, dx = \frac{e^{-j\pi f_X} - e^{j\pi f_X}}{-j2\pi f_X} = \frac{\sin(\pi f_X)}{\pi f_X} = \text{sinc}(f_X) $$

The key concept: Separability

When a 2D function can be written as $g(x,y) = g_1(x) \cdot g_2(y)$, its 2D Fourier transform factors:

$$ \mathcal{F}\{g_1(x) \cdot g_2(y)\} = G_1(f_X) \cdot G_2(f_Y) $$

This dramatically simplifies 2D transform calculations — you can do two 1D transforms instead of one 2D integral.

---

Part (b): Triangle to Sinc-Squared

$$ \mathcal{F}\{\Lambda(x) \cdot \Lambda(y)\} = \text{sinc}^2(f_X) \cdot \text{sinc}^2(f_Y) $$

The triangle function:

$$ \Lambda(x) = \begin{cases} 1 - |x| & |x| < 1 \\ 0 & |x| \geq 1 \end{cases} $$

It’s a tent-shaped function that rises linearly from 0 to 1 and back down.

Using the Convolution Theorem:

One of the most powerful results in Fourier analysis:

$$ \mathcal{F}\{f * g\} = F \cdot G $$

Convolution in the spatial domain becomes multiplication in the frequency domain.

Therefore:

$$ \mathcal{F}\{\Lambda(x)\} = \mathcal{F}\{\text{rect} * \text{rect}\} = \text{sinc}(f_X) \cdot \text{sinc}(f_X) = \text{sinc}^2(f_X) $$

Physical interpretation:

The sinc² intensity pattern (called the Airy pattern for circular apertures) describes the diffraction pattern from many optical elements.

---

Part (c): Constant to Delta

$$ \mathcal{F}\{1\} = \delta(f_X, f_Y) $$

What this means:

A constant signal (same value everywhere in space) transforms to a single spike at the origin in frequency space.

Physical intuition:

• A constant has no spatial variation — it’s completely “flat”
• The only frequency component is zero frequency (DC)
• Hence all the spectral energy is concentrated at $f_X = f_Y = 0$

---

Problem 3: Fourier-Bessel Transform (Goodman 2.6)

The Fourier-Bessel transform (also called the Hankel transform) is the natural transform for functions with circular symmetry. Instead of complex exponentials, it uses Bessel functions as its basis.

Definition:

$$ G(\rho) = 2\pi \int_0^{\infty} r \, g_R(r) \, J_0(2\pi r \rho) \, dr $$

where:

• $r = \sqrt{x^2 + y^2}$ is the radial coordinate
• $\rho = \sqrt{f_X^2 + f_Y^2}$ is the radial frequency
• $J_0$ is the zeroth-order Bessel function of the first kind

Why Bessel functions?

When you have circular symmetry and write the 2D Fourier transform in polar coordinates, the angular integration yields a Bessel function. It’s the natural “building block” for radially symmetric functions, just as $e^{j\omega t}$ is for 1D signals.

---

Part (a): Delta Ring

$$ \mathcal{B}\{\delta(r - r_0)\} = 2\pi r_0 J_0(2\pi r_0 \rho) $$

What is $\delta(r - r_0)$?

This is a delta ring — all the “mass” concentrated on a circle of radius $r_0$. Imagine an infinitely thin, infinitely bright ring.

The calculation:

Substituting into the Fourier-Bessel definition:

$$ G(\rho) = 2\pi \int_0^{\infty} r \, \delta(r - r_0) \, J_0(2\pi r \rho) \, dr $$

The sifting property picks out $r = r_0$:

$$ G(\rho) = 2\pi \cdot r_0 \cdot J_0(2\pi r_0 \rho) $$

Physical meaning:

A ring source produces an oscillating diffraction pattern described by $J_0$. The spacing of the oscillations depends on the ring radius $r_0$.

---

Part (b): Annular Region

For $g_R(r) = 1$ when $a \leq r \leq 1$ (and 0 otherwise):

$$ \mathcal{B}\{g_R(r)\} = \frac{J_1(2\pi\rho) - a J_1(2\pi a \rho)}{\rho} $$

What is this function?

An annulus or “washer” — a disk with a hole in the center. It’s 1 between radii $a$ and 1, and 0 elsewhere.

The key Bessel identity:

The solution uses the integration formula derived from:

$$ \frac{d}{dx}[x J_1(x)] = x J_0(x) $$

This is analogous to how $\frac{d}{dx}[\sin x] = \cos x$ helps with trigonometric integrals.

---

Part (c): Scaling Property

$$ \mathcal{B}\{g_R(ar)\} = \frac{1}{a^2} G\left(\frac{\rho}{a}\right) $$

The similarity theorem for radial functions:

Compressing the spatial function by factor $a$ (making it smaller):

• Expands the frequency spectrum by factor $a$ (spreads it out)
• Reduces the amplitude by factor $1/a^2$

Why $1/a^2$ instead of $1/a$?

In 1D, the scaling factor is $1/|a|$. In 2D with circular symmetry, you’re scaling in two dimensions, so it’s $1/a^2$.

---

Part (d): Gaussian Eigenfunction

$$ \mathcal{B}\{e^{-\pi r^2}\} = e^{-\pi \rho^2} $$

Derivation sketch:

1. The 2D Gaussian is separable: $e^{-\pi(x^2+y^2)} = e^{-\pi x^2} \cdot e^{-\pi y^2}$
2. For 1D: complete the square in the exponent of the Fourier integral
3. The resulting integral evaluates to another Gaussian
4. The 2D result follows by separability

---

Problem 4: Bessel Functions as Eigenfunctions

The claim:

$J_0(2\pi\rho_0 r)$ is an eigenfunction of any linear shift-invariant (LSI) system with circularly symmetric impulse response.

What does this mean?

When you input $J_0(2\pi\rho_0 r)$ to such a system, the output is the same function scaled by a constant:

$$ \text{output} = H(\rho_0) \cdot J_0(2\pi\rho_0 r) $$

The function passes through unchanged in shape — only its amplitude is modified.

Why is this true?

The key insight connects to Problem 3(a):

1. From 3(a): $\mathcal{B}\{\delta(r-r_0)\} = 2\pi r_0 J_0(2\pi r_0 \rho)$

2. The Fourier-Bessel transform is self-reciprocal (symmetric), so:

$$ \mathcal{B}\{J_0(2\pi\rho_0 r)\} = \frac{1}{2\pi\rho_0}\delta(\rho - \rho_0) $$

3. The Bessel function transforms to a delta ring in frequency space

4. For an LSI system, output = input × transfer function (in frequency domain):

$$ G_{\text{out}}(\rho) = H(\rho) \cdot \frac{1}{2\pi\rho_0}\delta(\rho - \rho_0) = \frac{H(\rho_0)}{2\pi\rho_0}\delta(\rho - \rho_0) $$

5. Transforming back: the output is $H(\rho_0) \cdot J_0(2\pi\rho_0 r)$

The eigenvalue:

$$ \lambda = H(\rho_0) = \mathcal{B}\{h_R(r)\}\big|_{\rho = \rho_0} $$

The eigenvalue is the transfer function evaluated at the spatial frequency $\rho_0$.

Physical interpretation:

Just as complex exponentials $e^{j2\pi fx}$ are eigenfunctions of 1D LTI systems, Bessel functions are eigenfunctions of circularly symmetric 2D systems. They’re the “natural modes” of such systems.

---

Problem 5: The Fourier Transform as a System

This problem asks us to analyze the Fourier transform itself as if it were a signal processing system.

---

Part (a): Is it Linear?

Yes, the Fourier transform is linear.

A system is linear if it satisfies superposition:

$$ \mathcal{F}\{a g_1 + b g_2\} = a \mathcal{F}\{g_1\} + b \mathcal{F}\{g_2\} $$

Proof:

Starting from the definition:

$$ \mathcal{F}\{ag_1 + bg_2\} = \iint [ag_1(x,y) + bg_2(x,y)] e^{-j2\pi(f_X x + f_Y y)} \, dx\,dy $$

By linearity of integration:

$$ = a \iint g_1 e^{-j2\pi(f_X x + f_Y y)} dx\,dy + b \iint g_2 e^{-j2\pi(f_X x + f_Y y)} dx\,dy $$$$ = a\mathcal{F}\{g_1\} + b\mathcal{F}\{g_2\} $$

---

Part (b): Does it Have a Transfer Function?

No, the Fourier transform system does not have a transfer function.

Why not?

A transfer function exists only for systems that are both:

1. Linear ✓ (proven above)
2. Shift-invariant ✗

Physical intuition:

The Fourier transform converts between fundamentally different representations (space ↔ frequency). A translation in one domain can’t map to a translation in the other because they represent different physical quantities.

Consequence:

Without shift-invariance, there’s no convolution relationship, and thus no transfer function or impulse response. Each input must be analyzed directly through the transform integral.

---

Homework 1 Summary

Problem	Core Concept
1	Delta function scaling and the Dirac comb for sampling
2	Fundamental transform pairs: rect↔sinc, separability
3	Fourier-Bessel transform for circular symmetry
4	Bessel functions as eigenfunctions of radially symmetric systems
5	Linearity vs. shift-invariance; when transfer functions exist

---

Homework 2: Computed Tomography and the Projection-Slice Theorem

Problem 1: X-ray Projections and the Projection-Slice Theorem

Consider the three-dimensional object shown in the figure below. A two-dimensional slice of this object is illuminated with a line source of x-rays, and a linear detector array is placed behind the object. The x-rays are attenuated by the material in the path of the beam; the denser the material, the less intense the beam. Assume the signal produced at a single detector in the array is equal to the integral of the density along the beam direction.

This problem introduces the mathematical foundation of Computed Tomography (CT) imaging.

---

The Physical Setup

Imagine a 3D object (like a human body). We take a 2D cross-sectional slice at some height $z_0$. This slice has a density distribution $f(x, y)$ — different tissues (bone, muscle, air) have different densities.

The x-ray imaging process:

1. A line source emits parallel x-ray beams along the $x$-direction
2. The beams pass through the object and are attenuated (weakened) by the material
3. Denser material absorbs more x-rays → weaker signal reaches the detector
4. A linear detector array on the other side measures what comes through

What is $f(x, y)$?

$f(x, y)$ is the density distribution of the cross-sectional slice — it’s what’s actually inside the object at each point $(x, y)$. Different tissues have different densities (bone is dense, muscle is medium, air is sparse). The function $f(x, y)$ assigns a density value to every point in the slice. It’s the unknown we want to reconstruct.

Symbol	What it is	Known or unknown?
$f(x,y)$	Internal density at point $(x,y)$	Unknown — can’t see inside
$p(y)$	Detector signal at position $y$	Measured — from part (a)

We measure $p(y)$. We want $f(x,y)$. The projection-slice theorem tells us how to invert this relationship.

---

Part (a): What does a single projection measure?

What is the one-dimensional signal received by the detector array for an object cross-section $f(x, y)$ located at $z_0$ when the beam is projected along the $x$-axis as in the figure? Assume the array is a continuous function of $y$.

Consider one detector element sitting at position $y$ on the detector array. The x-ray beam traveling to this detector passes through all points $(x, y)$ for $x \in (-\infty, \infty)$ at that fixed $y$ value.

The signal received is the sum (integral) of all the density values along that line:

$$ p(y) = \int_{-\infty}^{\infty} f(x, y) \, dx $$

Interpretation:

• For each $y$ position, we’re “collapsing” all the information along the $x$-direction into a single number
• The result $p(y)$ is a 1D function — it’s the projection or shadow of the 2D object onto the $y$-axis
• This operation is called the Radon transform

Example: If there’s a dense bone at some $(x_0, y_0)$, the projection $p(y_0)$ will have a higher value than if there were only soft tissue along that line.

Note on continuous vs. discrete detectors:

In the mathematical model, the detector is continuous — infinitely many beams with infinitesimal spacing. The problem even says “Assume the array is a continuous function of $y$.”

In a real CT scanner, detector elements are discrete (maybe 1000+ elements across), so there’s finite spacing. This introduces sampling considerations, but the theory treats it as continuous.

---

Spatial vs. Frequency Domain: The Interplay

Before diving into the theorem, it helps to understand the relationship between spatial and frequency representations.

Spatial domain — the “normal” view of an image or signal. Each point $(x, y)$ has a value $f(x, y)$ representing intensity, density, etc. This is what you see when you look at a photograph.

Frequency domain — a completely different representation of the same information. Each point $(\nu_x, \nu_y)$ tells you the amplitude and phase of a sinusoidal wave at that frequency.

Spatial Domain	Frequency Domain
Position $(x, y)$	Frequency $(\nu_x, \nu_y)$
Sharp edges, fine details	High frequencies (far from origin)
Smooth regions, gradual changes	Low frequencies (near origin)
Localized features	Spread across many frequencies
Periodic patterns	Concentrated at specific frequencies

Intuitive examples:

1. A sharp edge (like a boundary between black and white) requires many high-frequency sinusoids to represent — that’s why edges appear as energy far from the origin in the frequency domain.

2. A smooth gradient only needs low frequencies — the energy stays near the origin.

3. A periodic grating (stripes) appears as two bright dots in the frequency domain at $\pm$ the stripe frequency.

Why this matters for CT:

The projection-slice theorem exploits this duality. Instead of measuring $f(x,y)$ directly (which we can’t do — we can’t see inside the body), we measure projections $p(y)$. Each projection gives us a slice of $F(\nu_x, \nu_y)$. Collect enough slices → reconstruct the full frequency domain → inverse transform → recover the spatial image.

---

Part (b): The Projection-Slice Theorem

How is the one-dimensional Fourier transform of this distribution related to the two-dimensional Fourier transform of the cross-section $f(x,y)$? (Hint: Consider the two-dimensional Fourier transform along a particular line). This relationship is called the projection-slice theorem.

This is the key insight that makes CT reconstruction possible.

Question: How does the Fourier transform of the projection $p(y)$ relate to the 2D Fourier transform of $f(x,y)$?

The 2D Fourier transform of the object:

$$ F(\nu_x, \nu_y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) \, e^{-j2\pi(\nu_x x + \nu_y y)} \, dx \, dy $$

Breaking this down:

• $f(x, y)$ — the original 2D image/object (spatial domain)
• $F(\nu_x, \nu_y)$ — the transformed result (frequency domain)
• $\nu_x, \nu_y$ — spatial frequencies (cycles per unit length in x and y directions)
• $e^{-j2\pi(\nu_x x + \nu_y y)}$ — a 2D sinusoidal wave (explained below)

The double integral multiplies the image by each possible 2D wave and integrates over all space, measuring “how much” of that frequency is present. Think of it as decomposing the image into a sum of 2D sinusoidal gratings at different orientations and spacings.

Understanding the exponential term $e^{-j2\pi(\nu_x x + \nu_y y)}$:

Piece	What it is
$e$	Euler’s number (≈ 2.718)
$j$	Imaginary unit ($j = \sqrt{-1}$)
$2\pi$	Converts from “cycles” to “radians” (one cycle = $2\pi$ radians)
$\nu_x, \nu_y$	Spatial frequencies — cycles per unit length in each direction
$\nu_x x + \nu_y y$	Total phase at position $(x, y)$ — a dot product

Using Euler’s formula: $e^{-j\theta} = \cos(\theta) - j\sin(\theta)$

This is a 2D sinusoidal wave oscillating in the direction $(\nu_x, \nu_y)$. The wave is constant along lines where $\nu_x x + \nu_y y = \text{constant}$ and oscillates perpendicular to them. Why complex? Using $e^{j\theta}$ captures both amplitude and phase in one number, making the math cleaner.

$\nu_x$	$\nu_y$	Pattern
High	0	Vertical stripes
0	High	Horizontal stripes
High	High	Diagonal fine pattern
Low	Low	Slowly varying, smooth

Why do different frequencies produce different patterns?

• When $\nu_y = 0$: The phase becomes $\nu_x x$, which only depends on $x$. Lines of constant phase are vertical lines ($x = \text{const}$). The wave oscillates as you move horizontally (in $x$), creating vertical stripes.

• When $\nu_x = 0$: The phase becomes $\nu_y y$, which only depends on $y$. Lines of constant phase are horizontal lines ($y = \text{const}$). The wave oscillates as you move vertically (in $y$), creating horizontal stripes.

• When both are nonzero: Lines of constant phase are tilted at angle $\theta = \arctan(\nu_x / \nu_y)$ from horizontal. The wave oscillates perpendicular to this direction. If $\nu_x = \nu_y$, the lines are at 45°, and you get diagonal stripes. When both frequencies are high, the stripes are closely spaced.

• Low frequencies: Large wavelengths mean slow spatial variation — the pattern changes gradually across space, appearing smooth or like a gentle gradient.

The 1D Fourier transform of the projection:

$$ P(\nu_y) = \int_{-\infty}^{\infty} p(y) \, e^{-j2\pi \nu_y y} \, dy $$

Breaking this down:

• $p(y)$ — the 1D projection signal from the detector array
• $P(\nu_y)$ — the transformed result (frequency domain)
• $\nu_y$ — spatial frequency along the $y$-direction (how many cycles per unit length as you move in $y$; e.g., $\nu_y = 5$ cycles/cm means 5 oscillations per cm, wavelength = 0.2 cm)
• $e^{-j2\pi \nu_y y}$ — a 1D sinusoidal wave at frequency $\nu_y$

This integral decomposes the projection into its frequency components — how much of each sinusoidal frequency is present in the 1D signal.

Derivation: Substitute $p(y) = \int f(x,y) \, dx$ into the 1D transform:

$$ P(\nu_y) = \int_{-\infty}^{\infty} \left[ \int_{-\infty}^{\infty} f(x, y) \, dx \right] e^{-j2\pi \nu_y y} \, dy $$

Since $e^{-j2\pi \nu_y y}$ doesn’t depend on $x$, we can pull it inside the inner integral (constants pass through integrals). This gives us a single double integral:

$$ P(\nu_y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x, y) \, e^{-j2\pi \nu_y y} \, dx \, dy $$

Where does “$0 \cdot x$” come from? Compare this to the 2D Fourier transform:

$$ F(\nu_x, \nu_y) = \int \int f(x,y) \, e^{-j2\pi(\nu_x x + \nu_y y)} \, dx \, dy $$

If we set $\nu_x = 0$, we get $e^{-j2\pi(0 \cdot x + \nu_y y)} = e^{-j2\pi \nu_y y}$ — exactly what we have! So we write the “$0 \cdot x$” to make the pattern match obvious:

$$ P(\nu_y) = \int \int f(x, y) \, e^{-j2\pi(0 \cdot x + \nu_y \cdot y)} \, dx \, dy = F(0, \nu_y) $$

Notice that this is exactly $F(\nu_x, \nu_y)$ evaluated at $\nu_x = 0$:

$$ \boxed{P(\nu_y) = F(0, \nu_y)} $$

Visually:

• The object $f(x,y)$ lives in the spatial domain (x-y plane)
• Its 2D Fourier transform $F(\nu_x, \nu_y)$ lives in the frequency domain ($\nu_x$-$\nu_y$ plane)
• A projection along $x$ gives you the vertical line $\nu_x = 0$ in the frequency domain

---

Part (c): Reconstructing the Full Image

Locate the region of the Fourier transform plane $\mathcal{F}\{f(x, y)\}$ which can be recovered from projecting a line of x-rays down the x-axis from part (b). What region can be recovered from a $y$-axis projection? Can you describe a technique for recovering the entire cross-sectional image $f(x,y)$ from many projections? If you can, you have just invented a method of computer tomography.

What do we get from each projection?

Projection Direction	Spatial Result	Frequency Slice
Along $x$-axis	$p(y) = \int f \, dx$	Vertical line: $F(0, \nu_y)$
Along $y$-axis	$q(x) = \int f \, dy$	Horizontal line: $F(\nu_x, 0)$
At angle $\theta$	projection at angle $\theta$	Radial line at angle $\theta$ through origin

Connection to spatial frequency patterns:

A projection “sees” variations perpendicular to its direction, but is blind to variations parallel to it:

Pattern	Dominant frequencies	Captured by $x$-projection?	Captured by $y$-projection?
Vertical stripes	High $\nu_x$, low $\nu_y$	No (need $\nu_x \neq 0$)	Yes
Horizontal stripes	Low $\nu_x$, high $\nu_y$	Yes	No (need $\nu_y \neq 0$)
Diagonal/checkerboard	High $\nu_x$ and $\nu_y$	Partially	Partially
Smooth gradient	Low both	Yes	Yes

This is why you need projections from many angles — each one captures a different slice of the frequency content.

The reconstruction idea:

With just two projections (along $x$ and $y$), you only know two lines in the Fourier domain — that’s not enough to reconstruct the image.

But if you take many projections at different angles $\theta \in [0, \pi)$:

1. Each projection gives you one radial “spoke” through the origin in the Fourier domain
2. With enough projections, you fill in the entire 2D Fourier plane
3. Apply the inverse 2D Fourier transform to recover $f(x,y)$

This is Computed Tomography!

---

Problem 2: Equivalent Area and Equivalent Bandwidth (Goodman 2.5)

The “equivalent area” $\Delta_{xy}$ of a function $g(x, y)$ can be defined by

$$ \Delta_{xy} = \frac{\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} g(x, y) \, dx \, dy}{g(0, 0)}, $$

while the “equivalent bandwidth” $\Delta_{f_x f_y}$ of $g$ is defined in terms of its transform $G$ by

$$ \Delta_{f_x f_y} = \frac{\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G(f_x, f_y) \, df_x \, df_y}{G(0, 0)}. $$

Show that $\Delta_{xy} \cdot \Delta_{f_x f_y} = 1$.

This problem reveals a beautiful reciprocal relationship between a function’s spatial extent and its frequency content — a manifestation of the uncertainty principle.

---

The Definitions

Equivalent area of a function $g(x, y)$:

$$ \Delta_{xy} = \frac{\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} g(x, y) \, dx \, dy}{g(0, 0)} $$

Equivalent bandwidth of its Fourier transform $G(f_x, f_y)$:

$$ \Delta_{f_x f_y} = \frac{\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G(f_x, f_y) \, df_x \, df_y}{G(0, 0)} $$

Notation convention: Lowercase $g$ denotes the spatial domain function; uppercase $G$ denotes its Fourier transform (frequency domain). This is standard: the transform of $g$ is $G$, the transform of $f$ is $F$, etc.

---

What Do These Quantities Mean?

Equivalent area $\Delta_{xy}$:

• The numerator is the total “volume” under the function $g(x,y)$
• Dividing by $g(0,0)$ (the value at the origin) converts volume → area:
  • Think of $g(x,y)$ as a height above the $(x,y)$ plane
  • $\int \int g \, dx\,dy$ = volume under the surface $z = g(x,y)$
  • $g(0,0)$ = height at the origin
  • Dividing: volume / height = area

- Result: the area of a rectangle with height $g(0,0)$ that has the same volume as $g$ - **Intuition**: How "spread out" is the function in space? - **Note**: This definition assumes $g$ is centered at the origin. For symmetric functions like Gaussians, the origin is also the peak — but that's not guaranteed in general.

Equivalent bandwidth $\Delta_{f_x f_y}$:

• Same idea, but for the spectrum $G(f_x, f_y)$
• Intuition: How “spread out” is the function in frequency?

Quantity	Domain	Measures
$\Delta_{xy}$	Spatial	Effective spatial extent
$\Delta_{f_x f_y}$	Frequency	Effective frequency extent

Geometric interpretation:

The equivalent width $\Delta_x$ (or area $\Delta_{xy}$ in 2D) answers the question: What size rectangle, with height equal to the peak value, has the same area/volume as the original function?

---

The Proof

Goal: Show that $\Delta_{xy} \cdot \Delta_{f_x f_y} = 1$

Step 1: Use the Fourier transform at the origin

Look at the definition of $\Delta_{xy}$:

$$ \Delta_{xy} = \frac{\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} g(x,y) \, dx \, dy}{g(0,0)} $$

The numerator is the integral of $g$ over all space. We need to relate this to the Fourier transform $G$. The question is: does this integral appear anywhere in the Fourier transform?

Start with the definition of the 2D Fourier transform:

$$ G(f_x, f_y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} g(x, y) \, e^{-j2\pi(f_x x + f_y y)} \, dx \, dy $$

The trick is to evaluate this at the origin of the frequency domain ($f_x = 0$, $f_y = 0$), because this makes the exponential disappear:

$$ G(0, 0) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} g(x, y) \, e^{-j2\pi(0 \cdot x + 0 \cdot y)} \, dx \, dy $$

The exponent becomes zero:

$$ e^{-j2\pi(0 \cdot x + 0 \cdot y)} = e^{0} = 1 $$

So the exponential simply disappears, leaving:

$$ G(0, 0) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} g(x, y) \, dx \, dy $$

This is exactly the numerator of $\Delta_{xy}$! So we can write:

$$ \Delta_{xy} = \frac{G(0,0)}{g(0,0)} $$

Physical meaning: The value of the spectrum at zero frequency (the “DC component”) equals the total integral of the function over all space. This makes intuitive sense: the zero-frequency component represents the average or total “mass” of the signal, with no oscillation.

Step 2: Use the inverse Fourier transform at the origin

Now look at the definition of $\Delta_{f_x f_y}$:

$$ \Delta_{f_x f_y} = \frac{\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G(f_x, f_y) \, df_x \, df_y}{G(0,0)} $$

The numerator is the integral of $G$ over all frequencies. We need to relate this to the spatial function $g$. The question is: does this integral appear anywhere in the inverse Fourier transform?

Start with the definition of the inverse Fourier transform:

$$ g(x, y) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G(f_x, f_y) \, e^{j2\pi(f_x x + f_y y)} \, df_x \, df_y $$

The trick is to evaluate this at the origin of the spatial domain ($x = 0$, $y = 0$), because this makes the exponential disappear:

$$ g(0, 0) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G(f_x, f_y) \, e^{j2\pi(0 \cdot f_x + 0 \cdot f_y)} \, df_x \, df_y $$

The exponent becomes zero:

$$ e^{j2\pi(0 \cdot f_x + 0 \cdot f_y)} = e^{0} = 1 $$

So the exponential simply disappears, leaving:

$$ g(0, 0) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} G(f_x, f_y) \, df_x \, df_y $$

This is exactly the numerator of $\Delta_{f_x f_y}$! So we can write:

$$ \Delta_{f_x f_y} = \frac{g(0,0)}{G(0,0)} $$

Physical meaning: The value of the spatial function at the origin equals the total integral of the spectrum over all frequencies. This is the “dual” of Step 1 — the Fourier transform and its inverse have symmetric properties.

Step 3: Substitute into the definitions

From Step 1:

$$ \Delta_{xy} = \frac{\int \int g(x,y) \, dx \, dy}{g(0,0)} = \frac{G(0,0)}{g(0,0)} $$

From Step 2:

$$ \Delta_{f_x f_y} = \frac{\int \int G(f_x, f_y) \, df_x \, df_y}{G(0,0)} = \frac{g(0,0)}{G(0,0)} $$

Step 4: Multiply

$$ \Delta_{xy} \cdot \Delta_{f_x f_y} = \frac{G(0,0)}{g(0,0)} \cdot \frac{g(0,0)}{G(0,0)} = 1 $$$$ \boxed{\Delta_{xy} \cdot \Delta_{f_x f_y} = 1} $$

---

Physical Interpretation

If a function is…	Then its spectrum is…
Narrow in space (small $\Delta_{xy}$)	Wide in frequency (large $\Delta_{f_x f_y}$)
Wide in space (large $\Delta_{xy}$)	Narrow in frequency (small $\Delta_{f_x f_y}$)

Examples:

1. Gaussian: $g(x) = e^{-\pi x^2}$ transforms to $G(f) = e^{-\pi f^2}$

   • Both have the same shape — the Gaussian is “balanced” between domains
   • This is why Gaussians achieve the minimum uncertainty product

2. Rectangle: $\text{rect}(x)$ transforms to $\text{sinc}(f)$

   • Compact in space → infinite extent in frequency (though decaying)

3. Delta function: $\delta(x)$ transforms to $1$ (constant)

   • Infinitely narrow in space → infinitely wide (flat) in frequency

---

Connection to Optics

---

Problem 3: Periodic Functions and the Comb Function (Goodman 2.11/2.13)

$$ p(x, y) = g(x, y) * \left[ \text{comb}\left(\frac{x}{X}\right) \text{comb}\left(\frac{y}{Y}\right) \right] $$

defines a periodic function, with period $X$ in the $x$ direction and period $Y$ in the $y$ direction.

This problem explores how periodic functions arise from convolution with the Dirac comb, and how their Fourier transforms become discrete.

---

The Setup

We’re given a function defined as:

$$ p(x, y) = g(x, y) * \left[ \text{comb}\left(\frac{x}{X}\right) \text{comb}\left(\frac{y}{Y}\right) \right] $$

with period 𝑋 in the 𝑥 direction and period 𝑌 in the 𝑦 direction, and where $*$ denotes convolution.

What does this create?

The Dirac comb is an infinite train of delta functions:

$$ \text{comb}(x) = \sum_{n=-\infty}^{\infty} \delta(x - n) $$

Convolving any function $g$ with a delta function at position $x_0$ produces a copy of $g$ shifted to that position:

$$ g(x) * \delta(x - x_0) = g(x - x_0) $$

Therefore, convolving with the 2D Dirac comb creates infinite copies of $g$ at every point $(nX, mY)$ for all integers $n, m$. The result is a doubly-periodic function with period $X$ in $x$ and period $Y$ in $y$.

---

Part (a): Fourier Transform of a Periodic Function

Show that the Fourier transform of $p$ can be written

$$ P(f_x, f_y) = \sum_{n=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} G\left(\frac{n}{X}, \frac{m}{Y}\right) \delta\left(f_x - \frac{n}{X}, f_y - \frac{m}{Y}\right) $$

where $G$ is the Fourier transform of $g$.

Goal: Show that

$$ P(f_x, f_y) = \sum_{n=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} G\left(\frac{n}{X}, \frac{m}{Y}\right) \delta\left(f_x - \frac{n}{X}, f_y - \frac{m}{Y}\right) $$

where $G$ is the Fourier transform of $g$.

---

Step 1: Define the Dirac comb

$$ \text{comb}(x) = \sum_{n=-\infty}^{\infty} \delta(x - n) $$

Step 2: Scale the Dirac comb

Substituting $x/X$ into the Dirac comb definition:

$$ \text{comb}\left(\frac{x}{X}\right) = \sum_{n=-\infty}^{\infty} \delta\left(\frac{x}{X} - n\right) $$

Rewrite each delta’s argument: $\displaystyle\frac{x}{X} - n = \frac{x - nX}{X} = \frac{1}{X}(x - nX)$

Apply the scaling property $\delta(u/a) = |a| \delta(u)$ with $u = x - nX$ and $a = X$:

$$ \delta\left(\frac{x - nX}{X}\right) = |X| \delta(x - nX) $$

Therefore (assuming $X > 0$):

$$ \text{comb}\left(\frac{x}{X}\right) = \sum_{n=-\infty}^{\infty} X \delta(x - nX) = X \sum_{n=-\infty}^{\infty} \delta(x - nX) $$

Step 3: Fourier transform of the Dirac comb

The Dirac comb is its own Fourier transform: $\mathcal{F}\{\text{comb}(x)\} = \text{comb}(f)$.

The scaling theorem states: if $\mathcal{F}\{g(x)\} = G(f)$, then $\mathcal{F}\{g(ax)\} = \frac{1}{|a|} G\left(\frac{f}{a}\right)$.

Here $\text{comb}(x/X) = \text{comb}\left(\frac{1}{X} \cdot x\right)$, so $a = \frac{1}{X}$:

$$ \mathcal{F}\left\{\text{comb}\left(\frac{1}{X} \cdot x\right)\right\} = \frac{1}{|1/X|} \cdot \text{comb}\left(\frac{f_x}{1/X}\right) = |X| \cdot \text{comb}(Xf_x) $$

Expand the comb: $\text{comb}(Xf_x) = \sum_n \delta(Xf_x - n)$

Apply delta scaling to each term:

$$ \delta(Xf_x - n) = \delta\left(X\left(f_x - \frac{n}{X}\right)\right) = \frac{1}{|X|}\delta\left(f_x - \frac{n}{X}\right) $$

The $|X|$ and $\frac{1}{|X|}$ cancel:

$$ \mathcal{F}\left\{\text{comb}\left(\frac{x}{X}\right)\right\} = |X| \cdot \sum_n \frac{1}{|X|}\delta\left(f_x - \frac{n}{X}\right) = \sum_{n=-\infty}^{\infty} \delta\left(f_x - \frac{n}{X}\right) $$

Step 4: Extend to 2D

For the separable 2D Dirac comb:

$$ \mathcal{F}\left\{\text{comb}\left(\frac{x}{X}\right) \text{comb}\left(\frac{y}{Y}\right)\right\} = \sum_{n=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} \delta\left(f_x - \frac{n}{X}\right) \delta\left(f_y - \frac{m}{Y}\right) $$

Step 5: Apply the convolution theorem

Convolution in space becomes multiplication in frequency:

$$ P(f_x, f_y) = G(f_x, f_y) \cdot \sum_{n,m} \delta\left(f_x - \frac{n}{X}, f_y - \frac{m}{Y}\right) $$

Step 6: Apply the sifting property

Multiplying $G$ by a delta evaluates $G$ at the delta’s location:

$$ \boxed{P(f_x, f_y) = \sum_{n=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} G\left(\frac{n}{X}, \frac{m}{Y}\right) \delta\left(f_x - \frac{n}{X}, f_y - \frac{m}{Y}\right)} $$

---

Physical Interpretation

Spatial Domain	Frequency Domain
Periodic with period $(X, Y)$	Discrete with spacing $(1/X, 1/Y)$
Continuous (smooth tile)	Discrete (isolated spikes)
Infinite extent (repeats forever)	Finite number of significant harmonics

This is the 2D version of the Fourier series: periodic functions can be represented as sums of harmonics, and the transform reveals which harmonics are present.

---

Deep Dive: Delta Function Scaling

---

Part (b): Specific Example

Sketch the function $p(x, y)$ when

$$ g(x, y) = \text{rect}\left(\frac{2x}{X}\right) \text{rect}\left(\frac{2y}{Y}\right) $$

and find the corresponding Fourier transform $P(f_x, f_y)$.

Given: $g(x, y) = \text{rect}\left(\frac{2x}{X}\right) \text{rect}\left(\frac{2y}{Y}\right)$

This is a rectangle centered at the origin with width $X/2$ and height $Y/2$ — exactly one quarter of the unit cell area.

The periodic function $p(x, y)$:

Convolving with the Dirac comb creates copies at every $(nX, mY)$. The result is a grid of rectangles, each of size $\frac{X}{2} \times \frac{Y}{2}$, with centers separated by $X$ and $Y$.

Finding $G(f_x, f_y)$:

The rectangle function transforms to a sinc:

$$ \mathcal{F}\{\text{rect}(x)\} = \text{sinc}(f) $$

With scaling:

$$ \mathcal{F}\left\{\text{rect}\left(\frac{2x}{X}\right)\right\} = \frac{X}{2} \text{sinc}\left(\frac{X f_x}{2}\right) $$

For the separable 2D case:

$$ G(f_x, f_y) = \frac{X}{2} \text{sinc}\left(\frac{X f_x}{2}\right) \cdot \frac{Y}{2} \text{sinc}\left(\frac{Y f_y}{2}\right) = \frac{XY}{4} \text{sinc}\left(\frac{X f_x}{2}\right) \text{sinc}\left(\frac{Y f_y}{2}\right) $$

The transform $P(f_x, f_y)$:

Using the result from part (a):

$$ P(f_x, f_y) = \frac{XY}{4} \sum_{n,m} \text{sinc}\left(\frac{n}{2}\right) \text{sinc}\left(\frac{m}{2}\right) \delta\left(f_x - \frac{n}{X}, f_y - \frac{m}{Y}\right) $$

Which harmonics survive?

The sinc function has zeros at nonzero integers:

$$ \text{sinc}(k) = 0 \quad \text{for } k = \pm 1, \pm 2, \pm 3, \ldots $$

But $\text{sinc}(n/2)$ is zero only when $n/2$ is a nonzero integer, i.e., when $n$ is an even nonzero integer ($n = \pm 2, \pm 4, \ldots$).

$n$	$\text{sinc}(n/2)$
0	1
$\pm 1$	$2/\pi \approx 0.637$
$\pm 2$	0
$\pm 3$	$-2/(3\pi) \approx -0.212$
$\pm 4$	0

Only odd harmonics (and the DC term) have nonzero coefficients. This is characteristic of a 50% duty cycle pattern.

---

Problem 4: The Equivalent Object Theorem (Goodman 2.13/2.15)

The input to a certain imaging system is an object complex field distribution $U_o(x, y)$ of unlimited spatial frequency content, while the output of the system is an image field distribution $U_i(x, y)$. The imaging system can be assumed to act as a linear, invariant lowpass filter with a transfer function that is identically zero outside the region $|f_x| \leq B_X / 2, |f_y| \leq B_Y / 2$ in the frequency domain. Show that there exists an “equivalent” object $U_o'(x, y)$ consisting of a rectangular array of point sources that produces exactly the same image $U_i$ as does the true object $U_o$, and that the field distribution across the equivalent object can be written

$$ U_o'(x, y) = \sum_{n=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} \left[ \int_{\xi=-\infty}^{\infty} \int_{\eta=-\infty}^{\infty} U_o(\xi, \eta) \, \text{sinc}(n - B_X \xi) \, \text{sinc}(m - B_Y \eta) \, d\xi \, d\eta \right] \delta\left(x - \frac{n}{B_X}, y - \frac{m}{B_Y}\right). $$

This problem proves a remarkable result: for any bandlimited imaging system, there exists an equivalent object made entirely of point sources that produces exactly the same image.

---

The Setup

The imaging system:

• Input: object field $U_o(x, y)$ with unlimited spatial frequency content
• Output: image field $U_i(x, y)$
• The system acts as a lowpass filter with cutoff frequencies $B_X/2$ and $B_Y/2$
• Transfer function: $H(f_x, f_y) = 1$ inside the passband, $0$ outside

The claim: There exists an equivalent object $U_o'(x, y)$ consisting entirely of point sources (delta functions) on a rectangular grid that produces exactly the same image.

---

Why Does This Work?

This is analogous to how two audio signals that differ only in ultrasonic frequencies will sound identical through speakers that can’t reproduce ultrasonics.

---

The Derivation

Step 1: Express the image in terms of bandlimited content

The image is determined by convolving the object spectrum with the system’s transfer function:

$$ U_i(x, y) = \mathcal{F}^{-1}\{H(f_x, f_y) \cdot \mathcal{U}_o(f_x, f_y)\} $$

Since $H$ is zero outside the passband, only the portion of $\mathcal{U}_o$ inside the region $|f_x| \leq B_X/2$, $|f_y| \leq B_Y/2$ contributes to the image.

Define the bandlimited spectrum:

$$ \tilde{\mathcal{U}}_o(f_x, f_y) = \begin{cases} \mathcal{U}_o(f_x, f_y) & |f_x| \leq B_X/2 \text{ and } |f_y| \leq B_Y/2 \\ 0 & \text{otherwise} \end{cases} $$

Step 2: Represent the bandlimited spectrum as a Fourier series

Since $\tilde{\mathcal{U}}_o$ is nonzero only within a finite rectangle, we can expand it as a Fourier series over that region:

$$ \tilde{\mathcal{U}}_o(f_x, f_y) = \sum_{n=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} c_{nm} \, e^{j2\pi\left(\frac{n f_x}{B_X} + \frac{m f_y}{B_Y}\right)} $$

The coefficients are:

$$ c_{nm} = \frac{1}{B_X B_Y} \int_{-B_X/2}^{B_X/2} \int_{-B_Y/2}^{B_Y/2} \tilde{\mathcal{U}}_o(f_x, f_y) \, e^{-j2\pi\left(\frac{n f_x}{B_X} + \frac{m f_y}{B_Y}\right)} df_x \, df_y $$

Step 3: Relate coefficients to spatial integrals

This coefficient integral can be rewritten using the convolution theorem. The result involves the original object convolved with sinc functions:

$$ c_{nm} = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} U_o(\xi, \eta) \, \text{sinc}(n - B_X \xi) \, \text{sinc}(m - B_Y \eta) \, d\xi \, d\eta $$

Step 4: Construct the equivalent object

The inverse Fourier transform of $e^{j2\pi(f_x n/B_X + f_y m/B_Y)}$ is a delta function at position $(n/B_X, m/B_Y)$.

Therefore, the equivalent object is:

$$ \boxed{U_o'(x, y) = \sum_{n=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} \left[ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} U_o(\xi, \eta) \, \text{sinc}(n - B_X \xi) \, \text{sinc}(m - B_Y \eta) \, d\xi \, d\eta \right] \delta\left(x - \frac{n}{B_X}, y - \frac{m}{B_Y}\right)} $$

---

Physical Interpretation

Property	Equivalent Object
Structure	Grid of point sources (delta functions)
Grid spacing	$1/B_X$ in $x$, $1/B_Y$ in $y$
Amplitudes	Sinc-weighted integrals of $U_o$
Image produced	Identical to original object

Why does the spacing matter?

The grid spacing is $1/B_X$ and $1/B_Y$ — exactly the Nyquist sampling intervals for a signal bandlimited to $B_X/2$ and $B_Y/2$. This is no coincidence: the sampling theorem tells us that a bandlimited signal is completely determined by samples at the Nyquist rate.

---

Problem 5: Angular Spectrum of Apertures (Goodman 3.5)

Assuming unit-amplitude normally incident plane-wave illumination, find the angular spectrum of

This problem introduces the angular spectrum representation — a powerful method for analyzing diffraction by decomposing a wavefront into plane waves traveling in different directions.

---

Background: The Angular Spectrum

When a plane wave illuminates an aperture, the field immediately after the aperture is:

$$ U(x, y, 0) = t(x, y) $$

where $t(x, y)$ is the transmission function of the aperture.

The angular spectrum is simply the 2D Fourier transform of this field:

$$ A(f_x, f_y) = \mathcal{F}\{U(x, y, 0)\} = \mathcal{F}\{t(x, y)\} $$

---

Part (a): Circular Aperture

A circular aperture of diameter $d$.

The aperture: A circular opening of diameter $d$.

Transmission function:

$$ t(x, y) = \text{circ}\left(\frac{2r}{d}\right) $$

where $r = \sqrt{x^2 + y^2}$ and the circ function is defined as:

$$ \text{circ}(\rho) = \begin{cases} 1 & \rho \leq 1 \\ 0 & \rho > 1 \end{cases} $$

Finding the angular spectrum:

For circularly symmetric functions, the 2D Fourier transform reduces to the Fourier-Bessel (Hankel) transform:

$$ A(\rho_f) = 2\pi \int_0^{\infty} r \, t(r) \, J_0(2\pi r \rho_f) \, dr $$

where $\rho_f = \sqrt{f_x^2 + f_y^2}$ is the radial spatial frequency.

For the circ function, the standard transform pair (Goodman Table 2.1) gives:

$$ \mathcal{F}\left\{\text{circ}\left(\frac{r}{a}\right)\right\} = \pi a^2 \frac{2J_1(2\pi a \rho_f)}{2\pi a \rho_f} $$

With $a = d/2$ (the radius):

$$ A(f_x, f_y) = \frac{\pi d^2}{4} \cdot \frac{2J_1(\pi d \rho_f)}{\pi d \rho_f} $$

This can be written more compactly using the jinc function:

$$ \text{jinc}(x) = \frac{2J_1(\pi x)}{\pi x} $$

Therefore:

$$ \boxed{A(f_x, f_y) = \frac{\pi d^2}{4} \, \text{jinc}(d\rho_f)} $$

where $\rho_f = \sqrt{f_x^2 + f_y^2}$.

---

Properties of the Result

Physical interpretation:

The angular spectrum tells us how the light from the aperture is distributed among different propagation directions:

• Central peak: Most energy travels straight through (low spatial frequencies, small angles)
• Rings: Some energy diffracts to larger angles (higher spatial frequencies)
• Zeros: At certain angles, there’s complete destructive interference

The first zero occurs at $\rho_f = 1.22/d$, which corresponds to the angle $\theta = \arcsin(1.22\lambda/d)$ — the Rayleigh criterion for resolution.

---

Part (b): Circular Opaque Disk

A circular opaque disk of diameter $d$.

The obstacle: A circular opaque disk of diameter $d$.

Transmission function:

$$ t(x, y) = 1 - \text{circ}\left(\frac{2r}{d}\right) $$

This is the complement of the circular aperture — transparent everywhere except within the disk.

Finding the angular spectrum:

By linearity of the Fourier transform:

$$ A_{\text{disk}}(f_x, f_y) = \mathcal{F}\{1\} - \mathcal{F}\left\{\text{circ}\left(\frac{2r}{d}\right)\right\} $$

The Fourier transform of a constant is a delta function at the origin:

$$ \mathcal{F}\{1\} = \delta(f_x, f_y) $$

Therefore:

$$ \boxed{A_{\text{disk}}(f_x, f_y) = \delta(f_x, f_y) - \frac{\pi d^2}{4} \, \text{jinc}(d\rho_f)} $$

---

Babinet’s Principle

Physical interpretation:

Component	What it represents
$\delta(f_x, f_y)$	The unobstructed plane wave — light that would pass if nothing were there
$-\frac{\pi d^2}{4} \text{jinc}(d\rho_f)$	The diffracted wave — light removed by the disk

Remarkable consequence:

Except for the central delta function (the DC component, representing the undiffracted beam), the angular spectrum of the opaque disk has the same magnitude as that of the circular aperture — they differ only in sign!

This means:

• The diffraction patterns (intensity distributions) are identical away from the center
• An opaque disk and a circular aperture of the same size produce the same fringe pattern
• Only the on-axis intensity differs (the disk blocks the direct beam; the aperture transmits it)

---

Homework 2 Summary

Problem	Key Concept
1(a)	A projection is a line integral — it collapses one dimension
1(b)	Projection-slice theorem: 1D FT of projection = slice of 2D FT
1(c)	Many projections at different angles → fill Fourier space → inverse FT recovers image
2	Equivalent area × equivalent bandwidth = 1 (uncertainty principle)
3	Periodic functions have discrete spectra; Dirac comb convolution creates periodicity
4	Bandlimited systems can replace continuous objects with point source arrays
5	Angular spectrum method; circ ↔ jinc transform; Babinet’s principle

This elegant mathematical relationship — that projections in real space correspond to slices in Fourier space — is why CT scanners rotate around the patient, acquiring projections from many angles to reconstruct the internal structure. The area-bandwidth product of Problem 2 quantifies the fundamental tradeoff between spatial and spectral localization. The angular spectrum representation of Problem 6 provides the foundation for understanding diffraction from arbitrary apertures.

---

Homework 3: Fresnel and Fraunhofer Diffraction

What Is Diffraction?

When light encounters an obstacle or aperture, it doesn’t simply cast a sharp geometric shadow — it bends around edges and spreads out. This phenomenon is called diffraction, and it’s a direct consequence of the wave nature of light.

Two regimes describe diffraction depending on the observation distance:

Regime	Condition	Character
Fresnel (near-field)	Observation at moderate distance	Complex pattern; both amplitude and phase vary
Fraunhofer (far-field)	Observation at large distance ($z \gg d^2/\lambda$)	Pattern is the Fourier transform of the aperture

---

Problem 1: Cascaded Fresnel Propagation (Goodman 4-4/4-5)

Fresnel propagation over a sequence of successive distances $z_1, z_2, \ldots, z_n$ must be equivalent to Fresnel propagation over the single distance $z = z_1 + z_2 + \cdots + z_n$. Find a simple proof that this is the case. Repeat for propagation without the Fresnel approximation.

This problem establishes a fundamental consistency property of wave propagation: breaking a journey into stages doesn’t change the result. Light propagating from $A$ to $B$ directly must arrive in the same state as light propagating $A \to C_1 \to C_2 \to \cdots \to B$ through intermediate planes.

Background: Transfer Functions in Linear Systems

In linear systems theory (Goodman 2.3.2), the output of a system is related to its input through the transfer function $H(f)$:

$$ \text{Output spectrum} = H(f) \cdot \text{Input spectrum} $$

When two systems are cascaded (output of the first feeds into the second), their transfer functions multiply:

$$ H_{\text{total}}(f) = H_1(f) \cdot H_2(f) $$

Free-space propagation is a linear, shift-invariant operation on optical fields. Each propagation distance has its own transfer function, and successive propagation stages cascade.

---

Part (a): With the Fresnel Approximation

The Fresnel approximation (Goodman 4.2) assumes that propagation angles are small (paraxial regime). Under this approximation, the transfer function for propagation over distance $z$ is (Goodman 4.2.4):

$$ H(f_X, f_Y; z) = e^{j \frac{2\pi z}{\lambda}} \, e^{-j\pi\lambda z(f_X^2 + f_Y^2)} $$

Factor	Physical meaning
$e^{j 2\pi z/\lambda}$	On-axis phase accumulation — the phase picked up by a plane wave traveling straight along $z$
$e^{-j\pi\lambda z(f_X^2 + f_Y^2)}$	Quadratic phase — off-axis components acquire additional phase proportional to $f_X^2 + f_Y^2$

Cascading $n$ stages:

For $n$ successive propagation distances $z_1, z_2, \ldots, z_n$, the total transfer function is the product:

$$ H_{\text{total}} = \prod_{k=1}^{n} H(f_X, f_Y; z_k) = \prod_{k=1}^{n} e^{j \frac{2\pi z_k}{\lambda}} \, e^{-j\pi\lambda z_k(f_X^2 + f_Y^2)} $$

Since we’re multiplying exponentials, the exponents add:

$$ H_{\text{total}} = e^{j \frac{2\pi}{\lambda} \sum_{k=1}^{n} z_k} \, e^{-j\pi\lambda \left(\sum_{k=1}^{n} z_k\right)(f_X^2 + f_Y^2)} $$

Defining $z = z_1 + z_2 + \cdots + z_n$:

$$ H_{\text{total}} = e^{j \frac{2\pi z}{\lambda}} \, e^{-j\pi\lambda z(f_X^2 + f_Y^2)} = H(f_X, f_Y; z) $$

This is exactly the transfer function for a single propagation over the total distance $z$. The proof works because the Fresnel transfer function’s exponent is linear in $z$, so the exponents simply add when transfer functions are multiplied.

---

Part (b): Without the Fresnel Approximation (Exact)

The exact transfer function for free-space propagation (no paraxial approximation) is:

$$ H(f_X, f_Y; z) = e^{j \frac{2\pi z}{\lambda} \sqrt{1 - (\lambda f_X)^2 - (\lambda f_Y)^2}} $$

This comes directly from the Helmholtz equation (Goodman 3.3.1) via the angular spectrum propagation result (Goodman 3.10.2, 3.10.4). Each spatial frequency $(f_X, f_Y)$ corresponds to a plane wave with propagation direction $(\alpha, \beta, \gamma)$ where $\alpha = \lambda f_X$, $\beta = \lambda f_Y$, and $\gamma = \sqrt{1 - \alpha^2 - \beta^2}$. The exponential gives the phase accumulated by that plane wave over distance $z$.

The same argument applies:

The exponent is still linear in $z$:

$$ \prod_{k=1}^{n} e^{j \frac{2\pi z_k}{\lambda} \sqrt{1 - (\lambda f_X)^2 - (\lambda f_Y)^2}} = e^{j \frac{2\pi}{\lambda} \left(\sum_{k=1}^{n} z_k\right) \sqrt{1 - (\lambda f_X)^2 - (\lambda f_Y)^2}} = H(f_X, f_Y; z) $$

The proof is identical in structure. As long as the transfer function’s exponent is linear in $z$, cascading works by addition of exponents. This holds for both the exact and approximate transfer functions.

---

Problem 2: Double-Slit Fraunhofer Diffraction (Goodman 4-7/4-9)

Assuming unit-amplitude, normally incident plane-wave illumination:

(a) Find the intensity distribution in the Fraunhofer diffraction pattern of the double-slit aperture shown in the figure below.

(b) Sketch normalized cross sections of this pattern that appear along the $x$ and $y$ axes in the observation plane.

This is the classic Young’s double-slit experiment analyzed rigorously using Fourier optics. The result reveals how single-slit diffraction (sinc envelope) and two-slit interference (cosine fringes) combine to produce the observed pattern.

---

Background: Fraunhofer Diffraction

In the Fraunhofer (far-field) regime (Goodman 4.3), the diffracted field is proportional to the Fourier transform of the aperture transmittance:

$$ U(x, y) \propto \mathcal{F}\{t_A(\xi, \eta)\}\bigg|_{f_X = x/(\lambda z),\; f_Y = y/(\lambda z)} $$

The intensity is $I(x, y) = |U(x, y)|^2$.

---

The Convolution Approach

The aperture consists of two identical rectangular slits (each $X$ wide and $Y$ tall), separated by a center-to-center distance $\Delta$ along the $\eta$ axis.

Rather than integrating over this aperture directly, there is an elegant approach: express the aperture as a convolution of a single slit with a pair of delta functions:

$$ t_A(\xi, \eta) = \left[\operatorname{rect}\left(\frac{\xi}{X}\right) \operatorname{rect}\left(\frac{\eta}{Y}\right)\right] * \left[\delta(\xi)\left(\delta\left(\eta - \frac{\Delta}{2}\right) + \delta\left(\eta + \frac{\Delta}{2}\right)\right)\right] $$

Component	What it represents
$\operatorname{rect}(\xi/X)\operatorname{rect}(\eta/Y)$	A single rectangular slit centered at the origin
$\delta(\eta - \Delta/2) + \delta(\eta + \Delta/2)$	Two point sources at $\eta = \pm\Delta/2$
$*$ (convolution)	“Stamp” a copy of the rectangle at each delta location

---

Part (a): The Intensity Distribution

The Fourier transform of each factor is computed separately.

Rectangle → Sinc: The 1D rect → sinc transform is a standard pair from Goodman Table 2.1 (Section 2.1.6): $\mathcal{F}\{\operatorname{rect}(x/X)\} = X\operatorname{sinc}(Xf_X)$. The 2D version follows from separability (Goodman 2.1.4): because $\operatorname{rect}(\xi/X)\operatorname{rect}(\eta/Y)$ is a product of a function of $\xi$ alone times a function of $\eta$ alone, the 2D Fourier transform is just the product of the two 1D pairs from the table:

$$ \mathcal{F}\left\{\operatorname{rect}\left(\frac{\xi}{X}\right)\operatorname{rect}\left(\frac{\eta}{Y}\right)\right\} = XY\operatorname{sinc}(Xf_X)\operatorname{sinc}(Yf_Y) $$

Deltas → Cosine: The shift theorem gives $\mathcal{F}\{\delta(\eta - a)\} = e^{-j2\pi a f_Y}$. Applying this to each delta and simplifying with Euler’s formula (Goodman Appendix A):

$$ \mathcal{F}\left\{\delta(\xi)\left(\delta\left(\eta - \frac{\Delta}{2}\right) + \delta\left(\eta + \frac{\Delta}{2}\right)\right)\right\} = e^{-j\pi\Delta f_Y} + e^{j\pi\Delta f_Y} = 2\cos(\pi\Delta f_Y) $$

Combining: By the convolution theorem, convolution in the space domain becomes multiplication in the frequency domain:

$$ \mathcal{F}\{t_A\} = XY \operatorname{sinc}(X f_X)\operatorname{sinc}(Y f_Y) \cdot 2\cos(\pi\Delta f_Y) $$

The intensity (proportional to $|\mathcal{F}\{t_A\}|^2$) is:

$$ I(x, y) = \frac{4X^2 Y^2}{(\lambda z)^2} \operatorname{sinc}^2(X f_X) \operatorname{sinc}^2(Y f_Y) \cos^2(\pi\Delta f_Y) $$

where $f_X = x/(\lambda z)$ and $f_Y = y/(\lambda z)$.

Factor	Origin	Role in the pattern
$\operatorname{sinc}^2(X f_X)$	Diffraction by slit width $X$	Envelope along $x$
$\operatorname{sinc}^2(Y f_Y)$	Diffraction by slit height $Y$	Envelope along $y$
$\cos^2(\pi\Delta f_Y)$	Interference between two slits	Fringes along $y$

---

Part (b): Cross-Sectional Profiles

Along the $x$ axis ($y = 0$, so $f_Y = 0$):

Along the $x$ axis we set $y = 0$. In Fraunhofer diffraction, the spatial frequencies are:

$$ f_X = \frac{x}{\lambda z}, \qquad f_Y = \frac{y}{\lambda z} $$

So if $y = 0$, then $f_Y = \frac{0}{\lambda z} = 0$. This means we’re only seeing the zero-frequency (DC) behavior in the $y$ direction — no $y$-dependent structure is visible.

Substituting $f_Y = 0$ into the full intensity expression, the $y$-dependent factors evaluate to: $\operatorname{sinc}(Y \cdot 0) = \operatorname{sinc}(0) = 1$ (by definition, $\operatorname{sinc}(x) = \sin(\pi x)/(\pi x)$ approaches 1 as $x \to 0$) and $\cos(\pi\Delta \cdot 0) = \cos(0) = 1$. These factors collapse to 1 and drop out, leaving only the $x$-dependent term:

$$ I(x, 0) \propto \operatorname{sinc}^2\left(\frac{X}{\lambda z} x\right) $$

This is a pure sinc-squared envelope — the standard single-slit diffraction pattern. The first zeros occur at $x = \pm\lambda z / X$.

With $X/(\lambda z) = 10 \text{ m}^{-1}$, zeros are at $x = \pm 0.1 \text{ m}$.

Along the $y$ axis ($x = 0$, so $f_X = 0$):

$$ I(0, y) \propto \operatorname{sinc}^2\left(\frac{Y}{\lambda z} y\right) \cos^2\left(\frac{\pi\Delta}{\lambda z} y\right) $$

This is the product of:

• A sinc-squared envelope (from the slit height $Y$) with zeros at $y = \pm\lambda z / Y$
• Cosine-squared fringes (from the slit separation $\Delta$) with zeros at $y = \pm\lambda z / (2\Delta), \pm 3\lambda z / (2\Delta), \ldots$

With $Y/(\lambda z) = 1 \text{ m}^{-1}$ and $\Delta/(\lambda z) = 3/2 \text{ m}^{-1}$:

• The sinc envelope has zeros at $y = \pm 1 \text{ m}$
• The cosine fringes have zeros at $y = \pm 1/3, \pm 1, \pm 5/3, \ldots$ m

---

Problem 3: N-Slit Grating Diffraction (Goodman 4-8/4-10)

This problem generalizes the double slit to an $N$-slit grating — the foundation of spectroscopy and many optical instruments.

---

Part (a): The Grating Aperture

Sketch the aperture described by the amplitude transmittance function

$$ t_A(\xi, \eta) = \left\{\left[\operatorname{rect}\left(\frac{\xi}{X}\right)\operatorname{rect}\left(\frac{\eta}{Y}\right)\right] * \left[\frac{1}{\Delta}\operatorname{comb}\left(\frac{\eta}{\Delta}\right)\delta(\xi)\right]\right\} \operatorname{rect}\left(\frac{\eta}{N\Delta}\right) $$

where $N$ is an odd integer and $\Delta > Y$.

Building the aperture step by step:

Step	Operation	Result
1. Start	$\operatorname{rect}(\xi/X)\operatorname{rect}(\eta/Y)$	A single $X \times Y$ rectangular slit
2. Replicate	Convolve with $\frac{1}{\Delta}\operatorname{comb}(\eta/\Delta)\delta(\xi)$	Infinite copies along $\eta$, spaced $\Delta$ apart
3. Truncate	Multiply by $\operatorname{rect}(\eta/(N\Delta))$	Keep only $N$ slits (centered at origin)

Step 1: The unit cell. A single rectangle, $X$ wide and $Y$ tall. This is one slit of the grating.

Step 2: Infinite replication. The Dirac comb $\operatorname{comb}(\eta/\Delta)$ is a train of delta functions spaced $\Delta$ apart along $\eta$. Convolving with this creates identical copies of the rectangle at every integer multiple of $\Delta$. The $\delta(\xi)$ ensures no replication in the $\xi$ direction. The $1/\Delta$ normalization is a convention that simplifies the Fourier transform.

Step 3: Finite truncation. Multiplying by $\operatorname{rect}(\eta/(N\Delta))$ — a window of width $N\Delta$ centered at the origin — selects exactly $N$ slits. Since $N$ is odd, there is one slit centered at $\eta = 0$ and $(N-1)/2$ slits on each side.

The condition $\Delta > Y$ ensures the slits don’t overlap — there is opaque space between adjacent slits.

Result: $N$ identical rectangular slits ($X$ wide, $Y$ tall), equally spaced with center-to-center separation $\Delta$.

---

Part (b): Fraunhofer Diffraction Pattern

Find an expression for the intensity distribution in the Fraunhofer diffraction pattern of that aperture, assuming illumination by a normally incident plane wave.

The aperture $t_A$ has the structure $\{\text{rect} * \text{comb}\} \cdot \text{rect}$— a convolution inside, then a multiplication. The Fourier transform of each piece is computed separately, then combined using two theorems:

• Convolution theorem: convolution in space → multiplication in frequency
• Multiplication theorem: multiplication in space → convolution in frequency

The Fourier transform of each factor:

Rectangle: As in Problem 2, the separable rect pair gives (Goodman Table 2.1, 2.1.4):

$$ \mathcal{F}\left\{\operatorname{rect}\left(\frac{\xi}{X}\right)\operatorname{rect}\left(\frac{\eta}{Y}\right)\right\} = XY\operatorname{sinc}(Xf_X)\operatorname{sinc}(Yf_Y) $$

Comb: The comb function is a train of deltas: $\operatorname{comb}(\eta/\Delta) = \sum_{n=-\infty}^{\infty} \delta(\eta/\Delta - n)$. Applying the delta scaling property $\delta(\eta/\Delta - n) = |\Delta|\,\delta(\eta - n\Delta)$ gives:

$$ \frac{1}{\Delta}\operatorname{comb}\left(\frac{\eta}{\Delta}\right) = \sum_{n=-\infty}^{\infty} \delta(\eta - n\Delta) $$

The shift theorem gives $\mathcal{F}\{\delta(\eta - a)\} = e^{-j2\pi a f_Y}$. Each delta is at position $a = n\Delta$, so $\mathcal{F}\{\delta(\eta - n\Delta)\} = e^{-j2\pi n\Delta f_Y}$. Since the Fourier transform is linear, the transform of the sum is the sum of the transforms:

$$ \mathcal{F}\left\{\sum_{n} \delta(\eta - n\Delta)\right\} = \sum_{n} e^{-j2\pi n\Delta f_Y} $$

By the Poisson summation formula (Goodman 2.1.6), this infinite sum of complex exponentials equals another train of deltas:

$$ \sum_{n} e^{-j2\pi n\Delta f_Y} = \frac{1}{\Delta}\sum_{m=-\infty}^{\infty} \delta\left(f_Y - \frac{m}{\Delta}\right) = \operatorname{comb}(\Delta f_Y) $$

The result is a comb in the frequency domain with deltas at $f_Y = m/\Delta$ for integer $m$ — the grating orders.

Truncating rect: $\mathcal{F}\{\operatorname{rect}(\eta/(N\Delta))\} = N\Delta\operatorname{sinc}(N\Delta f_Y)$.

By the convolution theorem, the inner convolution (rect $*$ comb) becomes a product in the frequency domain: $XY\operatorname{sinc}(Xf_X)\operatorname{sinc}(Yf_Y) \cdot \operatorname{comb}(\Delta f_Y)$. This is the sinc envelope sampled at $f_Y = m/\Delta$.

By the multiplication theorem, the outer multiplication by $\operatorname{rect}(\eta/(N\Delta))$ becomes a convolution in the frequency domain with $N\Delta\operatorname{sinc}(N\Delta f_Y)$. Convolving the comb of deltas with the sinc replaces each delta at $f_Y = m/\Delta$ with a copy of the sinc, producing the ratio $\sin(N\pi\Delta f_Y)/\sin(\pi\Delta f_Y)$ — the array factor.

The result is:

$$ \mathcal{F}\{t_A\} = XY\operatorname{sinc}(X f_X)\operatorname{sinc}(Y f_Y) \cdot \frac{\sin(N\pi\Delta f_Y)}{\sin(\pi\Delta f_Y)} $$

The last factor is the array factor for $N$ equally-spaced sources. It arises from summing the contributions of $N$ slits with progressive phase shifts.

The intensity is:

$$ I(x, y) = \frac{X^2 Y^2}{(\lambda z)^2} \operatorname{sinc}^2(X f_X)\operatorname{sinc}^2(Y f_Y) \left[\frac{\sin(N\pi\Delta f_Y)}{\sin(\pi\Delta f_Y)}\right]^2 $$

Factor	Origin	Effect
$\operatorname{sinc}^2(X f_X)$	Single-slit width $X$	Envelope along $x$
$\operatorname{sinc}^2(Y f_Y)$	Single-slit height $Y$	Envelope along $y$
$[\sin(N\pi\Delta f_Y)/\sin(\pi\Delta f_Y)]^2$	$N$-slit interference	Sharp principal maxima at $f_Y = m/\Delta$

---

Part (c): Suppressing Even-Order Diffraction

What relationship between $Y$ and $\Delta$ can be expected to minimize the strength of the even-order diffraction components while leaving the zero-order component approximately unchanged?

The principal maxima of the grating (array factor) occur at spatial frequencies $f_Y = m/\Delta$, where $m$ is the diffraction order. Each principal maximum has an intensity modulated by the single-slit $\operatorname{sinc}^2$ envelope.

The sinc-squared envelope has zeros at $f_Y = n/Y$ for nonzero integer $n$. To suppress even orders ($m = \pm 2, \pm 4, \ldots$), we want the sinc zeros to coincide with the even-order maxima:

$$ f_Y = \frac{2m}{\Delta} = \frac{n}{Y} \quad \Longrightarrow \quad Y = \frac{n\Delta}{2m} $$

The simplest choice ($n = 1, m = 1$) gives:

$$ \boxed{Y = \frac{\Delta}{2}} $$

When each slit is exactly half as tall as the slit spacing, the sinc envelope has its first zero at $f_Y = 2/\Delta$ — which is where the second-order maximum sits. By periodicity of the sinc zeros, all even orders ($m = \pm 2, \pm 4, \ldots$) fall on zeros of the envelope and are suppressed.

The zero-order component ($m = 0$, at $f_Y = 0$) sits at the peak of the sinc envelope where $\operatorname{sinc}(0) = 1$, so it is unaffected.

---

Problem 4: Square Aperture with Central Obscuration (Goodman 4-9/4-11)

Find an expression for the intensity distribution in the Fraunhofer diffraction pattern of a square aperture with a square central obscuration. Assume unit-amplitude, normally incident plane-wave illumination.

This problem combines two ideas: the Fraunhofer diffraction of a square aperture and the principle of superposition (linearity) for handling obscured apertures.

---

Expressing the Aperture

The aperture is a large square (outer side $\ell_o$) with a smaller square hole blocked out (inner side $\ell_i$). We express the transmittance as the difference of two squares:

$$ t_A(\xi, \eta) = \operatorname{rect}\left(\frac{\xi}{\ell_o}\right)\operatorname{rect}\left(\frac{\eta}{\ell_o}\right) - \operatorname{rect}\left(\frac{\xi}{\ell_i}\right)\operatorname{rect}\left(\frac{\eta}{\ell_i}\right) $$

Term	What it represents
First rect product	Full square aperture (side $\ell_o$), no obscuration
Second rect product	The central square region (side $\ell_i$) that is blocked
Difference	The actual aperture — open where the outer square is, except blocked where the inner square is

---

The Fourier Transform

By linearity of the Fourier transform:

$$ \mathcal{F}\{t_A\} = \ell_o^2 \operatorname{sinc}(\ell_o f_X)\operatorname{sinc}(\ell_o f_Y) - \ell_i^2 \operatorname{sinc}(\ell_i f_X)\operatorname{sinc}(\ell_i f_Y) $$

Each rect pair transforms to a 2D sinc (as established in Homework 1 Problem 2a), scaled by the area of the corresponding square.

---

The Intensity Distribution

The intensity is the squared magnitude of the field:

$$ I(x, y) = \frac{1}{(\lambda z)^2} \left|\ell_o^2 \operatorname{sinc}(\ell_o f_X)\operatorname{sinc}(\ell_o f_Y) - \ell_i^2 \operatorname{sinc}(\ell_i f_X)\operatorname{sinc}(\ell_i f_Y)\right|^2 $$

where $f_X = x/(\lambda z)$ and $f_Y = y/(\lambda z)$.

---

Homework 3 Summary

Problem	Key Concept
1(a)	Fresnel transfer function exponent is linear in $z$ → cascading works by adding distances
1(b)	Exact (non-paraxial) transfer function also linear in $z$ → same cascading property
2(a)	Double-slit pattern = sinc envelope × cosine fringes; convolution theorem decomposes the aperture
2(b)	$x$-axis: pure sinc²; $y$-axis: sinc² envelope modulating cos² fringes
3(a–b)	$N$-slit grating: sinc envelope × array factor; principal maxima sharpen as $N$ grows
3(c)	$Y = \Delta/2$ (50% duty cycle) suppresses even diffraction orders
4	Obscured square aperture = difference of two sinc patterns; linearity of Fourier transform

The central theme of Homework 3 is Fraunhofer diffraction — the far-field regime where diffraction patterns are Fourier transforms of aperture functions. The convolution theorem provides an elegant decomposition for multi-element apertures (double slits, gratings), while linearity handles obscured apertures. Problem 1 establishes the consistency of Fresnel propagation, confirming that the transfer function approach gives physically meaningful results regardless of how we partition the propagation distance.

---

Homework 4: Phase Gratings, Fresnel Zones, and Self-Imaging

Homework 4 moves beyond simple amplitude apertures to explore phase gratings, on-axis Fresnel diffraction, and the remarkable self-imaging phenomenon. The unifying thread is the interplay between spatial structure and propagation: how periodic or structured objects encode information into diffraction orders, and how free-space propagation can reconstruct perfect images at specific distances without any lens.

What Is a Phase Grating?

In Homeworks 2 and 3, the apertures were amplitude objects — they blocked or transmitted light, modifying only its amplitude. A phase grating is different: it transmits all the light but modifies its phase as a function of position. Physically, this corresponds to a transparent plate with spatially varying thickness or refractive index.

Property	Amplitude grating	Phase grating
Transmittance magnitude	Varies ($0$ to $1$)	Constant ($= 1$)
Transmittance phase	Constant	Varies with position
Physical realization	Opaque/transparent stripes	Variable-thickness glass
Energy absorption	Yes — blocks some light	No — all light transmitted
Diffraction orders	Present	Present (often stronger than amplitude case)

The amplitude transmittance of a sinusoidal phase grating is:

$$ t_A(\xi, \eta) = e^{j(m/2)\sin(2\pi f_0 \xi)} $$

where $m$ is the modulation depth (peak-to-peak phase variation) and $f_0$ is the spatial frequency of the grating. Since $|t_A| = 1$ everywhere, no light is absorbed — all the incident energy is redistributed among the diffraction orders.

---

Problem 1: Resolving Power of a Sinusoidal Phase Grating (Goodman 4-13)

Two discrete spectral lines of a source are said to be “just resolved” by a diffraction grating if the peak of the $q$th-order diffraction component due to source wavelength $\lambda_1$ falls exactly on the first zero of the $q$th-order diffraction component due to source wavelength $\lambda_2$. The resolving power of the grating is defined as the ratio of the mean wavelength $\lambda$ to the minimum resolvable wavelength difference $\Delta\lambda$. Show that the resolving power of the sinusoidal phase grating is

$$ \frac{\lambda}{\Delta\lambda} = q \ell f_0 = qM $$

where $q$ is the diffraction order used in the measurement, $\ell$ is the width of the square grating, and $M$ is the number of spatial periods of the grating contained in the aperture. What phenomenon limits the use of arbitrarily high diffraction orders?

This problem connects the physical structure of a grating (its size and spatial frequency) to its ability to distinguish closely spaced wavelengths — the resolving power. The result is a cornerstone of spectroscopy.

---

Background: The Jacobi-Anger Expansion

The key to analyzing phase gratings is expanding the complex exponential transmittance into a sum of diffraction orders. The Jacobi-Anger identity provides this decomposition:

$$ e^{j(m/2)\sin(2\pi f_0 \xi)} = \sum_{q=-\infty}^{\infty} J_q(m/2) \, e^{j 2\pi q f_0 \xi} $$

where $J_q$ is the Bessel function of the first kind of order $q$.

Component	Meaning
$J_q(m/2)$	Amplitude of the $q$th diffraction order (depends on modulation depth $m$)
$e^{j 2\pi q f_0 \xi}$	Plane wave tilted at angle corresponding to frequency $q f_0$
$q = 0$	Undiffracted (zero-order) beam
$q = \pm 1$	First-order diffracted beams
$q = \pm 2, \pm 3, \ldots$	Higher-order diffracted beams

---

The Grating Transmittance

The full grating has finite extent — a square aperture of side $\ell$:

$$ t_A(\xi, \eta) = \operatorname{rect}\left(\frac{\xi}{\ell}\right) \operatorname{rect}\left(\frac{\eta}{\ell}\right) e^{j(m/2)\sin(2\pi f_0 \xi)} $$

The finite aperture multiplies the infinite phase grating. In the Fourier domain, each diffraction order (a delta function for an infinite grating) gets convolved with the Fourier transform of the rect aperture — a sinc function. This is exactly the same mechanism that gives $N$-slit gratings their finite peak widths (Homework 3 Problem 3).

The Fraunhofer diffraction pattern in the observation plane at distance $z$ is the Fourier transform of $t_A$, evaluated at $f_X = x/(\lambda z)$, $f_Y = y/(\lambda z)$. Using the Jacobi-Anger expansion:

$$ U(x, y) \propto \sum_{q=-\infty}^{\infty} J_q(m/2) \, \ell^2 \operatorname{sinc}\left[\ell\left(\frac{x}{\lambda z} - q f_0\right)\right] \operatorname{sinc}\left(\frac{\ell y}{\lambda z}\right) $$

Each order $q$ produces a sinc peak centered at position $x_q = q \lambda z f_0$ in the observation plane. The sinc envelope has its first zero at a distance $\lambda z / \ell$ from the peak center.

Quantity	Expression	Physical meaning
Peak position of order $q$	$x_q = q \lambda z f_0$	Where the $q$th-order beam lands
First zero of order $q$	$x_q + \lambda z / \ell$	Width set by finite aperture
Peak separation (orders $q$ and $q+1$)	$\lambda z f_0$	Set by grating spatial frequency

---

Deriving the Resolving Power

Two wavelengths $\lambda_1$ and $\lambda_2$ (with $\lambda_1 > \lambda_2$) produce their $q$th-order peaks at different positions because the peak location $x_q = q \lambda z f_0$ depends on $\lambda$. The Rayleigh criterion for “just resolved” states that the peak of one wavelength falls on the first zero of the other.

Peak of order $q$ at $\lambda_1$:

$$ x_q^{(\lambda_1)} = q \lambda_1 z f_0 $$

First zero of order $q$ at $\lambda_2$:

$$ x_q^{(\lambda_2), \text{zero}} = q \lambda_2 z f_0 + \frac{\lambda_2 z}{\ell} $$

Setting these equal (Rayleigh condition):

$$ q \lambda_1 z f_0 = q \lambda_2 z f_0 + \frac{\lambda_2 z}{\ell} $$

Canceling $z$ and rearranging:

$$ q f_0 (\lambda_1 - \lambda_2) = \frac{\lambda_2}{\ell} $$

With $\Delta\lambda = \lambda_1 - \lambda_2$ and $\lambda \approx \lambda_2$ (the two wavelengths are nearly equal):

$$ q f_0 \Delta\lambda = \frac{\lambda}{\ell} $$$$ \boxed{\frac{\lambda}{\Delta\lambda} = q f_0 \ell = qM} $$

where $M = f_0 \ell$ is the number of grating periods within the aperture.

---

What Limits the Use of High Diffraction Orders?

The $q$th-order diffracted beam propagates at an angle $\theta_q$ given by the grating equation:

$$ \sin\theta_q = q \lambda f_0 $$

For the beam to be a propagating (rather than evanescent) wave, we need $|\sin\theta_q| \leq 1$:

$$ q \lambda f_0 \leq 1 \quad \Longrightarrow \quad q \leq \frac{1}{\lambda f_0} $$

When $q\lambda f_0 \geq 1$, the diffracted wave becomes evanescent — it decays exponentially with distance from the grating and carries no energy to the far field. This is the evanescent wave phenomenon, and it imposes an upper limit on the usable diffraction order:

$$ q_{\max} = \left\lfloor \frac{1}{\lambda f_0} \right\rfloor $$

---

Problem 2: On-Axis Fresnel Diffraction (Goodman 4-19)

Find the intensity distribution on the aperture axis in the Fresnel diffraction patterns of apertures with the following transmittance functions (assume normally incident, unit-amplitude, plane-wave illumination):

(a) $t_A(\xi, \eta) = \operatorname{circ}\left(\sqrt{\xi^2 + \eta^2}\right)$

(b) $t_A(\xi, \eta) = \begin{cases} 1 & a \leq \sqrt{\xi^2 + \eta^2} < b \\ 0 & \text{otherwise} \end{cases}$

where $a < 1$, $b < 1$, and $a < b$.

This problem computes the intensity along the optical axis behind circular and annular apertures. Unlike Fraunhofer diffraction (which gives the far-field pattern), Fresnel diffraction captures the near-field behavior — and the on-axis intensity reveals a striking oscillatory pattern related to Fresnel zones.

---

Background: On-Axis Fresnel Diffraction Integral

For an observation point on the optical axis at distance $z$, the coordinates in the observation plane are $(x, y) = (0, 0)$. The Fresnel diffraction integral simplifies to:

$$ U(0, 0, z) = \frac{e^{jkz}}{j\lambda z} \iint t_A(\xi, \eta) \, e^{j\frac{\pi}{\lambda z}(\xi^2 + \eta^2)} \, d\xi \, d\eta $$

Factor	Meaning
$e^{jkz}$	On-axis propagation phase
$1/(j\lambda z)$	Amplitude scaling from the Huygens-Fresnel principle
$e^{j\frac{\pi}{\lambda z}(\xi^2 + \eta^2)}$	Quadratic (Fresnel) phase — the curvature of a spherical wavelet from each aperture point

For circularly symmetric apertures, converting to polar coordinates $(\rho, \theta)$ and integrating over $\theta$ eliminates the angular variable, leaving a one-dimensional radial integral.

---

Part (a): Circular Aperture

The transmittance $t_A = \operatorname{circ}(\sqrt{\xi^2 + \eta^2})$ describes a circular aperture of radius 1.

Converting to polar coordinates and applying $s = \rho^2$:

$$ U(0, 0, z) = \frac{e^{jkz}}{j\lambda z} \cdot 2\pi \int_0^1 \rho \, e^{j\frac{\pi}{\lambda z}\rho^2} \, d\rho = \frac{\pi e^{jkz}}{j\lambda z} \int_0^1 e^{j\frac{\pi}{\lambda z}s} \, ds $$

Evaluating the integral:

$$ \int_0^1 e^{j\frac{\pi}{\lambda z}s} \, ds = \frac{\lambda z}{j\pi} \left(e^{j\frac{\pi}{\lambda z}} - 1\right) $$

Combining:

$$ U(0, 0, z) = e^{jkz}\left(1 - e^{j\frac{\pi}{\lambda z}}\right) $$

The on-axis intensity is:

$$ I(0, 0, z) = \left|1 - e^{j\frac{\pi}{\lambda z}}\right|^2 = 4\sin^2\left(\frac{\pi}{2\lambda z}\right) $$

Feature	Value
Maximum intensity	$I = 4$ (when $\sin^2 = 1$, i.e., $z = 1/(2n\lambda)$ for odd $n$)
Zero intensity	$I = 0$ (when $\sin^2 = 0$, i.e., $z = 1/(2n\lambda)$ for even $n$)
Behavior as $z \to \infty$	$I \to 0$ (the argument $\pi/(2\lambda z) \to 0$)
Oscillatory period	Intensity oscillates with $1/z$, not with $z$

---

Part (b): Annular Aperture

The annular aperture transmits light only in the ring $a \leq \rho < b$. Following the same approach with $s = \rho^2$:

$$ U(0, 0, z) = \frac{\pi e^{jkz}}{j\lambda z} \int_{a^2}^{b^2} e^{j\frac{\pi}{\lambda z}s} \, ds = e^{jkz}\left(e^{j\frac{\pi a^2}{\lambda z}} - e^{j\frac{\pi b^2}{\lambda z}}\right) $$

The on-axis intensity is:

$$ I(0, 0, z) = \left|e^{j\frac{\pi a^2}{\lambda z}} - e^{j\frac{\pi b^2}{\lambda z}}\right|^2 = 4\sin^2\left(\frac{\pi(b^2 - a^2)}{2\lambda z}\right) $$

This has the same functional form as Part (a), with $b^2 - a^2$ (the area of the annular opening) replacing $1$ (the area of the full circular aperture).

Comparison	Circular aperture	Annular aperture
Transmitting area	$\pi \cdot 1^2 = \pi$	$\pi(b^2 - a^2)$
On-axis intensity	$4\sin^2\!\left(\frac{\pi}{2\lambda z}\right)$	$4\sin^2\!\left(\frac{\pi(b^2 - a^2)}{2\lambda z}\right)$
Oscillation “frequency” in $1/z$	$\pi/(2\lambda)$	$\pi(b^2 - a^2)/(2\lambda)$
Peak intensity	4	4

---

Problem 3: The Talbot Effect — Self-Imaging of Periodic Objects (Goodman 4-20)

Consider a general one-dimensional periodic object with an amplitude transmittance having an arbitrary periodic profile. Neglect the size of any bounding aperture, ignore the evanescent wave phenomenon, and assume that paraxial conditions hold. Show that at certain distances behind this object, perfect images of the amplitude transmittance are found. At what distances do these “self-images” appear?

This problem establishes one of the most remarkable results in wave optics: a periodic object illuminated by a plane wave produces perfect copies of itself at regular intervals downstream, with no lens required. This is the Talbot effect, discovered experimentally by Henry Fox Talbot in 1836.

---

Background: Fresnel Propagation of Spatial Frequency Components

A periodic object with period $L$ can be decomposed into a Fourier series:

$$ t(x) = \sum_{n=-\infty}^{\infty} c_n \, e^{j 2\pi n x / L} $$

Each term $e^{j 2\pi n x / L}$ is a plane wave component at spatial frequency $f_n = n/L$. Under Fresnel propagation over distance $z$, each spatial frequency component acquires a phase factor from the transfer function:

$$ H(f_n; z) = e^{-j\pi\lambda z f_n^2} = e^{-j\pi\lambda z n^2/L^2} $$

(The on-axis phase $e^{jkz}$ is a common factor that doesn’t affect the intensity or the relative structure of the field.)

The propagated field is therefore:

$$ U(x, z) = e^{jkz} \sum_{n=-\infty}^{\infty} c_n \, e^{j 2\pi n x / L} \, e^{-j\pi\lambda z n^2/L^2} $$

Component	Role
$c_n$	Fourier coefficient of the original grating (fixed)
$e^{j 2\pi n x / L}$	Spatial structure of the $n$th harmonic (fixed)
$e^{-j\pi\lambda z n^2/L^2}$	Propagation phase acquired by the $n$th harmonic (depends on $z$)

---

The Self-Imaging Condition

For the propagated field to reproduce the original transmittance (up to the common on-axis phase), we need:

$$ e^{-j\pi\lambda z n^2/L^2} = 1 \quad \text{for all integers } n $$

This requires:

$$ \frac{\pi\lambda z n^2}{L^2} = 2\pi p_n \quad \text{for some integer } p_n $$

Since this must hold for all $n$ simultaneously, and $n^2$ takes all perfect-square values $(0, 1, 4, 9, 16, \ldots)$, the condition simplifies to requiring it for $n = 1$ (since if it holds for $n = 1$, it automatically holds for all $n$ because $n^2 \cdot (\text{integer})$ is still an integer):

$$ \frac{\lambda z}{L^2} = 2m \quad \text{for positive integer } m $$

Solving for $z$:

$$ \boxed{z_m = \frac{2mL^2}{\lambda}, \quad m = 1, 2, 3, \ldots} $$

The fundamental self-imaging distance is the Talbot distance:

$$ z_T = \frac{2L^2}{\lambda} $$

Self-images appear at every integer multiple of $z_T$.

---

Problem 4: Self-Imaging with Concentric-Circle Frequency Spectra (Goodman 4-21)

A certain two-dimensional non-periodic object has the property that all of the frequency components of its amplitude transmittance fall on circles in the frequency plane, the radii of the circles being given by

$$ \rho_m = \sqrt{2ma} \quad m = 0, 1, 2, 3, \ldots $$

where $a$ is a constant. Assume uniform plane-wave illumination, neglect the finite size of the object and the evanescent wave phenomenon, and assume that paraxial conditions hold. Show that perfect images of the object form at periodic distances behind the object. Find the locations of these images.

This problem generalizes the Talbot effect from Problem 3 to a non-periodic, two-dimensional object. The key is that self-imaging doesn’t require periodicity per se — it requires that the propagation phases of all frequency components simultaneously return to their original values.

---

Background: Fresnel Propagation in 2D

Under Fresnel propagation, a 2D spatial frequency component at radial frequency $\rho$ in the frequency plane acquires the phase:

$$ H(\rho; z) = e^{-j\pi\lambda z \rho^2} $$

For a general object, $\rho$ takes a continuum of values and no single distance $z$ can reset all phases simultaneously (except $z = 0$). But this object has a special structure: its frequency content is restricted to discrete circles at $\rho_m = \sqrt{2ma}$.

---

Computing the Propagation Phase

For frequency components on the $m$th circle ($\rho_m = \sqrt{2ma}$):

$$ e^{-j\pi\lambda z \rho_m^2} = e^{-j\pi\lambda z \cdot 2ma} = e^{-j 2\pi m a\lambda z} $$

The $\rho_m^2 = 2ma$ structure is the crucial feature: squaring the radii produces values linear in $m$, which means the propagation phase is also linear in $m$:

$$ \phi_m = -2\pi m a\lambda z $$

If $\rho_m$ were…	Then $\rho_m^2 \propto$ …	Phase dependence on $m$	Self-imaging?
$\sqrt{2ma}$ (this problem)	$m$ (linear)	Linear in $m$	Yes — at integer multiples
$m \cdot \text{const}$ (periodic grating)	$m^2$ (quadratic)	Quadratic in $m$	Yes — but needs $n^2$ condition
Arbitrary	No pattern	No common period	Generally no

---

The Self-Imaging Condition

For a perfect image, all propagation phases must be integer multiples of $2\pi$:

$$ 2\pi m a\lambda z = 2\pi m p \quad \text{for some integer } p \text{ and all } m $$

Since $m$ appears as a common factor, this simplifies to a single condition:

$$ a\lambda z = p \quad \text{for positive integer } p $$$$ \boxed{z_p = \frac{p}{a\lambda}, \quad p = 1, 2, 3, \ldots} $$

Self-images appear periodically with fundamental period:

$$ \Delta z = \frac{1}{a\lambda} $$

---

Homework 4 Summary

Problem	Key Concept
1	Resolving power of phase grating: $\lambda/\Delta\lambda = qM$; evanescent waves limit usable orders
2(a)	On-axis Fresnel intensity behind circular aperture: $I = 4\sin^2(\pi/(2\lambda z))$ — oscillates with $1/z$
2(b)	Annular aperture: same form with $b^2 - a^2$ replacing 1 — fewer exposed Fresnel zones stretch oscillation
3	Talbot effect: periodic objects self-image at $z_m = 2mL^2/\lambda$
4	Non-periodic objects with $\rho_m = \sqrt{2ma}$ frequency structure also self-image, at $z_p = p/(a\lambda)$

The central theme of Homework 4 is the interplay between spatial structure and propagation phase. Problem 1 shows how a grating’s finite size limits its ability to distinguish wavelengths — resolving power depends on the total number of illuminated periods, whether the grating modulates amplitude or phase. Problem 2 reveals the oscillatory on-axis behavior behind circular apertures, a direct consequence of Fresnel zone interference. Problems 3 and 4 demonstrate self-imaging: when an object’s spatial frequency spectrum has the right discrete structure, free-space propagation periodically restores all phase relationships, producing perfect images without any lens.

---

Key Concepts Index

These mathematical tools form the foundation for analyzing optical systems, diffraction, and imaging:

• Dirac delta function — Impulse functions and sifting property
• Fourier transform — Spatial ↔ frequency domain conversion
• Convolution theorem — Multiplication ↔ convolution duality
• Fourier-Bessel transform — For circularly symmetric functions
• Bessel functions — Natural basis for radial symmetry
• Radon transform — Line integrals (projections)
• Projection-slice theorem — Foundation of CT imaging
• Computed tomography — Medical imaging from projections
• Angular spectrum — Plane wave decomposition of diffracted fields
• Airy pattern — Diffraction pattern of circular aperture (jinc function)
• Babinet’s principle — Complementary screens produce related diffraction patterns
• Nyquist sampling — Sampling rate for bandlimited signals
• Optical transfer function — Frequency-domain characterization of imaging system performance (MTF + PTF)
• Sinc function — Fourier transform of rectangle; ideal interpolation kernel
• Rayleigh criterion — Resolution limit for optical systems
• Scaling property — Narrow in space ↔ wide in frequency
• Fresnel diffraction — Near-field diffraction with quadratic phase transfer function
• Fraunhofer diffraction — Far-field diffraction as Fourier transform of aperture
• Transfer function — Frequency-domain characterization of propagation
• Diffraction grating — $N$-slit array; array factor produces sharp principal maxima
• Helmholtz equation — Time-independent wave equation governing monochromatic fields
• Talbot effect — Periodic objects self-image at multiples of $z_T = 2L^2/\lambda$
• Fresnel zones — Annular regions of alternating constructive/destructive interference
• Evanescent wave — Non-propagating field beyond the critical angle; limits grating resolution
• Grating equation — $\sin\theta_q = q\lambda f_0$; relates diffraction order to propagation angle
• Jacobi-Anger expansion — Decomposes sinusoidal phase modulation into Bessel-weighted harmonics