2
$\begingroup$

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

Above I have linked a data set containing a 2 seconds EEG sample I have recorded of someone "blinking". This data set was taken at a 220hz frequency. I am attempting to apply the function BandpassFilter[] function to isolate specific frequency bands in the EEG spectrum listed in the code.

data= Flatten@Import["C:\\Users\\YottaTech\\Desktop\\FacialData\\Duhe\\
John\\Blink\\1\\Blink.data.1.26.10.9.4.2167734.csv"];

(*These are the frequency band ranges (in Hz) I would like to filter out*)
delta = {.1, 3};
theta = {4, 7};
alpha = {8, 12};
beta = {12, 30};
gamma = {30,50};

(*Here is my initial attempt to filter out the ""Delta" frequency band*)
deltaFreq = BandpassFilter[data, {2*Pi*.1,2*Pi*3}];

Would anyone be able to show me where I am going wrong with this? Or show more a more powerful/efficient function to utilize to accomplish the task?

Improve this question
$\endgroup$
3
$\begingroup$

Here's an alternative to using BandpassFilter, which may or may not be any better.

You can use Fourier to extract only a few components. 

delta = {0.1,3}; (*in Hz*)
theta = {4, 7};
alpha = {8, 12};
beta = {12, 30};
gamma = {30, 50};

Take the elements you want, and create a filtered FFT, we're effectively creating a square bandpass filter. This will not handle frequencies above sample freq./2 without more care.

fftFilter[data_, sR_, band_] := Module[
  {dur, fRes, f1, f2, rSpan, f},
  dur = Length@data/sR;(*Total duration of data.*)
  fRes = 1/dur;(*Resolution of FFT is 1/duration.*)
  (*frequencies in terms of Fourier part.1=DC*)
  {f1, f2} = 1 + Quotient[#, fRes] & /@ band;

  rSpan = If[f1 == 1, Span @@ {2, -1}, Span @@ {1, -1}];
  f = Fourier[data, List /@ Range[f1, f2]];
  f = PadLeft[f, f2];
  f = Join[PadRight[f, (Length@data) - f2 + 1], 
    Conjugate@Reverse@(f[[rSpan]])];
  f = PadRight[f, Length@data];
  Re@InverseFourier[f]
  ]

I plot each of the bands separately and then combine all of the band information as a check. It seems to look ok. 

ListPlot[
 MapIndexed[# + 400*(First@#2 - 1) &,
  # - Mean@# & /@ (MapIndexed[
     fftFilter[data, 220, #] &, {delta, theta, alpha, beta, gamma}
     ])
  ]
 , Joined -> True
 , PlotLegends -> LineLegend[
   ColorData[97][#] & /@ Range[5]
   ,
   {"delta", "theta", "alpha", "beta", "gamma"}
   ]
 , PlotLabel -> "Band Data"
 ]

ListPlot[
 {
  data,
  Total[
   fftFilter[data, 220, #] & /@ {delta, alpha, beta, gamma, theta}
   ]
  }
 , PlotLegends -> LineLegend[
   ColorData[97][#] & /@ Range[2]
   ,
   {"Data", "\[CapitalSigma] bands"}
   ]
 , Joined -> True
 , PlotRange -> All
 , PlotLabel -> "Sum of Bands"
 ]

enter image description here

Improve this answer
$\endgroup$
  • $\begingroup$ So is the difference between red and blue the 2 Pi I missed? It would make sense since the red curve has more high frequency noise cut out. $\endgroup$ – Jason B. Feb 22 '16 at 19:09
  • $\begingroup$ @JasonB, I changed the pictures. The BandpassFilter solution isn't shown in these. But no, adding the factor of 2pi doesn't seem to reproduce what I've shown, it still shows more high frequency than seems reasonable to me. $\endgroup$ – N.J.Evans Feb 22 '16 at 19:11
  • $\begingroup$ What do you mean when you say it won't handle frequencies bellow sample frequency/2? Does this mean I can't filter frequencies bellow 110Hz? $\endgroup$ – William John-Pierre Duhe Feb 22 '16 at 19:31
  • $\begingroup$ Above sample frequency/2, not below. That's just the Nyquist limit. What I meant was if you ask the function I wrote for a frequency above that it's going to freak out because I didn't put any checks in. $\endgroup$ – N.J.Evans Feb 22 '16 at 19:34
  • 1
    $\begingroup$ Got cha, thanks. $\endgroup$ – William John-Pierre Duhe Feb 22 '16 at 19:35
6
$\begingroup$

If your sampling rate is 220 Hz, then you need to tell Mathematica that, otherwise the frequencies you are giving are meaningless. 

deltaFreq = BandpassFilter[data, delta, SampleRate -> 220];

ListLinePlot /@ {data, deltaFreq}

enter image description here

Here you can see the result of successive action of the BandpassFilter

Manipulate[
 Fold[(BandpassFilter[#1, #2, SampleRate -> 220] &), 
   data, {{0, 10000}, delta, theta, alpha, beta, gamma}[[;; n]]] // 
  ListLinePlot, {{n, 1}, 1, 6, 1}]

enter image description here

Improve this answer
$\endgroup$
  • 1
    $\begingroup$ delta is given in Hz so you need to multiply it by 2Pi $\endgroup$ – Simon Woods Feb 22 '16 at 17:24
  • $\begingroup$ @SimonWoods That makes sense, I'm away from the office so I can't fix it now but I did get the biggest mistake in the OP I think with the SampleRate option $\endgroup$ – Jason B. Feb 22 '16 at 18:19
  • $\begingroup$ I'm suspicious of the output of BandpassFilter, I feel like those peaks should be changing more significantly from band to band. $\endgroup$ – N.J.Evans Feb 22 '16 at 18:48
  • $\begingroup$ @N.J.Evans It might have been more interesting to apply the filters in the opposite order, or one at a time rather than cumulatively $\endgroup$ – Jason B. Feb 22 '16 at 19:04
  • 1
    $\begingroup$ I just noticed the fold and what you guys are talking about and am editing my response. $\endgroup$ – William John-Pierre Duhe Feb 22 '16 at 19:27
1
$\begingroup$ 
data = Flatten@
Import["C:\\Users\\YottaTech\\Desktop\\FacialData\\Duhe \
William\\Blink\\1\\Blink.data.1.26.10.26.49.6077426.csv"];

labels = {"Delta", "Theta", "Alpha", "Beta", "Gamma", "High"};
delta = 2*\[Pi]*{.1, 3};
theta = 2*\[Pi]*{4, 7};
alpha = 2*\[Pi]*{8, 12};
beta = 2*\[Pi]*{12, 30};
gamma = 2*\[Pi]*{30, 50};
high = 2*\[Pi]*{50, 80};

freqBands = {delta, theta, alpha, beta, gamma, high};

Manipulate[
Column[{
ListLinePlot[data, PlotRange -> Full, PlotLabel -> "Raw"], 
ListLinePlot[
BandpassFilter[data, freqBands[[n]], SampleRate -> 220], 
PlotLabel -> labels[[n]], PlotRange -> Full]
}],
{{n, 1}, 1, Length[freqBands], 1}]

After reviewing all the help I have arrived at a correctly working bandpass filter approach! Thank you guys very much for the help and I have included the working code above for reference.

Improve this answer
$\endgroup$

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.