Using BandpassFilter

Question

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?

Answer 1

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"
 ]

Answer 2

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}

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}]

Answer 3

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.