TREMSENToolbox.R

# -------------------------------- TREMSEN toolbox  ----------------------------------
# Version: 2.0
# Date: March 6, 2019
# Latest version available @ https://github.com/NIATS-UFU/TREMSEN-Toolbox.git

# Author: Prof. Adriano de Oliveira Andrade
# Contact email: adriano@ufu.br
# CV LATTES: http://lattes.cnpq.br/1229329519982110
# ORCID ID: http://orcid.org/0000-0002-5689-6606
# Researcher ID: http://www.researcherid.com/rid/D-9721-2012 
# GOOGLE ACADEMIC: https://scholar.google.com.br/citations?user=8nHKQHMAAAAJ&hl=pt-BR


# Adddress: Centre for Innovation and Technology Assessment in Health, 
#           Postgraduate Program in Electrical and Biomedical Engineering, 
#           Faculty of Electrical Engineering, Federal University of Uberlândia, 
#           Uberlândia, Brazil
# Webpage:  http://www.niats.feelt.ufu.br/

# Description: toolbox for processing data collected with the device TREMSEN
# ------------------------------------------------------------------------------------

# Install and load R packages
installAndLoadPackages <- function(requiredPackages)
{
  remainingPackages <- requiredPackages[!(requiredPackages %in% installed.packages()[,"Package"])];
  
  if(length(remainingPackages)) 
  {
    install.packages(remainingPackages);
  }
  for(packageName in requiredPackages)
  {
    library(packageName, character.only=TRUE, quietly=TRUE);
  }
}

# List of required packages
requiredPackages = c('readxl','mvnormtest', 'ggplot2','Hmisc',
'leaps','beanplot', 'moments', 'fBasics','lawstat','plotly',
   'robust','mclust', 'plyr','tsne', 'boot', 'reshape2',
 'pracma', 'seewave', 'psd', 'rlist', 'ggpubr', 'gridExtra', 
'grid', 'outliers', 'EMD','openxlsx',  'dygraphs','htmltools', 'tibble')

# Load and install packages if necessary
installAndLoadPackages(requiredPackages)


# Visualizacao de packages instalados
# search()


# LoadTREMSENFile ---------------------------------------------------------
# Load the data from the file `Filename` collected with TREMSEN 
# Input: 
#       Filename: file name, including path. This file should be a text file  
#                 generated by the data acquisition software of TREMSEN
# Output:
#       A dataframe with the loaded data
#
# Example of use:
#
#       testFilename <- file.choose()
#       df <- LoadTREMSENFile(testFilename)

LoadTREMSENFile <- function(Filename) {
  # Carregamento de arquivo do TREMSEN
  df = read.table(
    Filename,
    skip = 1,
    header = TRUE,
    sep = "\t",
    row.names = NULL,
    allowEscapes = TRUE
  )
  
  # foi necessário remover uma coluna pois o R não interpretou um caracter do arquivo original
  colnames(df) <- colnames(df)[2:ncol(df)] 
  df <- df[1:ncol(df) - 1]
  df[[1]] = as.numeric(df[[1]])
  
  
  # incluindo coluna que combina a resposta binária dos pulsos A e B
  X.PULSE <- combinePulseAB(df$X.PULSE.A., df$X.PULSE.B.)
  X.PULSE <- movavg(X.PULSE, n=10, type='m')
  
  indx1 <- which(X.PULSE >= 0.5)
  indx2 <- which(X.PULSE < 0.5)
  X.PULSE[indx1] <- 1
  X.PULSE[indx2] <- 0
  
  
  wnd <- as.numeric(vector(length = length(X.PULSE)))
  
  j <- 1
  FLAG <- FALSE
  
  for (i in 1:length(X.PULSE)) {
    if (X.PULSE[i] == 1) {
      FLAG <- TRUE
      wnd[i] <- as.numeric(c(paste(j, sep = "")))
    }
    else
    {
      if (FLAG == TRUE) {
        j <- j + 1
      }
      
      FLAG = FALSE
    }
  }
  
  X.PULSE.LABEL <- factor(wnd)
  
  df <- cbind(df, as.data.frame(X.PULSE))
  df <- cbind(df, as.data.frame(X.PULSE.LABEL))
  
  return(df)
}

