Построение и извлечение FFT фазы

Question

Вот код, который сравнивает построение БПФ фазы с двумя различными методами:

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

phase = np.pi / 4
f = 1
fs = f*20
dur=10
t = np.linspace(0, dur, num=fs*dur, endpoint=False)
y = np.cos(2 * np.pi * t + phase)
Y = scipy.fftpack.fftshift(scipy.fftpack.fft(y))
f = scipy.fftpack.fftshift(scipy.fftpack.fftfreq(len(t)))

p = np.angle(Y)
p[np.abs(Y) < 1] = 0

fig, ax = plt.subplots(2, 1)
ax[0].plot(t, y)
ax[1].plot(f*fs, p, label='from fft')
ax[1].phase_spectrum(y, fs, window=None, label='from phase_spectrum')
plt.legend()
plt.show()

вот результат:

Вот результат, когда номер сигнала периода не является целым числом:

У меня есть несколько вопросов:

• почему фазовый график с phase_spectrum или fft и угол так разные? Использование fft, а затем np.angle дает хорошие результаты, но как мы можем объяснить результат magnitude_spectrum?
• Здесь мы находимся в очень простом случае, когда у нас есть сигнал sine periodi c с N периодами Если у меня есть сигнал широкого диапазона, и я хочу извлечь фазу при f, как я могу это сделать? Например, здесь с обоими методами, представленными в примере, я не уверен, что могу извлечь точную фазу. С phase_spectrum, при f = 1 я не могу найти обратно pi / 4. И с помощью fft, а затем np.angle, чтобы извлечь хорошую фазу, мне нужно убедиться, что номер сигнала периода является целым числом.

Answer 1

Перед ответом просто небольшая заметка:
Удалите строку p[np.abs(Y) < 1] = 0. Большая часть вашего спектра имеет величину ниже 1, поэтому с этой линией ваш спектр выглядит в основном как плоская линия в нуле.

Теперь к ответу:
phase_spectrum делает три вещи, отличные от вас:

• Применяется разворачивание фазы.
  • Если вы хотите применить развертывание фазы в вашем коде, просто выполните np.unwrap(np.angle(Y)).
  • Если вы хотите, чтобы matplotlib отображал спектр без развертывания, используйте angle_spectrum вместо этого.
• Применяет оконную функцию к данным перед вычислением спектра.
  • Я знаю, что вы передали window=None, но по какой-то причине matplotlib решил, что window=None означает "используйте окно Hanning, пожалуйста" (см. документы ).
  • Если вы не хотите, чтобы matplotlib применял окно, одно из решений - передать window=lambda x: x.
    • В документах фактически предлагается передать window=matplotlib.mlab.window_none, но source для него всего лишь def window_none(x): return x.
• Он рассчитывает одностороннюю версию вашего спектра. docs говорят, что это значение по умолчанию, когда ввод действительный, а не сложный
  • Чтобы получить обычную двустороннюю версию, передайте sides='twosided' на phase_spectrum вызов.

---

Теперь о получении фазы с частотой f:

Для этого вы должны использовать фазу без развертывания .

Вы правы, что можете t непосредственно извлекать фазу однотонального сигнала, если у вас нет целого числа циклов. Это связано с тем, что частота сингла не попадает точно в верхнюю часть частотного бина в БПФ. Вы можете получить приближение с фазой ближайшего бина, хотя. Вы можете также выполнить интерполяцию спектра sin c, чтобы получить его значение на нужной частоте.

Если вы заботитесь только о фазе single частота f, тогда вы не должны использовать БПФ вообще. БПФ вычисляет фазу и амплитуду на всех частотах. Если вам нужна только одна частота, просто наберите Y_at_f = y @ np.exp(2j * np.pi * f * t) и получите эту фазу на np.angle(Y_at_f).

Answer 2

Вы можете извлечь фазу, на которую ссылается центр вашего окна данных, выполнив fftshift (круговое вращение на N / 2) перед FFT. Это связано с тем, что после fftshift функция atan2 () всегда связана с отношением нечетности к четности данных вокруг их центра (разложенных на нечетную функцию плюс четную функцию).

Итак, рассчитайте фазу сигнала в середине вашего окна во время его генерации, и используйте его вместо фазы в начале.