This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
############################################################## | |
# | |
# 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