时间戳

首页 > 解决方案 > Python rfft 算法

python-3.x - Python rfft 算法

问题描述

Numpy 具有rfft计算 FFT 或 Real 输入的功能。我想为此编写自己的代码。这是我的 DFT 和 FFT [ ref ] 的实际代码。


def DFT_slow(x):
    """Compute the discrete Fourier Transform of the 1D array x"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)

def FFT(x):
    """A recursive implementation of the 1D Cooley-Tukey FFT"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]

    if N % 2 > 0:
        raise ValueError("size of x must be a power of 2")
    elif N <= 32:  # this cutoff should be optimized
        return DFT_slow(x)
    else:
        X_even = FFT(x[::2])
        X_odd = FFT(x[1::2])
        factor = np.exp(-2j * np.pi * np.arange(N) / N)
        return np.concatenate([X_even + factor[:N / 2] * X_odd,
                               X_even + factor[N / 2:] * X_odd])

我对傅里叶变换不太了解，知道如何rfft从这些代码中计算？

标签： python-3.x

解决方案

rfft只是整个输入的正频率分量，fft因为fft实值信号的频率是共轭对称的。

def rDFT_slow(x):
    """Compute the discrete Fourier Transform of the 1D array of real data x"""
    x = np.real(np.asarray(x, dtype=float))
    N = x.shape[0]
    n = np.arange(N)
    k = np.arange(N//2+1).reshape((N//2+1,1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)

当您只计算一半的系数时，实现递归版本有一些更棘手的因素，但缓慢rDFT应该让您开始。

推荐阅读