nufft and aliasing errors

I am currently testing the nufft algortihm via the example:

1. function testnufft(N,r)
2. % N : number of nodes
3. % r : distorsion parameter 0 <= r <= 1
5. y = (0:N-1)'/N*2*pi; % equally spaced parametric absissas
6. t = y - r*sin(y); % uneven absissas
7. w = 1 - r*cos(y); % weights (dt/dy)
8. f = [0:N-1]'/2/pi; % equally spaced frequencies
9. %f = [0:fix(N/2)-1,-fix(N/2):-1]'/2/pi; % alternative frequencies
10. s = cos(t); % signal
11. S = nufft(s.*w,t,f)/N; % Fourier transform
13. subplot(211)
14. plot(abs(S))
15. subplot(212)
16. semilogy(abs(S))

For small r this seems to work, e.g., testnufft(16,1e-4);

If r increases, aliasing errors apper, e.g., testnufft(16,1e-2);

Incresing N does not help much, e.g. testnufft(128,1e-2);

Replacing the line: f = [0:N-1]'/2/pi; % equally spaced frequencies

by: f = [0:fix(N/2)-1,-fix(N/2):-1]'/2/pi; % alternative frequencies

helps a bit, e.g., now testnufft(128,1e-2); works.

But it does not work for more uneven data, e.g., testnufft(128,0.7); where aliasing errors appear.

Increasing the number of nodes does not help, e.g., testnufft(1024,1);

What did I miss in the use of nufft?

How to avoid this aliasing problem?

Answer 1

[Chris Turnes]

What you are observing here is the effects of non-uniform sampling; the further you get from uniform sampling, the less you can expect the evenly-spaced coefficients computed with nufft to resemble the "ground truth" you'd get from using uniform sampling.

Another way to consider this is to look at what happens to the transform itself if you try to "undo" it. That is, if the points are evenly spaced, you have the FFT, and if you reverse the transform you'll get back the original input (scaled by a factor of N). You can test this by explicitly building the DFT matrix in your code and seeing how close you get to an identity by doing the round trip transform:

N = 16;
y = (0:N-1)'/N*2*pi;    % equally spaced time points
f = (0:N-1)'/2/pi;      % equally spaced frequencies
Fideal = cospi(2*f(:)*y(:)') + 1j*sinpi(2*f(:)*y(:)');
norm(Fideal'*Fideal - N*eye(N))

The quantity Fideal' is a scaled inverse FFT, so it takes uniform frequency data back onto the uniform time grid; you expect the norm here to be small, because it means taking an FFT and then an inverse FFT gives you back something very close to the original data (there's some numerical round-off, so it won't be 0).

If you try to apply that inverse FFT from the uniformly-spaced frequency data that nufft gives you from your non-uniformly-sampled time grid, you'll find as you let the samples drift further and further from uniform you lose more and more of your ability to faithfully reconstruct the input because you're getting less and less of a faithful representation of the "true" spectrum:

r = 1e-8;
t = y - r*sin(y);       % uneven time points
Fact = cospi(2*f(:)*t(:)') + 1j*sinpi(2*f(:)*t(:)');
norm(Fideal'*Fact - N*eye(N))
% Compare this now as we increase r:
r = 1e-4;
t = y - r*sin(y);       % uneven time points
Fact = cospi(2*f(:)*t(:)') + 1j*sinpi(2*f(:)*t(:)');
norm(Fideal'*Fact - N*eye(N))
% Even more
r = 1e-2;
t = y - r*sin(y);       % uneven time points
Fact = cospi(2*f(:)*t(:)') + 1j*sinpi(2*f(:)*t(:)');
norm(Fideal'*Fact - N*eye(N))

It's just the mathematical nature of the transform and non-uniform sampling; when you stop sampling uniformly you start losing information about your data that is necessary for reconstruction. The more non-uniform you get, the less you can rely on the nufft results by themselves to give you a faithful estimate of the spectrum of the continuous-time signal you're sampling.

If what you're trying to do is analyze the spectrum of a non-uniformly sampled input, there are other tools out there that you might want to consider, such as plomb in the Signal Processing Toolbox.