This script has been written to pre-process and process uncalibrated triaxial accelerometer (Noraxon sensors) with sample frequency at 500 Hz.

Esports-players-kinematics-behaviour_accelerometry

# -*- coding: utf-8 -*-
"""
Created on Tue Aug  8 10:46:51 2023

@author: antoine dupuy

This script has been written to analyse uncalibrated triaxial accelerometer
(Noraxon sensors) with sample frequency at 500 Hz.
It uses a customed function to extract data from each participant csv file (FileOperation)
which will not be available here.
The accelerometers were located respectively on the right hand (between the distal
3rd metacarpal process, the proximal 2st metacarpal process, and the proximal
4th metacarpal process with the x-axis aligned in the continuation of the hand and pointing forward),
forearm (half-distance from the olecranon fossa and the imaginary line between
the ulnar and radial styloid process with the x-axis aligned in the continuation
of the forearm and pointing forward),
and arm (half-distance from the acromion and the lateral epicondyle with the x-axis
aligned in the continuation of the arm and pointing to the elbow).
All the following analysis is concerning the hand, the forearm, and the arm.
The script is reading .csv files. It is working with numpy dataframe.
In order, this script is: removing unwanted columns and rows, defining columns,
calculating the effective time of recording, removing outliers, calculating the percentage of unchanged values
imputing for missing values, subtracting the mean of the signal (because accelerometers are uncalibrated),
changing the unit (mG) into m.s-2, use an FFT to determine the filter cut-off frequency,
applying a 2nd order Butterworth bandpass (0.5-11 Hz) filter, integrating the signal to calculate velocity,
integrate again to calculate displacement, calculating position from displacement, accelerationa and velocity norms,
Root Mean Square (RMS) of acceleration and velocity, acceleration ratio between hand and forearm & forearm and arm,
number of zero crossings in velocity, cumulative distance and displacement ratio, and ellipse of confidence
of hand displacements.
Then exporting the results of all participants of the same folder to a .csv file.

"""

"import libraries"

import pandas as pd
import shutil
from scipy.integrate import cumulative_trapezoid
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import matplotlib.transforms as transforms
import scipy as sp
import numpy as np
import os

"definition of the different functions used for analysing data" #From https://gist.github.com/CarstenSchelp/b992645537660bda692f218b562d0712

def confidence_ellipse(x, y, ax, n_std=3.0, facecolor='none', **kwargs):
    """
    Created by Carsten Schelp
    Plot of the covariance confidence ellipse of *x* and *y*.

    Parameters
    ----------
    x, y : array-like, shape (n, )
        Input data.

    ax : matplotlib.axes.Axes
        The axes object to draw the ellipse into.

    n_std : float
        The number of standard deviations to determine the ellipse's radiuses.

    **kwargs
        Forwarded to `~matplotlib.patches.Ellipse`

    Returns
    -------
    matplotlib.patches.Ellipse
    """
    if x.size != y.size:
        raise ValueError("x and y must be the same size")

    cov = np.cov(x, y)
    pearson = cov[0, 1]/np.sqrt(cov[0, 0] * cov[1, 1])
    # Using a special case to obtain the eigenvalues of this
    # two-dimensional dataset.
    ell_radius_x = np.sqrt(1 + pearson)
    ell_radius_y = np.sqrt(1 - pearson) 
    ellipse = Ellipse((0, 0), width=ell_radius_x * 2, height=ell_radius_y * 2,
                      facecolor=facecolor, **kwargs)

    # Calculating the standard deviation of x from
    # the squareroot of the variance and multiplying
    # with the given number of standard deviations.
    scale_x = np.sqrt(cov[0, 0]) * n_std
    mean_x = np.mean(x)

    # calculating the standard deviation of y ...
    scale_y = np.sqrt(cov[1, 1]) * n_std
    mean_y = np.mean(y)

    transf = transforms.Affine2D() \
        .rotate_deg(45) \
        .scale(scale_x, scale_y) \
        .translate(mean_x, mean_y)

    ellipse.set_transform(transf + ax.transData)
    return ax.add_patch(ellipse)

