diff --git a/src/detrend.py b/src/detrend.py index 9de8111..89eb2d4 100644 --- a/src/detrend.py +++ b/src/detrend.py @@ -35,6 +35,7 @@ import pandas as pd import matplotlib.pyplot as plt from smoothingspline import spline +import emdecomp as emd from autoreg import ar_func import curvefit @@ -45,8 +46,17 @@ # differences def detrend(data, fit="spline", method="residual"): nullremoved_data = data.dropna() - yi = spline(data.dropna()) - residual(nullremoved_data, yi) + if fit=="spline": + yi = spline(nullremoved_data) + residual(nullremoved_data, yi) + + #Empirical Mode Decomposition (EMD) + elif fit=="emd": + imfs = emd.emd(nullremoved_data) + imfs_p = emd.phase_spectrum(imfs) + mis = emd.phase_mi(imfs_p) + stochastic_component, yi = emd.divide_signal(nullremoved_data, imfs, mis, cutoff=0.5) + residual(nullremoved_data, yi) # Detrends by finding ratio of original series data to curve data def residual(series, yi): diff --git a/src/emdecomp.py b/src/emdecomp.py new file mode 100644 index 0000000..1d345cb --- /dev/null +++ b/src/emdecomp.py @@ -0,0 +1,96 @@ +__copyright__ = """ + dplPy for tree ring width time series analyses + Copyright (C) 2022 OpenDendro + + 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. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +""" + +__license__ = "GNU GPLv3" + +#!/usr/bin/python +# -*- coding: utf-8 -*- + +# Date: 5/27/2022 +# Author: Grigoriy Lozhkin +# Title: emdecomp.py +# Description: This contains the spline method which fits +# a series to a spline curve. +# example usage (in other file): +# from emdecomp import emd +# yi = emd(series) + + +import numpy as np +from scipy.fft import fft +from PyEMD import EMD +from sktime.datasets import load_airline, load_shampoo_sales +from sklearn.feature_selection import mutual_info_regression +import matplotlib.pyplot as plt + +def emd(signal): + emd = EMD() + imfs = emd(signal.values) + return imfs + + +def phase_spectrum(imfs): + imfs_p = [] + for i, imf in enumerate(imfs): + trans = fft(imf) + imf_p = np.arctan(trans.imag / trans.real) + imfs_p.append(imf_p) + + return imfs_p + +def phase_mi(phases): + mis = [] + for i in range(len(phases)-1): + mis.append(mutual_info_regression(phases[i].reshape(-1, 1), phases[i+1])[0]) + + return np.array(mis) + +def divide_signal(signal, imfs, mis, cutoff=0.05): + + x = signal.index.to_numpy() + y = signal.to_numpy() + + cut_point = np.where(mis > cutoff)[0][0] + stochastic_component = np.sum(imfs[:cut_point], axis=0) + deterministic_component = np.sum(imfs[cut_point:], axis=0) + + t = [i for i in range(len(signal))] + + fig, axs = plt.subplots(3, 1, figsize=(15,8)) + axs[0].plot(t, signal.values) + axs[0].set_title('Original Signal') + + axs[1].plot(t, stochastic_component) + axs[1].set_title('Stochastic Component') + + axs[2].plot(t, deterministic_component) + axs[2].set_title('Deterministic Component') + + plt.show() + + plt.plot(x, y, "o", x, deterministic_component, "-") + plt.show() + return stochastic_component, deterministic_component + + + + + + + +