如何在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()

现在我想将这个波（或者说点所暗示的波）分解回原始的两个正弦波：

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

稍后我会让它返回重新组合的波形列表，但现在我遇到了一个我不太理解的异常情况。首先，我这样调用它，只是简单地传入一个数据集。

decompose_fft(sin3, threshold=0.0)

看起来很好，但我在y=0.2处出现了奇怪的线条，有人知道这是什么或是什么导致的吗？

编辑：

上面的问题已经在评论中被Mark回答了，谢谢！

- MetaStack

---

我的错，我认为你需要在复平面上绘制它们：https://www.tutorialspoint.com/How-to-Plot-Complex-Numbers-in-Python - Felipe Gutierrez

@FelipeGutierrez 我无法从中获得太多信息，但是无论如何，在复平面上简单绘制它并不能给我数据集 sin1 或 sin2（或告诉我它们的频率，以便我可以推导出它们），对吧？这是我尝试的图形，看起来很酷，但并没有告诉我太多信息。除非我当然不知道如何解释它：plt.scatter(fft3.real, fft3.imag, color='red') - MetaStack

我更新了我的答案，并附上了一些代码，可能会对你有所帮助。 - Mark Snyder

顺便提一下，np.fft.fftfreq() 的第二个参数应该与 np.arange() 中的第三个参数相同。 - Mark Snyder

添加参数ith: int = 0并修改代码，使得i > ith:，这个奇怪的行就消失了。 - Anibal Yeh

Answer 1

离散傅里叶变换可以给出复指数系数，当它们相加时，可以产生原始离散信号。特别地，第k个傅里叶系数可以给出在给定样本数内，具有k个周期的正弦波的幅度信息。

请注意，由于您的正弦波在1000个样本中没有整数周期，因此您实际上无法使用FFT检索原始正弦波。相反，您将得到许多不同正弦波的混合，包括约为0.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()

如果在一般情况下，您知道信号由某些正弦子集组成，其频率可能与信号长度不对齐，您可以通过零填充或扩展信号来识别频率。您可以在这里了解更多信息。如果信号完全是任意的，而您只想查看组件正弦波，则没有必要这样做。

- Mark Snyder

---

这太棒了！只有一个异常我不理解，我已经在问题中添加了一个编辑来解释我看到的内容。 - MetaStack

@LegitStack 这条线路并不是异常值，它是一个频率为0的正弦曲线 - 也就是说它是个常数！实际上，你应该发现这是数据的平均值。请记住，每个频率不等于0的正弦曲线都有一个平均值为0。如果没有直线来保持数据的平均值，你将无法对没有0平均值的数据进行傅里叶变换。 - Mark Snyder

Answer 2

离散傅里叶变换存在一些问题，这些问题并不像操作其连续对应物时那么显而易见。首先，您的输入数据的周期性应该与数据范围相匹配，因此如果使用以下内容，则会更加容易：

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次”（从0到4 pi的sin（x））和“每个信号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()

- Miłosz Wieczór

---

谢谢你，这真的很棒！我想知道你是否能帮我处理一个更一般的情况。我能够像您建议的那样完美地工作（实际上我不得不做 np.sin(i*x/2)）。但是当周期性与我的数据范围不匹配时怎么办？那时我该怎么办，因为我本质上无法访问 x？即使尝试对 sin1 和 sin2 进行不同的变异，例如 sin2 = np.sin(1.9*x)，也会破坏我能够推导出的图形，使它们与原始图形不匹配。你能帮我理解和处理更一般的情况吗？ - MetaStack

你大致上想要做的是将 sin(1.9x) 表达为一系列项 sin(x), sin(2x), sin(3x) 等的和（同时也包括相应的余弦项，因为余弦与正弦正交）。显然，这不能用单个项来完成，这就是为什么你会看到多个贡献。在我看来，讲义幻灯片中的“清晰”示例大多是针对明确定义的组件绘制的正弦和的总和，而不是基于实际信号的 DFT 分析。;) - Miłosz Wieczór

谢谢您的回答，这非常有教育意义。但不幸的是，它并不能解决一般情况。我只是为了创建 sin3（一个干净地表示一系列正弦波组合的数据点）而创建了 sin1 和 sin2。在这个演示作为跳板的真实世界应用中，我没有一开始就拥有 sin1 或 sin2，我只有类似于 sin3 的数据集，并且需要从仅给定该数据集的情况下推导出一组构成正弦波的子集：这是我试图解决的一般情况。 - MetaStack

问题在于DFT的“基本”频率单位是整个信号。因此，如果您可以恰好将一个周期放入正在转换的信号中，则这将对应于一个频率，依此类推。您想要的只有对于连续变换才有效果，因为它使用连续的频谱，并且任何“纯”正弦/余弦都将产生尖峰。从实际角度来看，我可以肯定，您的信号中拥有更多完整周期，您将拥有更好定义的单正弦组件-请尝试比较例如从0到4pi和从0到40 pi的范围。 - Miłosz Wieczór