Nufft for scalar function of 2 variables

Question

schlapp 2020-4-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/521022-nufft-for-scalar-function-of-2-variables

编辑： Chris Turnes 2020-4-28

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Chris Turnes 2020-4-28

编辑：Chris Turnes 2020-4-28

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 2020-4-27

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/521022-nufft-for-scalar-function-of-2-variables#answer_428656

编辑：schlapp 2020-4-27

在 MATLAB Online 中打开

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.