tpeplz.cpp

c经典算法(包括源码)

void toeplz(double r[], double x[], double y[], int n)
{
	int k,j,m,m1,m2;
	double sxn,sd,sgn1,shn,sgd,pt1,pt2,qt1,qt2,pp,qq;
    const int nmax = 100;
    double g[100], h[100];
    if (r[n] == 0.0) 
	{
		cout<<"levinson method fails: singular principal minor"<<endl;
		return;
	}
	x[1] = y[1] / r[n];
    if (n == 1)
	{
		return;
	}
    g[1] = r[n - 1] / r[n];
    h[1] = r[n + 1] / r[n];
    for (m = 1; m<=n; m++)
	{
        m1 = m + 1;
        sxn = -y[m1];
        sd = -r[n];
        for (j = 1; j<=m; j++)
		{
            sxn = sxn + r[n + m1 - j] * x[j];
            sd = sd + r[n + m1 - j] * g[m - j + 1];
        }
        if (sd == 0.0)
		{
			cout<<"levinson method fails: singular principal minor"<<endl;
			return;
		}
        x[m1] = sxn / sd;
        for (j = 1; j<=m; j++)
		{
            x[j] = x[j] - x[m1] * g[m - j + 1];
        }
        if (m1 == n)
		{
			return;
		}
        sgn1 = -r[n - m1];
        shn = -r[n + m1];
        sgd = -r[n];
        for (j = 1; j<=m; j++)
		{
            sgn1 = sgn1 + r[n + j - m1] * g[j];
            shn = shn + r[n + m1 - j] * h[j];
            sgd = sgd + r[n + j - m1] * h[m - j + 1];
        }
        if ((sd == 0.0) || (sgd == 0.0))
		{
			cout<<"levinson method fails: singular principal minor"<<endl;
			return;
		}
        g[m1] = sgn1 / sgd;
        h[m1] = shn / sd;
        k = m;
        m2 = (m + 1) / 2;
        pp = g[m1];
        qq = h[m1];
        for (j = 1; j<=m2; j++)
		{
            pt1 = g[j];
            pt2 = g[k];
            qt1 = h[j];
            qt2 = h[k];
            g[j] = pt1 - pp * qt2;
            g[k] = pt2 - pp * qt1;
            h[j] = qt1 - qq * pt2;
            h[k] = qt2 - qq * pt1;
            k = k - 1;
        }
    }
    cout<<"never get here"<<endl;
}