金融时间序列：python Matplotlib“specgram”y轴显示周期而不是频率

Question

python Matplotlib 的“频谱图”显示频率（y 轴）与时间（x 轴）的热图对时间序列分析很有用，但我希望 y 轴以周期（ = 1/频率），而不是频率。我仍然在问是否有人有一个完整的工作解决方案来实现这一目标？

紧随其后的 python 代码使用“specgram”生成作者的原始图，并且（当前已注释掉）与使用“mlab.specgram”提供的建议解决方案进行比较。这个建议的解决方案成功地从频率轻松转换为周期 = 1/频率，但没有为作者的示例生成可行的图。

from __future__ import division
from datetime import datetime
import numpy as np
from pandas import DataFrame, Series
import pandas.io.data as web
import pandas as pd
from pylab import plot,show,subplot,specgram
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt

################################################
# obtain data:
ticker = "SPY"
source = "google"
start_date = datetime(1999,1,1)
end_date = datetime(2012,1,1)
qt = web.DataReader(ticker, source, start_date, end_date)
qtC = qt.Close

################################################
data = qtC
fs = 1          # 1 sample / day
nfft = 128

# display the time-series data
fig = plt.figure()
ax1 = fig.add_subplot(311)
ax1.plot(range(len(data)),data)

#----------------

# Original version
##################

# specgram (NOT mlab.specgram) --> gives direct plot, but in Frequency space (want plot in Period, not freq).
ax2 = fig.add_subplot(212)
spec, freq, t = specgram(data, NFFT=nfft, Fs=fs, noverlap=0)

 #----------------
"""
# StackOverflow version (with minor changes to axis titles)
########################

# calcuate the spectrogram
spec, freq, t = mlab.specgram(data, NFFT=nfft, Fs=fs, noverlap=0)

# calculate the bin limits in time (x dir)
# note that there are n+1 fence posts
dt = t[1] - t[0]
t_edge = np.empty(len(t) + 1)
t_edge[:-1] = t - dt / 2.
# however, due to the way the spectrogram is calculates, the first and last bins 
# a bit different:
t_edge[0] = 0
t_edge[-1] = t_edge[0] + len(data) / fs

# calculate the frequency bin limits:
df = freq[1] - freq[0]
freq_edge = np.empty(len(freq) + 1)
freq_edge[:-1] = freq - df / 2.
freq_edge[-1] = freq_edge[-2] + df

# calculate the period bin limits, omit the zero frequency bin
p_edge = 1. / freq_edge[1:]

# we'll plot both
ax2 = fig.add_subplot(312)
ax2.pcolormesh(t_edge, freq_edge, spec)
ax2.set_ylim(0, fs/2)
ax2.set_ylabel('freq.[day^-1]')

ax3 = fig.add_subplot(313)
# note that the period has to be inverted both in the vector and the spectrum,
# as pcolormesh wants to have a positive difference between samples 
ax3.pcolormesh(t_edge, p_edge[::-1], spec[:0:-1])
#ax3.set_ylim(0, 100/fs)
ax3.set_ylim(0, nfft)
ax3.set_xlabel('t [days]')
ax3.set_ylabel('period [days]')
"""

Answer 1

如果你只是问如何以不同的方式显示频谱图，那么它实际上是相当简单的。

需要注意的是，有两个函数称为specgram：matplotlib.pyplot.specgram 和matplotlib.mlab.specgram。两者的区别在于前者绘制频谱图，而后者只计算一个（这就是我们想要的）。

唯一有点棘手的事情是计算颜色网格矩形的边缘位置。我们从specgram 得到以下信息：

• t：时间中心点
• freq: bins 的频率中心

对于时间维度，很容易通过中心计算 bin 限制：

• t_edge[n] = t[0] + (n - .5) * dt，其中dt是两个连续bin的时间差

频率也同样简单：

• f_edge[n] = freq[0] + (n - .5) * df

但我们想使用周期而不是频率。这使得第一个 bin 无法使用，我们将不得不丢弃 DC 分量。

一点代码：

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

# create some data: (fs = sampling frequency)
fs = 2000.
ts = np.arange(10000) / fs
sig = np.sin(500 * np.pi * ts)
sig[5000:8000] += np.sin(200 * np.pi * (ts[5000:8000] + 0.0005 * np.random.random(3000)))

# calcuate the spectrogram
spec, freq, t = mlab.specgram(sig, Fs=fs)

# calculate the bin limits in time (x dir)
# note that there are n+1 fence posts
dt = t[1] - t[0]
t_edge = np.empty(len(t) + 1)
t_edge[:-1] = t - dt / 2.
# however, due to the way the spectrogram is calculates, the first and last bins 
# a bit different:
t_edge[0] = 0
t_edge[-1] = t_edge[0] + len(sig) / fs

# calculate the frequency bin limits:
df = freq[1] - freq[0]
freq_edge = np.empty(len(freq) + 1)
freq_edge[:-1] = freq - df / 2.
freq_edge[-1] = freq_edge[-2] + df

# calculate the period bin limits, omit the zero frequency bin
p_edge = 1. / freq_edge[1:]

# we'll plot both
fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.pcolormesh(t_edge, freq_edge, spec)
ax1.set_ylim(0, fs/2)
ax1.set_ylabel('frequency [Hz]')

ax2 = fig.add_subplot(212)
# note that the period has to be inverted both in the vector and the spectrum,
# as pcolormesh wants to have a positive difference between samples 
ax2.pcolormesh(t_edge, p_edge[::-1], spec[:0:-1])
ax2.set_ylim(0, 100/fs)
ax2.set_xlabel('t [s]')
ax2.set_ylabel('period [s]')

这给出了：