Replies: 3 comments
-
Dear Kaya, This is very doable, using the idea that "shifting the frequency index is equivalent to post-multiplying by a phase in real space", for a type-2. Make sure you understand the math of that. It actually took me a while (90 mins) to write this 1D demo since I'm out of Python practice, but also to get all the phases correct. It's not trivial. But it is the right way to proceed also in 2D. Make sure you understand this 1D code first: # demo 1D NUFFT interpolation of real periodic func *without* length-doubling
# Barnett 8/13/25
import numpy as np
import finufft as fi
N = 5 # num samples
# a generic real test fun with bandlimit <= (N-1)/2 so interp exact...
fun = lambda t: 1.0 + np.sin(t+1) + np.sin(2*t-2) # bandlimit is 2
Nt = 100 # test targs
targs = np.random.rand(Nt)*2*np.pi
Nf = N//2 + 1 # num freq outputs for rfft; note rfft2 looks different :(
g = np.linspace(0,2*np.pi,N,endpoint=False) # sample grid
f = fun(g)
c = (1/N) * np.fft.rfft(f) # get coeffs 0,1,..Nf-1 (don't forget prefac)
assert c.size==Nf
# check the naive (double-length c array) NUFFT version:
cref = np.concatenate([np.conj(np.flip(c[1:])), c]) # reflect to 1-Nf...Nf-1 coeffs
ft = np.real(fi.nufft1d2(targs,cref,eps=1e-12,isign=1)) # f at targs (isign!)
# (taking Re here was just a formality; it is already real to eps_mach)
print("naive (reflected) 1d2 max err:", np.linalg.norm(fun(targs) - ft, np.inf))
# now demo avoid doubling the NUFFT length via freq shift and mult by phase:
c[1:] *= 2.0 # since each nonzero coeff appears twice in reflected array
N0 = Nf//2 # starting freq shift that FINUFFT interprets the c array
ftp = fi.nufft1d2(targs,c,eps=1e-12,isign=1) # f at targs but with phase error
# the key step: rephase (to account for shift), only then take Re (needed!)...
ft = np.real( ftp * (np.cos(N0*targs) + 1j*np.sin(N0*targs))) # guess 1j sign
print("unpadded 1d2 max err:", np.linalg.norm(fun(targs) - ft, np.inf))
The 2D version should be similar, but you can still only save a factor of two in the nufft2d2 input array size (hence FFT size), just as in 1D. It would need a rectangular (N/2 by N) array, where N is the number of regular sample points per dimension. Handling the origin and the coordinate lines will be a a bit tricky. The other difference is that rfft2 seems to spit out a full (Hermitian-symm) array (equivalent of Let me know if you're still stuck, after you parse and try out the above. PS I'm making this into a new tutorial. Best, Alex |
Beta Was this translation helpful? Give feedback.
-
Dear Alex, Thank you! If you are making this into a Python tutorial, perhaps I should share these tests if you needed a refresher on Python's FFT conventions. Yes previously, I was modifying the Hermitian array by multiplying most elements by 2 (but this depends on whether the array is even or odd) as shown here and padding all the positive frequencies to zero before sending it to Flatiron's nufft2 here. I'll share the 2D version without length doubling here once I get it tested and working. |
Beta Was this translation helpful? Give feedback.
-
https://finufft.readthedocs.io/en/latest/tutorial/realinterp1d.html
I don't have time to do the 2d version, but it looks like you'll be able
to, given the 1d version and your understanding... Let me know! (you can PR
to add to the tutorial).
Best, Alex
…On Wed, Aug 13, 2025 at 1:28 PM Kaya Unalmis ***@***.***> wrote:
Dear Alex,
Thank you! If you are making this into a Python tutorial, perhaps I should
share these tests
<https://github.com/PlasmaControl/DESC/blob/9341e97aa08e9ca3440a44aa2883411265644e9d/tests/test_interp_utils.py#L185>
if you needed a refresher on Python's FFT conventions.
Yes previously, I was modifying the Hermitian array by multiplying most
elements by 2 (but this depends on whether the array is even or odd) as
shown here
<https://github.com/PlasmaControl/DESC/blob/9341e97aa08e9ca3440a44aa2883411265644e9d/tests/test_interp_utils.py#L249>
and padding all the positive frequencies to zero
<https://github.com/PlasmaControl/DESC/blob/9341e97aa08e9ca3440a44aa2883411265644e9d/desc/integrals/_interp_utils.py#L422>
before sending it to Flatiron's nufft2 here
<https://github.com/PlasmaControl/DESC/blob/9341e97aa08e9ca3440a44aa2883411265644e9d/desc/integrals/_interp_utils.py#L419>
.
I'll share the 2D version without length doubling here once I get it
tested and working.
—
Reply to this email directly, view it on GitHub
<#720 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/ACNZRSVVJI4LPQPSFKJ2DNT3NNYVTAVCNFSM6AAAAACDVKKNYCVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTIMBZGU3DAMI>
.
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
|
Beta Was this translation helpful? Give feedback.
Uh oh!
There was an error while loading. Please reload this page.
Uh oh!
There was an error while loading. Please reload this page.
-
Does FINUFFT support type 2 nuffts where the input frequencies are purely positive along one dimension? The application is interpolation of real valued functions. Often we have Fourier series coefficients from, for example, numpy's real FFT
c = np.rfft2(f)
and we would like to evaluate these series at non uniform points with a nufftfq=nufft2_2d(c).real
. Right now we have to pad these coefficients with zeros at all the negative frequencies which doubles the size ofc
. I think the padding becoems more expensive when the precision parameter is chosen to be small.Beta Was this translation helpful? Give feedback.
All reactions