diff --git a/README.rst b/README.rst
index 531d61c..ccfbde0 100644
--- a/README.rst
+++ b/README.rst
@@ -59,6 +59,10 @@ where a = 0.2, b = 0.2 and c = 5.7.
 
 .. image:: https://raw.githubusercontent.com/capitanov/chaospy/master/img/Rossler_3D.png?sanitize=true
 
+Spectrum and auto-correlations::
+
+.. image:: https://raw.githubusercontent.com/capitanov/chaospy/master/img/Lorenz_Spectrum.png?sanitize=true
+
 Source code
 ~~~~~~~~~~~
 
diff --git a/img/Lorenz_Spectrum.png b/img/Lorenz_Spectrum.png
new file mode 100644
index 0000000..8249c8a
Binary files /dev/null and b/img/Lorenz_Spectrum.png differ
diff --git a/src/_attractors.py b/src/_attractors.py
deleted file mode 100644
index e01c747..0000000
--- a/src/_attractors.py
+++ /dev/null
@@ -1,468 +0,0 @@
-"""Chaotic attractors.
-
-Description   : This is the main module which collects attractors
-                together. You can choose: Lorenz, Rossler, Rikitake,
-                Nose-Hoover, Chua, Lottki, Duffing, Wang etc.
-
-    Attention!
-    You should set correct initial values of coordinates [Xo, Yo, Zo]
-    because of some attractors has an overflow effect.
-
-    ----------
-    Lorenz:
-    ----------
-    Lorenz attractor is ordinary differential equation (ODE) of
-    3rd order system.
-    In 1963, E. Lorenz developed a simplified mathematical model for
-    atmospheric convection.
-
-    Lorenz equations are:
-        dx/dt = sigma * (y - x)
-        dy/dt = rho * x - y - x * z
-        dz/dt = x * y - beta * z
-
-    where beta, rho and sigma - are Lorenz system parameters. Default
-    values are: beta = 8/3, rho = 28 and sigma = 10.
-
-    Wiki: If rho < 1 then there is only one equilibrium point.
-
-    ----------
-    Rikitake:
-    ----------
-    Rikitake system is ordinary differential equation (ODE) of
-    3rd order system.
-    Rikitake system attempts to explain the reversal of the Earth’s
-    magnetic field.
-
-    Rikitake equations are:
-        dx/dt = -mu * x + z * y
-        dy/dt = -mu * y + x * (z - a)
-        dz/dt = 1 - x * y
-
-    where a, mu - are Rikitake system parameters. Default values are
-    a = 5, mu = 2 or a = mu = 1.
-
-    ----------
-    Rossler:
-    ----------
-    Rossler attractor is the attractor for the Rössler system. Rossler
-    attractor is a system of three non-linear ordinary differential
-    equations.
-
-    Rossler equations are:
-        dx/dt = -(y + z)
-        dy/dt = x + a * y
-        dz/dt = b + z * (x - c)
-
-    where a, b and sigma - are Rossler system parameters. Default
-    values are: a = 0.2, b = 0.2 and c = 5.7.
-
-    1) a = 0.2, b = 0.2 and c = 5.7 (Standard model)
-    2) a = 0.1, b = 0.1 and c = 14 (another useful parameters)
-    3) a = 0.5, b = 1.0 and c = 3 (J. C. Sprott)
-
-    Wiki (Varying parameters);
-
-    Varying a:
-    b = 0.2 and c = 5.7 are fixed. Change a:
-
-    a <= 0      : Converges to the centrally located fixed point
-    a = 0.1     : Unit cycle of period 1
-    a = 0.2     : Standard parameter value selected by Rössler, chaotic
-    a = 0.3     : Chaotic attractor, significantly more Möbius strip-like
-                (folding over itself).
-    a = 0.35    : Similar to .3, but increasingly chaotic
-    a = 0.38    : Similar to .35, but increasingly chaotic
-
-    Varying b:
-    a = 0.2 and c = 5.7 are fixed. Change b:
-
-    If b approaches 0 the attractor approaches infinity, but if b would
-    be more than a and c, system becomes not a chaotic.
-
-    Varying c:
-    a = b = 0.1 are fixed. Change c:
-
-    c = 4       : period-1 orbit,
-    c = 6       : period-2 orbit,
-    c = 8.5     : period-4 orbit,
-    c = 8.7     : period-8 orbit,
-    c = 9       : sparse chaotic attractor,
-    c = 12      : period-3 orbit,
-    c = 12.6    : period-6 orbit,
-    c = 13      : sparse chaotic attractor,
-    c = 18      : filled-in chaotic attractor.
-
-    ----------
-    Nose–Hoover:
-    ----------
-    Nose–Hoover oscillator is ordinary differential equation (ODE) of
-    3rd order system.
-    Nose–Hoover system has only five terms and two quadratic
-    nonlinearities.
-
-    Nose–Hoover equations are:
-        dx/dt = y
-        dy/dt = y * z - x
-        dz/dt = 1 - y * y
-
-    Nose–Hoover system doesn't have any system parameters.
-
-    ----------
-    Wang attractor:
-    ----------
-    Wang system (improved Lorenz model) as classic chaotic attractor
-
-    Wang equations are:
-        dx/dt = x - y*z
-        dy/dt = x - y + x*z
-        dz/dt = -3*z + x*y
-
-    ----------
-    Duffing map:
-    ----------
-    It is a discrete-time dynamical system (2nd order)
-    The map depends on the two constants: a and b.
-    Z coordinate is a linear function.
-
-    Duffing equations are:
-    Eq. 1:
-        dx/dt = y
-        dy/dt = -a*y - x**3 + b * cos(z)
-        dz/dt = 1
-    where a = 0.1 and b = 11 (default parameters)
-
-    Eq. 2:
-        dx/dt = y
-        dy/dt = a*y - y**3 - b*x
-        dz/dt = 1
-    where a = 2.75 and b = 0.2 (default parameters)
-
-    ----------
-    Chua circuit:
-    ----------
-    Chua circuit. This is a simple electronic circuit that exhibits
-    classic chaotic behavior.
-
-    Chua equations are:
-    Eq. 1:
-        dx/dt = alpha * (y - x - ht)
-        dy/dt = x - y + z
-        dz/dt = -beta * y
-
-    where ht = mu1*x + 0.5*(mu0 - mu1)*(np.abs(x + 1) - np.abs(x - 1))
-    and alpha, beta, mu0 and mu1 - are Chua system parameters.
-
-    Default values are:
-    alpha = 15.6
-    beta = 28
-    mu0 = -1.143
-    mu1 = -0.714
-
-    Eq. 2:
-        dx/dt = 0.3*y + x - x**3
-        dy/dt = x + z
-        dz/dt = y
-
-    ----------
-    Lotka–Volterra:
-    ----------
-    The Lotka–Volterra equations, also known as the predator–prey
-    equations.
-
-    Chaotic Lotka-Volterra model require a careful tuning of
-    parameters and are even less likely to exhibit chaos as the number
-    of species increases.
-
-    Lotka–Volterra equations are:
-    Eq. 1:
-        dx/dt = x * (1 - x - 9*y)
-        dy/dt = -y * (1 - 6*x - y + 9*z)
-        dz/dt = z * (1 - 3*x - z)
-
-    Be careful! Init parameters of x, y, z should be set right.
-    For example, [x, y, z] = [0.6; 0.2; 0.01].
-
-------------------------------------------------------------------------
-
-GNU GENERAL PUBLIC LICENSE
-Version 3, 29 June 2007
-
-Copyright (c) 2019 Kapitanov Alexander
-
-This program is free software: you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation, either version 3 of the License, or
-(at your option) any later version.
-
-You should have received a copy of the GNU General Public License
-along with this program.  If not, see <http://www.gnu.org/licenses/>.
-
-THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
-APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
-HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT
-WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT
-NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
-FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND
-PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE
-DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR
-OR CORRECTION.
-
-------------------------------------------------------------------------
-"""
-
-
-# Authors       : Alexander Kapitanov
-# ...
-# Contacts      : <empty>
-# ...
-# Release Date  : 2019/05/31
-# License       : GNU GENERAL PUBLIC LICENSE
-
-
-import matplotlib.pyplot as plt
-import numpy as np
-from mpl_toolkits.mplot3d import axes3d  # noqa # pylint: disable=unused-import
-from scipy.fftpack import fft, fftshift
-from scipy.stats import gaussian_kde, kurtosis, skew
-from src.attractors.chua import chua
-from src.attractors.duffing import duffing
-from src.attractors.lorenz import lorenz
-from src.attractors.lotka_volterra import lotka_volterra
-from src.attractors.nose_hoover import nose_hoover
-from src.attractors.rikitake import rikitake
-from src.attractors.rossler import rossler
-from src.attractors.wang import wang
-
-# #####################################################################
-# Functions
-# #####################################################################
-
-# TODO: Reimplement this file! 2020/07/14
-
-
-def sel_chaos(chtype=None, xt=0, yt=0, zt=0, **kwargs):
-    """
-    Select chaotic system function by name
-
-    Parameters
-    ----------
-    chtype : str
-        Select function
-    xt, yt, zt : float
-        Input coordinates Z, Y, Z respectively
-    kwargs : float
-        Chaotic system parameters
-
-    """
-    ch_case = chtype.casefold()
-    if ch_case == "wang":
-        return wang(xt, yt, zt)
-    elif ch_case == "nose-hoover":
-        return nose_hoover(xt, yt, zt)
-    elif ch_case == "lotka":
-        return lotka_volterra(xt, yt, zt)
-    elif ch_case == "chua":
-        return chua(xt, yt, zt, **kwargs)
-    elif ch_case == "lorenz":
-        return lorenz(xt, yt, zt, **kwargs)
-    elif ch_case == "rossler":
-        return rossler(xt, yt, zt, **kwargs)
-    elif ch_case == "duffing":
-        return duffing(xt, yt, zt, **kwargs)
-    elif ch_case == "rikitake":
-        return rikitake(xt, yt, zt, **kwargs)
-    else:
-        raise Exception("Error: you should set the correct chaotic system")
-
-
-def sel_args(chtype=None, args=None):
-    """
-    Select arguments for chaotic system
-
-    Parameters
-    ----------
-    chtype : str
-        Select function
-    args : dict
-        Chaotic system arguments (as a dict)
-    """
-    ch_case = chtype.casefold()
-    if ch_case in ("lorenz", "rossler", "duffing", "chua", "rikitake"):
-        return args[ch_case]
-    return {}
-
-
-# #####################################################################
-# Input parameters
-# #####################################################################
-
-CHTYPE = "Wang"  # Chaotic system name
-CHARGS = {
-    "rikitake": {"a": 1, "mu": 1},
-    "duffing": {"a": 0.1, "b": 11},
-    "rossler": {"a": 0.2, "b": 0.2, "c": 5.7},
-    "lorenz": {"sigma": 10, "beta": 8 / 3, "rho": 28},
-    "chua": {"alpha": 0.1, "beta": 28, "mu0": -1.143, "mu1": -0.714},
-}  # Set arguments
-
-NW = 5000  # Number of ODE's dots
-DT = 100  # Step for equations
-NFFT = 2 ** 12  # Number of FFT dots
-
-
-# #####################################################################
-# Calculate attractor
-# #####################################################################
-def calc_chaos():
-    print("Start analyzing the next chaotic system: [%s]\n" % CHTYPE)
-
-    # Create zero arrays for coordinates
-    xt, yt, zt = np.zeros(NW), np.zeros(NW), np.zeros(NW)
-    # Set initial values for [X, Y, Z]
-    xt[0], yt[0], zt[0] = 1.0, 0.5, 1.0
-
-    for ii in range(NW - 1):
-        x_nt, y_nt, z_nt = sel_chaos(chtype=CHTYPE, xt=xt[ii], yt=yt[ii], zt=zt[ii], **sel_args(CHTYPE, CHARGS),)
-        xt[ii + 1] = xt[ii] + (x_nt / DT)
-        yt[ii + 1] = yt[ii] + (y_nt / DT)
-        zt[ii + 1] = zt[ii] + (z_nt / DT)
-
-    # #####################################################################
-    # Calculate standardized moments
-    # #####################################################################
-
-    ds_3d = ("X", "Y", "Z")
-    ch_3d = np.stack((xt, yt, zt))
-
-    nn, mm = ch_3d.shape[0], ch_3d.shape[1]
-    print("Calculate mean, variance, skewness, kurtosis and median for each " "coordinate of chaotic system:")
-
-    m15 = np.zeros([nn, mm])
-    for ii in range(nn):
-        m15[ii, 0] = np.mean(ch_3d[ii, :])  # Mean
-        m15[ii, 1] = np.var(ch_3d[ii, :])  # Variance
-        m15[ii, 2] = skew(ch_3d[ii, :])  # Skewness
-        m15[ii, 3] = kurtosis(ch_3d[ii, :])  # Kurtosis
-        m15[ii, 4] = np.median(ch_3d[ii, :])  # Median
-
-        print("%s axis: " % ds_3d[ii], end="")
-        print("M[t]:  {:+.4f},  ".format(m15[ii, 0]), end="")
-        print("D[t]:  {:+.4f},  ".format(m15[ii, 1]), end="")
-        print("A[t]:  {:+.4f},  ".format(m15[ii, 2]), end="")
-        print("E[t]:  {:+.4f},  ".format(m15[ii, 3]), end="")
-        print("Md[t]: {:+.4f}.  ".format(m15[ii, 4]))
-
-    # Find Probability density function
-    n_pdf = 3000
-    p_axi = np.zeros([nn, n_pdf])
-    d_kde = np.zeros([nn, n_pdf])
-    for ii in range(nn):
-        p_axi[ii] = np.linspace(ch_3d[ii, :].min(), ch_3d[ii, :].max(), n_pdf)
-        d_kde[ii] = gaussian_kde(ch_3d[ii, :]).evaluate(p_axi[ii, :])
-        d_kde[ii] /= d_kde[ii].max()
-
-    # #####################################################################
-    # Calculate Autocorrelation (Cross-correlation) and FFT
-    # #####################################################################
-
-    # TODO implement FFT and autocorrelation
-    d_acf = np.zeros([nn, mm])
-    d_fft = np.zeros([nn, NFFT], dtype=np.complex)
-    for ii in range(nn):
-        d_acf[ii] = np.correlate(ch_3d[ii, :], ch_3d[ii, :], "same")
-        d_acf[ii] /= d_acf[ii].max()
-        d_fft[ii] = fft(ch_3d[ii, 1:NFFT], NFFT)
-
-    d_fft = np.abs(fftshift(d_fft, axes=1))
-    d_fft /= np.max(d_fft)
-    d_fft = 20 * np.log10(d_fft)
-
-    plt.rc("font", size=8)  # controls default text sizes
-    plt.rc("axes", titlesize=10)  # fontsize of the axes title
-    plt.rc("axes", labelsize=10)  # fontsize of the x and y labels
-    plt.rc("legend", fontsize=8)  # legend fontsize
-
-    lim_xyz = [(np.min(ch_3d[ii]), np.max(ch_3d[ii])) for ii in range(3)]
-
-    # Plot 3D model
-    # fig1 = plt.figure('3D model of chaotic system')
-    # ax = fig1.gca(projection='3d')
-    # ax.plot(xt, yt, zt, 'o-', linewidth=0.1, markersize=0.3)
-    # ax.set_xlabel("X")
-    # ax.set_ylabel("Y")
-    # ax.set_zlabel("Z")
-    # ax.set_xlim(lim_xyz[0])
-    # ax.set_ylim(lim_xyz[1])
-    # ax.set_zlim(lim_xyz[2])
-    # ax.set_title("3D Chaotic Attractor")
-    # plt.tight_layout()
-
-    # Plot 2D coordinates of 3D model
-    fig2 = plt.figure("3D Coordinates")
-    plt.subplot(2, 2, 1)
-    plt.plot(yt, xt, linewidth=0.75)
-    plt.grid()
-    plt.xlabel("X")
-    plt.ylabel("Y")
-    plt.xlim(lim_xyz[1])
-    plt.ylim(lim_xyz[0])
-
-    plt.subplot(2, 2, 2)
-    plt.plot(yt, zt, linewidth=0.75)
-    plt.grid()
-    plt.xlabel("Z")
-    plt.ylabel("Y")
-    plt.xlim(lim_xyz[1])
-    plt.ylim(lim_xyz[2])
-
-    plt.subplot(2, 2, 3)
-    plt.plot(zt, xt, linewidth=0.75)
-    plt.grid()
-    plt.xlabel("X")
-    plt.ylabel("Z")
-    plt.xlim(lim_xyz[2])
-    plt.ylim(lim_xyz[0])
-
-    ax = fig2.add_subplot(2, 2, 4, projection="3d")
-    ax.plot(xt, yt, zt, linewidth=0.7)
-    ax.set_xlabel("X")
-    ax.set_ylabel("Y")
-    ax.set_zlabel("Z")
-    plt.tight_layout()
-
-    # Plot 2D data in timex
-    plt.figure("Coordinates evolution in time")
-    for ii in range(nn):
-        plt.subplot(2, 2, ii + 1)
-        plt.plot(ch_3d[ii], linewidth=0.75)
-        plt.grid()
-        plt.ylabel(ds_3d[ii])
-        plt.xlim([0, NW - 1])
-        plt.ylim(lim_xyz[ii])
-    plt.tight_layout()
-
-    # Plot Probability density function
-    plt.figure("Probability density function")
-    for ii in range(nn):
-        plt.plot(d_kde[ii], ".")
-        plt.xlim([0, n_pdf - 1])
-        plt.grid()
-    plt.tight_layout()
-
-    # Plot Autocorrelation and Spectrum
-    plt.figure("Autocorrelation and Spectrum")
-    for ii in range(nn):
-        plt.subplot(2, 3, ii + 1)
-        plt.plot(d_acf[ii], linewidth=0.75)
-        plt.grid()
-        plt.subplot(2, 3, ii + 1 + 3)
-        plt.plot(d_fft[ii], linewidth=0.75)
-        plt.grid()
-    plt.tight_layout()
-
-    plt.show()
-
-
-# Execute main function: Chaotic system
-calc_chaos()
diff --git a/src/dynamic_system.py b/src/dynamic_system.py
index 2e52513..c5a2a65 100644
--- a/src/dynamic_system.py
+++ b/src/dynamic_system.py
@@ -40,7 +40,6 @@
 
 from typing import Optional
 
