Can MATLAB do a 2-D FFT on a uniform parallelogram grid of measured points?

I have some measurements in slant coordinates, and I need to project them into the ground plane. When I just project the slant grid into ground coords, it ends up a uniform parallelogram grid. In order to form an image, I need to Fourier Transform my parallelogram grid of measured values into ground freq domain, the Inv FT that into ground time domain. The code below does this successfully, using a toy function f(x,y) = cos(3x-4y) in place of measurements. But in this code, I calculated non-fast Four Transform coefficients one at a time.

EDIT: I want to emphasize, this code does what I need, calculating the varialbe 'pfd' in Figure 41, I just need to know if some MATLAB function can somehow do it with a Fast Fourier Transform.

dim1 = 29; % Resolution of image we're creating
dim2 = 31;

% This block just creates xs & ys, coords of a random parallelogram of
% measured sample points.
theta1 = rand * pi/4;
theta2 = theta1 * (.5 + rand);
M = [2*(1+rand)/3 * [cos(theta1) sin(theta1)]; 2*(1+rand)/3 * [-sin(theta2) cos(theta2)]];
[Xs, Ys] = meshgrid(linspace(-3,3,300),linspace(-3,3,300));
gpt = M * [Xs(:)';Ys(:)']; %These points sprawl out all around the image area.
% Find the measured points that fall within the image.
keep = gpt(1,:)>-1/(2*dim1) & gpt(1,:)<1-1/(2*dim1) & gpt(2,:)>-1/(2*dim2) & gpt(2,:)<1-1/(2*dim2);
xs = gpt(1,keep);
ys = gpt(2,keep);
figure(1); scatter(xs,ys,.5); title('measured points within image')

% Create image pixels
[Xq, Yq] = meshgrid(linspace(0,1-1/(dim1),dim1),linspace(0,1-1/(dim2),dim2));
%Ratio of measured points within image to image pixels.
oversample = sum(keep) / (dim1 * dim2)
% 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]]';
% For this toy problem, use a closed form function in place of measured 
% data.  Calclate this function first for our image points.
td = cos(3*Xq-4*Yq);
% Calculate the frequency domain, from our td above and for our toy 
% measured data.
jfd = zeros(size(td));
nufd = zeros(size(td));
for j = 1:dim1
    for k = 1:dim2
        % Calculate the non-fast FFT, just for comparison
        jfd(k,j) = sum(sum(td          .* exp(-1i*2*pi*(kx(j)*Xq+ky(k)*Yq))));
        % Calculate same terms, but from our measured samples.
        pfd(k,j) = sum(cos(3*xs-4*ys) .* exp(-1i*2*pi*(kx(j)*xs+ky(k)*ys)));
    end
end
pfd = pfd / oversample;

figure(39); surf(Xq,Yq,td); title('orig function');
figure(40); surf(Xq,Yq,ifft2(jfd)); title('ifft2 of non-fast FT for comparison');
figure(41); surf(Xq,Yq,real(ifft2(pfd))); title('real of ifft2 of pfd');

So, is there any MATLAB function that do the same thing, calculate that pfd, but using an FFT?

I thought at one time that nufftn() did this, but I honestly couldn’t make any sense of the vague documentation or the comments in the code.

[Paul]

Hi Jerry,

If I'm following correctly, it seems that xs and ys define 1956 pairs of non-gridded points (xs and ys are both 1 x 1956). Is that the intent of the code? If so, my reading of nufftn is the "time" variables don't have to be equally spaced in each dimension, but they still have to be gridded, i.e., they can't be scattered points. In other words, we can have

t1 = sort(rand(1,N));
t2 = sort(rand(1,M));

and the function to be evaluated would be something like

