.. image:: https://mybinder.org/badge_logo.svg
   :target: https://mybinder.org/v2/gh/ehpor/hcipy/HEAD?labpath=doc/tutorial_notebooks/FourierTransforms/FourierTransforms.ipynb
   :alt: Open in Binder

Working with Fourier Transforms
===============================

This tutorial explores Fourier transforms in HCIPy, from basic usage to
advanced techniques. Fourier transforms are fundamental to optical
simulations—they let us propagate light between pupil and focal planes,
analyze spatial frequencies, filter images, calculate PSFs and MTFs, and
perform image registration via shifts, rotations, and distortions.

We’ll start with the relationship between spatial and frequency domain
grids, then move through each transform type with concrete examples.

.. code:: ipython3

    import hcipy
    import numpy as np
    import matplotlib.pyplot as plt

The relationship between spatial and frequency domain grids is set by
the sampling theorem. In HCIPy, we create a spatial grid first—it
defines where we sample our optical fields—and then derive the
corresponding Fourier grid from it. The grid spacing and total extent
are linked: a finer spatial sampling gives a wider frequency range, and
a larger spatial extent gives finer frequency resolution.

.. code:: ipython3

    # Create a spatial grid
    N = 128
    spatial_grid = hcipy.make_pupil_grid(N, 2.0)
    
    print("Spatial grid properties:")
    print(f"  Dimensions: {spatial_grid.dims}")
    print(f"  Spacing: {spatial_grid.delta}")
    print(f"  Total extent: {spatial_grid.delta[0] * spatial_grid.dims[0]:.4f}")


.. parsed-literal::

    Spatial grid properties:
      Dimensions: (128, 128)
      Spacing: [0.015625 0.015625]
      Total extent: 2.0000


The FFT grid has a reciprocal relationship to the spatial grid. The
maximum frequency we can represent—the Nyquist frequency—comes from the
spatial sampling: :math:`k_{\text{max}} = \pi / \Delta x`. The frequency
resolution comes from the total spatial extent:
:math:`\Delta k = 2\pi / L`. This means a larger aperture (bigger
:math:`L`) gives finer frequency sampling in the Fourier domain, while
finer spatial sampling extends the maximum frequency.

.. code:: ipython3

    # Create the corresponding FFT grid
    fft_grid = hcipy.make_fft_grid(spatial_grid)
    
    print("\nFFT grid properties:")
    print(f"  Dimensions: {fft_grid.dims}")
    print(f"  Spacing: {fft_grid.delta}")
    print(f"  Zero point: {fft_grid.zero}")
    
    # Calculate the relationship
    dx = spatial_grid.delta[0]
    N = spatial_grid.dims[0]
    L = dx * N
    df = 2 * np.pi / L  # Frequency spacing
    
    print(f"\nTheoretical relationships:")
    print(f"  Spatial extent L = {L:.4f}")
    print(f"  Frequency spacing df = 2π/L = {df:.4f}")
    print(f"  Maximum frequency = π/dx = {np.pi/dx:.4f}")


.. parsed-literal::

    
    FFT grid properties:
      Dimensions: (128, 128)
      Spacing: [3.14159265 3.14159265]
      Zero point: [-201.06192983 -201.06192983]
    
    Theoretical relationships:
      Spatial extent L = 2.0000
      Frequency spacing df = 2π/L = 3.1416
      Maximum frequency = π/dx = 201.0619


We can see that the FFT grid dimensions match the spatial grid, and the
spacing follows the :math:`2\pi / L` relationship. The zero-point is at
the correct location for an FFT-centered grid. This reciprocal
relationship is the foundation for all the Fourier transforms that
follow.

Let’s visualize the spatial and frequency grids side by side to make the
reciprocal relationship concrete.

