50Hz notch filter coefficients

edited August 2016 in OpenBCI_GUI
Dear all,

I tried to change notch filter coefficient for 50Hz noise. I got coefficient from Matlab [b a] = butter(2, [49 51]./(250/2), 'stop');
It's seem that same method as developer who developed openBCIGUI.

But, I still see very high magnitude of 50Hz . This make signal always be railed in visualizer.

Any suggestion please ....

image

Comments

  • Hi,

    I did the initial implementation of the filtering in the Processing GUI.  It was a long time ago.  It's a real pile of spaghetti, so you did a great job finding it.

    On first glance, it looks like you did it correctly.  I'll have to dig in a little deeper to see if I can find the problem.

    Chip
  • Hi chip,

    It might be other problems on my setup. You can see in other dicussion topics which were posted by me.

    Help me please. Thank you in advance! I hope to see EEG data soon! 

  • edited November 2014
    I recreated the filters in python:

    import matplotlib.pyplot as plt

    import numpy as np

    from scipy import signal


    # assumed sample rate of OpenBCI

    fs_Hz = 250.0


    # create the 60 Hz filter

    bp_stop_Hz = np.array([59.0, 61.0])

    b, a = signal.butter(2,bp_stop_Hz/(fs_Hz / 2.0), 'bandstop')

    # create the 50 Hz filter

    bp2_stop_Hz = np.array([49, 51.0])

    b2, a2 = signal.butter(2,bp2_stop_Hz/(fs_Hz / 2.0), 'bandstop')


    # compute the frequency response

    w, h = signal.freqz(b,a,1000)

    w, h2 = signal.freqz(b2,a2,1000)

    f = w * fs_Hz / (2*np.pi)             # convert from rad/sample to Hz


    And I got the following coefficients, which is exactly what you got:


    # 60 Hz Notch for fs = 250 Hz

    b = np.array([ 0.96508099, -0.24246832,  1.94539149, -0.24246832,  0.96508099])

    a = np.array([ 1.        , -0.24677826,  1.94417178, -0.23815838,  0.93138168])

    # 50 Hz Notch

    b2 = np.array([0.96508099, -1.19328255,  2.29902305, -1.19328255,  0.96508099])

    a2 = np.array([1.        , -1.21449348,  2.29780334, -1.17207163,  0.93138168])


    And I got the following coefficients, which is exactly what you got:

    image

    So far, so good.  Now to try it in the GUI itself.

    Chip



  • Thank you for proof that Chip!
  • I've got the Processing GUI revised to include a user-switchable notch filter...as soon as the software revision is discussed by the software team and accepted, you'll be able to switch between 60 Hz, 50 Hz, and none.

    I'm sorry that we North American authors were so provincial in our software design.  Hopefully, this'll make it easier for everyone else.  :)

    Chip
  • Hi Chip,

    I modified to be band pass 1 - 30 Hz. I got NaN as attached picture ?

    Any ideas?

    //butter(2,[1 30]/(250/2));  %bandpass filter
            b = new double[] { 
              8.6359264874329e-002, 0.0f, -1.72718529748659e-001, 0.0f, 8.6359264874329e-002
            };
            a = new double[] { 
              1.0f, -2.984389310913370e+000, 3.341365798136269e+000, -1.716417097831789e-001, 3.596733720384030e-001
            };
            filt_txt = "Bandpass 1-30Hz";
            short_txt = "1-30 Hz";
            break;

    image
  • P300,

    Can you post the entire function?  Or, better yet, do you have a Git where you keep your code?  Then you can point to your git and we can offer edits directly!

    The reason I need to see your code is that the snippet that you provided looks pretty good.  Also, in your screenshot, you are showing a NaN, but it doesn't look like you've actually selected your new filter.  It looks like you're still running the 1-50 Hz default filter.  I can tell because the button at the top right says "BP Filt, 1-50 Hz" and the title to the time-domain montage plot still says "1-50 Hz".

    So, either you have a bug that is preventing you from accessing your new filter, or you have introduced a bug that is affecting all of the filters (including the built-in ones).  Either way, we'd need to see more code.

    (The specific items I'll be looking for are: (1) did you add a new filter type to the existing list, or replace one of the existing ones?  (2) if you added a filter, did you make the necessary changes in the multiple different areas of the code...or did you miss one?)

    Chip
  • Oh, and the updated GUI with the 50 Hz notch filter has bee accepted and is now part of the master on the OpenBCI Github...


    Chip
  • edited November 2014

  • edited November 2014
    Sorry for confusing Chip. I really used 1-30Hz. I did not modified  title in that time I captured images. I modified as same as I did in notch 50Hz.

        //loop over all of the pre-defined filter types
        for (int Ifilt=0;Ifilt<n_filt;Ifilt++) {

          //define common notch filter
    //      b2 = new double[] { 
    //        9.650809863447347e-001, -2.424683201757643e-001, 1.945391494128786e+000, -2.424683201757643e-001, 9.650809863447347e-001
    //      };
    //      a2 = new double[] {    
    //        1.000000000000000e+000, -2.467782611297853e-001, 1.944171784691352e+000, -2.381583792217435e-001, 9.313816821269039e-001
    //      }; 
    //50Hz
          b2 = new double[] { 
            9.65080986344734e-001, -1.193282554333348e+000, 2.299023051351235e+000, -1.193282554333348e+000, 9.65080986344735e-001
            };
          a2 = new double[] {    
            1.000000000000000e+000, -1.214493479318978e+000, 2.297803341913800e+000, -1.172071629347717e+000, 9.31381682126903e-001}; 
          filtCoeff_notch[Ifilt] =  new FilterConstants(b2, a2, "Notch 60Hz", "60Hz");

          //define bandpass filter 
          switch (Ifilt) {
          case 0:
            //butter(2,[1 50]/(250/2));  %bandpass filter
    //        b = new double[] { 
    //          2.001387256580675e-001, 0.0f, -4.002774513161350e-001, 0.0f, 2.001387256580675e-001
    //        };
    //        a = new double[] { 
    //          1.0f, -2.355934631131582e+000, 1.941257088655214e+000, -7.847063755334187e-001, 1.999076052968340e-001
    //        };
    //butter(2,[1 30]/(250/2));  %bandpass filter
            b = new double[] { 
              8.6359264874329e-002, 0.0f, -1.72718529748659e-001, 0.0f, 8.6359264874329e-002
            };
            a = new double[] { 
              1.0f, -2.984389310913370e+000, 3.341365798136269e+000, -1.716417097831789e-001, 3.596733720384030e-001
            };
            filt_txt = "Bandpass 1-50Hz";
            short_txt = "1-50 Hz";
            break;
          case 1:
            //butter(2,[7 13]/(250/2));
            b = new double[] {  
              5.129268366104263e-003, 0.0f, -1.025853673220853e-002, 0.0f, 5.129268366104263e-003
            };
            a = new double[] { 
              1.0f, -3.678895469764040e+000, 5.179700413522124e+000, -3.305801890016702e+000, 8.079495914209149e-001
            };
            filt_txt = "Bandpass 7-13Hz";
            short_txt = "7-13 Hz";
            break;      
          case 2:
            //[b,a]=butter(2,[15 50]/(250/2)); %matlab command
  • edited November 2014
    For your comments.  Next time, I will post my code on Github .

    Thank you Chip
  • @chipaudette, I was looking for filter settings and configuration, which brought me here. As you are responsible for the filters I got a question for you: Is there a way to use the GUI to record filtered data rather than just the raw EMG? Also is it possible to record a new data set from playing back a [raw] recording that reflects the filters being applied in the GUI when you're looking at the waveforms you see in the time series (for example)?
    Thanks in advance.

  • wjcroftwjcroft Mount Shasta, CA

    @BCIStudent, hi.

    No, currently the GUI only records the raw data. This is so it can replay that recording using other filter settings. Each time EEG data is filtered (again), the fidelity of the signal is reduced. So having the GUI record filtered data is not a good idea. If you wish to filter EEG data, there are several libraries available, such as Brainflow, EEGLAB, Matlab, OpenViBE, etc.

    William

Sign In or Register to comment.