# combinePulseAB ---------------------------------------------------------
#   Combination of pulses A and B generated by TREMSEN. 
#   The combination is based on the application of the function
#   OR to the input pulses A and B
# 
# Input: 
#       PulseA: pulse A recorded in the data file of TREMSEN
#       PulseB: pulse B recorded in the data file of TREMSEN
# Output:
#       Binary pulse resulting from the combination of pulses A and B
#
# Example of use:
#
#       See example of use in the function LoadTREMSENFile

combinePulseAB <- function(PulseA, PulseB) {
  th <- 2

  ppA <- vector(length = length(PulseA))
  ppB <- vector(length = length(PulseB))
  
  ppA[which(PulseA < th)] <- 0
  ppA[which(PulseA > th)] <- 1
  
  ppB[which(PulseB > th)] <- 0
  ppB[which(PulseB < th)] <- 1
  
  pp <- as.numeric(ppA | ppB)
  
  return(pp)
}

# detrendTremsenData ---------------------------------------------------------
#   Linear detrend. The function uses the function detrend available in the Package 'pracma' 
#   (Practical Numerical Math Functions) for removing piecewise linear trends of the data
# 
# Input: 
#       df: dataframe resulting from LoadTREMSENFile
# Output:
#       Linear detrended dataframe
#
# Example of use:
#
#       df.detrended <- detrendTremsenData(df)

detrendTremsenData <- function(df, startColRange=2, endColRange=39){
  
  X <- data.matrix(df)
  dX <- detrend(X[, c(startColRange:endColRange)], 'linear')
  
  res <- (as.data.frame(dX))
  ss <- cbind(df["X.Time."], res)
  ss <- cbind(ss, df[c("X.PULSE.A.","X.PULSE.B.", "X.PULSE")])
  
  return(ss)
}

# nonLineardetrendTremsenData ---------------------------------------------------------
#   Nonlinear detrend. The function removes nonlinear trends from the data. 
#   The data are first smoothed (see Tukey's smoothers, in  the package stats, function smooth) 
#   and then the resulting signal is subtracted from a nonlinear trend.
#   The nonlinear trend is estimated by Local Polynomial Regression Fitting (see function loess)
#   in the package stats. The resulting signal is also linearly detrended.
#
# Input: 
#       df: dataframe resulting from LoadTREMSENFile
# Output:
#       Nonlinear detrended dataframe
#
# Example of use:
#
#       df.nonlineardetrended <- nonLineardetrendTremsenData(df)

nonLineardetrendTremsenData <- function(df, startColRange=2, endColRange=39) {
  
  df.smoothed <- smoothTremsenData(df)
  df.loess <- loessTremsenData(df.smoothed)
  
  res <- df
  res[startColRange:endColRange] <- df.smoothed[startColRange:endColRange] - df.loess[startColRange:endColRange]
  res[startColRange:endColRange] <- detrend(as.matrix(res[startColRange:endColRange]), 'linear')
  
  return(res)
}

# smoothTremsenData ---------------------------------------------------------
#   Smooth data based on the Tukey's (Running Median) Smoothing (see the package 
#   smooth in the package stats)
# 
# Input: 
#       df: dataframe resulting from LoadTREMSENFile
# Output:
#       smoothed data
#
# Example of use:
#
#       df.smoothed <- smoothTremsenData(df)

smoothTremsenData <- function(df, startColRange=2, endColRange=39) {
  
  X <- data.matrix(df)
  
  dx <- apply(X[,c(startColRange:endColRange)], 2, smooth,kind="3RS3R", twiceit = TRUE)
  
  res <- (as.data.frame(dx))
  ss <- cbind(df["X.Time."], res)
  ss <- cbind(ss, df[c("X.PULSE.A.","X.PULSE.B.","X.PULSE")])
  
  return(ss)
}

# loessTremsenData ---------------------------------------------------------
#       Estimate nonlinear trend
# Input: 
#       df: dataframe resulting from LoadTREMSENFile
# Output:
#       estimate nonlinear trend in a data
#
# Example of use:
#
#       df.nonlineardetrended <- nonLineardetrendTremsenData(df)

