2_FFT&DFT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "#  Feature extraction using Best Fourier Coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook was provided as part of [2_Machine_Learning_using_Best_Fourier_Coefficients.ipynb](2_Machine_Learning_using_Best_Fourier_Coefficients.ipynb) which contains some functions that take some time for execution. This notebook can give better insights in outputs and is easier to follow. For details, please go to:\n",
    "[2_Machine_Learning_using_Best_Fourier_Coefficients.ipynb](2_Machine_Learning_using_Best_Fourier_Coefficients.ipynb) and execute all cells.\n",
    "\n",
    "\n",
    "*Before usage of this notebook, please download folder `Files_for_notebooks.zip` from this link https://github.com/harislulic/ZeMA-machine-learning-tutorials/releases/tag/v0.1.2\n",
    "and store the files in the same folder which is location for this notebook.\n",
    "It is also necessary to install pip in your environment and using this, package PyDynamic:*\n",
    "\n",
    "pip install PyDynamic\n",
    "\n",
    "*In order to see interactive diagrams, write:* \n",
    "\n",
    "pip install ipywidgets\n",
    "\n",
    "jupyter nbextension enable --py widgetsnbextension   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### Importing the data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jugo01\\AppData\\Local\\Continuum\\anaconda3\\envs\\Zema\\lib\\site-packages\\PyDynamic\\identification\\fit_filter.py:34: DeprecationWarning: The module *identification* will be combined with the module *deconvolution* and renamed to *model_estimation* in the next major release 3.0. From then on you should only use the new module *model_estimation* instead.\n",
      "  DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "import h5py                                     # Importing the h5 package.\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import PyDynamic  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1.2.79'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from PyDynamic import __version__ as version\n",
    "version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    \n",
      "    Input matrices have dimensions: (2000, 6291), where 2000 represents number of measurements in time\n",
      "    and 6291 represents number of cycles.\n"
     ]
    }
   ],
   "source": [
    "filename = 'Sensor_data_2kHz.h5'                # Data filename.\n",
    "f = h5py.File(filename, 'r')                    # Importing the h5 file. \n",
    "\n",
    "#print(\"Keys: %s\" % f.keys())\n",
    "a_group_key = list(f.keys())[0]\n",
    "\n",
    "data = list(f[a_group_key])                     # Transforming data into list\n",
    "\n",
    "sensorADC=[]                                       # Initialising a list \"sensor\" and\n",
    "for i in range(11):                                # Filling it with data from all sensors \n",
    "    sensorADC.append(pd.DataFrame(data[i][:][:]))\n",
    "\n",
    "for i in range(11):                             \n",
    "    sensorADC[i]=sensorADC[i].iloc[:,:6291]           # Cuting the last cycle because it contains all zero elements.\n",
    "\n",
    "print(\"\"\"    \n",
    "    Input matrices have dimensions: %s, where %s represents number of measurements in time\n",
    "    and %s represents number of cycles.\"\"\" % (np.shape(sensorADC[0]),np.shape(sensorADC[0])[0],np.shape(sensorADC[0])[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Converting into SI units "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "offset=[0, 0, 0, 0, 0.00488591, 0.00488591, 0.00488591,  0.00488591, 1.36e-2, 1.5e-2, 1.09e-2]\n",
    "gain=[5.36e-9, 5.36e-9, 5.36e-9, 5.36e-9, 3.29e-4, 3.29e-4, 3.29e-4, 3.29e-4, 8.76e-5, 8.68e-5, 8.65e-5]\n",
    "b=[1, 1, 1, 1, 1, 1, 1, 1, 5.299641744, 5.299641744, 5.299641744]\n",
    "k=[250, 1, 10, 10, 1.25, 1, 30, 0.5, 2, 2, 2]\n",
    "units=['[Pa]', '[g]', '[g]', '[g]', '[kN]', '[bar]', '[mm/s]', '[A]', '[A]', '[A]', '[A]']\n",
    "\n",
    "sensor=[0]*len(sensorADC)\n",
    "\n",
    "for i in range(len(sensorADC)):\n",
    "    sensor[i]=((sensorADC[i]*gain[i])+offset[i])*b[i]*k[i]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "sensor=pd.read_csv(r'C:\\Users\\jugo01\\Desktop\\sensor_units.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### If you have problems with previous step, you can skip conversion into SI units by running next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sensor=sensorADC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading of train and test data\n",
    "*Note: see 2_Machine_Learning_using_Best_Fourier_Coefficients.ipynb. Data were split into train and test data for k=85% \n",
    "Based on this, target_train_vector and target_test_vector were provided. These vectors will be used for the Fourier transform into frequency domain.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import h5py\n",
    "\n",
    "train_test1= h5py.File(\"Train_test_data_split\",\"r\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_train_vector=train_test1[\"target_train_vector\"]\n",
    "target_test_vector=train_test1[\"target_test_vector\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Converting arrays into data frames:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_train_vector=pd.DataFrame(target_train_vector)\n",
    "target_test_vector=pd.DataFrame(target_test_vector)\n",
    "target=list(target_train_vector[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "So, after this step main data to work on are lists: \n",
    "\n",
    "\"sensor_train\" with their class labels \"train_target\"\n",
    " \n",
    "and \n",
    " \n",
    "\"sensor_test\" with their class labels \"test_target\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Traning data for one sensor has dimensions:  (2000, 5347) ,      ('sensor_train') \n",
      "and it's target vector has length:  (5347, 1) ,               ('target_train_vector') \n",
      "\n",
      "Testing data for one sensor has dimensions:  (2000, 944) ,      ('sensor_test') \n",
      "and it's target vector has length:  (944, 1) ,        ('target_test_vector') \n",
      "\n"
     ]
    }
   ],
   "source": [
    "sensor_train=[0]*11\n",
    "sensor_test=[0]*11\n",
    "\n",
    "for i in range(11):\n",
    "    sensor_train[i]=sensor[i].loc[:,target_train_vector.index]\n",
    "\n",
    "print(\"Traning data for one sensor has dimensions: \", sensor_train[10].shape,\",      ('sensor_train') \")\n",
    "print(\"and it's target vector has length: \", target_train_vector.shape,\",               ('target_train_vector') \\n\")\n",
    "\n",
    "for i in range(11):\n",
    "    sensor_test[i]=sensor[i].loc[:,target_test_vector.index]\n",
    "\n",
    "print(\"Testing data for one sensor has dimensions: \", sensor_test[10].shape,\",      ('sensor_test') \")\n",
    "print(\"and it's target vector has length: \", target_test_vector.shape,\",        ('target_test_vector') \\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can have a look at the data from one sensor after splitting for better understanding of structure for next steps. Number of rows is 2000 and each column is one random measurement cycle. Table shows only first five samples in time (five rows) for each cycle. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>...</th>\n",
       "      <th>5337</th>\n",
       "      <th>5338</th>\n",
       "      <th>5339</th>\n",
       "      <th>5340</th>\n",
       "      <th>5341</th>\n",
       "      <th>5342</th>\n",
       "      <th>5343</th>\n",
       "      <th>5344</th>\n",
       "      <th>5345</th>\n",
       "      <th>5346</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-0.052346</td>\n",
       "      <td>0.035353</td>\n",
       "      <td>0.150057</td>\n",
       "      <td>0.165182</td>\n",
       "      <td>0.007607</td>\n",
       "      <td>0.085452</td>\n",
       "      <td>-0.004497</td>\n",
       "      <td>-0.015737</td>\n",
       "      <td>0.004307</td>\n",
       "      <td>0.030425</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.175565</td>\n",
       "      <td>0.329240</td>\n",
       "      <td>0.210907</td>\n",
       "      <td>0.184818</td>\n",
       "      <td>0.050059</td>\n",
       "      <td>0.108024</td>\n",
       "      <td>0.142261</td>\n",
       "      <td>0.358908</td>\n",
       "      <td>0.161659</td>\n",
       "      <td>0.362495</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.013661</td>\n",
       "      <td>0.062294</td>\n",
       "      <td>0.073791</td>\n",
       "      <td>0.002720</td>\n",
       "      <td>0.045830</td>\n",
       "      <td>0.097229</td>\n",
       "      <td>0.054530</td>\n",
       "      <td>0.065786</td>\n",
       "      <td>0.283926</td>\n",
       "      <td>0.108324</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.218415</td>\n",
       "      <td>0.447456</td>\n",
       "      <td>0.349365</td>\n",
       "      <td>0.299143</td>\n",
       "      <td>0.048282</td>\n",
       "      <td>0.207124</td>\n",
       "      <td>0.431144</td>\n",
       "      <td>0.638964</td>\n",
       "      <td>0.378436</td>\n",
       "      <td>0.698642</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.088680</td>\n",
       "      <td>-0.157930</td>\n",
       "      <td>0.016507</td>\n",
       "      <td>0.020354</td>\n",
       "      <td>0.090973</td>\n",
       "      <td>-0.215526</td>\n",
       "      <td>0.061279</td>\n",
       "      <td>0.034248</td>\n",
       "      <td>0.184389</td>\n",
       "      <td>0.141985</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.369747</td>\n",
       "      <td>0.213356</td>\n",
       "      <td>0.078377</td>\n",
       "      <td>0.058016</td>\n",
       "      <td>-0.075154</td>\n",
       "      <td>0.073886</td>\n",
       "      <td>0.393118</td>\n",
       "      <td>0.545816</td>\n",
       "      <td>0.170193</td>\n",
       "      <td>0.216192</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.040834</td>\n",
       "      <td>-0.081864</td>\n",
       "      <td>0.264404</td>\n",
       "      <td>0.090820</td>\n",
       "      <td>0.189913</td>\n",
       "      <td>0.068431</td>\n",
       "      <td>-0.101957</td>\n",
       "      <td>0.105077</td>\n",
       "      <td>0.037419</td>\n",
       "      <td>0.232769</td>\n",
       "      <td>...</td>\n",
       "      <td>0.012711</td>\n",
       "      <td>0.306476</td>\n",
       "      <td>0.201671</td>\n",
       "      <td>0.298681</td>\n",
       "      <td>0.069645</td>\n",
       "      <td>0.011811</td>\n",
       "      <td>0.371483</td>\n",
       "      <td>0.673113</td>\n",
       "      <td>0.286716</td>\n",
       "      <td>0.490413</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.033578</td>\n",
       "      <td>-0.044107</td>\n",
       "      <td>0.056629</td>\n",
       "      <td>0.022157</td>\n",
       "      <td>0.091190</td>\n",
       "      <td>-0.215849</td>\n",
       "      <td>-0.147002</td>\n",
       "      <td>0.120792</td>\n",
       "      <td>0.155922</td>\n",
       "      <td>0.038234</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.026665</td>\n",
       "      <td>0.382832</td>\n",
       "      <td>0.299239</td>\n",
       "      <td>0.259937</td>\n",
       "      <td>0.135226</td>\n",
       "      <td>0.151655</td>\n",
       "      <td>0.500641</td>\n",
       "      <td>0.763004</td>\n",
       "      <td>0.247051</td>\n",
       "      <td>0.542166</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 5347 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       0         1         2         3         4         5         6     \\\n",
       "0 -0.052346  0.035353  0.150057  0.165182  0.007607  0.085452 -0.004497   \n",
       "1  0.013661  0.062294  0.073791  0.002720  0.045830  0.097229  0.054530   \n",
       "2  0.088680 -0.157930  0.016507  0.020354  0.090973 -0.215526  0.061279   \n",
       "3  0.040834 -0.081864  0.264404  0.090820  0.189913  0.068431 -0.101957   \n",
       "4  0.033578 -0.044107  0.056629  0.022157  0.091190 -0.215849 -0.147002   \n",
       "\n",
       "       7         8         9     ...      5337      5338      5339      5340  \\\n",
       "0 -0.015737  0.004307  0.030425  ... -0.175565  0.329240  0.210907  0.184818   \n",
       "1  0.065786  0.283926  0.108324  ... -0.218415  0.447456  0.349365  0.299143   \n",
       "2  0.034248  0.184389  0.141985  ... -0.369747  0.213356  0.078377  0.058016   \n",
       "3  0.105077  0.037419  0.232769  ...  0.012711  0.306476  0.201671  0.298681   \n",
       "4  0.120792  0.155922  0.038234  ... -0.026665  0.382832  0.299239  0.259937   \n",
       "\n",
       "       5341      5342      5343      5344      5345      5346  \n",
       "0  0.050059  0.108024  0.142261  0.358908  0.161659  0.362495  \n",
       "1  0.048282  0.207124  0.431144  0.638964  0.378436  0.698642  \n",
       "2 -0.075154  0.073886  0.393118  0.545816  0.170193  0.216192  \n",
       "3  0.069645  0.011811  0.371483  0.673113  0.286716  0.490413  \n",
       "4  0.135226  0.151655  0.500641  0.763004  0.247051  0.542166  \n",
       "\n",
       "[5 rows x 5347 columns]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sensor_train[0].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "###  Fast Fourier transform\n",
    "\n",
    "###### Steps:  \n",
    "    \n",
    "- transformation into frequency domain (FFT)\n",
    "- choose amplitudes with highest average absolute value (the top 10%)\n",
    "\n",
    "\n",
    "In this method of feature extraction, data is transformed into frequency domain using FFT function for discrete Fourier transform. More detail about FFT in [1_FFT_and_Reconstruction.ipynb](1_FFT_and_Reconstruction.ipynb)\n",
    "\n",
    "This step is an unsupervised extraction method (i.e. is done without knowledge of the cycle‘s group affiliation) and used is to reduce dimension for further steps.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### A function is created, which takes as input: \n",
    "- data from one sensor `sensor`,                                 \n",
    "- number of samples `n_of_samples`,                                    \n",
    "- percentage of data to choose `N`.\n",
    "\n",
    "Function does fast Fourier transform and chooses N% of sprectrum with highest average of absolute values for each sensor independently. Average of absolute values for one frequency is calculated through all cycles.                                   \n",
    "\n",
    "\n",
    "###### Function returns:\n",
    "- `freq_of_sorted_values` matrix sized [1, N% of features (amplitudes)] where elements are frequencies which are choosen and they are labels for second output from this function.\n",
    "- `sorted_values_matrix` sized [number of cycles, N% of features (amplitudes)] where row represents one cycle and columns are sorted by the average of absolute vales for each frequency (column)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chooseAndReturnOrdered(sensor, n_of_samples, N): \n",
    "    x_measurements=range(sensor.shape[0])                 # Number of measurements samples in time period.\n",
    "    x = np.true_divide(x_measurements, n_of_samples)      # Time values, used  as real time axis.\n",
    "    freq = np.fft.rfftfreq(x.size, 0.0005)                # Frequency axis, can be used for ploting in frequency domain.\n",
    "    fft_amplitudes = np.fft.rfft(sensor,n_of_samples,0)   # Ndarray of amplitudes after fourier transform.\n",
    "    fft_matrix = pd.DataFrame(fft_amplitudes)             # Transforming amplitudes into data frame (matrix)-\n",
    "                                                          # -where one column represents amplitudes of one-\n",
    "                                                          # -cycle.\n",
    "    fft_matrix=fft_matrix.transpose()                     # Transposing to matrix where rows are cycles.\n",
    "    n_rows, n_columns = np.shape(fft_matrix)\n",
    "\n",
    "    print(\"\\nNumber of cycles is: %s, and number of features is: %s\" % (n_rows, n_columns))\n",
    "    fft_matrix.columns = freq                    # Column labels are frequencies. \n",
    "    \n",
    "    # Calculating the average of absolute vales for each frequency (column).\n",
    "    absolute_average_values_from_columns=(np.abs(fft_matrix)).mean()\n",
    "    \n",
    "    # Sorting the fft_matrix by the average of absolute vales for each frequency (column).\n",
    "    fft_matrix=fft_matrix.reindex((np.abs(fft_matrix)).mean().sort_values(ascending=False).index, axis=1)\n",
    "    \n",
    "    # Taking first N percent columns from sorted fft_matrix. \n",
    "    sorted_values_matrix=fft_matrix.iloc[:,:round((N/100.0)*len(freq))]\n",
    "    \n",
    "    n_rows, n_columns = np.shape(sorted_values_matrix)\n",
    "    print(\"\\nNumber of cycles is: %s, and number of selected features is: %s\" % (n_rows, n_columns))\n",
    "    print(np.shape(sorted_values_matrix))\n",
    "    \n",
    "    # Informations about the selected frequencies are columns in sorted data frame. \n",
    "    freq_of_sorted_values=(pd.DataFrame(sorted_values_matrix.columns)).transpose()\n",
    "    print(\"\\nFirst 10 selected frequencies are:\\n\\n %s\" % freq_of_sorted_values.values[:,:10])\n",
    "    \n",
    "    sorted_values_matrix.columns=range(round((N/100.0)*len(freq))) # Resetting the column labels.\n",
    "    print(\"---------------------------------------------------------------------------------\\n\")\n",
    "    # Output \"sorted_values_matrix\" is data frame whose rows-\n",
    "    # -are cycles and columns are selected frequencies. For example,- \n",
    "    # -value at position (i,j) is amplitude for frequency j in cycle i.\n",
    "    \n",
    "    return freq_of_sorted_values, sorted_values_matrix;\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "###### Function execution\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Instead of executing the function, results of extracting 10% of highest amplitudes by FFT can be read in the next steps. Values were obtained by using factor of splitting data into train and test from the above.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import h5py\n",
    "\n",
    "amp_fft1= h5py.File(\"Sorted_vaules_from_all_sensors.hdf5\",\"r\")\n",
    "freq_fft1= h5py.File(\"Sorted_freq_from_all_sensors.hdf5\",\"r\")  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "freq_of_sorted_values=[0]*len(sensor_train)\n",
    "sorted_values_from_all_sensors=[0]*len(sensor_train)\n",
    "for i in range(len(sensor)):\n",
    "    freq_of_sorted_values[i]=freq_fft1[\"freq_of_sorted_values\"+str(i)]\n",
    "    sorted_values_from_all_sensors[i]=amp_fft1[\"sorted_values_from_all_sensors\"+str(i)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(sensor)):\n",
    "    freq_of_sorted_values[i]=pd.DataFrame( freq_of_sorted_values[i])\n",
    "    sorted_values_from_all_sensors[i]=pd.DataFrame( sorted_values_from_all_sensors[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "User is asked to define how many of features from frequency domain will be extracted in this step. Then, the function is executed for each sensor and extracted data is stored in 2 lists containing data frames mentioned above. Lists `freq_of_sorted_values` and `sorted_values_from_all_sensors` store function outputs and further selection of features is continued on list `sorted_values_from_all_sensors`. Informations about frequency are going to be used in for feature extraction from testing data, because these frequencies are pattern learned from training data and used for selecting from the testing data or some new data which need to be predicted. \n",
    "For all sensors, selected frequencies with most dominant amplitudes are listed as output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal and recommended percentage of features for this dataset is 10. \n",
      "\n",
      "Enter a percentage of features: 10\n",
      "\n",
      "\n",
      "\n",
      "Sensor number 0\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Number of cycles is: 5347, and number of features is: 1001\n",
      "\n",
      "Number of cycles is: 5347, and number of selected features is: 100\n",
      "(5347, 100)\n",
      "\n",
      "First 10 selected frequencies are:\n",
      "\n",
      " [[480.   0.  85. 640. 100. 120.   1.   2.   3. 121.]]\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Sensor number 1\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Number of cycles is: 5347, and number of features is: 1001\n",
      "\n",
      "Number of cycles is: 5347, and number of selected features is: 100\n",
      "(5347, 100)\n",
      "\n",
      "First 10 selected frequencies are:\n",
      "\n",
      " [[480.  80.  79. 479.   2. 481.  78.   0. 120.  77.]]\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Sensor number 2\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Number of cycles is: 5347, and number of features is: 1001\n",
      "\n",
      "Number of cycles is: 5347, and number of selected features is: 100\n",
      "(5347, 100)\n",
      "\n",
      "First 10 selected frequencies are:\n",
      "\n",
      " [[161. 162. 120. 163. 160. 406. 819. 826. 413. 820.]]\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Sensor number 3\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Number of cycles is: 5347, and number of features is: 1001\n",
      "\n",
      "Number of cycles is: 5347, and number of selected features is: 100\n",
      "(5347, 100)\n",
      "\n",
      "First 10 selected frequencies are:\n",
      "\n",
      " [[640. 819. 823. 824. 825. 818. 821. 820. 822. 826.]]\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Sensor number 4\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Number of cycles is: 5347, and number of features is: 1001\n",
      "\n",
      "Number of cycles is: 5347, and number of selected features is: 100\n",
      "(5347, 100)\n",
      "\n",
      "First 10 selected frequencies are:\n",
      "\n",
      " [[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]]\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Sensor number 5\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Number of cycles is: 5347, and number of features is: 1001\n",
      "\n",
      "Number of cycles is: 5347, and number of selected features is: 100\n",
      "(5347, 100)\n",
      "\n",
      "First 10 selected frequencies are:\n",
      "\n",
      " [[ 0.  1.  6.  5.  7.  2.  4.  8.  3. 13.]]\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Sensor number 6\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Number of cycles is: 5347, and number of features is: 1001\n",
      "\n",
      "Number of cycles is: 5347, and number of selected features is: 100\n",
      "(5347, 100)\n",
      "\n",
      "First 10 selected frequencies are:\n",
      "\n",
      " [[  0. 640. 999. 998. 996. 997. 993. 995. 994. 992.]]\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Sensor number 7\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Number of cycles is: 5347, and number of features is: 1001\n",
      "\n",
      "Number of cycles is: 5347, and number of selected features is: 100\n",
      "(5347, 100)\n",
      "\n",
      "First 10 selected frequencies are:\n",
      "\n",
      " [[  0. 640.   1.   2. 300.   3.   4.   5. 680.   6.]]\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Sensor number 8\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Number of cycles is: 5347, and number of features is: 1001\n",
      "\n",
      "Number of cycles is: 5347, and number of selected features is: 100\n",
      "(5347, 100)\n",
      "\n",
      "First 10 selected frequencies are:\n",
      "\n",
      " [[240. 880. 400. 239.   0. 241. 238. 540. 242.  60.]]\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Sensor number 9\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Number of cycles is: 5347, and number of features is: 1001\n",
      "\n",
      "Number of cycles is: 5347, and number of selected features is: 100\n",
      "(5347, 100)\n",
      "\n",
      "First 10 selected frequencies are:\n",
      "\n",
      " [[240. 880. 400.   0. 239. 241. 238. 540.  60. 242.]]\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Sensor number 10\n",
      "---------------------------------------------------------------------------------\n",
      "\n",
      "Number of cycles is: 5347, and number of features is: 1001\n",
      "\n",
      "Number of cycles is: 5347, and number of selected features is: 100\n",
      "(5347, 100)\n",
      "\n",
      "First 10 selected frequencies are:\n",
      "\n",
      " [[240. 880. 400.   0. 239. 241. 238. 540. 720.  60.]]\n",
      "---------------------------------------------------------------------------------\n",
      "\n"
     ]
    }
   ],
   "source": [
    "n_of_samples=np.shape(sensor_train[0])[0]\n",
    "\n",
    "N = int(input(\"Optimal and recommended percentage of features for this dataset is 10. \\n\\nEnter a percentage of features: \"))\n",
    "print(\"\\n\\n\")\n",
    "# Initialising the list woth 11 elements, which are data frames \"sorted_value_matrix\" from each sensor.\n",
    "freq_of_sorted_values=[0]*len(sensor_train)\n",
    "sorted_values_from_all_sensors=[0]*len(sensor_train)\n",
    "\n",
    "for i in range(len(sensor_train)):                     \n",
    "    print(\"Sensor number %s\" % i)\n",
    "    print(\"---------------------------------------------------------------------------------\")\n",
    "    freq_of_sorted_values[i],sorted_values_from_all_sensors[i]=chooseAndReturnOrdered(sensor_train[i], n_of_samples, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Features are complex numbers, because output from the Fourier transform is resulting with amplitudes and phase shifts. For methods used here, amplitudes are more important than phase shifts, and features that will be used are absolute values of amplitudes  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Example of frequency labels for features extracted from microphone with the best Fourier coefficients method._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example. We will take first sensor, which is microphone, as an example. The most dominant frequancy for microphone is 0 Hz, and then 480 Hz. That can be seen in `freq_of_sorted_values[0]`. First two columns in `sorted_values_from_all_sensors[0]` are  amplitudes through all measurement cycles for frequencies 0 and 480 Hz, respectively. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>...</th>\n",
       "      <th>90</th>\n",
       "      <th>91</th>\n",
       "      <th>92</th>\n",
       "      <th>93</th>\n",
       "      <th>94</th>\n",
       "      <th>95</th>\n",
       "      <th>96</th>\n",
       "      <th>97</th>\n",
       "      <th>98</th>\n",
       "      <th>99</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>480.0</td>\n",
       "      <td>85.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>640.0</td>\n",
       "      <td>100.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>120.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>...</td>\n",
       "      <td>149.0</td>\n",
       "      <td>25.0</td>\n",
       "      <td>838.0</td>\n",
       "      <td>86.0</td>\n",
       "      <td>128.0</td>\n",
       "      <td>839.0</td>\n",
       "      <td>407.0</td>\n",
       "      <td>68.0</td>\n",
       "      <td>759.0</td>\n",
       "      <td>42.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 100 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    0      1     2    3      4      5    6      7    8    9   ...     90  \\\n",
       "0  0.0  480.0  85.0  1.0  640.0  100.0  2.0  120.0  3.0  6.0  ...  149.0   \n",
       "\n",
       "     91     92    93     94     95     96    97     98    99  \n",
       "0  25.0  838.0  86.0  128.0  839.0  407.0  68.0  759.0  42.0  \n",
       "\n",
       "[1 rows x 100 columns]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "freq_of_sorted_values[0].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Discrete Fourier Transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to include uncertainties in features extraction using BFC, the following is going to be conducted:\n",
    "    \n",
    "- transformation into frequency domain (DFT) with uncertainty propagation through PyDynamic software [1_FFT_and_Reconstruction.ipynb](1_FFT_and_Reconstruction.ipynb)\n",
    "- extraction of amplitudes with highest average absolute value (the top 10%), with corresponding phases and uncertainties.\n",
    "\n",
    "The results of DFT are obtained by using function `Time2AmpPhase_multi`. Function is executed with uncertainties presented in sparse matrix in order to avoid memory errors during calculation for high number of cycles.The usage of sparse matrix is possible since most of covariance values are zeros. \n",
    "\n",
    "###### Function returns:\n",
    "- ` A` np.ndarray sized (M,N) where elements are amplitudes of time domain signals of length N in frequency domain for M cycles. \n",
    "- ` P` np.ndarray sized (M,N) where elements are phases of time domain signals of length N in frequency domain for M cycles. \n",
    "-  `UAP`np.ndarray sized (M, 3N) where elements are squared standard uncertainties of amplitudes and phases and their covariances for M cycles.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#this function is now part of PyDynamic\n",
    "from PyDynamic.uncertainty.propagate_DFT import GUM_DFT\n",
    "from PyDynamic.uncertainty.propagate_DFT import DFT2AmpPhase\n",
    "def Time2AmpPhase_multi(x, Ux, selector=None):\n",
    "    #\"\"\"Transformation from time domain to amplitude and phase for a set of M signals of the same type\n",
    "    # Parameters\n",
    "    # x: np.ndarray of shape (M,N)\n",
    "    # M time domain signals of length N\n",
    "    # Ux: np.ndarray of shape (M,)\n",
    "    # squared standard deviations representing noise variances of the signals x\n",
    "    # selector: np.ndarray of shape (L,), optional\n",
    "    # indices of amplitude and phase values that should be returned; default is 0:N-1\n",
    "    # Returns\n",
    "    \n",
    "    # A: np.ndarray of shape (M,N)\n",
    "    # amplitude values\n",
    "    # P: np.ndarray of shape (M,N)\n",
    "    # phase values\n",
    "    # UAP: np.ndarray of shape (M, 3N)\n",
    "    # diagonals of the covariance matrices: [diag(UPP), diag(UAA), diag(UPA)]\n",
    "    \n",
    "    M, nx = x.shape\n",
    "    assert(len(Ux)==M)\n",
    "    N = nx//2+1\n",
    "    if not isinstance(selector, np.ndarray):\n",
    "        selector = np.arange(nx//2+1)\n",
    "    ns = len(selector)\n",
    "\n",
    "    A = np.zeros((M,ns))\n",
    "    P = np.zeros_like(A)\n",
    "    UAP = np.zeros((M, 3*ns))\n",
    "    CxCos = None\n",
    "    CxSin = None\n",
    "    for m in range(M):\n",
    "        F, UF, CX = GUM_DFT(x[m,:], Ux[m], CxCos, CxSin, returnC=True)\n",
    "        CxCos = CX[\"CxCos\"]\n",
    "        CxSin = CX[\"CxSin\"]\n",
    "        A_m, P_m, UAP_m = DFT2AmpPhase(F, UF, keep_sparse=True)\n",
    "        A[m,:] = A_m[selector]\n",
    "        P[m,:] = P_m[selector]\n",
    "        UAP[m,:ns] = UAP_m.data[0][:N][selector]\n",
    "        UAP[m,ns:2*ns] = UAP_m.data[1][UAP_m.offsets[1]:2*N+UAP_m.offsets[1]][selector]\n",
    "        UAP[m, 2*ns:] = UAP_m.data[0][N:][selector]\n",
    "        \n",
    "    return A, P, UAP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Amplitudes are sorted through all cycles and  10% of the highest ones were extracted. Phases and uncertainties are sorted in the way to follow the sorting method for amplitudes. The values are corresponding to the values of 10% chosen amplitudes.\n",
    "###### Function `chooseAndReturnOrdered_with_uncertainty` returns:\n",
    "- `freq_of_sorted_values` matrix sized [1, N% of features (amplitudes)] where elements are frequencies which are choosen and they are labels for second output from this function.\n",
    "- `sorted_values_amp` sized [number of cycles, N% of features (amplitudes)] where row represents one cycle and columns are sorted by the average of absolute values of amplitudes for each frequency (column).\n",
    "- `sorted_values_phases` - sized [number of cycles, N% of features (phases)] where row represents one cycle and columns of phases are sorted by the absolute values of amplitudes  for each frequency (column). \n",
    "- `sorted_values_uncert_aa` - sized [number of cycles, N% of features (uncertainties)] where row represents one cycle and columns of squared standard uncertainties for amplitudes are sorted by the average of absolute vales for each frequency (column). \n",
    "- `sorted_values_uncert_ap` - sized [number of cycles, N% of features (uncertainties)] where row represents one cycle and columns of covariances for amplitudes and phases are sorted by the average of absolute vales for each frequency (column). \n",
    "- `sorted_values_uncert_pp`  - sized [number of cycles, N% of features (uncertainties)] where row represents one cycle and columns of squared standard uncertainties for phases are sorted by the average of absolute vales for each frequency (column). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from PyDynamic.uncertainty.propagate_DFT import GUM_DFT\n",
    "from PyDynamic.uncertainty.propagate_DFT import DFT2AmpPhase\n",
    "\n",
    "def chooseAndReturnOrdered_with_uncertainty(sensor, n_of_samples, N,sigma):\n",
    "    \n",
    "    x_measurements=range(sensor.shape[0])        \n",
    "    n_of_samples=np.shape(sensor)[0]                            # number of sampling points\n",
    "    x = np.true_divide(x_measurements, n_of_samples)            # time steps \n",
    "    freq=PyDynamic.uncertainty.propagate_DFT.GUM_DFTfreq(x.size, 0.0005)\n",
    "    ux=np.ones(sensor.shape[1])*sigma**2\n",
    "    a,b,c=Time2AmpPhase_multi(sensor.values.transpose(),ux)\n",
    "    a=pd.DataFrame(a) #amplitudes (M,N)\n",
    "    b=pd.DataFrame(b) #phases (M,N)\n",
    "    c=pd.DataFrame(c) #uncertainties (M,3*N)\n",
    "    a.columns = freq                    # Column labels are frequencies. \n",
    "    n_rows, n_columns=a.shape\n",
    "    print(\"\\nNumber of cycles is: %s, and number of features is: %s\" % (n_rows, n_columns))\n",
    "    # Calculating the average of absolute vales for each frequency (column).\n",
    "    absolute_average_values_from_columns=np.abs(a.mean())\n",
    "    # Sorting column indices in amplitudes for sorting phases and uncertainties\n",
    "    sorted_columns=np.argsort(absolute_average_values_from_columns)[::-1]\n",
    "    # Uncertaintites have one dimension larger than amplitudes and phases. # Columns indices(len(a):3*len(a) to follow the sorting)\n",
    "    sorted_columns_unc_ap=(sorted_columns+a.shape[1])\n",
    "    sorted_columns_unc_pp=(sorted_columns+a.shape[1]*2)\n",
    "    sorted_columns_unc=np.concatenate((sorted_columns,sorted_columns_unc_ap,sorted_columns_unc_pp))#sorted indices for uncertainties\n",
    "    # Reindexing all matrices based on columns.\n",
    "    a=a.reindex((np.abs(a)).mean().sort_values(ascending=False).index, axis=1)\n",
    "    b=b.reindex(columns=sorted_columns)\n",
    "    c=c.reindex(columns=sorted_columns_unc)\n",
    "    # Taking first N percent columns from sorted amplitudes,phases and ucertainties. \n",
    "    sorted_values_amp=a.iloc[:,:round((N/100.0)*len(freq))]\n",
    "    sorted_values_phases=b.iloc[:,:round((N/100.0)*len(freq))] \n",
    "    sorted_values_uncert_aa=c.iloc[:,:round((N/100.0)*len(freq))]\n",
    "    sorted_values_uncert_ap=c.iloc[:,len(freq):a.shape[1]+round((N/100.0)*len(freq))]\n",
    "    sorted_values_uncert_pp=c.iloc[:,2*len(freq):a.shape[1]*2+round((N/100.0)*len(freq))]                                             \n",
    "    n_rows, n_columns = np.shape(sorted_values_amp)\n",
    "    print(\"\\nNumber of cycles is: %s, and number of selected features is: %s\" % (n_rows, n_columns))\n",
    "    print(np.shape(sorted_values_amp))\n",
    "    \n",
    "    # Informations about the selected frequencies are columns in sorted data frame. \n",
    "    freq_of_sorted_values=(pd.DataFrame(sorted_values_amp.columns)).transpose()\n",
    "    print(\"\\nFirst 10 selected frequencies are:\\n\\n %s\" % freq_of_sorted_values.values[:,:10])\n",
    "    \n",
    "    # Resetting the column labels.\n",
    "    sorted_values_amp.columns=range(round((N/100.0)*len(freq)))\n",
    "    sorted_values_phases.columns=range(round((N/100.0)*len(freq)))\n",
    "    sorted_values_uncert_aa.columns=range(round((N/100.0)*len(freq)))\n",
    "    sorted_values_uncert_ap.columns=range(round((N/100.0)*len(freq)))\n",
    "    sorted_values_uncert_pp.columns=range(round((N/100.0)*len(freq)))\n",
    "\n",
    "    print(\"---------------------------------------------------------------------------------\\n\")\n",
    "    # Output \"sorted_values_matrix\" is data frame whose rows-\n",
    "    # -are cycles and columns are selected frequencies. For example,- \n",
    "    # -value at position (i,j) is amplitude for frequency j in cycle i.\n",
    "    return freq_of_sorted_values,sorted_values_amp,sorted_values_phases,sorted_values_uncert_aa,sorted_values_uncert_ap, sorted_values_uncert_pp\n",
    "    \n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Function execution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Instead of executing the function, results of extracting 10% of highest amplitudes by DFT can be read in the next steps. Values were obtained by using factor of splitting data into train and test from the above. Sigma value, representing white noise was assumed as 0.1*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import h5py\n",
    "\n",
    "amp_dft2= h5py.File(\"DFTSorted_vaules__from_all_sensors.hdf5\",\"r\")\n",
    "freq_dft2= h5py.File(\"DFTSorted_freq_from_all_sensors.hdf5\",\"r\") \n",
    "ph_dft2= h5py.File(\"DFTSorted_ph_from_all_sensors.hdf5\",\"r\")\n",
    "u_a_dft2= h5py.File(\"DFTSorted_uncer_from_all_sensors_a.hdf5\",\"r\")\n",
    "u_ap_dft2= h5py.File(\"DFTSorted_uncer_from_all_sensors_ap.hdf5\",\"r\")    \n",
    "u_pp_dft2= h5py.File(\"DFTSorted_uncer_from_all_sensors_pp.hdf5\",\"r\")    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "freq_of_sorted_values=[0]*len(sensor_train)\n",
    "sorted_values_amp_from_all_sensors=[0]*len(sensor_train)\n",
    "sorted_phases_from_all_sensors=[0]*len(sensor_train)\n",
    "sorted_uncer_from_all_sensors_a=[0]*len(sensor_train)\n",
    "sorted_uncer_from_all_sensors_ap=[0]*len(sensor_train)\n",
    "sorted_uncer_from_all_sensors_pp=[0]*len(sensor_train)\n",
    "for i in range(len(sensor)):\n",
    "    freq_of_sorted_values[i]=freq_dft2[\"freq_of_sorted_values\"+str(i)]\n",
    "    sorted_values_amp_from_all_sensors[i]=amp_dft2[\"sorted_values_amp_from_all_sensors\"+str(i)]\n",
    "    sorted_phases_from_all_sensors[i]=ph_dft2[\"sorted_phases_from_all_sensors\"+str(i)]\n",
    "    sorted_uncer_from_all_sensors_a[i]=u_a_dft2[\"sorted_uncer_from_all_sensors_a\"+str(i)]\n",
    "    sorted_uncer_from_all_sensors_ap[i]=u_ap_dft2[\"sorted_uncer_from_all_sensors_ap\"+str(i)]\n",
    "    sorted_uncer_from_all_sensors_pp[i]=u_pp_dft2[\"sorted_uncer_from_all_sensors_pp\"+str(i)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(sensor)):\n",
    "    freq_of_sorted_values[i]=pd.DataFrame(freq_of_sorted_values[i])\n",
    "    sorted_values_amp_from_all_sensors[i]=pd.DataFrame(sorted_values_amp_from_all_sensors[i])\n",
    "    sorted_phases_from_all_sensors[i]=pd.DataFrame(sorted_phases_from_all_sensors[i])\n",
    "    sorted_uncer_from_all_sensors_a[i]=pd.DataFrame(sorted_uncer_from_all_sensors_a[i])\n",
    "    sorted_uncer_from_all_sensors_ap[i]=pd.DataFrame(sorted_uncer_from_all_sensors_ap[i])\n",
    "    sorted_uncer_from_all_sensors_pp[i]=pd.DataFrame(sorted_uncer_from_all_sensors_pp[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "User is asked to define how many of features from frequency domain will be extracted in this step. Then, the function is executed for each sensor and extracted data is stored in 6 lists containing data frames mentioned above. Lists:\n",
    "- freq_of_sorted_values \n",
    "- sorted_values_amp_from_all_sensors \n",
    "- sorted_phases_from_all_sensors\n",
    "- sorted_uncer_from_all_sensors_a\n",
    "- sorted_uncer_from_all_sensors_ap\n",
    "- sorted_uncer_from_all_sensors_pp\n",
    "\n",
    "store function outputs and further selection of features is continued  through loop. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#function execution\n",
    "n_of_samples=np.shape(sensor_train[0])[0]\n",
    "\n",
    "N = int(input(\"Optimal and recommended percentage of features for this dataset is 10. \\n\\nEnter a percentage of features: \"))\n",
    "print(\"\\n\\n\")\n",
    "sigma=float(input(\"Assume standard deviation\"))\n",
    "# Initialising the list woth 11 elements, which are data frames \"sorted_value_matrix\" from each sensor.\n",
    "                    \n",
    "freq_of_sorted_values=[0]*len(sensor_train)\n",
    "sorted_values__amp_from_all_sensors=[0]*len(sensor_train)\n",
    "sorted_phases_from_all_sensors=[0]*len(sensor_train)\n",
    "sorted_uncer_from_all_sensors_a=[0]*len(sensor_train)\n",
    "sorted_uncer_from_all_sensors_ap=[0]*len(sensor_train)\n",
    "sorted_uncer_from_all_sensors_pp=[0]*len(sensor_train)\n",
    "\n",
    "    \n",
    "for i in range(len(sensor_train)):                     \n",
    "    print(\"Sensor number %s\" % i)\n",
    "    print(\"---------------------------------------------------------------------------------\")\n",
    "    freq_of_sorted_values[i],sorted_values_amp_from_all_sensors[i],sorted_phases_from_all_sensors[i],sorted_uncer_from_all_sensors_a[i],sorted_uncer_from_all_sensors_ap[i],sorted_uncer_from_all_sensors_pp[i]=chooseAndReturnOrdered_with_uncertainty(sensor_train[i], n_of_samples, N,sigma)\n",
    "    \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freq_of_sorted_values[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note: When amplitudes are small relative to the uncertainty associated with real and imaginary parts , the GUM uncertainty propagation becomes unreliable and a Monte Carlo method is recommended instead. Consequently, GUM2DFT does raise a warning to the user and recommends using a Monte Carlo method instead whenever an element of  is below a pre-defined threshold. The default threshold in GUM2DFT is 1.0, but may be adjusted for specific applications.[7]*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An overview of the results:\n",
    "\n",
    "*Note: Be aware of randomness of splitting data into train and test. Results shown here are for one random split. Functions FFT and DFT were executed for the same split of data.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>...</th>\n",
       "      <th>90</th>\n",
       "      <th>91</th>\n",
       "      <th>92</th>\n",
       "      <th>93</th>\n",
       "      <th>94</th>\n",
       "      <th>95</th>\n",
       "      <th>96</th>\n",
       "      <th>97</th>\n",
       "      <th>98</th>\n",
       "      <th>99</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>15.601349</td>\n",
       "      <td>42.066543</td>\n",
       "      <td>40.512832</td>\n",
       "      <td>19.052201</td>\n",
       "      <td>33.424006</td>\n",
       "      <td>25.445938</td>\n",
       "      <td>9.262591</td>\n",
       "      <td>30.340080</td>\n",
       "      <td>11.713543</td>\n",
       "      <td>18.716029</td>\n",
       "      <td>...</td>\n",
       "      <td>2.037884</td>\n",
       "      <td>6.216052</td>\n",
       "      <td>4.764562</td>\n",
       "      <td>2.592275</td>\n",
       "      <td>12.083291</td>\n",
       "      <td>2.941888</td>\n",
       "      <td>4.113036</td>\n",
       "      <td>7.766285</td>\n",
       "      <td>10.219041</td>\n",
       "      <td>2.441705</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>13.194672</td>\n",
       "      <td>53.251987</td>\n",
       "      <td>24.740147</td>\n",
       "      <td>8.909405</td>\n",
       "      <td>32.139716</td>\n",
       "      <td>51.545917</td>\n",
       "      <td>12.176153</td>\n",
       "      <td>31.096967</td>\n",
       "      <td>17.930797</td>\n",
       "      <td>8.678278</td>\n",
       "      <td>...</td>\n",
       "      <td>2.225403</td>\n",
       "      <td>12.098849</td>\n",
       "      <td>7.352120</td>\n",
       "      <td>7.171204</td>\n",
       "      <td>7.668650</td>\n",
       "      <td>10.721172</td>\n",
       "      <td>8.712233</td>\n",
       "      <td>3.294677</td>\n",
       "      <td>30.347054</td>\n",
       "      <td>14.614420</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 100 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          0          1          2          3          4          5   \\\n",
       "0  15.601349  42.066543  40.512832  19.052201  33.424006  25.445938   \n",
       "1  13.194672  53.251987  24.740147   8.909405  32.139716  51.545917   \n",
       "\n",
       "          6          7          8          9   ...        90         91  \\\n",
       "0   9.262591  30.340080  11.713543  18.716029  ...  2.037884   6.216052   \n",
       "1  12.176153  31.096967  17.930797   8.678278  ...  2.225403  12.098849   \n",
       "\n",
       "         92        93         94         95        96        97         98  \\\n",
       "0  4.764562  2.592275  12.083291   2.941888  4.113036  7.766285  10.219041   \n",
       "1  7.352120  7.171204   7.668650  10.721172  8.712233  3.294677  30.347054   \n",
       "\n",
       "          99  \n",
       "0   2.441705  \n",
       "1  14.614420  \n",
       "\n",
       "[2 rows x 100 columns]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted_values_amp_from_all_sensors[0].head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>...</th>\n",
       "      <th>90</th>\n",
       "      <th>91</th>\n",
       "      <th>92</th>\n",
       "      <th>93</th>\n",
       "      <th>94</th>\n",
       "      <th>95</th>\n",
       "      <th>96</th>\n",
       "      <th>97</th>\n",
       "      <th>98</th>\n",
       "      <th>99</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>20.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>...</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>20.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>...</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 100 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     0     1     2     3     4     5     6     7     8     9   ...    90  \\\n",
       "0  20.0  10.0  10.0  10.0  10.0  10.0  10.0  10.0  10.0  10.0  ...  10.0   \n",
       "1  20.0  10.0  10.0  10.0  10.0  10.0  10.0  10.0  10.0  10.0  ...  10.0   \n",
       "\n",
       "     91    92    93    94    95    96    97    98    99  \n",
       "0  10.0  10.0  10.0  10.0  10.0  10.0  10.0  10.0  10.0  \n",
       "1  10.0  10.0  10.0  10.0  10.0  10.0  10.0  10.0  10.0  \n",
       "\n",
       "[2 rows x 100 columns]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted_uncer_from_all_sensors_a[0].head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>...</th>\n",
       "      <th>90</th>\n",
       "      <th>91</th>\n",
       "      <th>92</th>\n",
       "      <th>93</th>\n",
       "      <th>94</th>\n",
       "      <th>95</th>\n",
       "      <th>96</th>\n",
       "      <th>97</th>\n",
       "      <th>98</th>\n",
       "      <th>99</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-1.387779e-17</td>\n",
       "      <td>-5.551115e-17</td>\n",
       "      <td>-2.775558e-17</td>\n",
       "      <td>1.831868e-15</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-1.110223e-16</td>\n",
       "      <td>-2.775558e-17</td>\n",
       "      <td>5.551115e-17</td>\n",
       "      <td>5.551115e-17</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.110223e-16</td>\n",
       "      <td>-2.220446e-16</td>\n",
       "      <td>-4.218847e-15</td>\n",
       "      <td>2.220446e-16</td>\n",
       "      <td>-4.440892e-16</td>\n",
       "      <td>-3.552714e-15</td>\n",
       "      <td>-2.220446e-15</td>\n",
       "      <td>3.330669e-16</td>\n",
       "      <td>2.220446e-16</td>\n",
       "      <td>-3.264056e-14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-4.336809e-19</td>\n",
       "      <td>1.387779e-17</td>\n",
       "      <td>5.551115e-17</td>\n",
       "      <td>1.998401e-15</td>\n",
       "      <td>6.938894e-18</td>\n",
       "      <td>-5.551115e-17</td>\n",
       "      <td>6.938894e-18</td>\n",
       "      <td>5.551115e-17</td>\n",
       "      <td>1.110223e-16</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-1.665335e-16</td>\n",
       "      <td>-1.831868e-15</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>-6.661338e-16</td>\n",
       "      <td>-9.436896e-16</td>\n",
       "      <td>-8.326673e-16</td>\n",
       "      <td>-4.440892e-16</td>\n",
       "      <td>-1.595946e-16</td>\n",
       "      <td>-6.938894e-15</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2 rows × 100 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    0             1             2             3             4             5   \\\n",
       "0  0.0 -1.387779e-17 -5.551115e-17 -2.775558e-17  1.831868e-15  0.000000e+00   \n",
       "1  0.0 -4.336809e-19  1.387779e-17  5.551115e-17  1.998401e-15  6.938894e-18   \n",
       "\n",
       "             6             7             8             9   ...            90  \\\n",
       "0 -1.110223e-16 -2.775558e-17  5.551115e-17  5.551115e-17  ... -1.110223e-16   \n",
       "1 -5.551115e-17  6.938894e-18  5.551115e-17  1.110223e-16  ...  0.000000e+00   \n",
       "\n",
       "             91            92            93            94            95  \\\n",
       "0 -2.220446e-16 -4.218847e-15  2.220446e-16 -4.440892e-16 -3.552714e-15   \n",
       "1 -1.665335e-16 -1.831868e-15  0.000000e+00 -6.661338e-16 -9.436896e-16   \n",
       "\n",
       "             96            97            98            99  \n",
       "0 -2.220446e-15  3.330669e-16  2.220446e-16 -3.264056e-14  \n",
       "1 -8.326673e-16 -4.440892e-16 -1.595946e-16 -6.938894e-15  \n",
       "\n",
       "[2 rows x 100 columns]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sorted_uncer_from_all_sensors_ap[0].head(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freq_of_sorted_values[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Additional: Transformation to time domain for all sensors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transformation from amplitude and phase to time domain is demonstrated with functions `Reconstruct_time_domain` (on the basis of PyDynamic´s function AmpPhase2Time) and `Reconstruct_time_domain_idft`(on the basis of PyDynamic´s functions AmpPhase2DFT, GUM_iDFT).\n",
    "\n",
    "First, zero arrays of amplitudes, phases and uncertainties are created (A, P, UAP). Then, function copies values of N% of arguments (amplitudes, phases, u_a, u_ap, u_pp) into column indices (frequencies) of initial zero arrays. \n",
    "\n",
    "Function `Reconstruct_time_domain`creates sparse matrix from the values contained in matrices of uncertainties for amplitudes, phases and their covariances, because UAP was previously created from the sparse matrix. For each cycle, function returns:  \n",
    "- `x` (np.ndarray ) – vector of time domain values and \n",
    "- `ux` -  (np.ndarray) – standard squared uncertainties of x.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from scipy.sparse import dia_matrix\n",
    "from PyDynamic.uncertainty.propagate_DFT import AmpPhase2Time\n",
    "\n",
    "def Reconstruct_time_domain(N,frequencies,amplitudes,phases,u_a,u_ap,u_pp):\n",
    "    \n",
    "    M,num=amplitudes.shape\n",
    "    length=int(amplitudes.shape[1])*N+1 \n",
    "    length_2=2*length\n",
    "    length_3=3*length \n",
    "    x=np.zeros((M, length_2-2))\n",
    "    #storing uncertainties in a list of arrays\n",
    "    ux=np.zeros((M, length_2-2))\n",
    "    #predefining A,P,UAP as zero arrays\n",
    "    A =np.zeros((M, length))\n",
    "    P= np.zeros((M, length))\n",
    "    UAP=np.zeros((M, 3*length))\n",
    "    assert(amplitudes.shape==phases.shape)\n",
    "    # Indices of columns with highest amplitudes in original matrix (resulted from DFT) are accessible from the\n",
    "    #sorted frequencies of all sensors\n",
    "    Index_amplitudes=frequencies[:,:N] \n",
    "    # Defining offsets for sparse matrix  \n",
    "    offset_UAP=[0,length,-length]\n",
    "     #indices(columns) of 10% highest amplitude values\n",
    "    Index_amplitudes=Index_amplitudes.astype(int)\n",
    "    col = np.array(Index_amplitudes[0])\n",
    "    #Values of 10% highest amplitudes(first N columns of input amplitudes) are copied in A,P,UAP in corresponding indices \n",
    "    #of columns. Other columns are zeros. \n",
    "    amp_col=np.arange(N)\n",
    "    A[:, [col]]= amplitudes[:, [amp_col]]\n",
    "    P[:, [col]]= phases[:, [amp_col]]\n",
    "    UAP[:,[col]]=u_a[:,[amp_col]]\n",
    "    UAP[:,[col+length]]=u_ap[:,[amp_col]]\n",
    "    UAP[:,[col+length_2]]=u_pp[:,[amp_col]]\n",
    "    for m in range(M): \n",
    "        # Defining diagonals for sparse matrix  \n",
    "        diag1=np.zeros(length_2)\n",
    "        diag2=np.zeros(length_2)\n",
    "        diag1[:length]=UAP[m,:length]\n",
    "        diag1[length:]=UAP[m,length_2:]\n",
    "        diag3=np.zeros(length_2)\n",
    "        diag2[offset_UAP[1]:length_2+offset_UAP[1]]=UAP[m][length:length_2] \n",
    "        diag3[offset_UAP[1]:length_2+offset_UAP[1]]=UAP[m][length:length_2]\n",
    "        diagonals =[diag1,diag2,diag3]\n",
    "        # Creating sparse matrix with three diagonals. Diag1 is the main diagonal.\n",
    "        Sparse_matr=dia_matrix((diagonals,offset_UAP),shape=((length_2, length_2)))\n",
    "        X,UX=AmpPhase2Time(A[m,:], P[m,:], Sparse_matr)\n",
    "        x[m,:] = X\n",
    "        ux[m,:]=np.diag(UX)\n",
    "  \n",
    "    return    x,ux\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function `Reconstruct_time_domain_idft` gradually performs transformation from amplitudes and phases to real and imaginary parts and then to time domain, taking into account squared standard uncertainties of amplitudes and phases.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "from PyDynamic.uncertainty.propagate_DFT import AmpPhase2DFT,GUM_iDFT\n",
    "import h5py\n",
    "def Reconstruct_time_domain_idft(N,frequencies,amplitudes,phases,u_a,u_pp):\n",
    "   \n",
    "    M,num=amplitudes.shape\n",
    "    length=int(amplitudes.shape[1])*N+1 #promijeniti index,generisati\n",
    "    length_2=2*length\n",
    "    x=np.zeros((M, length_2-2))\n",
    "    #storing uncertainties in a list of arrays\n",
    "    ux=np.zeros((M, length_2-2))\n",
    "    #predefining A,P,UAP as zero arrays\n",
    "    A =np.zeros((M, length))\n",
    "    P= np.zeros_like(A)\n",
    "    UAP=np.zeros((M, length_2)) #UAP contains squared standard uncertainties of amplitudes and phases\n",
    "    assert(amplitudes.shape==phases.shape)\n",
    "    # Indices of columns with highest amplitudes\n",
    "    Index_amplitudes=frequencies[:,:N] #promijeniti index,generisati\n",
    "    Index_amplitudes=Index_amplitudes.astype(int)\n",
    "    col = np.array(Index_amplitudes[0])\n",
    "    #indices(columns) of 10% highest amplitude values\n",
    "    #first N columns of input amplitudes\n",
    "    amp_col=np.arange(N)\n",
    "    A[:, [col]]= amplitudes[:, [amp_col]]\n",
    "    P[:, [col]]= phases[:, [amp_col]]\n",
    "    UAP[:,[col]]=u_a[:,[amp_col]]\n",
    "    UAP[:,[col+length]]=u_pp[:,[amp_col]]\n",
    "    for m in range(M):\n",
    "        F,UF=AmpPhase2DFT(A[m,:], P[m,:], UAP[m,:])\n",
    "        X,UX=GUM_iDFT(F, UF)\n",
    "        x[m,:] = X\n",
    "        ux[m,:]=np.diag(UX)\n",
    "     \n",
    "    return x,ux\n",
    "   \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Instead of executing the function, reconstructed time domain signals obtained from  extracting 10% of highest amplitudes by DFT can be read from hdf5 files in the next steps. Sigma value, representing white noise was assumed as 0.1*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reading data instead of execution \n",
    "import os\n",
    "import h5py\n",
    "\n",
    "x_time1= h5py.File(\"Reconstructed-time-signals10.hdf5\",\"r\")\n",
    "ux_time1= h5py.File(\"Reconstructed-uncert-time-signals10.hdf5\",\"r\")    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_time=[0]*len(sensor_train)\n",
    "ux_time=[0]*len(sensor_train)\n",
    "for i in range(len(sensor)):\n",
    "    x_time[i]=x_time1[\"x_time\"+str(i)]\n",
    "    ux_time[i]=ux_time1[\"ux_time\"+str(i)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Function execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#function execution\n",
    "x_time=[0]*len(sensor_train)\n",
    "ux_time=[0]*len(sensor_train)\n",
    "N=10 #percentage of amplitudes that were extracted from DFT results\n",
    "for i in range(len(sensor_train)):                     \n",
    "    print(\"Sensor number %s\" % i)\n",
    "    x_time[i],ux_time[i]=Reconstruct_time_domain(N,freq_of_sorted_values[i].values, sorted_values_amp_from_all_sensors[i].values,sorted_phases_from_all_sensors[i].values,sorted_uncer_from_all_sensors_a[i].values,sorted_uncer_from_all_sensors_ap[i].values,sorted_uncer_from_all_sensors_pp[i].values)\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#function execution\n",
    "x_time=[0]*len(sensor_train)\n",
    "ux_time=[0]*len(sensor_train)\n",
    "for i in range(len(sensor_train)):                     \n",
    "    print(\"Sensor number %s\" % i)\n",
    "    x_time[i],ux_time[i]=Reconstruct_time_domain_idft(N,freq_of_sorted_values[i].values, sorted_values_amp_from_all_sensors[i].values,sorted_phases_from_all_sensors[i].values,sorted_uncer_from_all_sensors_a[i].values,sorted_uncer_from_all_sensors_pp[i].values)\n",
    "       \n",
    "           "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualization of time domain signal through all cycles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f28f887d6e3d48e99726f13d3205248b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='sensor', options=(0, 1, 2, 3, 4, 5, 6, 7, 8, 9), value=0), Dropdow…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_sensor(sensor, cycle)>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook\n",
    "\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact, interact_manual\n",
    "units=['[Pa]', '[g]', '[g]', '[g]', '[kN]', '[bar]', '[mm/s]', '[A]', '[A]', '[A]', '[A]']\n",
    "labels1 = ['Microphone','Vibration plain bearing','Vibration piston rod','Vibration ball bearing', 'Axial force','Pressure','Velocity','Active current','Motor current phase 1','Motor current phase 2','Motor current phase 3']\n",
    "def plot_sensor(sensor,cycle):\n",
    "   \n",
    "    plt.figure(figsize=(15,12))\n",
    "    plt.plot(np.arange(0,1,0.0005),sensor_train[sensor].values.transpose()[cycle,:],label=\"Input time values\")\n",
    "    plt.ylabel(str(units[sensor]))\n",
    "    plt.xlabel(\"Time [s]\")\n",
    "    plt.title(str(labels1[sensor]))\n",
    "    plt.errorbar(np.arange(0,1,0.0005),x_time[sensor][cycle],yerr=np.sqrt((ux_time[sensor][cycle])),label=\"Reconstructed time values with DFT\", ecolor='orangered',\n",
    "            color='green')\n",
    "    # Adding legend to the plot    \n",
    "        \n",
    "\n",
    "interact(plot_sensor, sensor=range(10),cycle=range(sorted_values_amp_from_all_sensors[0].shape[0]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3d209b2c26514594869edb5f90c35092",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Dropdown(description='sensor', options=(0, 1, 2, 3, 4, 5, 6, 7, 8, 9), value=0), IntSlid…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.plot_sensor1(sensor, cycle)>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact, interact_manual\n",
    "units=['[Pa]', '[g]', '[g]', '[g]', '[kN]', '[bar]', '[mm/s]', '[A]', '[A]', '[A]', '[A]']\n",
    "labels1 = ['Microphone','Vibration plain bearing','Vibration piston rod','Vibration ball bearing', 'Axial force','Pressure','Velocity','Active current','Motor current phase 1','Motor current phase 2','Motor current phase 3']\n",
    "def plot_sensor1(sensor,cycle):\n",
    "    plt.figure(figsize=(15,12))\n",
    "    plt.plot(np.arange(0,1,0.0005),sensor_train[sensor].values.transpose()[cycle,:], label=\"Input time values\")\n",
    "    plt.ylabel(str(units[sensor]))\n",
    "    plt.xlabel(\"Time [s]\")\n",
    "    plt.title(str(labels1[sensor]))\n",
    "    plt.errorbar(np.arange(0,1,0.0005),x_time[sensor][cycle],yerr=np.sqrt((ux_time[sensor][cycle])),label=\"Reconstructed time values with DFT\", ecolor='orangered',\n",
    "            color='green')\n",
    "    # Adding legend to the plot    \n",
    "    plt.legend(loc='best', frameon=True)\n",
    "interact(plot_sensor1,sensor=range(10),cycle=widgets.IntSlider(min=0, max=sorted_values_amp_from_all_sensors[0].shape[0], step=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_test1.close()\n",
    "amp_fft1.close()\n",
    "freq_fft1.close()\n",
    "amp_dft2.close()\n",
    "freq_dft2.close()\n",
    "ph_dft2.close()\n",
    "u_a_dft2.close()\n",
    "u_ap_dft2.close()\n",
    "u_pp_dft2.close() \n",
    "x_time1.close()\n",
    "ux_time1.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References:\n",
    "\n",
    "[1]  PTB, ZeMA, - Deep dive into the ZeMA machine learning (ppt), January 2019\n",
    "\n",
    "[2]  https://www.nti-audio.com/en/support/know-how/fast-fourier-transform-fft\n",
    "\n",
    "[3]  http://www.sthda.com/english/wiki/correlation-test-between-two-variables-in-r\n",
    "\n",
    "[4]  https://en.wikipedia.org/wiki/Pearson_correlation_coefficient\n",
    "\n",
    "[5]  Edouard Duchesnay, Tommy Löfstedt, - Statistics and Machine Learning in Python, March 2018\n",
    "\n",
    "[6]  https://ipywidgets.readthedocs.io/en/latest/examples/Using%20Interact.html\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}