def count_sign_changes(values, threshold = 0.01): # Calculate the number of zero-crossings (values input = signal to analyse; threshold = the value to cross
                                                    # in the negative and positive side)
    if len(values) < 2:
        return 0

    sign_changes = 0
    current_sign = 1 if values[0] >= threshold else -1 if values[0] <= -threshold else 0
    i=1
    for value in values[1:]:
        new_sign = 1 if value >= threshold else -1 if value <= -threshold else 0

        if new_sign != 0 and new_sign != current_sign:
            #print (i,value)
            sign_changes += 1
            current_sign = new_sign
        i=i+1
    return sign_changes-1

##################################################
"time played"
time_played = (endingrow - beginningrow)/500 # Determine how much time participants played
print('time played = ', time_played, ' s')

"delta t"
fs = 500
dt = 1/fs
# print(dt)

df_output['time'] = time
#print (time)

"Columns identification"
Accel_1_Ax = df.iloc[:, 3]  # df['Ultium EMG.Internal Accel 1 Ax (mG)']
Accel_1_Ay = df.iloc[:, 4]  # df['Ultium EMG.Internal Accel 1 Ay (mG)']
Accel_1_Az = df.iloc[:, 5]  # df['Ultium EMG.Internal Accel 1 Az (mG)']
Accel_2_Ax = df.iloc[:, 6]  # df['Ultium EMG.Internal Accel 2 Ax (mG)']
Accel_2_Ay = df.iloc[:, 7]  # df['Ultium EMG.Internal Accel 2 Ay (mG)']
Accel_2_Az = df.iloc[:, 8]  # df['Ultium EMG.Internal Accel 2 Az (mG)']
Accel_3_Ax = df.iloc[:, 9]  # df['Ultium EMG.Internal Accel 3 Ax (mG)']
Accel_3_Ay = df.iloc[:, 10]  # df['Ultium EMG.Internal Accel 3 Ay (mG)']
Accel_3_Az = df.iloc[:, 11]  # df['Ultium EMG.Internal Accel 3 Az (mG)']
#print(Accel_1_Ay)
#plt.plot(time, Accel_1_Ay)
# plt.title('raw acceleration in the Medio-Lateral plane in function of time')
# plt.xlabel('time in seconds')
# plt.ylabel('acceleration in mG')

"Outliers management: Replacing data by missing values"
# Take the mean of the signal
mean_1_x = np.mean(Accel_1_Ax)
mean_1_y = np.mean(Accel_1_Ay)
mean_1_z = np.mean(Accel_1_Az)
mean_2_x = np.mean(Accel_2_Ax)
mean_2_y = np.mean(Accel_2_Ay)
mean_2_z = np.mean(Accel_2_Az)
mean_3_x = np.mean(Accel_3_Ax)
mean_3_y = np.mean(Accel_3_Ay)
mean_3_z = np.mean(Accel_3_Az)
#print(mean_1_x)

# Define a threshold
thh_1_x_h = mean_1_x + (3*np.std(Accel_1_Ax)) # 3 standard deviations threshold to remove accelerometers' falls and tasks' unrelated movements
thh_1_x_l = mean_1_x - (3*np.std(Accel_1_Ax))
thh_1_y_h = mean_1_y + (3*np.std(Accel_1_Ay))
thh_1_y_l = mean_1_y - (3*np.std(Accel_1_Ay))
thh_1_z_h = mean_1_z + (3*np.std(Accel_1_Az))
thh_1_z_l = mean_1_z - (3*np.std(Accel_1_Az))
thh_2_x_h = mean_2_x + (3*np.std(Accel_2_Ax))
thh_2_x_l = mean_2_x - (3*np.std(Accel_2_Ax))
thh_2_y_h = mean_2_y + (3*np.std(Accel_2_Ay))
thh_2_y_l = mean_2_y - (3*np.std(Accel_2_Ay))
thh_2_z_h = mean_2_z + (3*np.std(Accel_2_Az))
thh_2_z_l = mean_2_z - (3*np.std(Accel_2_Az))
thh_3_x_h = mean_3_x + (3*np.std(Accel_3_Ax))
thh_3_x_l = mean_3_x - (3*np.std(Accel_3_Ax))
thh_3_y_h = mean_3_y + (3*np.std(Accel_3_Ay))
thh_3_y_l = mean_3_y - (3*np.std(Accel_3_Ay))
thh_3_z_h = mean_3_z + (3*np.std(Accel_3_Az))
thh_3_z_l = mean_3_z - (3*np.std(Accel_3_Az))
#print(thh_1_x_h)
#print(thh_1_x_l)

