使用 Numpy 手动反转 FFT

Posted 2023-03-12

【中文标题】使用 Numpy 手动反转 FFT

【英文标题】：Manually inverting FFT using Numpy

【发布时间】：2016-03-02 14:52:55

【问题描述】：

我有一个小脚本用于计算方波的傅立叶变换，当我使用 numpy.fft.ifft() 反转 fft 时，它运行良好并正确返回方波。但是，在将谐波乘以我从numpy.fft.fft() 获得的相应系数之后，我无法通过手动相加谐波来反转变换。下面是我的脚本，我相信你会看到我的意图。

from numpy import zeros, concatenate, sin, pi, linspace
from numpy.fft import fft, fftfreq, ifft
import numpy as np
import matplotlib.pyplot as plt

N = 1024 # samples
T = 1 # period
dt = T/N # sampling period
fs = 1/dt # sampling frequency
t = linspace(0, T, N) # time points
functime = .... # square wave

funcfft = fft(functime) # fft
fftcoeffs = np.abs(funcfft)/N # coefficients, divide by N to get actual coeff.s(I believe?)
freqs = fftfreq(N, dt) # frequencies

plt.plot(freqs, fftcoeffs) # gives me reasonable output
plt.show()

FF = ifft(funcfft)
plt.plot(t, FF) # plots exactly the same function as functime defined above
plt.show()

到目前为止一切都很好。现在我的问题是，如果我在上面的脚本之后运行下面的脚本，我不应该收敛到原始函数吗？：

FFF = zeros(N)
for i in range(300):
    FFF += fftcoeffs[i]*sin(2*pi*freqs[i]*t)
plt.plot(t, FFF)
plt.show()

假设range(300) 足以收敛。现在当我这样做时，FFF 与我原来的功能不同。我认为如果我将各个频率的谐波乘以它们对应的系数，我认为这些系数存储在 fftcoeffs 中，然后我会收敛到原始函数。我做错了什么？

更新：根据 DanielSank 的建议，我更新了我的 for 循环如下，不幸的是它没有给我想要的结果：

freqs2 = np.abs(freqs)
freqs2 = np.sort(freqs2)
for k in range(300):
    FFF += fftcoeffs[k]*exp(2j*pi*freqs2[k]*t/N)

我不确定我是否在此处执行“按绝对值排序 fftfreq”部分。

【问题讨论】：

您没有正确排序频率 :D 另外，一般来说，当您在任何在线论坛上发布代码时，请包含一个自包含的工作示例 (SCWE)。这被认为是标准做法。不这样做会非常非常频繁地使本来会给您答案的人忽略该帖子。

【参考方案1】：

术语

此问题中的任何内容均不针对 fast Fourier transform (FFT)。 FFT 是一种计算discrete Fourier transform (DFT) 的特殊算法，所以我会说“DFT”而不是“FFT”。

在这个答案中，m 表示离散时间索引，k 表示离散频率索引。

什么是 DFT？

这里有几个问题，所有这些问题都源于对 DFT 工作原理的误解。 取自numpy.fft 模块文档字符串，numpy 将离散傅里叶变换定义为

A_k = \sum_m=0^n-1 a_m \exp[-2 \pi i (m k / n)]

LaTeX 表示法表示离散傅里叶变换是复指数的线性组合 exp[2 pi i m k / n] 其中n 是点总数，m 是离散时间索引。 在您的符号中，这将是 exp[2 pi i m k / N]，因为您使用 N 来表示总点数。

使用exp，而不是sine

请注意，DFT 使用指数；这些是不是 sine 函数。 如果要从离散傅立叶变换系数构建时域信号，则必须使用 DFT 本身使用的相同函数！ 因此，你可以试试这个：

FFF = zeros(N)
for i in range(300):
    FFF += fftcoeffs[i]*np.exp(2*pi*i*freqs[i]*t)
plt.plot(t, FFF)
plt.show()

