Question: Can FINUFFT support spatial-domain matched field processing?

[xichaoqiang]

Hi and thanks for the great work on FINUFFT!

I'm working on a seismic imaging problem using matched field processing (MFP), and I'm wondering whether FINUFFT can be used or adapted for spatial-domain beamforming.

The imaging formula I’m using is:

$$ F(x, y) = \sum_j c_j \cdot \exp\left(i \cdot \omega \cdot \frac{\sqrt{(x - x_j)^2 + (y - y_j)^2}}{v}\right) $$

Where:

Where:

• $(x_j, y_j)$ are the receiver positions,
• $(c_j)$ is the complex signal at each sensor (at a given frequency),
• $(x, y)$ is a candidate source location,
• $\omega = 2\pi f$ is the angular frequency,
• and $v$ is the wave speed.

This phase term is nonlinear in $(x, y)$ because it depends on the Euclidean distance between each receiver and candidate grid point.

My question:

Since this is not a linear phase expression (i.e., not of the form $\exp(i k \cdot x )$, it cannot be directly written as a standard NUFFT transform.

Is there any way to approximate or accelerate this kind of spatial-domain imaging using FINUFFT (or a generalization of it)?

Any ideas or references would be greatly appreciated.

Thanks again for your excellent work!

[xichaoqiang]

Complexity Analysis

Let:

• $M_x$: number of grid points in x-direction
• $M_y$: number of grid points in y-direction
• $N$: number of receivers

Then the total number of evaluations is:

$\mathcal{O}(M_x \cdot M_y \cdot N)$

This is quite expensive for high-resolution grids.

[ahbarnett]

Hi xichaoqiang,
If your expression did not have the sqrt, then it would take the form of a paraxial (Fresnel) approximation, which can be massaged into a NUFFT form: see https://github.com/ahbarnett/fresnaq and associated paper of mine.
However, the sqrt prevents this. Instead your form is closest to a Helmholtz kernel interaction in the plane,
and I actually don't understand physically why your kernel is $e^{ikr}$ where r = distance, and k = omega/v.
If your interaction was 2D waves, it would be the Hankel function $H_0^{(1)}(kr)$, which you note has an extra decay factor $r^{-1/2}$, and the tool to apply such a kernel between many sources and targets is the FMM: https://github.com/flatironinstitute/fmm2d
If the max $kr$ is not too large ($<10^3$) then you might be able to apply some KIFMM to your strange kernel,
eg https://github.com/bempp/kifmm

I'm moving this to a Discussion.

Best, Alex

[xichaoqiang]

Thanks for your kindly reply.
In seismic research, it is common to transform data from the time-space domain ( (t, x) ) to the frequency-wavenumber or frequency-velocity domain ( (f, k) ) or ( (f, v) ). This transformation typically involves two main steps:

1. Temporal Fourier Transform
   The first step is to transform the time-domain signals at each receiver into the frequency domain using the Fast Fourier Transform (FFT), which is computationally efficient:

   $$ \tilde{u}(f, x) = \int u(t, x) e^{-i 2\pi f t} , dt $$

2. Slant Stack / Beamforming (Transformation to Wavenumber or Velocity Domain)
   To account for irregular spatial sampling, a method such as slant stacking is applied. This process maps the spatially distributed data into the frequency-wavenumber ( (f, k) ) or frequency-velocity ( (f, v) ) domain by summing along linear moveout paths:

$$ U(f, p) = \sum_{j=1}^{N} \tilde{u}(f, x_j) , e^{-i 2\pi f x_j \frac{1}{v}} $$

where $(v)$ is the velocity , $( x_j )$ is the receiver location, and $( N )$ is the number of receivers.
This step typically requires complex-valued exponential computations and summation over all sensors for each velocity (or wavenumber), resulting in high computational complexity—especially for dense velocity grids or wide frequency bands.

To accelerate this step, I recently adopted the FINUFFT library, specifically the finufft1d3 routine, which is optimized for type-3 NUFFT (nonuniform source and nonuniform target). This allowed me to significantly improve performance while maintaining accuracy.

[xichaoqiang]

To further localize seismic sources in physical space, the transformed frequency-domain data can be projected from the frequency-slowness domain ((f, p)) back into the spatial domain ((x, y)), resulting in a process known as Matched Field Processing (MFP).

This projection is based on the principle of phase alignment: at each trial location ((x, y)), we compute the expected phase delays to each receiver based on the distance and wave velocity. The match between the observed and modeled wavefields provides an energy (or coherence) measure for source localization.

Mathematically, the spatial-domain beamforming or MFP power at position ((x, y)) can be expressed as:

$$ P(x, y) = \left| \sum_{j=1}^N \tilde{u}(f, x_j, y_j) , e^{-i \omega \tau_j(x, y)} \right|^2 $$

where:

• $(\tilde{u}(f, x_j, y_j))$ is the frequency-domain data at receiver (j),
• $(\omega = 2\pi f)$ is the angular frequency,
• $(\tau_j(x, y) = \frac{d_j(x, y)}{v})$ is the travel time from trial location ((x, y)) to receiver (j),
• $(d_j(x, y))$ is the distance between ((x, y)) and receiver (j),
• and $(v)$ is the assumed propagation velocity.

[xichaoqiang]

The library you recommended, https://github.com/flatironinstitute/fmm2d, the Helmholtz FMM is very close to what I need. The hfmm2d_st_c_g function appears to be functioning properly. It likely just needs the Hankel function to be replaced with an exponential function.

[xichaoqiang]

https://github.com/flatironinstitute/fmm3d, is also close to what I need. The hfmm3d_t_c_g function appears to be functioning properly. The denominator term is not needed.

[ahbarnett]

Unfortunately you can't just change the kernel function in FMM codes... they are specific to the kernel, because it solves a PDE. (The only exception is KIFMM's, which might be your answer).