Ploteador de pulsos para NAS (Nacho Acquisition System)

library (ggplot2)
library(reshape2)

Partimos de un archivo CSV que contiene una línea por cada disparo registrado.

MyData <- read.csv(file="capture.csv", header=FALSE, sep=",")
num_shots <- nrow(MyData)
samples_per_shot <- 500

Para este ejemplo hay 2000 disparos con 500 muestras por disparo. Incluimos el n´úmero de disparo en el conjunto original de datos para luego poder aplanarlos usando ese número como identificador. Aprovechamos también para poner nombres a las cosas.

MyData$shot_num=factor(1:num_shots)
colnames(MyData)[colnames(MyData)=="V1"] <- "data1"
colnames(MyData)[colnames(MyData)=="V2"] <- "data2"
colnames(MyData)[colnames(MyData)=="V3"] <- "data3"
colnames(MyData)[colnames(MyData)=="V4"] <- "data4"
colnames(MyData)[colnames(MyData)=="V5"] <- "time_stamp_1"
colnames(MyData)[colnames(MyData)=="V6"] <- "data5"
colnames(MyData)[colnames(MyData)=="V7"] <- "time_stamp_2"
colnames(MyData)[colnames(MyData)=="shot_num"] <- "shot_num"
colnames(MyData)[8:507] <- paste("S",1:samples_per_shot,sep="")
melted_data <- melt(MyData[,c(508,8:507)], id.vars = "shot_num")
colnames(melted_data)[colnames(melted_data)=="shot_num"] <- "shot_num"
colnames(melted_data)[colnames(melted_data)=="variable"] <- "sample_num"
colnames(melted_data)[colnames(melted_data)=="value"] <- "sample_value"
melted_data$sample_num <- as.numeric(melted_data$sample_num)

Tras reordenar los datos nos quedamos con melted_data que contiene tres variables: el número de disparo shot_num, el número de muestra de cada disparo sample_num y el valor de la muestra sample_value.

Ahora calculamos los valores m´áximos para cada disparo y hacemos un gráfico de cajas para ver si podemos eliminar cosas que no tengan pinta de tener un pulso.

Probamos a quedarnos con los disparos cuyo valor máximo esté dentro del segundo y tercer cuartil.

library(stats)
filter_limits <- quantile(maxs_per_shot$max_value)[c(2,3)]
low_limit = filter_limits[1] * 1.2
high_limit = filter_limits[2] * 1.1
filtered_shot_numbers <- maxs_per_shot[
  maxs_per_shot$max_value > low_limit  & 
  maxs_per_shot$max_value < high_limit, "shot_num"]
filtered_melted_data <- melted_data[melted_data$shot_num %in% filtered_shot_numbers[1:30],]

A mano he añadido unos coeficientes para modificar el máximo y el mínimo, de forma que se quedaran fuera cosas que aparentemente parecen ruido, aunque lo suyo sería hacer distintos grupos y analizarlos por separado no sea que sean cosas interesantes. Lo mismo con los pulsos cuyos máximos están en el cuarto cuartil; los primeros que he visto son ruido, pero lo mismo hay cosas interesantes.

Por último ploteamos las cosas en una matriz para vez los pulsos a vista de pájaro. He incluido dos líneas horizontales (en gris), con los límites inferior y superior del filtro que he usado para discriminar los disparos. Representamos tres gráficas con los primeros 60 disparos que pasaron el filtro.

library (ggplot2)
shot_subset <- melted_data[melted_data$shot_num %in% filtered_shot_numbers[1:20],]
ggplot(shot_subset, aes(x=sample_num, y=sample_value)) +
  geom_point(colour="blue", size=0) +
  geom_line(colour="blue", alpha=0.2) +
  xlab("Sample number") +
  ylab("Sample value") +
  geom_hline(yintercept=low_limit, color="gray") +
  geom_hline(yintercept=high_limit, color="gray") +
  facet_wrap(~shot_num)

shot_subset <- melted_data[melted_data$shot_num %in% filtered_shot_numbers[21:40],]
ggplot(shot_subset, aes(x=sample_num, y=sample_value)) +
  geom_point(colour="blue", size=0) +
  geom_line(colour="blue", alpha=0.2) +
  xlab("Sample number") +
  ylab("Sample value") +
  geom_hline(yintercept=low_limit, color="gray") +
  geom_hline(yintercept=high_limit, color="gray") +
  facet_wrap(~shot_num)

