> For the complete documentation index, see [llms.txt](https://zavalab.gitbook.io/statschemebook/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://zavalab.gitbook.io/statschemebook/chapter-6-statistical-data-analysis-and-learning/01-signal-processing/06-2d-fourier-transforms.md).

# 2D Fourier Transforms

* FTs for 2D signals are used to analyze and perform conv operations.
* Discrete FT of 2D field is:

  $$
  \hat{f}*{k\_1,k\_2}=\frac{1}{N}\sum*{j\_1=0}^{N-1}\sum\_{j\_2=0}^{N-1}f\_{j\_1,j\_2}e^{-i2\pi\left(\frac{j\_1k\_1}{N}+\frac{j\_2k\_2}{N}\right)},; k\_1,k\_2=0,...,N-1.
  $$
* 2D FT returns a 2D field in the frequency domain.
* Each freq comp contains a wave (sheet); original field is a superposition of sheets.
* As in 1D, can use convolution theorem to obtain $$f\*g=\mathcal{F}^{-1}{\mathcal{F}{f}\cdot \mathcal{F}{g}}$$.

## Example: 2D Signals as Superposition of 2D Waves

{% hint style="danger" icon="flask" %}

* 2D field (e.g., an image) is superposition of 2D fields (each with a different freq).
* 2D field with a given freq is sheet with wave pattern (higher freqs have finer waves).
* Consider 2D function:

  $$
  f(x\_1,x\_2)=\sin ; 2\pi \omega\_1(x\_1+x\_2)-\sin ; 2\pi \omega\_2(x\_1+x\_2)
  $$

  with frequencies $$\omega\_1=20$$ and $$\omega\_2=60$$.
* Below we plot this function in 3D and we can see that this forms a wavy field (sheet) and show two waves that form function.
* We also visualize function as 2D projection. This is what you would see if you observed wavy sheet from the top.
  {% endhint %}

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-57f42a5b1383d636cd24a46578518525c1b6f8e2%2Fch6_2dsinetotal.png?alt=media" alt="" width="489"><figcaption><p>2D function 𝑓(𝑥₁, 𝑥₂) shown in 3D and as a projection in 2D (top panels). Low-frequency component (middle panels) and high-frequency component (bottom panels).</p></figcaption></figure>

#### Python implementation

{% hint style="success" %}
Notebook [`ch6_fourier_transform_2D.ipynb`](https://colab.research.google.com/github/vzavala/StatsBook/blob/master/Python/ch6_fourier_transform_2D.ipynb) — open it in Colab and run the cells alongside this section.

MATLAB live script: [`ch6_fourier_transform_2D.mlx`](https://github.com/vzavala/StatsBook/blob/master/Matlab-Live/ch6_fourier_transform_2D.mlx)
{% endhint %}

A 2D signal as a superposition of 2D waves (Victor Z and Ugochukwu Ikegwu, UW-Madison, 2024): $$f(x\_1,x\_2)=\sin 2\pi\omega\_1(x\_1+x\_2)-\sin 2\pi\omega\_2(x\_1+x\_2)$$ with $$\omega\_1=20$$, $$\omega\_2=60$$ (`omega_1`, `omega_2`). The grid is fine enough (8 samples per period of the 60 Hz wave) to avoid aliasing.

**Note:** This notebook reproduces Example 6.20 (Section 6.3.5) of *Statistics for Chemical Engineers*: a 2D signal built from a low-frequency and a high-frequency wave.

```python
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 0.25, 125)
X1, X2 = np.meshgrid(x, x)
omega_1, omega_2 = 20, 60
f1 = np.sin(2*np.pi*omega_1*(X1 + X2))
f2 = np.sin(2*np.pi*omega_2*(X1 + X2))
f = f1 - f2

fig = plt.figure(figsize=(11, 13))
fig.subplots_adjust(left=0.03, right=0.97, bottom=0.04, top=0.97, wspace=0.15, hspace=0.3)  # fixed spacing (3D axes)
for row, (Z, lab) in enumerate(((f, 'f'), (f1, 'f_1'), (f2, 'f_2'))):
    ax = fig.add_subplot(3, 2, 2*row + 1, projection='3d')
    ax.plot_surface(X1, X2, Z, cmap='gray', linewidth=0)
    ax.view_init(elev=30., azim=10)
    ax.set_xlabel('$x_1$'); ax.set_ylabel('$x_2$'); ax.set_title(f'${lab}(x_1,x_2)$')
    ax.set_xticklabels([]); ax.set_yticklabels([]); ax.set_zticklabels([])
    ax2 = fig.add_subplot(3, 2, 2*row + 2)
    ax2.imshow(Z, cmap='gray', origin='lower', extent=[0, 0.25, 0, 0.25])
    ax2.set_xlabel('$x_1$'); ax2.set_ylabel('$x_2$')
plt.show()
```

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-3ba48a4d17735fa9f66ccc54a8d01eeb80a546d6%2Fch6_fourier_transform_2D-03-0.png?alt=media" alt="" width="520"><figcaption></figcaption></figure>

**Note:** As in the book, any image is a superposition of waves of different frequencies; identifying them lets us design filters.

***

## Example: 2D FT Applied to Filtering

{% hint style="danger" icon="flask" %}

* We analyze frequency content of a 2D uniform average filter (defined by a 2D step function) and of a 2D Gaussian filter (defined by a 2D Gaussian kernel).
  {% endhint %}

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-a084623b54708b3ac081fb595273fa94e303e8cf%2Fch6_2dfourierfilter.png?alt=media" alt="" width="489"><figcaption><p>Left panels: Uniform filter, Fourier transform in 2D, and Fourier transform in 3D. Right panels: Gaussian filter, Fourier transform in 2D, and Fourier transform in 3D.</p></figcaption></figure>

#### Python implementation

{% hint style="success" %}
Notebook [`ch6_fourier_transform_filters_2D.ipynb`](https://colab.research.google.com/github/vzavala/StatsBook/blob/master/Python/ch6_fourier_transform_filters_2D.ipynb) — open it in Colab and run the cells alongside this section.

MATLAB live script: [`ch6_fourier_transform_filters_2D.mlx`](https://github.com/vzavala/StatsBook/blob/master/Matlab-Live/ch6_fourier_transform_filters_2D.mlx)
{% endhint %}

2D Fourier transforms of a uniform (box) filter and of a Gaussian filter (Victor Z and Ugochukwu Ikegwu, UW-Madison, 2024). We show each filter $$g$$ (`g`), the magnitude $$\log(1+|\hat{g}(\omega\_1,\omega\_2)|)$$ (centered with `fftshift`), and the same as a 3D surface.

**Note:** This notebook reproduces the filters of Example 6.21 (Section 6.3.5) of *Statistics for Chemical Engineers*: 2D uniform and Gaussian filters and their Fourier transforms.

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft2, fftshift
from scipy.ndimage import gaussian_filter

L = 50
g_box = np.zeros((L, L))
g_box[19:30, 19:30] = 1                     # 2D step (uniform) filter
g_gauss = np.zeros((L, L))
g_gauss[L // 2, L // 2] = 1
g_gauss = gaussian_filter(g_gauss, sigma=1)  # 2D Gaussian filter

fig = plt.figure(figsize=(11, 13))
fig.subplots_adjust(left=0.03, right=0.97, bottom=0.04, top=0.97, wspace=0.15, hspace=0.3)  # fixed spacing (3D axes)
for col, (g, title) in enumerate(((g_box, 'uniform filter'), (g_gauss, 'Gaussian filter'))):
    g_hat = fftshift(np.log(1 + np.abs(fft2(g))))
    ax = fig.add_subplot(3, 2, col + 1)
    ax.imshow(g[10:40, 10:40], cmap='gray_r')
    ax.set_title(title); ax.set_xlabel('$x_1$'); ax.set_ylabel('$x_2$')
    ax = fig.add_subplot(3, 2, col + 3)
    ax.imshow(g_hat, cmap='gray_r')
    ax.set_xlabel(r'$\omega_1$'); ax.set_ylabel(r'$\omega_2$')
    ax = fig.add_subplot(3, 2, col + 5, projection='3d')
    W1, W2 = np.meshgrid(np.arange(L), np.arange(L))
    ax.plot_surface(W1, W2, g_hat, cmap='gray_r', linewidth=0)
    ax.set_xlabel(r'$\omega_1$'); ax.set_ylabel(r'$\omega_2$'); ax.set_title(r'$|\hat{g}|$')
    ax.set_xticklabels([]); ax.set_yticklabels([]); ax.set_zticklabels([])
plt.show()
```

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-63e61ffb965915d99aa2af165617946a6c92ce69%2Fch6_fourier_transform_filters_2D-03-0.png?alt=media" alt="" width="520"><figcaption></figcaption></figure>

**Note:** As in the book, the Gaussian filter keeps only the low frequencies (the bright region in the middle of its spectrum), which is why it preserves only coarse features.

***

{% hint style="danger" icon="flask" %}

* We now illustrate how to use FTs to perform convolutions.
* We take original image $$f$$ and obtain its 2D FT $$\hat{f}$$.
* Construct Gaussian filter with $$\sigma=1$$ ($${g}\_{1}$$). We obtain its FT $$\hat{g}\_1=\mathcal{F}{g\_1}$$.
* Construct Gaussian filter with $$\sigma=10$$ ($${g}\_{2}$$). We obtain its FT $$\hat{g}\_2=\mathcal{F}{g\_2}$$.
* Value of $$\sigma$$ for first filter is smaller (filter is designed to only eliminate high freqs).
* Obtain filtered images in freq domain: $$\hat{f}\_1=\hat{f}\cdot \hat{g}\_1$$ and $$\hat{f}\_2=\hat{f}\cdot \hat{g}\_2$$.
* Obtain filtered images in orig domain via inverse FT: $$f\_1=\mathcal{F}^{-1}{\hat{f}\_1}$$, $$f\_2=\mathcal{F}^{-1}{\hat{f}\_2}$$.
  {% endhint %}

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-34c6123e1bde81004f539a15c64157fd2267fa60%2Fch6_2d_fourier_image.png?alt=media" alt="" width="489"><figcaption><p>Top panels: Original and filtered images. Bottom panels: Frequency spectrum of original and filtered images.</p></figcaption></figure>

#### Python implementation

{% hint style="success" %}
Notebook [`ch6_image_filters_2D.ipynb`](https://colab.research.google.com/github/vzavala/StatsBook/blob/master/Python/ch6_image_filters_2D.ipynb) — open it in Colab and run the cells alongside this section.

MATLAB live script: [`ch6_image_filters_2D.mlx`](https://github.com/vzavala/StatsBook/blob/master/Matlab-Live/ch6_image_filters_2D.mlx)
{% endhint %}

Filtering an image with the convolution theorem (Victor Z and Ugochukwu Ikegwu, UW-Madison, 2024).

| Slides                                           | Python             | Meaning                                   |
| ------------------------------------------------ | ------------------ | ----------------------------------------- |
| $$f$$, $$\hat{f}=\mathcal{F}{f}$$                | `f`, `f_hat`       | image and its 2D FT                       |
| $$g\_1$$, $$g\_2$$ ($$\sigma=1$$, $$\sigma=10$$) | `g1`, `g2`         | Gaussian filters (same size as the image) |
| $$\hat{f}\_i=\hat{f}\cdot\hat{g}\_i$$            | `f1_hat`, `f2_hat` | filtered images in the frequency domain   |
| $$f\_i=\mathcal{F}^{-1}{\hat{f}\_i}$$            | `f1`, `f2`         | filtered images                           |

**Note:** This notebook reproduces Examples 6.19 and 6.21 (Section 6.3.5) of *Statistics for Chemical Engineers*: Sobel and Gaussian filters applied to an image, directly and with the convolution theorem.

```python
import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft2, ifft2, fftshift, ifftshift
from scipy.ndimage import gaussian_filter

f = plt.imread(path+'Data/droplet1g.jpg').astype(float)
if f.ndim == 3:
    f = f[..., 0]

def gaussian_kernel(shape, sigma):
    """Gaussian filter g as an image: a centred delta blurred with sigma (unit sum)."""
    g = np.zeros(shape)
    g[shape[0] // 2, shape[1] // 2] = 1.0
    return gaussian_filter(g, sigma=sigma)

f_hat = fft2(f)
results = []
for sigma in (1, 10):
    g = gaussian_kernel(f.shape, sigma)
    g_hat = fft2(ifftshift(g))                 # kernel centred at the origin
    fi_hat = f_hat * g_hat                     # convolution theorem
    fi = np.real(ifft2(fi_hat))
    results.append((sigma, fi, fi_hat))

def spectrum(z_hat):
    return fftshift(np.log(1 + np.abs(z_hat)))

fig, axs = plt.subplots(2, 3, figsize=(13, 8))
axs[0, 0].imshow(f, cmap='gray'); axs[0, 0].set_title('$f$')
axs[1, 0].imshow(spectrum(f_hat), cmap='gray'); axs[1, 0].set_title(r'$\hat{f}$')
for k, (sigma, fi, fi_hat) in enumerate(results, start=1):
    axs[0, k].imshow(fi, cmap='gray'); axs[0, k].set_title(rf'$f_{k}$ ($\sigma={sigma}$)')
    axs[1, k].imshow(spectrum(fi_hat), cmap='gray'); axs[1, k].set_title(rf'$\hat{{f}}_{k}$')
for ax in axs.ravel():
    ax.axis('off')
plt.tight_layout()
plt.show()
```

<figure><img src="https://1753452229-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2FakKE106umUCSWFtTnl31%2Fuploads%2Fgit-blob-cdc8a49bee2c65ee11a9dc85eebc8f5bbd7eaedc%2Fch6_image_filters_2D-03-0.png?alt=media" alt="" width="520"><figcaption></figcaption></figure>

Check: the result equals a direct (spatial) Gaussian convolution.

```python
direct = gaussian_filter(f, sigma=1, mode='wrap')
print(f"max |f_1 (FFT) - f_1 (direct)| = {np.max(np.abs(results[0][1] - direct)):.2e}")
```

```
max |f_1 (FFT) - f_1 (direct)| = 2.84e-13
```

**Note:** As in the book, the Sobel filter reveals edges and the Gaussian filter smooths the image; a small part of the spectrum is enough to capture its coarse features.