# Replace outliers (above or below 12 sd) by 9999
Accel_1_Ax[(Accel_1_Ax > thh_1_x_h) | (Accel_1_Ax < thh_1_x_l)] = 9999
Accel_1_Ay[(Accel_1_Ay > thh_1_y_h) | (Accel_1_Ay < thh_1_y_l)] = 9999
Accel_1_Az[(Accel_1_Az > thh_1_z_h) | (Accel_1_Az < thh_1_z_l)] = 9999
Accel_2_Ax[(Accel_2_Ax > thh_2_x_h) | (Accel_2_Ax < thh_2_x_l)] = 9999
Accel_2_Ay[(Accel_2_Ay > thh_2_y_h) | (Accel_2_Ay < thh_2_y_l)] = 9999
Accel_2_Az[(Accel_2_Az > thh_2_z_h) | (Accel_2_Az < thh_2_z_l)] = 9999
Accel_3_Ax[(Accel_3_Ax > thh_3_x_h) | (Accel_3_Ax < thh_3_x_l)] = 9999
Accel_3_Ay[(Accel_3_Ay > thh_3_y_h) | (Accel_3_Ay < thh_3_y_l)] = 9999
Accel_3_Az[(Accel_3_Az > thh_3_z_h) | (Accel_3_Az < thh_3_z_l)] = 9999
#plt.plot(time, Accel_1_Ax)
#plt.axis([150, 175, -100, 10000])

# Percentage of data unchanged
total_data_1_Ax = len(Accel_1_Ax)
total_data_1_Ay = len(Accel_1_Ay)
total_data_1_Az = len(Accel_1_Az)
total_data_2_Ax = len(Accel_2_Ax)
total_data_2_Ay = len(Accel_2_Ay)
total_data_2_Az = len(Accel_2_Az)
total_data_3_Ax = len(Accel_3_Ax)
total_data_3_Ay = len(Accel_3_Ay)
total_data_3_Az = len(Accel_3_Az)
#print(total_data_1_Ax)

total_removed_1_Ax = np.count_nonzero(Accel_1_Ax == 9999)
total_removed_1_Ay = np.count_nonzero(Accel_1_Ay == 9999)
total_removed_1_Az = np.count_nonzero(Accel_1_Az == 9999)
total_removed_2_Ax = np.count_nonzero(Accel_2_Ax == 9999)
total_removed_2_Ay = np.count_nonzero(Accel_2_Ay == 9999)
total_removed_2_Az = np.count_nonzero(Accel_2_Az == 9999)
total_removed_3_Ax = np.count_nonzero(Accel_3_Ax == 9999)
total_removed_3_Ay = np.count_nonzero(Accel_3_Ay == 9999)
total_removed_3_Az = np.count_nonzero(Accel_3_Az == 9999)
#print(total_removed_1_Ay)

percentage_data_unchanged_1_Ax = 100-(total_removed_1_Ax*100)/total_data_1_Ax
percentage_data_unchanged_1_Ay = 100-(total_removed_1_Ay*100)/total_data_1_Ay
percentage_data_unchanged_1_Az = 100-(total_removed_1_Az*100)/total_data_1_Az
percentage_data_unchanged_2_Ax = 100-(total_removed_2_Ax*100)/total_data_2_Ax
percentage_data_unchanged_2_Ay = 100-(total_removed_2_Ay*100)/total_data_2_Ay
percentage_data_unchanged_2_Az = 100-(total_removed_2_Az*100)/total_data_2_Az
percentage_data_unchanged_3_Ax = 100-(total_removed_3_Ax*100)/total_data_3_Ax
percentage_data_unchanged_3_Ay = 100-(total_removed_3_Ay*100)/total_data_3_Ay
percentage_data_unchanged_3_Az = 100-(total_removed_3_Az*100)/total_data_3_Az
#print(percentage_data_unchanged_1_Ay)

