python - 在 x 轴上移动频谱图

Question

我的信号从-1 秒开始。

但是，在绘制频谱图时，第一个 bin 边缘从 0 开始（中点在 0.25） 如何更改它以便在 x 轴上绘制时准确表示我的数据？

使用xextent=(time[0] + 0.125, time[-1])似乎可以解决这个问题。但是，我不确定是什么变量决定了 bin 宽度，因此担心这可能会随着不同采样率、点数等的其他数据集而改变。

示例代码：

import numpy as np
import matplotlib.pyplot as plt
import signal

time = np.linspace(-1, 16, 65536)
signal = np.sin(2 * np.pi * time)

fig = plt.figure()

ax_top = fig.add_subplot(211)
ax_spec = fig.add_subplot(212)

ax_top.plot(time, signal)
ax_top.set_xlim(-1, 10)
ax_spec.set_xlim(-1, 10)

Pxx, freqs, bins, cax = ax_spec.specgram(signal, NFFT=2048, Fs=4096, noverlap=2048 / 2, mode='magnitude', pad_to=2048 * 16)

ax_spec.set_xticks(np.arange(time[0], time[-1], 1))
ax_top.set_xticks(np.arange(time[0], time[-1], 1))

plt.show()

标签： pythonmatplotlibsignal-processingspectrogram

Answer 1

我将稍微简化您的示例，减少手动设置，以确保轴正确开箱即用：

import numpy as np
import matplotlib.pyplot as plt

time = np.linspace(-1, 16, 65536)
signal = np.sin(2 * np.pi * time)

fig, ax = plt.subplots(constrained_layout=True)
# In your case, use this instead:
#fig, (ax_top, ax_spec) = plt.subplots(2, 1, constrained_layout=True)

ax.plot(time, 1024 * signal + 1024, c='r')
Pxx, freqs, bins, cax = ax.specgram(signal, NFFT=2048, Fs=4096, noverlap=2048/2, mode='magnitude', pad_to=2048*16)
plt.show()

目标是使信号和图像的时间轴对齐。

频谱图代码的重要部分在matplotlib.mlab._spectral_helper. 频谱图中的每一列都来自一个由matplotlib.mlab._stride_windows. 时间窗口是NFFT 样本宽度，因此n第 th 个时间箱的中心位于0.5 * (time[0] + time[NFFT - 1]) + n * (NFFT - noverlap) * dt。您可以在 中看到此计算mlab._spectral_helper。唯一的区别dt是假定为1 / Fs，并且假定起点为零。

回到 的代码specgram，您可以看到xextent只有在您不手动设置时才会用必要的半像素填充它。总而言之，这意味着您的图像只是移动了-time[0]. 您可以使用 模拟相同的班次xextent。

在我展示如何做到这一点之前，请记住您的值Fs不正确：

>>> 1 / np.diff(time).mean()
3855.0

发生这种情况是因为Fs应该是(len(time) - 1) / (time[-1] - time[0])，不管其他任何事情。您可以从两个样本案例中轻松直观地理解这一点。当您手动将半像素填充应用于范围的边缘时，这将很有用。下面的填充代码直接取自以下情况xextent = None：

time = np.linspace(-1, 16, 65536)
signal = np.sin(2 * np.pi * time)
fs = (len(time) - 1) / (time[-1] - time[0])
half_pixel = 512 / fs   # (NFFT-noverlap) / Fs / 2
half_bin = 1024 / fs  # NFFT / Fs / 2
xextent = (time[0] - half_pixel + half_bin, time[-1] + half_pixel - half_bin)

fig, ax = plt.subplots(constrained_layout=True)

ax.plot(time, (fs / 4) * (signal + 1), c='r')
Pxx, freqs, bins, cax = ax.specgram(signal, NFFT=2048, Fs=fs, noverlap=2048 / 2, mode='magnitude', pad_to=2048 * 16, xextent=xextent)
plt.show()

果然，图像完全位于信号图的中心，信号的末端伸出图像边缘正好半个时间窗口：