loessTremsenData <- function(df, startColRange=2, endColRange=39) {
  
  X <- data.matrix(df)
  
  loessData <- function(yamp,t) {
    dat <- data.frame(x = t, y = yamp)
    yp = predict(loess(y ~ x, dat, span = 0.1))
    return(yp)
  }
  
  dx <- apply(X[,c(startColRange:endColRange)], 2, loessData, t = X[,c(1)])
  
  res <- (as.data.frame(dx))
  ss <- cbind(df["X.Time."], res)
  ss <- cbind(ss, df[c("X.PULSE.A.","X.PULSE.B.","X.PULSE")])
  
  return(ss)
}

# windowTremsenData ---------------------------------------------------------
#       Apply a rectangular window to the input signal. The input signal is multiplied by PULSE (0=LOW and 1=HIGH)
# Input: 
#       df: dataframe resulting from LoadTREMSENFile
# Output:
#       windowed data
#
# Example of use:
#
#       df.windowed <- windowTremsenData(df)

windowTremsenData <- function(df, startColRange=2, endColRange=39) {
  
  M <- apply(df[,c(startColRange:endColRange)], 2, 
             wnd <- function(x,y) { return(x*y) }, 
             y = df$X.PULSE)
  
  dxx <- df
  dxx[,c(startColRange:endColRange)] <- M
  return(dxx)
}

# psdTremsenData ---------------------------------------------------------
#             Estimate the power spectrum of the input data
# Input: 
#       df: dataframe resulting from LoadTREMSENFile
# Output:
#       power spectrum based on the 'Adaptive sine multitaper power spectral density estimation'. See the function
#       pspectrum in the package psd
#
# Example of use:
#
#       pp <- psdTremsenData(df.nonlineardetrended) 

psdTremsenData <- function(df, startColRange=2, endColRange=39) {
  
  psf <- function(vec, fs) {
    
    # library(psd)
    sss <- pspectrum(vec, verbose = FALSE, niter = 10, 
                     AR = TRUE, x.frqsamp = fs, plot = FALSE) 
    
    return (data.frame("freq" = sss$freq, "spec" = sss$spec))
  }
  
  
  X <- data.matrix(df)
  
  fs <- 1 / (df$X.Time.[2] - df$X.Time.[1])
  
  Nwindows <- nlevels(df$X.PULSE.LABEL)-1
  
  dx <- list()
  
  for (i in 1:Nwindows){
    
    indx <- which(df$X.PULSE.LABEL==i)
    
    dx[[i]] <- apply(X[indx, c(startColRange:endColRange)], 2, psf, fs=fs)
  }
  
  return(dx)
}

# featExtractFromTremenDataSet ---------------------------------------------------------
#             Estimate features from the input data
# Input: 
#       df: dataframe resulting from LoadTREMSENFile
#       w: length of the window, from which features are estimated
#       s: spot step. Windows of length w are positioned on each spot, starting from 1, being then incremented by s, 
#          to the last spot, which should be less or equal to (total - ww), where total is the number of samples of
#          the time-series
#
# Output:
#       set of estimated features for each window
#
# Example of use:
#
#       df.featTremsenData <- featExtractFromTremenDataSet(df.nonlineardetrended,w=50,s=10, method = "rms")

