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=201nufft2 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.
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.
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:
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
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
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”.
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):
This helped very much indeed. I wasn’t really aware that an NFFT is more like an inverse FFT rather than a FFT (that thought didn’t even cross my mind, and it didn’t seem terribly obvious from the stuff I read).
That makes this quote make much more sense to me as well
– I started digging into it, hopefully I can make some progress.
