样条之最小二乘算法求多项式

核心代码：

  1 // 使用最小二乘算法求多项式
  2 void YcLeastSquaresFitSpline::CalculateMultinomialValues(const void* valuesPtr, int stride, int n, int m, float* a) const
  3 {
  4     memset(a, 0, sizeof(float)*m);
  5
  6     float xStep = 1.0f/(n - 1);
  7
  8     int i,j,k;
  9     float z,p,c,g,q,d1,d2,s[20],t[20],b[20];
 10     for (i=0; i<=m-1; i++)
 11     {
 12         a[i]=0.0f;
 13     }
 14
 15     if (m>n)
 16     {
 17         m=n;
 18     }
 19     if (m>20)
 20     {
 21         m=20;
 22     }
 23
 24     z=0.0f;
 25     for (i=0; i<=n-1; i++)
 26     {
 27         z=z+xStep*i/(1.0f*n);
 28     }
 29
 30     b[0]=1.0f;
 31     d1=1.0f*n;
 32     p=0.0f;
 33     c=0.0f;
 34     for (i=0; i<=n-1; i++)
 35     {
 36         p=p+(xStep*i-z);
 37         c=c+YfGetFloatValue(valuesPtr, stride, i);
 38     }
 39     c=c/d1;
 40     p=p/d1;
 41     a[0]=c*b[0];
 42     if (m>1)
 43     {
 44         t[1]=1.0f;
 45         t[0]=-p;
 46         d2=0.0f;
 47         c=0.0f;
 48         g=0.0f;
 49
 50         for (i=0; i<=n-1; i++)
 51         {
 52             q=xStep*i-z-p;
 53             d2=d2+q*q;
 54             c=c+YfGetFloatValue(valuesPtr, stride, i)*q;
 55             g=g+(xStep*i-z)*q*q;
 56         }
 57         c=c/d2;
 58         p=g/d2;
 59         q=d2/d1;
 60         d1=d2;
 61         a[1]=c*t[1];
 62         a[0]=c*t[0]+a[0];
 63     }
 64
 65     for (j=2; j<=m-1; j++)
 66     {
 67         s[j]=t[j-1];
 68         s[j-1]=-p*t[j-1]+t[j-2];
 69         if (j>=3)
 70         {
 71             for (k=j-2; k>=1; k--)
 72             {
 73                 s[k]=-p*t[k]+t[k-1]-q*b[k];
 74             }
 75         }
 76
 77         s[0]=-p*t[0]-q*b[0];
 78         d2=0.0f;
 79         c=0.0f;
 80         g=0.0f;
 81         for (i=0; i<=n-1; i++)
 82         {
 83             q=s[j];
 84             for (k=j-1; k>=0; k--)
 85             {
 86                 q=q*(xStep*i-z)+s[k];
 87             }
 88             d2=d2+q*q;
 89             c=c+YfGetFloatValue(valuesPtr, stride, i)*q;
 90             g=g+(xStep*i-z)*q*q;
 91         }
 92
 93         c=c/d2;
 94         p=g/d2;
 95         q=d2/d1;
 96         d1=d2;
 97         a[j]=c*s[j];
 98         t[j]=s[j];
 99         for (k=j-1; k>=0; k--)
100         {
101             a[k]=c*s[k]+a[k];
102             b[k]=t[k];
103             t[k]=s[k];
104         }
105     }
106 }

切图：

相关软件的下载地址为：http://files.cnblogs.com/WhyEngine/TestSpline.zip