featExtractFromTremenDataSet <- function(df, w, s, method="rms", startColRange=2, endColRange=39) {
  
  methods <- c("rms", "mav", "peak", "mavfd", "mavfdn","mavsd", "mavsdn")
  
  meth <- pmatch(method, methods)
  
  if (is.na(meth)) 
    stop("invalid feature extraction method")
  else
    selectedMethod <- methods[meth]
  
  slideFunct <- function(data, ww, ss) {
    total <- length(data)
    spots <- seq(from = 1,
                 to = (total - ww),
                 by = ss)
    result <- vector(length = length(spots))
    
    for (i in 1:length(spots)) {
      
      # Estima a raíz do valor médio quadrado
      if (selectedMethod == "rms")
        result[i] <- sqrt(mean(data[spots[i]:(spots[i] + ww-1)]^2))
      
      # Estima o valor absoluto médio de um vetor.
      if (selectedMethod == "mav")
        result[i] <- sum(abs(data[spots[i]:(spots[i] + ww-1)]))/ww
      
      # Estima o valor máximo de um vetor, considera somente valores positivos da janela
      if (selectedMethod == "peak")
        result[i] <- max(data[spots[i]:(spots[i] + ww-1)])
      
      # Estima a média do valor absoluto da primeira diferença
      if (selectedMethod == "mavfd")
        result[i] <- sum(abs(diff(data[spots[i]:(spots[i] + ww-1)]) ))/(ww-1)
      
      # Estima a média do valor absoluto da primeira diferença do sinal normalizado
      if (selectedMethod == "mavfdn") {
        dd <- data[spots[i]:(spots[i] + ww-1)]
        dd <- (dd-mean(dd)) / std(dd)
        result[i] <- sum(abs(diff(dd))) / (ww-1)
      }
      
      # Estima a média do valor absoluto da segunda diferença 
      if (selectedMethod == "mavsd")
        result[i] <- sum(abs(diff(data[spots[i]:(spots[i] + ww-1)], lag = 2) )) / (ww-2)
      
      # Estima a média do valor absoluto da segunda diferença do sinal normalizado
      if (selectedMethod == "mavsdn") {
        dd <- data[spots[i]:(spots[i] + ww-1)]
        dd <- (dd - mean(dd)) / std(dd)
        result[i] <- sum(abs(diff(dd,lag = 2))) / (ww-2)
      }

    }
    
    return(result)
  }
  
  X <- data.matrix(df)
  dx <- apply(X[,c(startColRange:endColRange)], 2, slideFunct, ww = w, ss = s)
  
  dt <- df$X.Time.[2] - df$X.Time.[1]
  
  tt <- dt * (seq(from = 1, to = (length(df$X.Time.) - w), by = s) + w/2)
  
  PulseA <- approx(df$X.Time., df$X.PULSE.A., xout = tt, method = "linear")
  PulseB <- approx(df$X.Time., df$X.PULSE.B., xout = tt, method = "linear")
  
  res <- (as.data.frame(dx))
  ss <- cbind(data.frame(X.Time.= tt), res)
  ss <- cbind(ss,data.frame(X.PULSE.A.= PulseA$y))
  ss <- cbind(ss,data.frame(X.PULSE.B.= PulseB$y))
  
  # incluindo coluna que combina a resposta binária dos pulsos A e B
  pp <- combinePulseAB(ss$X.PULSE.A., ss$X.PULSE.B.)
  ss <- cbind(ss,data.frame(X.PULSE = pp))
  
  return(ss)
}

# getStatisticsFromWindowedTremenDataSet ---------------------------------------------------------
#       Estimate statistics from a set of features
# Input: 
#       df: dataframe resulting from LoadTREMSENFile
#       f: statistic
#
# Output:
#       estimated statistic
#
# Example of use:
#
#       medianRMS <- getStatisticsFromWindowedTremenDataSet(df.featTremsenData, f=median)

getStatisticsFromWindowedTremenDataSet <- function(df, f = median, startColRange=2, endColRange=39) {
  
  X <- df
  
  xx <- X$X.PULSE
  
  # condicao em que o pulso começa em nível alto
  if (xx[1] == 1) 
    xx[1] <- 0 
  
  ss <- abs(diff(xx))
  Nwindows <- sum(ss == 1)/2 
  
  boundwnd <- which(ss == 1)
  boundwnd[seq(1,length(boundwnd),2)] <- boundwnd[seq(1,length(boundwnd),2)] + 1
  
  result <- vector()
  
  windx <- 1
  
  for (n in seq(1,length(boundwnd),2)) 
  {
    indxo <- boundwnd[n]
    indxf <- boundwnd[n+1]
    name <- paste('window:', windx, sep='')
    windx <- windx + 1
    tmp <- list (apply(X[indxo:indxf,startColRange:endColRange], 2, FUN=f))
    result[[name]] <- tmp
  }
  
  return(result)
}

