Why n_trans result is shifted for cufinufft

[sijieh23]

Hi, I tried to perform a few NUFFT for interpolation simultaneously using the type-2 transform. Here is my code

import cupy as cp
import cupyx.scipy.fft as cpft
import cufinufft as cnft

d = 3
m = 1024
n = 10

x = cp.linspace(0, 2 * pi, m, endpoint=False)

# generate random data
u = cp.zeros((d, m))
for _ in range(10):
    for j in range(d):
        a = cp.random.randint(-20, 20)
        u[j] += cp.sin(a * x)

# FFT
u_hat = cpft.fftshift(cpft.fft(u, axis=1, norm='forward'))

# define query points
xq = cp.random.uniform(0, x[-1], n)
uq = cp.zeros((d, n))

# interpolation using type-2 NUFFT
plan = cnft.Plan(2, m, d, isign=1, eps=1e-15, dtype=np.complex128, modeord=0)
plan.setpts(xq)
uq[:] = plan.execute(u_hat).real

What I observed is that the interpolation result, uq, is the interpolation of u shifted by one index, i.e., uq[1] is the interpolation of u[0], and uq[2] is u[1]'s, etc. I was expecting uq[0] to be the interpolation of u[0], etc.

Why is this, or am I missing something?

Thanks!

[ahbarnett]

This was actually amusing to debug. You fftshifted along the n_trans axis as well as the freq axis!
Check https://numpy.org/devdocs/reference/generated/numpy.fft.fftshift.html

Note you had several other issues that made a non-working code. Here's the CPU code I hacked from yours:

import numpy as cp
import numpy as np
import scipy.fft as cpft
import finufft as finufft
pi = np.pi


d = 3        # ntrans, not dim!
m = 1024
n = 10        # query pts

x = cp.linspace(0, 2 * pi, m, endpoint=False)

# generate random data
u = cp.zeros((d, m))
for _ in range(10):
    for j in range(d):
        a = cp.random.randint(-20, 20)
        u[j] += cp.sin(a * x)

# FFT
u_hat = cpft.fftshift(cpft.fft(u, axis=1, norm='forward'),axes=(1,))

# define query points
xq = 0.0*cp.random.uniform(0, 2*pi, size=n) + x[1]   # hack to make 2nd x pt
uq = cp.zeros((d, n))

# interpolation using type-2 NUFFT
# (m,) singleton tuple must need trailing comma! :
plan = finufft.Plan(2, (m,), n_trans=d, isign=1, eps=1e-15, dtype=np.complex128, modeord=0)
plan.setpts(xq)
uq = plan.execute(u_hat).real

uq
u

[sijieh23]

What a stupid mistake I made! Thanks for going through the trouble.

Two more questions:

1. Does n_trans always assume to be in the first axis? (In this case, 3)? In other words, it will execute 3 independent NUFFTs, right?
2. If I were to select modeord=1, I don't need the fftshift, right?

[ahbarnett]

The answers are yes and yes. For 1 keep in mind the "independent" NUFFTs are faster than if doing many separate calls. Enjoy, Alex

…

On Mon, May 5, 2025 at 6:36 PM SJH ***@***.***> wrote: What a stupid mistake I made! Thanks for going through the trouble. Two more questions: 1. Does n_trans always assume to be in the first axis? (In this case, 3)? In other words, it will execute 3 independent NUFFTs, right? 2. If I were to select modeord=1, I don't need the fftshift, right? — Reply to this email directly, view it on GitHub <#669 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACNZRSW3SWG4HV2BBI5CBFD247RW3AVCNFSM6AAAAAB4EILSGCVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTGMBUGQ4TKNQ> . 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