Feature Extraction

Due to the stochastic and transient nature of the EMG signal, the raw signal provides little information. As such, hand-crafted features are extracted from the underlying signal to increase its information density before being passed to the machine learning model. These features are extracted from windows (i.e., a predefined number of samples) of EMG. The selection of features (or feature sets) is essential for robust control systems.

This module aims to simplify the feature extraction stage for developers. We have included a number of features from the literature that we have robustly tested and validated. The following code snippet exemplifies how to extract features from raw EMG:

Example Code

import numpy as np
from libemg.feature_extractor import FeatureExtractor
from libemg.utils import get_windows

fe = FeatureExtractor()

# Load data from a csv file:
data = np.loadtxt("emg_data.csv", delimiter=",")
# Split the raw EMG into windows
windows = get_windows(data, 50, 25)

# Extract a list of features
feature_list = ['MAV', 'SSC', 'ZC', 'WL']
features_1 = fe.extract_features(feature_list, windows)

# Extract a predefined feature group
features_2 = fe.extract_feature_group('HTD', windows)

Note that some features have parameters that must be tuned for a particular problem and hardware implementation.

Feature Performance

Each of the ~50 features that were implemented were tested individually with a linear discriminant analysis (LDA) classifier on the 3DCdataset. Note: some of these features are designed to improve robustness to factors (i.e., power line interference, limb position effect, contraction intensity variability), and as such, don’t achieve high accuracy on their own for this gesture recognition task. Their value shouldn’t be discounted when used in a rounded feature set for real-world use. The code associated with Figure 1 can be found in the example code tab below.

alt text

Figure 1: Individual Accuracy of Each Feature on the 3DCDataset

Example Code

import os
import sys
import numpy as np
import pickle
import matplotlib.pyplot as plt
from libemg.datasets import _3DCDataset
from libemg.emg_predictor import EMGClassifier
from libemg.feature_extractor import FeatureExtractor
from libemg.utils import make_regex
from libemg.data_handler import OfflineDataHandler
from libemg.offline_metrics import OfflineMetrics
from libemg.filtering import Filter

if __name__ == "__main__":
    # get the 3DC Dataset using library handle - this downloads the dataset
    dataset = _3DCDataset(save_dir='example_data',
                          redownload=False)
    # take the downloaded dataset and load it as an offlinedatahandler
    odh = dataset.prepare_data(format=OfflineDataHandler)

    # Perform an analysis where we test all the features available of the library individually for within-subject
    # classification accuracy.

    # setup out model type and output metrics
    model = "LDA"
    om = OfflineMetrics()
    metrics = ['CA']

    # get the subject list
    subject_list = list(range(0,22))

    # initialize our feature extractor
    fe = FeatureExtractor()
    feature_list = fe.get_feature_list()

    # get the variable ready for where we save the results
    results = np.zeros((len(feature_list), len(subject_list)))

    for s in subject_list:
        subject_data = odh.isolate_data("subjects",[s])
        subject_train = subject_data.isolate_data("sets",[0])
        subject_test  = subject_data.isolate_data("sets",[1])

        # apply a standardization on the raw data (x - mean)/std
        filter = Filter(sampling_frequency=1000)
        filter_dic = {
            "name": "standardize",
            "data": subject_train
        }
        filter.install_filters(filter_dic)
        filter.filter(subject_train)
        filter.filter(subject_test)

        # from the standardized data, perform windowing
        train_windows, train_metadata = subject_train.parse_windows(200,100)
        test_windows, test_metadata = subject_test.parse_windows(200,100)

        # for each feature in the feature list
        for i, f in enumerate(feature_list):
            train_features = fe.extract_features([f], train_windows)
            test_features  = fe.extract_features([f], test_windows)

            # get the dataset ready for the classifier
            data_set = {}
            data_set['training_features'] = train_features
            data_set['training_labels'] = train_metadata["classes"]

            # setup the classifier
            classifier = EMGClassifier(model)
            classifier.fit(feature_dictionary=data_set.copy())

            # running the classifier analyzes the test data we already passed it
            preds,probs = classifier.run(test_features, test_metadata['classes'])
            # get the CA: classification accuracy offline metric and add it to the results
            results[i,s] = om.extract_offline_metrics(metrics, test_metadata['classes'], preds)[metrics[0]] * 100
            print(f"S{s} {f}: {results[i,s]}%")
    # the feature accuracy is represented by the mean accuracy across the subjects
    mean_feature_accuracy = results.mean(axis=1)
    std_feature_accuracy  = results.std(axis=1)


    plt.bar(feature_list, mean_feature_accuracy, yerr=std_feature_accuracy)
    plt.grid()
    plt.xlabel("Features")
    plt.ylabel("Accuracy")
    plt.show()

