2
$\begingroup$

Let's say I have a 2 second data set taken at 220Hz sample rate and I would like to filter out the frequency bands associated with the EEG Spectrum: $$\begin{align} \Delta:& [1,3]\text{ Hz}\\ \theta:& [4,7]\text{ Hz}\\ \alpha_1:& [8,9]\text{ Hz}\\ \alpha_2:& [10,12]\text{ Hz}\\ \beta_1:& [13,17]\text{ Hz}\\ \beta_2:& [18,30]\text{ Hz}\\ \gamma_1:& [31,40]\text{ Hz}\\ \gamma_2:& [41,50]\text{ Hz} \end{align} $$ What would be the most simple approach to do this ?

$\endgroup$
2
$\begingroup$

What would be the most simple approach to do this ?

I'm a GNU Radio person. So the simplest approach was this GNU Radio companion-designed flowgraph:

Bandpass flow graph

Generating this flowgraph file generates a python program that uses GNU Radio to do the signal processing, on as many CPU cores as there are.

You can then run it as

$ ./bandpass_filters.py -f {input file containing float32 samples, one after the other}

To decrease the complexity of all these bandpasses, I first used a low pass to reduce the sampling rate of the signal by half – that's Ok, since you are only interested in frequencies below half of half of your input sampling rate.

I tried this with a 200MB "dummy" file, creating eight 100MB output files – this, to and from a temp directory, took a whole 13 seconds. I'm pretty hopeful it's fast enough! I benchmarked this a bit, and it seems the source is able to push through about 20 million samples per second; the slow part is the writing to the 8 files.

Note that this is a quick and easy solution – a proper, speed-optimized solution would probably use

  1. decimation, because you're output sampling rate is still 110 Hz, which contains a maximum of 9 Hz (that's a waste of processing power, and hence, speed)
  2. another pre-decimation by 2 through low-pass filtering to 27.5 Hz for the four filters below that rate

However, really, at your modest amount of samples (16GB = 4 Billion 32bit floats?) this isn't really necessary, if you wanted to grab a coffee, anyways.

That flowgraph was really nice and easy to design. Here I am, shamelessly plugging GNU Radio's Guided tutorials if you want to learn to design such signal processing flow graphs yourself; in fact, it's even more fun to use GNU Radio in live flowgraphs, i.e. in systems where the in- or output (or both) are physical devices (microphones, ecg sensors, sonars, radio frontends). Exactly the same flow graph can work with a sound card, if you replace the "file source" by an "audio source" (both come with GNU Radio) and specify a sound card-compatible sampling rate via the -r flag.

Now, since you probably know better than me how to work with Mathematica, here's what I'd recommend you implement:

  1. (first, one time) You use the GNU Radio companion to generate a python program that does the signal processing for you.
  2. You use Mathematica's BinaryWrite() to write your raw samples to a file as 32 bit floating point numbers, e.g. BinaryWrite("/tmp/samples.dat.f32", your_sample_vector, "Real32").
  3. You use Mathematica's abilities to call executables on your machine to execute that python script. Now, I'm no Mathematica expert, so this is all Google-based-guessing: RunProcess({"/path/of/python/executable", "/path/of/generated/python/program", "-f", "/tmp/samples.dat.f32"})
  4. You read (BinaryRead) the resulting /tmp/samples.dat.f32_1_3.dat (and …_4_7.dat etc) in Mathematica.

Optionally, you really consider doing less workflow in Mathematica and use the sheer mass of signal processing things that come with GNU Radio or its module ecosystem.

$\endgroup$
1
$\begingroup$

The sampling frequency seems right (220 Hz > 2*50 Hz) so there's no loss of information.

What you have to do is design a bandpass filter for each frequency band and apply it to your signal to extract each band. Since you seem to have "unlimited" computing amount, you can easily apply high-order filters and get your signals.

Since you are using mathematica, I suggest this function : https://reference.wolfram.com/language/ref/BandpassFilter.html Which allows to specify frequencies and other parameters easily without having to design the whole thing.

For example, let's say you have a signal $S$ in mathematica, which is defined by an array of numbers {12,21,...,36,...}. And you want to extract the $\alpha_2$ waves. You'll have to use the function as follow (n is the filter's kernel length)

BandPassFilter[S, {w1,w2},n]

w1 and w2 being defined as : $$ w_1 = 2*\pi*10\\ w_2 = 2*\pi*12 $$

If you are free to use any software, I also suggest you switch to MatLab, which is more appropriate than mathematica (it's my opinion) to process data as you want to do, and widely used in the signal engineering community.

$\endgroup$
  • $\begingroup$ Thanks for the information! I have looked at the BandpassFilters as well as using LowpassFilters and BandstopFilters. I was wondering if you could perhaps ellaborate a little bit on how to actually design one or each of these functions in the native Mathematica code. $\endgroup$ – William John-Pierre Duhe Feb 22 '16 at 15:47
  • $\begingroup$ I also do have several VERY powerful computation stations and access to clusters if need be. $\endgroup$ – William John-Pierre Duhe Feb 22 '16 at 15:48
  • $\begingroup$ I'll add this in my answer $\endgroup$ – MaximGi Feb 22 '16 at 15:53
  • $\begingroup$ By the way, you won't need super calculators to do this (unless you have GB of data, but for a 220 Hz sampling, I doubt it). I was talking about "unlimited" computation because you're using a desktop computer and not implementing this on a lighter system like an arduino/raspberry-pi/whatever... $\endgroup$ – MaximGi Feb 22 '16 at 16:05
  • 1
    $\begingroup$ 16GB ? Assuming these are 16bit samples each, that's a mere 8 Billion samples … meeh, will take a minute or two on a modern workstation with enough cores to dedicate one core to every filter, unless you "overspecify" your filters. Really no need to bring out the big guns, @MaximGi. Notice that I think that processing effort would be in the range of the writing speed of your storage. $\endgroup$ – Marcus Müller May 22 '16 at 22:33

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.