Help needed with FFT and NFFT

Hi, I think I could use some help understanding nonuniform discrete Fourier transforms (NFFT). Either I’m not using the tools correctly or I’m misunderstanding how they work in general.

I read that when my signal has nonuniform time steps I should use a NFFT (type II, I think) instead of a FFT. Since my results looked all kind of weird, I tried with a small example, shown below. This is similar to what I want to use it on: multiple periods of sine waves with uneven time steps. The goal is to detect additional harmonics (of course there are none in this example).

x = range(0, 360; length=200)
x´ = x[1:end-1]  # skip 360° so it's periodic
c = cosd.(x´)

using FFTW
rfft(c) .|> abs |> plot

using FastTransforms
nufft2(complex.(c), x´, eps()) .|> abs |> plot

using NFFT
k = x´ ./ x[end] ./ 2  # nfft wants k to be between -0.5 and 0.5
nfft(k, complex.(c)) .|> abs |> plot

I don’t understand what’s going on with nfft, while the results so nufft2 make somewhat sense to me. However, if I use length=201 nufft2 is just zeros.

When I add the line

x´ = x´ .+ randn(length(x´))

before running cosd to simulate the uneven time steps, rfft still looks great, nfft looks better-ish than before, while nufft2 is all over the place.
I could even use 10randn and the rfft results still look good – shouldn’t I just always use the normal FFT? Breaking the FFT isn’t that hard (quadratically spaced time steps will do it) so it would be good to know how to work with NFFTs nonetheless.

If helpful I can also write a bit more about my actual application.

1 Like

This maps k over [0, 0.5). Mapping it over [-0.5, 0.5) produces good results.

1 Like

Ah yes, that must have been the remnants of some experiment I tried.

So the results I’m seeing here look okay (I need to get better at reading FFT from complex values…). However, if I now go to length=201 the results look off again, which doesn’t happen with rfft. Is this expected?

In the NFFT.jl package’s reference paper linked, it is specified that the size N, of the NDFT in the equidistant sampling domain, is restricted to even numbers.

1 Like

That issue is resolved. We now also support odd N. I need to update the arXiv paper…

7 Likes

I did a quick comparison, with less points so it clearer. (length refers here to c and , not x).
fft nfft

And one where the are shifted byrandn(length(x´)).
fft nfft randn
(Here I calculated k = (x´ .- minimum(x´)) ./ (maximum(x´) - minimum(x´)) .- 0.5)

Are these the correct results?
To my eyes in all cases fft looks better. What would be the propper use case for NFFT then?

I think here you may be mapping onto [-0.5, 0.5].
As indicated, if you map over [-0.5, 0.5) you get better NFFT results, which overlay the FFT perfectly:

1 Like

Ah, thanks for spelling it out for me. How do you make that happen though, could you share your code?
Having k = -0.5:0.03447285262855501:0.4997127262280953 for example didn’t make it any better.

The code for the plot above is provided here below:

using Plots; gr()
n = 30
x = 360/n * (iseven(n) ? (-n÷2:(n÷2-1)) : (-(n-1)÷2:(n-1)÷2))
c = cosd.(x)

using FFTW
using Plots; gr()
fft(c) |> fftshift .|> abs |> plot

using NFFT
k = x / 360       # k over [-0.5, 0.5) if even, (-0.5, 0.5) if odd
nfft(k, complex.(c)) .|> abs |> x -> plot!(x; ls=:dash)

EDIT:
The symmetry of support x above has been improved for the n-odd case

2 Likes

That looks much better, thank you.

Is this the general rule to form k? If so, would it make sense to have nfft calculate this automatically? If not, how would I do this with my uneven timesteps? I don’t really understand why this is correct and the other k aren’t.

The documentation says
image
and in the examples uses

k = range(-0.4, stop=0.4, length=J)
k = rand(2, J) .- 0.5
k = Float32.(LinRange(-0.5,0.5,64))

in different places. From there I didn’t know any better than “it doesn’t really matter” which clearly isn’t true. If all of your answers are basic stuff I’m not aware of, where could I read up on this?

No, there is no general rule for forming k and it makes no sense to let the NFFT automatically chose k in an equidistant way.

The point here really is that the NFFT and the FFT are two different things. Only in case that the nodes are taken equidistantly, they are the same. Otherwise not.

In the examples you show we just use some exemplary set of nodes, which you need to adapt to the actual nodes you have in your application.

yes, shifting nodes in frequency space has a huge impact.

Well, one needs some deeper experience with Fourier mathematics in order to use the NFFT correctly. Since you aim at calculating the Fourier coefficients from a time signal you likely need to investigate the inverse NFFT, which is different from the adjoint NFFT. Probably this documentation Tools · NFFT
or this paper
https://downloads.hindawi.com/journals/ijbi/2007/024727.pdf
are helpful. In MRI the NFFT is regularly used, usually under the name “gridding”.

3 Likes

Just so there is no misunderstanding, I meant the process of scaling of the time steps I have so they are [-0.5, 0.5). It’s clear that I can’t “invent” an equidistant grid if I want to take advantage of NFFT.

Well, that is unfortunate. I took a look at the NFFT paper (https://arxiv.org/pdf/2208.00049.pdf) and got immediately lost, but I’ll take a look at the other resources, too. I understand the basic priciple behind the normal fourier transform and can use a FFT, but that’s about where it ends.
I was hoping for a solution where could give a function the time steps in addition to the signal and get a more correct result than using a normal FFT on the signal (or at least equally straight-forward). If that isn’t possible and it takes getting deep in the mathematical weeds I probably can’t justify spending the time, sadly.

Ok, this scaling is no magic. Within the FFT (Fast Fourier transform - Wikipedia) one has the term kn/N in the exponential. The k/N leads to sampling points being located between [0, 1). If we wrap this at 0.5 (open interval) we get the interval [-1/2, 1/2).

Not totally your question but this really helped me understand the NFFT, so figured I may as well share somewhere. With an FFT, you’re rewriting your function sampled at N grid points into a series sum of N sines and cosines. If you were to inverse FFT, you’re evaluating those sines and cosines on the original grid (getting exactly your function back). But with the NFFT you can evaluate that series at any arbitrary points. This gets you something like interpolation (since its basically assuming theres no power beyond the nyquist of the original grid):

using NFFT, Distributions, Random, PyPlot
N_coarse = 32
N_fine = 256
grid_coarse = range(-0.5, 0.5, length=N_coarse+1)[1:end-1]
grid_fine = sort(rand(Uniform(-0.5, 0.5), N_fine))

f = randn(N_coarse)
f_resampled = real.(nfft(grid_fine, nfft_adjoint(grid_coarse, N_coarse, f)) / N_coarse)

plot(grid_coarse, f, "o-", label="original")
plot(grid_fine, f_resampled, ".-", label="resampled")
legend()
ylim(-2, 2)

image