z = cos(t1) + sin(t2');

Same thing for the frequency variables, i.e., they define a non-uniformly spaced grid (in general).

I'd try running some code here, but Answers doesn't seem to executing code right now (at least not for me).

[Jerry Guern]

Thanks for the reply. xs and ys will always have the same length, but that length will be different in each run, because they're on a randomly oriented grid . But it is always a uniform grid, albeit a slanted parallelogram grid. If you run the code, Figure 1 will show you what I mean. I tried to use nufftn() for this, but I couldn't make any sense of the documentation or code comments.

Answer 1

[Chris Turnes]

Yes, nufftn can compute those two loops:

jfd2 = reshape(nufftn(td, { Yq(:,1), Xq(1,:) },    { ky kx }),   dim2, dim1);
pfd2 = reshape(nufftn(cos(3*xs-4*ys), [ys(:) xs(:)], { ky kx }), dim2, dim1) / oversample;

If you compute norm(jfd - jfd2) and norm(pfd - pfd2) you'll find they match to within around 1e-12, 1e-13 in value.

[Paul]

@Chris Turnes

If I'm understanding pfd2 correctly ....

The 1D vectors xs and ys define pairs of points that are not gridded in x-y space. The sample points to nufftn uses [ ] to indicated that the sample points are given in a list. The query points input uses {} to indicate to nufftn that the query points are gridded. Is that correct?

Is the reshape really needed for pfd2?

Do you think some additional examples at nufftn would be helpful?

[Chris Turnes]

@Paul - yes, your understanding is correct. When sample/query points are specified as an array, such as t = [ys(:) xs(:)], each row of the array t specifies a coordinate set for a single location. When they are specified as a cell array of vectors, then they are specifying grid vectors for each dimension separately (though the grid itself may be irregularly spaced). If you have gridded data the cell syntax is preferable, as nufftn might be able to better optimize some code paths in that case.

The reshape is necessary because nufftn always outputs a vector so that it can have consistent output among equivalent calls. For example,

x = cos(3*xs-4*ys);
t = [ys(:) xs(:)]
pfd2 = nufftn(x, t, { ky, kx }) 

will yield the same shape of result as

[Ky, Kx] = ndgrid(ky, kx);
pfd3 = nufftn(x, t, [Ky(:) Kx(:)]);

[Paul]

The doc page nufftn sates: "When f is specified as a cell array of D vectors, the length of each dimension of the output Y is equal to the length of the corresponding vector in the cell array."

I'd try a quick experiment here, but can't currently run code on Answers.

[Chris Turnes]

@Paul - thanks, let me look into that. The behavior and the documentation don't appear to be in agreement, so one or the other should be fixed.

[Jerry Guern]

@Chris Turnes Thanks for the help! I saw your Answers, and I'm testing your solution right now. Btw, if someone is going to edit the nufftn() documentation, it would be good if they added a lot to it, like a full working example that shows how to use the output. The current documentation just shows you get an output from nufftn() but doesn't tell you it's structure or contents. Why is the output 1xN for 2-D inputs? How specifically must it be reshape()'d? Does it need an fftshift()? I have these answers now (except the first one) thanks to your code example, but they'd be helpful to others in the documentation.

[Jerry Guern]

@Chris Turnes Also, nothing in the documentation suggested in any way that I needed to input the wave numbers ky & kx to nufftn().

[Chris Turnes]

@Jerry Guern I'll reach out to our doc writer and see if we can find ways to make the nufftn documentation more helpful. It is a real challenge, because the formula for the non-uniform DFT in N-D is a difficult to write in a very succinct way, especially when you have mixes of sample and query points that are gridded and scattered.

I'll try to use the same naming conventions as in your script, and hopefully this will help. For a 2-D problem, the non-uniform DFT summation has the formula

In this formula, h is the input array to nufftn, G is its output, each pair (x(n), y(n)) is a sample point that gives the coordinates of h(n), and each pair (kx(m), ky(m)) is a query point that gives the frequency coordinates of G(m). For nufftn to be able to compute this sum, you have to tell it what all the sample and query points are.

Since in your code you're summing the terms over all of the xs and ys values, you know that xs and ys form your sample points. Since in the exponential term xs multiplies kx and ys multiplies ky, you know that kx and ky form your query points. As a result, to use nufftn, you need to form the matrix [ys(:) xs(:)] to give your scattered sample points and the set of vectors { ky, kx } to give your gridded query points. You can use vectors for the query points because your frequencies form a separable grid, which is not the case for your sample points.

[Jerry Guern]

Hey, @Chris Turnes, I cannibalized your code for this problem, and I got it to work for my other problem with random samples! https://www.mathworks.com/matlabcentral/answers/2160255-how-to-use-nufftn-for-interpolation

[Paul]

@Chris Turnes

For me, the doc page for nufftn would benefit greatly from additional examples illustrating the different options [] and {} for the sample points and query points that should/could be used depending on how the sample point and query point data is arranged. Thanks for all of your responses on nufftn threads. Very helpful.

@Jerry Guern

Did nufftn improve execution speed (or any other aspects) of your code, other than not having to write the loop yourself?

[Jerry Guern]

I'm glad you thought to ask this. I made my term-by-term code a more efficient to test. With dim1=291, dim2=311, and oversample coming out to 3.91, my term-by-term code finished in 336 seconds, while nufftn() finished in 246.6 seconds. So it's faster but tbh I expected more of a speedup than that, given that I'm doing a pre-Gauss non-fast FT, and nufftn() is purportedly doing a modified FFT.

@Chris Turnes are these times and this speedup within expectations? I really expected nufftn() running a modified FFT to blow away my naive term-by-term FT.

Here is my complete code that calculated these times. It includes the lines you gave me above.

% Resolution of ground plane image
dim1 = 291;
dim2 = 311;

% This block just creates xs & ys, coords of a random parallelogram of
% measured sample points.
theta1 = rand * pi/4;
theta2 = theta1 * (.5 + rand);
M = [2*(1+rand)/3 * [cos(theta1) sin(theta1)]; 2*(1+rand)/3 * [-sin(theta2) cos(theta2)]];
[Xs, Ys] = meshgrid(linspace(-3,3,3000),linspace(-3,3,3000));
gpt = M * [Xs(:)';Ys(:)'];
% Find the measured points that fall within the image.
keep = gpt(1,:)>-1/(2*dim1) & gpt(1,:)<1-1/(2*dim1) & gpt(2,:)>-1/(2*dim2) & gpt(2,:)<1-1/(2*dim2);
xs = gpt(1,keep);
ys = gpt(2,keep);
figure(1); scatter(xs,ys,.5); title('measured points in image')

% Create image pixels
[Xq, Yq] = meshgrid(linspace(0,1-1/(dim1),dim1),linspace(0,1-1/(dim2),dim2));
%Ratio of measured points within image to image pixels.
oversample = sum(keep) / (dim1 * dim2)
% 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]]';
% For this toy problem, use a closed form function in place of measured 
% data.  Calclate this function first for our image points.
td = cos(3*Xq-4*Yq);

