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

[razgzy]

Thanks a lot for your detailed reply.
I tried the MathWork's nufft because I need 2d type 1 nufft from nonuniform space domain to uniform frequency domain. But it seems that standard softwares implement nonuniform frequency domain to uniform space domain. I'd appreciate it if you could tell me where i can find a 2d type 1 nufft algorithm for nonuniform space domain implemented in python GPU.
Best wish

[razgzy]

t1 = np.arange(0,300,1)
t2 =np.arange(500.5,700.5,1)
x = np.concatenate([t1,t2]).astype(np.float64)
c = 2*np.sin(0.1*np.pi*x) + np.sin(0.02*np.pi*x)
c = c.astype(np.complex128)
f = np.arange(0,500,1).astype(np.float64) / 500 
Y = finufft.nufft1d3(x, c, f * 2 * np.pi)
plt.figure()
plt.plot(f, np.abs(Y))
plt.show()

I tried the nufft1d3 and got the similar results shown in matlab's example.
It seems I misunderstood the type 1 nfft. What i need is 2d type3 nfft. But cufinufft hasn't implement it yet.

[ahbarnett]

So, are you ok now? It seems to be a matter of reading the docs and matching to your needs... Best, A

…

On Wed, Jul 12, 2023 at 10:50 PM razgzy ***@***.***> wrote: t1 = np.arange(0,300,1) t2 =np.arange(500.5,700.5,1) x = np.concatenate([t1,t2]).astype(np.float64) c = 2*np.sin(0.1*np.pi*x) + np.sin(0.02*np.pi*x) c = c.astype(np.complex128) f = np.arange(0,500,1).astype(np.float64) / 500 Y = finufft.nufft1d3(x, c, f * 2 * np.pi) plt.figure() plt.plot(f, np.abs(Y)) plt.show() I tried the nufft1d3 and got the similar results shown in matlab's example. It seems I misunderstood the type 1 nfft — Reply to this email directly, view it on GitHub <#316 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNZRSUPSC4W4X4UTQEG7YDXP5PA7ANCNFSM6AAAAAA2IHUJ74> . You are receiving this because you commented.Message ID: ***@***.*** com>

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