How to use nufftn() for interpolation?

The documentation for nufftn() is very unclear, and neither of the Q/A's I found here give me enough info.

I put together a toy problem that illustrates what I'm attempting to do, that I think I should be able to do with nufftn(). The code below generates 100 scattered sample points x,y and 100 corresponding values of a smooth function v(x,y). I'd like to use a non-uniform discreet fft to convert those 100 scattered samples into a nice square 10x10 grid of interpolated values. I thought nufftn() would give me a 10x10 of frequency coefficients, but apparently not. The t doesn't come out the grid of gaussian values I expected. So, what must I do to make this function do what I was trying to do? Thanks!

EDIT : I rewrote my code as recommended in the Answer below, but something is still off:

x = rand(100,1);
y = rand(100,1);        % Random sampling pts in 2D
v = exp(-x.^2-y.^2);    % Values at sample points
t = {x,y};
[Xq, Yq] = meshgrid(linspace(0, 1, 10), linspace(0, 1, 10));
f = {Xq(:), Yq(:)};
c = nufftn(v,t,f);
p = reshape(c,[10,10]); % I thought is would be 10x10 freq values...
surf(Xq,Yq,p);
% If this code did what I'm trying, p would be approx exp(-(.05:.1:.95).^2-(.05:.1:.95)'.^2)

Here's the result:

Error using nufftn (line 107)

The total number of sample points must match the number of elements in the first input argument.

Error in ttt (line 8)

c = nufftn(v,t,f);

But if I change the curly brackets in Lines 5&7 to square:

t=[x,y];
[Xq, Yq] = meshgrid(linspace(0, 1, 10), linspace(0, 1, 10));
f = [Xq(:), Yq(:)];

...then I get big complex numbers for c and p. The abs(p) does have the Gaussian-ish shape I expected, but the values are complex and way off from the values of v.

I also tried oversampling, by changing the 100's in my first two lines to 10000, and it ran okay, but I get the same result of huge complex values in p with gaussian-ish magnitudes.

So, this is still just not quite doing what I wanted. Thanks in advance for help.

My old code in my original post:

x = rand(100,1);
y = rand(100,1);        % Random sampling pts in 2D
v = exp(-x.^2-y.^2);    % Values at sample points
c = nufftn(v,[x y]);    % ???
d = reshape(c,[10,10]); % I thought is would be 10x10 freq values...
t = ifft2(d);           % ...but apparently not.

Answer 1

[Jerry Guern]

Okay, so it turns out that

1) My method is non-ideal for this specific problem, but that's what I needed to find out.

2) nufftn() can totally do a successful FFT on random points, which was the goal of my Question

3) The nufftn() documentation available at the time I posted this doesn't really tell the story, and I had to get help from @Chris Turnes at Mathworks to get this to work.

The code below displays a smooth function and my attempt to estimate it from random samples using nufftn(). The range and overall shape of the estimate comes out right on, but it's very noisy, becaue the random samples aren't uniform and therefore aren't great for estimating Fourier Transform coefficients. And it takes major oversampling to make up for the randomness. But with some serious smoothing, this estimate would be not-too-bad, and this approach might be viable-ish for certain problems.

dim1 = 21; % Dimensions of output grid. Make these odd.
dim2 = 23;
nsample = 200000; % Number of random samples
oversample = nsample / (dim1 *dim2);

xs = rand(nsample,1) - 1/(2*dim1);  % Random sampling pts in 2D
ys = rand(nsample,1) - 1/(2*dim2);
v = exp(-xs.^2-ys.^2);    % Values at sample points

% Wave numbers for Fourier Transform
kx = [0 [1:floor(dim1/2)] [-floor(dim1/2) : -1]];
ky = [0 [1:floor(dim2/2)] [-floor(dim2/2) : -1]]';

pfd2 = reshape(nufftn(v, [ys(:) xs(:)], { ky kx }), dim2, dim1) / oversample;
estimate = ifft2(pfd2); % Our attempt at estimating v
[Xq, Yq] = meshgrid(linspace(0,1-1/(dim1),dim1),linspace(0,1-1/(dim2),dim2));

figure(50); % Display our estimate of v
surf(Xq,Yq,estimate); 
title_str = sprintf('ifft2 of nufftn()\n%d random data pts\noversampling = %.2f',nsample,oversample);
title(title_str);

figure(51);
surf(Xq,Yq,exp(-Xq.^2-Yq.^2));  % Display the actual v
title('The actual function, for comparison');

