# INSTRUCTIONS# 0. Create a folder named "emg" and another named "cycles". The EMG and Cycles files within must have exactly the same name.# 1. There must be n EMG files inside the folder named "emg". These data must have time in the first column and the name of the muscles in the first row, e.g. # # ## time ME MA FL RF VM VL ST # # ## 1 0.014 0.201416 -6.445313 22.65930 -0.100708 -0.906372 7.351685 -1.309204 # # ## 2 0.015 -2.316284 -0.100708 24.16992 1.812744 -1.913452 -4.531860 2.920532 # # ## 3 0.016 -7.351685 -7.150269 23.46497 0.704956 -5.337524 3.424072 -0.604248 # # ## 4 0.017 -5.538940 -3.222656 27.49329 5.236816 -4.330444 -1.611328 0.503540 # # ## 5 0.018 -10.675049 -5.740356 23.16284 -0.704956 2.014160 1.007080 -2.719116 # # ## 6 0.019 -12.487793 -3.927612 19.94019 2.014160 -5.136108 -0.805664 0.000000 # # 2. Subsequently, there should be n stride cycle files associated with the EMG table, with equivalent time data, in the following format: # ## V1 V2 # # ## 1 1.414 2.074 # # ## 2 2.448 3.115 # # ## 3 3.488 4.141 # # ## 4 4.515 5.168 # # ## 5 5.549 6.216 # # ## 6 6.596 7.249 # # Each row represents a cycle (stride)# To create V1 and V2, explicitly consider the following (based on the previous data):# # V1 # V2# # 1 1.414 (heel strike) 2.074 (midstance) Takeoff is not considered. In the next row of V1, consider the next heel strike.# To calculate V2, you can obtain the duty factor (%), then calculate the total cycle duration, multiply the cycle time by the duty factor, and finally add this result to V1. If you're still unsure about the value of V1, note that it is the first value of the step. # 3. Within the RStudio environment, go to the lower right panel and navigate to "File". Find the folder where your EMG and Cycles files are located. Click on the gear next to "Rename" and select "Set as Working Directory". You should do this every time you restart the program or want to change folders. # 4. The author of this routine suggests that the data be in seconds.# 5. You can modify the plot time in the variable plot_time.# Questions for Leonardo Lagos H library(musclesyneRgies) # ----------------------------------------------------------------- #--TRAE ARCHIVOS DESDE CARPETAS EMG Y CYCLES Y CREA UN raw_data_from_files-- #-------------------------------------------------------- #Get current working directory data_path <- getwd() data_path <- paste0(data_path, .Platform$file.sep) raw_data_from_files <- rawdata( path_cycles = paste0(data_path, "/cycles/"), path_emg = paste0(data_path, "/emg/"), header_cycles = TRUE ) # ----------------------------------------------------------------- #--RECORTA CICLOS DE LOS EXTREMOS--- #______________________________________________________________ # Say you recorded more cycles than those you want to consider for the analysis # You can subset the raw data (here we keep only 3 cycles, starting from the first) # RAW_DATA_subset <- lapply( # raw_data_from_files, # function(x) { # subsetEMG(x, # cy_max = 3, # cy_start = 1 # ) # } # ) # ----------------------------------------------------------------- #--GRAFICA DATOS BRUTOS- #-------------------------------------------------------- # Raw EMG can be plotted with the following (the first three seconds are plot by default) # Now also in dark mode if you fancy it # pp <- plot_rawEMG(raw_data_from_files[[1]], # trial = names(raw_data_from_files)[1], # plot_time= 1.5, # row_number = 4, # col_number = 4, # dark_mode = TRUE, # line_col = "tomato3" # ) #trace(plot_rawEMG, edit = TRUE) #editar fila Y # ----------------------------------------------------------------- #--APLICACIÓN FILTRADO SIN MODIFICACIÓN DE PARAMETROS-- #-------------------------------------------------------- # The raw EMG data set then needs to be filtered # If you don't want to subset the data set, just filter it as it is # Here we filter the whole data set with the default parameters for locomotion: # - Demean EMG # - High-pass IIR Butterworth 4th order filter (cut-off frequency 50 Hz) # - Full-wave rectification (default) # - Low-pass IIR Butterworth 4th order filter (cut-off frequency 20 Hz) # - Minimum subtraction # - Amplitude normalisation # trace(filtEMG, edit = TRUE) # #filtered_EMG <- lapply(raw_data_from_files, function(x) filtEMG(x)) # # # # ----------------------------------------------------------------- # #--GRAFICA CON FILTRO ANTERIOR-- # #-------------------------------------------------------- # pp <- plot_rawEMG(filtered_EMG[[1]], # trial = names(filtered_EMG)[1], # plot_time= 1.5, # row_number = 4, # col_number = 4, # dark_mode = TRUE, # ) # ----------------------------------------------------------------- #--APLICACIÓN FILTRADO "CON" POSIBILIDAD MODIFICACIÓN DE PARAMETROS-- #-------------------------------------------------------- # If you decide to change filtering parameters, just give them as arguments: another_filtered_EMG <- lapply( raw_data_from_files, function(x) { filtEMG(x, demean = TRUE, rectif = "fullwave", HPf = 30, #https://www.nature.com/articles/s41598-020-65257-w HPo = 4, LPf = 6, #https://www.nature.com/articles/s41598-020-65257-w, #https://www.frontiersin.org/articles/10.3389/fnhum.2015.00048/full LPo = 4, min_sub = TRUE, ampl_norm = TRUE ) } ) # # # trace(plot_rawEMG, edit = TRUE) #--# ----------------------------------------------------------------- #--GRAFICA CON PARAMETROS MODIFICABLES # #-------------------------------------------------------- # pp <- plot_rawEMG(another_filtered_EMG[[1]], # trial = names(another_filtered_EMG)[1], # plot_time= 1.5, # row_number = 4, # col_number = 4, # dark_mode = TRUE, # line_col = "tomato3" # ) # ----------------------------------------------------------------- #--NOMALIZACIÓN CON SEPARACIÓN CICLO EN 100 Y 100 PUNTOS-- #-------------------------------------------------------- # Now the filtered EMG needs some time normalisation so that cycles will be comparable # Here we time-normalise the filtered EMG, including only three cycles and trimming first # and last to remove unwanted filtering effects # Each cycle is divided into two parts, each normalised to a length of 100 points #trace(norm_EMG, edit = TRUE) norm_EMG <- lapply( another_filtered_EMG, function(x) { normEMG(x, trim = FALSE, cy_max = 2, #indica el numero de ciclos, en este caso 2 ciclos, con divisiones de 200 y 200 cycle_div = c(100, 100) ) } ) # ----------------------------------------------------------------- #--GRAFICA CON NORMALIZACION 100 Y 100 PUNTOS-- #-------------------------------------------------------- # The filtered and time-normalised EMG can be plotted with the following #trace(plot_meanEMG, edit = TRUE) pp <- plot_meanEMG( norm_EMG[[1]], trial = names(norm_EMG)[1], row_number = 4, col_number = 4, dark_mode = FALSE, line_size = 0.8, line_col = "tomato3" ) # ----------------------------------------------------------------- #--NORMALIZACIÓN CON DATA DISTRIBUIDA 60% Y 40%-- #-------------------------------------------------------- # If this cycle division does not work for you, it can be changed # But please remember to have the same amount of columns in the cycle times as the number # of phases you want your cycles to be divided into # Here we divide each cycle with a ratio of 60%-40% and keep only two cycles (first and last # are still trimmed, so to have two cycles you must start with at least four available) # another_norm_EMG <- lapply( # filtered_EMG, # function(x) { # normEMG(x, # trim = TRUE, # cy_max = 2, # cycle_div = c(120, 80) # ) # } # ) # ----------------------------------------------------------------- #--GRAFICA CON NORMALIZACION 60%-40% PUNTOS-- #-------------------------------------------------------- # The filtered and time-normalised EMG can be plotted with the following # pp <- plot_meanEMG( # another_norm_EMG[[1]], # trial = names(another_norm_EMG)[1], # row_number = 4, # col_number = 4, # dark_mode = TRUE, # line_size = 0.8, # line_col = "tomato3" # ) # --------------------------------------------------------------------------- # ----------------------------------------------------------------- #--APLICACIÓN NMF PARA EXTRACCION SINERGIA-- #-------------------------------------------------------- # At this stage, synergies can be extracted # This is the core function to extract synergies via NMF SYNS <- lapply(norm_EMG, synsNMF) # ----------------------------------------------------------------- #--GRAFICA CON SINERGIAS EXTRAIDAS #-------------------------------------------------------- # The extracted synergies can be plotted with the following #trace(plot_syn_trials, edit = TRUE) pp <- plot_syn_trials( SYNS[[1]], max_syns = max(unlist(lapply(SYNS, function(x) x$syns))), trial = names(SYNS)[1], dark_mode = TRUE, line_size = 0.8, line_col = "tomato1", sd_col = "tomato4", #limits= c(0.4,0.6), ) # Load data from text file on desktop file_path <- "~/Desktop/mi_carpeta/6.-(H)Motor_Primitive.txt" # Replace with your file path act_pattern <- read.table(file_path, header = TRUE) # Extract activation signal columns act_cols <- 7:ncol(act_pattern) act_signal <- as.matrix(act_pattern[, act_cols]) # Calculate FWHM and CoA for the first cycle first_cycle_signal <- act_signal[1:which.max(act_pattern$time), 1] act_sub_FWHM <- FWHM(first_cycle_signal) act_sub_CoA <- CoA(first_cycle_signal) # Plot activation signal, FWHM, and CoA hm <- min(first_cycle_signal) + (max(first_cycle_signal) - min(first_cycle_signal)) / 2 hm_plot <- first_cycle_signal hm_plot[which(hm_plot > hm)] <- hm hm_plot[which(hm_plot < hm)] <- NA plot(first_cycle_signal, type = "l", xlab = "Time", ylab = "Amplitude") lines(hm_plot, lwd = 3, col = 2) # FWHM (horizontal, in red) abline(v = act_sub_CoA, lwd = 3, col = 4) # CoA (vertical, in blue) # Calculate HFD and H for full activation signal nonlin_HFD <- HFD(act_signal)$Higuchi nonlin_H <- Hurst(act_signal, min_win = max(act_pattern$time))$Hurst message("Higuchi's fractal dimension: ", round(nonlin_HFD, 3)) message("Hurst exponent: ", round(nonlin_H, 3)) library(ggplot2) # Crear data frame con ángulos y valores df <- data.frame( angle = seq(0, 2*pi, length.out = 100), value = seq(0, 100, length.out = 101)[-1] ) # Crear gráfica circular con líneas desde el centro ggplot(df, aes(x = 0, y = 0)) + geom_segment(aes(xend = sin(angle)*value/100, yend = cos(angle)*value/100)) + coord_equal() + theme_void() + labs(title = "Gráfico circular con líneas desde el centro", x = "", y = "") + theme(plot.title = element_text(size = 20, hjust = 0.5)) library(ggplot2) # Crear data frame con ángulos y valores df <- data.frame( angle = seq(0, 2*pi, length.out = 300), value = c(20, 32, 87) ) # Crear gráfica circular con líneas desde el centro ggplot(df, aes(x = 0, y = 0)) + geom_segment(aes(xend = sin(angle)*value/100, yend = cos(angle)*value/100), size = 3) + coord_equal() + theme_void() + labs(title = "Gráfico circular con líneas desde el centro", x = "", y = "") + theme(plot.title = element_text(size = 20, hjust = 0.5)) # Now synergies don't have a functional order and need classification # Let's load the built-in data set to have some more trial to classify # (clustering cannot be done on only one trial and having just a few, # say less than 10, won't help) data("SYNS") # Classify with k-means and produce a plot that shows how the clustering went with: # - Full width at half maximum on the x-axis # - Centre of activity on the y-axis # (both referred to the activation patterns of the classified muscle synergies) SYNS_classified <- classify_kmeans(SYNS) # Classified synergies can be finally plotted with pp <- plot_classified_syns( SYNS_classified, dark_mode = TRUE, line_col = "tomato1", sd_col = "tomato4", condition = "TW" ) # "TW" = Treadmill Walking, change with your own # ----------------------------------------------------------------------- # A 2D UMAP plot of the classified synergies can be obtained with pp <- plot_classified_syns_UMAP( SYNS_classified, condition = "TW" ) # From now on, it's all about the analysis # For example, one can measure the full width at half maximum (FWHM) # of the activation patterns or their centre of activity (CoA) # Load a typical activation pattern of 30 cycles (from locomotion) data("act_pattern") # Reduce activation pattern to the first cycle act_sub <- act_pattern$signal[1:which(act_pattern$time == max(act_pattern$time))[1]] # Calculate FWHM of the first cycle act_sub_FWHM <- FWHM(act_sub) # Calculate CoA of the first cycle act_sub_CoA <- CoA(act_sub) # Half maximum (for the plots) hm <- min(act_sub) + (max(act_sub) - min(act_sub)) / 2 hm_plot <- act_sub hm_plot[which(hm_plot > hm)] <- hm hm_plot[which(hm_plot < hm)] <- NA # Plots plot(act_sub, ty = "l", xlab = "Time", ylab = "Amplitude") lines(hm_plot, lwd = 3, col = 2) # FWHM (horizontal, in red) graphics::abline(v = act_sub_CoA, lwd = 3, col = 4) # CoA (vertical, in blue) # Or perhaps one might want to investigate the nonlinear behaviour of a long activation pattern act <- act_pattern$signal # Calculate the local complexity or Higuchi's fractal dimension (HFD) nonlin_HFD <- HFD(act)$Higuchi # Calculate the global complexity or Hurst exponent (H) nonlin_H <- Hurst(act, min_win = max(act_pattern$time))$Hurst message("Higuchi's fractal dimension: ", round(nonlin_HFD, 3)) ## Higuchi's fractal dimension: 1.047 message("Hurst exponent: ", round(nonlin_H, 3)) ## Hurst exponent: 0.338