# resampleTremsenData ---------------------------------------------------------
#       resample data set
# Input: 
#       df: dataframe resulting from LoadTREMSENFile
#       fs: new sampling frequency in Hz
#
# Output:
#       resampled data set
#
# Example of use:
#
#       df.resampled <- resampleTremsenData(df.nonlineardetrended,200)

resampleTremsenData <- function(df, fs) {
  
  reseampleData <- function(y,t,xx) {
    yp <- spline(t, y, xout = xx, method = "fmm")
    return(yp$y)
  }
  
  reseampleDataLinear <- function(y,t,xx) {
    yp <- approx(t, y, xout = xx, method = "linear")
    return(yp$y)
  }
  
  
  X <- data.matrix(df)
  
  tnew <- seq(from = df$X.Time.[1], to = df$X.Time.[length(df$X.Time.)], by = 1/fs)
  
  dx <- apply(X[,c(2:39)], 2, reseampleData, t=df$X.Time., xx = tnew)
  dx1 <- apply(X[,c(40:41)], 2, reseampleDataLinear, t=df$X.Time., xx = tnew)
  
  res <- (as.data.frame(dx))
  res1 <- (as.data.frame(dx1))
  
  ss <- cbind(data.frame(X.Time.= tnew), res)
  ss <- cbind(ss, res1)
  
  # incluindo coluna que combina a resposta binária dos pulsos A e B
  pp <- combinePulseAB(ss$X.PULSE.A., ss$X.PULSE.B.)
  ss <- cbind(ss,data.frame(X.PULSE = pp))
  
  return(ss)
}

# plotPSTremsenDataSet ---------------------------------------------------------
#             plot the power spectrum of a data set
# Input: 
#       pp: power spectrum, as estimated by psdTremsenData
#       printplot: boolean flag (if true, then it will plot the power spectrumm, otherwise a figure handle is generated)
#
# Output:
#       figure displaying the power spectrum of each time series in the data set
#
# Example of use:
#
#  pp <- psdTremsenData(df.nonlineardetrended) ## It is a good practice to remove trends prior to use this function
#  g1<-plotPSTremsenDataSet(pp[[1]],printplot = FALSE)
#  print(g1)


plotPSTremsenDataSet <- function(pp, printplot = TRUE) {
  
  # Lookup table for changing the label of the graphs:
  lookupTable <- c(
    
    X.G1.X. = "G1X",
    X.G1.Y. = "G1Y",
    X.G1.Z. = "G1Z",
    X.G2.X. = "G2X",
    X.G2.Y. = "G2Y",
    X.G2.Z. = "G2Z",
    X.G3.X. = "G3X",
    X.G3.Y. = "G3Y",
    X.G3.Z. = "G3Z",
    X.G4.X. = "G4X",
    X.G4.Y. = "G4Y",
    X.G4.Z. = "G4Z",
    
    X.A1.X. = "A1X",
    X.A1.Y. = "A1Y",
    X.A1.Z. = "A1Z",
    X.A2.X. = "A2X",
    X.A2.Y. = "A2Y",
    X.A2.Z. = "A2Z",
    X.A3.X. = "A3X",
    X.A3.Y. = "A3Y",
    X.A3.Z. = "A3Z",
    X.A4.X. = "A4X",
    X.A4.Y. = "A4Y",
    X.A4.Z. = "A4Z",
    
    X.M1.X. = "M1X",
    X.M1.Y. = "M1Y",
    X.M1.Z. = "M1Z",
    X.M2.X. = "M2X",
    X.M2.Y. = "M2Y",
    X.M2.Z. = "M2Z",
    X.M3.X. = "M3X",
    X.M3.Y. = "M3Y",
    X.M3.Z. = "M3Z",
    X.M4.X. = "M4X",
    X.M4.Y. = "M4Y",
    X.M4.Z. = "M4Z",
    
    X.PULSE = "PULSE",
    X.EMG1. = "EMG1",
    X.EMG2. = "EMG2"
    
  )

  #pp should be estimated with psdTremsenData
  dat <- melt(pp, id = c("freq","spec"))
  facs <- factor(dat$L1)
  
  c1 <- c("X.G1.X.","X.G1.Y.","X.G1.Z.","X.G2.X.","X.G2.Y.","X.G2.Z.")
  c2 <- c("X.A1.X.","X.A1.Y.","X.A1.Z.","X.A2.X.","X.A2.Y.","X.A2.Z.")
  c3 <- c("X.M1.X.","X.M1.Y.","X.M1.Z.","X.M2.X.","X.M2.Y.","X.M2.Z.")
  
  g <- ggplot() 
  g <- g + geom_line(data=dat[which(facs == c1),], alpha = 0.8, aes(x = freq, y = spec , group=L1), size=1)

  g <- g + geom_line(data=dat[which(facs == c2),], alpha = 0.8, aes(x = freq, y = spec , group=L1), size=1)
 
  g <- g + geom_line(data=dat[which(facs == c3),], alpha = 0.8, aes(x = freq, y = spec , group=L1), size=1)
  
  
  g <- g + scale_y_continuous(trans = "log10")
  
  g <- g + facet_grid(L1 ~ ., scales = "free_y", labeller = labeller(L1=lookupTable)) 
  
  g <- g + geom_vline(xintercept = c(2, 5, 10, 15), colour = "red")
  
  g <- g + xlab("frequency (Hz)") + ylab("energy")
  
  g <- g + theme_bw(12)
  
  if(printplot == TRUE) 
    print(g)
  
  return(list(g))
  
}

