The problem with using a normal band-pass filter for EEG analysis is that the output is an oscillating signal, and you'd have to do peak-detection to find the amplitude. Given the short duration of transient features in the EEG stream, it would be much better to have a continuous (DC) signal as output from the filter.

For example, we have this wavelet in the input data that we want to detect, corresponding to about 5 wavelengths at 12Hz in 256Hz sampled data:

Passing it through a band-pass filter (centered on 12Hz with 4Hz width at -3dB), we get the following signal as output, which is not adequate for our purposes:

The input signal is multipled by the output of a complex oscillator
(**e ^{iwt}**) and the result of the multiplication is
passed through a low-pass filter before the modulus (absolute value)
of the complex number is taken. The oscillator runs at the centre
frequency of the band to be detected (

Applying this technique to the example wavelet above, with the same filter width and centre frequency, we now get this output:

A nice smooth output like this, which reflects the amplitude envelope of the input wavelet rather than its oscillations, is much easier to use in further processing.

Take an example input sine wave:

2*A*cos(a*t)

which can also be written:

A *
(e^{iat}+e^{-iat})

Apply the complex oscillator:

A * e^{iwt} *
(e^{iat} + e^{-iat})

giving:

A * (e^{i(w+a)t} +
e^{i(w-a)t})

Apply the low-pass filter F(w):

A * (F(w+a) * e^{i(w+a)t} + F(w-a) *
e^{i(w-a)t})

Assuming the low-pass filter is efficient and appropriate to the
frequencies being tested [*], then **F(w+a)**
will be close to 0, so we can drop that term, giving:

A * F(w-a) * e^{i(w-a)t}

Taking the modulus gives:

A * F(w-a)

So, we have successfully extracted the amplitude of the sine wave
(**A**) as a continuous value with no oscillations to worry
about.

This assumes perfect unchanging infinite-length waves, as do most
filter analysis proofs!! In real life **A** is a function of
time **A(t)**, and this amplitude envelope interacts with the
envelope of the low-pass filter by convolution. So as a rule of
thumb, the length of the measured amplitude envelope will be roughly
the length of the input amplitude envelope plus the length of the
filter envelope (i.e. the length of the filter impulse response).

[* Note: for this technique to work properly
the low-pass filter must give **F(w+a)** close to zero for all
positive '**a**'. So for a centre frequency of 15Hz, the low-pass
filter must give near zero response at 15Hz and above. If
the **F(w+a)** term isn't near 0, then that term remains and this
introduces an oscillating interference to the output magnitude. A
practical filter's response never quite goes to zero, so you need to
bear in mind that there will be a small oscillating error in your
output magnitude corresponding to whatever small amount of signal is
leaked through your filter. Check the response of your filter at and
above the centre frequency to see whether it is acceptable for your
application, and if it isn't, then narrow your filter width or select
a higher order filter (which also has its trade-offs).]

The stages can be executed exactly as described above. For example, here is some example code adapted from EEGMIR, which uses IIR filters generated with fidlib:

// 'input' is current input sample // 'lp0' and 'lp1' are working buffers for low-pass filters // 'osc' is a working buffer for generating the complex oscillator // 'output' is output of filter double val = input; sincos_step(osc); double re = lp_func(lp0, val * osc[0]); double im = lp_func(lp1, val * osc[1]); double output = hypot(re, im);

The `sincos_step()` function steps the complex oscillator
forward one angular increment corresponding to the sample period. Of
course you could just evaluate **sin(w*t)** and **cos(w*t)**
functions but it is more efficient to do it as a complex
multiplication with a constant value (**e ^{ia}**), as is
done in EEGMIR. The low-pass filter

Compared to a simple IIR band-pass filter, twice the calculations
are required (two filters) plus the oscillator and `hypot()`
calculations (two multiplies and a square root).

The low-pass filter should be selected bearing in mind the normal trade-off of sharpness in frequency response against sharpness in time response, i.e. a very sharp frequency cut-off gives very poor time-response, so transient features would go undetected. For best recognition of transient features, simple Bessel IIR filters give a nice quick and localised time-response, as well as a good locality of frequency response, although the response is domed rather than square.

Rather than executing the stages as described above, it is possible to pre-multiply the complex oscillator into the FIR filter. This means that the phase is fixed instead of changing, but as we are extracting the modulus, we don't care about the phase (and anyway we could compensate for the fixed phase if that was necessary). So we take the FIR, and pre-multiply it by the complex oscillator output resulting in two FIR filters, one giving the real part of the result, and the other giving the imaginary part. The modulus of this complex output value can be taken as before to give the result.

As with IIR, the low-pass FIR filter should be selected to give good transient detection. Window functions are well suited to this. (Window functions are all low-pass filters). Using a Blackman window, for example, the code to set up the FIR filters could be as follows (to make the example simple, parameters are hardcoded):

// Length of filter, which determines impulse response length // and filter width (in frequency). #define LEN 109 // Number of waves within filter, which determines center // frequency. #define WAVES 5 double filt[LEN], filt_re[LEN], filt_im[LEN]; void setup() { // Basic Blackman window; see fidlib source for more flexible // parameterised Blackman creation code int wid = (LEN+1)/2; memset(filt, 0, sizeof(filt)); for (int a = 0; a<wid; a++) { filt[wid-a] = filt[wid+a] = 0.42 + 0.5 * cos(M_PI * a / wid) + 0.08 * cos(M_PI * 2.0 * a / wid); } // Multiply with complex oscillator for (int a = 0; a<LEN; a++) { double theta = 2 * M_PI * WAVES * a / LEN; filt_re[a] = filt[a] * sin(theta); filt_im[a] = filt[a] * cos(theta); } }

The code to apply it might look like this:

double calculate(double *input) { double tot_re = 0, tot_im = 0; for (int a = 0; a<LEN; a++) { tot_re += filt_re[a] * input[a]; tot_im += filt_im[a] * input[a]; } return hypot(tot_re, tot_im); }

Compared to using an FIR band-pass filter, twice the calculations
are required (two FIR filters), plus the square root in
the `hypot()`. However, there is an additional optimisation
possible when using FIR filters with this technique. Because the
output of the filter calculation is continuous (DC) rather than
fluctuating (AC), it is not necessary to calculate every output sample
in order to do peak-detection. So it is possible to sample the output
at a lower frequency, and still be sure to detect all the features.
For example for EEG data at 256Hz, perhaps the FIR filters could be
applied once every 8 or 16 samples, resulting in output at 32Hz or
16Hz. This means that using this technique, the filters require much
less CPU power than a band-pass FIR.

For fast integer FIR calculations (e.g. in a DSP chip), the only
tricky part of this method is calculating the square root in integers.
There are some code samples
here.
You need to allow 16 iterations to calculate the square root of a
32-bit number. However, if a squared signal value is acceptable for
later calculations, then the square root can be skipped, e.g. using
**re*re+im*im** directly instead of calculating
**(re ^{2}+im^{2})^{0.5}**.