data_unchanged_hand = (percentage_data_unchanged_1_Ax + percentage_data_unchanged_1_Ay + percentage_data_unchanged_1_Az)/3
data_unchanged_forearm = (percentage_data_unchanged_2_Ax + percentage_data_unchanged_2_Ay + percentage_data_unchanged_2_Az)/3
data_unchanged_arm = (percentage_data_unchanged_3_Ax + percentage_data_unchanged_3_Ay + percentage_data_unchanged_3_Az)/3
print(data_unchanged_hand, '%')
print(data_unchanged_forearm, '%')
print(data_unchanged_arm, '%')


"Imputation for missing values"
# Replace 9999 by Not a Number (nan)
Accel_1_Ax.replace(9999, float('nan'), inplace=True)
Accel_1_Ay.replace(9999, float('nan'), inplace=True)
Accel_1_Az.replace(9999, float('nan'), inplace=True)
Accel_2_Ax.replace(9999, float('nan'), inplace=True)
Accel_2_Ay.replace(9999, float('nan'), inplace=True)
Accel_2_Az.replace(9999, float('nan'), inplace=True)
Accel_3_Ax.replace(9999, float('nan'), inplace=True)
Accel_3_Ay.replace(9999, float('nan'), inplace=True)
Accel_3_Az.replace(9999, float('nan'), inplace=True)
# #print(Accel_1_Ax[57000:57500])

df_output['Accel_1_Ax'] = Accel_1_Ax
df_output['Accel_1_Ay'] = Accel_1_Ay
df_output['Accel_1_Az'] = Accel_1_Az
df_output['Accel_2_Ax'] = Accel_2_Ax
df_output['Accel_2_Ay'] = Accel_2_Ay
df_output['Accel_2_Az'] = Accel_2_Az
df_output['Accel_3_Ax'] = Accel_3_Ax
df_output['Accel_3_Ay'] = Accel_3_Ay
df_output['Accel_3_Az'] = Accel_3_Az

Accel_1_Ax = Accel_1_Ax.interpolate(method='linear').ffill().bfill()
Accel_1_Ay = Accel_1_Ay.interpolate(method='linear').ffill().bfill()
Accel_1_Az = Accel_1_Az.interpolate(method='linear').ffill().bfill()
Accel_2_Ax = Accel_2_Ax.interpolate(method='linear').ffill().bfill()
Accel_2_Ay = Accel_2_Ay.interpolate(method='linear').ffill().bfill()
Accel_2_Az = Accel_2_Az.interpolate(method='linear').ffill().bfill()
Accel_3_Ax = Accel_3_Ax.interpolate(method='linear').ffill().bfill()
Accel_3_Ay = Accel_3_Ay.interpolate(method='linear').ffill().bfill()
Accel_3_Az = Accel_3_Az.interpolate(method='linear').ffill().bfill()

# Replace by the mean of the representative window
# Accel_1_Ax.fillna(mean_1_x, inplace=True)
# Accel_1_Ay.fillna(mean_1_y, inplace=True)
# Accel_1_Az.fillna(mean_1_z, inplace=True)
# Accel_2_Ax.fillna(mean_2_x, inplace=True)
# Accel_2_Ay.fillna(mean_2_y, inplace=True)
# Accel_2_Az.fillna(mean_2_z, inplace=True)
# Accel_3_Ax.fillna(mean_3_x, inplace=True)
# Accel_3_Ay.fillna(mean_3_y, inplace=True)
# Accel_3_Az.fillna(mean_3_z, inplace=True)
#plt.plot(time, Accel_1_Ax)


"mean accelerations all axis + subtracting the mean"
# Remove the mean acceleration of the entire signal in all axes (because accelerometers are uncalibrated)
Accel_1_Ax = Accel_1_Ax - mean_1_x
Accel_1_Ay = Accel_1_Ay - mean_1_y
Accel_1_Az = Accel_1_Az - mean_1_z
Accel_2_Ax = Accel_2_Ax - mean_2_x
Accel_2_Ay = Accel_2_Ay - mean_2_y
Accel_2_Az = Accel_2_Az - mean_2_z
Accel_3_Ax = Accel_3_Ax - mean_3_x
Accel_3_Ay = Accel_3_Ay - mean_3_y
Accel_3_Az = Accel_3_Az - mean_3_z
# print(mean_Accel_1_Ax)