-import numpy as np
 import pandas as pd
 from src.utils.calculator import Calculator
 from src.utils.drawer import PlotDrawer
@@ -61,7 +60,6 @@ class DynamicSystem:
     calculator: Calculator()
     settings: Settings()
 
-
     Examples
     --------
     # TODO: Add examples
@@ -102,10 +100,10 @@ def initialize(self, input_args: Optional[tuple] = None, show_logs: bool = False
         # Update main calculator
         self.calculator = Calculator()
 
-    def collect_statistics(self, points: np.ndarray):
+    def collect_statistics(self):
         math_dict = {}
-        _min_max = self.calculator.check_min_max(points)
-        _moments = self.calculator.check_moments(points)
+        _min_max = self.calculator.check_min_max()
+        _moments = self.calculator.check_moments()
         math_dict.update({"Min": _min_max[0]})
         math_dict.update({"Max": _min_max[1]})
         math_dict.update(_moments)
@@ -115,22 +113,24 @@ def collect_statistics(self, points: np.ndarray):
     def run(self):
         # Get vector of coordinates
         _points = self.model.get_coordinates()
+        self.calculator.coordinates = _points
 
         # Calculate
-        stats = self.collect_statistics(_points)
+        stats = self.collect_statistics()
         if self.settings.show_logs:
             print(f"[INFO]: Show statistics:\n{stats}\n")
 
-        self.calculator.check_probability(_points)
-        self.calculator.calculate_fft(_points)
+        self.calculator.check_probability()
+        spectrums = self.calculator.calculate_spectrum()
+        correlations = self.calculator.calculate_correlation()
 
         # Draw results
         if self.settings.show_plots or self.settings.save_plots:
             self.drawer.coordinates = _points
-
+            self.drawer.show_spectrum_and_correlation(spectrums, correlations)
             # self.drawer.show_time_plots()
             # self.drawer.show_3d_plots()
-            self.drawer.make_3d_plot_gif(50)
+            # self.drawer.make_3d_plot_gif(50)
             # self.drawer.show_all_plots()
 
 
@@ -139,13 +139,13 @@ def run(self):
         "--init_point",
         "1 -1 2",
         "--points",
-        "2000",
+        "3000",
         "--step",
         "50",
         "--save_plots",
         # "--show_plots",
         "--add_2d_gif",
-        "rossler",
+        "lorenz",
     )
 
     dynamic_system = DynamicSystem(input_args=command_line, show_log=True)
