Smoothing algorithm

Question

I wrote this code to smooth the curve. He takes 5 points next to the point and adds them and averages.

/* 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(ij).x; a.at(i).y+=tmp.at(ij).y; } a.at(i).x/=NO_OF_NEIGHBOURS; a.at(i).y/=NO_OF_NEIGHBOURS; } } } 

But I get very high values ​​for each point, and not similar values ​​for the previous point. The form is maximized, what happens in this algorithm?

Answer 1

in the next 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 neighbor, you add a.at (i) x and y respectively to neighboring 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 2

What you see here is the bass implementation of a finite impulse response (FIR) filter that implements the boxcar window function . Thinking about the problem in terms of DSP, you need to filter your incoming vector with equal FIR coefficients NO_OF_NEIGHBOURS , each of which has a value of 1/NO_OF_NEIGHBOURS . It is usually better to use the established algorithm, rather than reinvent the wheel.

Here is a rather messy implementation that I quickly developed, and the filters are doubled. You can easily change this to filter your data type. The demo shows the filtering of several cycles of the upward saw function (0, .25, .5,1) for demonstration purposes only. 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; } } 

Increase the number of filter coefficients using this line to see a smoother result. There is no smoothing with one filter coefficient.

 boxFIR box(1); 

The code is flexible enough so you can even change the shape of the window if you want. Do this by changing the coefficients defined in the constructor.

Note: this will give a slightly different result of your implementation, since it is a causal filter (it depends only on the current selection and previous samples). Your implementation is not causal, as it expects the future in future samples in order to get the average value, and therefore you need conditional operators for the situation when you are approaching the end of the vector. If you want to get a result similar to what you are trying to do with your filter using this algorithm, run your vector using this algorithm in the reverse order (this works fine while the window function is symmetrical). Thus, you can get a similar conclusion without the unpleasant conditional part of the algorithm.

Answer 3

Filtering is good for smoothing memory. This is the return pass for 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; } 

For more on zero phase FIR filters in MATLAB: http://www.mathworks.com/help/signal/ref/filtfilt.html

Answer 4

The current point value is used twice: once because you use += and once if y==0 . Thus, you build the amount, for example, 6 points, but only divide by 5. This problem occurs both in the case of IF and in ELSE. Also: you should check that the vector is long enough, otherwise your ELSE case will read with negative indices.

This is not a problem in itself, but just a thought: do you think that you are using an algorithm that applies only to each point twice: you can save the temporary value xy (initialized so that it is identical to the first point), then as you visit each point, you simply add a new point and subtract the oldest point, if it is larger than your NEIGHBOURS . You keep this “current amount” updated for each point and keep this value divided by NEIGHBOURS -number in the new point.

Answer 5

You make an addition with the point itself, when you need to take neighboring points - just shift the 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 } 

Answer 6

This works great for me:

 for (i = 0; i < lenInput; i++) { float x = 0; for (int j = -neighbours; j <= neighbours; j++) { x += input[(i + j <= 0) || (i + j >= lenInput) ? i : i + j]; } output[i] = x / (neighbours * 2 + 1); }