最小二乘法三次多项式曲线拟合 算法 C++ 实现，该怎么处理

/******************************************

//参考 《常用算法程序集 （C语言描述 第三版)》

 

//最小二乘法

//x[n] y[n] 已知输入

//n输入点个数

//a[m] 返回m-1次拟合多项式的m个系数

//m  拟合多项式的项数,即拟合多项式的最高次为m-1

//dt[3] dt[0]返回拟合多项式与各数据点误差的平方和，

        dt[1]返回拟合多项式与各数据点误差的绝对值之和

        dt[2]返回拟合多项式与各数据点误差的绝对值的最大值

 

//

//拟合多项式的输出

//Y(x) = a0 + a1(x-X) + a2(x-X)^2 + …… am(x-X)^m

// 其中X为已知点x的平均值

******************************************/

 

#include "math.h"

void pir1(x,y,n,a,m,dt)

int n,m

double x[],y[],a[],dt[]

{ 

  int i,j,k

  double z,p,c,g,q,d1,d2,s[20],t[20],b[20]

  for (i=0 i<=m-1 i++) a[i]=0.0

  if (m>n) m=n

  if (m>20) m=20

  z=0.0

  for (i=0 i<=n-1 i++) z=z+x[i]/(1.0*n)

  b[0]=1.0 d1=1.0*n p=0.0 c=0.0

  for (i=0 i<=n-1 i++)

  { p=p+(x[i]-z) c=c+y[i]}

  c=c/d1 p=p/d1

  a[0]=c*b[0]

  if (m>1)

  { t[1]=1.0 t[0]=-p

  d2=0.0 c=0.0 g=0.0

  for (i=0 i<=n-1 i++)

  { q=x[i]-z-p d2=d2+q*q

  c=c+y[i]*q

  g=g+(x[i]-z)*q*q

  }

  c=c/d2 p=g/d2 q=d2/d1

  d1=d2

  a[1]=c*t[1] a[0]=c*t[0]+a[0]

  }

  for (j=2 j<=m-1 j++)

  { s[j]=t[j-1]

  s[j-1]=-p*t[j-1]+t[j-2]

  if (j>=3)

    for (k=j-2 k>=1 k--)

      s[k]=-p*t[k]+t[k-1]-q*b[k]

    s[0]=-p*t[0]-q*b[0]

    谨拍d2=0.0 c=0.0 g=0.0

    for (i=0 i<=n-1 i++)

    { q=s[j]

    for (k=j-1 k>=0 k--)

      q=q*(x[i]-z)+s[k]

   唤晌首 d2=d2+q*q c=c+y[i]*q

    g=g+(x[i]-z)*q*q

    }

    c=c/d2 p=g/d2 q=d2/d1

    d1=d2

    a[j]=c*s[j] t[j]=s[j]

    for (k=j-1 k>=0 k--)

    { a[k]=c*s[k]+a[k]

    b[k]=t[k] t[k]=s[k]

    }

  }

  dt[0]=0.0 dt[1]=0.0 dt[2]=0.0

  for (i=0 i<=n-1 和数i++)

  { q=a[m-1]

  for (k=m-2 k>=0 k--)

    q=a[k]+q*(x[i]-z)

  p=q-y[i]

  if (fabs(p)>dt[2]) dt[2]=fabs(p)

  dt[0]=dt[0]+p*p

  dt[1]=dt[1]+fabs(p)

  }

  return

}

任意波形的生成 (geneartion of arbitrary waveform) 在商业，军事等领域都有着重要的应用，诸如空间光通信 (free-space optics communication)， 高速信号处理 (high-speed signal processing)，雷达 (radar) 等。在任意波形生成后， 如何评估生成的任意波形 成为另外一个重要的话题。

假设有一组实验数据，已知他们之间的函数关系：y=f(x)，通过这些信息，需要确定函数中的一些参数项。例如，f 是一个线型函数 f(x)=k*x+b，那么参数 k 和 b 就是需要确春州纤定的值。如果这些参数用 p 表示的话，那么就需要找到一组 p 值使得如下公式中的 S 函数最小：

这种算法被称之为 最小二乘拟合 (least-square fitting)。scipy 中的子函数库 optimize 已迹雹经提供实现最小二乘拟合算法的函数 leastsq 。下面是 leastsq 函扒仿数导入的方式：

scipy.optimize.leastsq 使用方法

在 Python科学计算——Numpy.genfromtxt 一文中，使用 numpy.genfromtxt 对数字示波器采集的三角波数据导入进行了介绍，今天，就以 4GHz三角波 波形的拟合为案例介绍任意波形的拟合方法。

在 Python科学计算——如何构建模型？ 一文中，讨论了如何构建三角波模型。在标准三角波波形的基础上添加了 横向，纵向的平移和伸缩特征参数 ，最后添加了 噪声参数 模拟了三角波幅度参差不齐的随机性特征。但在波形拟合时，并不是所有的特征参数都要纳入考量，例如，噪声参数应是 波形生成系统 的固有特征，正因为它的存在使得产生的波形存在瑕疵，因此，在进行波形拟合并评估时，不应将噪声参数纳入考量，最终模型如下：

在调用 scipy.optimize.leastsq 函数时，需要构建误差函数：

有时候，为了使图片有更好的效果，需要对数据进行一些处理：

leastsq 调用方式如下：

合理的设置 p0 可以减少程序运行时间，因此，可以在运行一次程序后，用拟合后的相应数据对 p0 进行修正。

在对波形进行拟合后，调用 pylab 对拟合前后的数据进行可视化：

均方根误差 (root mean square error) 是一个很好的评判标准，它是观测值与真值偏差的平方和观测次数n比值的平方根，在实际测量中，观测次数n总是有限的，真值只能用最可信赖（最佳）值来代替.方根误差对一组测量中的特大或特小误差反映非常敏感，所以，均方根误差能够很好地反映出测量的精密度。

RMSE 用程序实现如下：

拟合效果，模型参数输出：

leastsq 函数适用于任何波形的拟合，下面就来介绍一些常用的其他波形：

---