c++opencvsignal-processingdigital-filter

FIR filter design in C++ using coefficients from MATLAB, Filter not giving correct results


I am working on a project in OpenCV, as I am new to this field, I am facing a few problems with it. I have to design a 128th order FIR bandpass filter for which I have used the coefficients calculated through MATLAB. I found the following code online.

const int taps = 129;
double buffer[129] = {0.0};
int offset = 0;
double input;

double coefficients[] =

{
    0.0005,0.0004,0.0002,0.0000,-0.0001,-0.0000,0.0001,0.0004,0.0007,0.0011,
    0.0015,0.0017,0.0018,0.0016,0.0011,0.0006,0.0001,-0.0003,-0.0002,0.0003,
    0.0011,0.0022,0.0032,0.0038,0.0038,0.0031,0.0015,-0.0005,-0.0026,-0.0044,
    -0.0053,-0.0050,-0.0037,-0.0016,0.0007,0.0023,0.0027,0.0012,-0.0022,-0.0072,
    -0.0127,-0.0178,-0.0212,-0.0220,-0.0199,-0.0153,-0.0093,-0.0036,-0.0001,-0.0003,
    -0.0053,-0.0147,-0.0274,-0.0408,-0.0519,-0.0573,-0.0545,-0.0419,-0.0197,0.0102,
    0.0442,0.0779,0.1064,0.1254,0.1321,0.1254,0.1064,0.0779,0.0442,0.0102,-0.0197,
    -0.0419,-0.0545,-0.0573,-0.0519,-0.0408,-0.0274,-0.0147,-0.0053,-0.0003,-0.0001,-0.0036,
    -0.0093,-0.0153,-0.0199,-0.0220,-0.0212,-0.0178,-0.0127,-0.0072,-0.0022,0.0012,0.0027,
    0.0023,0.0007,-0.0016,-0.0037,-0.0050,-0.0053,-0.0044,-0.0026,-0.0005,0.0015,0.0031,
    0.0038,0.0038,0.0032,0.0022,0.0011,0.0003,-0.0002,-0.0003,0.0001,0.0006,0.0011,
    0.0016,0.0018,0.0017,0.0015,0.0011,0.0007, 0.0004,0.0001,-0.0000,-0.0001,0.0000,
    0.0002,0.0004,0.0005
};


double filter( double input)

{
    double output = 0;


    for(int i= taps-1; i>0 ; i-- )
    {
       buffer[i] = buffer[i-1];
    }

       buffer[0] = input;


    for(int j = 0; j<taps; j++ )
    {

        output += (coefficients[j]*buffer[j]);

    }
    return output; 
}

I try to use the filter as described in the code snippet below. filteredTrajectory is the vector of vector, where I want to store the filtered values. YCoordinates is a vector of vector containing the data to be filtered. The dimension of YCoordinates is [1000x266]

// Calling FIR filter function filter
// filtering all 266 values across 1000 frames.

const int numFrames = 1000;

//featuresPrevious.size() = 266 (predefined fixed value)

vector<vector<double>> filteredTrajectory;
filteredTrajectory.resize(numFrames-1);

   for (int i = 0; i<numFrames-1; i++)

    {

       filteredTrajectory[i].resize(featuresPrevious.size());

       for (int j = 0; j<featuresPrevious.size(); j++)

        {

            filteredTrajectory[i][j] = filter (Ycoordinates[i][j]); 

     }

    }

When I run this code, the filtered values do not match with those of MATLAB. However if I filter only one value throughout all 1000 frames, I get matching results with MATLAB (see code below).

// Filtering value at index Ycoordinates[i][0] across all frames 

filteredTrajectory.resize(numFrames-1);

   for (int i = 0; i<numFrames-1; i++)

    {

     filteredTrajectory[i].resize(1);

     filteredTrajectory[i][0] = filter (Ycoordinates[i][0]);

} 

What could be the possible flaws in the code which provide wrong results when filtering all trajectories inside the loop but matching results when filtered individually.


Solution

  • Matlab's filter is a one dimensional filter. When provided an input which is a two dimensional matrix, it filters each column independently (the default if the dim parameter is not set, or set to 1; setting dim to 2 would do the same but on each row), essentially resetting the FIR state in between each processed column (or row).

    As such, to achieve the correct result you must also reset the filter state which you store in buffer between each set of values to be filtered by setting resetting the buffer to zeros:

    void reset_buffer()
    {
      for(int i=0; i<taps ; i++ )
      {
        buffer[i] = 0.0;
      }
    }
    

    Also, the filter must see each values belonging to the same set sequentially (i.e. not interleaving the values to the filtering function), so the order of the loop is important. Depending on whether you want to filter along the column or row dimension of your data, you can use either of the following function:

    // Filter along columns, similar to Matlab's filter(coefficients, 1, data, [], 1)
    void filterDim1()
    {
      for (int i = 0; i<numFrames-1; i++)
      {
        filteredTrajectory[i].resize(featuresPrevious.size());
      }
      for (int j = 0; j<featuresPrevious.size(); j++)
      {
        for (int i = 0; i<numFrames-1; i++)
        {
          filteredTrajectory[i][j] = filter (Ycoordinates[i][j]); 
        }
        reset_buffer();
      }
    }
    // Filter along rows, similar to Matlab's filter(coefficients, 1, data, [], 2)
    void filterDim2()
    {
      for (int i = 0; i<numFrames-1; i++)
      {
        filteredTrajectory[i].resize(featuresPrevious.size());
        for (int j = 0; j<featuresPrevious.size(); j++)
        {
          filteredTrajectory[i][j] = filter (Ycoordinates[i][j]); 
        }
        reset_buffer();
      }
    }