Features

Let \(x_{i}\) represents the signal in segment i and \(N\) denotes the number of samples in the timeseries.

Let \(fj\) be the frequency of the spectrum at frequency bin \(j\), \(P_{j}\) is the EMG power spectrum at frequency bin \(j\), and \(M\) is the length of the frequency bin.

Note: every feature and feature group is associated with an abbreviation (e.g., MAV). This is how we refer to a specific feature in our library and is how you should interface with it.

Mean Absolute Value (MAV) ^[1]

The average of the absolute value of the EMG signal. This is one of the most commonly used features for EMG.

\( \text{MAV} = \frac{1}{N}\sum_{i=1}^{N} |x_{i}| \)

Zero Crossings (ZC) ^[1]

The number of times that the amplitude of the EMG signal crosses a zero amplitude threshold.

\( \text{ZC} = \sum_{i=1}^{N-1}[\text{sgn}(x_{i} \times x_{i+1}) \cap |x_{i} - x_{i+1}| \ge \text{threshold}] \)

\( \text{sgn(}x\text{)} = \left\{\begin{array}{lr} 1, & \text{if } x \ge \text{threshold} \\ 0, & \text{otherwise } \end{array}\right\} \)

Slope Sign Change (SSC) ^[1]

The number of times that the slope of the EMG signal changes (i.e., the number of times that it changes between positive and negative).

\( \text{SSC} = \sum_{i=2}^{N-1}[f[(x_{i} - x_{i-1}) \times (x_{i} - x_{i+1})]] \\ \)

\( f(x) = \left\{\begin{array}{lr} 1, & \text{if } x \ge \text{threshold} \\ 0, & \text{otherwise } \end{array}\right\} \)

Waveform Length (WL) ^[1]

The cumulative length of the EMG waveform over the passed in window. This feature is used to measure the overall complexity of the signal (i.e., higher WL means more complex).

\( \text{WL} = \sum_{i=1}^{N-1}|x_{i+1} - x_{i}| \)

L-score (LS) ^[2]

A feature that uses the Legendre Shift Polynomial.

\( \text{LS} = \sum_{i=1}^{N} \text{LSP} * B \)

Maximum Fractal Length (MFL) ^[3]

A nonlinear information feature for the EMG waveform that is relatively invariant to contraction intensity.

\( \text{MFL} = log_{10} \sum_{i=1}^{N-1} | x_{i+1} - x_{i} | \)

Mean Squared Ratio (MSR) ^[4]

Determines the total amount of activity in a window.

\( \text{MSR} = \frac{1}{N}\sum_{i=1}^{N}(x_{i})^{\frac{1}{2}} \)

Willison Amplitude (WAMP) ^[1]

The number of times that there is a difference between amplitudes of two seperate samples exceeding a pre-defined threshold. This feature is related to the firing of motor unit action potentials and muscle contraction force.

\( \text{WAMP} = \sum_{i=1}^{N-1}[f(|x_{n}-x_{n+1})] \\ \)

\( f(x) = \left\{\begin{array}{lr} 1, & \text{if } x \ge \text{threshold} \\ 0, & \text{otherwise } \end{array}\right\} \)

Root Mean Square (RMS) ^[1]

Models the signal as an amplitude modulated Gaussian random process that relates to constant force and non-fatiguing contraction.

\( \text{RMS} = \sqrt{\frac{1}{N}\sum_{i=1}^{N}x_{i}^{2}} \)

Integral of Absolute Value (IAV) ^[1]

The integral of the absolute value. Also known as IEMG.

\( \text{IAV} = \sum_{i=i}^{N} |x_{i}| \)

Difference Absolute Standard Deviation Value (DASDV) ^[1]

Standard deviation of the wavelength.

\( \text{DASDV}=\sqrt{\frac{1}{N-1}\sum_{i=1}^{N-1} (x_{i+1} - x_{i})^2} \)