02_12_2023 Implementation of empirical mode decomposition.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9kEVwMCf8o1g"
   },
   "source": [
    "A simpel sine wave is an IMF. An IMF have a zero mean.\n",
    "\n",
    "*   List item\n",
    "*   List item\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 34
    },
    "id": "gEks1f569sb5",
    "outputId": "907f6c3c-fa11-4688-b277-e870d48db76b"
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "import scipy.io\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#setting\n",
    "minus_step = 9\n",
    "frequency = 1\n",
    "total_time = 1\n",
    "res = 0.01\n",
    "\n",
    "time = np.arange(0, total_time, res)\n",
    "amplitude = np.sin(2 * np.pi * frequency * total_time \n",
    "                   * (time/time[-1]))\n",
    "np.mean(amplitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3-LpuzKQ_D5M"
   },
   "source": [
    "The difference between minima and maxima of IMF is at most by 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 315
    },
    "id": "Pz_9GQmUA0Fp",
    "outputId": "cf8e50e8-ec79-4758-9e1d-e28a256bd213"
   },
   "outputs": [],
   "source": [
    "nodata_flag = np.nan\n",
    "envelope_max = np.zeros(len(amplitude))\n",
    "envelope_min = np.zeros(len(amplitude))\n",
    "envelope_max[:] = nodata_flag\n",
    "envelope_min[:] = nodata_flag\n",
    "\n",
    "minima_total = 0\n",
    "maxima_total = 0\n",
    "graph_stat_history = 'up'\n",
    "graph_stat_recent = 'up'\n",
    "for i in range(len(amplitude)):\n",
    "  if i == 1:\n",
    "    if amplitude[i] > amplitude[i-1]:\n",
    "      graph_stat_recent = 'up'\n",
    "      graph_stat_history = 'up'\n",
    "    elif amplitude[i] < amplitude[i-1]:\n",
    "      graph_stat_recent = 'down'\n",
    "      graph_stat_history = 'down'\n",
    "    else:\n",
    "      graph_stat_recent = graph_stat_history\n",
    "    continue\n",
    "  if amplitude[i] > amplitude[i-1]:\n",
    "    graph_stat_recent = 'up'\n",
    "  elif amplitude[i] < amplitude[i-1]:\n",
    "    graph_stat_recent = 'down'\n",
    "  else:\n",
    "    graph_stat_recent = graph_stat_history\n",
    "\t\n",
    "  if graph_stat_recent != graph_stat_history:\n",
    "    if graph_stat_recent == 'up':\n",
    "      envelope_min[i] = amplitude[i-1]\n",
    "      minima_total += 1\n",
    "      graph_stat_history = 'up'\n",
    "    elif graph_stat_recent == 'down':\n",
    "      envelope_max[i] = amplitude[i-1]\n",
    "      maxima_total += 1\n",
    "      graph_stat_history = 'down'\n",
    "\n",
    "print('Total minima', minima_total)\n",
    "print('Total maxima', maxima_total)\n",
    "plt.title('A Sine Wave (Red dot = maxima, Blue dot = minima)')\n",
    "plt.plot(time, amplitude)\n",
    "plt.plot(time, envelope_min, 'bo')\n",
    "plt.plot(time, envelope_max, 'ro')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZzyVMbuKCus4"
   },
   "source": [
    "In this article, I will show you the explanation of EMD algorithm that easy to understand (using practical code and many graph plot). First, let's decide our toy test signal that will be decomposed using EMD. In this article, I will merge 1 Hz and 4 Hz sine wave. It will be like this.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 809
    },
    "id": "rLxksZHVH78v",
    "outputId": "bb7bf7a8-314c-4d67-ba39-0c0c42e00641"
   },
   "outputs": [],
   "source": [
    "frequency_1 = 1\n",
    "frequency_2 = 4\n",
    "total_time_test = 5\n",
    "res_test = 0.01\n",
    "\n",
    "time_test = np.arange(0, total_time_test, res_test)\n",
    "amplitude_1 = np.sin(2 * np.pi * frequency_1 * total_time_test \\\n",
    "                     * (time_test/time_test[-1]))\n",
    "amplitude_2 = np.sin(2 * np.pi * frequency_2 * total_time_test \\\n",
    "                     * (time_test/time_test[-1]))\n",
    "amplitude_fusion = amplitude_1 + amplitude_2\n",
    "plt.plot(time_test, amplitude_1)\n",
    "plt.title('%s Hz Wave'%(frequency_1))\n",
    "plt.show()\n",
    "plt.close()\n",
    "plt.plot(time_test, amplitude_2)\n",
    "plt.title('%s Hz Wave'%(frequency_2))\n",
    "plt.show()\n",
    "plt.close()\n",
    "plt.plot(time_test, amplitude_fusion)\n",
    "plt.title('A Fusion Wave')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fbzjPq8vKgn2"
   },
   "source": [
    "**Step 1, Find the Minima and Maxima**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 332
    },
    "id": "Cu19PGvdLIlT",
    "outputId": "d7fec245-7ce7-4849-f882-8e6487a88f61"
   },
   "outputs": [],
   "source": [
    "nodata_flag = np.nan\n",
    "envelope_max_test = np.zeros(len(amplitude_fusion))\n",
    "envelope_min_test = np.zeros(len(amplitude_fusion))\n",
    "envelope_max_test[:] = nodata_flag\n",
    "envelope_min_test[:] = nodata_flag\n",
    "\n",
    "minima_total = 0\n",
    "maxima_total = 0\n",
    "graph_stat_history = 'up'\n",
    "graph_stat_recent = 'up'\n",
    "for i in range(len(amplitude_fusion)):\n",
    "  if i == 1:\n",
    "    if amplitude_fusion[i] > amplitude_fusion[i-1]:\n",
    "      graph_stat_recent = 'up'\n",
    "      graph_stat_history = 'up'\n",
    "    elif amplitude_fusion[i] < amplitude_fusion[i-1]:\n",
    "      graph_stat_recent = 'down'\n",
    "      graph_stat_history = 'down'\n",
    "    else:\n",
    "      graph_stat_recent = graph_stat_history\n",
    "    continue\n",
    "  if amplitude_fusion[i] > amplitude_fusion[i-1]:\n",
    "    graph_stat_recent = 'up'\n",
    "  elif amplitude_fusion[i] < amplitude_fusion[i-1]:\n",
    "    graph_stat_recent = 'down'\n",
    "  else:\n",
    "    graph_stat_recent = graph_stat_history\n",
    "\t\n",
    "  if graph_stat_recent != graph_stat_history:\n",
    "    if graph_stat_recent == 'up':\n",
    "      envelope_min_test[i] = amplitude_fusion[i-1]\n",
    "      minima_total += 1\n",
    "      graph_stat_history = 'up'\n",
    "    elif graph_stat_recent == 'down':\n",
    "      envelope_max_test[i] = amplitude_fusion[i-1]\n",
    "      maxima_total += 1\n",
    "      graph_stat_history = 'down'\n",
    "\n",
    "print(maxima_total)\n",
    "print(minima_total)\n",
    "plt.title('A Fusion Wave')\n",
    "plt.plot(time_test, amplitude_fusion)\n",
    "plt.plot(time_test, envelope_min_test, 'bo')\n",
    "plt.plot(time_test, envelope_max_test, 'ro')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "f-4Kk3VSY5f4"
   },
   "source": [
    "**Step 2, create an envelope of minima and maxima from the array of minima and maxima**. Interpolate the value of minima and maxima to create an envelope of minima and maxima using a cubic spline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 298
    },
    "id": "zcSDOkVQcM3-",
    "outputId": "b9802b2e-1196-4585-c796-139fe94e3a4c"
   },
   "outputs": [],
   "source": [
    "from scipy.interpolate import CubicSpline\n",
    "\n",
    "all_minima_point = []\n",
    "all_minima_point_time = []\n",
    "all_maxima_point = []\n",
    "all_maxima_point_time = []\n",
    "for i in range(len(envelope_min_test)):\n",
    "  if not np.isnan(envelope_min_test[i]):\n",
    "    all_minima_point.append(envelope_min_test[i])\n",
    "    all_minima_point_time.append(time_test[i])\n",
    "\n",
    "for i in range(len(envelope_max_test)):\n",
    "  if not np.isnan(envelope_max_test[i]):\n",
    "    all_maxima_point.append(envelope_max_test[i])\n",
    "    all_maxima_point_time.append(time_test[i])\n",
    "\n",
    "all_minima_point = np.asarray(all_minima_point)\n",
    "all_minima_point_time = np.asarray(all_minima_point_time)\n",
    "all_maxima_point = np.asarray(all_maxima_point)\n",
    "all_maxima_point_time = np.asarray(all_maxima_point_time)\n",
    "\n",
    "#print(all_maxima_point, all_minima_point)\n",
    "cs_minima = CubicSpline(all_minima_point_time, \\\n",
    "                        all_minima_point)\n",
    "cs_maxima = CubicSpline(all_maxima_point_time, \\\n",
    "                        all_maxima_point)\n",
    "#print(cs_minima, cs_maxima)\n",
    "envelope_minima_cs = cs_minima(time_test)\n",
    "envelope_maxima_cs = cs_maxima(time_test)\n",
    "#print(envelope_minima_cs)\n",
    "\n",
    "plt.title('A Fusion Wave')\n",
    "plt.plot(time_test, amplitude_fusion)\n",
    "plt.plot(time_test, envelope_min_test, 'bo')\n",
    "plt.plot(time_test, envelope_max_test, 'ro')\n",
    "plt.plot(time_test, envelope_minima_cs)\n",
    "plt.plot(time_test, envelope_maxima_cs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WbbjScXT2Ejf"
   },
   "source": [
    "**Step 4, from the envelope of minima and maxima, get the middle value**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 298
    },
    "id": "sQn6rsxH3ILJ",
    "outputId": "665a9f75-e03b-441f-fe6c-aed3cb22627a"
   },
   "outputs": [],
   "source": [
    "mid_value = np.zeros(len(envelope_maxima_cs))\n",
    "for i in range(len(mid_value)):\n",
    "  mid_value[i] = (envelope_maxima_cs[i] + \\\n",
    "                  envelope_minima_cs[i]) / 2\n",
    "\n",
    "plt.title('A Fusion Wave')\n",
    "plt.plot(time_test, amplitude_fusion)\n",
    "plt.plot(time_test, envelope_min_test, 'bo')\n",
    "plt.plot(time_test, envelope_max_test, 'ro')\n",
    "plt.plot(time_test, envelope_minima_cs)\n",
    "plt.plot(time_test, envelope_maxima_cs)\n",
    "plt.plot(time_test, mid_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-XttMZ0H9H_5"
   },
   "source": [
    "**Step 5, decrease the value of the real toy test signal by the middle value of the envelope**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 298
    },
    "id": "ar61Ph0f-jMh",
    "outputId": "6a8c08ea-6790-451b-f70c-98539a88e816"
   },
   "outputs": [],
   "source": [
    "extracted_val = amplitude_fusion - mid_value\n",
    "plt.title('An Extracted Wave')\n",
    "plt.plot(time_test, extracted_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "knsAq-XmFO1X"
   },
   "source": [
    "If you look closely, this is our 4 Hz sine wave successfully extracted. But looks like there is some error like a sudden leap on both sides of the wave. **This has come from the error of the interpolation of the cubic spline. The interpolation of the cubic spline becomes so bad at the beginning and the end of the signal.** The solution is we must define a boundary condition. What I mean about boundary condition is **we must define the envelope of minima and maxima at the beginning and the end of the signal before doing the cubic spline**. In this article, the envelope of minima and maxima of the beginning and the end of the signal is the same as with the nearest maxima and minima from them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 545
    },
    "id": "eYpWfj9xUHDB",
    "outputId": "5260d48e-2526-4885-e78c-400b552ceff8"
   },
   "outputs": [],
   "source": [
    "nodata_flag = np.nan\n",
    "envelope_max_test = np.zeros(len(amplitude_fusion))\n",
    "envelope_min_test = np.zeros(len(amplitude_fusion))\n",
    "envelope_max_test[:] = nodata_flag\n",
    "envelope_min_test[:] = nodata_flag\n",
    "\n",
    "minima_total = 0\n",
    "maxima_total = 0\n",
    "graph_stat_history = 'up'\n",
    "graph_stat_recent = 'up'\n",
    "for i in range(len(amplitude_fusion)):\n",
    "  if i == 1:\n",
    "    if amplitude_fusion[i] > amplitude_fusion[i-1]:\n",
    "      graph_stat_recent = 'up'\n",
    "      graph_stat_history = 'up'\n",
    "    elif amplitude_fusion[i] < amplitude_fusion[i-1]:\n",
    "      graph_stat_recent = 'down'\n",
    "      graph_stat_history = 'down'\n",
    "    else:\n",
    "      graph_stat_recent = graph_stat_history\n",
    "    continue\n",
    "  if amplitude_fusion[i] > amplitude_fusion[i-1]:\n",
    "    graph_stat_recent = 'up'\n",
    "  elif amplitude_fusion[i] < amplitude_fusion[i-1]:\n",
    "    graph_stat_recent = 'down'\n",
    "  else:\n",
    "    graph_stat_recent = graph_stat_history\n",
    "\t\n",
    "  if graph_stat_recent != graph_stat_history:\n",
    "    if graph_stat_recent == 'up':\n",
    "      envelope_min_test[i] = amplitude_fusion[i-1]\n",
    "      minima_total += 1\n",
    "      graph_stat_history = 'up'\n",
    "    elif graph_stat_recent == 'down':\n",
    "      envelope_max_test[i] = amplitude_fusion[i-1]\n",
    "      maxima_total += 1\n",
    "      graph_stat_history = 'down'\n",
    "\n",
    "#boundary condition of minima\n",
    "for i in range(len(envelope_min_test)):\n",
    "  if not np.isnan(envelope_min_test[i]):\n",
    "    envelope_min_test[0] = envelope_min_test[i]\n",
    "    break\n",
    "envelope_min_test_reverse = envelope_min_test[::-1]\n",
    "for i in range(len(envelope_min_test_reverse)):\n",
    "  if not np.isnan(envelope_min_test_reverse[i]):\n",
    "    envelope_min_test[-1] = envelope_min_test_reverse[i]\n",
    "    break\n",
    "\n",
    "#boundary condition of maxima\n",
    "for i in range(len(envelope_max_test)):\n",
    "  if not np.isnan(envelope_max_test[i]):\n",
    "    envelope_max_test[0] = envelope_max_test[i]\n",
    "    break\n",
    "envelope_max_test_reverse = envelope_max_test[::-1]\n",
    "for i in range(len(envelope_max_test_reverse)):\n",
    "  if not np.isnan(envelope_max_test_reverse[i]):\n",
    "    envelope_max_test[-1] = envelope_max_test_reverse[i]\n",
    "    break\n",
    "\n",
    "all_minima_point = []\n",
    "all_minima_point_time = []\n",
    "all_maxima_point = []\n",
    "all_maxima_point_time = []\n",
    "for i in range(len(envelope_min_test)):\n",
    "  if not np.isnan(envelope_min_test[i]):\n",
    "    all_minima_point.append(envelope_min_test[i])\n",
    "    all_minima_point_time.append(time_test[i])\n",
    "\n",
    "for i in range(len(envelope_max_test)):\n",
    "  if not np.isnan(envelope_max_test[i]):\n",
    "    all_maxima_point.append(envelope_max_test[i])\n",
    "    all_maxima_point_time.append(time_test[i])\n",
    "\n",
    "all_minima_point = np.asarray(all_minima_point)\n",
    "all_minima_point_time = np.asarray(all_minima_point_time)\n",
    "all_maxima_point = np.asarray(all_maxima_point)\n",
    "all_maxima_point_time = np.asarray(all_maxima_point_time)\n",
    "\n",
    "#print(all_maxima_point, all_minima_point)\n",
    "cs_minima = CubicSpline(all_minima_point_time, \n",
    "                        all_minima_point)\n",
    "cs_maxima = CubicSpline(all_maxima_point_time, \n",
    "                        all_maxima_point)\n",
    "#print(cs_minima, cs_maxima)\n",
    "envelope_minima_cs = cs_minima(time_test)\n",
    "envelope_maxima_cs = cs_maxima(time_test)\n",
    "#print(envelope_minima_cs)\n",
    "\n",
    "mid_value = np.zeros(len(envelope_maxima_cs))\n",
    "for i in range(len(mid_value)):\n",
    "  mid_value[i] = (envelope_maxima_cs[i] + \n",
    "                  envelope_minima_cs[i])/2\n",
    "\n",
    "plt.title('A Fusion Wave')\n",
    "plt.plot(time_test, amplitude_fusion)\n",
    "plt.plot(time_test, envelope_min_test, 'bo')\n",
    "plt.plot(time_test, envelope_max_test, 'ro')\n",
    "plt.plot(time_test, envelope_minima_cs)\n",
    "plt.plot(time_test, envelope_maxima_cs)\n",
    "plt.plot(time_test, mid_value)\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "extracted_val = amplitude_fusion - mid_value\n",
    "plt.title('An Extracted Wave')\n",
    "plt.plot(time_test, extracted_val)\n",
    "plt.show()\n",
    "plt.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lrujikRxCYXD"
   },
   "source": [
    "Ouch, this time the error is affected the entire extracted signal. But the sudden leap has vanished.\n",
    "\n",
    "**Step 6, Check if this extracted signal is an IMF**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 68
    },
    "id": "rppS5mN5CtmW",
    "outputId": "aae545b2-8ad9-43c7-f2f4-2257d2a942a6"
   },
   "outputs": [],
   "source": [
    "np.mean(extracted_val)\n",
    "\n",
    "minima_total = 0\n",
    "maxima_total = 0\n",
    "graph_stat_history = 'up'\n",
    "graph_stat_recent = 'up'\n",
    "for i in range(len(extracted_val)):\n",
    "  if i == 1:\n",
    "    if extracted_val[i] > extracted_val[i-1]:\n",
    "      graph_stat_recent = 'up'\n",
    "      graph_stat_history = 'up'\n",
    "    elif extracted_val[i] < extracted_val[i-1]:\n",
    "      graph_stat_recent = 'down'\n",
    "      graph_stat_history = 'down'\n",
    "    else:\n",
    "      graph_stat_recent = graph_stat_history\n",
    "    continue\n",
    "  if extracted_val[i] > extracted_val[i-1]:\n",
    "    graph_stat_recent = 'up'\n",
    "  elif extracted_val[i] < extracted_val[i-1]:\n",
    "    graph_stat_recent = 'down'\n",
    "  else:\n",
    "    graph_stat_recent = graph_stat_history\n",
    "\t\n",
    "  if graph_stat_recent != graph_stat_history:\n",
    "    if graph_stat_recent == 'up':\n",
    "      minima_total += 1\n",
    "      graph_stat_history = 'up'\n",
    "    elif graph_stat_recent == 'down':\n",
    "      maxima_total += 1\n",
    "      graph_stat_history = 'down'\n",
    "\n",
    "print('mean:', np.mean(extracted_val))\n",
    "print('Total minima', minima_total)\n",
    "print('Total maxima', maxima_total)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "z_3WJHPdGI7I"
   },
   "source": [
    "The mean is nearly zero, I think we can ignore it and round it to zero. The number of Minima and Maxima also fulfill the requirement, **we get an IMF**. If you are not satisfied with the mean of the extracted signal, you can process it 1 more time following step 2 - step 6. But **this time, instead of inputting a real toy test signal, we use that extracted signal**. This the route of the algorithm when in the first or second or so on triest, the extracted signal not satisfied IMF condition. To make it more simple, let's just admit the signal is an IMF.\n",
    "\n",
    "**Step 7, decrease the original toy test signal with this IMF.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 281
    },
    "id": "9xDe65wsFI7M",
    "outputId": "96858ddc-b0ab-461d-c18c-138de548171f"
   },
   "outputs": [],
   "source": [
    "amplitude_fusion = amplitude_fusion - extracted_val\n",
    "\n",
    "plt.title('A Fusion Wave')\n",
    "plt.plot(time_test, amplitude_fusion)\n",
    "plt.show()\n",
    "plt.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "joEZiwgTFvNn"
   },
   "source": [
    "Woah, we get our 1 Hz signal (but as we thought, the beginning and the end is wrong). It looks like we success to decompose our fused signal. But, the condition to stop in the algorithm of EMD is when our residuum signal is just a constant, monotonic, or just have 1 extremum (just have 1 minima or maxima). **Our signal which we have reduced by the IMF before is called residuum**. Let's continue the process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 860
    },
    "id": "uxg8YhhgIK9n",
    "outputId": "78068d84-33a6-4542-b193-e57f0306b3af"
   },
   "outputs": [],
   "source": [
    "envelope_max_test = np.zeros(len(amplitude_fusion))\n",
    "envelope_min_test = np.zeros(len(amplitude_fusion))\n",
    "envelope_max_test[:] = nodata_flag\n",
    "envelope_min_test[:] = nodata_flag\n",
    "\n",
    "minima_total = 0\n",
    "maxima_total = 0\n",
    "graph_stat_history = 'up'\n",
    "graph_stat_recent = 'up'\n",
    "for i in range(len(amplitude_fusion)):\n",
    "  if i == 1:\n",
    "    if amplitude_fusion[i] > amplitude_fusion[i-1]:\n",
    "      graph_stat_recent = 'up'\n",
    "      graph_stat_history = 'up'\n",
    "    elif amplitude_fusion[i] < amplitude_fusion[i-1]:\n",
    "      graph_stat_recent = 'down'\n",
    "      graph_stat_history = 'down'\n",
    "    else:\n",
    "      graph_stat_recent = graph_stat_history\n",
    "    continue\n",
    "  if amplitude_fusion[i] > amplitude_fusion[i-1]:\n",
    "    graph_stat_recent = 'up'\n",
    "  elif amplitude_fusion[i] < amplitude_fusion[i-1]:\n",
    "    graph_stat_recent = 'down'\n",
    "  else:\n",
    "    graph_stat_recent = graph_stat_history\n",
    "\t\n",
    "  if graph_stat_recent != graph_stat_history:\n",
    "    if graph_stat_recent == 'up':\n",
    "      envelope_min_test[i] = amplitude_fusion[i-1]\n",
    "      minima_total += 1\n",
    "      graph_stat_history = 'up'\n",
    "    elif graph_stat_recent == 'down':\n",
    "      envelope_max_test[i] = amplitude_fusion[i-1]\n",
    "      maxima_total += 1\n",
    "      graph_stat_history = 'down'\n",
    "\n",
    "all_minima_point = []\n",
    "all_minima_point_time = []\n",
    "all_maxima_point = []\n",
    "all_maxima_point_time = []\n",
    "for i in range(len(envelope_min_test)):\n",
    "  if not np.isnan(envelope_min_test[i]):\n",
    "    all_minima_point.append(envelope_min_test[i])\n",
    "    all_minima_point_time.append(time_test[i])\n",
    "\n",
    "for i in range(len(envelope_max_test)):\n",
    "  if not np.isnan(envelope_max_test[i]):\n",
    "    all_maxima_point.append(envelope_max_test[i])\n",
    "    all_maxima_point_time.append(time_test[i])\n",
    "\n",
    "all_minima_point = np.asarray(all_minima_point)\n",
    "all_minima_point_time = np.asarray(all_minima_point_time)\n",
    "all_maxima_point = np.asarray(all_maxima_point)\n",
    "all_maxima_point_time = np.asarray(all_maxima_point_time)\n",
    "\n",
    "#boundary condition of minima\n",
    "for i in range(len(envelope_min_test)):\n",
    "  if not np.isnan(envelope_min_test[i]):\n",
    "    envelope_min_test[0] = envelope_min_test[i]\n",
    "    break\n",
    "envelope_min_test_reverse = envelope_min_test[::-1]\n",
    "for i in range(len(envelope_min_test_reverse)):\n",
    "  if not np.isnan(envelope_min_test_reverse[i]):\n",
    "    envelope_min_test[-1] = envelope_min_test_reverse[i]\n",
    "    break\n",
    "\n",
    "#boundary condition of maxima\n",
    "for i in range(len(envelope_max_test)):\n",
    "  if not np.isnan(envelope_max_test[i]):\n",
    "    envelope_max_test[0] = envelope_max_test[i]\n",
    "    break\n",
    "envelope_max_test_reverse = envelope_max_test[::-1]\n",
    "for i in range(len(envelope_max_test_reverse)):\n",
    "  if not np.isnan(envelope_max_test_reverse[i]):\n",
    "    envelope_max_test[-1] = envelope_max_test_reverse[i]\n",
    "    break\n",
    "\n",
    "#print(all_maxima_point, all_minima_point)\n",
    "cs_minima = CubicSpline(all_minima_point_time, \n",
    "                        all_minima_point)\n",
    "cs_maxima = CubicSpline(all_maxima_point_time, \n",
    "                        all_maxima_point)\n",
    "#print(cs_minima, cs_maxima)\n",
    "envelope_minima_cs = cs_minima(time_test)\n",
    "envelope_maxima_cs = cs_maxima(time_test)\n",
    "#print(envelope_minima_cs)\n",
    "\n",
    "mid_value = np.zeros(len(envelope_maxima_cs))\n",
    "for i in range(len(mid_value)):\n",
    "  mid_value[i] = (envelope_maxima_cs[i] + \n",
    "                  envelope_minima_cs[i]) / 2\n",
    "\n",
    "plt.title('A Fusion Wave')\n",
    "plt.plot(time_test, amplitude_fusion)\n",
    "plt.plot(time_test, envelope_min_test, 'bo')\n",
    "plt.plot(time_test, envelope_max_test, 'ro')\n",
    "plt.plot(time_test, envelope_minima_cs)\n",
    "plt.plot(time_test, envelope_maxima_cs)\n",
    "plt.plot(time_test, mid_value)\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "extracted_val_1 = amplitude_fusion - mid_value\n",
    "plt.title('A Fusion Wave')\n",
    "plt.plot(time_test, extracted_val_1)\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "print('mean:', np.mean(extracted_val))\n",
    "print('Total minima', minima_total)\n",
    "print('Total maxima', maxima_total)\n",
    "\n",
    "amplitude_fusion_1 = amplitude_fusion - extracted_val_1\n",
    "\n",
    "plt.title('A Residuum')\n",
    "plt.plot(time_test, amplitude_fusion_1)\n",
    "plt.show()\n",
    "plt.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nGbZds2egQRZ"
   },
   "source": [
    "It looks like the requirement of the residuum not yet fulfilled. Let's processed it again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 860
    },
    "id": "8iFIOq5tgnww",
    "outputId": "b46e5395-41aa-4aaf-ef6b-5dfdeff30c54"
   },
   "outputs": [],
   "source": [
    "envelope_max_test = np.zeros(len(amplitude_fusion_1))\n",
    "envelope_min_test = np.zeros(len(amplitude_fusion_1))\n",
    "envelope_max_test[:] = nodata_flag\n",
    "envelope_min_test[:] = nodata_flag\n",
    "\n",
    "minima_total = 0\n",
    "maxima_total = 0\n",
    "graph_stat_history = 'up'\n",
    "graph_stat_recent = 'up'\n",
    "for i in range(len(amplitude_fusion_1)):\n",
    "  if i == 1:\n",
    "    if amplitude_fusion_1[i] > amplitude_fusion_1[i-1]:\n",
    "      graph_stat_recent = 'up'\n",
    "      graph_stat_history = 'up'\n",
    "    elif amplitude_fusion_1[i] < amplitude_fusion_1[i-1]:\n",
    "      graph_stat_recent = 'down'\n",
    "      graph_stat_history = 'down'\n",
    "    else:\n",
    "      graph_stat_recent = graph_stat_history\n",
    "    continue\n",
    "  if amplitude_fusion_1[i] > amplitude_fusion_1[i-1]:\n",
    "    graph_stat_recent = 'up'\n",
    "  elif amplitude_fusion_1[i] < amplitude_fusion_1[i-1]:\n",
    "    graph_stat_recent = 'down'\n",
    "  else:\n",
    "    graph_stat_recent = graph_stat_history\n",
    "\t\n",
    "  if graph_stat_recent != graph_stat_history:\n",
    "    if graph_stat_recent == 'up':\n",
    "      envelope_min_test[i] = amplitude_fusion_1[i-1]\n",
    "      minima_total += 1\n",
    "      graph_stat_history = 'up'\n",
    "    elif graph_stat_recent == 'down':\n",
    "      envelope_max_test[i] = amplitude_fusion_1[i-1]\n",
    "      maxima_total += 1\n",
    "      graph_stat_history = 'down'\n",
    "\n",
    "all_minima_point = []\n",
    "all_minima_point_time = []\n",
    "all_maxima_point = []\n",
    "all_maxima_point_time = []\n",
    "for i in range(len(envelope_min_test)):\n",
    "  if not np.isnan(envelope_min_test[i]):\n",
    "    all_minima_point.append(envelope_min_test[i])\n",
    "    all_minima_point_time.append(time_test[i])\n",
    "\n",
    "for i in range(len(envelope_max_test)):\n",
    "  if not np.isnan(envelope_max_test[i]):\n",
    "    all_maxima_point.append(envelope_max_test[i])\n",
    "    all_maxima_point_time.append(time_test[i])\n",
    "\n",
    "all_minima_point = np.asarray(all_minima_point)\n",
    "all_minima_point_time = np.asarray(all_minima_point_time)\n",
    "all_maxima_point = np.asarray(all_maxima_point)\n",
    "all_maxima_point_time = np.asarray(all_maxima_point_time)\n",
    "\n",
    "#boundary condition of minima\n",
    "for i in range(len(envelope_min_test)):\n",
    "  if not np.isnan(envelope_min_test[i]):\n",
    "    envelope_min_test[0] = envelope_min_test[i]\n",
    "    break\n",
    "envelope_min_test_reverse = envelope_min_test[::-1]\n",
    "for i in range(len(envelope_min_test_reverse)):\n",
    "  if not np.isnan(envelope_min_test_reverse[i]):\n",
    "    envelope_min_test[-1] = envelope_min_test_reverse[i]\n",
    "    break\n",
    "\n",
    "#boundary condition of maxima\n",
    "for i in range(len(envelope_max_test)):\n",
    "  if not np.isnan(envelope_max_test[i]):\n",
    "    envelope_max_test[0] = envelope_max_test[i]\n",
    "    break\n",
    "envelope_max_test_reverse = envelope_max_test[::-1]\n",
    "for i in range(len(envelope_max_test_reverse)):\n",
    "  if not np.isnan(envelope_max_test_reverse[i]):\n",
    "    envelope_max_test[-1] = envelope_max_test_reverse[i]\n",
    "    break\n",
    "\n",
    "#print(all_maxima_point, all_minima_point)\n",
    "cs_minima = CubicSpline(all_minima_point_time, all_minima_point)\n",
    "cs_maxima = CubicSpline(all_maxima_point_time, all_maxima_point)\n",
    "#print(cs_minima, cs_maxima)\n",
    "envelope_minima_cs = cs_minima(time_test)\n",
    "envelope_maxima_cs = cs_maxima(time_test)\n",
    "#print(envelope_minima_cs)\n",
    "\n",
    "mid_value = np.zeros(len(envelope_maxima_cs))\n",
    "for i in range(len(mid_value)):\n",
    "  mid_value[i] = (envelope_maxima_cs[i] + envelope_minima_cs[i]) / 2\n",
    "\n",
    "plt.title('A Fusion Wave')\n",
    "plt.plot(time_test, amplitude_fusion_1)\n",
    "plt.plot(time_test, envelope_min_test, 'bo')\n",
    "plt.plot(time_test, envelope_max_test, 'ro')\n",
    "plt.plot(time_test, envelope_minima_cs)\n",
    "plt.plot(time_test, envelope_maxima_cs)\n",
    "plt.plot(time_test, mid_value)\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "extracted_val_final = amplitude_fusion_1 - mid_value\n",
    "plt.title('A Fusion Wave')\n",
    "plt.plot(time_test, extracted_val_final)\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "print('mean:', np.mean(extracted_val))\n",
    "print('Total minima', minima_total)\n",
    "print('Total maxima', maxima_total)\n",
    "\n",
    "amplitude_fusion_final = amplitude_fusion_1 - extracted_val_final\n",
    "\n",
    "plt.title('A Residuum')\n",
    "plt.plot(time_test, amplitude_fusion_final)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5IFQVIvTh3OC"
   },
   "source": [
    "And at last, we find the residuum (monotonic). The last extracted looks like just noise because of the accumulative error of boundary. **This is can be prevented by raising the threshold of the mean of the signal**. When the mean of the signal is not fulfilling the threshold (the extracted signal is not an IMF), it will try again following step 2 - 6 until the extracted signal fulfill the threshold. The more convergence the threshold to zero, the more likely this noise vanishes. So that this article is not too long, in the EMD algorithm demonstration earlier, the threshold we used was not too strict. Because in the end what we want to achieve here is that readers understand step by step the EMD algorithm easily. See you."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}