diff --git a/src/utils/calculator.py b/src/utils/calculator.py
index 0fdd01c..5e079de 100644
--- a/src/utils/calculator.py
+++ b/src/utils/calculator.py
@@ -42,40 +42,52 @@
 from typing import Tuple
 
 import numpy as np
-from scipy.stats import kurtosis, skew
+from scipy.fftpack import fft, fftshift
+from scipy.stats import gaussian_kde, kurtosis, skew
 
 
 class Calculator:
-    """Main class for calculate math parameters.
+    """Main class for calculate math parameters: FFTs, Auto-Correlation, KDE (Prob) etc.
 
     See Also:
     -----
 
     """
 
-    def __init__(self):
-        pass
+    def __init__(self, kde_dots: int = 1000, fft_dots: int = 4096):
+        self.kde_dots = kde_dots
+        self.fft_dots = fft_dots
+        self._coordinates = None
 
-    @staticmethod
-    def check_min_max(coordinates: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
-        """Calculate minimum and maximum for X, Y, Z coordinates.
+    def __len__(self):
+        return len(self.coordinates)
+
+    @property
+    def coordinates(self) -> np.ndarray:
+        return self._coordinates
+
+    @coordinates.setter
+    def coordinates(self, value: np.ndarray):
+        """3D coordinates.
 
         Parameters
         ----------
-        coordinates: np.array
+        value : np.ndarray
             Numpy 3D array of dynamic system coordinates.
+
         """