# plotTremsenDataset ---------------------------------------------------------
#       plot time-series in the data set
# Input: 
#       dflist: list of dataframes to be plotted
#       indxs: index of signals to be plotted from the data set
#       printplot: boolean flag (if true, then it will plot the power spectrumm, otherwise a figure handle is generated)
#
# Output:
#       figure displaying the data set
#
# Example of use:
#
# gg1 <- plotTremsenDataset(list(df.nonlineardetrended,df.resampled ),indxs=c(1,2),printplot=FALSE, droplvs = c("X.PULSE.A.","X.PULSE.B."), alphavec = c(1, 0.3))
# gg2 <- plotTremsenDataset(list(df.nonlineardetrended,df.resampled ),indxs=c(3,4),printplot=FALSE, droplvs = c("X.PULSE.A.","X.PULSE.B."), alphavec = c(1, 0.3))
# gg3 <- plotTremsenDataset(list(df.nonlineardetrended,df.resampled ),indxs=c(5,6),printplot=FALSE, droplvs = c("X.PULSE.A.","X.PULSE.B."), alphavec = c(1, 0.3))
# grid.arrange(gg1[[1]], gg2[[1]], gg3[[1]], nrow = 1, ncol=3, top = "Data set")

plotTremsenDataset <- function(dflist, indxs = c(1:13), printplot = TRUE, droplvs = "", alphavec = NULL) {
  
  c1 <- c("X.Time.","X.G1.X.","X.G1.Y.","X.G1.Z.", "X.PULSE","X.PULSE.A.","X.PULSE.B.")
  c2 <- c("X.Time.","X.G2.X.","X.G2.Y.","X.G2.Z.", "X.PULSE","X.PULSE.A.","X.PULSE.B.")
  c3 <- c("X.Time.","X.A1.X.","X.A1.Y.","X.A1.Z.","X.PULSE","X.PULSE.A.","X.PULSE.B.")
  c4 <- c("X.Time.","X.A2.X.","X.A2.Y.","X.A2.Z.","X.PULSE","X.PULSE.A.","X.PULSE.B.")
  c5 <- c("X.Time.","X.M1.X.","X.M1.Y.","X.M1.Z.","X.PULSE","X.PULSE.A.","X.PULSE.B.")
  c6 <- c("X.Time.","X.M2.X.","X.M2.Y.","X.M2.Z.","X.PULSE","X.PULSE.A.","X.PULSE.B.")
  
  c7 <- c("X.Time.","X.G3.X.","X.G3.Y.","X.G3.Z.", "X.PULSE","X.PULSE.A.","X.PULSE.B.")
  c8 <- c("X.Time.","X.G4.X.","X.G4.Y.","X.G4.Z.", "X.PULSE","X.PULSE.A.","X.PULSE.B.")
  c9 <- c("X.Time.","X.A3.X.","X.A3.Y.","X.A3.Z.","X.PULSE","X.PULSE.A.","X.PULSE.B.")
  c10 <- c("X.Time.","X.A4.X.","X.A4.Y.","X.A4.Z.","X.PULSE","X.PULSE.A.","X.PULSE.B.")
  c11 <- c("X.Time.","X.M3.X.","X.M3.Y.","X.M3.Z.","X.PULSE","X.PULSE.A.","X.PULSE.B.")
  c12 <- c("X.Time.","X.M4.X.","X.M4.Y.","X.M4.Z.","X.PULSE","X.PULSE.A.","X.PULSE.B.")
  
  c13 <- c("X.Time.","X.EMG1.","X.EMG2.","X.PULSE","X.PULSE.A.","X.PULSE.B.")
  
  hh <- list(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13)
  
  #seleciona um subconjunto
  hh <- list.subset(hh,indxs)
  
  heads <- levels(factor(unlist(hh)))
  
  heads <- levels(droplevels(factor(unlist(hh)), droplvs))
  
  if (isempty(alphavec) || length(dflist)!=length(alphavec)){
    aa <- rep(0.7,length(dflist)) #alpha
  }else{
    aa <- alphavec
  }
  
  g <- ggplot() 
  for(indx in 1:length(dflist)){
    g <- g + geom_line(data = melt(dflist[[indx]][heads], id = "X.Time."), 
                       aes(x = X.Time., y = value), alpha = aa[indx], 
                       color = indx)
    
    g <- g + labs(y = "")
    g <- g + labs(x = "time (s)")
    
  }
  
  # Lookup table for changing the label of the graphs:
  lookupTable <- c(
    
    X.G1.X. = "G1X\n(dps)",
    X.G1.Y. = "G1Y\n(dps)",
    X.G1.Z. = "G1Z\n(dps)",
    X.G2.X. = "G2X\n(dps)",
    X.G2.Y. = "G2Y\n(dps)",
    X.G2.Z. = "G2Z\n(dps)",
    X.G3.X. = "G3X\n(dps)",
    X.G3.Y. = "G3Y\n(dps)",
    X.G3.Z. = "G3Z\n(dps)",
    X.G4.X. = "G4X\n(dps)",
    X.G4.Y. = "G4Y\n(dps)",
    X.G4.Z. = "G4Z\n(dps)",
    
    X.A1.X. = "A1X\n(g)",
    X.A1.Y. = "A1Y\n(g)",
    X.A1.Z. = "A1Z\n(g)",
    X.A2.X. = "A2X\n(g)",
    X.A2.Y. = "A2Y\n(g)",
    X.A2.Z. = "A2Z\n(g)",
    X.A3.X. = "A3X\n(g)",
    X.A3.Y. = "A3Y\n(g)",
    X.A3.Z. = "A3Z\n(g)",
    X.A4.X. = "A4X\n(g)",
    X.A4.Y. = "A4Y\n(g)",
    X.A4.Z. = "A4Z\n(g)",
    
    X.M1.X. = "M1X\n(gauss)",
    X.M1.Y. = "M1Y\n(gauss)",
    X.M1.Z. = "M1Z\n(gauss)",
    X.M2.X. = "M2X\n(gauss)",
    X.M2.Y. = "M2Y\n(gauss)",
    X.M2.Z. = "M2Z\n(gauss)",
    X.M3.X. = "M3X\n(gauss)",
    X.M3.Y. = "M3Y\n(gauss)",
    X.M3.Z. = "M3Z\n(gauss)",
    X.M4.X. = "M4X\n(gauss)",
    X.M4.Y. = "M4Y\n(gauss)",
    X.M4.Z. = "M4Z\n(gauss)",
    
    X.PULSE = "PULSE",
    X.EMG1. = "EMG1",
    X.EMG2. = "EMG2"
    
  )
  
  g <- g + facet_grid(variable ~ ., scales = "free_y", labeller = labeller(variable = lookupTable))

  g <- g + theme_bw(18)
  #g+theme(axis.text=element_text(size=18), axis.title=element_text(size=18,face="bold"))
  if (printplot == TRUE) 
    print(g)
  
  return(list(g))
}


