基于b样条插值的最小二乘正弦拟合c代码

正弦拟合c代码以及正弦拟合算法原理一、正弦拟合基本原理

二、完整c代码

#include

#include

#include

#define n 201

#define pi 3.1416

using namespace std;

double byangtiao(double x[], double y[], int count);//样条插值

double sinnihe(double x[], double y[], int count);//正弦拟合

int gssd(double *a,double *b,double *x,int count,double eps,int iter);//线性方程组求解

double *pu;

double *pv;

int main()

{

double x[8]={21.098901098901099, 59.340659340659343 ,93.626373626373621, 135.82417582417582 ,195.16483516483515, 250.54945054945054 ,303.29670329670330 ,340.2197821978020};

double

y[8]={3036.1252441406250 ,3036.1218261718750 ,3036.1184082031250 ,3036.1103515625000 ,3036.1135253906250 ,3036.1235351562500, 3036.1252441406250 ,3036.1252441406250};

for(int i=0;i<9;i++)

x[i]=x[i]*pi/180;

byangtiao(x, y, (sizeof x) / (sizeof x[0]));

sinnihe(pu, pv, ((sizeof x) / (sizeof x[0]) - 1)*n);

delete[]pu; delete[]pv;

return 0;

}

double byangtiao(double a[], double b[], int count)

{

double *x; double*y;int i=0;int j=0;

x = new double[count + 2];

y = new double[count + 2];

x[0] = a[0]; y[0] = b[0];

x[count + 1] = a[count - 1]; y[count + 1] = b[count - 1];

for ( i = 0; i

{

x[i + 1] = a[i]; y[i + 1] = b[i];

}

pu = new double[((count-1)*n)];

pv = new double[((count-1)*n)];

double t[n] = { 0 };

int k = 0;

for ( i = 0; i

{

t[i] =(i)*0.005;

}

for ( i = 0; i < count - 1; i++)

{

for ( j = 0; j < n; j++)

{

pu[k] = x[i] * 1 / 6 * (-t[j] * -t[j] * -t[j] +3* t[j] * t[j] - 3 * t[j] + 1) + x[i + 1] *

1 / 6 * (3 * t[j] * t[j] * t[j] - 6 * t[j] * t[j] + 4) + x[i + 2] * 1 / 6 *

(-3 * t[j] * t[j] * t[j] + 3 * t[j] * t[j] + 3 * t[j] + 1) + x[i + 3] *

1 / 6 * t[j] * t[j] * t[j];

pv[k] = y[i] * 1 / 6 * (-t[j] * -t[j] * -t[j] + 3*t[j] * t[j] - 3 * t[j] + 1) + y[i + 1] *

1 / 6 * (3 * t[j] * t[j] * t[j] - 6 * t[j] * t[j] + 4) + y[i + 2] * 1 / 6 *

(-3 * t[j] * t[j] * t[j] + 3 * t[j] * t[j] + 3 * t[j] + 1) + y[i + 3] *

1 / 6 * t[j] * t[j] * t[j];

k++;

}

}

delete[]x; delete[]y;

return 0;

}

int gssd(double *a,double *b,double *x,int count,double eps,int iter)

{

int i,j,piter=0;

double fr=1.0,t,*x2;

x2=(double*)malloc(n*sizeof(double));

if(x2==NULL)

{

return(0);

}

for(i=0;i

if(fabs(a[i*count+i])<1e-12)

{

return(0);

}

for(i=0;i

x[i]=0.0;

while((fr>eps)&&(piter

{

fr=0.0;

for(i=0;i

{

t=0.0;

for(j=0;j

t=t+a[i*count+j]*x2[j];

for(j=i+1;j

t=t+a[i*count+j]*x[j];

x2[i]=(b[i]-t)/a[i*count+i];

t=fabs(x2[i]-x[i])/(1.0+fabs(x[i]));

if(t>fr)

fr=t;

}

for(i=0;i

x[i]=x2[i];

piter++;

}

return(piter);

}

double sinnihe(double x[], double y[], int count) {

int k = 1; int i, j;

double A[9] = { 0 };

double B[3] = { 0 };

for (i = 0; i < count; i++)

{

A[0] = A[0]+1;

A[1] = A[1] + sin(pu[i]);

A[2] = A[2] + cos(pu[i]);

A[3] = A[3] + sin(pu[i]);

A[4] = A[4] + sin(pu[i])*sin(pu[i]);

A[5] = A[5] + sin(pu[i])*cos(pu[i]);

A[6] = A[6] + cos(pu[i]);

A[7] = A[7] + sin(pu[i])*cos(pu[i]);

A[8] = A[8] + cos(pu[i])*cos(pu[i]);

B[0] = B[0]+pv[i];

B[1] = B[1] + pv[i] * sin(pu[i]);

B[2] = B[2] + pv[i] * cos(pu[i]);

}

for(i=0;i<9;i++)

cout<

for(i=0;i<3;i++)

cout<

double w = 0,h=0,d=0,C[3]={0};

gssd(A,B,C,3,0.000001,30);

d = C[0];

w = atan(C[2] / C[1]);

h = C[1] / cos(w);

cout << " 幅度值为：" << h <

return 0;

}