python - 如何计算nfft

我试图了解如何使用Jake Vanderplas' nfft module的nfft方法。不幸的是，该示例不是很好的说明，因为我尝试仅基于样本的输入列表（[(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)

我试图在一个示例中计算f_k，其中样本之间的间隔约为10 ms，抖动为1或2 ms。

实施文档：

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).

我卡住的地方：

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)

最佳答案

该示例定义了变量f_k，该变量作为nfft的f_hat参数传递。
根据定义

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

给定的f_hat表示指定采样时刻的时域信号。在您的情况下，这仅对应于sample_values。

x的另一个参数nfft是这些样本的实际时刻。您还需要分别提供：

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)

关于python - 如何计算nfft，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/57435660/