使用 Numpy 手动反转 FFT

Posted

技术标签:

【中文标题】使用 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 Nq = N - pb=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 结果也反映了复数共轭。您需要所有这些共轭来获得逆。

以上是关于使用 Numpy 手动反转 FFT的主要内容,如果未能解决你的问题,请参考以下文章

Numpy fft 冻结更长的样本

FFTW 从 numpy.fft 产生不同的结果

Python NumPy 将 FFT 转换为文件

Numpy FFT错误 - 带信封的梳状滤波器

Matlab 中的 FFT 和 numpy / scipy 给出不同的结果

无法正确生成 Numpy FFT