tic
% Calculate the frequency domain, from our td above and for our toy 
% measured data.
jfd = zeros(size(td));
nufd = zeros(size(td));
fs = cos(3*xs-4*ys);
for j = 1:dim1
    for k = 1:dim2
        % Calculate the non-fast FFT, just for comparison
        jfd(k,j) = sum(sum(td          .* exp(-1i*2*pi*(kx(j)*Xq+ky(k)*Yq))));
        % Calculate same terms, but from our measured samples.
        pfd(k,j) = sum(fs .* exp(-1i*2*pi*(kx(j)*xs+ky(k)*ys)));
    end
end
pfd = pfd / oversample;
toc

tic
jfd2 = reshape(nufftn(td, { Yq(:,1), Xq(1,:) },    { ky kx }),   dim2, dim1);
pfd2 = reshape(nufftn(cos(3*xs-4*ys), [ys(:) xs(:)], { ky kx }), dim2, dim1) / oversample;
toc

figure(39); surf(Xq,Yq,td); title('orig function');
figure(40); 
surf(Xq,Yq,ifft2(jfd)); 
title('ifft2 of non-fast FT for comparison');
figure(41); surf(Xq,Yq,real(ifft2(pfd))); title('real of ifft2 of pfd');
figure(50); 
surf(Xq,Yq,real(ifft2(jfd2))); 
title('ifft2 of nufft for comparison');
figure(51); surf(Xq,Yq,real(ifft2(pfd2))); title('real of ifft2 of nufft');

[Jerry Guern]

@Chris Turnes These run times really are surprisingly slow. Are we sure it's running an FFT internally? I can't use the MATLAB version if it's really this slow. Another issue is, this code only works if dim1 and dim2 are odd. I know that's a fft2() thing, but I don't kow how to modify my inputs to nufftn() to accomodate it.

[Chris Turnes]

@Jerry Guern - this is because nufftn is having a hard time recognizing that the grids ky and kx are uniform and is using a slow algorithm as a result. You can give nufftn a little nudge to help it, namely, by putting kx and ky in monotonically increasing order:

tic
jfd3 = ifftshift(reshape(nufftn(td, { Yq(:,1), Xq(1,:) },    { fftshift(ky) fftshift(kx) }),   dim2, dim1));
pfd3 = ifftshift(reshape(nufftn(cos(3*xs-4*ys), [ys(:) xs(:)], { fftshift(ky) fftshift(kx) }), dim2, dim1)) / oversample;
toc

When I run that on MATLAB online, it cuts the runtime from 170 seconds to 0.13 s.

Please know that this performance issue is on our radar and we're working to address it.

[Jerry Guern]

@Chris Turnes Well, that dramatically improved my run time, but now it gives a totally wrong solution. This is a serious bug if I'm understanding correctly. The fast algorithm assumes k's in the wrong order?

[Jerry Guern]

@Chris Turnes

[Chris Turnes]

@Jerry Guern -- the algorithm gives the correct answers for the inputs it is given. However, if nufftn sees that the points are in uniform order, it can use faster algorithms. The detection of uniform points requires them to be monotonically increasing or decreasing though; in your case, you have something like [0:N -N:-1], which has a discontinuity, so it doesn't detect the points as uniform.

When I ran your code and used fftshift to put the points into the order [-N:N] instead (which also requires undoing this permutation on the output), I got floating point round-off differences on the order of 1e-9 and 1e-5, respectively, which the order of difference I'd expect to see based on the algorithm choices. I've run it below just to demonstrate that, at least in the example from above, this is giving results with reasonable round-off. If you're seeing significantly different results, can you share the exact code you're running?

% Resolution of ground plane image
dim1 = 291;
dim2 = 311;
% This block just creates xs & ys, coords of a random parallelogram of
% measured sample points.
theta1 = rand * pi/4;
theta2 = theta1 * (.5 + rand);
M = [2*(1+rand)/3 * [cos(theta1) sin(theta1)]; 2*(1+rand)/3 * [-sin(theta2) cos(theta2)]];
[Xs, Ys] = meshgrid(linspace(-3,3,3000),linspace(-3,3,3000));
gpt = M * [Xs(:)';Ys(:)'];
% Find the measured points that fall within the image.
keep = gpt(1,:)>-1/(2*dim1) & gpt(1,:)<1-1/(2*dim1) & gpt(2,:)>-1/(2*dim2) & gpt(2,:)<1-1/(2*dim2);
xs = gpt(1,keep);
ys = gpt(2,keep);
% Create image pixels
[Xq, Yq] = meshgrid(linspace(0,1-1/(dim1),dim1),linspace(0,1-1/(dim2),dim2));
%Ratio of measured points within image to image pixels.
oversample = sum(keep) / (dim1 * dim2)

