Turn 2D data into a diagram instantly using python's matplotlib The second in the series, "Illustrate in an instant using python's matplotlib". Anyway, for those who want to easily ** find an outstanding cycle with one-dimensional time series data as shown below **.

Using python,

-[Spectrum analysis](https://ja.wikipedia.org/wiki/%E5%91%A8%E6%B3%A2%E6%95%B0%E3%82%B9%E3%83%9A%E3 Perform calculation of% 82% AF% E3% 83% 88% E3% 83% AB) --Illustration of calculation results

Do all at once.

1. Code

Set the following classes in advance.

`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 :Time axis value
         - f :Data value
         - time_unit :Time axis unit
         - Po :Power spectral density
        """
        assert t.size == f.size  #Make sure that the length of the time axis and the length of the data are the same
        assert np.unique(np.diff(t)).size == 1  #Make sure all time intervals are constant
        self.time_unit = time_unit   #Unit of time
        T = (t[1] - t[0]) * t.size
        self.period = 1.0 / (np.arange(t.size / 2)[1:] / T)

        #Calculate power spectral density
        f = f - np.average(f)         #Zero the average.
        F = fftpack.fft(f)                          #Fast Fourier transform
        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):
        #Plot the power spectral density over time on the horizontal axis
        fig, ax = plt.subplots(figsize=(fsizex, fsizey))   #Specifying the size of the 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:   #Draw a line where the power spectral density value is greater than the threshold and describe the period value.
            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

As you can see, FFT (Discrete Fourier Transform) in chronological order using scipy's fft function % 83% 95% E3% 83% BC% E3% 83% AA% E3% 82% A8% E5% A4% 89% E6% 8F% 9B) is performed to calculate the power spectral density. The period at which the power spectral density shows a larger value than the surroundings is the period that is predominant in the time series.

2. Verification

Create your own time series data and use it to verify the above code.

#coding:utf-8

import numpy as np
import matplotlib.pyplot as plt
from spectra import Spectra

#Creating time series data
#Monthly data for 30 years(That is, the number of data is 360)
#Big year(=12 months)10 years gradual in addition to the cycle(=120 months)A small amount of noise at each time
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

#Drawing the original time series
plt.figure(figsize=(20, 6))
plt.plot(t, f)
plt.xlim(0, N)
plt.xlabel('Month')
plt.show()

#Drawing of outstanding cycles
spectra = Spectra(t, f, 'Month')
plt = spectra.draw_with_time()
plt.show()

3. Result

The original time series data diagram looks like the one at the top of this post. And the diagram of the outstanding cycle is as follows.

[ 120. 12.] Month components are dominant!

As you can see, it ’s an outstanding cycle.

12 months --120 months

It extracted wonderfully. If the values at both ends of the original time series data are significantly different, [Window function](https://ja.wikipedia.org/wiki/%E7%AA%93%E9%96%A2%] Process the data into a form suitable for spectral analysis, such as by applying E6% 95% B0).

References

-[Practical climate data analysis using UNIX / Windows / Macintosh-Three points that writers of reports and papers on climate science, meteorology, ocean science, etc. should know](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 & keywords =% 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

Instantly illustrate the predominant period in time series data using spectrum analysis

1. Code

spectra.py

2. Verification

3. Result

`spectra.py`