Permalink
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
220 lines (129 sloc) 6.7 KB
##############################################################################################################
#
# Project: Music as Medicine, Pilot two
# URL: https://pavlov.tech/2018/06/29/syncing-bobbi-music-as-medicine/
# Purpose: First exploration of Music as Medicine ECG data
# Date: 2018-07-03
# Author: Robin van Emden, http://pavlov.tech
# Affiliation: Jheronimus Academy of Data Science, http://jads.nl
#
##############################################################################################################
# A first fast and frugal exploratory Bobbi MaM import/export script
library(here)
library(tuneR)
library(data.table)
library(signal)
source("qrs_detect.R")
source("peak_detect.R")
source("utils.R")
setwd(here::here("R"))
################################### import files #############################################################
# music was precut to between beeps with audacity - length: 0 to 4139.022 seconds
# here, we load resulting WAV between 17 to 24 minutes - that is, the first musical intermezzo
# load and plot music
music <- readWave("../music/music-as-medicine-cut-to-beeps.wav", from = 17, to = 24, units = "minutes")
plot(music)
# load ECG data
subject_list <- list()
files = list.files(path = "../ecg/", pattern="*.csv")
for (i in 1:length(files))
# for now, we work with two data frames seperately
assign(paste0("subject_",i), fread(paste0("../ecg/",files[i]), sep=";"))
# for future analyses:
subject_list[[gsub("*.csv$", "", files[i])]] <-
fread(paste0("../ecg/",files[i]), sep=";")
# delete empty column
subject_2$V14 = NULL
subject_3$V14 = NULL
# add an index
subject_2$idx <- seq.int(nrow(subject_2))
subject_3$idx <- seq.int(nrow(subject_3))
################################### filter ecg ###############################################################
# invert polarity participant 3
subject_3$heartrate <- -1* subject_3$heartrate
# basic high pass filtering
highpass <- butter(2, 0.001, type="high")
subject_2$heartrate <- filter(highpass, x=subject_2$heartrate)
subject_3$heartrate <- filter(highpass, x=subject_3$heartrate)
# minimal low pass filtering
lowpass <- butter(2, 0.4, type="low")
subject_2$heartrate <- filter(lowpass, x=subject_2$heartrate)
subject_3$heartrate <- filter(lowpass, x=subject_3$heartrate)
subject_2 <- peak_detect_wavelet(subject_2, 10)
subject_3 <- peak_detect_wavelet(subject_3, 10)
################################### cut to beeps #############################################################
# find the start and end beeps
s2start = min(which(subject_2$mx > 4000))
s3start = min(which(subject_3$mx < -4000))
s3end = max(which(subject_3$mx > -1000))
# remove all past first and before last beeps
subject_2 <- subject_2[idx>=s2start]
subject_3 <- subject_3[idx>=s3start & idx<=s3end]
# start ms at zero
subject_2$ms <- subject_2$ms - subject_2$ms[1]
subject_3$ms <- subject_3$ms - subject_3$ms[1]
# add a "seconds" column
subject_2$seconds <- subject_2$ms/1000
subject_3$seconds <- subject_3$ms/1000
################################### scale data to audio beeps ###############################################
# scale ECG data to *exactly* length of audio recording between beeps
#
# nice and clean way:
#
# subject_2$seconds <- lin_map(subject_3$seconds, 0, max_audio)
# subject_3$seconds <- lin_map(subject_3$seconds, 0, max_audio)
#
# but for current analysis, we miss some data
# so here, we get a scale factor from s3 and apply to s3 and s2
# get max subject_3
max_s3 <- max(subject_3$second)
max_audio <- 4139.022
# calculate scale factor
scale_factor <- max_audio/max_s3
# apply scale factor to both subjects
subject_3$seconds <- subject_3$seconds * scale_factor
subject_2$seconds <- subject_2$seconds * scale_factor
################################### save between beeps ECG data #############################################
fwrite(subject_2[,c("seconds","heartrate")], "../export/subject_2_hr_between_beeps.csv")
fwrite(subject_3[,c("seconds","heartrate")], "../export/subject_3_hr_between_beeps.csv")
fwrite(subject_3[,c("seconds","mx")], "../export/subject_3_mx_between_beeps.csv")
################################### plot first musical intermezzo ###########################################
start <- 17*1000*60
end <- 24*1000*60
subject_2 <- subject_2[ms>start & ms<end]
subject_3 <- subject_3[ms>start & ms<end]
subject_2 <- qrs_detect(subject_2,120,30)
subject_3 <- qrs_detect(subject_3,120,30)
subject_2 <- delete_close_peaks(subject_2)
subject_3 <- delete_close_peaks(subject_3)
# set start ms and seconds to zero
subject_2$ms <- subject_2$ms - subject_2$ms[1]
subject_3$ms <- subject_3$ms - subject_3$ms[1]
subject_2$seconds <- subject_2$seconds - subject_2$seconds[1]
subject_3$seconds <- subject_3$seconds - subject_3$seconds[1]
# plot overview of HR and MX to check if all seems right
plot(subject_2$seconds/60, subject_2$mx, type = "l")
plot(subject_3$seconds/60, subject_3$mx, type = "l")
plot(subject_2$seconds/60, subject_2$heartrate, type = "l")
points(subject_2[is_peak == TRUE]$seconds/60,subject_2[is_peak == TRUE]$heartrate, col = "red")
plot(subject_3$seconds/60, subject_3$heartrate, type = "l")
points(subject_3[is_peak == TRUE]$seconds/60,subject_3[is_peak == TRUE]$heartrate, col = "red")
################################### ECG close up plot #######################################################
start <- 2*1000*60
end <- 2.3*1000*60
plot(subject_2[ms>start & ms<end]$seconds, subject_2[ms>start & ms<end]$heartrate, type = "l")
points(subject_2[ms>start & ms<end & is_peak == TRUE]$seconds,
subject_2[ms>start & ms<end & is_peak == TRUE]$heartrate, col = "red")
plot(subject_3[ms>start & ms<end]$seconds, subject_3[ms>start & ms<end]$heartrate, type = "l")
points(subject_3[ms>start & ms<end & is_peak == TRUE]$seconds,
subject_3[ms>start & ms<end & is_peak == TRUE]$heartrate, col = "red")
################################### save first musical intermezzo ###########################################
fwrite(subject_2[,c("seconds","heartrate")], "../export/subject_2_hr_between_17_24.csv")
fwrite(subject_3[,c("seconds","heartrate")], "../export/subject_3_hr_between_17_24.csv")
fwrite(subject_2[is_peak==TRUE][,c("seconds","heartrate")], "../export/subject_2_peaks_between_17_24.csv")
fwrite(subject_3[is_peak==TRUE][,c("seconds","heartrate")], "../export/subject_3_peaks_between_17_24.csv")
################################### PySpike export ##########################################################
s2 <- subject_2[is_peak == TRUE][, "seconds"]
s3 <- subject_3[is_peak == TRUE][, "seconds"]
spikes <- t(as.matrix (Cbind(s2,s3)))
write.table(spikes, file = "../export/spikes.txt", row.names = FALSE, col.names = FALSE, na = "")