pynufft/pynufft

Help with example of non-uniform sample points

Closed this issue · 1 comments

LibrEars commented

Hi @jyhmiinlin,

thank you for the library! I am currently trying to use it to calculate noise spectral densities of current measurements. This could be added to the examples since it compares to the rFFT from numpy.

I am wondering how to use pynufft with nonuniform sample points that would be of type NUFFT-II according to wikipedia.

Here is my sample code to reproduce the issue

# Example uniform (i.e. integer) frequencies f_{k}=k but nonuniform sample points p_{n}. From wikipedpa: This corresponds to evaluating a Fourier series at nonequispaced points. It is also known as adjoint NDFT. https://en.wikipedia.org/wiki/Non-uniform_discrete_Fourier_transform
import numpy as np
import matplotlib.pyplot as plt
from pynufft import NUDFT

w,h = plt.rcParams['figure.figsize'] # figure size matplotlib to scale figure with subfigures

# %% Generate sample data
m = 1000 # dimension

# Some dummy signal frequency
f1 = 15

# uniform data
uniform_time_stamps = np.linspace(0, 20, m)  # time vlaues for uniform distributed data
uniform_current = np.sin(2*np.pi*uniform_time_stamps*f1) + np.random.randn(m)/10 # uniform distributed sample points p_{n} (a current in Ampere, with some pick-up noise at 15Hz)

# non-uniform data
time_stamps = np.linspace(0, 20, m) + np.random.randn(m)/10  # time vlaues for non-uniform distributed data
current = np.sin(2*np.pi*time_stamps*f1) + np.random.randn(m)/10 # non-uniform distributed sample points p_{n} (a current in Ampere, with some pick-up noise at 2Hz)

# %% First calculate rFFT for comparison
sf = 1/m # sample freq for normalization 1/sqrt(Hz)
n = uniform_time_stamps.size
t = uniform_time_stamps[-1]
d = t / n # sample spacing

# Calculate frequency range and concat to table
freqRange = np.fft.rfftfreq(n, d=d)
   
# Calculate real FFT  / normalised by the bandwith 1/sqrt(n) with norm="ortho" AND /np.sqrt(sf)
FFT = np.abs(np.fft.rfft(uniform_current, norm="ortho") / np.sqrt(sf))

# %% Try to calculate non-uniform discrete Fourier transform
# Initiate the NufftObj object
NufftObj = NUFFT()

# Select grid parameter
Nd = (1000,) # number of sample points
Kd = (2*1000,) # at least two times?
Jd = (6,)

x = time_stamps.reshape(1000,1)
y = current.reshape(1000,1)

NufftObj.plan(x, Nd, Kd, Jd)

# is this right?:
f = NufftObj.adjoint(y)  # calculate frequencies from NufftObj -> should yield something from 0 Hz to 25 Hz?
y_f = NufftObj.forward(y) # Calculate spectral density


# %% Plot
fig, ax = plt.subplots(1,3, figsize=(3*w,h), layout="constrained")

# Data
ax[0].set_title("Data")
ax[0].plot(time_stamps, current, "k-", label="data", alpha=0.3)
ax[0].set_xlabel("Time")
ax[0].set_ylabel("Curren /A")

# FFT
ax[1].set_title("FFT")
ax[1].plot(freqRange, FFT, "g-", label="Real part spectral density")
ax[1].set_xlabel("Frequency Hz")
ax[1].set_ylabel("Spectral density in A /Hz$^{1/2}$")

# NUFFT
ax[2].set_title("NUFFT")
ax[2].plot(f, y_f.real, "g-", label="Real part spectral density")
ax[2].set_xlabel("Frequency Hz")
ax[2].set_ylabel("Spectral density in A /Hz$^{1/2}$")

plt.show()