emd分解.txt

EMD分解C语言算法

///////////////////////////////////：
 /*
 插值算法1：
 
        输入：
				n：整型变量，输入参数，节点个数
				x[]: 大小为n的一维实数组，存放样条节点：切x1<x2<x3……
				y[]: 大小为n的一维实数组，输入参数，存放样条节点上的函数值
				yp1：在第一个节点处的函数值的一阶导数
				ypn：在第n个节点处的函数值的一阶导数
		输出：
		        y2[]: 插值点的二阶导数
		          
 */  
//////////////////////////////////
void CCaculate::spline1(double x[], double y[], int n, 
					   double yp1,double ypn, double y2[])
{
     double *uk,aaa,sig,p,bbb,ccc,qn,un;
	 uk=new double [n];
	int i,k;
    if( yp1 > 9.9e+29 )
	{
        y2[1] = 0;
        uk[1] = 0;
	}
    else
	{
        y2[1] = -0.5;
        aaa = (y[2] - y[1]) / (x[2] - x[1]);
        uk[1] = (3.0 / (x[2] - x[1])) * (aaa - yp1);
    }
    for (i = 2; i<=n-1; i++)
	{
        sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
        p = sig * y2[i - 1] + 2.0;
        y2[i] = (sig - 1.0) / p;
        aaa = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
        bbb = (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
        ccc = x[i + 1] - x[i - 1];
        uk[i] = (6.0 * (aaa - bbb) / ccc - sig * uk[i - 1]) / p;
    }
    if (ypn > 9.9e+29 )
	{
        qn = 0.0;
        un = 0.0;
	}
    else
	{
        qn = 0.5;
        aaa = ypn - (y[n] - y[n - 1]) / (x[n] - x[n - 1]);
        un = (3.0/ (x[n] - x[n - 1])) * aaa;
    }
    y2[n] = (un - qn * uk[n - 1]) / (qn * y2[n - 1] + 1.0);
    for (k = n - 1;k>=1;k--)
        y2[k] = y2[k] * y2[k + 1] + uk[k];
	delete [] uk;


}
///////////////////////////////////：
 /*
 插值算法2：
 
        输入：
				n：整型变量，输入参数，节点个数
				xa[]: 大小为n的一维实数组，存放样条节点：切x1<x2<x3……
				ya[]: 大小为n的一维实数组，输入参数，存放样条节点上的函数值
				y2a[]: 插值算法1中所求的输出：插值点的二阶导数
				x：     需求的样条点
		输出：
		        y：输出所求的样条插值结果

注意：splint1:只调用一次，一旦调用后，保存输出，对任意的x，插值函数的值通过
      多次调用splint2（需要几次就调用几次）
         		          
 */  
//////////////////////////////////
void CCaculate::splint2(double xa[], double ya[], double y2a[], 
						int n, double& x, double& y)
{
	int klo,khi,k;
	double h,a,b,aaa,bbb;
    klo = 1;
    khi = n;
    loop:   if (khi - klo > 1 )
	{
		k = (khi + klo) / 2;
		if (xa[k] > x)
			khi = k;
		else
		{
			klo = k;
		}
           goto loop;
	}
    h = xa[khi] - xa[klo];
    if (h == 0 )
	{
	//	AfxMessageBox("输入的xa[]数组有问题");
		return;
	}
	a = (xa[khi] - x) / h;
	b = (x - xa[klo]) / h;
	aaa = a * ya[klo] + b * ya[khi];
	bbb = (a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi];
	y = aaa + bbb * h*h /6.0;

}
///////////////////////////////////：
 /*
 插值算法1：
 
        输入：
				n：整型变量，输入参数，节点个数
				x[]: 大小为n的一维实数组，存放样条节点：切x1<x2<x3……
				y[]: 大小为n的一维实数组，输入参数，存放样条节点上的函数值
				yp1：在第一个节点处的函数值的一阶导数
				ypn：在第n个节点处的函数值的一阶导数
		输出：
		        y2[]: 插值点的二阶导数
		          
 */  
//////////////////////////////////
void CCaculate::spline1(double x[], double y[], int n, 
					   double yp1,double ypn, double y2[])
{
     double *uk,aaa,sig,p,bbb,ccc,qn,un;
	 uk=new double [n];
	int i,k;
    if( yp1 > 9.9e+29 )
	{
        y2[1] = 0;
        uk[1] = 0;
	}
    else
	{
        y2[1] = -0.5;
        aaa = (y[2] - y[1]) / (x[2] - x[1]);
        uk[1] = (3.0 / (x[2] - x[1])) * (aaa - yp1);
    }
    for (i = 2; i<=n-1; i++)
	{
        sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
        p = sig * y2[i - 1] + 2.0;
        y2[i] = (sig - 1.0) / p;
        aaa = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
        bbb = (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
        ccc = x[i + 1] - x[i - 1];
        uk[i] = (6.0 * (aaa - bbb) / ccc - sig * uk[i - 1]) / p;
    }
    if (ypn > 9.9e+29 )
	{
        qn = 0.0;
        un = 0.0;
	}
    else
	{
        qn = 0.5;
        aaa = ypn - (y[n] - y[n - 1]) / (x[n] - x[n - 1]);
        un = (3.0/ (x[n] - x[n - 1])) * aaa;
    }
    y2[n] = (un - qn * uk[n - 1]) / (qn * y2[n - 1] + 1.0);
    for (k = n - 1;k>=1;k--)
        y2[k] = y2[k] * y2[k + 1] + uk[k];
	delete [] uk;


}
///////////////////////////////////：
 /*
 插值算法2：
 
        输入：
				n：整型变量，输入参数，节点个数
				xa[]: 大小为n的一维实数组，存放样条节点：切x1<x2<x3……
				ya[]: 大小为n的一维实数组，输入参数，存放样条节点上的函数值
				y2a[]: 插值算法1中所求的输出：插值点的二阶导数
				x：     需求的样条点
		输出：
		        y：输出所求的样条插值结果

注意：splint1:只调用一次，一旦调用后，保存输出，对任意的x，插值函数的值通过
      多次调用splint2（需要几次就调用几次）
         		          
 */  
//////////////////////////////////
void CCaculate::splint2(double xa[], double ya[], double y2a[], 
						int n, double& x, double& y)
{
	int klo,khi,k;
	double h,a,b,aaa,bbb;
    klo = 1;
    khi = n;
    loop:   if (khi - klo > 1 )
	{
		k = (khi + klo) / 2;
		if (xa[k] > x)
			khi = k;
		else
		{
			klo = k;
		}
           goto loop;
	}
    h = xa[khi] - xa[klo];
    if (h == 0 )
	{
	//	AfxMessageBox("输入的xa[]数组有问题");
		return;
	}
	a = (xa[khi] - x) / h;
	b = (x - xa[klo]) / h;
	aaa = a * ya[klo] + b * ya[khi];
	bbb = (a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi];
	y = aaa + bbb * h*h /6.0;

}