"mg to g" # mg is a "milli graviton" or a 1000th of the gravitational constant
Accel_1_Ax = Accel_1_Ax/1000
Accel_1_Ay = Accel_1_Ay/1000
Accel_1_Az = Accel_1_Az/1000
Accel_2_Ax = Accel_2_Ax/1000
Accel_2_Ay = Accel_2_Ay/1000
Accel_2_Az = Accel_2_Az/1000
Accel_3_Ax = Accel_3_Ax/1000
Accel_3_Ay = Accel_3_Ay/1000
Accel_3_Az = Accel_3_Az/1000
# plt.plot(time, Accel_1_Ay)

"g to m.s-2"
g = 9.80665
Accel_1_Ax = Accel_1_Ax*g
Accel_1_Ay = Accel_1_Ay*g
Accel_1_Az = Accel_1_Az*g
Accel_2_Ax = Accel_2_Ax*g
Accel_2_Ay = Accel_2_Ay*g
Accel_2_Az = Accel_2_Az*g
Accel_3_Ax = Accel_3_Ax*g
Accel_3_Ay = Accel_3_Ay*g
Accel_3_Az = Accel_3_Az*g
#plt.plot(time, Accel_1_Ay)

"FFT to determine cut-off frequency"
N = len(Accel_1_Ax)
FFT_Accel_1_Ax = sp.fftpack.fft(Accel_1_Ax.values)
xf = sp.fftpack.fftfreq(N, dt)[0:N//2]
#plt.plot(xf, 2.0/N * np.abs(Accel_1_Ax[0:N//2]))
#plt.plot(xf, FFT_Accel_1_Ax[0:N//2]) # Uncomment to see the FFT of the signal
# plt.axis([0, 2, 0, 4000]) # Uncomment to show the FFT with origin at 0 and to choose the window you want to display (first two digits = x axis)
# plt.title('Fast Fourier Transform (FFT) of hand acceleration in the Medio-Lateral plane (zoom in)')
# plt.xlabel('frequencies (Hz)')
# plt.ylabel('amplitude of the spectral components')

"butterworth filter - acceleration"
# Apply a 2nd order Butterworth filter with bandpass (0.5-11 Hz) cutoff frequency, 500 Hz frequency sampling
sos = sp.signal.butter(2, [0.5, 11], 'bandpass', fs=500, output='sos')
Accel_1_Ax_f = sp.signal.sosfilt(sos, Accel_1_Ax)
Accel_1_Ay_f = sp.signal.sosfilt(sos, Accel_1_Ay)
Accel_1_Az_f = sp.signal.sosfilt(sos, Accel_1_Az)
Accel_2_Ax_f = sp.signal.sosfilt(sos, Accel_2_Ax)
Accel_2_Ay_f = sp.signal.sosfilt(sos, Accel_2_Ay)
Accel_2_Az_f = sp.signal.sosfilt(sos, Accel_2_Az)
Accel_3_Ax_f = sp.signal.sosfilt(sos, Accel_3_Ax)
Accel_3_Ay_f = sp.signal.sosfilt(sos, Accel_3_Ay)
Accel_3_Az_f = sp.signal.sosfilt(sos, Accel_3_Az)
# print(Accel_1_Ax_f)
# plt.plot(time, Accel_1_Ay_f)

df_output['Accel_1_Ax_f'] = Accel_1_Ax_f
df_output['Accel_1_Ay_f'] = Accel_1_Ay_f
df_output['Accel_1_Az_f'] = Accel_1_Az_f

"integration (velocities)"
#Veloc_1_Vx_p = cumtrapz_example(Accel_1_Ax_f,time)
Veloc_1_Vx_p = cumulative_trapezoid(Accel_1_Ax_f, time, initial=0)
Veloc_1_Vy_p = cumulative_trapezoid(Accel_1_Ay_f, time, initial=0)
Veloc_1_Vz_p = cumulative_trapezoid(Accel_1_Az_f, time, initial=0)
Veloc_2_Vx_p = cumulative_trapezoid(Accel_2_Ax_f, time, initial=0)
Veloc_2_Vy_p = cumulative_trapezoid(Accel_2_Ay_f, time, initial=0)
Veloc_2_Vz_p = cumulative_trapezoid(Accel_2_Az_f, time, initial=0)
Veloc_3_Vx_p = cumulative_trapezoid(Accel_3_Ax_f, time, initial=0)
Veloc_3_Vy_p = cumulative_trapezoid(Accel_3_Ay_f, time, initial=0)
Veloc_3_Vz_p = cumulative_trapezoid(Accel_3_Az_f, time, initial=0)
# print(Veloc_1_Vy_p)

df_output['Veloc_1_Vx_p'] = Veloc_1_Vx_p
df_output['Veloc_1_Vy_p'] = Veloc_1_Vy_p
df_output['Veloc_1_Vz_p'] = Veloc_1_Vz_p
df_output['Veloc_2_Vx_p'] = Veloc_2_Vx_p
df_output['Veloc_2_Vy_p'] = Veloc_2_Vy_p
df_output['Veloc_2_Vz_p'] = Veloc_2_Vz_p
df_output['Veloc_3_Vx_p'] = Veloc_3_Vx_p
df_output['Veloc_3_Vy_p'] = Veloc_3_Vy_p
df_output['Veloc_3_Vz_p'] = Veloc_3_Vz_p
# print(Veloc_1_Vy_p)
# plt.plot(time, Veloc_1_Vy_p)

"double integration (displacement)"
# displacement signal
Trav_1_Px = cumulative_trapezoid(Veloc_1_Vx_p, time, initial=0)
Trav_1_Py = cumulative_trapezoid(Veloc_1_Vy_p, time, initial=0)
Trav_1_Pz = cumulative_trapezoid(Veloc_1_Vz_p, time, initial=0)
Trav_2_Px = cumulative_trapezoid(Veloc_2_Vx_p, time, initial=0)
Trav_2_Py = cumulative_trapezoid(Veloc_2_Vy_p, time, initial=0)
Trav_2_Pz = cumulative_trapezoid(Veloc_2_Vz_p, time, initial=0)
Trav_3_Px = cumulative_trapezoid(Veloc_3_Vx_p, time, initial=0)
Trav_3_Py = cumulative_trapezoid(Veloc_3_Vy_p, time, initial=0)
Trav_3_Pz = cumulative_trapezoid(Veloc_3_Vz_p, time, initial=0)
# plt.plot(time, Trav_1_Py)
# plt.plot(Trav_1_Py, Trav_1_Px)

df_output['Trav_1_Px'] = Trav_1_Px
df_output['Trav_1_Py'] = Trav_1_Py
df_output['Trav_1_Pz'] = Trav_1_Pz
df_output['Trav_2_Px'] = Trav_2_Px
df_output['Trav_2_Py'] = Trav_2_Py
df_output['Trav_2_Pz'] = Trav_2_Pz
df_output['Trav_3_Px'] = Trav_3_Px
df_output['Trav_3_Py'] = Trav_3_Py
df_output['Trav_3_Pz'] = Trav_3_Pz

"Displacement to Position"
Pos_1_Px = (df_output['Trav_1_Px'].astype(float) -
            df_output['Trav_1_Px'].shift(-1).astype(float))*-1
Pos_1_Py = (df_output['Trav_1_Py'].astype(float) -
            df_output['Trav_1_Py'].shift(-1).astype(float))*-1
Pos_1_Pz = (df_output['Trav_1_Pz'].astype(float) -
            df_output['Trav_1_Pz'].shift(-1).astype(float))*-1
Pos_2_Px = (df_output['Trav_2_Px'].astype(float) -
            df_output['Trav_2_Px'].shift(-1).astype(float))*-1
Pos_2_Py = (df_output['Trav_2_Py'].astype(float) -
            df_output['Trav_2_Py'].shift(-1).astype(float))*-1
Pos_2_Pz = (df_output['Trav_2_Pz'].astype(float) -
            df_output['Trav_2_Pz'].shift(-1).astype(float))*-1
Pos_3_Px = (df_output['Trav_3_Px'].astype(float) -
            df_output['Trav_3_Px'].shift(-1).astype(float))*-1
Pos_3_Py = (df_output['Trav_3_Py'].astype(float) -
            df_output['Trav_3_Py'].shift(-1).astype(float))*-1
Pos_3_Pz = (df_output['Trav_3_Pz'].astype(float) -
            df_output['Trav_3_Pz'].shift(-1).astype(float))*-1

df_output['Pos_1_Px'] = Pos_1_Px
df_output['Pos_1_Py'] = Pos_1_Py
df_output['Pos_1_Pz'] = Pos_1_Pz
df_output['Pos_2_Px'] = Pos_2_Px
df_output['Pos_2_Py'] = Pos_2_Py
df_output['Pos_2_Pz'] = Pos_2_Pz
df_output['Pos_3_Px'] = Pos_3_Px
df_output['Pos_3_Py'] = Pos_3_Py
df_output['Pos_3_Pz'] = Pos_3_Pz

"zero line plot"
#plt.axhline(y=0, color='r', linestyle='-') # Display an horizontal line at whatever "y" you want

"acceleration norm"
# Normalisation of every point of the signal
A_hand = np.sqrt(Accel_1_Ax_f**2)+np.sqrt(Accel_1_Ay_f**2)+np.sqrt(Accel_1_Az_f**2)
A_forearm = np.sqrt(Accel_2_Ax_f**2) + \
    np.sqrt(Accel_2_Ay_f**2)+np.sqrt(Accel_2_Az_f**2)
A_arm = np.sqrt(Accel_3_Ax_f**2)+np.sqrt(Accel_3_Ay_f**2) + \
    np.sqrt(Accel_3_Az_f**2)
#print('Acceleration normalization: ', A_hand)
# plt.plot(time, A_hand)

df_output['A_hand'] = A_hand

"velocity norm"
V_hand = np.sqrt(Veloc_1_Vx_p**2)+np.sqrt(Veloc_1_Vy_p**2) + \
    np.sqrt(Veloc_1_Vz_p**2)
V_forearm = np.sqrt(Veloc_2_Vx_p**2) + \
    np.sqrt(Veloc_2_Vy_p**2)+np.sqrt(Veloc_2_Vz_p**2)
V_arm = np.sqrt(Veloc_3_Vx_p**2)+np.sqrt(Veloc_3_Vy_p**2) + \
    np.sqrt(Veloc_3_Vz_p**2)
# print(V_hand)
# plt.plot(time, V_hand)

"RMS accelerations"
rms_Accel_1_Ax_f = np.sqrt(np.mean(Accel_1_Ax_f**2))
rms_Accel_1_Ay_f = np.sqrt(np.mean(Accel_1_Ay_f**2))
rms_Accel_1_Az_f = np.sqrt(np.mean(Accel_1_Az_f**2))
rms_Accel_2_Ax_f = np.sqrt(np.mean(Accel_2_Ax_f**2))
rms_Accel_2_Ay_f = np.sqrt(np.mean(Accel_2_Ay_f**2))
rms_Accel_2_Az_f = np.sqrt(np.mean(Accel_2_Az_f**2))
rms_Accel_3_Ax_f = np.sqrt(np.mean(Accel_3_Ax_f**2))
rms_Accel_3_Ay_f = np.sqrt(np.mean(Accel_3_Ay_f**2))
rms_Accel_3_Az_f = np.sqrt(np.mean(Accel_3_Az_f**2))
# print('rms hand acc x: ', rms_Accel_1_Ax_f)

"RMS accelerations all axes"
# RMS of hand, forearm and arm acceleration
rms_acc_hand = rms_Accel_1_Ax_f + rms_Accel_1_Ay_f + rms_Accel_1_Az_f
rms_acc_forearm = rms_Accel_2_Ax_f + rms_Accel_2_Ay_f + rms_Accel_2_Az_f
rms_acc_arm = rms_Accel_3_Ax_f + rms_Accel_3_Ay_f + rms_Accel_3_Az_f
print('RMS acc hand: ', rms_acc_hand, ' m.s-2')
print('RMS acc forearm: ', rms_acc_forearm, ' m.s-2')
print('RMS acc arm: ', rms_acc_arm, ' m.s-2')

"acceleration ratio"
# Ratio of acceleration between hand and forearm & forearm and arm
A_hand_forearm_ratio = rms_acc_hand/rms_acc_forearm
A_forearm_arm_ratio = rms_acc_forearm/rms_acc_arm
print('Acceleration ratio between hand and forearm: ', A_hand_forearm_ratio)
print('Acceleration ratio between forearm and arm: ', A_forearm_arm_ratio)

"Zero crossings in velocity"
# Number of hand direction changes
zc_Veloc_1_Vx = count_sign_changes(Veloc_1_Vx_p) # cf. definition of count_sign_changes function
print(zc_Veloc_1_Vx, 'zero crossings veloc in x')
zc_Veloc_1_Vy = count_sign_changes(Veloc_1_Vy_p)
print(zc_Veloc_1_Vy, 'zero crossings veloc in y')
zc_Veloc_1_Vz = count_sign_changes(Veloc_1_Vz_p)
print(zc_Veloc_1_Vz, 'zero crossings veloc in z')
zc_Veloc_hand = zc_Veloc_1_Vx + zc_Veloc_1_Vy + zc_Veloc_1_Vz
print(zc_Veloc_hand)

"Cumulative distance"
# Sum of the absolute position in all axes
travelled_dist_hand = np.sum(np.abs(Pos_1_Px)) + \
    np.sum(np.abs(Pos_1_Py)) + np.sum(np.abs(Pos_1_Pz))
print('travelled distance: ', travelled_dist_hand, ' m')
travelled_dist_forearm = np.sum(np.abs(Pos_2_Px)) + \
    np.sum(np.abs(Pos_2_Py)) + np.sum(np.abs(Pos_2_Pz))
print('travelled distance: ', travelled_dist_forearm, ' m')
travelled_dist_arm = np.sum(np.abs(Pos_3_Px)) + \
    np.sum(np.abs(Pos_3_Py)) + np.sum(np.abs(Pos_3_Pz))
print('travelled distance: ', travelled_dist_arm, ' m')

"Ellipse area (95% confidence interval)"
# Cf. definition of confidence_ellipse function
n_std = 2 # Approximately 95% of trajectories (3 sd = approx. 99.7% of trajectories)
cov = np.cov(Trav_1_Px, Trav_1_Py)
pearson = cov[0, 1]/np.sqrt(cov[0, 0] * cov[1, 1])
ell_radius_x = np.sqrt(1 + pearson)
ell_radius_y = np.sqrt(1 - pearson)

scaled_ell_x_radius = np.sqrt(cov[0, 0]) * n_std # Scaling the semi-major axis to be able to calculate the exact width of the ellipse of confidence
scaled_ell_y_radius = np.sqrt(cov[1, 1]) * n_std # Scaling the semi-minor axis to be able to calculate the exact width of the ellipse of confidence
width_short = scaled_ell_x_radius*2
width_long = scaled_ell_y_radius*2

print('short radius: ', scaled_ell_x_radius, ' m ; ', 'long radius: ', scaled_ell_y_radius, ' m')
print('width of the ellipse short : ', width_short, ' m ; ', 'and long: ', width_long, ' m')
ell_area_of_use = np.pi*scaled_ell_x_radius*scaled_ell_y_radius # ellipse area formula: pi*a*b (a = length of the semi-major axis; b = length of the semi-minor axis)
print('ellipse area of use: ', ell_area_of_use, 'm2')


"Ellipse plot"
# ax = plt.gca()     # plt.plot(Trav_1_Py, Trav_1_Px) to visualize the trajectories contained into the ellipse
# confidence_ellipse(Trav_1_Py, Trav_1_Px, ax, n_std=2, edgecolor='red')


# -------------------------------------------------------------------------------------------