-        return np.min(coordinates, axis=0), np.max(coordinates, axis=0)
+        self._coordinates = value
 
-    @staticmethod
-    def check_moments(coordinates: np.ndarray, is_common: bool = False) -> dict:
+    def check_min_max(self) -> Tuple[np.ndarray, np.ndarray]:
+        """Calculate minimum and maximum for X, Y, Z coordinates.
+        """
+        return np.min(self.coordinates, axis=0), np.max(self.coordinates, axis=0)
+
+    def check_moments(self, is_common: bool = False) -> dict:
         """Calculate stochastic parameters: mean, variance, skewness, kurtosis etc.
 
         Parameters
         ----------
-        coordinates: np.array
-            Numpy 3D array of dynamic system coordinates.
-
         is_common : bool
             If False - method returns moments for each coordinate. Otherwise
             returns moments over all ndarray. Similar for axis or axes along
@@ -84,21 +96,62 @@ def check_moments(coordinates: np.ndarray, is_common: bool = False) -> dict:
         """
         axis = None if is_common else 0
         return {
-            "Mean": np.mean(coordinates, axis=axis),
-            "Variance": np.var(coordinates, axis=axis),
-            "Skewness": skew(coordinates, axis=axis),
-            "Kurtosis": kurtosis(coordinates, axis=axis),
-            "Median": np.median(coordinates, axis=axis),
+            "Mean": np.mean(self.coordinates, axis=axis),
+            "Variance": np.var(self.coordinates, axis=axis),
+            "Skewness": skew(self.coordinates, axis=axis),
+            "Kurtosis": kurtosis(self.coordinates, axis=axis),
+            "Median": np.median(self.coordinates, axis=axis),
         }
 
