I have implemented an N-dimensional Non-Uniform Discrete Fourier Transform (in this case it's specifically an inverse NUDFT) using PyTorch. My goal with this implementation is to have a function which is given an N-dimensional array of values and a list of non-grid-aligned, non-uniform points to interpolate its value at. For example, if I had the following pseudocode:

image = load_image("some_image.jpg")
frequency_space = fft_n(image)
target_points = random(32, 2)

output = nuidft(frequency_space, target_points)

I would expect output to contain plausible interpolated values for 32 randomly sampled, non-grid-aligned points in the input image.

Currently, I am doing this interpolation by evaluating the Fourier series of the input at the points, but my current implementation is egregiously time- and space-hungry, as it evaluates the IDFT equation for every target point on every one of the points in the frequency space representation of the signal, i.e.:

frequency_space = dft_n(x)

target_points = rand(batch_size, n_points, n_dims)

results = empty_like(*x.shape, *target_points.shape)
for coordinates, frequency in frequency_space, :
    for target_point in target_points:
        # evaluate IDFT complex sinusoids at target_point
        # store value in results array.

# take the n-dim mean (double sum for 2D input, triple sum for 3D input, etc)
output = results.mean(dims=range(n_dims))

This approach is also given as the first method for uniform to uniform IDFT in this person's blog post I found. (specifically the "Direct Evaluation of IDFT" section). Unfortunately, the later methods discussed in the same post are both specific to uniform synthesis (like scaling an image up or down uniformly by a constant factor), which is not what I'm trying to do here.

Now, I know that in the uniform IDFT, you can drastically reduce the time complexity by leveraging the transform's separability, i.e.:

idft_1d(idft_1d(x, dim=0), dim=1) == idft_2d(x)

So instead of doing $N^2$ operations for a square input of size $N$, you do $2N$ 1D DFTs, which can be cheaply computed using the FFT algorithm.

In my research, I've stumbled across the concept of an N-dimensional Non-Uniform Discrete Fourier Transform, which comes in three variants:

  • The nonuniform discrete Fourier transform of type I (NUDFT-I) uses uniform sample points ${\displaystyle p_{n}=n/N}$ but nonuniform (i.e. non-integer) frequencies ${\displaystyle f_{k}}$. This corresponds to evaluating a generalized Fourier series at equispaced points. It is also known as NDFT or forward NDFT
  • The nonuniform discrete Fourier transform of type II (NUDFT-II) uses uniform (i.e. integer) frequencies ${\displaystyle f_{k}=k}$ but nonuniform sample points ${\displaystyle p_{n}}$. This corresponds to evaluating a Fourier series at nonequispaced points. It is also known as adjoint NDFT.
  • The nonuniform discrete Fourier transform of type III (NUDFT-III) uses both nonuniform sample points ${\displaystyle p_{n}}$ and nonuniform frequencies ${\displaystyle f_{k}}$. This corresponds to evaluating a generalized Fourier series at nonequispaced points. It is also known as NNDFT.

But I am having trouble wrapping my head around which, if any, of these variants might be helpful to me. Thus, my question is: am I missing any obvious shortcuts? Am I doomed to evaluate every sinusoid once for every sample point? Or is there some separability magic I'm missing here than can reduce my complexity to something nicer, something with a $log$ in it instead of just a bunch of exponents?

Apologies if the question is vague, but if I've managed to communicate what I'm searching for adequately, I would massively appreciate being pointed in the right direction.


1 Answer 1


There are python modules available for the NUFFT, the NUFFT function can be used for this sort-of inverse operation by exchanging the frequencies and times, and multiplying the times by negative 1. Note that this requires the NUFFT moduThis can be shown from the definition of the NUDFT:
$X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i p_n f_k}, \quad 0 \leq k \leq N-1,$
you can see that exchanging p (the sample points) with f (the frequencies) and inverting one of them, we get
$X_k = \sum_{n=0}^{N-1} x_n e^{2\pi i f_n p_k}, \quad 0 \leq k \leq N-1,$
which is equivalent to the definition of the Uniform IDFT given below after swapping $X_k$ and $x_n$ in the uniform case.
$x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k \cdot e^{i2\pi \tfrac{k}{N} n}$
I am not familiar enough with image processing or any of these particular libraries to say how to extend this to the two-dimensional case, but I don't believe it should be terribly complicated.

Several NUFFT modules exist and I am not sure which one would be the best for your application.

References: https://en.wikipedia.org/wiki/Discrete_Fourier_transform


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.