Nufft for scalar function of 2 variables

2 次查看（过去 30 天）
schlapp
schlapp 2020-4-27
编辑： 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 个评论
Chris Turnes
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

schlapp
schlapp 2020-4-27
编辑：schlapp 2020-4-27
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.

类别

在 Help Center 和 File Exchange 中查找有关 Fast Fourier Transforms 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by