8

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

通过 FFT 进行正弦波分解

很酷。他使用 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并获取以下fft内容sin3：

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

fft 天真地绘制

好的，不知道该怎么办。

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

plt.plot(sin1)
plt.show()

sin1 数据绘制

plt.plot(sin2)
plt.show()

sin2 数据绘制

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

我知道这与编程无关，而与数学有关，但有人可以在这里给我提示吗？

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

使用 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 回答，谢谢！

4

2 回答 2

4

离散傅里叶变换为您提供复指数的系数，当它们相加时，会产生原始的离散信号。特别是，第 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()

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

于 2020-01-14T00:47:33.787 回答
3

离散傅里叶变换存在一些问题，这些问题在使用其连续对应物时并不明显。一方面，您的输入的周期性应该与您的数据范围相匹配，因此如果您使用它会容易得多：

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)

由于在 FFTsin中直接进入虚部，您可以尝试仅绘制虚部：

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

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

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

for i in range(1,100):
  plt.plot(x, fft3.imag[i] * np.sin(i*x)/100)
plt.show()
于 2020-01-14T00:43:33.470 回答