How do I invert a non uniform fft using the nufft function?

I was trying to get the nufft to work, unfortunately, I don't understand how to invert it correctly. I have some scaling issues. I have only been working on symmetric domains.

To compute the inverse, I used the same reasoning of applying the fft twice (with a flip on the vector). And use the following expression,

% Sine
Us_nufft =   fac * nufft(Us, xnodes, k);
Us_t1 = J2 * fac * nufft( Us_nufft, -k, xnodes);

where fac is the scaling factor that I'm trying to find, xnodes is the space grid, k is the spectrum (of same length of xnodes), and J2 is another factor I tried using, but is set to 1 for this example.

Using the regular fft, I found that I had to scale the transform with (L / N), where L is the length of the domain in space, and N is the number of points. So I wanted to use that same scaling factor, working with the nufft on an equispaced grid. I wanted to have it work with equispaced grid at first, but I already encountered scaling problems.

By varying the number of points and the length of the domain, the scaling on the invert changed, so I thought I should still put a factor somewhere that should depend on N and L. I tried several different factors but they didn't work.

I also plotted the Fourier transforms and it appears that fac is correct, by comparing the Fourier transform of a Gaussian and the analytical version of the Fourier transform

as can be seen in the first plot. As factor I used the following code, and I don't know why using these J3 and J made it so that the curves aligned much better...

% scaling factor for equispaced
J = k(end) - k(1);
J3 = k(end-1) - k(1);
fac =  (L * J)/(N*J3);
J2 = 1;

I'll just mention here that the clear spikes in the cosine and the sine transform convinced me that the spectrum I chose was appropriate. To compare the inverse curves I plotted Us_t1, Us (the sine evaluated at xnodes) and Us_ifft against each other, where Us_ifft is obtained by computing the ifft on the transformed vector (on equispaced points it should not be problematic).

Us_ifft = fftshift(ifft(fftshift(Us_nufft)));

It appears that taking the ifft is completely off, I don't understand why, and the nufft version is close but a little off... And I don't know how to correct that.

In the plot, M is the maximum of the nufft-transformed vector. To obtain the two first plots, I used

N = 400;
dim = [-3*pi 3*pi];
dx = (dim(2) - dim(1)) / N;
xnodes = (dim(1):dx:dim(2)-dx);

When I wanted to use non-uniform grid points, things got even weirder. Most importantly, the scaling didn't seem uniform. Upon testing different mappings, it would seem that computing the inverse that way would amplify in denser regions and "muffle" in less dense regions. For this example, I chose the following nodes and dimensions :

N = 400;
dim = [-pi pi];
dx = (dim(2) - dim(1)) / N;

xnodes = linspace(dim(1), dim(2), N);
map1 =@(x) sign(x).* x.^2;
xnodes = map1(xnodes);

And the plot of the sine became :

The Fourier transform plots are also wrong :

I don't understand why it is the case, I don't know if I should multiply by the step-size vector or something (but doing that would use a vector of size N-1), and multiplying it pointwise with the nufft, wouldn't be possible. Should I consider increasing my vector using boundary conditions to do that ? Is it even the right thing to do ?

Thank you for your time !

Here's the full code :

% Testing of the nufft
%% Create space points
N = 400;
dim = [-pi pi];
dx = (dim(2) - dim(1)) / N;
%xnodes = (dim(1):dx:dim(2)-dx);
xnodes = linspace(dim(1), dim(2), N);
map1 =@(x) sign(x).* x.^2;

xnodes = map1(xnodes);
%dx_int = testch(2:end) - testch(1:end-1);
dx_int = (1 / N);

% The length of the domain. 
L = xnodes(end) - xnodes(1);
%% frequencies
% Fix a spectrum
k = linspace(-8, 8, N);


%% Attempt of scaling factor for the Fourier transform
% scaling factor for equispaced
J = k(end) - k(1);
J3 = k(end-1) - k(1);
fac =  (L * J)/(N*J3);
%J2 = J3;
% Nufft, cheby not mapping


% For nufft
%J2 = k(end-1) - k(1);
%J = k(end) - k(1);
%fac = (L) /(N);
%J2 = J3 / L;
J2 = 1;

%% Create the function vectors
% Sine function
Us = sin(xnodes);

% cosine
Uc = cos(xnodes);

