1
$\begingroup$

Data: https://drive.google.com/file/d/0B7G3oUuZG1TwYWhMVUhSQzNwTEk/view?usp=sharing

Here I have 20 datasets collected from custom built EMG sensors monitoring a users face. I am particularly just interested in the 3rd and 6th rows of the data ( These are the sensors monitoring for Blinks ). I am interested in automating the process of detecting where the blink features occur, clipping the data in a window around the blink, and storing that for analysis.

Here is the full data set:

ListLinePlot[blinkData, PlotRange -> Full]

graph1

blinkData = 
Flatten[Part[Import["FilePath"], {3,6}]];
Manipulate[
blinkSmoothDistribution1D = 
SmoothKernelDistribution[blinkData[[n ;; m]]];
blinkEntropy = 
NIntegrate[
With[{f = PDF[blinkSmoothDistribution1D, x]}, 
If[f > 0, f Log[f], 0]], {x, -\[Infinity], \[Infinity]}, 
MaxRecursion -> 20, WorkingPrecision -> 12];
ListLinePlot[blinkData[[n ;; m]], PlotRange -> Full, 
PlotLegends -> {"Entropy: " <> ToString[blinkEntropy]}]
, {n, 1, m - 1, 1}, {m, 100, Length[blinkData], 1}]

graph2

Here I have manually found the window that the feature exists in and calculated the differential entropy for that window, as this is the feature of the data I am most interested in. I am interested to know if there is a way to write a function to find this window where the differential entropy is actually the lowest - that corresponds to when the user was blinking.

Mvariance[x_, p_] := (
Maouty = {};
Stopv = Length[x] - p + 1;
If[p > Length[x], 
Print["moving window size must NOT be greater than ", Length[x]], 
Print["mvariance is "]];
Do[AppendTo[Maouty, Variance[x[[i ;; i + p - 1]]]], {i, Stopv}];
Return[Maouty]  )

variance = Mvariance[blinkData, 1000]

graph3

Per Bill's suggestion I have calculated the moving variance of a the data with a window size of 1000 and I can clearly see where the blink emerges. Now how do I extract the range of this spike and correspondingly clip the other data set appropriately?

$\endgroup$
  • 1
    $\begingroup$ Is what you want to detect that big spike of stuff at around samples 12000-13000? If so, how about using a moving variance calculation? $\endgroup$ – bill s Feb 29 '16 at 18:49
  • $\begingroup$ Shown above I calculated the moving variance, but how do I automate this to extract the range of this spike from the corresponding data set? $\endgroup$ – William John-Pierre Duhe Feb 29 '16 at 19:04
  • $\begingroup$ Here's one starting point for work: Module[{var, first, last}, var = MovingMap[Variance, Differences@#, 1000]; first = First@FirstPosition[var, x_ /; x > 0.000025, Length@# + 1]; last = first - 1 + First@FirstPosition[var[[first ;;]], x_ /; x <= 0.000025, Length@var - first + 1]; ListLinePlot[#[[first + 250 ;; last + 750]], PlotRange -> All]] &@ Import["Blink-3 #20.csv"][[3]] $\endgroup$ – kirma Feb 29 '16 at 19:39
  • $\begingroup$ Just to clarify, above code attempts to extract the first episode from a recording. This works by computing moving map of variances of sample-to-sample differences, and checking where variance threshold is exceeded, and where it drops below it again. $\endgroup$ – kirma Feb 29 '16 at 19:46
  • 1
    $\begingroup$ Another tool for processing variance data is the recently added FindPeaks function. $\endgroup$ – kirma Feb 29 '16 at 19:52
3
$\begingroup$

I used to work on this long time ago, we used matlab back in the day.

This is how I remember we did it.

First, import the data

blinkData = 
  Table[Import[
    "Downloads/Blink-3/Blink-3 #" <> ToString[i] <> ".csv"], {i, 2, 
    20}];

Then, you listen to the difference of channel 6 and 3, this increased the signal to noise ratio in the few channels I looked at. I use a lowpass filter to get rid of the noise. (This could be improved since you know the frequencies you are looking at, so you don't care about gamma or beta waves). On that filtered data you look for the maximum, that will be your centre for the ROI.

Show[ListLinePlot[
  LowpassFilter[#[[6]] - #[[3]], .0050] & /@ blinkData, 
  PlotRange -> Full], 
 ListPlot[{First@
      First@Position[LowpassFilter[#[[6]] - #[[3]], .0050], 
        Max[LowpassFilter[#[[6]] - #[[3]], .0050]]], 
     Max[LowpassFilter[#[[6]] - #[[3]], .0050]]} & /@ blinkData, 
  PlotStyle -> {Red, Large}, Joined -> False]]

enter image description here

then your new data is a window around those maxima

newData = 
 MapIndexed[blinkData[[First[#2]]][[All, #1 - 1000 ;; # + 1000]] &, 
  First@First@
      Position[LowpassFilter[#[[6]] - #[[3]], .0050], 
       Max[LowpassFilter[#[[6]] - #[[3]], .0050]]] & /@ blinkData]

If somebody can make it look nicer, please edit. That First@First@Position looks horrible but is already late and it just worked.

EDIT:

Sorry, forgot to add the plot of the clipped data and the mean of it, which was what we were interested on at the end of the day.

enter image description here

And finally you get an average spike with the new data, whose duration we used to "clean" the raw data of the eye movements.

ListLinePlot[
 Mean[LowpassFilter[#[[6]] - #[[3]], 0.0050] & /@ newData], 
 PlotRange -> Full]

Hope this helps.

enter image description here

$\endgroup$
  • $\begingroup$ I was just about to post my own answer but yours is more elegant it seems! I will go ahead and post my methodology either way because I use a few different techniques but over all I don't think I could have asked for a better answer. It sounds like you have a lot of experience in this particular field. If you are free to Skype some time and share a little bit about your experience that would be great, you may think the projects I am working on are pretty cool too. Just shoot me a message. $\endgroup$ – William John-Pierre Duhe Mar 1 '16 at 17:47
  • $\begingroup$ well, I did most of the matlab code for a neuroscience lab for a few months, implemented algorithms and generally helped with coding problems for the guys of the lab. Then I went back to physics for my masters and phd... just drop me a line in skype or gmail, both have the same username as here. $\endgroup$ – tsuresuregusa Mar 1 '16 at 19:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.