How to use nufftn() for interpolation?
Accepted Answer
More Answers (2)
- 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.
- 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.
0 votes
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 entriesTo 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 nufftnAt 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 plotThe 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!
12 Comments
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.
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.
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.- Transform these uniform samples into the Fourier domain.
- Take a non-uniform DFT to go back into the spatial domain at non-uniform locations.
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:
(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.
Categories
Find more on Data Exploration in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
