Data Ingestion

I requested the data from University of Alberta. Alejandro Ley, MBA, Research Assistant, Human Neurophysiology Laboratory, University of Alberta has provided me a binary EMG file with 40 sweeps. The file was unpacked and transformed to data.frame format using Readbin.

# edit this path!
    fl = "/Users/jheintz/Box/02_Data_Analytics_Projects/000_KIN_PhD/01_Classes/KIN457_MotorLearningControl/00_Paper/02_EEG_Data/02_UniversityofAlberta/SampleRC.bin"
    sampledata <- file(fl, "rb")

# Read the header
    header <- readBin(sampledata, "raw", 381)

# The file is 840381 bytes. Minus the header is 
# 840000 bytes. We're reading 2 byte integers so that's 420000 samples over all the channels.
    nsamples <- 420000
    nchannels <- 6
    channeldata <- readBin(sampledata, "int", n=nsamples, size=2, endian="big") 

# Slicing the channels into separate arrays,
# 70000x6 matrix if that's better for you.
    Trigger     <- channeldata[seq(length=nsamples/nchannels, from=1, by=nchannels)]
    channel2    <- channeldata[seq(length=nsamples/nchannels, from=2, by=nchannels)]
    channel3    <- channeldata[seq(length=nsamples/nchannels, from=3, by=nchannels)]
    channel4    <- channeldata[seq(length=nsamples/nchannels, from=4, by=nchannels)]
    StimCurrent <- channeldata[seq(length=nsamples/nchannels, from=5, by=nchannels)]
    EMG         <- channeldata[seq(length=nsamples/nchannels, from=6, by=nchannels)]

# Introducing Time Scale, sampling rate = 5000Hz
    int = 1/5000
    t = seq(0, int * (70000-1), int)*1000

# merging channels, trigger, EMG, and time into one dataframe
    df = data.frame(t, Trigger, channel2, channel3, channel4, StimCurrent, EMG)

# Taking out the average for EMG, StimuCurrent, Channel 2-4
    df[,3:7] = apply(df[,3:7], 2, function(x){ x - mean(x)})
    head(df[, -c(3,4,5)])

Plotting the entire recording

The recording shows different stimulation intensities. The sweeps are not sorted by stimulation intensity.

# Plotting all sweeps
    ggplot(data = df, aes(x = t)) + 
          geom_line(aes(y = StimCurrent/100), col = "blue", size = 0.2) +
          geom_line(aes(y = Trigger/1100), col = "red", size = 0.2) + 
          geom_line(aes(y = EMG/200, col = factor(sweep)), col = "black", size = 0.2) +
          xlab("ms") +
          ylab("mv") + 
          ggtitle("40 Sweeps EMG") + 
            theme(plot.title = element_text(size = 8, 
                                            hjust = 0.5
            ))

Splitting into seperate sweeps for plotting and analysis

Split recording at trigger signal and limiting each sweep to 60 ms.

# Parameters, 
    s = 0; r = 0; start = 1
    df$sweep = as.integer(1:length(df$t))
# looping trough entire dataset and identify trigger. trigger identified if signal is > 10
for (i in 1:length(df$Trigger)){
    # if trigger signal is greater 10, and trigger signal doesn't belong to the same signal 
    # (trigger 1ms (5 x 0.2 ms)) or last record of data set. 
    if ((df$Trigger[i] > 10 & r == 0) | i == length(df$Trigger)){
    # Trigger signal is 3 to 4 sequential numbers, 
        r = 1 # reset if trigger went to 0
    # annotating sweep segment with sequential numbering 
        s  = s + 1 # sweep number
    # print(paste(start, i, s))
        df$sweep[start:i] = s-1
    # calculate duration in ms for each sweep
        l = length(c(start:i))
    # print(paste(l, length(seq(0, by = 0.2, length.out = l))))
       df$t[start:i] = seq(0, by = 0.2, length.out = l)
    # store last data point position as starting point for following sweep
        start = i
  }
  if(r != 0 & df$Trigger[i] < 10) {
    # reset of trigger marker
    r = 0
  }
}

# Deleting the pre-signals before very first stimulus, and trigger
    df =  df[!df$sweep == 0, ]
    head(df)

Plotting sweeps

# function to plot
f.plot = function(temp) {
    ggplot(data = temp, aes(x = t)) + 
        geom_line(aes(y = StimCurrent/100), col = "blue", size = 0.2) +
        geom_line(aes(y = Trigger/1100), col = "red", size = 0.2) + 
        geom_line(aes(y = EMG/200, col = factor(sweep)), col = "black", size = 0.2) +
        xlim(0,60) + ylim(min(df$EMG/200), max(df$EMG/200)) +
        xlab("ms") + ylab("mv") + 
        ggtitle(paste0("Sweep = ", temp$sweep[1])) + theme(plot.title = element_text(size = 8, 
                hjust = 0.5, margin = margin(t = 10, b = -10) 
        ))
}
 
# generating plots for all sweeps
    p.hflx = list()
    for (s in 1:length(unique(df$sweep))){
      temp = subset(df, sweep == s & t < 60)
      p.hflx[[s]] = f.plot(temp)
    }
    
# store plots in list
    n = length(p.hflx)
    nCol = floor(sqrt(n))
    do.call("grid.arrange", c(p.hflx, ncol = nCol))

Identifying H-reflex(max) and M-wave(max) for Recruitment Curve

# time segments in which maximums will be identified
    t.stim = c(seq(0, 1, by = 0.2))
    t.mwave = c(seq(10, 20, by = 0.2))
    t.hrflx = c(seq(30, 45, by = 0.2))

# looping through sweeps and identifying max values conditioned on time segment
    df.recruit = data.frame()
    for (s in 1:length(unique(df$sweep))){
        temp = subset(df, sweep == s & t < 60)
        
        df.temp = data.frame(
            sweep = unique(temp$sweep), 
            Stim.max = max(temp[temp$t %in% t.stim, ]$StimCurrent), 
            Hrflx.max = max(temp[temp$t %in% t.hrflx, ]$EMG), 
            Mwave.max = max(temp[temp$t %in% t.mwave, ]$EMG)
        )
        df.recruit = rbind(df.recruit, df.temp)
    }

# order after stimulation instensity and show head of dataframe
    df.recruit = df.recruit[order(df.recruit$Stim.max), ]
    head(df.recruit)

# Plot recruitment curve
    ggplot(data = df.recruit, aes(x = Stim.max)) +
        geom_point(aes( y = Hrflx.max), size = 1, color = "salmon") + 
        geom_smooth(aes (y = Hrflx.max), 
                    method = stats::loess, 
                    formula = y ~ x,
                    se = FALSE, color = 'salmon') +
        geom_point(aes( y = Mwave.max), size = 1, color = "dodgerblue")+
        geom_smooth(aes (y = Mwave.max), 
                    method = stats::loess, 
                    formula = y ~ x,
                    se = FALSE, color = 'dodgerblue') +
        xlab("Stimulation") + 
        ylab("EMG Signals") +
        ggtitle("Recruitment Curves  (H-reflex, M-wave)")

Sweeps sorted by stimulus intensity