[Jerry Guern]

I haven't checked yet, but I would bet money that the lumpiness of this result would mirror the clumpiness of the randome sample points.

Answer 2

[Chris Turnes]

There are a few things that are probably worth explaining to answer your question. The key points are really:

• when using non-uniform transforms, you have to be careful about managing the scaling of the grid, and
• nufft is not an inverse of fft. The more non-uniform the sampling modalities get, the further and further you get from being able to get something close to the theoretical result of the signal from its original samples by using nufft and fft/ifft back to back.

That said, it might be helpful to fully explain things by beginning with a 1-D version. The kernel you're considering is the exponential function

This is a separable kernel, so we can simplify things by just reducing to 1 dimension:

The generalized m-point non-uniform DFT is given by

where the f(i) are whatever frequencies at which you choose to evaluate the transform. The default query points when you don't supply them are to take m = n and f = (0:(n-1))/n. So let's evaluate that:

% Set to a known random state.
rng('default');
n = 10;
x = rand(n,1);
v = exp(-x.^2);
V = nufft(v, x);
fdefault = (0:(n-1))/n;
Vref = exp(-1j*2*pi*fdefault(:)*x') * v;
norm(V-Vref)

ans = 3.3471e-14

So now let's go back to the original question. From the documentation for nufftn:

"When no query points are specified, the transform is computed at N evenly spaced query points in each dimension, where N = ceil(numel(X).^(1/D)) and D is the number of columns in t."

This means we expect a 10 x 10 grid that's an analog of what we had in 1-D, so let's see how it compares:

rng('default');
n = 10;
x = rand(n^2,1);
y = rand(n^2,1);        % Random sampling pts in 2D
v = exp(-x.^2-y.^2);    % Values at sample points
V = nufftn(v,[x y]);
fdefault = (0:(n-1))/n;
% Create the frequency grids for the row dimension and column dimension:
f1 = kron(ones(n,1), fdefault(:));
f2 = kron(fdefault(:), ones(n,1));
Vref = exp(-1j*2*pi*(f1*x(:)' + f2*y(:)'))*v;
norm(V-Vref)

ans = 1.2920e-13

So from this you should be able to tell what your frequency grid is -- it is the 10 x 10 grid with nodes (0:(n-1))/n in each dimension.

Now let's revisit your goal, and again to simplify things let's go to 1-D. What would happen if x was uniformly spaced?

% Set to a known random state.
rng('default');
n = 10;
x = 0:0.1:0.9;
v = exp(-x.^2);
V = nufft(v, x);
vr = ifft(V);

You'd think since we've uniformly sampled the function that this should give us back v right? Not exactly:

max(abs(vr - v), [], 'all')

ans = 3.8703

The issue here is grid scaling. nufft(g, x) is not the same thing as fft(g) -- if we wanted it to be equivalent, we'd have to scale x

norm(nufft(v,x*10) - fft(v))

ans = 0

or the query points

norm(nufft(v,x,0:9) - fft(v))

ans = 0

This is because of the grid scaling in the FFT. Typically, we think of the FFT as either taking sample points 0:(n-1) and query points (0:(n-1))/n or vice versa. You can consider it either way, or any other scaling such that the product of the ranges of values for the two equals just below n. So given that, let's consider the 1-D case again, and let's imagine you have ideal, uniformly-spaced samples:

n = 10;
x = (0:(n-1))'/n;
v = exp(-x.^2);
vr = ifft(nufft(v, n*x));
vt = exp(-((0:(n-1))/n).^2)';
norm(vr-vt)

ans = 4.3356e-16

As we move from uniform samples to non-uniform samples, you're going to start seeing the effects of what non-uniformity in sampling does. Let's move the samples just a tiny bit:

rng('default');
x = (0:(n-1))'/n + 1e-6*randn(n,1);
v = exp(-x.^2);
vr = ifft(nufft(v, n*x));
norm(vr-vt)

ans = 1.3277e-04

We've really not moved the samples by much, but already we're a good bit away from the "idealized" answer. Jitter the samples more and you'll see the effect take off:

rng('default');
x = (0:(n-1))'/n + 1e-3*randn(n,1);
v = exp(-x.^2);
vr = ifft(nufft(v, n*x));
norm(vr-vt)

ans = 0.1338

And if you get to the equivalent of the original example, then you've really lost most everything:

rng('default');
x = randn(n,1);
v = exp(-x.^2);
vr = ifft(nufft(v, n*x));
norm(vr-vt)

ans = 2.3695

This happens because the more non-uniform the samples get, the further and further the non-uniform DFT matrix gets from being the inverse of the (inverse) Fourier transform matrix.

So, in summary:

• Part of the reason the results are so different from what you expect is because you need to apply some grid scaling.
• Even with the grid scaling, though, the non-uniformity is severe enough that you're unlikely to see something that really looks like an exponential.

I'm not sure what your larger goal is, but typically folks use iterative methods to try to do this type of interpolation onto uniform samples with a degree of accuracy. You can also try interpolating onto a coarser grid and using "density compensation" for how your non-uniform samples are scattered.

Answer 3

[Umar]

Hi @Jerry Guern ,

After going through documentation provided at the link below,

https://www.mathworks.com/help/matlab/ref/double.nufftn.html

Let me break down the steps required and clarify the misunderstandings:

The nufftn( )function computes the NUDFT based on provided sample points and optionally, query points. Your goal is to convert your scattered samples into a structured grid, which requires careful setup of both the input data and the target grid. You have 100 random sample points (x, y) and corresponding values v. Make sure that x and y are properly formatted as matrices or vectors that correspond to the shape of your input data. The vector v should be a column vector of size 100, containing the values at each sample point. So, you need to specify the sample points correctly when calling nufftn( ). Instead of passing [x y], you should provide a cell array where each dimension's sample points are specified separately. For example:

     t = {x, y};  % Cell array containing x and y as separate entries

To get a grid output, you need to define query points that correspond to the desired grid structure. For a 10x10 grid, create evenly spaced points within the range of your original data.For instance:

     [Xq, Yq] = meshgrid(linspace(0, 1, 10), linspace(0, 1, 10));  % Create query 
     points
     f = {Xq(:), Yq(:)};  % Reshape them into vectors for nufftn

At this point, you can call nufftn() using both your sample points and your query points:

     c = nufftn(v, t, f);  % Compute NUDFT at query points

This will give you frequency coefficients at your specified grid locations. Since you will receive coefficients in an array format corresponding to your query points, you may not need to reshape it into a square grid directly unless required for further processing. If you want to visualize it as a grid or perform an inverse transform (like in your original code), ensure that you handle dimensions appropriately. The resulting variable c will contain your interpolated values at the specified grid points. If you want to visualize this output or compare it with expected Gaussian values:

     t = reshape(c, [10, 10]);  % Reshape if necessary for visualization
     surf(Xq, Yq, t);            % Visualize using surface plot

The expectation that the output resembles Gaussian values is valid; however, remember that interpolation will depend heavily on how well your scattered samples represent the underlying function across the defined domain.

Hope this helps.

If further issues arise or if additional clarifications are needed on specific aspects of this process, feel free to ask!

[Jerry Guern]

Thanks, I appreciate the response. I recoded it as you recommended, but now it just errors out. I edited my post to add my complete updated code. Any idea what's wrong?

[Jerry Guern]

I tried a couple of different changes and included those in my long edit. It's closer now, but something is still off.

[Umar]

Hi @Jerry Guern ,

I agree with @Chris Turnes comments. His feedback effectively underscores the critical importance of grid scaling in non-uniform transforms and highlights how small variations in sampling can lead to significant impacts on reconstruction accuracy. These insights are valuable for anyone working with Fourier transforms in non-uniform contexts, providing a clear framework for approaching potential pitfalls in signal processing tasks. If nufftn() does not yield the desired results, an alternative approach is to use the griddata() function, which is specifically designed for interpolating scattered data onto a grid. Here’s how you can implement this:

% Generate scattered sample points
numPoints = 100;
x = rand(numPoints, 1) * 10; % Random x-coordinates
y = rand(numPoints, 1) * 10; % Random y-coordinates

% Define a smooth function v(x,y)
v = exp(-((x-5).^2 + (y-5).^2)); % Gaussian function

% Create a grid for interpolation
gridSize = 10;
[xGrid, yGrid] = meshgrid(linspace(0, 10, gridSize), linspace(0, 10, gridSize));

% Interpolate using griddata
VGrid = griddata(x, y, v, xGrid, yGrid, 'cubic');

% Display the results
figure;
surf(xGrid, yGrid, VGrid);
title('Interpolated Values on 10x10 Grid using griddata');
xlabel('X-axis');
ylabel('Y-axis');
zlabel('Interpolated Values');

please see attached.

In nutshell, the nufftn() function can be used for nonuniform Fourier transforms, but it may not directly provide the interpolated grid you are looking for. The griddata() function serves as an effective alternative for interpolating scattered data onto a structured grid. By following the provided code examples, you should be able to achieve the desired interpolation of your scattered sample points into a 10x10 grid.

Hope this helps.

If you have further questions or need additional assistance, feel free to ask.

[Jerry Guern]

@Umar and @Chris Turnes I really appreciate the effort you've both put into writing all of that, but I need to understand how to use specifically nufftn(). That was the sole purpose of my toy problem, so using some other function to solve it defeats the purpose. I did read the nufftn() docs page too, but it's so vague and has no clear examples, so it was frustratingly unhelpful.

The problem with my result from nufftn() isn't just poor sampling, there's something I'm not getting about it even processes the inputs. If I crank the sampling waaaay up (see code below), I get the same weird results: a complex p with an abs() that's sort of the same shape as v but way off in magnitude, phase, and shape. I understand I'd get better results with griddata(), but I need to know how to use nufftn() to get the best results on this it can, and I don't think I'm even setting is up correctly to do that.

Thanks again for your help and patience.

x = rand(20000,1);       % We are now 200x oversampled
y = rand(20000,1);        % Random sampling pts in 2D
v = exp(-x.^2-y.^2);    % Values at sample points
t = [x,y];
[Xq, Yq] = meshgrid(linspace(0,1,10),linspace(0,1,10));
f = [Xq(:), Yq(:)];
c = nufftn(v,t,f);
p = fft2(reshape(c,[10,10]));
figure(1); surf(Xq,Yq,abs(p));   % If I was using nufftn() right, this would be near v,
figure(2); surf(Xq,Yq,real(p));  % ...this would also be near v,
figure(3); surf(Xq,Yq,imag(p));  % ...and this would be small.

[Umar]

Hi @Jerry Guern ,

You mentioned, “but I need to understand how to use specifically nufftn(). That was the sole purpose of my toy problem, so using some other function to solve it defeats the purpose. I did read the nufftn() docs page too, but it's so vague and has no clear examples, so it was frustratingly unhelpful. The problem with my result from nufftn() isn't just poor sampling, there's something I'm not getting about it even processes the inputs. If I crank the sampling waaaay up (see code below), I get the same weird results: a complex p with an abs() that's sort of the same shape as v but way off in magnitude, phase, and shape. I understand I'd get better results with griddata(), but I need to know how to use nufftn() to get the best results on this it can, and I don't think I'm even setting is up correctly to do that. “

Let me break down the function and your code to ensure you are using it correctly. First, try to understand the parameters of nufftn( ).

Input Array (X): This is the data you want to transform. It can be a vector, matrix, or multidimensional array.

Sample Points (t): These are the points at which your data is sampled. They can be provided as a vector, matrix, or a cell array of vectors, depending on the dimensionality of your input data.

Query Points (f): These are the points at which you want to evaluate the Fourier transform. Similar to sample points, they can be specified in various formats.

Now after taking a closer look at your provided code snippet:

x = rand(20000,1);       % Random x-coordinates
y = rand(20000,1);       % Random y-coordinates
v = exp(-x.^2 - y.^2);   % Values at sample points
t = [x, y];              % Sample points in 2D
[Xq, Yq] = meshgrid(linspace(0, 1, 10), linspace(0, 1, 10)); % 
Query grid
f = [Xq(:), Yq(:)];      % Query points in 2D
c = nufftn(v, t, f);     % Compute NUDFT
p = fft2(reshape(c, [10, 10])); % Reshape and compute 2D FFT

Sample Points (t): You are correctly defining t as a 2D array containing your random sample points. Ensure that the dimensions of t match the number of elements in v. In your case, both t and v have 20,000 elements, which is correct.

Query Points (f): The query points f are defined using meshgrid, which is appropriate for evaluating the Fourier transform over a grid. The reshaping of c into a 10x10 matrix for the fft2 function is also valid.

Output Interpretation: The output c from nufftn() represents the Fourier coefficients at the specified query points. The subsequent fft2 operation is applied to these coefficients, which may not yield results that directly resemble your original values v. This is because nufftn() computes the NUDFT, which is inherently different from a standard Fourier transform.

My recommendations for improvement would be

Check the Range of Sample Points: Ensure that the random values in x and y are within the range of your query points. If the sample points are too far from the query points, the results may not align well.

Increase Query Points: Instead of using a fixed grid of 10x10, consider increasing the resolution of your query points to better capture the variations in your data. For example, you could use linspace(0, 1, 100) instead of 10.

Visualize the Output: When visualizing the results, consider plotting the original values v against the transformed values to better understand the discrepancies. This can help identify if the issue lies in the sampling or the transformation process.

Here is an adjusted version of your code with higher resolution for query points:

x = rand(20000, 1);       % Random x-coordinates
y = rand(20000, 1);       % Random y-coordinates
v = exp(-x.^2 - y.^2);    % Values at sample points
t = [x, y];               % Sample points in 2D
[Xq, Yq] = meshgrid(linspace(0, 1, 100), linspace(0, 1, 100)); % 
Higher 
resolution query grid
f = [Xq(:), Yq(:)];       % Query points in 2D
c = nufftn(v, t, f);      % Compute NUDFT
p = fft2(reshape(c, [100, 100])); % Reshape and compute 2D FFT
figure(1); surf(Xq, Yq, abs(p));   % Visualize magnitude
figure(2); surf(Xq, Yq, real(p));  % Visualize real part
figure(3); surf(Xq, Yq, imag(p));  % Visualize imaginary part

Please see attached.

Make sure that your sample and query points are appropriately defined and by increasing the resolution of your query grid, you should be able to achieve better results with nufftn(). Remember, the nature of the NUDFT means that the output may not directly resemble the original sampled values, but with careful adjustments, you can obtain meaningful insights from your data. If you continue to experience issues, consider revisiting the documentation or exploring additional examples to deepen your understanding of the function's capabilities.

[Jerry Guern]

Thanks again for your time and help, @Umar At this point, I suspect the function just doesn't work as advertised. Your plots above show that the results is off by a factor of 10^7 and shaped nothing like v(x,y). That's not a mere discrepancy or sampling problem. @Chris Turnes I think this is a bug to report.

[Chris Turnes]

@Jerry Guern -- there isn't a bug in nufft / nufftn here -- they are computing what they advertise to compute, which is the summation formula for sample points t and query points f. You can verify this by calling the functions and comparing to a direct calculation of that formula -- the results agree to within the expected precision range. The difficulty here is in setting up the formulas to do what you want them to do.

Using nufft/nufftn to do interpolation is a very tricky problem, and it comes with a lot of pitfalls. It's not the inverse of a Fourier transform, so you can't expect it to give you back the formulaic samples you're looking for.

To dive into this a bit more, let's continue to focus on a 1-D case, because that's going to be easier to reason about and everything extends to 2-D and beyond fairly easily. Let's start with a far easier problem, which is to consider the opposite of what you're trying to do: suppose we have uniform samples of your function and we want to use nufft to interpolate onto a non-uniform grid. Before we even begin, I want to point out something that is very important -- doing interpolation involving the Fourier domain is necessarily going to involve a lot of assumptions about periodicity, because that's fundamental to Fourier domain representations. It's crucial to keep that in mind, because that is where a lot of mistakes are made.

That said, let's start with a set of uniform samples, and I'm going to pick an odd sample size because when you have FFTs of even lengths, some strange and hard-to-explain things happen at the Nyquist frequency. So let's take an example:

rng('default');
N = 11;
% Sort so that it is ordered in plots
x = sort(rand(N,1));
% Uniform samples on [0 1)
xu = (0:(N-1))'/N;
v = exp(-xu.^2);

In theory, to do interpolation, you're trying to:

• Transform these uniform samples into the Fourier domain.
• Take a non-uniform DFT to go back into the spatial domain at non-uniform locations.

So what happens if we try to directly do this?

% Transform to integer frequencies 0:(N-1), since the uniform time points
% are in [0 1)
V = fft(v);
% Minus sign for sample points because we're computing the adjoint
% transform, scale by N to account for Fourier matrix definitions.
cn = nufft(V, -(0:(N-1)), x) / N;
plot(x, exp(-x.^2), x, abs(cn));

This isn't a terrible interpolation -- it generally follows the shape of the theoretical thing we're chasing -- but one thing we can observe is that the results are complex, where we'd want them to be real:

cn.'

ans = 
  Columns 1 through 9

   0.9608 + 0.0433i   0.7816 + 0.0195i   0.8730 - 0.1085i   0.9353 + 0.0324i   0.7497 - 0.0013i   0.6523 + 0.0171i   0.5183 + 0.0308i   0.4677 + 0.0352i   0.4074 - 0.0548i

  Columns 10 through 11

   0.8265 - 0.5011i   0.9430 - 0.4496i

Let's consider that the FFT is giving us back a result with the DC frequency as the first element, and for real-valued inputs, the second half of the output is the flipped, conjugated copy of the first half:

[V(2:6) flip(conj(V(7:end)))]

ans = 
   0.0976 - 1.2275i   0.0976 - 1.2275i
   0.2582 - 0.5076i   0.2582 - 0.5076i
   0.2867 - 0.2782i   0.2867 - 0.2782i
   0.2958 - 0.1459i   0.2958 - 0.1459i
   0.2990 - 0.0459i   0.2990 - 0.0459i

When we are working with regular, uniform transforms, this symmetry is what gives us back real values if we invert the transform, because the exponentials for these frequency components are conjugate symmetric also:

[exp(1j*2*pi*xu) exp(1j*2*pi*xu*(N-1))]

ans = 
   1.0000 + 0.0000i   1.0000 + 0.0000i
   0.8413 + 0.5406i   0.8413 - 0.5406i
   0.4154 + 0.9096i   0.4154 - 0.9096i
  -0.1423 + 0.9898i  -0.1423 - 0.9898i
  -0.6549 + 0.7557i  -0.6549 - 0.7557i
  -0.9595 + 0.2817i  -0.9595 - 0.2817i
  -0.9595 - 0.2817i  -0.9595 + 0.2817i
  -0.6549 - 0.7557i  -0.6549 + 0.7557i
  -0.1423 - 0.9898i  -0.1423 + 0.9898i
   0.4154 - 0.9096i   0.4154 + 0.9096i
   0.8413 - 0.5406i   0.8413 + 0.5406i

But this doesn't happen with non-uniform spatial sampling points:

[exp(1j*2*pi*x) exp(1j*2*pi*x*(N-1))]

ans = 
   0.8180 + 0.5752i   0.9881 - 0.1539i
   0.6982 + 0.7159i  -0.1245 + 0.9922i
   0.5484 + 0.8362i  -0.8878 - 0.4603i
  -0.1781 + 0.9840i   0.2180 - 0.9759i
  -0.9569 - 0.2903i  -0.9809 + 0.1947i
  -0.6737 - 0.7390i  -0.4461 + 0.8950i
   0.3956 - 0.9184i   0.6017 + 0.7987i
   0.8299 - 0.5580i   0.9345 + 0.3559i
   0.8555 - 0.5178i   0.6671 + 0.7449i
   0.9646 - 0.2638i  -0.8908 - 0.4544i
   0.9758 - 0.2188i  -0.5934 - 0.8049i

This is where the periodicity becomes crucial. Since the DFT represents a periodic function, we can observe that the "integer" frequency (N-1) is the periodic repetition of the "integer" frequency -1:

[exp(1j*2*pi*xu*(N-1)) exp(1j*2*pi*xu*-1)]

ans = 
   1.0000 + 0.0000i   1.0000 + 0.0000i
   0.8413 - 0.5406i   0.8413 - 0.5406i
   0.4154 - 0.9096i   0.4154 - 0.9096i
  -0.1423 - 0.9898i  -0.1423 - 0.9898i
  -0.6549 - 0.7557i  -0.6549 - 0.7557i
  -0.9595 - 0.2817i  -0.9595 - 0.2817i
  -0.9595 + 0.2817i  -0.9595 + 0.2817i
  -0.6549 + 0.7557i  -0.6549 + 0.7557i
  -0.1423 + 0.9898i  -0.1423 + 0.9898i
   0.4154 + 0.9096i   0.4154 + 0.9096i
   0.8413 + 0.5406i   0.8413 + 0.5406i

For uniform samples spaced at 1/N, we can actually only represent integral frequencies up to N/2 (the Nyquist rate), and afterwards we are producing the periodic repetitions of the the "negative" frequency components -N/2 to -1. If we keep this in mind, then our "sample" points going into NUFFT are actually more accurately specified as [0:(N-1)/2, -(N-1)/2:-1]. And then an interesting thing happens with our non-uniform samples:

[exp(1j*2*pi*x) exp(-1j*2*pi*x)]

ans = 
   0.8180 + 0.5752i   0.8180 - 0.5752i
   0.6982 + 0.7159i   0.6982 - 0.7159i
   0.5484 + 0.8362i   0.5484 - 0.8362i
  -0.1781 + 0.9840i  -0.1781 - 0.9840i
  -0.9569 - 0.2903i  -0.9569 + 0.2903i
  -0.6737 - 0.7390i  -0.6737 + 0.7390i
   0.3956 - 0.9184i   0.3956 + 0.9184i
   0.8299 - 0.5580i   0.8299 + 0.5580i
   0.8555 - 0.5178i   0.8555 + 0.5178i
   0.9646 - 0.2638i   0.9646 + 0.2638i
   0.9758 - 0.2188i   0.9758 + 0.2188i

Now, the first and last columns of the non-uniform Fourier matrix are back to being conjugate symmetric, and we see:

cn = nufft(V, -[0:(N-1)/2 -(N-1)/2:-1], x) / N

cn = 
   0.9775 + 0.0000i
   0.9372 + 0.0000i
   0.9433 + 0.0000i
   0.9192 - 0.0000i
   0.7429 + 0.0000i
   0.6748 + 0.0000i
   0.5221 + 0.0000i
   0.4254 + 0.0000i
   0.4555 - 0.0000i
   0.7339 + 0.0000i
   0.7873 + 0.0000i

plot(x, exp(-x.^2), x, abs(cn));

This is a much better result, and it's (basically) real-valued. This pretty closely matches the exponential, except at the end where we get a jump up. That jump comes again from periodicity -- we've assumed by interpolating in the Fourier domain that we're dealing with a periodic function, and the last sample point is beyond (n-1)/n, so we're getting an interpolation between the end of the sampling range and the value of the function at the beginning of the range, which is 1. We're jumping back up because the periodic version of our uniformly-sampled function would have a large jump between the end of its range and the beginning, so this is interpolating that and smoothing that jump out.

The conclusion we can draw from this is that, yes, we can use interpolation with nufft to go from uniform samples to non-uniform samples, and things largely work as long as we take the right care.

So what about the other way? This is a much harder problem, because unlike the case where you start with uniform samples, you don't actually know the exact periodicity since you're starting with scattered samples. But we can use the easier problem to inform us about it. The DFT and non-uniform DFT are linear operators, so we can represent everything we've done here with matrices and analyze the problem with linear algebraic tools:

F = fft(eye(N));
A = nufft(eye(N), -[0:(N-1)/2 -(N-1)/2:-1], x) / N;
% Should be small, since this is the same operation just represented as
% matrices:
norm(A*F*v - cn)

ans = 1.0015e-15

In the easy "forward" problem -- uniform-to-non-uniform -- we are trying to compute when we are given x. If we consider the "backward" problem -- non-uniform-to-uniform -- we are effectively trying to solve that same linear system for x when we are given y. So let's try solving that system:

vn = exp(-x.^2);
% Use lsqminnorm to find the least-squares solution
vu_estimate = lsqminnorm(A*F, vn);
plot(xu, exp(-xu.^2), xu, abs(vu_estimate))

This isn't a great solution, but it's in part because the matrix is ill-conditioned. What if oversample though?

n = 10*N;
x = sort(rand(n,1));
vn = exp(-x.^2);
A = nufft(eye(N), -[0:(N-1)/2 -(N-1)/2:-1], x) / N;
% Use lsqminnorm to find the least-squares solution
vu_estimate = lsqminnorm(A*F, vn);
plot(xu, exp(-xu.^2), xu, abs(vu_estimate))

Now we have a very reasonable solution, especially given that we have implied periodicity, so we again get these edge effects.

This brings us to the most crucial element, though: unlike the DFT, the adjoint of the nufft is not anything close to the inverse of the system. Put another way, if we want to interpolate from non-uniform samples to uniform samples, we can't do it by computing (which is effectively what you're trying to do by computing a NUFFT and then taking an IFFT); we instead have to compute , where the dagger is the pseudo-inverse. And is not particularly close to unless the sample points are very close to uniform.

