Python Type1 nufft can not handle nonuniform signal out of (-3pi, 3pi) ?

[razgzy]

t1 = np.arange(0,300,1)
t2 =np.arange(500.5,700.5,1)
t = np.concatenate([t1,t2])
S = 2*np.sin(0.1*np.pi*t) + np.sin(0.02*np.pi*t)
Y = finufft.nufft1d1(t, S, 500)
f = np.arange(0,500,1)/500
plt.figure()
plt.plot(f, Y)
plt.show()

I tried to test the example given in Matlab:
https://www.mathworks.com/help/matlab/ref/double.nufft.html
But I failed with RuntimeError: FINUFFT spreader nonuniform point out of range [-3pi,3pi]^d in type 1 or 2

[ahbarnett]

Dear razgzy, Looks like the error message is self-explanatory. See: https://finufft.readthedocs.io/en/latest/math.html Our design choice is simply because fmod by 2pi is too slow, so we leave that up to the user. However, a more major problem is that the math definitions of TheMathWorks nufft are strange and do not match standard software (FINUFFT, NFFT3, etc). I don't believe you have the conversion right. You will have to write out the formulae for TMW and our libraries. You will also need to decide if you want type 1 or type 3. This test matlab code shows type 1 and type 3 compared in the two libraries: ``` % compare FINUFFT to The MathWorks (TMW) nufft introduced in 2020a. % Barnett 7/19/22 clear % choose params... isign = -1; % sign of imaginary unit in exponential: matches TMW eps = 1e-9; % requested accuracy o.debug = 0; % choose 1 for timing breakdown text output FFTW_ESTIMATE = bitshift(1,6); FFTW_MEASURE = bitshift(1,0); o.fftw = FFTW_ESTIMATE; % or see fftw3.h o.upsampfac=0; % 0 (auto), 2.0 (default), or 1.25 (low-RAM, small-FFT) o.chkbnds=0; % a few percent faster M = 1e6; % # of NU pts (in all dims) N = 1e5; % # of modes (approx total, used in all dims) maxNumCompThreads(8); % for 8-core laptop % ----------------------------------------- 1D % type 1 x = pi*(2*rand(M,1)-1); c = randn(M,1)+1i*randn(M,1); nt = floor(0.37*N); % pick output mode index to test fe = sum(c.*exp(1i*isign*nt*x)); % exact ans of1 = 1+floor(N/2); % mode index offset (CMCL ord) %of1 = 1; % mode index offset (FFT ordering) fprintf('do FINUFFT 1d type 1... (N=%d, M=%d, tol=%.3g)\n',N,M,eps) f = finufft1d1(x,c,isign,eps,N,o); tic f = finufft1d1(x,c,isign,eps,N,o); fprintf('\t%.3g sec\trel err in F[%d] is %.3g\n',toc,nt,abs((fe-f(nt+of1))/max(f))) fprintf('do Mathworks nufft 1d type 1 (explicit equispaced freq pts)...\n') tic freqtargs = -N/2:N/2-1; % default modeord fm = nufft(c,(1/2/pi)*x,freqtargs); fprintf('\t%.3g sec\trel err in F[%d] is %.3g\n',toc,nt,abs((fe-fm(nt+of1))/max(f))) fprintf('\trel l2 diff btw two output vecs: %.3g\n',norm(f-fm)/norm(f)) % type 3 bandwidth=1000; % small to help TMW; usually this is N to make type 3 hardness sim to type 1 s = bandwidth*rand(N,1); % freq targets for type 3, our convention j = floor(0.37*N); st = s(j); % index to test fe = sum(c.*exp(1i*isign*st*x)); % exact ans fprintf('\ndo FINUFFT 1d type 3... (N,M,tol as above; freq bandwidth=%.3g)\n',bandwidth) tic f = finufft1d3(x,c,isign,eps,s,o); fprintf('\t%.3g sec\trel err in F[%d] is %.3g\n',toc,nt,abs((fe-f(j))/max(f))) fprintf('do Mathworks nufft 1d type 3 (arb freq query pts)...\n') tic fm = nufft(c,(1/2/pi)*x,s); fprintf('\t%.3g sec\trel err in F[%d] is %.3g\n',toc,nt,abs((fe-fm(j))/max(f))) fprintf('\trel l2 diff btw two output vecs: %.3g\n',norm(f-fm)/norm(f)) ``` Sample output (note we're at least 100x faster...): ```

> addpath ../finufft/matlab/ > benchmark_nufft_1d_types1and3

do FINUFFT 1d type 1... (N=100000, M=1000000, tol=1e-09) 0.0243 sec rel err in F[37000] is 2.73e-10 do Mathworks nufft 1d type 1 (explicit equispaced freq pts)... 2.39 sec rel err in F[37000] is 8.48e-09 rel l2 diff btw two output vecs: 3.84e-08 do FINUFFT 1d type 3... (N,M,tol as above; freq bandwidth=1e+03) 0.0483 sec rel err in F[37000] is 7.76e-10 do Mathworks nufft 1d type 3 (arb freq query pts)... [I got bored of waiting...] ``` Hope this helps. Best, Alex PS since not a bug, I'll move to Discussion.

…

On Wed, Jul 12, 2023 at 1:06 AM razgzy ***@***.***> wrote: t1 = np.arange(0,300,1) t2 =np.arange(500.5,700.5,1) t = np.concatenate([t1,t2]) S = 2*np.sin(0.1*np.pi*t) + np.sin(0.02*np.pi*t) Y = finufft.nufft1d1(t, S, 500) f = np.arange(0,500,1)/500 plt.figure() plt.plot(f, Y) plt.show() I tried to test the example given in Matlab: https://www.mathworks.com/help/matlab/ref/double.nufft.html But I failed with RuntimeError: FINUFFT spreader nonuniform point out of range [-3pi,3pi]^d in type 1 or 2 — Reply to this email directly, view it on GitHub <#315>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNZRSWJPGQHR7G4GXS2BODXPYWD7ANCNFSM6AAAAAA2G5TTZQ> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

-- *-------------------------------------------------------------------~^`^~._.~' |\ Alex Barnett Center for Computational Mathematics, Flatiron Institute | \ http://users.flatironinstitute.org/~ahb 646-876-5942