-    def check_probability(self, coordinates: np.ndarray):
-        pass
-        # TODO: Implement this method!
+    def check_probability(self):
+        """Check probability for each chaotic coordinates.
+
+        """
+        p_axi = np.zeros([3, self.kde_dots])
+        d_kde = np.zeros([3, self.kde_dots])
+        for ii in range(3):
+            p_axi[ii] = np.linspace(self.coordinates[ii, :].min(), self.coordinates[ii, :].max(), self.kde_dots)
+            d_kde[ii] = gaussian_kde(self.coordinates[ii, :]).evaluate(p_axi[ii, :])
+            d_kde[ii] /= d_kde[ii].max()
+        return d_kde
+
+    def calculate_spectrum(self):
+        """Calculate FFT (in dB) for input 3D coordinates. You can set number of FFT points into the object instance.
 
-    def calculate_fft(self, coordinates: np.ndarray):
-        pass
-        # TODO: Implement this method!
+        """
+        spectrum = fft(self.coordinates, self.fft_dots, axis=0)
+        spectrum = np.abs(fftshift(spectrum, axes=0))
+        # spectrum = np.abs(spectrum)
+        spectrum /= np.max(spectrum)
+        spec_log = 20 * np.log10(spectrum + np.finfo(np.float32).eps)
+        return spec_log
+
+    def calculate_correlation(self):
+        """Calculate auto correlation function for chaotic coordinates.
+
+        """
+        nn, mm = 3, len(self.coordinates)
+        auto_corr = np.zeros([mm, nn])
+        for ii in range(nn):
+            auto_corr[:, ii] = np.correlate(self.coordinates[:, ii], self.coordinates[:, ii], "same")
+        return auto_corr
 
 
 if __name__ == "__main__":
