График БПФ как набор синусоидальных волн в python?

Question

Я видел, как кто-то делал это в презентации, но мне трудно воспроизвести то, что он смог сделать. Вот слайд из его презентации:

Довольно круто. Он разложил набор данных с использованием БПФ, а затем нанес на график соответствующие синусоидальные волны, указанные в БПФ.

Поэтому, чтобы воссоздать то, что он сделал, я создал серию точек, которые соответствуют комбинации из двух синусоид:

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

Хорошо, не уверен, что с этим делать.

Я хочу перейти отсюда к наборам данных, которые выглядят как 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()

позже я ' заставлю его вернуть рекомбинированный список волн, но сейчас я получаю аномальное ie, я не совсем понимаю. Прежде всего, я называю это так, просто передавая набор данных.

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)

Поскольку в БПФ sin переходит непосредственно в мнимый компонент, вы можете попробовать построить только мнимую часть:

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

То, что вы должны увидеть, будет иметь пики с центром в x=2 и x=4 которые соответствуют исходным синусоидальным компонентам, которые имели частоты «2 на сигнал» (sin (x) от 0 до 4 пи) и «4 на сигнал» (sin (2x) от 0 до 4 пи).

Чтобы построить все отдельные компоненты, вы можете go с помощью:

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

Answer 2

Дискретное преобразование Фурье дает вам коэффициенты комплексных экспонент, которые при суммировании дают исходный дискретный сигнал. В частности, k-й коэффициент Фурье дает вам информацию об амплитуде синусоиды, которая имеет k циклов в течение заданного числа выборок.

Обратите внимание, что поскольку ваши синусы не имеют целых чисел циклов в 1000 выборок, вы на самом деле не сможете получить свои первоначальные синусоиды с использованием БПФ. Вместо этого вы получите смесь из множества различных синусоид, включая постоянный компонент ~ .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, вы также можете выбрать исключение синусоид низкой мощности и посмотреть, как это повлияет на окончательную реконструкцию.

РЕДАКТИРОВАТЬ: в приведенном выше коде есть небольшая ловушка , хотя это не так. Он скрывает внутреннюю симметрию ДПФ для реальных сигналов и наносит на график каждую из синусоид дважды по половине их истинной амплитуды. Этот код более производительный и отображает синусоиды с их правильной амплитудой:

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

Если в общем случае вы знаете, что сигнал состоит из некоторого подмножества синусоид с частотами, которые могут не совпадать с По длине сигнала вы можете определить частоты с помощью дополнения нулями или расширения вашего сигнала. Вы можете узнать больше об этом здесь . Если сигналы абсолютно произвольны, и вам просто интересно посмотреть на составляющие синусоиды, в этом нет необходимости.