oversample = 2.3694

% 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]]';
% For this toy problem, use a closed form function in place of measured 
% data.  Calclate this function first for our image points.
td = cos(3*Xq-4*Yq);
tic
% Calculate the frequency domain, from our td above and for our toy 
% measured data.
jfd = zeros(size(td));
nufd = zeros(size(td));
fs = cos(3*xs-4*ys);
for j = 1:dim1
    for k = 1:dim2
        % Calculate the non-fast FFT, just for comparison
        jfd(k,j) = sum(sum(td          .* exp(-1i*2*pi*(kx(j)*Xq+ky(k)*Yq))));
        % Calculate same terms, but from our measured samples.
        pfd(k,j) = sum(fs .* exp(-1i*2*pi*(kx(j)*xs+ky(k)*ys)));
    end
end
pfd = pfd / oversample;
toc

Elapsed time is 159.395508 seconds.

tic
jfd2 = ifftshift(reshape(nufftn(td, { Yq(:,1), Xq(1,:) },    { fftshift(ky) fftshift(kx) }),   dim2, dim1));
pfd2 = ifftshift(reshape(nufftn(cos(3*xs-4*ys), [ys(:) xs(:)], { fftshift(ky) fftshift(kx) }), dim2, dim1) / oversample);
toc

Elapsed time is 0.098690 seconds.

norm(jfd - jfd2, 'fro')

ans = 3.2676e-09

norm(pfd - pfd2, 'fro')

ans = 6.7795e-05

[Jerry Guern]

@Chris Turnes Yes, that is indeed much faster, and my timing results are similar to yours. But the result from nufftn() is completely wrong now. Look at the figures the code produces. It's shaped nothing like the cos(3*Xq-4*Yq) surface.

[Chris Turnes]

@Jerry Guern - that's what I'm trying to understand. The actual values of pfd and pfd2 from my script above are only different by very small amounts (~7e-5). So I'm not sure how the plot you generated looks substantially different from the same plot generated with pfd given that its values are on the order of 1.

When I do a surface plot of pfd and pfd2 as computed with the fftshifted version side by side, they don't look different:

subplot(2,2,1); surf(Xq,Yq,td); title('orig function');
subplot(2,2,2); surf(Xq,Yq,ifft2(jfd)); title('ifft2 of non-fast FT for comparison');
subplot(2,2,3); surf(Xq,Yq,real(ifft2(pfd))); title('real of ifft2 of pfd');
subplot(2,2,4); surf(Xq,Yq,real(ifft2(pfd2))); title('real of ifft2 of pfd2');

My point is that if it doesn't match the result you're expecting it's not because of the output of nufftn -- it's producing something equivalent to the loop calculation.

[Jerry Guern]

@Chris Turnes Then what am I doing wrong here?

I tidied up all the code, took out the stuff we no longer need, and picked params that get the point across but have a more reasonable run time.

Figure 1 is the true function. Figure 2 is my term-by-term version, slow but correct. Figure 3 is my first version with nufftn, also slow but correst. Figure 4 is the same thing but with the k values in order, very fast but the result is a complete mess.

I'm not going to be able to figure it out myself, because I'm just working with different cannibalized versions of the code you gave me, which I don't fully understand.

So, what must I fix in my code that generates Figure 4, to make it fast but get a correct result?

I feel like we're sooo close to a success here.

% Resolution of desired image
dim1 = 191;
dim2 = 211;

% This block just creates xs & ys, coords of a random parallelogram of
% measured sample points.
theta1 = rand * pi/4;
theta2 = theta1 * (.5 + rand);
M = [2*(1+rand)/3 * [cos(theta1) sin(theta1)]; 2*(1+rand)/3 * [-sin(theta2) cos(theta2)]];
[Xs, Ys] = meshgrid(linspace(-3,3,2000),linspace(-3,3,2000));
gpt = M * [Xs(:)';Ys(:)'];
% Find the measured points that fall within the image.
keep = gpt(1,:)>-1/(2*dim1) & gpt(1,:)<1-1/(2*dim1) & gpt(2,:)>-1/(2*dim2) & gpt(2,:)<1-1/(2*dim2);
xs = gpt(1,keep);
ys = gpt(2,keep);

% Create image pixels
[Xq, Yq] = meshgrid(linspace(0,1-1/(dim1),dim1),linspace(0,1-1/(dim2),dim2));

%Ratio of measured points within image to image pixels.
oversample = sum(keep) / (dim1 * dim2)

% 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]]';
% Use a closed form function in place of measured data.
td = cos(3*Xq-4*Yq); % Image grid points
fs = cos(3*xs-4*ys); % Sampled data points
figure(1); surf(Xq,Yq,td); title('orig function');

