Python,Matplotlib:specgram数据数组值与specgram图不匹配

5
我正在使用 matplotlib.pyplot.specgram 和 matplotlib.pyplot.pcolormesh 制作地震信号的频谱图。 使用 pcolormesh 的原因是我需要在频谱图数据数组上进行算术运算,然后重新绘制结果频谱图(对于三分量地震图 - 东、北和垂直 - 我需要计算水平频谱幅度并将垂直频谱除以水平频谱)。使用频谱图数组数据比单个振幅频谱更容易完成这项工作。
但是我发现,在进行算术操作后,频谱图的绘制结果有意外值。进一步调查后发现,使用 pyplot.specgram 方法制作的频谱图与使用 pyplot.pcolormesh 制作及其返回的数据数组中的值不同。两个绘图/数组应该包含相同的值,但我无法弄清楚为什么它们不同。
例如: 绘制的图形为:
plt.subplot(513)
PxN, freqsN, binsN, imN = plt.specgram(trN.data, NFFT = 20000, noverlap = 0, Fs = trN.stats.sampling_rate, detrend = 'mean', mode = 'magnitude')
plt.title('North')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.clim(0, 150)
plt.colorbar()
#np.savetxt('PxN.txt', PxN)

看起来与情节不同

plt.subplot(514)
plt.pcolormesh(binsZ, freqsZ, PxN)
plt.clim(0,150)
plt.colorbar()

尽管“PxN”数据数组(即每个段的频谱图数据值)是由第一种方法生成并在第二种方法中重复使用,但仍然存在这种情况。您是否知道为什么会发生这种情况?顺便说一句,我意识到我的NFFT值不是平方数,但在我的编码阶段这并不重要。此外,我不知道“imN”数组(pyplot.specgram的第四个返回变量)是什么以及它的用途...

它们有什么不同?你能举个例子吗?你能提供一个trN的示例,让其他人能够重现你的问题吗? - tmdavison
1个回答

18

首先,让我们举一个你所描述的例子,以便其他人能够理解。

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(1)

# Brownian noise sequence
x = np.random.normal(0, 1, 10000).cumsum()

fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(8, 10))

values, ybins, xbins, im = ax1.specgram(x, cmap='gist_earth')
ax1.set(title='Specgram')
fig.colorbar(im, ax=ax1)

mesh = ax2.pcolormesh(xbins, ybins, values, cmap='gist_earth')
ax2.axis('tight')
ax2.set(title='Raw Plot of Returned Values')
fig.colorbar(mesh, ax=ax2)

plt.show()

输入图像描述

幅度差异

您会立即注意到绘制值的幅度差异。

默认情况下,plt.specgram 不会绘制其返回的 "原始" 值。相反,它将它们缩放为分贝(换句话说,它绘制振幅的 10 * log10)。如果您不希望缩放物体,请指定 scale="linear"。但是,对于查看频率组成,对数刻度将是最有意义的。

考虑到这一点,让我们模仿 specgram 的操作:

plotted = 10 * np.log10(values)

fig, ax = plt.subplots()
mesh = ax.pcolormesh(xbins, ybins, plotted, cmap='gist_earth')

ax.axis('tight')
ax.set(title='Plot of $10 * log_{10}(values)$')
fig.colorbar(mesh)

plt.show()

在此输入图片描述

使用对数色彩刻度

或者,我们可以在图像上使用对数标准化,并获得类似的结果,但更清晰地表明颜色值处于对数刻度上:

from matplotlib.colors import LogNorm

fig, ax = plt.subplots()
mesh = ax.pcolormesh(xbins, ybins, values, cmap='gist_earth', norm=LogNorm())

ax.axis('tight')
ax.set(title='Log Normalized Plot of Values')
fig.colorbar(mesh)

plt.show()

输入图像描述

imshowpcolormesh

需要注意的是,我们展示的示例没有应用任何插值方法,而原始的 specgram 图有。 specgram 使用 imshow,而我们使用 pcolormesh 进行绘图。对于这种情况(网格间距相等),我们可以选择使用任意一种。

在这种情况下,imshowpcolormesh 都是非常好的选择。然而,如果您要处理大型数组,imshow 的性能会显著提高。因此,即使您不想进行插值(例如,使用 interpolation='nearest' 关闭插值),也可能考虑使用它。

例如:

extent = [xbins.min(), xbins.max(), ybins.min(), ybins.max()]

fig, ax = plt.subplots()
mesh = ax.imshow(values, extent=extent, origin='lower', aspect='auto',
                 cmap='gist_earth', norm=LogNorm())

ax.axis('tight')
ax.set(title='Log Normalized Plot of Values')
fig.colorbar(mesh)

plt.show()

在此输入图像描述


太棒了,非常感谢,这解释了一切,现在我明白了。 - Nukolas
我觉得这个答案应该有数百个赞,它为我节省了很多时间。 - eric
需要为Fs添加值:specgram(x, Fs=1, cmap='gist_earth') - Naomi Fridman

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接