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.

[Paul]

Hi @Chris Turnes

A. "when you stop sampling uniformly you start losing information about your data that is necessary for reconstruction."

If the number of query points is the same as the number of sample points, then Fact is square and invertible, isn't it, i.e., reconstructiong should be possible. Why is the code comparing Fideal to Fact? (I realize I'm resurrecting an old question here ....)

But, if the number of query points and sample points are not the same, then Fact is not square and reconstruction would be ..... difficult (at best).

B. "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."

Given A and B, I'm unclear on the utility of the NUDFT for non-uniformly sampled sample points (I can see some utility with uniformly sampled sample points and non-uniform query points). Do you know of any examples that show the utility of the NUDFT with non-uniform sampling of a continuous-time signal?

[Chris Turnes]

If the number of query points is the same as the number of sample points, then Fact is square and invertible, isn't it, i.e., reconstructiong should be possible.

Square, definitely -- invertible, from a practical sense, no. Even if it is invertible, the inverse is not going to be the Fourier matrix, so it would require solving a linear system. In practice, though, the further from uniform the sampling grid is, the more ill-conditioned the matrix is going to be and the harder it's going to be to get a reasonable solution. You can see my comment in this other nufft-related question for an example.

Do you know of any examples that show the utility of the NUDFT with non-uniform sampling of a continuous-time signal?

Yes, there are many. The key thing to remember is that the NUDFT is not an inverse of a Fourier transform; that is, when you have non-uniform samples you can't take a single NUFFT/NUDFT to simply recover perfect uniform frequency samples of the underlying function. In other words, if Fis the DFT matrix, and A is the NUDFT matrix for some set of samples, we don't have the property , so you can't expect to start with non-uniform data, take a NUDFT, then take an inverse DFT and get a nice uniform representation.

The utility of the NUFFT is that it is a fast way of applying the matrix A. So for a given problem, if you have a set of non-uniformly spaced samples y, you can view it through the following lens: suppose the underlying function is bandlimited and satisfies the requirements you typically see of a Fourier transform. Then if we had uniform samples x of the signal, we could take its Fourier transform into the frequency domain and then a non-uniform Fourier transform out of the frequency domain back onto (non-uniform) samples in the spatial domain and get a decent interpolation. That is, . Typically applications that have such samples then view this as an inverse problem and are trying to solve the system in some least-squares sense, and typically with some kind of regularization or constraint. The ability to apply A very fast allows iterative solvers to solve these systems quite efficiently.

One area I have some limited experience with this in is 2-D/3-D MRI, where the samples acquired from the MRI machine are typically representative of Fourier-domain data and are measured along geometric contours that don't correspond to a uniform grid of samples. To reconstruct the spatial volume, this problem is sometimes phrased as I have above and least-squares solvers are used to generate the spatial volume that would produce those measurements. Without the NUFFT, such calculations would be prohibitively expensive.

[Paul]

Hi Chris,

Based on your response, I see now why your comment above was comparing Fact to Fideal. But, to be clear, I wasn't talking about going from non-uniform samples in the time domain back to uniform samples in the time domain. Rather, I was pointing out that if numel(f) == numel(t) then the F matrix is square and (always?) invertible, and therefore we can, in theory, recover the original non-uniform samples in the time domain from the NUDFT samples in the frequency domain.

rng(100);
N = 16;
y = (0:N-1)'/N*2*pi;    % equally spaced time points
f = (0:N-1)'/2/pi;      % equally spaced frequencies
% Even more
r = 1e-2;
t = y - r*sin(y);       % uneven time points
f = f + 1e-2*rand(size(f)); % uneven query points
%Fact = cospi(2*f(:)*t(:)') + 1j*sinpi(2*f(:)*t(:)');
Fact = cospi(2*f(:)*t(:)') - 1j*sinpi(2*f(:)*t(:)'); % changed sign to match nufft

cond(Fact)

ans = 1.1363

x = rand(N,1);
X = nufft(x,t,f);
norm(Fact\X - x)

ans = 2.4722e-15

In this case, I can (kind of, but not yet completely) see the NUDFT as some sort of frequency domain representation of the original signal, analagous to the DFT.

But, when numel(t) ~= numel(f) the original samples (x) are not reconstructible from the NUDFT (X). If I can't reconstruct x from X, what exactly does X mean?

[Vilnis Liepins]

Hi Paul,

But, when numel(t) ~= numel(f) the original samples (x) are not reconstructible from the NUDFT (X).

Not really, you can do the same as in Matlab fft:

X = nufft(x,t,f);

• If x is a vector and the length of x is less than length f, then pad x with trailing zeros and vector t with coresponding sample points (could be uniform) to length f.
• If x is a vector and the length of x is greater than length f, then truncate x and t to length f.

[Chris Turnes]

@Paul - Yes, to be clear, I think there are two separate thoughts here. The first, related to the original post, is that when you have non-uniform time samples and you take a non-uniform FFT from them, then you're getting a distorted sampling of the ideal spectrum that you don't have. That is, assuming there's some underlying continuous, bandlimited periodic function and you had sampled above the Nyquist rate at locations , the DFT would give you coefficients representing exact samples of the analog spectrum. If you have non-uniform locations t, if the locations are pretty close to a uniform grid, then the non-uniform DFT gives you interpolated samples that are somewhat close to being samples of the analog spectrum, but the more irregular the sampling gets, the more distorted those interpolated frequency spectrum values get from being samples of the ideal analog curve.

The second is that, yes, so long as your sampling gives you a well-conditioned enough problem, you could effectively solve a linear system represented by the NUDFT to recover your non-uniform samples. Through that lens, you could view the whole problem as, given uniform Fourier coefficients, use the NUDFT to interpolate a set of non-uniform time samples, which would look like solving the system where y represents the uniform Fourier coefficients and x represents non-uniform time samples.

I'll note that since you can always apply an inverse Fourier transform, this also exactly corresponds to the following problem: given uniform time samples, use the Fourier transform to interpolate onto non-uniform time samples, which then means you're solving the system . However, the more irregular your samples get, the worse and worse the conditioning of that system gets. So while it may be theoretically invertible, in practice things can start to "blow up" if the conditioning of the system is poor.

Of course, you can also view the opposite: given non-uniform samples x at locations t can we recover a uniformly-sampled version of the same underlying signal, which is then

I think the point most people get tripped up on with the non-uniform DFT is that if you're trying to use it to interpolate samples of a function between uniform and non-uniform grids, you pretty much always have to solve a linear system. Applying a NUFFT (+ FFT/IFFT) alone usually won't give you something that's a very good interpolation by itself.