Several simple smoothing algorithms for data points __ Numerical calculation

Recently in the process of writing some data processing. It is often necessary to smooth the data. FIR Filter directly or IIR filter has a startup problem, the filter is completed after the total data Qiatouquwei. So I went to find some simple ways to deal with the data smoothing.

Several data smoothing algorithms based on least squares are found in an old version of the Math Handbook. It is written in C code, tested a bit, the effect is also OK. Here is a simple record, as a note for yourself.

The principle of the algorithm is simple, with 5.3 times smoothing as an example. The adjacent 5 data points can be fitted out with a 3-time curve, and then the data values of the corresponding positions on the 3 curves are used as filtering results. Simply put, the Savitzky-golay filter. But the Savitzky-golay filter does not specifically consider the boundaries of several data points, and this algorithm also deliberately put the edge of several points of data fitting results are deduced.

Not much said, the following code. The first is the linear fit smoothing Process code. They are three-point linear smoothing, five-point linear smoothing and seven-point linear smoothing respectively.

void LinearSmooth3 (double in[], double out[], int N) {int i;
        if (N < 3) {for (i = 0; I <= N-1; i++) {out[i] = In[i];

        } else {Out[0] = (5.0 * IN[0] + 2.0 * in[1]-in[2])/6.0;
        for (i = 1; I <= N-2; i++) {Out[i] = (In[i-1] + in[i] + in[i + 1])/3.0;
    } Out[n-1] = (5.0 * IN[N-1] + 2.0 * in[n-2]-in[n-3])/6.0;
    } void LinearSmooth5 (double in[], double out[], int N) {int i;
        if (N < 5) {for (i = 0; I <= N-1; i++) {out[i] = In[i];
        } else {Out[0] = (3.0 * In[0] + 2.0 * in[1] + in[2]-in[4])/5.0;
        OUT[1] = (4.0 * IN[0] + 3.0 * In[1] + 2 * in[2] + in[3])/10.0;  for (i = 2; I <= N-3; i++) {Out[i] = (In[i-2] + in[i-1] + in[i] + in[i + 1] + In[i + 2])
        /5.0; } Out[n-2] = (4.0 * IN[N-1] + 3.0 * In[n-2] + 2 * in[n-3] + in[n-4])/10.0;
    Out[n-1] = (3.0 * In[n-1] + 2.0 * in[n-2] + in[n-3]-in[n-5])/5.0;
    } void LinearSmooth7 (double in[], double out[], int N) {int i;
        if (N < 7) {for (i = 0; I <= N-1; i++) {out[i] = In[i]; } else {out[0] = (13.0 * In[0] + 10.0 * IN[1] + 7.0 * IN[2] + 4.0 * In[3] + in[4

        ]-2.0 * IN[5]-5.0 * in[6])/28.0;

        OUT[1] = (5.0 * IN[0] + 4.0 * IN[1] + 3 * in[2] + 2 * in[3] + in[4]-in[6])/14.0; OUT[2] = (7.0 * In[0] + 6.0 * in [1] + 5.0 * IN[2] + 4.0 * IN[3] + 3.0 * IN[4] + 2.0 * in[5] + in[6]

        /28.0;  for (i = 3; I <= N-4 i++) {Out[i] = (in[i-3] + in[i-2] + in[i-1] + in[i] + in[i + 1] +
        In[i + 2] + In[i + 3])/7.0; } Out[n-3] = (7.0 * In[n-1] + 6.0 * in [N- 2] + 5.0 * IN[N-3] + 4.0 * IN[N-4] + 3.0 * IN[N-5] + 2.0 * in[n-6] + in[n-7])/28.0; Out[n-2] = (5.0 * IN[N-1] + 4.0 * IN[N-2] + 3.0 * In[n-3] + 2.0 * in[n-4] + in[n

        -5]-in[n-7])/14.0;  Out[n-1] = (13.0 * in[n-1] + 10.0 * In[n-2] + 7.0 * In[n-3] + 4 * in[n-4] + in[n-5]-2
    * IN[N-6]-5 * in[n-7])/28.0;
 }
}

Then we use the two-time function to fit the smoothing.

void QuadraticSmooth5 (double in[], double out[], int N) {int i;
        if (N < 5) {for (i = 0; I <= N-1; i++) {out[i] = In[i];
        } else {out[0] = (31.0 * In[0] + 9.0 * IN[1]-3.0 * IN[2]-5.0 * IN[3] + 3.0 * in[4])/35.0;
        OUT[1] = (9.0 * In[0] + 13.0 * in[1] + * IN[2] + 6.0 * IN[3]-5.0 *in[4])/35.0;
                      for (i = 2; I <= N-3 i++) {Out[i] = (-3.0 * (in[i-2) + In[i + 2]) +
        12.0 * (In[i-1] + in[i + 1]) + * In[i])/35.0; } Out[n-2] = (9.0 * In[n-1] + 13.0 * In[n-2] + 12.0 * In[n-3] + 6.0 * IN[N-4]-5.0 * in[n-5])/35.
        0;
    Out[n-1] = (31.0 * in[n-1] + 9.0 * In[n-2]-3.0 * IN[N-3]-5.0 * IN[N-4] + 3.0 * in[n-5])/35.0;
    } void QuadraticSmooth7 (double in[], double out[], int N) {int i;
     if (N < 7) {for (i = 0; I <= N-1; i++) {       Out[i] = In[i]; } else {out[0] = (32.0 * in[0] + 15.0 * IN[1] + 3.0 * IN[2]-4.0 * IN[3]-6.0

        * IN[4]-3.0 * IN[5] + 5.0 * in[6])/42.0;

        OUT[1] = (5.0 * IN[0] + 4.0 * IN[1] + 3.0 * IN[2] + 2.0 * in[3] + in[4]-in[6])/14.0; OUT[2] = (1.0 * IN[0] + 3.0 * in [1] + 4.0 * IN[2] + 4.0 * IN[3] + 3.0 * IN[4] + 1.0 * IN[5]-2.0 *
        [6])/14.0;
                       for (i = 3; I <= N-4 i++) {Out[i] = ( -2.0 * (in[i-3) + In[i + 3]) +
        3.0 * (In[i-2] + in[i + 2]) + 6.0 * (In[i-1] + in[i + 1]) + 7.0 * In[i])/21.0; } Out[n-3] = (1.0 * IN[N-1] + 3.0 * in [N-2] + 4.0 * IN[N-3] + 4.0 * In[n-4] + 3.

        0 * IN[N-5] + 1.0 * IN[N-6]-2.0 * in[n-7])/14.0; Out[n-2] = (5.0 * IN[N-1] + 4.0 * IN[N-2] + 3.0 * In[n-3] + 2.0 * in[N-4] + in[n-5]-in[n-7])/14.0; Out[n-1] = (32.0 * in[n-1] + 15.0 * In[n-2] + 3.0 * In[n-3]-4.0 * IN[N-4]-6.0 * in[n
    -5]-3.0 * IN[N-6] + 5.0 * in[n-7])/42.0;
 }
}

Finally, three times function fitting smoothing.

/** * 5.3 times Smooth * * */void CubicSmooth5 (double in[], double out[], int N) {int i;
    if (N < 5) {for (i = 0; I <= N-1; i++) out[i] = In[i];
        else {out[0] = (69.0 * in[0] + 4.0 * IN[1]-6.0 * IN[2] + 4.0 * in[3]-in[4])/70.0;
        OUT[1] = (2.0 * in[0] + 27.0 * IN[1] + 12.0 * IN[2]-8.0 * IN[3] + 2.0 * in[4])/35.0;  for (i = 2; I <= N-3; i++) {Out[i] = ( -3.0 * (in[i-2) + In[i + 2]) + 12.0 * (in[i-1) + in[i
        + 1]) + 17.0 * In[i])/35.0; 
        } Out[n-2] = (2.0 * in[n-5]-8.0 * IN[N-4] + 12.0 * In[n-3] + 27.0 * In[n-2] + 2.0 * in[n-1])/35.0;
    Out[n-1] = (-in[n-5] + 4.0 * IN[N-4]-6.0 * in[n-3] + 4.0 * In[n-2] + 69.0 * in[n-1])/70.0;
} return;
    } void CubicSmooth7 (double in[], double out[], int N) {int i;
    if (N < 7) {for (i = 0; I <= N-1; i++) {out[i] = In[i];    } else {out[0] = (39.0 * in[0] + 8.0 * IN[1]-4.0 * IN[2]-4.0 * In[3] + 1
        .0 * In[4] + 4.0 * IN[5]-2.0 * in[6])/42.0; OUT[1] = (8.0 * In[0] + 19.0 * IN[1] + 16.0 * IN[2] + 6.0 * IN[3]-4.0 * in[4]-7.0* in[5] + 4.0 * in
        [6])/42.0; OUT[2] = ( -4.0 * In[0] + 16.0 * in [1] + 19.0 * IN[2] + 12.0 * IN[3] + 2.0 * IN[4]-4.0 * IN[5] + 1.0
        * In[6])/42.0;
                       for (i = 3; I <= N-4 i++) {Out[i] = ( -2.0 * (in[i-3) + In[i + 3]) +
        3.0 * (In[i-2] + in[i + 2]) + 6.0 * (In[i-1] + in[i + 1]) + 7.0 * In[i])/21.0; } out[n-3] = ( -4.0 * in[n-1] + 16.0 * in [N-2] + 19.0 * In[n-3] + 12.0 * In[n-4]
        + 2.0 * IN[N-5]-4.0 * IN[N-6] + 1.0 * in[n-7])/42.0; Out[n-2] = (8.0 * In[n-1] + 19.0 * In[n-2] + 16.0 * In[n-3] + 6.0 * IN[N-4]-4.0 * IN[N-5]-7.0 * IN[N-6] + 4.0 * in[n-7])/42.0;  Out[n-1] = (39.0 * in[n-1] + 8.0 * IN[N-2]-4.0 * IN[N-3]-4.0 * IN[N-4] + 1.0 * in[n-
    5] + 4.0 * IN[N-6]-2.0 * in[n-7])/42.0;
 }
}

The code above is a simple test that can be used with ease.

I've been asked to get a 9-point linear smoothing and 11-point linear smoothing coefficients. I simply calculate, the factor is listed in the back. But I don't think the two sets of coefficients mean much, so it's better to design an FIR filter or a SG filter.

9 Point Linear Smoothing

Yy[0]= (17*y[0] + 14*y[1] + 11*y[2] + 8*y[3] + 5*y[4] + 2*y[5]-y[6]-4*y[7]-7*y[8])/45;
Yy[1]= (56*y[0] + 47*y[1] + 38*y[2] + 29*y[3] + 20*y[4] + 11*y[5] + 2*y[6]-7*y[7]-16*y[8])/180;
Yy[2]= (22*y[0] + 19*y[1] + 16*y[2] + 13*y[3] + 10*y[4] + 7*y[5] + 4*y[6] + y[7]-2*y[8])/90;
Yy[3]= (32*y[0] + 29*y[1] + 26*y[2] + 23*y[3] + 20*y[4] + 17*y[5] + 14*y[6] + 11*y[7] + 8*y[8])/180;

Yy[4]= (Y[0] + y[1] + y[2] + y[3] + y[4] + y[5] + y[6] + y[7] + y[8])/9;

Yy[5]= (8*y[0] + 11*y[1] + 14*y[2] + 17*y[3] + 20*y[4] + 23*y[5] + 26*y[6] + 29*y[7] + 32*y[8])/180;
Yy[6]= ( -2*y[0] + y[1] + 4*y[2] + 7*y[3] + 10*y[4] + 13*y[5] + 16*y[6] + 19*y[7] + 22*y[8])/90;
Yy[7]= ( -16*y[0]-7*y[1] + 2*y[2] + 11*y[3] + 20*y[4] + 29*y[5] + 38*y[6] + 47*y[7] + 56*y[8])/180;
Yy[8]= ( -7*y[0]-4*y[1]-y[2] + 2*y[3] + 5*y[4] + 8*y[5] + 11*y[6] + 14*y[7] + 17*y[8])/45;

11 Point Linear Smoothing

yy[0]=0.3181818181818182*
Y[0] + 0.2727272727272727*y[1] + 0.22727272727272727*y[2] +
0.18181818181818182*Y[3] + 0.13636363636363635*y[4] +
0.09090909090909091*Y[5] + 0.045454545454545456*y[6] +
1.6152909428804953*^-17*Y[7]-0.045454545454545456*y[8]-
0.09090909090909091*Y[9]-0.13636363636363635*Y[10]

yy[1]=0.27272727272727276*
Y[0] + 0.23636363636363633*y[1] + 0.19999999999999998*y[2] +
0.16363636363636364*Y[3] + 0.12727272727272726*y[4] +
0.09090909090909091*Y[5] + 0.05454545454545455*y[6] +
0.018181818181818202*Y[7]-0.01818181818181818*y[8]-
0.054545454545454536*Y[9]-0.09090909090909088*Y[10]

yy[2]=0.2272727272727273*
Y[0] + 0.19999999999999996*y[1] + 0.17272727272727267*y[2] +
0.14545454545454542*Y[3] + 0.11818181818181815*y[4] +
0.09090909090909091*Y[5] + 0.06363636363636363*y[6] +
0.03636363636363638*Y[7] + 0.009090909090909094*y[8]-
0.01818181818181816*Y[9]-0.0454545454545454*Y[10]

yy[3]=0.18181818181818188*
Y[0] + 0.16363636363636355*y[1] + 0.1454545454545454*y[2] +
0.12727272727272723*Y[3] + 0.10909090909090904*y[4] +
0.0909090909090909*Y[5] + 0.07272727272727272*y[6] +
0.054545454545454564*Y[7] + 0.036363636363636376*y[8] +
0.01818181818181823*Y[9] + 8.326672684688674*^-17*y[10]

yy[4]=0.13636363636363644*
Y[0] + 0.12727272727272718*y[1] + 0.11818181818181811*y[2] +
0.10909090909090903*Y[3] + 0.09999999999999995*y[4] +
0.0909090909090909*Y[5] + 0.08181818181818182*y[6] +
0.07272727272727275*Y[7] + 0.06363636363636364*y[8] +
0.05454545454545459*Y[9] + 0.04545454545454555*y[10]

Yy[5]=0.090909090909091*
Y[0] + 0.09090909090909077*y[1] + 0.09090909090909083*y[2] +
0.09090909090909083*Y[3] + 0.09090909090909084*y[4] +
0.0909090909090909*Y[5] + 0.09090909090909091*y[6] +
0.09090909090909094*Y[7] + 0.09090909090909093*y[8] +
0.090909090909091*Y[9] + 0.09090909090909102*y[10]

yy[6]=0.04545454545454558*
Y[0] + 0.0545454545454544*y[1] + 0.06363636363636352*y[2] +
0.07272727272727264*Y[3] + 0.08181818181818173*y[4] +
0.0909090909090909*Y[5] + 0.1*y[6] + 0.10909090909090911*y[7] +
0.11818181818181821*Y[8] + 0.12727272727272737*y[9] +
0.13636363636363652*Y[10]

Yy[7]=1.1102230246251565*^-\
16*y[0] + 0.01818181818181802*y[1] + 0.03636363636363624*y[2] +
0.05454545454545445*Y[3] + 0.07272727272727264*y[4] +
0.0909090909090909*Y[5] + 0.10909090909090909*y[6] +
0.12727272727272732*Y[7] + 0.1454545454545455*y[8] +
0.16363636363636372*Y[9] + 0.181818181818182*y[10]

yy[8]=-0.0454545454545453*
Y[0]-0.018181818181818354*y[1] + 0.009090909090908955*y[2] +
0.03636363636363624*Y[3] + 0.06363636363636355*y[4] +
0.09090909090909088*Y[5] + 0.11818181818181818*y[6] +
0.1454545454545455*Y[7] + 0.17272727272727273*y[8] +
0.2000000000000001*Y[9] + 0.22727272727272746*y[10]

yy[9]=-0.09090909090909072*
Y[0]-0.05454545454545473*y[1]-0.018181818181818327*y[2] +
0.01818181818181805*Y[3] + 0.05454545454545444*y[4] +
0.09090909090909088*Y[5] + 0.12727272727272726*y[6] +
0.16363636363636366*Y[7] + 0.2*y[8] + 0.23636363636363647*y[9] +
0.27272727272727293*Y[10]

yy[10]=-0.1363636363636362*
Y[0]-0.09090909090909116*y[1]-0.04545454545454561*y[2]-
1.6653345369377348*^-16*Y[3] + 0.04545454545454533*y[4] +
0.09090909090909088*Y[5] + 0.13636363636363635*y[6] +
0.18181818181818188*Y[7] + 0.2272727272727273*y[8] +
0.27272727272727293*Y[9] + 0.3181818181818184*y[10]