Is the FilterDesign's Butterworth filter wrong in BCI2000

Forum for software developers to discuss BCI2000 software development
Locked
shjlynn
Posts: 22
Joined: 27 Dec 2010, 08:28

Is the FilterDesign's Butterworth filter wrong in BCI2000

Post by shjlynn » 05 Jul 2012, 19:53

hello,

Have anyone used the Butterworth Bandpass filter in the FilterDesign.cpp??
I've tried it, but did not get the right results. I looked into it, and found that the Poles and Zeros caculated by the program is wrong. I can't figure out why.

for example, I want to bandpass the signal between 5.25-6.75Hz, my Signal SamplingRate=1000Hz, bandpass filter order=5. So I write the code like this:

TransferFunction bp =
FilterDesign::Butterworth()
.Order( 5 )
.Bandpass( 5.25/1000, 6.75/1000 )
.TransferFunction();

Is any parameter setting mistake? thank you!

Haijun

mellinger
Posts: 1052
Joined: 12 Feb 2003, 11:06

Re: Is the FilterDesign's Butterworth filter wrong in BCI200

Post by mellinger » 06 Jul 2012, 08:48

How do you know that the poles and zeros are wrong? In case you want to compare the designed filter with one designed by Matlab, you need to understand that Matlab computes filter _coefficients_, while BCI2000 filter design computes a transfer function in terms of poles and zeros.

For more information, see
http://www.bci2000.org/wiki/index.php/P ... efficients
http://www.bci2000.org/wiki/index.php/P ... nformation

Regards,
Juergen

shjlynn
Posts: 22
Joined: 27 Dec 2010, 08:28

Re: Is the FilterDesign's Butterworth filter wrong in BCI200

Post by shjlynn » 06 Jul 2012, 13:18

when I program the code like what I showed, I observe the all the poles in complex plane, they were all in right half plane. that's obviously wrong.

so is there any mistake in setting the parameters?

Haijun
mellinger wrote:How do you know that the poles and zeros are wrong? In case you want to compare the designed filter with one designed by Matlab, you need to understand that Matlab computes filter _coefficients_, while BCI2000 filter design computes a transfer function in terms of poles and zeros.

For more information, see
http://www.bci2000.org/wiki/index.php/P ... efficients
http://www.bci2000.org/wiki/index.php/P ... nformation

Regards,
Juergen

shjlynn
Posts: 22
Joined: 27 Dec 2010, 08:28

Re: Is the FilterDesign's Butterworth filter wrong in BCI200

Post by shjlynn » 06 Jul 2012, 17:58

sorry, I've figured out the problem. My understanding make a mistake, the code itself is right.

But still have a quick question, when it comes to butterworth bandpass, what's the value that be set in " bp.Evaluate( ??? ) "

typedef Ratpoly<FilterDesign::Complex> TransferFunction;
TransferFunction tf( 1.0 );
Real gain = 1.0;
TransferFunction bp =
FilterDesign::Butterworth()
.Order( mBandpassOrder )
.Bandpass( LowEdge, HighEdge )
.TransferFunction();
tf *= bp;

gain /= abs( bp.Evaluate( ??? ) );
shjlynn wrote:when I program the code like what I showed, I observe the all the poles in complex plane, they were all in right half plane. that's obviously wrong.

so is there any mistake in setting the parameters?

Haijun
mellinger wrote:How do you know that the poles and zeros are wrong? In case you want to compare the designed filter with one designed by Matlab, you need to understand that Matlab computes filter _coefficients_, while BCI2000 filter design computes a transfer function in terms of poles and zeros.

For more information, see
http://www.bci2000.org/wiki/index.php/P ... efficients
http://www.bci2000.org/wiki/index.php/P ... nformation

Regards,
Juergen

mellinger
Posts: 1052
Joined: 12 Feb 2003, 11:06

Re: Is the FilterDesign's Butterworth filter wrong in BCI200

Post by mellinger » 09 Jul 2012, 07:33

For a bandpass, you will want to normalize the gain in the center of the passband, so you evaluate at a position corresponding to that normalized frequency f_c:
z = exp( +/- i * 2 * pi * f_c )

For a highpass, one evaluates at f_c = 0, i.e. z = 1.
For a lowpass, one evaluates at f_c = 0.5, i.e. z = -1.

Regards,
Juergen

shjlynn
Posts: 22
Joined: 27 Dec 2010, 08:28

Re: Is the FilterDesign's Butterworth filter wrong in BCI200

Post by shjlynn » 09 Jul 2012, 12:33

Thank you,Juergen.
So, if my passband is 5.25-6.75Hz, Signal SamplingRate=1000Hz.
Then it will be: z = exp( +/- i * 2 * pi * 0.006) ???

Haijun
mellinger wrote:For a bandpass, you will want to normalize the gain in the center of the passband, so you evaluate at a position corresponding to that normalized frequency f_c:
z = exp( +/- i * 2 * pi * f_c )

For a highpass, one evaluates at f_c = 0, i.e. z = 1.
For a lowpass, one evaluates at f_c = 0.5, i.e. z = -1.

Regards,
Juergen

mellinger
Posts: 1052
Joined: 12 Feb 2003, 11:06

Re: Is the FilterDesign's Butterworth filter wrong in BCI200

Post by mellinger » 09 Jul 2012, 13:34

Right.

griffin.milsap
Posts: 58
Joined: 08 Jun 2009, 12:42

Re: Is the FilterDesign's Butterworth filter wrong in BCI200

Post by griffin.milsap » 11 Jul 2012, 14:44

Just in case anyone else is having trouble with this, here are my thoughts. I've been playing with this FilterDesign toolbox for the last few days, and I had a rather limited grasp on this whole concept when I started. Hopefully this clarifies a few points (and I'm not just spouting incorrect information)

