Финансовые временные ряды: спектрограмма Python Matplotlib по оси Y, отображающая период вместо частоты

Отображение «спектрограммы» 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]')
"""

person TonyMorland    schedule 16.07.2014    source источник


Ответы (1)


Если вы только спрашиваете, как отображать спектрограмму по-другому, то это на самом деле довольно просто.

Следует отметить, что есть две функции с именем specgram: matplotlib.pyplot.specgram и matplotlib.mlab.specgram. Разница между ними в том, что первый рисует спектрограмму, а второй вычисляет только ее (и это то, что нам нужно).

Единственная небольшая сложность заключается в вычислении положения краев прямоугольника цветовой сетки. Получаем из specgram следующее:

  • t: центральные точки во времени
  • freq: частотные центры бинов

Для измерения времени легко вычислить пределы бинов по центрам:

  • t_edge[n] = t[0] + (n - .5) * dt, где dt — разница во времени двух последовательных бинов

Это было бы так же просто для частот:

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

но мы хотим использовать период вместо частоты. Это делает первый бункер непригодным для использования, и нам придется выбросить компонент постоянного тока.

Немного кода:

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]')

Это дает:

введите здесь описание изображения

person DrV    schedule 16.07.2014
comment
Преобразование частоты в период точно соответствует требованиям, но предлагаемое отображение цветовой сетки не работает, как показано в приложенном примере, в котором теперь показаны исходный код и рисунок автора по сравнению с предлагаемым решением. - person TonyMorland; 17.07.2014
comment
@TonyMorland: Боюсь, я не совсем понимаю вас. Что значит не работает в данном контексте? По крайней мере, set_ylim выглядит немного странно, так как ось должна быть в единицах времени. Или, может быть, проблема в том, что график должен быть линейным, но только метки оси Y должны отражать период? - person DrV; 17.07.2014
comment
ТКС DrV. я не могу вставить график, чтобы показать вам, но если вы попробуете версии кода с комментариями и без комментариев, которые теперь вставлены под исходным вопросом, вы увидите проблему. - person TonyMorland; 17.07.2014