Nufft for scalar function of 2 variables
3 views (last 30 days)
Show older comments
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 Comment
Chris Turnes
on 28 Apr 2020
Edited: Chris Turnes
on 28 Apr 2020
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.
See Also
Categories
Find more on Fourier Analysis and Filtering in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!