在Python中将FFT绘制为一组正弦波？

Question

我在演示中看到有人这样做，但我很难重现他能够做到的事情。这是他演讲中的一张幻灯片：

很酷。他使用 FFT 分解数据集，然后绘制 FFT 指定的适当正弦波。

因此，为了重现他所做的事情，我创建了一系列对应于 2 个正弦波组合的点：

import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

x = np.arange(0, 10, 0.01)
x2 = np.arange(0, 20, 0.02)
sin1 = np.sin(x)
sin2 = np.sin(x2)
x2 /= 2
sin3 = sin1 + sin2
plt.plot(x, sin3)
plt.show()

现在我想将这个波（或者更确切地说，这些点暗示的波）分解回原来的 2 个正弦波：

# goal: sin3 -> sin1, sin2
# sin3 
array([ 0.00000000e+00,  2.99985000e-02,  ... 3.68998236e-01])
# sin1 
array([ 0.        ,  0.00999983,  0.01999867,  ... -0.53560333])
# sin2 
array([ 0.        ,  0.01999867,  0.03998933, ... 0.90460157])

我首先导入并numpy获取：fftsin3

import numpy as np
fft3 = np.fft.fft(sin3)

好吧，据我所知，这就是了。现在我有一个包含复数的数组：

array([ 2.13316069e+02+0.00000000e+00j,  3.36520138e+02+4.05677438e+01j,...])

如果我天真地绘制它，我会看到：

plt.plot(fft3)
plt.show()

好吧，不知道该怎么办。

我想从这里获取看起来像 sin1 和 sin2 的数据集：

plt.plot(sin1)
plt.show()

plt.plot(sin2)
plt.show()

我了解数据集中复数的实部和虚部fft3，我只是不确定如何处理它们来从中导出数据集sin1。sin2

我知道这与编程关系不大，与数学关系更大，但有人可以在这里给我提示吗？

编辑：更新马克·斯奈德的答案：

使用马克的代码，我能够得到我所期望的结果，并最终得到了这个方法：

def decompose_fft(data: list, threshold: float = 0.0):
    fft3 = np.fft.fft(data)
    x = np.arange(0, 10, 10 / len(data))
    freqs = np.fft.fftfreq(len(x), .01)
    recomb = np.zeros((len(x),))
    for i in range(len(fft3)):
        if abs(fft3[i]) / len(x) > threshold:
            sinewave = (
                1 
                / len(x) 
                * (
                    fft3[i].real 
                    * np.cos(freqs[i] * 2 * np.pi * x) 
                    - fft3[i].imag 
                    * np.sin(freqs[i] * 2 * np.pi * x)))
            recomb += sinewave
            plt.plot(x, sinewave)
    plt.show()

    plt.plot(x, recomb, x, data)
    plt.show()

稍后我会让它返回重新组合的波浪列表，但现在我遇到了一个我不太明白的异常情况。首先，我这样称呼它，只需传入一个数据集。

decompose_fft(sin3, threshold=0.0)

但看起来不错，但我在“y=0.2有谁知道这可能是什么或导致它的原因是什么？”处得到了这条奇怪的线？

编辑：

上述问题已由Mark在评论中解答，谢谢！

Answer 1

离散傅里叶变换给出复指数的系数，当它们相加时，产生原始的离散信号。特别是，第 k 个傅立叶系数为您提供有关在给定数量的样本上具有 k 个周期的正弦曲线幅度的信息。

请注意，由于您的正弦在 1000 个样本中没有整数个周期，因此您实际上无法使用 FFT 检索原始正弦波。相反，您将得到许多不同正弦曲线的混合，包括 ~.4 的恒定分量。

您可以使用以下代码绘制各个正弦波分量并观察它们的总和是原始信号：

freqs = np.fft.fftfreq(len(x),.01)
threshold = 0.0
recomb = np.zeros((len(x),))
for i in range(len(fft3)):
    if abs(fft3[i])/(len(x)) > threshold:
        recomb += 1/(len(x))*(fft3[i].real*np.cos(freqs[i]*2*np.pi*x)-fft3[i].imag*np.sin(freqs[i]*2*np.pi*x))
        plt.plot(x,1/(len(x))*(fft3[i].real*np.cos(freqs[i]*2*np.pi*x)-fft3[i].imag*np.sin(freqs[i]*2*np.pi*x)))
plt.show()

plt.plot(x,recomb,x,sin3)
plt.show()

通过更改threshold，您还可以选择排除低功率的正弦曲线，并查看它如何影响最终的重建。

编辑：上面的代码中有一点陷阱，尽管它没有错误。它隐藏了真实信号 DFT 的固有对称性，并以真实幅度的一半绘制每个正弦曲线两次。此代码性能更高，并以正确的幅度绘制正弦曲线：

import cmath

freqs = np.fft.fftfreq(len(x),.01)
threshold = 0.0
recomb = np.zeros((len(x),))
middle = len(x)//2 + 1
for i in range(middle):
    if abs(fft3[i])/(len(x)) > threshold:
        if i == 0:
            coeff = 2
        else:
            coeff = 1
        sinusoid = 1/(len(x)*coeff/2)*(abs(fft3[i])*np.cos(freqs[i]*2*np.pi*x+cmath.phase(fft3[i])))
        recomb += sinusoid
        plt.plot(x,sinusoid)
plt.show()

plt.plot(x,recomb,x,sin3)
plt.show()

如果在一般情况下，您知道信号由某些正弦波子集组成，其频率可能与信号长度不正确对齐，则您可以通过零填充或扩展信号来识别频率。您可以在这里了解更多信息。如果信号完全是任意的，并且您只是对查看正弦波分量感兴趣，则无需这样做。