Algorithm for smoothing

What it looks like you have here is a bass-ackwards implementation of a finite impulse response (FIR) filter that implements a boxcar window function. Thinking about the problem in terms of DSP, you need to filter your incoming vector with NO_OF_NEIGHBOURS equal FIR coefficients that each have a value of 1/NO_OF_NEIGHBOURS. It is normally best to use an established algorithm rather than reinvent the wheel.

Here is a pretty scruffy implementation that I hammered out quickly that filters doubles. You can easily modify this to filter your data type. The demo shows filtering of a few cycles of a rising saw function (0,.25,.5,1) just for demonstration purposes. It compiles, so you can play with it.

#include <iostream>
#include <vector>

using namespace std;

class boxFIR
{
    int numCoeffs; //MUST be > 0
    vector<double> b; //Filter coefficients
    vector<double> m; //Filter memories

public:
    boxFIR(int _numCoeffs) :
    numCoeffs(_numCoeffs)
    {
        if (numCoeffs<1)
            numCoeffs = 1; //Must be > 0 or bad stuff happens

        double val = 1./numCoeffs;
        for (int ii=0; ii<numCoeffs; ++ii) {
            b.push_back(val);
            m.push_back(0.);
        }
    }    

    void filter(vector<double> &a)
    {
        double output;

        for (int nn=0; nn<a.size(); ++nn)
        {
            //Apply smoothing filter to signal
            output = 0;
            m[0] = a[nn];
            for (int ii=0; ii<numCoeffs; ++ii) {
                output+=b[ii]*m[ii];
            }

            //Reshuffle memories
            for (int ii = numCoeffs-1; ii!=0; --ii) {
                m[ii] = m[ii-1];
            }                        
            a[nn] = output;
        }
    }


};

int main(int argc, const char * argv[])
{
    boxFIR box(1); //If this is 1, then no filtering happens, use bigger ints for more smoothing

    //Make a rising saw function for demo
    vector<double> a;
    a.push_back(0.); a.push_back(0.25); a.push_back(0.5); a.push_back(0.75); a.push_back(1.);
    a.push_back(0.); a.push_back(0.25); a.push_back(0.5); a.push_back(0.75); a.push_back(1.);
    a.push_back(0.); a.push_back(0.25); a.push_back(0.5); a.push_back(0.75); a.push_back(1.);
    a.push_back(0.); a.push_back(0.25); a.push_back(0.5); a.push_back(0.75); a.push_back(1.);

    box.filter(a);

    for (int nn=0; nn<a.size(); ++nn)
    {
        cout << a[nn] << endl;
    }
}

Up the number of filter coefficients using this line to see a progressively more smoothed output. With just 1 filter coefficient, there is no smoothing.

boxFIR box(1);

The code is flexible enough that you can even change the window shape if you like. Do this by modifying the coefficients defined in the constructor.

Note: This will give a slightly different output to your implementation as this is a causal filter (only depends on current sample and previous samples). Your implementation is not causal as it looks ahead in time at future samples to make the average, and that is why you need the conditional statements for the situation where you are near the end of your vector. If you want output like what you are attempting to do with your filter using this algorithm, run the your vector through this algorithm in reverse (This works fine so long as the window function is symmetrical). That way you can get similar output without the nasty conditional part of algorithm.


in following block:

            for(int j=0;j<NO_OF_NEIGHBOURS;j++)
            {
                a.at(i).x=a.at(i).x+a.at(i+j).x;
                a.at(i).y=a.at(i).y+a.at(i+j).y;
            }

for each neighbour you add a.at(i)'s x and y respectively to neighbour values.

i understand correctly, it should be something like this.

            for(int j=0;j<NO_OF_NEIGHBOURS;j++)
            {
                a.at(i).x += a.at(i+j+1).x
                a.at(i).y += a.at(i+j+1).y
            }

Filtering is good for 'memory' smoothing. This is the reverse pass for the learnvst's answer, to prevent phase distortion:

for (int i = a.size(); i > 0; --i)
{
    // Apply smoothing filter to signal
    output = 0;
    m[m.size() - 1] = a[i - 1];

    for (int j = numCoeffs; j > 0; --j) 
        output += b[j - 1] * m[j - 1];

    // Reshuffle memories
    for (int j = 0; j != numCoeffs; ++j) 
        m[j] = m[j + 1];

    a[i - 1] = output;
}

More about zero-phase distortion FIR filter in MATLAB: http://www.mathworks.com/help/signal/ref/filtfilt.html

Tags:

Algorithm

C++

Signal Processing

Filter

Recent Posts