但是，这也会以一种可能让您困惑的方式失败。

别名

最后一块拼图与称为混叠的效果有关。 假设您对信号exp[2 pi i (N + p) m / N] 进行 DFT。 如果你计算它，你会发现所有的A_k 都是零，除了A_p。 事实上，如果你参加了exp[2 pi i p m / N] 的 DFT，你会得到相同的。 您可以看到频率大于N 的任何复指数都显示为好像它是较低的频率。 特别是，频率为q + b N 的任何复指数，其中b 是任何整数，看起来它的频率为q。

现在假设我们有一个时域信号cos(2 pi p m / N)。 这等于

(1/2)[ (exp(2 pi i p m / N) + exp(-2 pi i p m / N) ].

负频率对 DFT 产生了有趣的影响。 频率-p可以写成(N-p) + N。 这具有q + b N 与q = N - p 和b=1 的形式。 所以，那个负频率-p 最终看起来像N - p！

numpy 函数fftfreq 知道这一点。 看看fftfreq 的输出，您会看到它从零开始，运行到采样频率的一半（称为奈奎斯特频率），然后变为负数！ 这是为了帮助您处理我们刚刚描述的锯齿效应。

所有这一切的结果是，如果您想通过对最低频率傅立叶分量求和来逼近一个函数，您确实不想从fftfreq 中获取最低的几个元素。 相反，您希望按绝对值对 fftfreq 进行排序，然后将复指数与这些频率相加。

还可以查看np.fft.hfft。 该函数旨在处理实值函数以及与之相关的别名问题。

代码

由于这是一个很难口头讨论的问题，这里有一个脚本可以完全满足您的需求。 请注意，我将 cmets 置于 那些 cmets 描述的代码块之后。 确保您已安装 matplotlib（在您的虚拟环境中......您使用虚拟环境，对吗？）。 如果您有任何问题，请发表评论。

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt


pi = np.pi


def square_function(N, square_width):
    """Generate a square signal.

    Args:
        N (int): Total number of points in the signal.
        square_width (int): Number of "high" points.

    Returns (ndarray):
        A square signal which looks like this:

              |____________________
              |<-- square_width -->
              |                    ______________
              |
              |^                   ^            ^
        index |0             square_width      N-1

        In other words, the output has [0:N]=1 and [N:]=0.
    """
    signal = np.zeros(N)
    signal[0:square_width] = 1
    return signal


def check_num_coefficients_ok(N, num_coefficients):
    """Make sure we're not trying to add more coefficients than we have."""
    limit = None
    if N % 2 == 0 and num_coefficients > N // 2:
        limit = N/2
    elif N % 2 == 1 and num_coefficients > (N - 1)/2:
        limit = (N - 1)/2
    if limit is not None:
        raise ValueError(
            "num_coefficients is  but should not be larger than ".format(
                num_coefficients, limit))


def test(N, square_width, num_coefficients):
    """Test partial (i.e. filtered) Fourier reconstruction of a square signal.

    Args:
        N (int): Number of time (and frequency) points. We support both even
            and odd N.
        square_width (int): Number of "high" points in the time domain signal.
            This number must be less than or equal to N.
        num_coefficients (int): Number of frequencies, in addition to the dc
            term, to use in Fourier reconstruction. This is the number of
            positive frequencies _and_ the number of negative frequencies.
            Therefore, if N is odd, this number cannot be larger than
            (N - 1)/2, and if N is even this number cannot be larger than
            N/2.
    """
    if square_width > N:
        raise ValueError("square_width cannot be larger than N")
    check_num_coefficients_ok(N, num_coefficients)

    time_axis = np.linspace(0, N-1, N)

    signal = square_function(N, square_width)
    ft = np.fft.fft(signal)

    reconstructed_signal = np.zeros(N)
    reconstructed_signal += ft[0] * np.ones(N)
    # Adding the dc term explicitly makes the looping easier in the next step.

    for k in range(num_coefficients):
        k += 1  # Bump by one since we already took care of the dc term.
        if k == N-k:
            reconstructed_signal += ft[k] * np.exp(
                1.0j*2 * pi * (k) * time_axis / N)
        # This catches the case where N is even and ensures we don't double-
        # count the frequency k=N/2.

        else:
            reconstructed_signal += ft[k] * np.exp(
                1.0j*2 * pi * (k) * time_axis / N)
            reconstructed_signal += ft[N-k] * np.exp(
                1.0j*2 * pi * (N-k) * time_axis / N)
        # In this case we're just adding a frequency component and it's
        # "partner" at minus the frequency

    reconstructed_signal = reconstructed_signal / N
    # Normalize by the number of points in the signal. numpy's discete Fourier
    # transform convention puts the (1/N) normalization factor in the inverse
    # transform, so we have to do it here.

    plt.plot(time_axis, signal,
             'b.', markersize=20,
             label='original')
    plt.plot(time_axis, reconstructed_signal.real,
             'r-', linewidth=3,
             label='reconstructed')
    # The imaginary part is zero anyway. We take the real part to
    # avoid matplotlib warnings.

    plt.grid()
    plt.legend(loc='upper right')

【讨论】：

感谢您的回复。首先，我真的不知道如何定义i 又名-1**0.5。显然-1**0.5 没有定义它。但无论如何，我做了freqs2 = np.abs(freqs)，然后freqs2 = np.sort(freqs2)，然后按照你的建议更新了我的for循环，使用exp而不是sin。我使用k 作为for 循环变量，而不是i，这样它就不会与想象中的i 混淆。我猜你忘了在exp() 中包含/N，我尝试了使用和不使用它，但我在所有尝试中都得到了更糟糕的结果。我怀疑我是否完全理解您的建议。

@Deniz 很难理解代码的书面描述。您可以编辑原始帖子以显示您尝试过的内容吗？

我刚刚看到您对这篇文章的更新并试用了您的脚本。它完美无缺。就我而言，一个周期的正方形（实际上是矩形）信号不像对称的 +a/- 一种正方形信号，所以我必须根据我的情况进行调整，但您的代码仍然可以帮助我理解程序和 DFT 理论。非常感谢您帮助我丹尼尔。我选择这个答案作为正确答案。

@Deniz 很高兴能帮上忙。我发现关于信号处理的文献不能令人满意，所以我就这些主题写了自己的一套笔记。如果您愿意，我很乐意给您一份副本。如果您查看我的物理堆栈交换用户资料并点击我的专业网站的链接，您可以联系我。

@Deniz 你不能@我，因为你在评论我的帖子。用户会自动在自己的帖子中收到有关 cmets 的通知。另外，我想提请您注意，我为您提供了一个完整的工作代码示例，它对您有所帮助。请在您的问题中也这样做。举个例子，问题更容易理解；描述代码几乎是无用的。此外，编写一个最小的示例可以让您在 99% 的时间里解决自己的问题。

【参考方案2】：

最大的问题是您只使用了 FFT 结果的幅度 (np.abs)，因此丢弃了所有相位信息。

如果您保留 FFT 的复杂相位结果，并在正弦波重新合成中使用该相位信息，您手动添加的谐波将更接近原始结果。

可能的提示：您可能必须在 FFT 之前对波形进行 fftshift 以获得有用的相位结果，具体取决于方波的频率与 FFT 孔径宽度的关系。

【讨论】：

感谢您的回答。使用我上面给出的相同 for 循环，唯一的例外是使用 funcfft 而不是 fftcoeffs，我得到了更糟糕的结果。使用我当前的代码，它给我的完全错误，但它是一个正方形，当我不使用 fftcoeffs 时，它给我的是指数而不是正方形。

为了从复指数中得到一个真实的结果，所有的频率都必须用复共轭来镜像。请注意，完整的复数 FFT 结果也反映了复数共轭。您需要所有这些共轭来获得逆。