title: bandpass butterworth filter implementation in C++
body: <p>I am implementing an image analysis algorithm using openCV and c++, but I found out openCV doesnt have any function for Butterworth Bandpass filter officially. 
  in my project I have to pass a time series of pixels into the Butterworth 5 order filter and the function will return the filtered time series pixels. Butterworth(pixelseries,order, frequency), if you have any idea to help me of how to start please let me know. Thank you</p>
  
  <p><strong>EDIT :</strong>
  after getting help, finally I come up with the following code. which can calculate the Numerator Coefficients and Denominator Coefficients, but the problem is that some of the numbers is not as same as matlab results. here is my code:</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;stdio.h&gt;
  #include &lt;vector&gt;
  #include &lt;math.h&gt;
  
  using namespace std;
  
  #define N 10 //The number of images which construct a time series for each pixel
  #define PI 3.14159
  
  double *ComputeLP( int FilterOrder )
  {
      double *NumCoeffs;
      int m;
      int i;
  
      NumCoeffs = (double *)calloc( FilterOrder+1, sizeof(double) );
      if( NumCoeffs == NULL ) return( NULL );
  
      NumCoeffs[0] = 1;
      NumCoeffs[1] = FilterOrder;
      m = FilterOrder/2;
      for( i=2; i &lt;= m; ++i)
      {
          NumCoeffs[i] =(double) (FilterOrder-i+1)*NumCoeffs[i-1]/i;
          NumCoeffs[FilterOrder-i]= NumCoeffs[i];
      }
      NumCoeffs[FilterOrder-1] = FilterOrder;
      NumCoeffs[FilterOrder] = 1;
  
      return NumCoeffs;
  }
  
  double *ComputeHP( int FilterOrder )
  {
      double *NumCoeffs;
      int i;
  
      NumCoeffs = ComputeLP(FilterOrder);
      if(NumCoeffs == NULL ) return( NULL );
  
      for( i = 0; i &lt;= FilterOrder; ++i)
          if( i % 2 ) NumCoeffs[i] = -NumCoeffs[i];
  
      return NumCoeffs;
  }
  
  double *TrinomialMultiply( int FilterOrder, double *b, double *c )
  {
      int i, j;
      double *RetVal;
  
      RetVal = (double *)calloc( 4 * FilterOrder, sizeof(double) );
      if( RetVal == NULL ) return( NULL );
  
      RetVal[2] = c[0];
      RetVal[3] = c[1];
      RetVal[0] = b[0];
      RetVal[1] = b[1];
  
      for( i = 1; i &lt; FilterOrder; ++i )
      {
          RetVal[2*(2*i+1)]   += c[2*i] * RetVal[2*(2*i-1)]   - c[2*i+1] * RetVal[2*(2*i-1)+1];
          RetVal[2*(2*i+1)+1] += c[2*i] * RetVal[2*(2*i-1)+1] + c[2*i+1] * RetVal[2*(2*i-1)];
  
          for( j = 2*i; j &gt; 1; --j )
          {
              RetVal[2*j]   += b[2*i] * RetVal[2*(j-1)]   - b[2*i+1] * RetVal[2*(j-1)+1] +
                  c[2*i] * RetVal[2*(j-2)]   - c[2*i+1] * RetVal[2*(j-2)+1];
              RetVal[2*j+1] += b[2*i] * RetVal[2*(j-1)+1] + b[2*i+1] * RetVal[2*(j-1)] +
                  c[2*i] * RetVal[2*(j-2)+1] + c[2*i+1] * RetVal[2*(j-2)];
          }
  
          RetVal[2] += b[2*i] * RetVal[0] - b[2*i+1] * RetVal[1] + c[2*i];
          RetVal[3] += b[2*i] * RetVal[1] + b[2*i+1] * RetVal[0] + c[2*i+1];
          RetVal[0] += b[2*i];
          RetVal[1] += b[2*i+1];
      }
  
      return RetVal;
  }
  
  double *ComputeNumCoeffs(int FilterOrder)
  {
      double *TCoeffs;
      double *NumCoeffs;
      int i;
  
      NumCoeffs = (double *)calloc( 2*FilterOrder+1, sizeof(double) );
      if( NumCoeffs == NULL ) return( NULL );
  
      TCoeffs = ComputeHP(FilterOrder);
      if( TCoeffs == NULL ) return( NULL );
  
      for( i = 0; i &lt; FilterOrder; ++i)
      {
          NumCoeffs[2*i] = TCoeffs[i];
          NumCoeffs[2*i+1] = 0.0;
      }
      NumCoeffs[2*FilterOrder] = TCoeffs[FilterOrder];
  
      free(TCoeffs);
  
      return NumCoeffs;
  }
  
  double *ComputeDenCoeffs( int FilterOrder, double Lcutoff, double Ucutoff )
  {
      int k;            // loop variables
      double theta;     // PI * (Ucutoff - Lcutoff) / 2.0
      double cp;        // cosine of phi
      double st;        // sine of theta
      double ct;        // cosine of theta
      double s2t;       // sine of 2*theta
      double c2t;       // cosine 0f 2*theta
      double *RCoeffs;     // z^-2 coefficients
      double *TCoeffs;     // z^-1 coefficients
      double *DenomCoeffs;     // dk coefficients
      double PoleAngle;      // pole angle
      double SinPoleAngle;     // sine of pole angle
      double CosPoleAngle;     // cosine of pole angle
      double a;         // workspace variables
  
      cp = cos(PI * (Ucutoff + Lcutoff) / 2.0);
      theta = PI * (Ucutoff - Lcutoff) / 2.0;
      st = sin(theta);
      ct = cos(theta);
      s2t = 2.0*st*ct;        // sine of 2*theta
      c2t = 2.0*ct*ct - 1.0;  // cosine of 2*theta
  
      RCoeffs = (double *)calloc( 2 * FilterOrder, sizeof(double) );
      TCoeffs = (double *)calloc( 2 * FilterOrder, sizeof(double) );
  
      for( k = 0; k &lt; FilterOrder; ++k )
      {
          PoleAngle = PI * (double)(2*k+1)/(double)(2*FilterOrder);
          SinPoleAngle = sin(PoleAngle);
          CosPoleAngle = cos(PoleAngle);
          a = 1.0 + s2t*SinPoleAngle;
          RCoeffs[2*k] = c2t/a;
          RCoeffs[2*k+1] = s2t*CosPoleAngle/a;
          TCoeffs[2*k] = -2.0*cp*(ct+st*SinPoleAngle)/a;
          TCoeffs[2*k+1] = -2.0*cp*st*CosPoleAngle/a;
      }
  
      DenomCoeffs = TrinomialMultiply(FilterOrder, TCoeffs, RCoeffs );
      free(TCoeffs);
      free(RCoeffs);
  
      DenomCoeffs[1] = DenomCoeffs[0];
      DenomCoeffs[0] = 1.0;
      for( k = 3; k &lt;= 2*FilterOrder; ++k )
          DenomCoeffs[k] = DenomCoeffs[2*k-2];
  
  
      return DenomCoeffs;
  }
  
  void filter(int ord, double *a, double *b, int np, double *x, double *y)
  {
      int i,j;
      y[0]=b[0] * x[0];
      for (i=1;i&lt;ord+1;i++)
      {
          y[i]=0.0;
          for (j=0;j&lt;i+1;j++)
              y[i]=y[i]+b[j]*x[i-j];
          for (j=0;j&lt;i;j++)
              y[i]=y[i]-a[j+1]*y[i-j-1];
      }
      for (i=ord+1;i&lt;np+1;i++)
      {
          y[i]=0.0;
          for (j=0;j&lt;ord+1;j++)
              y[i]=y[i]+b[j]*x[i-j];
          for (j=0;j&lt;ord;j++)
              y[i]=y[i]-a[j+1]*y[i-j-1];
      }
  }
  
  
  
  
  int main(int argc, char *argv[])
  {
      //Frequency bands is a vector of values - Lower Frequency Band and Higher Frequency Band
  
      //First value is lower cutoff and second value is higher cutoff
      double FrequencyBands[2] = {0.25,0.375};//these values are as a ratio of f/fs, where fs is sampling rate, and f is cutoff frequency
      //and therefore should lie in the range [0 1]
      //Filter Order
  
      int FiltOrd = 5;
  
      //Pixel Time Series
      /*int PixelTimeSeries[N];
      int outputSeries[N];
      */
      //Create the variables for the numerator and denominator coefficients
      double *DenC = 0;
      double *NumC = 0;
      //Pass Numerator Coefficients and Denominator Coefficients arrays into function, will return the same
  
      NumC = ComputeNumCoeffs(FiltOrd);
      for(int k = 0; k&lt;11; k++)
      {
          printf("NumC is: %lf\n", NumC[k]);
      }
      //is A in matlab function and the numbers are correct
      DenC = ComputeDenCoeffs(FiltOrd, FrequencyBands[0], FrequencyBands[1]);
      for(int k = 0; k&lt;11; k++)
      {
          printf("DenC is: %lf\n", DenC[k]);
      }
      double y[5];
      double x[5]={1,2,3,4,5};
      filter(5, DenC, NumC, 5, x, y);    
      return 1;
  }
  </code></pre>
  
  <p>I get this resutls for my code :</p>
  
  <blockquote>
    <p>B= 1,0,-5,0,10,0,-10,0,5,0,-1
    A= 1.000000000000000,  -4.945988709743181,  13.556489496973796,  -24.700711850327743,
    32.994881546824828,  -33.180726698160655,  25.546126213403539,  -14.802008410165968,
    6.285430089797051,  -1.772929809750849,  0.277753012228403</p>
  </blockquote>
  
  <p>but if I want to test the coefficinets in same frequency band in MATLAB, I get the following results:</p>
  
  <blockquote>
  <pre><code>&gt;&gt; [B, A]=butter(5, [0.25,0.375])
  </code></pre>
    
    <p>B =   0.0002,         0,   -0.0008,         0,    0.0016,         0,   -0.0016,         0,  0.0008,         0,   -0.0002</p>
    
    <p>A =   1.0000,   -4.9460,   13.5565,  -24.7007,   32.9948,  -33.1806,   25.5461,  -14.8020,    6.2854,   -1.7729,    0.2778</p>
  </blockquote>
  
  <p>I have test this website :http://www.exstrom.com/journal/sigproc/ code, but the result is equal as mine, not matlab. anybody knows why? or how can I get the same result as matlab toolbox?</p>