% Guassian
gauss =@(x) exp(-x.^2);

% The analytical expression for the Fourier transform.
tr_gauss =@(x) sqrt(pi) * exp(-(pi * x).^2);
Ug = gauss(xnodes);
%% Compute the Fourier transforms using the nufft, and the inverses as well.
% Sine
Us_nufft =   fac * nufft(Us, xnodes, k);
Us_t1 = J2 * fac * nufft( Us_nufft, -k, xnodes);
Us_ifft = fftshift(ifft(fftshift(Us_nufft)));

% Cosine
Uc_nufft = fac * nufft(Uc, xnodes, k);
Uc_t1 = J2 * fac * nufft(Uc_nufft, -k, xnodes);
Uc_ifft = fftshift(ifft(fftshift(Uc_nufft)));

% Gaussian
Ug_nufft =  fac * nufft(Ug, xnodes, k);
Ug_t1 =  J2* fac *nufft(Ug_nufft, k, xnodes);
Ug_ifft = fftshift(ifft(fftshift(Ug_nufft)));

% Analytical fourier transform of Gaussian
fftg = tr_gauss(k);

%% Plot the results

%Sine

fig3 = figure('units', 'normalized', 'Outerposition', [0.3 0.2 0.4 0.6]);
ax2 = gca;
plot(xnodes, real(Us_t1))
hold on;
plot(xnodes, Us)
hold on;
plot(xnodes, real(Us_ifft))

% Add the max value to the plot
Ms = max(real(Us_t1)); 

% Labels and legend
ax2.YLabel.String = '"sin"';
ax2.XLabel.String = ['x (mapping nodes), M = ' num2str(Ms)];
title(['Inverse transform of the fft of the sine function on a modified space grid of ' num2str(N) ' points']);
legend(ax2, 'nufft', 'analytical', 'ifft', 'location', 'NorthEast');


% Cosine

fig32 = figure('units', 'normalized', 'Outerposition', [0.3 0.2 0.4 0.6]);
ax22 = gca;
plot(xnodes, real(Uc_t1))
hold on;
plot(xnodes, Uc)
hold on;
plot(xnodes, real(Uc_ifft))

% The max value
Mc = max(real(Uc_t1));

% labels and legend
ax22.YLabel.String = '"cos"';
ax22.XLabel.String = ['x (mapping nodes), M = ' num2str(Mc)];
title(['Inverse transform of the fft of the cosine function on a modified space grid of ' num2str(N) ' points']);
legend(ax22, 'nufft', 'analytical', 'ifft', 'location', 'NorthEast');


% Gaussian
fig4 = figure('units', 'normalized', 'Outerposition', [0.3 0.2 0.4 0.6]);
ax3 = gca;
plot(xnodes, real(Ug_t1))
hold on;
plot(xnodes, Ug)
hold on;
plot(xnodes, real(Ug_ifft))

% The max value
M = max(real(Ug_t1));

% Labels and legend
ax3.YLabel.String = '"gaussian"';
ax3.XLabel.String = ['x (mapping nodes), M = ' num2str(M)];
title(['Inverse transform of the fft of the gaussian function on a modified space grid of ' num2str(N) ' points']);
legend(ax3, 'nufft', 'analytical', 'ifft', 'location', 'NorthEast');


%% Plot on the fourier side
fig5 = figure('units', 'normalized', 'Outerposition', [0.15 0.05 0.7 0.9]);
subplot(3,1,1)
ax4 = gca;
plot(k, abs(Us_nufft))
subplot(3,1,2)
ax5 = gca;
plot(k, abs(Uc_nufft))
subplot(3,1,3)
ax6 = gca;
plot(k, abs(Ug_nufft))
hold on;

plot(k, fftg)
ax4.YLabel.String = '"fft(cos)"';
ax5.YLabel.String = '"fft(sin)"';
ax6.YLabel.String = '"fft(gaussian)"';
ax4.XLabel.String = ['k'];
ax5.XLabel.String = ['k'];
ax6.XLabel.String = ['k'];
sgtitle(['abs of the nufft of different functions on a modified space grid of ' num2str(N) ' points']);
legend(ax6, 'numerical', 'analytical', 'location', 'NorthEast');

Answer 1

[Michael Malmberg]

