Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
##############################################################
#
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.