# Multiple plot function ---------------------------------------------------------
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols: Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

#' plotMultiPanelData ---------------------------------------------------------
#'
#' @param df1 -> data frame
#'
#' @return
#' @export
#'
#' @examples
#' Filename <- file.choose() # select file
#' df <- LoadTREMSENFile(Filename) # load tremsenfile
#' plotMultiPanelData(df) # plot tremsen data

plotMultiPanelData <- function(df1)
{
  colnames(df1)[1] <- "time"
  
  browsable(
    tagList(list(
      tags$div(
        style = 'width:33%;display:block;float:left;',
        dygraph(data.frame(time=df1$time, G1x=df1$X.G1.X., G2x=df1$X.G2.X.), group = "ensync", height = 200, width = "100%") %>%
          dyLegend(show="always")%>%dyOptions(colors = c("rgb(255,0,0)", "rgb(155,0,0)"), colorSaturation=c(0.5, 0.1)),
        dygraph(data.frame(time=df1$time, G1y=df1$X.G1.Y., G2y=df1$X.G2.Y.), group = "ensync", height = 200, width = "100%") %>%
          dyLegend(show="always")%>%dyOptions(colors = c("rgb(0,255,0)", "rgb(0,155,0)")),
        dygraph(data.frame(time=df1$time, G1z=df1$X.G1.Z., G2z=df1$X.G2.Z.), group = "ensync", height = 200, width = "100%") %>%
          dyLegend(show="always")%>%dyOptions(colors = c("rgb(0,0,255)", "rgb(0,0,155)"))%>%dyRangeSelector()
      ),
      tags$div(
        style = 'width:33%;display:block;float:left;',
        dygraph(data.frame(time=df1$time, A1x=df1$X.A1.X., A2x=df1$X.A2.X.), group = "ensync", height = 200, width = "100%") %>%
          dyLegend(show="always")%>%dyOptions(colors = c("rgb(255,0,0)", "rgb(155,0,0)")),
        dygraph(data.frame(time=df1$time, A1y=df1$X.A1.Y., A2y=df1$X.A2.Y.), group = "ensync", height = 200, width = "100%") %>%
          dyLegend(show="always")%>%dyOptions(colors =  c("rgb(0,255,0)", "rgb(0,155,0)")),
        dygraph(data.frame(time=df1$time, A1z=df1$X.A1.Z., A2z=df1$X.A2.Z.), group = "ensync", height = 200, width = "100%") %>%
          dyLegend(show="always")%>%dyOptions(colors = c("rgb(0,0,255)", "rgb(0,0,155)"))%>%dyRangeSelector()
      ),
      tags$div(
        style = 'width:33%;display:block;float:left;',
        dygraph(data.frame(time=df1$time, M1x=df1$X.M1.X., M2x=df1$X.M2.X.), group = "ensync", height = 200, width = "100%") %>%
          dyLegend(show="always")%>%dyOptions(colors = c("rgb(255,0,0)", "rgb(155,0,0)")),
        dygraph(data.frame(time=df1$time, M1y=df1$X.M1.Y., M2y=df1$X.M2.Y.), group = "ensync", height = 200, width = "100%") %>%
          dyLegend(show="always")%>%dyOptions(colors =  c("rgb(0,255,0)", "rgb(0,155,0)")),
        dygraph(data.frame(time=df1$time, M1z=df1$X.M1.Z., M2z=df1$X.M2.Z., df1$X.PULSE), group = "ensync", height = 200, width = "100%") %>%
          dyLegend(show="always")%>%dyOptions(colors = c("rgb(0,0,255)", "rgb(0,0,155)", "rgb(0,0,0)"))%>%dyRangeSelector()
        
      )
    )
    ))
}