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

  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 }