绘制和提取 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 结果呢？
• 这里我们处于一个非常简单的情况，我们有 N 个周期的正弦周期信号如果我有一个宽泛信号并且我想提取 f 处的相位我该怎么办？例如，对于示例中介绍的两种方法，我不确定是否可以提取精确的相位。对于phase_spectrum，在 f = 1 时我无法找到 pi/4。使用 fft 和 np.angle，为了提取良好的相位，我需要确保周期信号数是整数。

Answer 1

在回答之前，先做一个小注释：
删除该p[np.abs(Y) < 1] = 0行。大多数频谱的幅度都低于 1，这就是为什么这条线的频谱看起来大多像一条位于 0 处的平坦线。

现在回答：
phase_spectrum做了三件事与你不同：

• 它应用相位展开。
  • 如果您想在代码中应用相位展开，只需执行np.unwrap(np.angle(Y)).
  • 如果您希望 matplotlib 绘制频谱而不展开，请angle_spectrum改用。
• 它在计算频谱之前对数据应用窗口函数。
  • 我知道你通过了一个window=None，但是，出于某种原因，matplotlib 认为这window=None意味着“请使用汉宁窗口”（请参阅​​文档）。
  • 如果您不希望 matplotlib 应用窗口，一种解决方案是传递window=lambda x: x.
    • 该文档实际上建议通过window=matplotlib.mlab.window_none，但它的来源只是一个def window_none(x): return x.
• 它计算频谱的单侧版本。文档说，只要输入是真实的而不是复杂的，这就是默认值。
  • 要获得正常的双面版本，请传递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)。