python - 绘制和提取 fft 相位

Question

下面是比较 fft 相位图与 2 种不同方法的代码:

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结果？
• 这是一个非常简单的情况，其中 we have sine periodic signal with N periods如果我有一个宽泛信号并且我想提取 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意思是“请使用汉宁窗”(参见docs)。
  • 如果您不希望 matplotlib 应用窗口，一种解决方案是传递 window=lambda x: x 。
    • docs实际上建议通过window=matplotlib.mlab.window_none ，但是source因为它只是一个 def window_none(x): return x .
• 它计算频谱的单侧版本。 docs只要输入是真实的而不是复杂的，就说这是默认值。
  • 要获取普通双面版本，请传递 sides='twosided'到phase_spectrum打电话。

<小时/>

现在，关于获取频率 f 的相位:

为此，您必须使用阶段而不展开。

你是对的，如果没有整数个周期，则无法直接提取单音信号的相位。这是因为信号的频率并不正好落在 FFT 中任何频率仓的顶部。不过，您可以通过最近的 bin 的相位获得近似值。您还可以对频谱进行正弦插值，以获得所需频率处的值。

如果您只关心单个频率的相位f ，那么你根本不应该使用 FFT。 FFT 计算所有频率的相位和幅度。如果你只关心单个频率，就做 Y_at_f = y @ np.exp(2j * np.pi * f * t)并通过 np.angle(Y_at_f) 获得该阶段.