-    pass
+    calc = Calculator()
+
+    num_dots = 10000
+    np.random.seed(42)
+    calc.coordinates = np.random.randn(num_dots, 3)
+
+    dd_kde = calc.check_probability()
+    import matplotlib.pyplot as plt
+
+    print(dd_kde.shape)
+    plt.figure("Probability density function")
+    for idx in range(dd_kde.shape[0]):
+        plt.plot(dd_kde[idx], ".")
+        # plt.xlim([0, dd_kde.shape[1] - 1])
+        plt.grid()
+    plt.tight_layout()
+    plt.show()
diff --git a/src/utils/drawer.py b/src/utils/drawer.py
index 1a84532..ed7f03c 100644
--- a/src/utils/drawer.py
+++ b/src/utils/drawer.py
@@ -56,7 +56,7 @@ class PlotDrawer:
 
     """
 
-    _plot_axis = ((1, 0), (1, 2), (2, 0))
+    _plot_axis = ((0, 1), (2, 1), (0, 2))
     _plot_labels = {0: "X", 1: "Y", 2: "Z"}
 
     def __init__(
@@ -98,10 +98,53 @@ def coordinates(self, value: np.ndarray):
     def min_max_axis(self):
         return np.vstack([np.min(self.coordinates, axis=0), np.max(self.coordinates, axis=0)]).T
 
+    # TODO: Add plots for KDE
+    # # Plot Probability density function
+    # plt.figure("Probability density function")
+    # for ii in range(nn):
+    #     plt.plot(d_kde[ii], ".")
+    #     plt.xlim([0, n_pdf - 1])
+    #     plt.grid()
+    # plt.tight_layout()
+
+    def show_spectrum_and_correlation(self, spectrums: np.ndarray, correlations: np.ndarray):
+        """Plot 3D coordinates as time series."""
+        _ = plt.figure("Autocorrelation and Spectrum", figsize=(8, 6), dpi=100)
+
+        x_corr = np.linspace(-len(self.coordinates) // 2, len(self.coordinates) // 2, len(self.coordinates))
+        x_ffts = np.linspace(-0.5, 0.5, len(spectrums))
+
+        plt.suptitle(f"{self.model_name} attractor", x=0.1)
+        for ii, axis in enumerate(self._plot_labels.values()):
+            plt.subplot(3, 3, ii + 1)
+            plt.title("Time plots", y=1.0) if ii % 2 == 1 else None
+            plt.plot(self.coordinates[:, ii], linewidth=0.75)
+            plt.grid(True)
+            plt.ylabel(axis)
+            plt.xlim([0, len(self.coordinates) - 1])
+            plt.subplot(3, 3, ii + 4)
+            plt.title("Spectrum plots", y=1.0) if ii % 2 == 1 else None
+            plt.plot(x_ffts, spectrums[:, ii], linewidth=0.75)
+            plt.grid(True)
+            plt.ylabel(axis)
+            plt.xlim([np.min(x_ffts), np.max(x_ffts)])
+            plt.subplot(3, 3, ii + 7)
+            plt.title("Correlation plots", y=1.0) if ii % 2 == 1 else None
+            plt.plot(x_corr, correlations[:, ii], linewidth=0.75)
+            plt.grid(True)
+            plt.ylabel(axis)
+            plt.xlim([np.min(x_corr), np.max(x_corr)])
+        plt.tight_layout()
+
+        if self.save_plots:
+            plt.savefig(f"{self.model_name}_spectrum_correlations.png")
+        if self.show_plots:
+            plt.show()
+
     def show_time_plots(self):
         """Plot 3D coordinates as time series."""
         _ = plt.figure("Coordinates evolution in time", figsize=(8, 6), dpi=100)
-        for ii, axis in enumerate(["X", "Y", "Z"]):
+        for ii, axis in enumerate(self._plot_labels.values()):
             plt.subplot(3, 1, ii + 1)
             plt.plot(self.coordinates[:, ii], linewidth=0.75)
             plt.grid(True)
@@ -116,7 +159,7 @@ def show_time_plots(self):
             plt.show()
 
     def _axis_defaults_3d(self, plots):
-        plots.set_title(f"{self.model_name} attractor")
+        plots.set_title(f"{self.model_name} model")
         plots.set_xlabel("X")
         plots.set_ylabel("Y")
         plots.set_zlabel("Z")
@@ -151,12 +194,11 @@ def show_3d_plots(self):
             plt.grid()
             plt.xlabel(self._plot_labels[xx])
             plt.ylabel(self._plot_labels[yy])
-            # TODO: 2020/07/26: Set limits!
             plt.xlim(self.min_max_axis[xx])
             plt.ylim(self.min_max_axis[yy])
 
         ax = fig.add_subplot(2, 2, 4, projection="3d")
-        ax.plot(self.coordinates[:, 0], self.coordinates[:, 1], self.coordinates[:, 2], linewidth=0.7)
+        ax.plot(self.coordinates[:, 0], self.coordinates[:, 1], self.coordinates[:, 2], linewidth=0.75)
         self._axis_defaults_3d(ax)
         plt.tight_layout()
 
@@ -189,7 +231,7 @@ def make_3d_plot_gif(self, step_size: int = 10):
         self._color_map = plt.cm.get_cmap("hsv", step_dots)
 
         ani = animation.FuncAnimation(
-            fig, self.update_coordinates, step_dots, fargs=(step_size,), interval=100, blit=False, repeat=True
+            fig, self._update_coordinates, step_dots, fargs=(step_size,), interval=100, blit=False, repeat=True
         )
 
         if self.save_plots:
@@ -197,7 +239,8 @@ def make_3d_plot_gif(self, step_size: int = 10):
         if self.show_plots:
             plt.show()
 
-    def update_coordinates(self, num, step):
+    def _update_coordinates(self, num, step):
+        """Update plots for making gif"""
         self._plot_list[0].set_data(self.coordinates[0 : 1 + num * step, 0], self.coordinates[0 : 1 + num * step, 1])
         self._plot_list[0].set_3d_properties(self.coordinates[0 : 1 + num * step, 2])
         self._plot_list[0].set_color(self._color_map(num))
@@ -208,31 +251,55 @@ def update_coordinates(self, num, step):
                 )
                 self._plot_list[ii + 1].set_color(self._color_map(num))
 
-    def show_all_plots(self):
-        """Cannot show all plots while 'show_plots' is True.
-        Note: After closing plots you cannot reopen them!
-        """
-        if not self.show_plots:
-            plt.show()
+    # def show_all_plots(self):
+    #     """Cannot show all plots while 'show_plots' is True.
+    #     Note: After closing plots you cannot reopen them!
+    #     """
+    #     if not self.show_plots:
+    #         plt.show()
+    #
+    # @staticmethod
+    # def close_all_plots():
+    #     plt.clf()
+    #     plt.cla()
+    #     plt.close()
+
+
+def random_circle(num_points: int = 100) -> np.ndarray:
+    r"""Create random circle for 3D plots.
+
+    Parameters
+    ----------
+    num_points : int
+        Number of points to draw
+
+    Returns
+    -------
+
+    arr: np.ndarray
+        Numpy 3D array with shapes [3, num_points]
+    """
 
-    @staticmethod
-    def close_all_plots():
-        plt.clf()
-        plt.cla()
-        plt.close()
+    angle = 3
+    ratio = 0.020
+
+    theta = np.linspace(0, angle * np.pi, num_points)
+    xyz = np.vstack([np.cos(theta), np.sin(theta), np.sin(np.linspace(0, 2.05 * np.pi, num_points))]).T
+    xyz += ratio * np.cumsum(np.random.randn(num_points, 3), axis=1)
+    return xyz
 
 
 if __name__ == "__main__":
     np.random.seed(42)
-    points = np.cumsum(np.random.randn(50, 3), axis=1)
+
+    circle_points = random_circle(num_points=200)
 
     drawer = PlotDrawer(show_plots=True, add_2d_gif=True)
-    drawer.coordinates = points
+    drawer.coordinates = circle_points
     drawer.model_name = "Chaotic"
 
     # print(drawer.min_max_axis)
-    drawer.make_3d_plot_gif()
-    # drawer.make_3d_plot_gif(show_2d_plots=True)
+    drawer.make_3d_plot_gif(step_size=20)
     # drawer.show_time_plots()
     # drawer.show_3d_plots()
     # drawer.show_all_plots()