Nufft for scalar function of 2 variables

I need some help applying the new nufftn function from matlab to calculate Fourier transforms of scalar functions of two variables. Basically I want to transform wavefunctions from 2D-momentum (kx,ky) to 2D-position (x,y) space. For each grid (kx,ky,x,y) I use the same nonuniform grid of a scaled Gauss-Legendre quadrature representing values from -infinity to infinity. The given function values evaluated at my nonuniform sample grid (kx,ky) are given as a 2d matrix f.

I have written my own function gqft before. Its basically an implementation of the Gauss-Legendre quadrature to calculate the Fourier integral. The functions I transform are either purely even or odd, hence I use either the cos, or the sin part of the imaginary exponential. Here is the cos-Version shown for brevity:

function [output] = gqft(f,grid,weights)
% Fourier Transformation with Gausssian quadrature
Npp=length(grid);
output = zeros(Npp,Npp);

parfor ii = 1:Npp
    for jj = 1:Npp
        kernel = f.*cos((grid(jj)*grid + grid(ii)*grid'));% f.*sin((grid(jj)*grid + grid(ii)*grid'));
        output(ii,jj)=sum(weights(Nmin:Nmax)'.*weights(Nmin:Nmax).*kernel(Nmin:Nmax,Nmin:Nmax),'all')/2/pi;
    end
end

%output = sum(f.*cos(grid.*reshape(grid,1,1,[]) + grid'.*reshape(grid,1,1,1,[])).*reshape(weights,1,1,[]).*reshape(weights,1,1,1,[]),[3,4]);
% (requires too much memory)
end

This gives me the desired result, however it is very slow and scales horribly with increased grid size due to the double for-loops. The commented line is a vectorized approach however not usable because of exessive memory requirement (the integrand is basically a 4d matrix).

Recently my university got access to the new matlab version, introducing the functions nufft, nufftn. I would like to use them instead of my own function for performance reasons. To test it, I applied both to the example function

f = exp(-(kx^2 + ky^2 + kx*ky))

Using my own function

ftilde = gqft(f,kx,weights)

yields the correct result, with some errors at the grid borders (close to infinity). However if I restrict the evaluation to +-5 in each direction, the result coincides with the analytical expression.

Using

ftilde = reshape(nufftn(f,{kx,kx'},{kx,kx'}),length(kx),length(kx))

gives me plain zeros. Cropping again everything outside of +-5 however still gives me not the desired result. The function is way too narrow and the function values are way too high.

Using

ftilde = reshape(nufftn(f,{k5,k5'},{k5,k5'}),length(k5),length(k5)); k5 = [-5:5/200:5];

with a uniform grid on +-5 also does not give me the result that I need.

Soo to conclude I either need some advice on how to get my own code to run a lot faster or on how to use nufft in a correct way.

[Chris Turnes]

I think likely the issue here is the scaling of the nodes. The nufftn function applies the formula listed on its documentation page. Typically in FFTs, one domain is "scaled" to have its grid lie within [0 1) while the other repsents some sample indices. Here, you are using the same input and output locations. I'm not sure what the appropriate scaling for your data is, but you'll likely find better results using nufftn if you can define the proper relationship between the input and output domains in terms of their discrete sampling representation.

Answer 1

[schlapp]

I understood now that using

k = [-d:d/Ngrid:d];
kscld = k/sqrt(2*pi);
ftilde = reshape(nufftn(f,{kscld,kscld'},{kscld,kscld'}),length(k),length(k))/2/pi*(d/Ngrid)^2;

yields in fact the same result as my function, however here I am using now a uniform grid. Utilmately I would like to use the same nonuniform grid as in my function.

With the nonuniform grid I seem to get only plain zeros as results.