[Umar]

Hi @Jerry,

I understand your concerns about the significant deviation in results, which has led to your suspicion of a potential bug. However, as Chris pointed out, it appears that the nufft/nufftn functions are indeed performing as intended based on their underlying summation formula. The complexity of setting up these formulas correctly cannot be understated, especially when transitioning from uniform to non-uniform sampling. Chris's breakdown of the interpolation process highlights critical considerations regarding periodicity and its implications on our results. The examples provided effectively illustrate how careful attention to sampling points can significantly influence output accuracy.

Thank you once again @Chris Tunes for your valuable contributions to this discussion.

I will look forward to hearing @Jerry’s thoughts on this proposal.

[Jerry Guern]

@Chris Turnes and @Umar I appreciate that you've put a lot of time into this, but right now I am not a millimeter closer to understanding how to fix my code. You've convinced me I'm using nufftn incorrectly, please show me how to use it correctly. I beg you, no more paragraphs or 1-D examples. If you would just correct my little example code block so that it does a non-uniform DFFT on the samples and produces a grid that approximates the sample values, that would be helpful to me. Thank you both again for all the time you've put into this. I hope it ends up being productive time.

[Chris Turnes]

@Jerry Guern There's not a way to "correct" your sample code to get a better interpolation -- nufftn by itself is not a sufficient tool for this problem. Or rather, there's no call to nufftn that by itself will give you a good approximation of the uniform sample values. As I was aiming to demonstrate in my previous reply, you need to combine it with some kind of linear solver to get a more reasonable result.

