[PYTHON] Illustrez instantanément la période prédominante dans les données de séries chronologiques à l'aide de l'analyse spectrale
Transformez instantanément les données 2D en diagramme à l'aide de matplotlib de python Le deuxième de la série "Illustrer en un instant avec matplotlib de python" suit. Quoi qu'il en soit, pour ceux qui veulent facilement ** trouver un cycle exceptionnel avec des données chronologiques unidimensionnelles comme indiqué ci-dessous **.
En utilisant python,
- [Analyse spectrale](https://ja.wikipedia.org/wiki/%E5%91%A8%E6%B3%A2%E6%95%B0%E3%82%B9%E3%83%9A%E3 Effectuer le calcul de% 82% AF% E3% 83% 88% E3% 83% AB) --Illustration des résultats des calculs
Faites tout en même temps.
1. Code
Définissez les classes suivantes à l'avance.
spectra.py
# coding:utf-8
from scipy import fftpack
import numpy as np
import matplotlib.pyplot as plt
class Spectra(object):
def __init__(self, t, f, time_unit):
"""
- t :Valeur de l'axe des temps
- f :Valeur des données
- time_unit :Unité de l'axe des temps
- Po :Densité du spectre de puissance
"""
assert t.size == f.size #Assurez-vous que la longueur de l'axe des temps et la longueur des données sont identiques
assert np.unique(np.diff(t)).size == 1 #Assurez-vous que tous les intervalles de temps sont constants
self.time_unit = time_unit #Unité de temps
T = (t[1] - t[0]) * t.size
self.period = 1.0 / (np.arange(t.size / 2)[1:] / T)
#Calculer la densité spectrale de puissance
f = f - np.average(f) #Zéro la moyenne.
F = fftpack.fft(f) #Transformée de Fourier Rapide
self.Po = np.abs(F[1:(t.size // 2)]) ** 2 / T
def draw_with_time(self, fsizex=8, fsizey=6, print_flg=True, threshold=1.0):
#Tracez la densité spectrale de puissance au fil du temps sur l'axe horizontal
fig, ax = plt.subplots(figsize=(fsizex, fsizey)) #Spécifier la taille de la figure
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_xlabel(self.time_unit)
ax.set_ylabel("Power Spectrum Density")
ax.plot(self.period, self.Po)
if print_flg: #Tracez une ligne où la valeur de densité spectrale de puissance est supérieure au seuil et décrivez la valeur de période.
dominant_periods = self.period[self.Po > threshold]
print(dominant_periods, self.time_unit +
' components are dominant!')
for dominant_period in dominant_periods:
plt.axvline(x=dominant_period, linewidth=0.5, color='k')
ax.text(dominant_period, threshold,
str(round(dominant_period, 3)))
return plt
Comme vous pouvez le voir, [FFT (transformée de Fourier discrète)](https://ja.wikipedia.org/wiki/%E9%AB%98%E9%80%9F%E3 % 83% 95% E3% 83% BC% E3% 83% AA% E3% 82% A8% E5% A4% 89% E6% 8F% 9B) est effectuée pour calculer la densité spectrale de puissance. La période à laquelle la densité spectrale de puissance montre une valeur plus grande que l'environnement est la période prédominante dans la série temporelle.
2. Vérification
Créez vos propres données de séries chronologiques et utilisez-les pour vérifier le code ci-dessus.
#coding:utf-8
import numpy as np
import matplotlib.pyplot as plt
from spectra import Spectra
#Création de données de séries chronologiques
#Données mensuelles sur 30 ans(Autrement dit, le nombre de données est de 360)
#Grande année(=12 mois)10 ans progressifs en plus du cycle(=120 mois)Un peu de bruit à chaque fois
N = 360
t = np.arange(0, N)
td = t * np.pi / 6.0
f = np.sin(td) + 35.0 + 0.2 * np.sin(td * 0.1) + np.random.randn(N) * 0.1
#Dessin de la série chronologique originale
plt.figure(figsize=(20, 6))
plt.plot(t, f)
plt.xlim(0, N)
plt.xlabel('Month')
plt.show()
#Dessin de cycles exceptionnels
spectra = Spectra(t, f, 'Month')
plt = spectra.draw_with_time()
plt.show()
3. Résultat
Le diagramme de données de série chronologique d'origine ressemble à celui en haut de cet article. Et le diagramme du cycle exceptionnel est le suivant.
[ 120. 12.] Month components are dominant!
Comme vous pouvez le voir, c ’est un cycle exceptionnel.
- 12 mois --120 mois
Il a extrait à merveille. Si les valeurs aux deux extrémités des données de la série chronologique d'origine sont sensiblement différentes, [Window function](https://ja.wikipedia.org/wiki/%E7%AA%93%E9%96%A2%] Traitez les données sous une forme adaptée à l'analyse spectrale, par exemple en appliquant E6% 95% B0).
Les références
- [Analyse pratique des données climatiques sous UNIX / Windows / Macintosh-Trois points que les rédacteurs de rapports et d'articles sur les sciences du climat, la météorologie, les sciences de l'océan, etc. devraient connaître](https://www.amazon. co.jp/UNIX-Windows-Macintosh%E3%82%92%E4%BD%BF%E3%81%A3%E3%81%9F%E5%AE%9F%E8%B7%B5-%E6%B0 % 97% E5% 80% 99% E3% 83% 87% E3% 83% BC% E3% 82% BF% E8% A7% A3% E6% 9E% 90-% E6% B0% 97% E5% 80% 99% E5% AD% A6-% E6% B0% 97% E8% B1% A1% E5% AD% A6-% E6% B5% B7% E6% B4% 8B% E5% AD% A6% E3% 81% AA% E3% 81% A9% E3% 81% AE% E5% A0% B1% E5% 91% 8A% E6% 9B% B8-% E8% AB% 96% E6% 96% 87% E3% 82% 92 % E6% 9B% B8% E3% 81% 8F% E4% BA% BA% E3% 81% 8C% E7% 9F% A5% E3% 81% A3% E3% 81% A6% E3% 81% 8A% E3 % 81% 8D% E3% 81% 9F% E3% 81% 843% E3% 81% A4% E3% 81% AE% E3% 83% 9D% E3% 82% A4% E3% 83% B3% E3% 83 % 88-% E6% 9D% BE% E5% B1% B1 / dp / 4772241221 / ref = sr_1_1? Ie = UTF8 & qid = 14482492333 & sr = 8-1 & mots-clés =% E6% B0% 97% E5% 80% 99% E3% 83 % 87% E3% 83% BC% E3% 82% BF% E8% A7% A3% E6% 9E% 90)
- https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.fftpack.fft.html
Recommended Posts