read_accelerator_data.R

##############################################################
#
#  read_totem_data.R
#
#  Import script
#
#  Parses raw Totem Sensor accelerometer data in R
#  And looks for peaks (steps, cycles, repeats) 
#
#  Robin van Emden - robin@pavlov.tech - 2015
#
##############################################################

##############################################################
#
#  Example accelerometer data:
#
#  Timestamp,X,Y,Z,Temperature
#  4646,12332,12212,-1888,23
#  4662,13188,10128,312,23
#  4679,20108,14820,1264,23
#  ...
#
##############################################################

## configuration

#  working directory
wdir <- "C:/Users/robin/Desktop/read_totem_data"
#  sensor data file
sdf <- "total-clean.csv"
#  set working directory
setwd(wdir)
#  sensor data directory
sddir <- "sensor_data"
#  build path to sensor file
sdfpath <- paste0(sddir, "/", sdf)
#  dataframe cache directory
dfcdir <- "df_cache"
#  build path to cache file
dfpath <- paste0(dfcdir,"/",sdf,".Rda")


## check if conversion has been cached and load data

if (!file.exists(dfpath) || TRUE) {

  #  load file and separate by ,
  data <- read.table(sdfpath, sep = "," , header = T)
  # set fractional seconds
  options(digits.secs=3)
  #  add timestamp
  data$timestamp_iso <- ISOdatetime(2015,10,27,16,23,0) + data$Timestamp/1e3
  #  save to cache
  save(data,file=dfpath)
  ## cache file does exist, load dataframe
  
} else { load(file=dfpath)  }


##############################################################

## select subset

## all
title <- "All"
start_value <-   '2015-10-27 16:23:04.645'
end_value <-     '2015-10-27 17:48:00.000' 

## overwrite with subselection

#title <- "Walking to gym"
#start_value <-   '2015-10-27 16:28:00.000'
#end_value <-     '2015-10-27 16:34:00.000' 

#title <- "Rowing with sensor in pocket"
#start_value <-   '2015-10-27 16:42:00.000'
#end_value <-     '2015-10-27 16:55:00.000' 

# clean the data 

data <-na.omit(data)

# selection of the data

data_subset <- data[data$timestamp_iso >= as.POSIXct(start_value ) & data$timestamp_iso <= as.POSIXct(end_value ),]

##############################################################

## XYZ plot

## Running median smoothing

X <- runmed(data_subset$X,31)
Y <- runmed(data_subset$Y,31)
Z <- runmed(data_subset$Z,31)

numbered_time <- data_subset$timestamp_iso

plot(numbered_time, X, type="l", ylab="X,Y,Z", xlab="t", main=paste0(title," - Raw XYZ"))
par(new=TRUE)
plot(numbered_time, Y, type="l", col="green", axes=FALSE, main=NULL, xlab = "" , ylab="" )
par(new=TRUE)
plot(numbered_time, Z, type="l", col="red", axes=FALSE, main=NULL, xlab = "" , ylab="" )

##############################################################



# find peaks

library(wmtsa)
library(zoo)

PeakCycle <- function(data_subset=as.vector(sunspots), SearchFrac=0.02){ 
  
  # the SearchFrac parameter just controls how much to look to either side 
  # of wavCWTPeaks()'s estimated maxima for a bigger value
  # see dRange
  
  Wave <- wavCWT(data_subset)
  WaveTree <- wavCWTTree(Wave)
  WavePeaks <- wavCWTPeaks(WaveTree, snr.min=5)
  WavePeaks_Times <- attr(WavePeaks, which="peaks")[,"iendtime"]
  
  NewPeakTimes <- c()
  dRange <- round(SearchFrac*length(data_subset))
  for(i in 1:length(WavePeaks_Times)){
    NewRange <- max(c(WavePeaks_Times[i]-dRange, 1)):min(c(WavePeaks_Times[i]+dRange, length(data_subset)))
    NewPeakTimes[i] <- which.max(data_subset[NewRange])+NewRange[1]-1
  }
  
  return(matrix(c(NewPeakTimes, data_subset[NewPeakTimes]), ncol=2, dimnames=list(NULL, c("PeakIndices", "Peaks"))))
}

# Sum acceleration

data_subset$total_acceleration <- sqrt(data_subset$X^2+data_subset$Y^2+data_subset$Z^2)

# Run median

data_subset$total_acceleration <- runmed(data_subset$total_acceleration ,81)

# Find peaks

peaks = unique(PeakCycle(as.vector(data_subset$total_acceleration),0.001))

# Filter peaks

peaks <- subset(peaks, peaks[,2]>17000)

# Count peaks

count_peaks <- length(peaks[,1]);

# Create an X-Axis

X_axis <- 1:length(data_subset$total_acceleration)

# AUC approximated by looking at a lot of trapezium figures 
# each time bound between x_i, x_{i+1}, y{i+1} and y_i 
# Using the rollmean of the zoo package
# Very rough estimate of energy expenditure

id <- order(data_subset$total_acceleration)
AUC <- sum(diff(data_subset$total_acceleration[id])*rollmean(X_axis[id],2))
energy = round(AUC / 1000)

# Plot total force and peaks

plot(X_axis,data_subset$total_acceleration, type="l",col="grey", xlab="Series", ylab="Movement", main=paste0(title," - Accelerometer plot"))
points(peaks, col="blue", pch=20)
mtext(paste0("Peaks =", count_peaks, " - Energy expenditure = ",energy), side=3)