How to get conditioned Euler angles, angular velocity, and linear acceleration from MetaBase output?

Comments

• dryheat122

  June 2023 edited June 2023

  I got this figured out with considerable help from ChatGPT. Here's the code for anyone trying to do the same thing. For reasons I'm not sure about, for my test data this produces angles in the range -8564 to 10,306. I'm only using the result in plots so it doesn't really matter to me. If you wanted it in the range -360 to 360 you'd need to scale the angles.

  import csv
  import numpy as np
  import pandas as pd
  from scipy.integrate import cumtrapz
  import matplotlib.pyplot as plt
  from scipy.signal import butter,filtfilt
  
  def lowpass(data, cutoff, fs, order):
      normal_cutoff = cutoff / nyq
      # Get the filter coefficients 
      b, a = butter(order, normal_cutoff, btype='low', analog=False)
      y = filtfilt(b, a, data)
      return y
  
  # Read the files and populate the acceleration and angular_velocity arrays
  avelo = []
  filename = "accel-gyro_MetaWear_2023-06-09T13.59.19.303_EABF46360571_Gyroscope_100.000Hz_1.7.3.csv"
  with open(filename) as csvfile:
      csvreader = csv.reader(csvfile)
      next(csvreader) # skip header
      for row in csvreader:
          avelo.append([row[3],row[4],row[5]])
  avelo = np.array(avelo,dtype=float)
  
  accel = []
  filename = "accel-gyro_MetaWear_2023-06-09T13.59.19.303_EABF46360571_Accelerometer_100.000Hz_1.7.3.csv"
  with open(filename) as csvfile:
      csvreader = csv.reader(csvfile)
      next(csvreader) # skip header
      for row in csvreader:
          accel.append([row[3],row[4],row[5]])
  accel = np.array(accel,dtype=float)
  
  # ensure arrays are the same length and create index for plots
  arows,foo = avelo.shape
  brows,foo = accel.shape
  if arows > brows: # then remove rows from avelo
      avelo = avelo[0:brows]
      rows = brows
  elif brows > arows: # then remove rows from accel
      accel = accel[0:arows]
      rows = arows
  else:
      rows = arows
  
  idx = np.arange(0, rows, dtype=int) # index for plots
  
  # Integrate angular velocity to get Euler angles
  euler_angles = cumtrapz(avelo, axis=0, initial=0.0)
  
  # combine the dataframes into one
  dfs = {
      'accel': pd.DataFrame(accel, columns=list('xyz')),
      'angle': pd.DataFrame(euler_angles, columns=list('xyz')),
      'avelo': pd.DataFrame(avelo, columns=list('xyz'))
  }
  df = pd.concat(dfs, axis=1)
  
  # configure and apply lowpass filter to remove noise
  fs = 100        # sample rate, Hz
  cutoff = 3      # desired cutoff frequency of the filter
  order = 4       # desired order of the filter
  nyq = 0.5 * fs  # Nyquist Frequency
  
  parameters = ["accel","angle","avelo"]
  for p in parameters:
      for var in ('x','y','z'):
          df[p,var] =  lowpass(df[p,var], cutoff, fs, order)
  
  # Plot the three series
  ax1 = plt.subplot(3, 1, 1)
  ax1.plot(idx, df['angle']['x'], label='x')
  ax1.plot(idx, df['angle']['y'], label='y')
  ax1.plot(idx, df['angle']['z'], label='z')
  ax1.set_title('angle')
  ax1.legend(loc="best")
  ax2 = plt.subplot(3, 1, 2)
  ax2.plot(idx, df['avelo']['x'], label='x')
  ax2.plot(idx, df['avelo']['y'], label='y')
  ax2.plot(idx, df['avelo']['z'], label='z')
  ax2.set_title('avelo')
  ax2.legend(loc="best")
  ax3 = plt.subplot(3, 1, 3)
  ax3.plot(idx, df['accel']['x'], label='x')
  ax3.plot(idx, df['accel']['y'], label='y')
  ax3.plot(idx, df['accel']['z'], label='z')
  ax3.legend(loc="best")
  ax3.set_title('accel')
  # plt.subplots_adjust(hspace=10)
  plt.tight_layout()
  plt.show()