This is my best stab at modifying your original code to give a more reasonable reconstruction. You can play around with the values of M (the number of non-uniform 2-D sample points) and Nx/Ny (the x- and y- interpolated grid sizes) and see how it affects the results. It is using lsqr and function handles to solve the inverse problem.

rng('default')
% --- Random sampling pts in 2D
M = 200;
t = rand(M, 2);     
% Number of uniform points we want in each dimension.
% Best to stick to an odd number, as things become very difficult with even
% numbers because the Nyquist frequency occupies a point in the spectrum.
Nx = 11;
Ny = 11;
% Shifted frequency points.  Combining this with fftshift/ifftshift to get
% better performance out of nufftn.
fx = (0:(Nx-1)) - floor(Nx/2);
fy = (0:(Ny-1)) - floor(Ny/2);
% Evaluate the function at the non-uniform points.
v = exp(-t(:,1).^2-t(:,2).^2);
% Set up the parameters for lsqr to solve the system.
% We'll use the "transformFcn" defined below for the forward/adjoint.
tol = 1e-12;
maxIt = 1000;
vu = lsqr(@(q, dir) transformFcn(q, dir, t, { fx fy }, [Nx Ny]), v, tol, maxIt);

lsqr converged at iteration 135 to a solution with relative residual 0.058.

