Algorithm for smoothing

Question

I wrote this code for smoothing of a curve . It takes 5 points next to a point and adds them and averages it .

/* Smoothing */
void smoothing(vector<Point2D> &a)
{
    //How many neighbours to smooth
    int NO_OF_NEIGHBOURS=10;
    vector<Point2D> tmp=a;
    for(int i=0;i<a.size();i++)
    {

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

        }
        else
        {
            for(int j=1;j<NO_OF_NEIGHBOURS;j++)
            {
                a.at(i).x+=tmp.at(i-j).x;
                a.at(i).y+=tmp.at(i-j).y;
            }
            a.at(i).x/=NO_OF_NEIGHBOURS;
            a.at(i).y/=NO_OF_NEIGHBOURS;
        }

    }

}

But i get very high values for each point, instead of the similar values to the previous point . The shape is maximized a lot , what is going wrong in this algorithm ?

Answer 1

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.

Answer 2

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
            }

Answer 3

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

Answer 4

The current-value of the point is used twice: once because you use += and once if y==0. So you are building the sum of eg 6 points but only dividing by 5. This problem is in both the IF and ELSE case. Also: you should check that the vector is long enough otherwise your ELSE-case will read at negative indices.

Following is not a problem in itself but just a thought: Have you considered to use an algorithm that only touches every point twice?: You can store a temporary x-y-value (initialized to be identical to the first point), then as you visit each point you just add the new point in and subtract the very-oldest point if it is further than your NEIGHBOURS back. You keep this "running sum" updated for every point and store this value divided by the NEIGHBOURS-number into the new point.

Answer 5

You make addition with point itself when you need to take neighbor points - just offset index by 1:

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
 }