demo_signal_feature_extraction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Digital EMG Signal Feature Extraction"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Import packages."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "import numpy as np\n",
    "from filter import BandpassFilter1D, NotchFilter1D\n",
    "from processing import MeanShift1D, Detrend1D, Resample1D, Normalize1D\n",
    "from segment import SegmentND\n",
    "from feature_extraction import *\n",
    "%matplotlib inline"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "EMG signal processing functions"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [],
   "source": [
    "# process signal of each channel\n",
    "def process_signal1d(x, raw_fs=1000, low_fs=1, high_fs=120, notch_fs=60, Q=20, window_size=250, step_size=50, target_fs=512):\n",
    "    \"\"\"\n",
    "    @param x: signal of a single channel\n",
    "    @param raw_fs: original sampling rate\n",
    "    @param low_fs: low cutoff frequency\n",
    "    @param high_fs: high cutoff frequency\n",
    "    @param notch_fs: notch cutoff frequency\n",
    "    @param Q: Q factor\n",
    "    @param window_size: windows size for detrending\n",
    "    @param step_size: step size for detrending\n",
    "    @param target_fs: target sampling rate for resampling step\n",
    "    \"\"\"\n",
    "    # mean-correct signal\n",
    "    x_processed = MeanShift1D.apply(x)\n",
    "    # filtering noise\n",
    "    x_processed = BandpassFilter1D.apply(x_processed, low_fs, high_fs, order=4, fs=raw_fs)\n",
    "    x_processed = NotchFilter1D.apply(x_processed, notch_fs, Q=Q, fs=raw_fs)\n",
    "    # detrend\n",
    "    x_processed = Detrend1D.apply(x_processed, detrend_type='locreg', window_size=window_size, step_size=step_size)\n",
    "    # resample\n",
    "    x_processed = Resample1D.apply(x_processed, raw_fs, target_fs)\n",
    "    # rectify\n",
    "    x_processed = abs(x_processed)\n",
    "    # normalize\n",
    "    x_processed = Normalize1D.apply(x_processed, norm_type='min_max')\n",
    "    return x_processed"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [],
   "source": [
    "# process multi-channel signal\n",
    "def process_signalnd(x, raw_fs=1000, low_fs=1, high_fs=120, notch_fs=60, Q=20, window_size=250, step_size=50, target_fs=512):\n",
    "    \"\"\"\n",
    "    @param x: signal of a single channel\n",
    "    @param raw_fs: original sampling rate\n",
    "    @param low_fs: low cutoff frequency\n",
    "    @param high_fs: high cutoff frequency\n",
    "    @param notch_fs: notch cutoff frequency\n",
    "    @param Q: Q factor\n",
    "    @param window_size: windows size for detrending\n",
    "    @param step_size: step size for detrending\n",
    "    @param target_fs: target sampling rate for resampling step\n",
    "    \"\"\"\n",
    "    num_channels = x.shape[1]\n",
    "    x_processed = np.array([])\n",
    "    for i in range(num_channels):\n",
    "        # process each channel\n",
    "        channel_processed = process_signal1d(x[:, i], raw_fs, low_fs, high_fs, notch_fs, Q, window_size, step_size, target_fs)\n",
    "        channel_processed = np.expand_dims(channel_processed, axis=1)\n",
    "        if i == 0:\n",
    "            x_processed = channel_processed\n",
    "            continue\n",
    "        x_processed = np.hstack((x_processed, channel_processed))\n",
    "    return x_processed"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Load multi-channel signal from file."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EMG samples:\n",
      "      CH1    CH2    CH3    CH4\n",
      "0  727.0  725.0  724.0  726.0\n",
      "1  736.0  733.0  737.0  735.0\n",
      "2  737.0  735.0  743.0  739.0\n",
      "3  738.0  738.0  749.0  741.0\n",
      "4  739.0  739.0  749.0  736.0\n",
      "5  737.0  735.0  740.0  730.0\n",
      "6  734.0  730.0  735.0  727.0\n",
      "7  732.0  729.0  730.0  725.0\n",
      "8  731.0  725.0  726.0  723.0\n",
      "9  724.0  720.0  720.0  714.0\n",
      "\n",
      "Number of samples:\n",
      " 10000\n",
      "\n",
      "Number of channels:\n",
      " 4\n"
     ]
    }
   ],
   "source": [
    "# load emg signal from csv file\n",
    "emg_raw = pd.read_csv('./data/example/emg2.csv')\n",
    "\n",
    "# print first ten EMG samples\n",
    "print('EMG samples:\\n', emg_raw.head(10))\n",
    "\n",
    "# print number of samples\n",
    "print('\\nNumber of samples:\\n', len(emg_raw))\n",
    "\n",
    "# print number of channels\n",
    "print('\\nNumber of channels:\\n', len(emg_raw.columns))\n",
    "\n",
    "# get emg raw data\n",
    "emg_raw = emg_raw.values"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Clean raw signal."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [],
   "source": [
    "# Clean emg raw signal\n",
    "raw_fs = 1000\n",
    "target_fs = 512\n",
    "low_fs = 10\n",
    "high_fs = 120\n",
    "notch_fs = 60\n",
    "Q = 20\n",
    "windows = 512\n",
    "steps = 50\n",
    "emg_processed = process_signalnd(emg_raw, raw_fs, low_fs, high_fs, notch_fs, Q, windows, steps, target_fs)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "Segment multi-channel signal."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Number of segments:\n",
      " 144\n"
     ]
    }
   ],
   "source": [
    "window_size = 512\n",
    "step_size = 32\n",
    "segments = SegmentND.apply(emg_processed, window_size, step_size)\n",
    "\n",
    "# print number of segments\n",
    "print('\\nNumber of segments:\\n', len(segments))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Part 1: EMG Feature Extraction\n",
    "\n",
    "With the aim to facilitate the classification process, the raw signal should be transformed into relevant data structure or feature vector to highlight the important data for pattern recognition. There are three types of features in EMG control system: 1) time domain, 2) frequency domain, 3) time-frequency domain."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, channel length): (512, 4)\n"
     ]
    }
   ],
   "source": [
    "# extract feature from one segment of signal\n",
    "segment = segments[0]\n",
    "\n",
    "# print shape of segment\n",
    "# (number of channels, number of samples per channel)\n",
    "print('Shape (# of channels, channel length): {}'.format(segment.shape))\n",
    "\n",
    "window_size = 10    # window length to get feature"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### 1.1. Time Domain Features\n",
    "\n",
    "The time domain features are the most popular in EMG pattern recognition, because they are easy and quick to calculate as well as they do not require transformation. The time domain features usually are computed based on signal amplitude, and resultant values."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Maximum Peak Value (MPK)\n",
    "\n",
    "It is the maximum absolute vlaue of the signal window.\n",
    "\n",
    "$$\\text{MPK}_{k}(x)=\\text{max}{|x_{i}|}$$\n",
    "\n",
    "where $x_{i}$ is value of signal window $k$"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, feature length): (52, 4)\n"
     ]
    }
   ],
   "source": [
    "peaks = MaxPeak.apply(segment, window_size)\n",
    "# print feature shape\n",
    "print('Shape (# of channels, feature length): {}'.format(peaks.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Mean Value (MN)\n",
    "\n",
    "When the analysis time is not too short. the mean is essentially a constant. Therefore, any shift in value of the mean are indictive of changes in potential.\n",
    "\n",
    "$$\\text{MN}_{k}(x)=\\frac{1}{N}\\sum_{i=1}^{N}{x_{i}}$$\n",
    "\n",
    "where $x_{i}$ is value of signal window $k$, $N$ is the length of the signal window."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, feature length): (52, 4)\n"
     ]
    }
   ],
   "source": [
    "means = Mean.apply(segment, window_size)\n",
    "# print feature shape\n",
    "print('Shape (# of channels, feature length): {}'.format(means.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Variance (VAR)\n",
    "\n",
    "It is a statistical feature.\n",
    "\n",
    "$$\\text{VAR}_{k}(x)=\\sqrt{\\frac{1}{N}\\sum_{i=1}^{N}{(x_{i}-\\bar{x})^{2}}}$$\n",
    "\n",
    "where $x_{i}$ is value of signal window $k$, $\\bar{x}$ is mean value of signal window, $N$ is the length of the signal window."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, feature length): (52, 4)\n"
     ]
    }
   ],
   "source": [
    "vars = Variance.apply(segment, window_size)\n",
    "# print feature shape\n",
    "print('Shape (# of channels, feature length): {}'.format(vars.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Standard Deviation (STD)\n",
    "\n",
    "Is is a statistical feature.\n",
    "\n",
    "$$\\text{STD}_{k}(x)=\\sqrt{\\frac{\\sum_{i=1}{N}{(x_{i} - \\bar{x})}^{2}}{N-1}}$$\n",
    "\n",
    "where $x_{i}$ is value of signal window $k$, $\\bar{x}$ is mean value of signal window, $N$ is the length of the signal window."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, feature length): (52, 4)\n"
     ]
    }
   ],
   "source": [
    "stds = StandardDeviation.apply(segment, window_size)\n",
    "# print feature shape\n",
    "print('Shape (# of channels, feature length): {}'.format(stds.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Skewness (S)\n",
    "\n",
    "It is statistical feature. It measures the degree of deviation from the symmetry of Normal or Gaussian distribution.\n",
    "\n",
    "$$\\text{S}_{k_{mc}}(x)=\\frac{\\sum_{i=1}^{N}{\\frac{(x_i-\\bar{x})^{3}}{N}}}{(\\sum_{i=1}^{N}{\\frac{(x_{i} - \\bar{x})}{N})^{3/2}}$$\n",
    "\n",
    "where $x_{i}$ is value of signal window $k$, $\\bar{x}$ is mean value of signal window, $N$ is the length of signal window, $S_{k_{mc}}$ is the moment coefficient of skewness."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, feature length): (52, 4)\n"
     ]
    }
   ],
   "source": [
    "skews = Skewness.apply(segment, window_size)\n",
    "# print feature shape\n",
    "print('Shape (# of channels, feature length): {}'.format(skews.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Kurtosis (K)\n",
    "\n",
    "It is a statistical feature. It measures the peakedness or flatness of a distribution.\n",
    "\n",
    "$$\\text{K}_{k_{mc}}(x)=\\frac{\\sum_{i=1}{N}{\\frac{(x_{i} - \\bar{x})^{4}}{N}}}{\\sum_{i=1}^{N}(\\frac{(x_{i} - \\bar{x})^{2}}{N})^{2}} - 3$$\n",
    "\n",
    "where $x_{i}$ is value of signal window $k$, $\\bar{x}$ is mean value of signal window, $N$ is the length of signal window, $K_{k_{mc}}$ is the moment coefficient of kurtosis."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, feature length): (52, 4)\n"
     ]
    }
   ],
   "source": [
    "kurts = Kurtosis.apply(segment, window_size)\n",
    "# print feature shape\n",
    "print('Shape (# of channels, feature length): {}'.format(kurts.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Root Mean Square (RMS)\n",
    "\n",
    "It is modeled as amplitude modulated Gaussian random process whose RMS is related to the constant force and non-fatiguing contraction.\n",
    "\n",
    "$$\\text{RMS}_{k}=\\sqrt{\\frac{1}{N}\\sum_{i=1}^{N}{x_{i}^{2}}}$$\n",
    "\n",
    "where $x_{i}$ is value of signal segment $k$, $N$ is the length of the signal segment."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, feature length): (52, 4)\n"
     ]
    }
   ],
   "source": [
    "rmss = RootMeanSquare.apply(segment, window_size)\n",
    "# print feature shape\n",
    "print('Shape (# of channels, feature length): {}'.format(rmss.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Waveform Length (WL)\n",
    "\n",
    "It is the cumulative length of the waveform over the segment. The resultant values indicate a measure of waveform amplitude, frequency and duration all within a single parameter.\n",
    "\n",
    "$$\\text{WL}_{k}=\\sum_{i=1}^{N-1}{|x_{i+1}-x_{i}|}$$\n",
    "\n",
    "where $x_{i}$ and $x_{i+1}$ is consecutive values of signal segment $k$, $N$ is the length of the signal segment."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, feature length): (52, 4)\n"
     ]
    }
   ],
   "source": [
    "wls = WaveformLength.apply(segment, window_size)\n",
    "# print feature shape\n",
    "print('Shape (# of channels, feature length): {}'.format(wls.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### William Amplitude (WAMP)\n",
    "\n",
    "It calculates the number of times that the absolute value of the difference between EMG signal amplitude of two consecutive samples ($x_{i}$ and $x_{i+1}$) exceeds a predetermined threshold $\\epsilon$.\n",
    "\n",
    "$$\\text{WAMP}_{k}=\\sum_{i=1}^{N-1}{f(x_{i}-x_{i+1})}$$\n",
    "\n",
    "$$f(x)=\\begin{cases}{1}&&{x>\\epsilon}\\\\{0}&&{\\text{otherwise}\\end{cases}$$\n",
    "\n",
    "where $x_{i}$ and $x_{i+1}$ is consecutive values of signal segment $k$, $N$ is the length of the signal segment, $\\epsilon$ is threshold value."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, feature length): (52, 4)\n"
     ]
    }
   ],
   "source": [
    "eps = 0.3\n",
    "wamps = WillisonAmplitude.apply(segment, window_size, eps)\n",
    "# print feature shape\n",
    "print('Shape (# of channels, feature length): {}'.format(wamps.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### 1.2. Frequency Domain Features\n",
    "\n",
    "The frequency domain features are based on estimated Power Spectrum Density (PSD) of signal and are computed by Fast Fourier Transforms (FFT). In comparison with time domain features, frequency domain features requires more computation."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Power Spectrum Density (PSD)\n",
    "\n",
    "The feature shows at what point frequency variations are strong and weak.\n",
    "\n",
    "$$\\text{PSD}={\\Bigg|\\sum_{i=0}^{N-1}{x_{i}e^{\\frac{-j2{\\pi}ki}{N}}}\\Bigg|}^{2}$$\n",
    "\n",
    "where $k=0,1,...,N-1$, $N$ is the length if the signal data, $x_{i}$ represents the discrete samples of signal data."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Frequency Median (FMD)\n",
    "\n",
    "The frequency median splits the power spectrum density into two equal parts.\n",
    "\n",
    "$$\\text{FMD}=\\frac{1}{2}\\sum_{i=1}^{N}{\\text{PSD}_{i}}$$\n",
    "\n",
    "where $N$ is the length of the power spectrum density, and $\\text{PSD}_{i}$ is the $i^{th}$ line of power spectrum density."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, feature length): (257, 4)\n"
     ]
    }
   ],
   "source": [
    "min_fs = 2\n",
    "fs = 512\n",
    "window_size = 512\n",
    "average = 'median'\n",
    "fmds = PSD.apply(segment, window_size, fs, average)\n",
    "# print feature shape\n",
    "print('Shape (# of channels, feature length): {}'.format(fmds.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Frequency Mean (FMN)\n",
    "\n",
    "$$\\text{FMN}=\\frac{\\sum_{i=1}^{M}{f_{i}\\text{PSD}_{i}}}{\\sum_{i=1}^{M}{\\text{PSD}_{i}}}$$\n",
    "\n",
    "where $N$ is the length of the power spectrum density, $f_{i}=(i\\times\\text{sampling_rate})$ and $\\text{PSD}_{i}$ is the $i^{th}$ line of power spectrum density."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, feature length): (257, 4)\n"
     ]
    }
   ],
   "source": [
    "fs = 512\n",
    "window_size = 512\n",
    "average = 'mean'\n",
    "fmns = PSD.apply(segment, window_size, fs, average)\n",
    "# print feature shape\n",
    "print('Shape (# of channels, feature length): {}'.format(fmns.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### 1.3. Time-Frequency domain features\n",
    "\n",
    "Time-frequency domain representation can localize the energy of the signal both in time and in frequency, allowing a more accurate description of the physical phenomenon, but these features generally requires a transformation that could be computationally heavy."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Short TIme Fourier Transform (STFT).\n",
    "\n",
    "As the extension of Fourier transform method by devising the input signal into segments, by doing each signal in the each window can be assumed to be stationary.\n",
    "\n",
    "$$\\text{STFT}_{x}(t, w)=\\int{W^{*}(\\tau-t)x(\\tau)e^{-jw\\tau}}d\\tau$$\n",
    "\n",
    "where $W$ is the window function, $*$ is complex conjugate, $\\tau$ represents time, $w$ stands for frequency."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, frequency, stft per frequency): (64, 49, 4)\n"
     ]
    }
   ],
   "source": [
    "stfts = STFT.apply(segment,np.arange(1, 128))\n",
    "# return stfts\n",
    "# print feature shape\n",
    "print('Shape (# of channels, frequency, stft per frequency): {}'.format(stfts.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "##### Wavelet Transform (WT)\n",
    "\n",
    "It is a transform where a signal is integrated with a shifted and scaled mother wavelet function.\n",
    "\n",
    "$$\\text{WT}_{x}(a,b)=\\int{{x(t)}{(\\frac{1}{\\sqrt{a}})\\Psi^{*}{(\\frac{t-b}{a})}}dt}$$\n",
    "\n",
    "where $x(t)$ is the input signal, $\\Psi^{*}$ is the complex conjugate of the mother wavelet transform, $\\Psi(\\frac{t-b}{a})$ is the shifted and scaled version of the wavelet at time $b$ and scale $a$."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape (# of channels, # of scale, cwt per scale): (50, 512, 4)\n"
     ]
    }
   ],
   "source": [
    "window_size = 128\n",
    "widths = 50\n",
    "# cwts = CWT.apply(segment, window_size, widths)\n",
    "L = segment.shape[0]\n",
    "num_channels = segment.shape[1]\n",
    "if isinstance(widths, int):\n",
    "    widths = np.arange(1, widths + 1)\n",
    "cwts = np.array([])\n",
    "for i in range(num_channels):\n",
    "    cwt = signal.cwt(segment[:, i], signal.ricker, widths, dtype=np.float64)\n",
    "    if i == 0:\n",
    "        cwts = cwt\n",
    "        continue\n",
    "    cwts = np.dstack((cwts, cwt))\n",
    "# print feature shape\n",
    "print('Shape (# of channels, # of scale, cwt per scale): {}'.format(cwts.shape))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Part 3: Practices\n",
    "\n",
    "### Practice 1:\n",
    "\n",
    "Try to change some parameters from feature extraction functions and find the differences.\n",
    "\n",
    "### Practice 2:\n",
    "\n",
    "Try to extract the features list of signal segments. Save the features into a list."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}