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.
