Realtime filtering EEG signal (Bandpass Butterworth filter)

[CodeWithMichal]

Hi,

My issue is related to already existing ones: 30 27 I was trying to use approach suggested there but with no success

I want to apply Butterworth bandpass filter to realtime EEG signal. It consists of several wave bands, for example: -Theta (4-8Hz) -Alpha (8-12Hz) And others ranging from 2 to 50 Hz. The raw EEG signal ranges from -100microVolts to 100microVolts. My sampling rate is 250hz

First I tried to just create an filter and process incoming signal sample by sample:

            var filterDef = new BandPassFilter(4f/250f, 8f/ 250f, 6);

but the Process method always returns NaN. Then I switched to approach suggested in 27

            var filterDef = new BandPassFilter(4f/250f, 8f/ 250f, 6);
            var sosFilter = DesignFilter.TfToSos(filterDef.Tf);
            var chain = new FilterChain(sosFilter);

and using that I am getting results but while it should be something in range [-50,50] microVolts I am getting the values of billions.

I know it is possible to filter these bands because I am using already existing software which does it. Unfortunately the only information I can take from this software is filter definition, nothing more is provided. Below you can find images showing filter definition and realtime fitlering, where the output is delayed not more than 400ms

I know I am not providing much information but maybe you have some idea what else could be done to accomplish it? The work you do with this library is amazing - I havent found anything comparable, which could be easily used in C#

Best regards, Michal

[ar1st0crat]

Hi!

Thanks for the compliment ) It took for me 3 issues on the github to finally look carefully in my code and realize that actually we have here classes for filtering with double precision :-D. While coding filtering operations I was thinking more about audio processing, so default versions of functions operate on data with single precision, which is usually sufficient in "audio cases". But actually NWaves was intended to be a universal DSP lib. And it is (hopefully)). Sometimes I just forget about all useful features in this lib:

var tf = new Filters.Butterworth.BandPassFilter(4f/250, 8f/ 250, 5).Tf;

// it's ok: TransferFunction has always been with double precision, so use it here:
var filter = new Filters.Base64.IirFilter64(tf);

// now the filter carries out all its operations with double precision:
var filtered = _signal.Samples.Select(s => filter.Process(s));

// filter.Process() accepts one sample of type 'double' and returns 'double'
// so you have Enumerable<double> - you can LINQ-materialize it whenever you want

// maybe you have a loop, something like this:

// while (double_sample)
//     filteredSample = filter.Process(double_sample)
//     you can downcast filteredSample to float if you need

The bad news is that I didn't add FilterChain implementation for doubles (I'm going to add it in the next version). Starting at some order of the filter, it'll be more efficient to switch to SOS. It's easy to derive the correspoding class (simply replace floats with doubles and ...Filter with ...Filter64 in FilterChain class):

    public class FilterChain64 : IFilter64, IOnlineFilter64
    {
        /// <summary>
        /// List of filters in the chain
        /// </summary>
        private readonly List<IOnlineFilter64> _filters;

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="filters"></param>
        public FilterChain(IEnumerable<IOnlineFilter64> filters = null)
        {
            _filters = filters?.ToList() ?? new List<IOnlineFilter64>();
        }

        /// <summary>
        /// Constructor from collection of transfer functions (e.g., SOS sections).
        /// This constructor will create IIR (!) filters.
        /// </summary>
        /// <param name="tfs"></param>
        public FilterChain(IEnumerable<TransferFunction> tfs)
        {
            _filters = new List<IOnlineFilter64>();

            foreach (var tf in tfs)
            {
                _filters.Add(new IirFilter64(tf));
            }
        }

        /// <summary>
        /// Add filter to the chain
        /// </summary>
        /// <param name="filter"></param>
        public void Add(IOnlineFilter64 filter) => _filters.Add(filter);

        /// <summary>
        /// Insert filter at specified index into the chain
        /// </summary>
        /// <param name="idx"></param>
        /// <param name="filter"></param>
        public void Insert(int idx, IOnlineFilter64 filter) => _filters.Insert(idx, filter);

        /// <summary>
        /// Remove filter at specified index from the chain
        /// </summary>
        /// <param name="idx"></param>
        public void RemoveAt(int idx) => _filters.RemoveAt(idx);

        /// <summary>
        /// Process sample by the chain of filters
        /// </summary>
        /// <param name="input"></param>
        /// <returns></returns>
        public double Process(double input)
        {
            var sample = input;

            foreach (var filter in _filters)
            {
                sample = filter.Process(sample);
            }

            return sample;
        }

        /// <summary>
        /// Reset state of all filters
        /// </summary>
        public void Reset()
        {
            foreach (var filter in _filters)
            {
                filter.Reset();
            }
        }

        /// <summary>
        /// Offline filtering
        /// </summary>
        /// <param name="signal"></param>
        /// <param name="method"></param>
        /// <returns></returns>
        public double[] ApplyTo(double[] signal, FilteringMethod method = FilteringMethod.Auto)
        {
            return signal.Select(s => Process(s)).ToArray();
        }
    }

Regards, Tim

[ar1st0crat]

You can also try lower filter orders (like 3 and 4). Maybe the frequency selective properties of lower-order filters will be OK for your task. Order=5 means that there will be 10 coeffs in TF numerator and 10 coeffs in TF denominator.

[CodeWithMichal]

Hi Tim,

Its good to hear that! Thank you for explaining in details. I definitely need to deep dive into theory and stop using it like a black box. Could you leave this issue open for next 2-3 days? I will come back there with results and let you know if it worked for me.

Best regards, Michał

[CodeWithMichal]

Hi again Tim,

Indeed, using order 3-4 seems to solve my problem. However anything above order 8 gives again overranged data. I tried to use SOS and FilterChain64 but even that didnt help.

Could you shortly explain me what are the consequences of switching to SOS?

btw the data I am getting from device are not scaled initially to microVolts and values are quite big (range of billions). Is there any difference whether I filter original values or scaled values?

[CodeWithMichal]

Hi Tim, I am closing it for now since your suggestion helped. Orders 3-6 work fine, what is the most important for me :)

[ar1st0crat]

Hi!

Glad it helped ) Anyway, I'd like to answer your previous questions )))

I figured out what you can do with higher order filters. Take a look at this notebook with explanations regarding SOS. And in the last cell you'll find necessary C# code.

https://colab.research.google.com/drive/1A_y7sTt3qJoQMyhSv-tOT_pv6-pk4a8d?usp=sharing

Regards, Tim