tarabyte
tarabyte

Reputation: 19172

How to compute nfft

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)

Upvotes: 3

Views: 3199

Answers (1)

SleuthEye
SleuthEye

Reputation: 14577

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)

Upvotes: 4

Related Questions