tic
% Calculate the frequency domain from our toy measured data.
nufd = zeros(size(td));
for j = 1:dim1
    for k = 1:dim2
        pfd1(k,j) = sum(fs .* exp(-1i*2*pi*(kx(j)*xs+ky(k)*ys)));
    end
end
pfd1 = pfd1 / oversample;
uu = toc
figure(2); surf(Xq,Yq,real(ifft2(pfd1))); 
title_str = sprintf('real of ifft2 of term-by-term naive FT\nelapsed time = %f',uu);
title(title_str);

tic
nufft_out1 = reshape(nufftn(fs, [ys(:) xs(:)], { ky kx }), dim2, dim1) / oversample;
uu = toc
figure(3); surf(Xq,Yq,real(ifft2(nufft_out1))); 
title_str = sprintf('real of ifft2 of nufft ver1\nelapsed time = %f',uu);
title(title_str);

tic
nufft_out2 = reshape(nufftn(fs, [ys(:) xs(:)], { sort(ky) sort(kx) }), dim2, dim1) / oversample;
uu = toc
figure(4); surf(Xq,Yq,real(ifft2(nufft_out2)));
title_str = sprintf('real of ifft2 of nufft with sorted k values\nelapsed time = %f',uu);
title(title_str);

[Jerry Guern]

@Chris Turnes btw, the code also still only works for odd dim1 and dim2, so there's something else going on too that I don't know how to fix.

[Paul]

I believe that in the last case we have to ifftshift before taking ifft2. See below.

% Resolution of desired image
dim1 = 19; %191;
dim2 = 21; %211;
% This block just creates xs & ys, coords of a random parallelogram of
% measured sample points.
theta1 = rand * pi/4;
theta2 = theta1 * (.5 + rand);
M = [2*(1+rand)/3 * [cos(theta1) sin(theta1)]; 2*(1+rand)/3 * [-sin(theta2) cos(theta2)]];
[Xs, Ys] = meshgrid(linspace(-3,3,2000),linspace(-3,3,2000));
gpt = M * [Xs(:)';Ys(:)'];
% Find the measured points that fall within the image.
keep = gpt(1,:)>-1/(2*dim1) & gpt(1,:)<1-1/(2*dim1) & gpt(2,:)>-1/(2*dim2) & gpt(2,:)<1-1/(2*dim2);
xs = gpt(1,keep);
ys = gpt(2,keep);
% Create image pixels
[Xq, Yq] = meshgrid(linspace(0,1-1/(dim1),dim1),linspace(0,1-1/(dim2),dim2));
%Ratio of measured points within image to image pixels.
oversample = sum(keep) / (dim1 * dim2)

oversample = 410.7970

% 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]]';
% Use a closed form function in place of measured data.
td = cos(3*Xq-4*Yq); % Image grid points
fs = cos(3*xs-4*ys); % Sampled data points
figure(1); surf(Xq,Yq,td); title('orig function');zlim([-1.5,1.5])

tic
% Calculate the frequency domain from our toy measured data.
nufd = zeros(size(td));
for j = 1:dim1
    for k = 1:dim2
        pfd1(k,j) = sum(fs .* exp(-1i*2*pi*(kx(j)*xs+ky(k)*ys)));
    end
end
pfd1 = pfd1 / oversample;
uu = toc

uu = 0.3626

figure(2); surf(Xq,Yq,real(ifft2(pfd1))); 
title_str = sprintf('real of ifft2 of term-by-term naive FT\nelapsed time = %f',uu);
title(title_str);

tic
nufft_out1 = reshape(nufftn(fs, [ys(:) xs(:)], { ky kx }), dim2, dim1) / oversample;
uu = toc

uu = 0.4512

figure(3); surf(Xq,Yq,real(ifft2(nufft_out1))); 
title_str = sprintf('real of ifft2 of nufft ver1\nelapsed time = %f',uu);
title(title_str);

tic
nufft_out2 = reshape(nufftn(fs, [ys(:) xs(:)], { sort(ky) sort(kx) }), dim2, dim1) / oversample;
uu = toc

uu = 0.0765

%figure(4); surf(Xq,Yq,real(ifft2(nufft_out2)));
figure(4); surf(Xq,Yq,real(ifft2(ifftshift(nufft_out2))));
title_str = sprintf('real of ifft2 of nufft with sorted k values\nelapsed time = %f',uu);
title(title_str);

When dim1 is odd, kx will have the numel(kx) = dim1

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

kx = 1×11

0 1 2 3 4 5 -5 -4 -3 -2 -1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

But, when kx is even, we get the same result

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

kx = 1×11

0 1 2 3 4 5 -5 -4 -3 -2 -1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

I'm going to speculate that when dim1 is even we really want

kx = [0, 1:dim1/2-1, -dim1/2, -dim1/2+1:-1]

kx = 1×10

0 1 2 3 4 -5 -4 -3 -2 -1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

I put the negative sign on the Nyquist point so that when we use sort(kx) in case 4 we get a frequency vector consistent with what's assumed by ifftshift when N is even

clearvars