Lets say you have a signal sampled at 1000Hz with 10 muV of noise and a 500 muV 100Hz Sine wave. We can design a filter for the 100Hz component of this signal by creating a transfer function, but we first need to look at the Z plane.

Image

In order to make sense of this space, all of the frequencies must be normalized to the sampling rate. A normalized frequency of 0 Hz is (0/1000) = 0.0, and a normalized frequency of 1000Hz is (1000/1000) = 1.0. According to the nyquist criterion, frequencies above half of the sampling rate (500Hz) will show artifacts from aliasing, so we are only interested in frequencies up to 500Hz in this signal (500/1000) = 0.5.

Frequencies map to locations in the Z plane by using this transformation: z = exp( i * 2Pi * f ) This means that a frequency of 0Hz - 0.0 (DC) exists at z = exp( i*2*Pi*0 ) = exp( 0 + 0i ) = (1.0 + i0.0), and a frequency of 500Hz - 0.5 (nyquist) exists at z = exp( i * 2Pi * 0.5 ) = exp( i * Pi ) = (-1.0 + 0.0i). It just so happens that if you interpolate this number from 0Hz to the nyquist frequency (500Hz in this case), you'll find that this traces along the top half of the unit circle from (1.0,0.0) to (-1.0,0.0).

The transfer function created by FilterDesign is a fraction of complex polynomials. Everywhere the numerator of this fraction is equal to zero, the resulting evaluation of this transfer function is zero. As such, these locations are called zeros. Everywhere the denominator of this fraction is equal to zero, the resulting evaluation of this transfer function is infinite (divide by zero). These locations are called poles. If you were to imagine a large sheet of tent fabric on this plane, propped up by tent-poles at the pole locations, and nailed to the ground at the zero locations, you would have a nice visual representation of your transfer function in the z-plane. Where this "tent-cloth" intersects the top half of the unit circle, you will see the frequency response of the filter you've designed. You can plot it in a figure called a bode plot and see what frequencies are suppressed (stopped) and which are passed.

These filters do not necessarily maintain the amplitude of the original signal, however. If you were to band-pass filter our example signal from 90-110 Hz, you would not see 500muV of output 100Hz activity. About all we can do is attempt to normalize the gain, but even this will not accurately output the original amplitude of this signal.

In the case of a low-pass filter, we want to normalize the gain such that the gain on the lowest frequency within the pass band (DC) is unity or close to it. As such, we want to Evaluate() our transfer function at (1.0 + 0.0i), or simply 1.0. We divide our unity gain by the result of this evaluation such that when we multiply by the gain later, we get a unity gain at this frequency.

For a high-pass filter, we want to apply the same idea, but we have to move our evaluation point to somewhere in the pass band. The nyquist frequency should be in the pass band of a high pass filter, so we can evaluate it there (-1.0 + 0.0i).

A Band-stop filter is interesting, in that it should pass both low and high frequencies but stop some frequencies in between. In this case, you could likely evaluate the transfer function at either -1.0 or 1.0.

A band-pass filter, however, will stop frequencies at both -1.0 and 1.0 (depending on your design parameters), so you will want to evaluate the gain at the peak of your pass-band, which should be the center frequency. In this instance, to isolate the 100Hz component, we will want a band-pass filter with a cuton at 90Hz and a cutoff of 110Hz. This puts our center frequency, f_c = 90+((110-90)/2.0) = 100Hz = 0.1. We can then evaluate the transfer function at exp( complex( 0, 2*pi*f_c ) ) and divide by the magnitude of the complex result. It just so happens that abs() is overloaded to deal with complex numbers in just this way.

In short, here is some helpful source code for dealing with FilterDesign.

Code: Select all

    if( type == "Butterworth" || type == "Chebyshev" )
    {
      FilterDesign::Butterworth* filt = NULL;
      if( type == "Butterworth" )
        filt = new FilterDesign::Butterworth();
      else
      {
        FilterDesign::Chebyshev* chebyFilt = new FilterDesign::Chebyshev();
        chebyFilt->Ripple_dB( ripple );
        filt = chebyFilt;
      }

      filt->Order( order );
      
      float cuton_freq = ( float )cuton / ( float )input.SamplingRate();
      float cutoff_freq = ( float )cutoff / ( float )input.SamplingRate();
      FilterDesign::Complex z( 0, 0 );
      
      if( character == "Lowpass" )
      {
        filt->Lowpass( cutoff_freq );
        z = 1.0f;
      }
      else if( character == "Highpass" )
      {
        filt->Highpass( cutoff_freq );
        z = -1.0f;
      }
      else if( character == "Bandpass" )
      {
        filt->Bandpass( cuton_freq, cutoff_freq );
        double f_c = cuton_freq + ( ( cutoff_freq - cuton_freq ) / 2.0f );
        z = exp( FilterDesign::Complex( 0, 2 * M_PI * f_c ) );
      }
      else if( character == "Bandstop" )
      {
        filt->Bandstop( cuton_freq, cutoff_freq );
        z = 1.0f;
      }
      else
        bcierr << "Unknown filter characteristic \"" << character 
               << "\" on row " << row << " of " << matrix->Name() << endl;
               
      transferFunction *=  filt->TransferFunction();
      FilterDesign::Complex complexGain = filt->TransferFunction().Evaluate( z );
      gain /= abs( complexGain );
    }
You guys probably knew most of that, but I'm still trying to understand it myself and the pure mathematical descriptions didn't really make as much sense to me. I find it easier to understand what's going on intuitively by thinking about it this way.

Hope that helps,

-Griff

Locked

Who is online

Users browsing this forum: No registered users and 1 guest