How to compute nfft

2.8k Views Asked by At

I'm trying to understand how to use the nfft method of Jake Vanderplas' nfft module. The example unfortunately isn't very illustrative as I try to parametrize everything based on just an input list of samples ([(time0, signal0), (time1, signal1), ...]):

import numpy as np
from nfft import nfft

# define evaluation points
x = -0.5 + np.random.rand(1000)

# define Fourier coefficients
N = 10000
k = - N // 2 + np.arange(N)
f_k = np.random.randn(N)

# non-equispaced fast Fourier transform
f = nfft(x, f_k)

I'm trying to compute f_k in an example where the samples are about 10 ms apart with 1 or 2 ms jitter in that interval.

The implementation documentation:

def nfft(x, f_hat, sigma=3, tol=1E-8, m=None, kernel='gaussian',
         use_fft=True, truncated=True):
    """Compute the non-equispaced fast Fourier transform

    f_j = \sum_{-N/2 \le k < N/2} \hat{f}_k \exp(-2 \pi i k x_j)

    Parameters
    ----------
    x : array_like, shape=(M,)
        The locations of the data points. Each value in x should lie
        in the range [-1/2, 1/2).
    f_hat : array_like, shape=(N,)
        The amplitudes at each wave number k = range(-N/2, N/2).

Where I'm stuck:

import numpy as np
from nfft import nfft

def compute_nfft(sample_instants, sample_values):
    """
    :param sample_instants: `numpy.ndarray` of sample times in milliseconds
    :param sample_values: `numpy.ndarray` of samples values
    :return: Horizontal and vertical plot components as `numpy.ndarray`s
    """
    N = len(sample_instants)
    T = sample_instants[-1] - sample_instants[0]
    x = np.linspace(0.0, 1.0 / (2.0 * T), N // 2)
    y = 2.0 / N * np.abs(y[0:N // 2])
    y = nfft(x, y)
    return (x, y)
1

There are 1 best solutions below

2
SleuthEye On

The example defines a variable f_k which is passed as nfft's f_hat argument. According to the definition

f_j = \sum_{-N/2 \le k < N/2} \hat{f}_k \exp(-2 \pi i k x_j)

given, f_hat represents the time-domain signal at the specified sampling instants. In your case this simply corresponds to sample_values.

The other argument x of nfft are the actual time instants of those samples. You'd need to also provide those separately:

def compute_nfft(sample_instants, sample_values):
    N = len(sample_instants)
    T = sample_instants[-1] - sample_instants[0]
    x = np.linspace(0.0, 1.0 / (2.0 * T), N // 2)
    y = nfft(sample_instants, sample_values)
    y = 2.0 / N * np.abs(y[0:N // 2])
    return (x, y)