绘制和提取FFT相位

Question

这是一个比较两种不同方法的fft相位绘图的代码：

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，为了提取正确的相位，我需要确保信号周期数是整数。

- rem

Answer 1

在回答问题之前，先简要说明一下：
删除 p[np.abs(Y) < 1] = 0 这一行。大部分频谱的幅值都在1以下，这也是为什么加了这一行后你的频谱看起来像是一个基本平坦的直线。

现在开始回答：
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。
• 它计算的是单面频谱。当输入为实数而非复数时，默认情况下文档说这是默认值。
  • 要获得正常的双面版本，请在 phase_spectrum 调用中传递 sides='twosided'。

---

现在，关于如何获取频率为 f 的相位：

为此，你必须使用未经解缠的相位。

如果信号的频率不恰好落在FFT中的任何频率点上，则无法直接提取单音信号的相位，这是正确的。你可以获得最近频率点的相位值作为近似值。你还可以通过sinc插值来获取你想要的频率处的频谱值。

如果你只关心单一频率f的相位，那么就完全不需要使用FFT。FFT计算所有频率的相位和幅度。如果您只关心单个频率，请使用Y_at_f = y @ np.exp(2j * np.pi * f * t)来计算该频率下的数值，并用np.angle(Y_at_f)来获得相位值。

- Maxpxt

Answer 2

在进行FFT之前，您可以通过执行fftshift（N/2的循环旋转）来提取相位与数据窗口中心相关。这是因为，在fftshift之后，atan2（）始终与其中心周围数据的奇偶比有关（分解为奇函数加上偶函数）。

因此，在生成窗口期间计算信号中间的相位，并使用该相位代替开头的相位。

- hotpaw2