几个简单的数据点平滑处理算法

最近在写一些数据处理的程序。经常需要对数据进行平滑处理。直接用FIR滤波器或IIR滤波器都有一个启动问题，滤波完成后总要对数据掐头去尾。因此去找了些简单的数据平滑处理的方法。

在一本老版本的《数学手册》中找到了几个基于最小二乘法的数据平滑算法。将其写成了C 代码，测试了一下，效果还可以。这里简单的记录一下，算是给自己做个笔记。

算法的原理很简单，以五点三次平滑为例。取相邻的5个数据点，可以拟合出一条3次曲线来，然后用3次曲线上相应的位置的数据值作为滤波后结果。简单的说就是 Savitzky-Golay 滤波器 。只不过Savitzky-Golay 滤波器并不特殊考虑边界的几个数据点，而这个算法还特意把边上的几个点的数据拟合结果给推导了出来。

不多说了，下面贴代码。首先是线性拟合平滑处理的代码. 分别为三点线性平滑、五点线性平滑和七点线性平滑。

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 ; } }

 来自CODE的代码片

linearSmooth.c

然后是利用二次函数拟合平滑。

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 ] + 12 * 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 ]) + 17 * 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 * in [ 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 ; } }

 来自CODE的代码片

quadraticSmooth.c

最后是三次函数拟合平滑。

/** * 五点三次平滑 * */ 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 ; } }

 来自CODE的代码片

cubicSmooth.c

上面的代码经过了简单的测试，可以放心使用。

最近还有人向我要9点线性平滑和11点线性平滑的系数。我简单算了算，把系数列在了后面。不过我觉得这两组系数意义不大，与其这样做还不如设计一个FIR 滤波器或者 SG 滤波器。

9点线性平滑

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点线性平滑

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]