shot_subset <- melted_data[melted_data$shot_num %in% filtered_shot_numbers[41:60],]
ggplot(shot_subset, aes(x=sample_num, y=sample_value)) +
  geom_point(colour="blue", size=0) +
  geom_line(colour="blue", alpha=0.2) +
  xlab("Sample number") +
  ylab("Sample value") +
  geom_hline(yintercept=low_limit, color="gray") +
  geom_hline(yintercept=high_limit, color="gray") +
  facet_wrap(~shot_num)

Otra forma curiosa de verlo es apilando todos los 379 disparos que pasaron el filtro impuesto.

filtered_melted_data <- melted_data[melted_data$shot_num %in% filtered_shot_numbers,]
filtered_melted_data$shot_num <- paste("Shot:",filtered_melted_data$shot_num)
ggplot(filtered_melted_data, aes(x=sample_num, y=sample_value)) +
  geom_point(size=0, alpha=0.1) +
  geom_line(alpha=0.1) +
  xlab("Sample number") +
  ylab("Sample value") +
  geom_hline(yintercept=low_limit, color="gray") +
  geom_hline(yintercept=high_limit, color="gray") +
  theme(legend.position="none")

LS0tCnRpdGxlOiAiUGxvdGVhZG9yIGRlIHB1bHNvcyBwYXJhIE5BUyAoTmFjaG8gQWNxdWlzaXRpb24gU3lzdGVtKSIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQKICBwZGZfZG9jdW1lbnQ6IGRlZmF1bHQKLS0tCmBgYHtyfQpsaWJyYXJ5IChnZ3Bsb3QyKQpsaWJyYXJ5KHJlc2hhcGUyKQoKYGBgCgoKUGFydGltb3MgZGUgdW4gYXJjaGl2byBDU1YgcXVlIGNvbnRpZW5lIHVuYSBsw61uZWEgcG9yIGNhZGEgZGlzcGFybyByZWdpc3RyYWRvLgoKYGBge3J9Ck15RGF0YSA8LSByZWFkLmNzdihmaWxlPSJjYXB0dXJlLmNzdiIsIGhlYWRlcj1GQUxTRSwgc2VwPSIsIikKbnVtX3Nob3RzIDwtIG5yb3coTXlEYXRhKQpzYW1wbGVzX3Blcl9zaG90IDwtIDUwMApgYGAKClBhcmEgZXN0ZSBlamVtcGxvIGhheSBgciBudW1fc2hvdHNgIGRpc3Bhcm9zIGNvbiBgciBzYW1wbGVzX3Blcl9zaG90YCBtdWVzdHJhcyBwb3IgZGlzcGFyby4gSW5jbHVpbW9zIGVsIG7CtMO6bWVybyBkZSBkaXNwYXJvIGVuIGVsIGNvbmp1bnRvIG9yaWdpbmFsIGRlIGRhdG9zIHBhcmEgbHVlZ28gcG9kZXIgYXBsYW5hcmxvcyB1c2FuZG8gZXNlIG7Dum1lcm8gY29tbyBpZGVudGlmaWNhZG9yLiBBcHJvdmVjaGFtb3MgdGFtYmnDqW4gcGFyYSBwb25lciBub21icmVzIGEgbGFzIGNvc2FzLgoKYGBge3J9Ck15RGF0YSRzaG90X251bT1mYWN0b3IoMTpudW1fc2hvdHMpCmNvbG5hbWVzKE15RGF0YSlbY29sbmFtZXMoTXlEYXRhKT09IlYxIl0gPC0gImRhdGExIgpjb2xuYW1lcyhNeURhdGEpW2NvbG5hbWVzKE15RGF0YSk9PSJWMiJdIDwtICJkYXRhMiIKY29sbmFtZXMoTXlEYXRhKVtjb2xuYW1lcyhNeURhdGEpPT0iVjMiXSA8LSAiZGF0YTMiCmNvbG5hbWVzKE15RGF0YSlbY29sbmFtZXMoTXlEYXRhKT09IlY0Il0gPC0gImRhdGE0Igpjb2xuYW1lcyhNeURhdGEpW2NvbG5hbWVzKE15RGF0YSk9PSJWNSJdIDwtICJ0aW1lX3N0YW1wXzEiCmNvbG5hbWVzKE15RGF0YSlbY29sbmFtZXMoTXlEYXRhKT09IlY2Il0gPC0gImRhdGE1Igpjb2xuYW1lcyhNeURhdGEpW2NvbG5hbWVzKE15RGF0YSk9PSJWNyJdIDwtICJ0aW1lX3N0YW1wXzIiCmNvbG5hbWVzKE15RGF0YSlbY29sbmFtZXMoTXlEYXRhKT09InNob3RfbnVtIl0gPC0gInNob3RfbnVtIgpjb2xuYW1lcyhNeURhdGEpWzg6NTA3XSA8LSBwYXN0ZSgiUyIsMTpzYW1wbGVzX3Blcl9zaG90LHNlcD0iIikKCm1lbHRlZF9kYXRhIDwtIG1lbHQoTXlEYXRhWyxjKDUwOCw4OjUwNyldLCBpZC52YXJzID0gInNob3RfbnVtIikKY29sbmFtZXMobWVsdGVkX2RhdGEpW2NvbG5hbWVzKG1lbHRlZF9kYXRhKT09InNob3RfbnVtIl0gPC0gInNob3RfbnVtIgpjb2xuYW1lcyhtZWx0ZWRfZGF0YSlbY29sbmFtZXMobWVsdGVkX2RhdGEpPT0idmFyaWFibGUiXSA8LSAic2FtcGxlX251bSIKY29sbmFtZXMobWVsdGVkX2RhdGEpW2NvbG5hbWVzKG1lbHRlZF9kYXRhKT09InZhbHVlIl0gPC0gInNhbXBsZV92YWx1ZSIKbWVsdGVkX2RhdGEkc2FtcGxlX251bSA8LSBhcy5udW1lcmljKG1lbHRlZF9kYXRhJHNhbXBsZV9udW0pCmBgYAoKVHJhcyByZW9yZGVuYXIgbG9zIGRhdG9zIG5vcyBxdWVkYW1vcyBjb24gYG1lbHRlZF9kYXRhYCBxdWUgY29udGllbmUgdHJlcyB2YXJpYWJsZXM6IGVsIG7Dum1lcm8gZGUgZGlzcGFybyBgc2hvdF9udW1gLCBlbCBuw7ptZXJvIGRlIG11ZXN0cmEgZGUgY2FkYSBkaXNwYXJvIGBzYW1wbGVfbnVtYCB5IGVsIHZhbG9yIGRlIGxhIG11ZXN0cmEgYHNhbXBsZV92YWx1ZWAuCgpBaG9yYSBjYWxjdWxhbW9zIGxvcyB2YWxvcmVzIG3CtMOheGltb3MgcGFyYSBjYWRhIGRpc3Bhcm8geSBoYWNlbW9zIHVuIGdyw6FmaWNvIGRlIGNhamFzIHBhcmEgdmVyIHNpIHBvZGVtb3MgZWxpbWluYXIgY29zYXMgcXVlIG5vIHRlbmdhbiBwaW50YSBkZSB0ZW5lciB1biBwdWxzby4KCmBgYHtyfQptYXhzX3Blcl9zaG90IDwtIGFnZ3JlZ2F0ZShtZWx0ZWRfZGF0YSRzYW1wbGVfdmFsdWUsIGJ5ID0gbGlzdChtZWx0ZWRfZGF0YSRzaG90X251bSksIG1heCkKY29sbmFtZXMobWF4c19wZXJfc2hvdClbMV0gPC0gInNob3RfbnVtIgpjb2xuYW1lcyhtYXhzX3Blcl9zaG90KVsyXSA8LSAibWF4X3ZhbHVlIgptYXggPC0gZ2dwbG90KG1heHNfcGVyX3Nob3QsIGFlcyh4PSJtYXgiLCB5PW1heF92YWx1ZSkpICsgCiAgZ2VvbV9ib3hwbG90KCkKbWF4CmBgYAoKUHJvYmFtb3MgYSBxdWVkYXJub3MgY29uIGxvcyBkaXNwYXJvcyBjdXlvIHZhbG9yIG3DoXhpbW8gZXN0w6kgZGVudHJvIGRlbCBzZWd1bmRvIHkgdGVyY2VyIGN1YXJ0aWwuCgpgYGB7cn0KbGlicmFyeShzdGF0cykKCmZpbHRlcl9saW1pdHMgPC0gcXVhbnRpbGUobWF4c19wZXJfc2hvdCRtYXhfdmFsdWUpW2MoMiwzKV0KbG93X2xpbWl0ID0gZmlsdGVyX2xpbWl0c1sxXSAqIDEuMgpoaWdoX2xpbWl0ID0gZmlsdGVyX2xpbWl0c1syXSAqIDEuMQpmaWx0ZXJlZF9zaG90X251bWJlcnMgPC0gbWF4c19wZXJfc2hvdFsKICBtYXhzX3Blcl9zaG90JG1heF92YWx1ZSA+IGxvd19saW1pdCAgJiAKICBtYXhzX3Blcl9zaG90JG1heF92YWx1ZSA8IGhpZ2hfbGltaXQsICJzaG90X251bSJdCmBgYAoKQSBtYW5vIGhlIGHDsWFkaWRvIHVub3MgY29lZmljaWVudGVzIHBhcmEgbW9kaWZpY2FyIGVsIG3DoXhpbW8geSBlbCBtw61uaW1vLCBkZSBmb3JtYSBxdWUgc2UgcXVlZGFyYW4gZnVlcmEgY29zYXMgcXVlIGFwYXJlbnRlbWVudGUgcGFyZWNlbiBydWlkbywgYXVucXVlIGxvIHN1eW8gc2Vyw61hIGhhY2VyIGRpc3RpbnRvcyBncnVwb3MgeSBhbmFsaXphcmxvcyBwb3Igc2VwYXJhZG8gbm8gc2VhIHF1ZSBzZWFuIGNvc2FzIGludGVyZXNhbnRlcy4gTG8gbWlzbW8gY29uIGxvcyBwdWxzb3MgY3V5b3MgbcOheGltb3MgZXN0w6FuIGVuIGVsIGN1YXJ0byBjdWFydGlsOyBsb3MgcHJpbWVyb3MgcXVlIGhlIHZpc3RvIHNvbiBydWlkbywgcGVybyBsbyBtaXNtbyBoYXkgY29zYXMgaW50ZXJlc2FudGVzLgoKUG9yIMO6bHRpbW8gcGxvdGVhbW9zIGxhcyBjb3NhcyBlbiB1bmEgbWF0cml6IHBhcmEgdmV6IGxvcyBwdWxzb3MgYSB2aXN0YSBkZSBww6FqYXJvLiBIZSBpbmNsdWlkbyBkb3MgbMOtbmVhcyBob3Jpem9udGFsZXMgKGVuIGdyaXMpLCBjb24gbG9zIGzDrW1pdGVzIGluZmVyaW9yIHkgc3VwZXJpb3IgZGVsIGZpbHRybyBxdWUgaGUgdXNhZG8gcGFyYSBkaXNjcmltaW5hciBsb3MgZGlzcGFyb3MuIFJlcHJlc2VudGFtb3MgdHJlcyBncsOhZmljYXMgY29uIGxvcyBwcmltZXJvcyA2MCBkaXNwYXJvcyBxdWUgcGFzYXJvbiBlbCBmaWx0cm8uCgpgYGB7cn0Kc2hvdF9zdWJzZXQgPC0gbWVsdGVkX2RhdGFbbWVsdGVkX2RhdGEkc2hvdF9udW0gJWluJSBmaWx0ZXJlZF9zaG90X251bWJlcnNbMToyMF0sXQoKZ2dwbG90KHNob3Rfc3Vic2V0LCBhZXMoeD1zYW1wbGVfbnVtLCB5PXNhbXBsZV92YWx1ZSkpICsKICBnZW9tX3BvaW50KGNvbG91cj0iYmx1ZSIsIHNpemU9MCkgKwogIGdlb21fbGluZShjb2xvdXI9ImJsdWUiLCBhbHBoYT0wLjIpICsKICB4bGFiKCJTYW1wbGUgbnVtYmVyIikgKwogIHlsYWIoIlNhbXBsZSB2YWx1ZSIpICsKICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQ9bG93X2xpbWl0LCBjb2xvcj0iZ3JheSIpICsKICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQ9aGlnaF9saW1pdCwgY29sb3I9ImdyYXkiKSArCiAgZmFjZXRfd3JhcCh+c2hvdF9udW0pCmBgYAoKYGBge3J9CnNob3Rfc3Vic2V0IDwtIG1lbHRlZF9kYXRhW21lbHRlZF9kYXRhJHNob3RfbnVtICVpbiUgZmlsdGVyZWRfc2hvdF9udW1iZXJzWzIxOjQwXSxdCgpnZ3Bsb3Qoc2hvdF9zdWJzZXQsIGFlcyh4PXNhbXBsZV9udW0sIHk9c2FtcGxlX3ZhbHVlKSkgKwogIGdlb21fcG9pbnQoY29sb3VyPSJibHVlIiwgc2l6ZT0wKSArCiAgZ2VvbV9saW5lKGNvbG91cj0iYmx1ZSIsIGFscGhhPTAuMikgKwogIHhsYWIoIlNhbXBsZSBudW1iZXIiKSArCiAgeWxhYigiU2FtcGxlIHZhbHVlIikgKwogIGdlb21faGxpbmUoeWludGVyY2VwdD1sb3dfbGltaXQsIGNvbG9yPSJncmF5IikgKwogIGdlb21faGxpbmUoeWludGVyY2VwdD1oaWdoX2xpbWl0LCBjb2xvcj0iZ3JheSIpICsKICBmYWNldF93cmFwKH5zaG90X251bSkKYGBgCgpgYGB7cn0Kc2hvdF9zdWJzZXQgPC0gbWVsdGVkX2RhdGFbbWVsdGVkX2RhdGEkc2hvdF9udW0gJWluJSBmaWx0ZXJlZF9zaG90X251bWJlcnNbNDE6NjBdLF0KCmdncGxvdChzaG90X3N1YnNldCwgYWVzKHg9c2FtcGxlX251bSwgeT1zYW1wbGVfdmFsdWUpKSArCiAgZ2VvbV9wb2ludChjb2xvdXI9ImJsdWUiLCBzaXplPTApICsKICBnZW9tX2xpbmUoY29sb3VyPSJibHVlIiwgYWxwaGE9MC4yKSArCiAgeGxhYigiU2FtcGxlIG51bWJlciIpICsKICB5bGFiKCJTYW1wbGUgdmFsdWUiKSArCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0PWxvd19saW1pdCwgY29sb3I9ImdyYXkiKSArCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0PWhpZ2hfbGltaXQsIGNvbG9yPSJncmF5IikgKwogIGZhY2V0X3dyYXAofnNob3RfbnVtKQpgYGAKCk90cmEgZm9ybWEgY3VyaW9zYSBkZSB2ZXJsbyBlcyBhcGlsYW5kbyB0b2RvcyBsb3MgYHIgbGVuZ3RoKGZpbHRlcmVkX3Nob3RfbnVtYmVycylgIGRpc3Bhcm9zIHF1ZSBwYXNhcm9uIGVsIGZpbHRybyBpbXB1ZXN0by4KCmBgYHtyfQpmaWx0ZXJlZF9tZWx0ZWRfZGF0YSA8LSBtZWx0ZWRfZGF0YVttZWx0ZWRfZGF0YSRzaG90X251bSAlaW4lIGZpbHRlcmVkX3Nob3RfbnVtYmVycyxdCmZpbHRlcmVkX21lbHRlZF9kYXRhJHNob3RfbnVtIDwtIHBhc3RlKCJTaG90OiIsZmlsdGVyZWRfbWVsdGVkX2RhdGEkc2hvdF9udW0pCmdncGxvdChmaWx0ZXJlZF9tZWx0ZWRfZGF0YSwgYWVzKHg9c2FtcGxlX251bSwgeT1zYW1wbGVfdmFsdWUpKSArCiAgZ2VvbV9wb2ludChzaXplPTAsIGFscGhhPTAuMSkgKwogIGdlb21fbGluZShhbHBoYT0wLjEpICsKICB4bGFiKCJTYW1wbGUgbnVtYmVyIikgKwogIHlsYWIoIlNhbXBsZSB2YWx1ZSIpICsKICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQ9bG93X2xpbWl0LCBjb2xvcj0iZ3JheSIpICsKICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQ9aGlnaF9saW1pdCwgY29sb3I9ImdyYXkiKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uPSJub25lIikKYGBg