Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?

share|improve this question
    
I saw this question before. Wasn't it migrated to Signal Processing? – Dr. belisarius Feb 22 at 16:07
    
    
I am more interested in high efficiency ways of designing functions and kernels to accomplish the goal though. I will have 100's of GB's of data to analyze in the near future and have 16GB right now. – William John-Pierre Duhe Feb 22 at 16:10
1  
See also related question here: mathematica.stackexchange.com/q/42027/58 – Thomas Feb 22 at 17:02
up vote 3 down vote accepted

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

share|improve this answer
    
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. – JasonB Feb 22 at 19:09
    
@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. – N.J.Evans Feb 22 at 19:11
    
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? – William John-Pierre Duhe Feb 22 at 19:31
    
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. – N.J.Evans Feb 22 at 19:34
1  
Got cha, thanks. – William John-Pierre Duhe Feb 22 at 19:35

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

share|improve this answer
1  
delta is given in Hz so you need to multiply it by 2Pi – Simon Woods Feb 22 at 17:24
    
@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 – JasonB Feb 22 at 18:19
    
I'm suspicious of the output of BandpassFilter, I feel like those peaks should be changing more significantly from band to band. – N.J.Evans Feb 22 at 18:48
    
@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 – JasonB Feb 22 at 19:04
1  
I just noticed the fold and what you guys are talking about and am editing my response. – William John-Pierre Duhe Feb 22 at 19:27
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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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