将 FFT 绘制为 python 中的一组正弦波？

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 并获取 sin3 的 fft:

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 数据集来自它。

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

编辑：更新 Mark Snyder 的回答：

使用 Mark 的代码，我能够得到我所期望的，并以这种方法结束：

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

稍后我会将其设为 return 重组的波浪列表，但现在我遇到了一个我不太明白的异常现象。首先我这样称呼它，简单地传入一个数据集。

decompose_fft(sin3, threshold=0.0)

但看起来不错，但我在 y=0.2 看到了这条奇怪的线有谁知道这可能是什么或是什么原因造成的？

编辑：

以上问题已在评论区马克回答，谢谢！

Answer 1

离散傅里叶变换存在一些问题，这些问题在使用连续变换时并不是很明显。一方面，输入的周期应与数据的范围相匹配，因此如果使用：

会容易得多

x = np.linspace(0, 4*np.pi, 200)

然后您可以按照您最初的想法：

sin1 = np.sin(x)
sin2 = np.sin(2*x)
sin3 = sin1 + sin2
fft3 = np.fft.fft(sin3)

由于 FFT sin 直接进入虚部，您可以尝试只绘制虚部：

plt.plot(fft3.imag)
plt.show()

您应该看到以 x=2 和 x=4 为中心的峰值，对应于原始正弦分量，其频率为“每个信号 2”（sin(x) 从 0 到4 pi) 和“每个信号 4”（从 0 到 4 pi 的 sin(2x)）。

要绘制所有单独的组件，您可以使用：

for i in range(1,100):
  plt.plot(x, fft3.imag[i] * np.sin(i*x)/100)
plt.show()

Answer 2

离散傅立叶变换为您提供了复指数的系数，当这些系数加在一起时，会产生原始的离散信号。特别是，第 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 的固有对称性，并以其真实幅度的一半绘制每个正弦曲线两次。此代码性能更高，并以正确的幅度绘制正弦曲线：

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

如果在一般情况下您知道信号由正弦波的某些子集组成，其频率可能与信号的长度不正确对齐，您可以通过补零来识别频率或延长你的信号。您可以了解更多相关信息 here。如果信号完全是任意的，而您只对查看分量正弦曲线感兴趣，则没有必要。

将 FFT 绘制为 python 中的一组正弦波？

Plot FFT as a set of sine waves in python?

python

numpy

fft

sine-wave