.. code:: ipython3

    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Plot spatial grid (1D slice)
    spatial_1d = spatial_grid.x.reshape(spatial_grid.dims)[N//2, :]
    axes[0].plot(spatial_1d, np.zeros(N), 'o', markersize=2, label='Grid points')
    axes[0].axvline(x=0, color='red', linestyle='--', alpha=0.5, label='Center')
    axes[0].set_xlabel('x')
    axes[0].set_ylabel('(zero)')
    axes[0].set_title('Spatial Domain Grid (1D slice)')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # Plot frequency grid points (1D slice)
    fft_1d = fft_grid.x.reshape(fft_grid.dims)[N//2, :]
    axes[1].plot(fft_1d, np.zeros(N), 'o', markersize=2, label='Grid points')
    axes[1].axvline(x=0, color='red', linestyle='--', alpha=0.5, label='Center')
    axes[1].set_xlabel('Spatial frequency (k)')
    axes[1].set_ylabel('(zero)')
    axes[1].set_title('Frequency Domain Grid (1D slice)')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()



.. image:: output_8_0.png


The spatial grid is centered at zero with 128 points spanning 2 units.
The frequency grid is also centered at zero but spans
:math:`[-\pi/\Delta x, \pi/\Delta x]`, also with 128 points. The density
of samples in each domain is inversely related—this is exactly what the
sampling theorem predicts.

We can control the Fourier domain sampling independently from the
spatial domain using two parameters: ``q`` (oversampling) and ``fov``
(field of view). The ``q`` parameter multiplies the number of samples in
the Fourier domain—higher ``q`` means finer frequency sampling at the
cost of a larger grid. The ``fov`` parameter scales the extent of the
Fourier domain—smaller ``fov`` crops the field of view, producing fewer
samples and faster computation, but risks aliasing if the field is too
small. Together they let us trade off resolution, field size, and speed.

Let’s demonstrate with a circular aperture.

.. code:: ipython3

    # Create an aperture
    aperture = hcipy.make_circular_aperture(0.5)(spatial_grid)
    
    # Different FFT grid configurations
    configs = [
        (1, 1, 'Critical sampling'),
        (2, 1, '2x oversampling'),
        (1, 0.5, 'Half FOV'),
        (2, 0.5, '2x oversampling, half FOV')
    ]

Now we’ll create FFT grids and propagators for each configuration and
visualize the resulting focal-plane field.

.. code:: ipython3

    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    axes = axes.flatten()
    
    for idx, (q, fov, title) in enumerate(configs):
        # Create FFT grid
        fft_grid_config = hcipy.make_fft_grid(spatial_grid, q=q, fov=fov)
    
        # Create FFT
        fft = hcipy.FastFourierTransform(spatial_grid, q=q, fov=fov)
    
        # Compute FFT
        wf = hcipy.Wavefront(aperture)
        focal_field = fft.forward(wf.electric_field)
    
        # Plot
        im = hcipy.imshow_field(np.log10(np.abs(focal_field)**2 / (np.abs(focal_field)**2).max()), vmin=-10, grid=fft_grid_config, ax=axes[idx])
        axes[idx].set_title(f'{title}\nq={q}, fov={fov}\nGrid: {fft_grid_config.dims}')
        plt.colorbar(im, ax=axes[idx])
    
    plt.tight_layout()
    plt.show()



.. image:: output_13_0.png


The effect is clear: at critical sampling (``q=1, fov=1``), the focal
plane just fits the Airy pattern. With ``q=2``, we get twice as many
samples across the same field—a smoother rendering of the PSF. Reducing
``fov`` to 0.5 halves the field shown, which is fine if we only care
about the core of the PSF and want faster computation. In physical
terms, ``q`` controls the sampling density in the focal plane while
``fov`` controls the total extent.

HCIPy provides several Fourier transform implementations.
``FastFourierTransform`` uses the FFT algorithm—O(N log N) and very
fast—but requires the output grid to be a regularly-spaced FFT grid
derived from the input grid. ``MatrixFourierTransform`` uses direct
matrix multiplication—O(N²) and slower—but accepts any arbitrary output
grid, giving complete flexibility. The choice is speed versus
generality.

Let’s compare them head-to-head.

.. code:: ipython3

    from hcipy.fourier import FastFourierTransform, MatrixFourierTransform
    import time
    
    # Create test field
    grid_256 = hcipy.make_pupil_grid(256)
    aperture_256 = hcipy.make_circular_aperture(1)(grid_256)
    wf_256 = hcipy.Wavefront(aperture_256)

First, the FFT on its native grid. We’ll time 10 forward transforms and
take the average.

.. code:: ipython3

    # Test 1: FFT (native grid)
    fft_grid_native = hcipy.make_fft_grid(grid_256)
    fft = hcipy.FastFourierTransform(grid_256)
    
    start = time.time()
    for _ in range(10):
        result_fft = fft.forward(wf_256.electric_field)
    fft_time = (time.time() - start) / 10

Now the MFT with a custom output grid—arbitrary size and spacing. The
flexibility comes at a computational cost.

.. code:: ipython3

    # Test 2: MFT (arbitrary output grid)
    custom_output_grid = hcipy.make_pupil_grid(300, 15.0)
    mft = hcipy.MatrixFourierTransform(grid_256, custom_output_grid)
    
    start = time.time()
    for _ in range(10):
        result_mft = mft.forward(wf_256.electric_field)
    mft_time = (time.time() - start) / 10

Let’s look at the timing comparison and the resulting fields side by
side.

.. code:: ipython3

    print(f"FFT time: {fft_time*1000:.2f} ms")
    print(f"MFT time: {mft_time*1000:.2f} ms")
    print(f"MFT/FFT time ratio: {mft_time/fft_time:.1f}x")
    
    # Visualize
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    hcipy.imshow_field(np.log10(np.abs(result_fft)**2 / (np.abs(result_fft)**2).max()),
                 vmin=-10, grid=fft_grid_native, ax=axes[0])
    axes[0].set_title(f'FFT (native grid)\n{fft_grid_native.dims} pixels')
    
    hcipy.imshow_field(np.log10(np.abs(result_mft)**2 / (np.abs(result_mft)**2).max()),
                 vmin=-10, grid=custom_output_grid, ax=axes[1])
    axes[1].set_title(f'MFT (custom grid)\n{custom_output_grid.dims} pixels')
    
    plt.tight_layout()
    plt.show()


.. parsed-literal::

    FFT time: 0.52 ms
    MFT time: 24.15 ms
    MFT/FFT time ratio: 46.4x



.. image:: output_22_1.png


The FFT is dramatically faster—often 10-100× depending on grid size. Use
it whenever you can work with the native FFT grid. The MFT shines when
you need a non-uniform or oddly-shaped output grid, such as zooming into
a small region of the focal plane without computing the full field, or
when working with irregular aperture shapes that don’t fit an
FFT-friendly grid.

Both ``FastFourierTransform`` and ``MatrixFourierTransform`` support
``.forward()`` for the transform from spatial to frequency domain and
``.backward()`` for the inverse transform. HCIPy handles the
normalization automatically—no need to multiply by pixel area or grid
spacing manually. The ``.backward()`` method undoes ``.forward()`` to
within numerical precision.

``FourierFilter`` lets us apply a filter function in the frequency
domain. We provide a function that operates on the internal FFT grid,
and ``FourierFilter`` handles the forward transform, filtering, and
inverse transform automatically. The ``q`` parameter controls
oversampling for the internal FFT grid, which helps avoid aliasing
artifacts when the filter has sharp features.

.. code:: ipython3

    from hcipy.fourier import FourierFilter
    
    # Low-pass filter
    def low_pass_filter(grid):
        k = np.sqrt(grid.x**2 + grid.y**2)
        cutoff = 20.0
        return hcipy.Field(np.exp(-(k / cutoff)**2), grid)

We’ll create a test field consisting of a smooth Gaussian signal plus
high-frequency sinusoidal noise.

.. code:: ipython3

    # Create a field with noise
    grid_128 = hcipy.make_pupil_grid(128)
    x = grid_128.x
    y = grid_128.y
    
    # Pattern with high-frequency noise
    signal = np.exp(-(x**2 + y**2) / 0.1)  # Gaussian
    noise = 0.1 * np.sin(2 * np.pi * 10 * x) * np.sin(2 * np.pi * 10 * y)
    noisy_field = hcipy.Field(signal + noise, grid_128)

.. code:: ipython3

    # Apply filter
    filt = hcipy.FourierFilter(grid_128, low_pass_filter, q=2)
    filtered_field = filt.forward(noisy_field)
    
    # Visualize
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    hcipy.imshow_field(noisy_field, ax=axes[0])
    axes[0].set_title('Noisy Input')
    
    fft_grid_internal = filt.internal_grid
    filter_field = low_pass_filter(fft_grid_internal)
    hcipy.imshow_field(filter_field, ax=axes[1])
    axes[1].set_title('Low-Pass Filter\n(Fourier Domain)')
    
    hcipy.imshow_field(filtered_field.real, ax=axes[2])
    axes[2].set_title('Filtered Output')
    
    plt.tight_layout()
    plt.show()



.. image:: output_29_0.png


The low-pass filter attenuates the high-frequency noise while preserving
the low-frequency Gaussian signal. In the Fourier domain panel, we can
see the filter profile—a Gaussian centered at zero frequency. After
filtering, the output is much smoother than the input, with most of the
10-cycle-period ripple removed. The residual noise comes from
frequencies below the cutoff.

Shifts, shears, and rotations in the spatial domain correspond to linear
phase ramps, shears, and rotations in the Fourier domain — these are
fundamental properties of the Fourier transform, with the shift theorem
covering translations. HCIPy provides dedicated classes that apply these
transformations efficiently via the frequency domain: ``FourierShift``,
``FourierShear``, and ``FourierRotation``. These are useful for image
registration, sub-pixel alignment, and distortion correction.

Let’s demonstrate with a simple test pattern.

.. code:: ipython3

    # Test pattern: letter F
    grid_2d = hcipy.make_pupil_grid(128, 2.0)
    
    def f_pattern(grid):
        pattern = grid.zeros()
    
        vertical_bar = hcipy.make_rectangular_aperture((0.15, 1.15), center=(-0.3, 0))
        top_bar = hcipy.make_rectangular_aperture((0.5, 0.15), center=(0, 0.5))
        middle_bar = hcipy.make_rectangular_aperture((0.5, 0.15), center=(0, 0))
    
        pattern[vertical_bar(grid) > 0.5] = 1
        pattern[top_bar(grid) > 0.5] = 1
        pattern[middle_bar(grid) > 0.5] = 1
    
        return pattern
    
    test_pattern = f_pattern(grid_2d)

Now let’s apply a shift, shear, and rotation to the test pattern.

.. code:: ipython3

    # Apply shift
    shifter = hcipy.FourierShift(grid_2d, [0.2, 0.3])
    shifted = shifter.forward(test_pattern)
    
    # Apply shear
    shearer = hcipy.FourierShear(grid_2d, shear=0.3, shear_dim=0)
    sheared = shearer.forward(test_pattern)
    
    # Apply rotation
    rotator = hcipy.FourierRotation(grid_2d, angle=np.pi/6)
    rotated = rotator.forward(test_pattern)

Each transform is created once and applied via ``.forward()``.

.. code:: ipython3

    # Visualize
    fig, axes = plt.subplots(2, 2, figsize=(12, 12))
    
    hcipy.imshow_field(test_pattern, ax=axes[0,0], cmap='gray')
    axes[0,0].set_title('Original F Pattern')
    
    hcipy.imshow_field(shifted.real, ax=axes[0,1], cmap='gray')
    axes[0,1].set_title('Shifted by (0.2, 0.3)')
    
    hcipy.imshow_field(sheared.real, ax=axes[1,0], cmap='gray')
    axes[1,0].set_title('Sheared (x-direction)')
    
    hcipy.imshow_field(rotated.real, ax=axes[1,1], cmap='gray')
    axes[1,1].set_title('Rotated by 30°')
    
    plt.tight_layout()
    plt.show()



.. image:: output_36_0.png


Each operation modifies the F pattern as expected: the shift moves it to
a new position, the shear skews it along one axis, and the rotation
turns it by 30 degrees. The slight artifacts at the edges come from the
periodic boundary conditions inherent in FFT-based transforms—padding
the input grid helps reduce these in practice.

As a practical example, let’s calculate the Modulation Transfer Function
(MTF) of a square-aperture optical system. The MTF describes how
contrast is transferred from the object to the image as a function of
spatial frequency—it’s one of the most important metrics for
characterizing an optical system’s resolution.

.. code:: ipython3

    # Create optical system
    pupil_grid = hcipy.make_pupil_grid(256, 2.0)
    
    # Square aperture
    aperture = hcipy.make_rectangular_aperture([1.0, 1.0])(pupil_grid)
    
    # FFT for PSF
    focal_grid = hcipy.make_fft_grid(pupil_grid, q=2)
    fft = hcipy.FastFourierTransform(pupil_grid, q=2)
    
    # Calculate PSF
    wf = hcipy.Wavefront(aperture)
    psf = np.abs(fft.forward(wf.electric_field))**2
    psf /= psf.max()

The MTF is the normalized magnitude of the Fourier transform of the PSF.
We’ll compute it by taking another Fourier transform of the PSF itself.

.. code:: ipython3

    # Calculate MTF
    mtf_grid = hcipy.make_fft_grid(focal_grid)
    fft_for_mtf = hcipy.FastFourierTransform(focal_grid)
    mtf_complex = fft_for_mtf.forward(psf)
    mtf = np.abs(mtf_complex) / np.abs(mtf_complex).max()

Now let’s visualize the pupil, PSF, and MTF side by side.

.. code:: ipython3

    # Plot results
    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
    
    # Pupil
    hcipy.imshow_field(aperture, ax=axes[0], cmap='gray')
    axes[0].set_title('Square Aperture (Pupil)')
    
    # PSF
    hcipy.imshow_field(np.log10(psf), vmin=-5, grid=focal_grid, ax=axes[1])
    axes[1].set_title('Point Spread Function (log scale)')
    
    # MTF
    hcipy.imshow_field(mtf, grid=mtf_grid, ax=axes[2])
    axes[2].set_title('Modulation Transfer Function')
    
    plt.tight_layout()
    plt.show()


.. parsed-literal::

    /tmp/ipykernel_3515/4071187771.py:9: RuntimeWarning: divide by zero encountered in log10
      hcipy.imshow_field(np.log10(psf), vmin=-5, grid=focal_grid, ax=axes[1])



.. image:: output_43_1.png


The MTF shows how contrast is transferred at different spatial
frequencies. For our square aperture, the MTF drops linearly from 1 at
zero frequency to 0 at the cutoff frequency—a triangular shape
characteristic of a sinc² PSF (the Fourier transform of a square). The
MTF tells us the finest detail the optical system can resolve and how
contrast degrades at smaller spatial scales.

When performing many transforms with different configurations, we can
use ``make_fourier_transform()`` to automatically select the best
algorithm based on the requested field of view and oversampling. This is
convenient for writing grid-agnostic simulation code that adapts to the
situation.

.. code:: ipython3

    from hcipy.fourier import make_fourier_transform
    import time
    
    # Large grid
    grid_large = hcipy.make_pupil_grid(1024)
    field = grid_large.zeros(dtype='complex')

With the full field of view (``fov=1``), the output grid is compatible
with the FFT, so ``make_fourier_transform`` selects the fast algorithm.

.. code:: ipython3

    # Full FOV: FFT
    fft_transform = hcipy.make_fourier_transform(grid_large, fov=1)
    print(f"Full FOV (fov=1): {type(fft_transform).__name__}")
    print(f"  Output grid: {fft_transform.output_grid.dims}")


.. parsed-literal::

    Full FOV (fov=1): FastFourierTransform
      Output grid: (1024, 1024)


With a small, highly-sampled region (``fov=1/32, q=4``), the output grid
no longer matches an FFT-friendly layout, so the MFT is automatically
selected.

.. code:: ipython3

    # Zoomed region: MFT
    mft_transform = hcipy.make_fourier_transform(grid_large, fov=1/32, q=4)
    print(f"\nZoomed FOV (fov=1/32, q=4): {type(mft_transform).__name__}")
    print(f"  Output grid: {mft_transform.output_grid.dims}")


.. parsed-literal::

    
    Zoomed FOV (fov=1/32, q=4): MatrixFourierTransform
      Output grid: (128, 128)


``make_fourier_transform`` automatically picks the right tool for the
job. When the full field of view is requested, the grid is
FFT-compatible and the fast algorithm wins. When zooming into a small
region with high oversampling, the MFT is selected because the
non-standard grid can’t use the FFT efficiently. This auto-selection
makes it easy to write grid-agnostic simulation pipelines.