% Resolution of desired image
dim1 = 18; %191;
dim2 = 20; %211;
% This block just creates xs & ys, coords of a random parallelogram of
% measured sample points.
theta1 = rand * pi/4;
theta2 = theta1 * (.5 + rand);
M = [2*(1+rand)/3 * [cos(theta1) sin(theta1)]; 2*(1+rand)/3 * [-sin(theta2) cos(theta2)]];
[Xs, Ys] = meshgrid(linspace(-3,3,2000),linspace(-3,3,2000));
gpt = M * [Xs(:)';Ys(:)'];
% Find the measured points that fall within the image.
keep = gpt(1,:)>-1/(2*dim1) & gpt(1,:)<1-1/(2*dim1) & gpt(2,:)>-1/(2*dim2) & gpt(2,:)<1-1/(2*dim2);
xs = gpt(1,keep);
ys = gpt(2,keep);
% Create image pixels
[Xq, Yq] = meshgrid(linspace(0,1-1/(dim1),dim1),linspace(0,1-1/(dim2),dim2));
%Ratio of measured points within image to image pixels.
oversample = sum(keep) / (dim1 * dim2)

oversample = 383.3556

% 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]]';

kx = [0, 1:dim1/2-1, -dim1/2, -dim1/2+1:-1];
ky = [0, 1:dim2/2-1, -dim2/2, -dim2/2+1:-1]';


% Use a closed form function in place of measured data.
td = cos(3*Xq-4*Yq); % Image grid points
fs = cos(3*xs-4*ys); % Sampled data points
figure; surf(Xq,Yq,td); title('orig function');zlim([-1.5,1.5])

tic
% Calculate the frequency domain from our toy measured data.
nufd = zeros(size(td));
for j = 1:dim1
    for k = 1:dim2
        pfd1(k,j) = sum(fs .* exp(-1i*2*pi*(kx(j)*xs+ky(k)*ys)));
    end
end
pfd1 = pfd1 / oversample;
uu = toc

uu = 0.2969

figure; surf(Xq,Yq,real(ifft2(pfd1))); 
title_str = sprintf('real of ifft2 of term-by-term naive FT\nelapsed time = %f',uu);
title(title_str);

tic
nufft_out1 = reshape(nufftn(fs, [ys(:) xs(:)], { ky kx }), dim2, dim1) / oversample;
uu = toc

uu = 0.3502

figure; surf(Xq,Yq,real(ifft2(nufft_out1))); 
title_str = sprintf('real of ifft2 of nufft ver1\nelapsed time = %f',uu);
title(title_str);

tic
nufft_out2 = reshape(nufftn(fs, [ys(:) xs(:)], { sort(ky) sort(kx) }), dim2, dim1) / oversample;
uu = toc

uu = 0.0480

%figure(4); surf(Xq,Yq,real(ifft2(nufft_out2)));
figure; surf(Xq,Yq,real(ifft2(ifftshift(nufft_out2))));
title_str = sprintf('real of ifft2 of nufft with sorted k values\nelapsed time = %f',uu);
title(title_str);

[Paul]

@Chris Turnes

I have additional concerns about the correctness of the doc page at nufftn and therefore have more questions about nufftn implementation.

Consider a small set of sample points. Define some query points.

x = 0:5; y = 0:7; f1 = 0:10; f2 = 0:11; f = {f1,f2};

Grid the sample points using ndgrid

[X1,Y1] = ndgrid(x,y);

and compute the input array.

Z1 = X1.^2 + Y1;

According to the doc page, if using {} on the sample points: "When t is specified as a cell array of D vectors, the length of each vector must equal the length of the corresponding dimension of X."

Let's check

[size(Z1,1) == numel(x),size(Z1,2) == numel(y)]

ans = 1x2 logical array
   1   1

We're good according to the doc page, so we have

nz1a = nufftn(Z1,{x,y},f);

Or, we can string out the sample points pairwise (we are still using the input array as a matrix)

nz1b = nufftn(Z1,[X1(:),Y1(:)],f);

The two results are the same.

isequal(nz1a,nz1b)

ans = logical
   1

And we see that nz1a and nz1b are both returned as column vectors, contra to the doc page when the query points are specified with {} as we previously discussed upstream

[size(nz1a),size(nz1b)]

ans = 1×4

132 1 132 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Is it true that nz1a and nz1b are returned in the same order as one would get if we ndgridded f1 and f2 and then strung out the grid arrays pairwise, i.e., [F1,F2]=ndgrid(f1,f2), [F1(:), F2(:)] maps to nz1a(:) ?

Now use meshgrid to grid the sample points

[X2,Y2] = meshgrid(x,y);
Z2 = X2.^2 + Y2;

Now we don't satisfy the doc page criteria

[size(Z2,1) == numel(x),size(Z2,2) == numel(y)]

ans = 1x2 logical array
   0   0

but nufftn is happy to process the input array anyway instead of throwing an error when using {x,y} for the sample points

nz2a = nufftn(Z2,{x,y},f);

Using pairs for the sample points with [] is programmatically correct because all that's required for [] is that "the number of rows of t must equal the number of elements in X."

nz2b = nufftn(Z2,[X2(:),Y2(:)],f);

