Matlab smooth函数原理

由于项目上要用到平滑一维数组数据，参考Matlab  smooth函数转成c++代码

// x,g均为数组，具体内容略
plot(x,g);hold on,plot(x,smooth(g,50),'r');
z1 = (g1-smooth(g1,50)');  figure,plot(x,z1,'.-')

蓝色为平滑前，红色为平滑后

为了要找到缺陷，即灰度值突变很大地方，可以平滑前后相减，注意这里平滑窗宽尽量选大，选择原则是较小甚至不影响缺陷突变的地方

    平滑前后相减

                     

举例

例如g=[1 2 1 5 1 1 1 3 1 1 1]

smooth(g,5)'

ans =

1.0000    1.3333    2.0000    2.0000    1.8000    2.2000    1.4000    1.4000    1.4000    1.0000    1.0000

smooth(g,10)'

ans =

    1.0000    1.3333    2.0000    1.7143    1.7778    1.7778    1.6667    1.2857    1.4000    1.0000    1.0000

以窗宽10为例(因为除的奇数，所以最大为9)  (matlab 下标以1开始,而vc从0开始)

令b=(k-1)/2;

当i小于b时    

Y[0]=g[0] =1

Y[1]=(g[0]+g[1]+g[2])/3=4/3=1.3333

Y[2]=(g[0]+g[1]+g[2]+g[3]+g[4])/5=2

Y[3]=(g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6])/7=1.7143 

Y[3]= (g[3+(-3)]+g[3+(-2)]+g[3+(-1)]+g[3+(0)]+g[3+1]+g[3+2]+g[3+3])/(2*3+1)

那么从j=-i：i    (这里k是变动的)

Y[i]=0;

   For( j=-i;j<i+1;i++)

{

      Y[i]+=g[i+j]

}

Y[i]=y[i]/(2*i+1);

当i大于等于时 ，同时length(g)-i>b  (这里k是固定的)

如果n为偶数，那么应该是k=n-1；然后范围是-(k-1)/2: (k-1)/2   b=(k-1)/2

Y[4]=(g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6] +g[7]+g[8])/9=1.7778

Y[5]=(g[1]+g[2]+g[3]+g[4]+g[5]+g[6] +g[7]+g[8] +g[9])/9=1.7778

Y[6]=(g[2]+g[3]+g[4]+g[5]+g[6]+g[7]+g[8] +g[9]+g[10])/9=1.6667

Y[6]= (g[6+(-4)]+g[6+(-3)]+g[6+(-2)]+g[6+(-1)]+g[6]+g[7]+g[8] +g[9]+g[10])/9…..

Y[i]=(g[i+(-b)]+g[i+(-b+1)]+g[i+(-b+2)]+g[i+(-b+3)]+g[i+(-b+4)]+g[i+(-b+5)]+g[i+(-b+6)] +g[i+(-b+7)]+g[i+(-b+8)])/k

-b：b即

If(i>b)||(i=b)

{

   Y[i]=0;

   For( j=-b;j<b+1;b++)

{

      Y[i]+=g[i+j]

}

Y[i]=y[i]/k;}

当i大于b  同时 length(g)-i<=b   从后往前推

Y[7]=( g[4]+ g[5]+g[6] +g[7]+g[8] +g[9]+g[10])/7=1.2857

Y[8]=( g[6] +g[7]+g[8] +g[9]+g[10])/5=1.4

Y[9]=( g[8] +g[9]+g[10])/3=1

Y[10]=(g[10] =1

当i=7时，i1=10-7=3

那么Y[7]= (g[7+(-3)]+g[7+(-2)]+g[7+(-1)]+g[7+(0)]+g[7+1]+g[7+2]+g[7+3])/(2*3+1)

当i=8时，i1=2

那么Y[8]= (g[8+(-2)]+g[8+(-1)]+g[8+(0)]+g[8+1]+g[8+2])/(2*2+1)

令 i1=(n-1)-i,那么从-i1：i1

Y[i]=0;

   For( j=-i1;j<i1+1;i1++)

{

      Y[i]+=g[i+j]

}

Y[i]=y[i]/(2*i1+1);

 

C++代码如下，nn为数组元素个数；

            int b=(k-1)/2;
            vector<double> Z(nn);
            if(nn>k)
            {
                for(int i=0;i<nn;i++)
                {

                    if(i<b)
                    {
                        Z[i]=0;
                        for(int j=-i;j<i+1;j++)
                        {
                            Z[i]+=D[i+j];
                        }
                        Z[i]=Z[i]/(2*i+1);
                    }
                        
                  else  if(((i>b)||(i=b))&((nn-i)>b))
                    {
                        Z[i]=0;
                        for( int j=-b;j<b+1;j++)
                        {
                            Z[i]+=D[i+j];
                        }
                        Z[i]=Z[i]/k;
                    }
                    else
                    {
                        Z[i]=0;
                        int i1=(nn-1)-i;
                        for( int j=-i1;j<i1+1;j++)
                        {
                            Z[i]+=D[i+j];
                        }
                        Z[i]=Z[i]/(2*i1+1);
                    }
                }
            }

  平滑前

 平滑后

 平滑前后相减