So this is a partial answer. I have not run your problem specifically. However, I was working on my own code in matlab 2021a and tried to use nufftn to do a 2D inverse transform.

I had taken an image and used nufftn to sample onto some radial coordinates, and then used nufftn again to transform it back to the original cartesian coordinates. It appears that all I had to do to make it work was apply a density compensation matrix to the frequency space data, in this manner (pseudo code). I work in spatial frequency, so if you are dealing with time based signals, then the analogy still follows, just in different words

yf = nufftn(original_image , spatial_coordinates , frequency_coordinates);

new_image = nufftn(yf .* densityCompensationMatrix, frequency_coordinates, -spatial_coordinates).

This worked for me. In actual implementation, I couldn't do a negative spatial coordinates directly without it taking super long (due to matlab's nufftn not dealing as effiiently with gridded coordinates going from positive to negative (see the isuniform function within nufftn)). So instead I just did:

new_image = fliplr(flipud(nufftn(yf .* densityCompensationMatrix, frequency_coordinates, spatial_coordinates)).

This worked for me because my spatial coordinates were symmetric as is. Hopefully this helps!

In general, density compensation just reflects the relative linear distance / area / volume (for 1D/2D/3D) in frequency space afforded to each sample, normalized to the linear distance/area/volume for if the samples were completely uniformly spaced. So it sounds like you were on the right track. I haven't analyzed mine completely for perfect quantitiative accuracy, but qualitatively so far it looks right.

[Chris Turnes]

In actual implementation, I couldn't do a negative spatial coordinates directly without it taking super long (due to matlab's nufftn not dealing as effiiently with gridded coordinates going from positive to negative (see the isuniform function within nufftn)).

This performance issue has been addressed in R2022b, as well as other improvements to try to detect uniformity that can be exploited. You should now observe the function's performance to be roughly equivalent whether or not the spatial coordinates are negated:

x = randn(100,1);
t = 0:99;
f = rand(100,1);
% warmup
for k = 1:5, nufft(x, t, f); end
tic; for k = 1:100, nufft(x, t, f); end, toc

Elapsed time is 0.050349 seconds.

tic; for k = 1:100, nufft(x, -t, f); end, toc

Elapsed time is 0.054934 seconds.

[Ajin Joy]

Can you give an example for how to choose the units or range of 'spatial/cartesian coordinates' and 'frequency/raidal coordinates' for the image?

[Michael Malmberg]

I guess it depends on your application. If you don't really have any real-world units you are comparing to, then you could probably just assign the range of the image coordinates arbitrarily. E.g. if you have a 256 by 256 image, you can make the "spatial coordinates" be defined with meshgrid as just -0.5:1/256:0.5-1/256. The choise of going rom -0.5 to 0.5 is arbitrary. What is important though is that you make the fourier coordinates consistent with the spatial coordinates, using their inverse relationships. I deal with MRI stuff, so this site helps things make sense for me:

k-space: FOV - Questions and Answers ​in MRI (mriquestions.com)

The radial coordinates I was choosing were just covering the same general area of frequency space with radial lines, with each line centered around 0. There are plenty of other nonlinear coordinates that you could choose though that I'm sure nufft would handle just fine.

[Michael Malmberg]

@Chris Turnes That is excellent news. Thanks for the update!

Answer 2

[James Burnet]

This is ultimately a linear algebra question.

The NUFFT is just an algorithm that approximates the NUDFT. If you construct the NUDFT matrix, F, you can see that F’*F ~= I (this would just be the DFT), which makes it obvious why you would not get the same answer.

What can be done though, is apply the QR decomposition to F.

[Q,R] = qr(F)
Q*R = F
Q’*Q = I
F*inv(R) = Q
% meaning 
(F*inv(R))’*(F*inv(R)) = I = inv(R)’*F’*F*inv(R)

Replacing F and F’ with the nufft like you already have shown, we can then see that

inv(R)’*nufft(nufft(inv(R)*x)) = x

This approach is not unique to nufft, but allows for orthonornalizing any transform. This could all be done with linear algebra instead obviously, without any need for nufft, but depending on the dimensions of the problem, it can be more time efficient.

It’s not necessary to generate the actual NUDFT matrices to accomplish this. The quick and dirty way would be to generate an m by m matrix A, use nufft on A to get an n by m matrix B, and multiply B*inv(A) to get an approximate for F.