nz1b and nz2b are equal

isequal(nz1b,nz2b)

ans = logical
   1

Conceptually, I can undestand that result IF in those cases nufftn maps [X1(:),Y1(:)] to Z1(:) and maps [X2(:),Y2(:)] to Z2(:), because those maps are the same, albeit in a different order.

But what I can't understand is how nz2a also returns the same result

isequal(nz1a,nz2a)

ans = logical
   1

In the calls to nufffn for the nz1a and nz2a cases, the sample point and query point input arguments are the same, but Z1 and Z2 are not

isequal(Z1,Z2)

ans = logical
   0

Z1 and Z2 have the same elements, but in a different order

isequal(sort(Z1(:)),sort(Z2(:)))

ans = logical
   1

But nufftn returns the same results by somehow mapping the same x,y pairs to the same values in Z1 and Z2 that are in different locations in Z1 and Z2?

At first I thought that maybe nufftn is comparing the lengths of x and y to the dimensions of the input array to get the ndgrid vs. meshgrid sorted out, but that would be ambiguous if the input array is square.

[Jerry Guern]

@Chris Turnes Okay, my code now seems to do just what I wanted. I'll stay tuned to see the resolution of @Paul's concern, but I think I'm good. I'm going to paste my latest code below, because this version is maximum clarity and addresses all the concerns I had going in and the concerns discovered along the way. When they revamp the documentation, I think giving a working example this complete would be very helpful to future nufftn() users.

% Resolution of desired image
dimy = 111;
dimx = 122;

% This block just creates xs & ys, coords of a random parallelogram of
% toy measured sample points that fall within the image area.
theta1 = rand * pi/4;
theta2 = theta1 * (.5 + rand);
M = [2*(1+rand)/3 * [cos(theta1) sin(theta1)]; 2*(1+rand)/3 * [-sin(theta2) cos(theta2)]];
[Xs, Ys] = meshgrid(linspace(-3,3,2000),linspace(-3,3,2000));
gpt = M * [Xs(:)';Ys(:)'];
% Find the measured points that fall within the image.
keep = gpt(1,:)>-1/(2*dimx) & gpt(1,:)<1-1/(2*dimx) & gpt(2,:)>-1/(2*dimy) & gpt(2,:)<1-1/(2*dimy);
xs = gpt(1,keep);
ys = gpt(2,keep);

% Create image pixels
[Xq, Yq] = meshgrid(linspace(0,1-1/(dimx),dimx),linspace(0,1-1/(dimy),dimy));
%Ratio of measured points within image to image pixels.
oversample = sum(keep) / (dimx * dimy)
% Wave numbers for Fourier Transform
kx = [ceil(-dimx/2) : ceil(dimx/2)-1];
ky = [ceil(-dimy/2) : ceil(dimy/2)-1];
% Use a closed form function in place of measured data for this toy problem
td = cos(3*Xq-4*Yq); % Image grid points
fs = cos(3*xs-4*ys); % Sampled data points
figure(1); surf(Xq,Yq,td); title('orig function');

% Run the nufftn() and display results
tic
nufft_out2 = reshape(nufftn(fs, [ys(:) xs(:)], { ky kx }), dimy, dimx) / oversample;
uu = toc
figure(4); surf(Xq,Yq,real(ifft2(ifftshift(nufft_out2))));
title_str = sprintf('real of ifft2 of ifftshift of nufft\nelapsed time = %f',uu);
title(title_str);

[Chris Turnes]

@Paul:

Is it true that nz1a and nz1b are returned in the same order as one would get if we ndgridded f1 and f2 and then strung out the grid arrays pairwise, i.e., [F1,F2]=ndgrid(f1,f2), [F1(:), F2(:)] maps to nz1a(:) ?

Yes, but let me modify your example slightly to be fair:

x = 0:5; y = 0:7; f1 = (0:10)/11; f2 = (0:11)/12; f = {f1,f2};
[X1,Y1] = ndgrid(x,y);
Z1 = X1.^2 + Y1;
nz1a = nufftn(Z1,{x,y},f);
nz1b = nufftn(Z1,[X1(:),Y1(:)],f);
[F1,F2] = ndgrid(f1, f2);
nz1c = nufftn(Z1,[X1(:),Y1(:)],[F1(:) F2(:)]);
max(abs(nz1a - nz1c))

ans = 8.5265e-13

Now we don't satisfy the doc page criteria but nufftn is happy to process the input array anyway instead of throwing an error when using {x,y} for the sample points

Yes, this is essentially the same discontinuity between behavior and documentation but for the sample points rather than the query points. The current behavior is basically always vectorizing input and output whereas the doc suggests different shapes for input/output. This is an issue we need to address.

Using pairs for the sample points with [] is programmatically correct because all that's required for [] is that "the number of rows of t must equal the number of elements in X." nz1b and nz2b are equal. Conceptually, I can undestand that result IF in those cases nufftn maps [X1(:),Y1(:)] to Z1(:) and maps [X2(:),Y2(:)] to Z2(:), because those maps are the same, albeit in a different order. But what I can't understand is how nz2a also returns the same result.