vu = reshape(vu, [Nx Ny]);
% Check against reference answer:
tx = (0:(Nx-1))/Nx;
ty = (0:(Ny-1))/Ny;
surf(tx,ty,abs(vu));
hold on;
[txx, tyy] = ndgrid(tx, ty);
vu_actual = exp(-tx(:).^2 - ty(:)'.^2);
plot3(txx, tyy, vu_actual, '*r')
% For better view of the result:
caz = 96.7;
cav = 31.7;
view(caz,cav);

function b = transformFcn(q, direction, t, f, N)
% Forward direction = "notransp" = fft --> nufftn
% Adjoint direction = "transp" = nufftn --> fft
if direction == "notransp"
    b = forward(q, t, f, N);
else
    b = adjoint(q, t, f, N);
end
end

function b = forward(a, t, f, N)
% FFT from uniform samples (0:(N-1))/N --> 0:(N-1)
b = fftn(reshape(a, N));
% fft shift to get a centered domain:
b = fftshift(b);
% NUFFT to go from (-N/2):(N/2-1) --> t
b = nufftn(b, f, -t) / prod(N);
end

function b = adjoint(a, t, f, N)
% Adjoint NUFFT to go from t --> (-N/2):(N/2-1)
b = nufftn(a, t, f) / prod(N);
% ifft shift to get back to a centered domain
b = ifftshift(reshape(b, N));
% iFFT from uniform samples 0:(N-1) --> (0:(N-1))/N
b = ifftn(b) * prod(N);
b = b(:);
end

[Jerry Guern]

Chris, I posted a related Question that's more open-ended and probably within your wheelhouse. Would you mind taking a look? I'm trying to do an FFT on a set of sample points that are uniform but on a slanted paralellogram grid. I thought this was something nufftn() could do, but you convinced me it's not. Thanks in advance. https://www.mathworks.com/matlabcentral/answers/2167093-can-matlab-do-a-2-d-fft-on-a-uniform-parallelogram-grid-of-measured-points

[Chris Turnes]

@Jerry Guern -- I've answered on your other post, and indeed nufftn can help you compute what you want. Your other question is a bit different from the one you've posed here, and is well-suited to what nufftn can do.