绘制和提取 fft 阶段答案

【问题标题】：Plotting and extracting fft phase绘制和提取 fft 阶段
【发布时间】：2023-03-30 03:02:01
【问题描述】：

这是一个将 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 会产生很好的结果，但是我们如何解释幅度谱的结果呢？
这是一个非常简单的情况，we have sine periodic signal with N periods 如果我有一个宽带信号并且我想在 f 处提取相位，我该怎么做？例如，这里使用示例中介绍的两种方法，我不确定我是否可以提取精确的相位。使用 phase_spectrum，在 f = 1 时，我无法找到 pi/4。然后使用 fft 和 np.angle，为了提取好的相位，我需要确保周期的信号数是整数。

【问题讨论】：

标签： python numpy fft phase

【解决方案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 的相位的近似值。您还可以对频谱进行 sinc 插值，以获得所需频率的值。

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

【讨论】：

【解决方案2】：

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

因此，计算窗口中间的信号在生成过程中的相位，并使用它来代替开始时的相位。

【讨论】：