Your particular example is just leading you astray. The reason why for this example everything is always exactly equal is that you've chosen integer values for both the sample and query points, which means every single exponential evaluation takes the form for integer m and n, which is always 1. So the reason they match is because every element of the output is identical and is equal to the DC value:

x = 0:5; y = 0:7; f1 = (0:10); f2 = (0:11); f = {f1,f2};
[X1,Y1] = ndgrid(x,y);
Z1 = X1.^2 + Y1;
nz1a = nufftn(Z1,{x,y},f);
numel(unique(nz1a))

ans = 1

sum(Z1(:)) == unique(nz1a)

ans = logical
   1

If you use fractional frequencies, things will make more sense:

x = 0:5; y = 0:7; f1 = (0:10)/11; f2 = (0:11)/12; f = {f1,f2};
[X1,Y1] = ndgrid(x,y);
[X2,Y2] = meshgrid(x,y);
Z1 = X1.^2 + Y1;
Z2 = X2.^2 + Y2;
nz1a = nufftn(Z1, {x y}, f);
nz1b = nufftn(Z1, [X1(:) Y1(:)], f);
max(abs(nz1b - nz1a))

ans = 2.2388e-12

% We have to carefully manage the shape to make things line up.
nz2a = nufftn(Z2', {x y}, f);
max(abs(nz1a - nz2a))

ans = 0

% Alternatively, flip x and y, and f2 and f1, then
% reshape, transpose, and flatten the output.
nz2a = reshape(nufftn(Z2, {y x}, {f2 f1}), [numel(f2) numel(f1)]).';
nz2a = nz2a(:);
max(abs(nz1a - nz2a))

ans = 0

% Naively passing in won't work
nz2a_wrong = nufftn(Z2, {x y}, f);
max(abs(nz1a - nz2a_wrong))

ans = 263.7607

% But passing in the point list works, because we are specifying the
% correct coordinates for each sample.
nz2b = nufftn(Z2, [X2(:) Y2(:)], f);
max(abs(nz1a - nz2b))

ans = 2.3020e-12

Generally, because of the confusion that can result with meshgrid, ndgrid is simpler to use.

[Paul]

I was wondering if my example was the problem. Thanks for pointing that out. I was more focused on the fact that Z1 and Z2 had different dimensions but nufftn didn't seem to care.

In summary, to make sure I'm clear:

If the input sample points are specified as a cell array {} of D vectors, then the input array must be either

a) an N-D array as would be formed from operating on the ndgrid of those D vectors, or

b) the same as (a) and then reshaped into a vector

x = 0:5; y = 0:7; z = 0:.1:1; f1 = (0:10)/11; f2 = (0:11)/12; f3 = (0:4)/5; f = {f1,f2,f3};
[X1,Y1,Z1] = ndgrid(x,y,z);
A1 = X1.^2 + Y1 + 5*Z1;
out1 = nufftn(A1,{x,y,z},f);
out2 = nufftn(A1(:),{x,y,z},f);
out3 = nufftn(A1(:).',{x,y,z},f);
isequal(out1,out2,out3)

ans = logical
   1

I'll still suggest that if the sample points are specified with {} and the input array is a matrix or N-D array, then nufftn should throw an error if the dimensions of the input array don't line up with the lengths of the sample point vectors.

If the query points are specified as a cell array {} of D vectors, then the output will be a column vector formed (approximately) as if we are first forming the N-D array from an ndgrid of query points and then reshaping into a column vector

[F1,F2,F3] = ndgrid(f1,f2,f3);
out1 = nufftn(A1,{x,y,z},f);
out2 = nufftn(A1,{x,y,z},[F1(:),F2(:),F3(:)]);
isequal(out1,out2)

ans = logical
   0

norm(out1-out2)

ans = 7.2127e-11

out1 and out2 are slightly different because nufftn goes down different code paths? (similarly for nz1a and nz1c in your first example?)

"Generally, because of the confusion that can result with meshgrid, ndgrid is simpler to use."

I think I'd stay away from meshgrid altogether when dealing with nufftn.

[Chris Turnes]

Yes, your understanding is correct. Thanks for pointing out the inconsistencies between the doc and behavior, and I'll definitely pass on your feedback.

Answer 2

[Matt J]

I need to Fourier Transform my parallelogram grid of measured values into ground freq domain, the Inv FT that into ground time domain.

I don't know what "ground time/freq" domain means.

I think if you just pretend the samples are Cartesian and take the FFT and IFFT, you will obtain the non-Fourier domain result sampled on your original parallelgoram. You can then just use griddata or scatteredInterpolant to reample that, if needed.

[Jerry Guern]

I'm creating an image in ground plane coordinates. When I said ground time/freq domain, the time domain would be the values of the ground plane image I'm trying to create, and the freq domain is the fft of that. I said ground to contrast it with the slant plane, the coordionates I'm projecting into ground plane. I don't know what you meant by what you said. Could you show me in code?

[Matt J]

I mean,

slanted =  ifft2(  some_transform(  fft2(td)  ) );

ground = griddata(___,slanted,___)