NUFFT for zeropadded data

Hi all,

I am trying to use the NUFFT package from NYU (https://cims.nyu.edu/cmcl/nufft/nufft.html) with upsampled data from a given function, i.e. I zero-pad the data wich I apply the FFT to. Then I want to evaluate this function with the NUFFT at randomized data locations. I do not manage to get the scaling of the input to the NUFFT correct. I do however get correct results when the target locations are a uniform grid. Please see the code attached below. I know there are other implementations of NUFFT, but I would like to use this one.

Thanks,

Fredrik

%%
% EXAMPLE SCRIPT FOR 2D NUFFT TYPE II, USING CODE WRITTEN BY GREENGARD-LEE
%
% USE TO COMPUTE f(xj,yj), WITH (xj,yj) NONUNIFORM, BY GIVEN FOURIER 
%COEFFICIENTS ON A UNIFORM GRID.
%
% NOTE:
%
%   i, BY GREENGARD-LEE CONVENTION, FFT(X) = fft(X)/length(X),  
% WHERE fft IS MATLABS FFT:
%
%                     N
%       X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
%                    n=1
%
% THE SAME RELATION APPLIES TO IFFT.

clear all
close all
clc
%%
% Testfunction that is sampled to get uniform Fourier coefficients Fkl. Choose periodic.
scale = 0.3;
 f = @(x,y) sin(2*pi*x).*cos(4*pi*y).*exp(-(sqrt((x-17/97).^2 + (y+3/577).^2)/scale).^2);
% f = @(x,y) exp(sin(x+y));
% f = @(x,y) exp(-(sqrt(x.^2 + y.^2)/scale).^2);
% f = @(x,y) sin(2*pi*x).*cos(4*pi*y).*exp(sin(8*pi*x).*sin(3*pi*x));

%%
% The Domain [a,b]x[c,d].
L = 12;
a = -L/2;
b = L/2;
c = a;
d = b;
Lx = (b-a);
Ly = (d-c);
%% Set grid sizes
% Sampling factors
sf = 1; 
M = 2^5; % Number of uniform sampling points in real space
N = M^2; % Number of scattered evaluation points in real space
Msf = M * sf; % Number of Fourier modes
%%
% Points (xj,yj) to evaluate fj = f(xj,yj) at with NUFFT. XX, YY are 
%representations on the unform grid in R^2
% from which the uniformly spaced Fkl are obtained.


xj = a + Lx*rand(N,1);          %xj in [a b]
yj = c + Ly*rand(N,1);          %yj in [a b]

fjexact = f(xj,yj);               %Correct value at (xj,yj)

xx = linspace(a,b,M+1)';
xx(end) = [];
yy = xx;
[XX, YY] = meshgrid(xx,yy);
fXYexact = f(XX,YY);            %Uniform sampling of f on [a,b]x[c,d]
%% FFT
Fk = fftshift(fft2(ifftshift(fXYexact),Msf,Msf)); %Shift in accordance with Greengard_lee convention
%% NUFFT
eps = 1e-14;
iflag = +1;
[fnufft,ier] = nufft2d2(M^2,XX(:)*2*pi/Lx/sf+pi/sf,YY(:)*2*pi/Ly/sf+pi/sf,iflag,eps,M*sf,M*sf,Fk.'/(M*sf)^2);
sftemp = Msf^2/N;
[fnufftnonuniform,ier] = nufft2d2(N,xj*2*pi/Lx/sftemp+pi/sftemp,yj*2*pi/Ly/sftemp+pi/sftemp,iflag,eps,M*sf,M*sf,Fk.'/(M*sf)^2);
fnufft = real(reshape(fnufft,size(XX)));
fnufftnonuniform = real(fnufftnonuniform);

%%
nufffterror = abs(fXYexact-fftshift(fnufft));
nufftnonuniformerror = abs(fjexact - fftshift(fnufftnonuniform));
disp(' ')
disp(['Absolute max error uniform: ' num2str(max(nufffterror(:))) ])
disp(['Absolute max error random: ' num2str(max(nufftnonuniformerror(:))) ])
disp(' ')

figure(1)
plot3(xj,yj,fjexact,'ro')
hold on
plot3(xj,yj,ifftshift(fnufftnonuniform),'bx')
hold off

figure(2)
plot3(xj,yj,log10(abs(fjexact - fftshift(fnufftnonuniform))),'k^')