mchgenbg.c

speech signal process tools

		fprintf (stderr, "measure did not decrease\n");
	/* do interpolation by 1/2 */
	/* try to find a decreased measure */
	    rttflag = iinterp ();
	    if (rttflag == 1) {
		cnt++;
		if (moflag)
		    fprintf (stderr, "no decrease in measure during interpolation\n");
	    }
	/* iinterp() found a decreased measure */
	    if (cvtest (R_measure, Rtemp_measure))
		return (5);	/* conclude iteration */
	    return (0);
	}
	if (Rnew_measure < R_measure) {
	    if (cvtest (R_measure, Rnew_measure)) {
		Rnew_R ();
		return (5);	/* conclude iteration */
	    }
	    Rnew_R ();
	}
	Rnew_R0 ();		/* use Rnew for next iteration */
	return (rtflag);
    }
    fprintf (stderr, "trouble ckRnew\n");
    exit (1);
 /* NOTREACHED */
}
cvtest(R_measure,Rnew_measure)
/* convergence test on measure change */
register double R_measure,Rnew_measure;
{
    register double temp;
/*
		"Shankar's fix"
	Note this change hase been made to avoid division by zero.
	If given sample covariance matrix is toeplitz, division by zero
	could occur here also.
*/
    if (Rnew_measure > cvfact) {
	temp = (R_measure-Rnew_measure) / R_measure;
	temp = (temp >= 0) ? temp : -temp;
    }
    else
	temp = Rnew_measure;
    if (temp < cvfact) {
	if (moflag)
	    fprintf (stderr,
		    "Stopping due to convergence, rel_ratio=%g cvfact=%g\n",
		    temp, cvfact);
	return (1);
    }
    return (0);
}

inv_det(tR,itR,siginv,isiginv,pdet)
	/* return inverse and natural log of determinant */
double *tR,*itR,*siginv,*isiginv,*pdet;
{
    int     rtflag;

    if (scflag) {
	if (rflag) {
	    rtflag = ltoepinv (tR, siginv, pdet, qp);
	    return (rtflag);
	}
	if (cflag) {
	    rtflag = lctoepinv (tR, itR, siginv, isiginv, pdet, qp);
	    return (rtflag);
	}
    }
    if (mdflag || mcflag) {
	makesigma (tR, itR, sigtemp, isigtemp, Pi);
	if (rflag) {
	/* test for positive definitness with cholesky but use more
	   accurate orthogonalizaton technique to get siginv if positive
	   def */


	/* 
	 dmmove(sigtemp,sig2temp,qp);
	 rtflag=lsyminv(sig2temp,siginv,sig1temp,qp,pdet);
	 */

	    rtflag = lsyminv (sigtemp, siginv, sig1temp, qp, pdet);
	    if (rtflag)
		return (rtflag);

	/* 
	 rtflag=ldmsyminv(sigtemp,siginv,pdet,qp,0);
	 */
	    return (rtflag);
	}
	if (cflag) {
	    fprintf (stderr, "no ln|R| available\n");
	    exit (1);
	/* NOTREACHED */
	/* 
	 rtflag=lhdmsym(siginv,isiginv,siginv,isiginv,pdet,qp);
	 return(rtflag);
	 */
	}
    }
 /* JTB: lint thinks this is reached, anyway */
    return (rtflag);
}

double
domeasure(siginv,isiginv,det)
	/* det is the natural logrithm of the determinant of R */
register double det,*siginv,*isiginv;
{
    register double measure;
    double  hdmtrace (), dmtrace ();

    if (rflag) {
	measure = dmtrace (siginv, S, qp) - qp;
	measure = measure + det - logdetS;
    }
    if (cflag) {
	measure = hdmtrace (siginv, isiginv, S, iS, qp) - qp;
	measure = measure + det - logdetS;
    }

    return (measure);
}

do_algorithm()
/* make an Rnew[] from the R0[] using the algorithm */
/* R0 is known to be positive definite */
{
    register int    rtflag;
    double  dum;

    doBi (siginv, Pi, Bi, qp);	/* make Bi */

 /* sig1temp=siginv*S*siginv */
 /* sigtemp=siginv*S*siginv - siginv */

    if (anderson && (initflag == 2 || initflag == 3 || it_cnt > 1) && rflag) {
				/* do not allow anderson on first
				   iteration unless initflag=2 or 3, do
				   not allow anderson in complex case */
	dmadd (sigtemp, sig1temp, sig2temp, qp);
				/* make 2*siginv*S*siginv-siginv */
	doAij (sig2temp, siginv, Aij);/* make Aij */
    }
    else {
	if (rflag)
	    doAij (siginv, siginv, Aij);/* make Aij */
	if (cflag)
	    docAij (siginv, isiginv, siginv, isiginv, Aij);/* make Aij */
    }

 /* 
  rtflag=ldmsyminv(Aij,Ainv,&dum,tNi,0);
  if(rtflag==1)
  { fprintf(stderr,"Aij singular\n"); return(3); }
  doRnew(Ainv,Bi);
  */

 /* solves Aij*x=y, Aij pos definite */
 /* tBi=x, Bi=y, ttBi=tttBi=temp vector, M and Ainv are temp matrices */

    if (posdefsol (Aij, M, Ainv, tBi, Bi, ttBi, tttBi, tNi)) {
	if (moflag)
	    fprintf (stderr, "Aij npd from posdefsol\n");
	if (anderson) {		/* try burg algorithm */
	    fprintf (stderr, "trying burg algorithm\n");
	    anderson = 0;
	    rtflag = do_algorithm ();
	    return (rtflag);
	}
	else {
	    rtflag = ldmsyminv (Aij, Ainv, &dum, tNi, 0);
	    if (rtflag == 1) {
		fprintf (stderr, "Aij singular from ldmsyminv\n");
		return (3);
	    }
	    doRnew (Ainv, Bi);
	    return (0);
	}
    }

    else
	DooRnew (tBi, Rnew, R0);/* breakup tBi[] into Rnew[] and iRnew[] 
				*/
    return (0);

}

DooRnew(tBi,Rnew,R0)/* breakup tBi[] into Rnew[] and iRnew[] */
register double tBi[],Rnew[],R0[];
{
    register int    i, fNi;
    register double dum;

    fNi = Ni;

    if (anderson && (initflag == 2 || initflag == 3 || it_cnt > 1) && rflag) {
	for (i = 0; i < fNi; i++) {
	    dum = tBi[i] + R0[i];
	    Rnew[i] = dum;
	}
    }
    else {
	dvmove (tBi, Rnew, fNi);
	if (cflag)
	    dvmove (tBi + fNi, iRnew + 1, tNi - fNi);
    }
    return;
}

doAspec(siginv,matrix,Aij)
register double matrix[],siginv[],Aij[];
{
    register double At;
    register int    tab;
    register short *p;

    p = Aijtab;
    while (*p != -2) {
	At = 0;
	while ((tab = (*p)) != -1) {
	    p++;
	    At += siginv[tab] * matrix[*p];
	    p++;
	}
	p++;
	*(Aij + *p) = At;
	p++;
	*(Aij + *p) = At;
	p++;
    }
    return;
}

doAij(matrix,siginv,Aij)/* make Aij real case */
register double matrix[],siginv[],Aij[];
/* Aij=trace Pi*siginv*Pj*matrix */
/* Aij= (Pi)mn*(siginv)np*(Pj)pq*(matrix)qm with Einstein summation convent. */
{
    register int    i, j, iM, jM, fNi, ftNi;
    register double At;
    double  dooAij ();

 /* 
  cntAij++;
  if(cntAij>3)return;
  */
    if (iiflag == 2) {		/* unroll table has been made */
	if (qdim == qdimsave && pdim == pdimsave) {
				/* unroll table has been made and it is
				   for the current dimension */
	    doAspec (siginv, matrix, Aij);
	    return;
	}			/* do fast Aij */
	else {
	/* unroll table is not for current dimension */
	    freedp (&Aijtab);	/* free storage */
	    iiflag = 0;		/* set to count for new table size */
	}
    }

 /* setup for unroll */
 /* first time (iiflag==0) just count for size of unroll table, second
    time (iiflag==1) assign indices to table */
    ii = 0;
    pAijtab = Aijtab;

    fNi = Ni;
    ftNi = tNi;

 /* do Pi*siginv*Pj*matrix */
 /* see comments on Aij in Aijxxx */

    for (i = iM = 0; i < fNi; i++, iM += ftNi) {
	make (Pi, i);
	for (j = i, jM = iM + i; j < fNi; j++, jM += ftNi) {
	    if (i == j) {
		At = dooAij (siginv, Pi, matrix, Pi);
		Aij[iM + j] = At;
	    }
	    else {
		make (Pj, j);
		At = dooAij (siginv, Pi, matrix, Pj);
		Aij[jM] = Aij[iM + j] = At;
	    }
	    setunroll (jM, iM + j);
	}
    }
    finishunroll ();
    return;
}

docAij(matrix,imatrix,siginv,isiginv,Aij)/* do complex case Aij */
register double matrix[],imatrix[],siginv[],isiginv[],Aij[];
/* Aij=trace Pi*siginv*Pj*matrix */
/* Aij= (Pi)mn*(siginv)np*(Pj)pq*(matrix)qm with Einstein summation convent. */
{
    register int    i, j, iM, jM, fNi, ftNi;
    int     NixtNi;
    register double At;
    double  doocAij ();

    fNi = Ni;
    ftNi = tNi;

 /* do Pi*siginv*Pj*matrix - Pi*isiginv*Pj*imatrix */
 /* see comments on Aij */

    for (i = 0, iM = 0; i < fNi; i++, iM += ftNi) {
	make (Pi, i);
	for (j = i, jM = iM + i; j < fNi; j++, jM += ftNi) {
	    if (i == j) {
		At = doocAij (siginv, isiginv, Pi, matrix, imatrix, Pi, 0);
		Aij[iM + j] = At;
	    }
	    else {
		make (Pj, j);
		At = doocAij (siginv, isiginv, Pi, matrix, imatrix, Pj, 0);
		Aij[jM] = Aij[iM + j] = At;
	    }
	}
    }

    NixtNi = fNi * ftNi;

 /* do Pi*siginv*iPj*imatrix + Pi*isiginv*iPj*matrix */

    for (i = 0, iM = 0; i < fNi; i++, iM += ftNi) {
	make (Pi, i);
	for (j = fNi, jM = NixtNi; j < ftNi; j++, jM += ftNi) {
	    make (Pj, j - fNi + 1);
	    At = doocAij (siginv, isiginv, Pi, matrix, imatrix, Pj, 1);
	    Aij[jM + i] = Aij[iM + j] = At;
	}
    }

 /* do iPi*siginv*iPj*matrix - iPi*isiginv*iPj*imatrix */

    for (i = fNi, iM = NixtNi; i < ftNi; i++, iM += ftNi) {
	make (Pi, i - fNi + 1);
	for (j = i, jM = iM; j < ftNi; j++, jM += ftNi) {
	    if (i == j) {
		At = doocAij (siginv, isiginv, Pi, matrix, imatrix, Pi, 2);
		Aij[iM + j] = At;
	    }
	    else {
		make (Pj, j - fNi + 1);
		At = doocAij (siginv, isiginv, Pi, matrix, imatrix, Pj, 2);
		Aij[jM + i] = Aij[iM + j] = At;
	    }
	}
    }
    return;
}

finishunroll()
{
    if (iiflag == 1) {		/* making table */
    /* signal end of table */
	*pAijtab = -2;
	pAijtab++;
    }
    ii++;
    if (iiflag == 0) {		/* not yet made table */
	if (allocsp (&Aijtab, ii))
	    fprintf (stderr, "cannot allocate for Aijtab\n");
	else {
	    iiflag = 1;		/* set to make table */
	    qdimsave = qdim;	/* save dimensions of table */
	    pdimsave = pdim;
	}
    }
    else if (iiflag == 1)	/* if makeing table, set table finished */
	iiflag = 2;		/* table finished */
    return;
}

setunroll(i1,i2)/* set up addresses for unroll or just count for allocation */
register int i1,i2;
{
    if (iiflag == 1) {		/* making table */
	*pAijtab = -1;		/* signal new element */
	pAijtab++;
	*pAijtab = i1;
	pAijtab++;
	*pAijtab = i2;
	pAijtab++;
    }
    else {
    /* doing count only */
	ii++;
	ii++;
	ii++;
    }
    return;
}

double
dooAij(siginv,Pi,matrix,Pj)
register double *matrix,*siginv;
register short *Pj,*Pi;

/* Aij=trace Pi*siginv*Pj*matrix */
/* matrix is siginv or 2*siginv*S*siginv-siginv */
/* Aij= (Pi)mn*(siginv)np*(Pj)pq*(matrix)qm with Einstein summation convent. */
{
    register int    m, n, nM, p, q, qM, qpr;
    register short *pPj;
    register double At;

    qpr = qp;
    At = 0;
    for (m = 0; m < qpr; m++) {
	for (n = 0, nM = 0; n < qpr; n++, Pi++, nM += qpr) {
	    if (!*Pi)
		continue;
	    pPj = Pj;
	    for (p = 0; p < qpr; p++) {
		for (q = 0, qM = m; q < qpr; q++, qM += qpr, pPj++) {
		    if (!*pPj)
			continue;
		    At += *(siginv + nM + p) * *(matrix + qM);
		    if (iiflag) {/* making table */
			*pAijtab = nM + p;
			pAijtab++;
			*pAijtab = qM;
			pAijtab++;
		    }
		    else {	/* just counting */
			ii++;
			ii++;
		    }
		}
	    }
	}
    }
    return (At);
}

double
doocAij(siginv,isiginv,Pi,matrix,imatrix,Pj,flag)/* complex case */
register double *matrix,*imatrix,*siginv,*isiginv;
register short *Pj,*Pi;
register int flag;/* indicates the terms in the complex expansion */

/* Aij=trace Pi*siginv*Pj*matrix */
/* matrix is siginv or 2*siginv*S*siginv-siginv */
/* Aij= (Pi)mn*(siginv)np*(Pj)pq*(matrix)qm with Einstein summation convent. */
{
    register int    m, n, nM, p, q, qM, qpr;
    register short *pPj;
    register double At;
    double  doooA ();

    qpr = qp;
    At = 0;
    for (m = 0; m < qpr; m++) {
	for (n = 0, nM = 0; n < qpr; n++, Pi++, nM += qpr) {
	    if (!*Pi)
		continue;
	    pPj = Pj;
	    for (p = 0; p < qpr; p++) {
		for (q = 0, qM = m; q < qpr; q++, qM += qpr, pPj++) {
		    if (!*pPj)
			continue;
		    At += doooA (siginv, isiginv, matrix, imatrix, nM + p,
			    qM, flag, n, m, p, q);
		}
	    }
	}
    }
    return (At);
}

double /* do multiplication of matrices, see comments in Aijxxx */
doooA(siginv,isiginv,matrix,imatrix,i1,i2,flag,n,m,p,q)
register double siginv[],isiginv[],matrix[],imatrix[];
register int i1,i2,flag,n,m,p,q;
{
    register double term;

    if (!flag) {
	term = *(siginv + i1) * *(matrix + i2) - *(isiginv + i1) * *(imatrix + i2);
	return (term);
    }
    else if (flag == 1) {
	term = (siginv[i1] * imatrix[i2] + isiginv[i1] * matrix[i2]);
	if (q > p)
	    return (-term);
	else
	    return (term);
    }
    else {		/*(flag == 2)*/
    /* fact1=(n>m)?1:-1; fact2=(q>p)?1:-1; */
    /* here complex Pi are positive one in upper triangle */

	term = siginv[i1] * matrix[i2] - isiginv[i1] * imatrix[i2];
	if (((n > m) && (q > p)) || ((m > n) && (p > q)))
	    return (-term);
	else
	    return (term);
    }
}
/*
Aijxxx see also the discussion in doBi

Aij=trace[Pi*siginv*M*Pj] where M=2*siginv*S*siginv-siginv or =siginv

siginv*S*siginv is hermitian
M is hermitian since it is composed of hermitian matrices

Aij is hermitian,

Aij=trace[Pi*siginv*Pj*M]=trace[Pj*M*Pi*siginv]=trace[(M*Pj)'*(siginv*Pi)']=
trace[siginv*Pi*M*Pj]complex conjugate=trace[Pj*siginv*Pi*M]complex conjugate=
Aji complex conjugate, Q.E.D.

For M=siginv, Aij is symmetric and real
For the real case, Aij is symmtric and real

The Pi and Pj are symmetric, the iPi and iPj are anti-symmetric

trace [ Pi*siginv*Pj*matrix ]=

trace [ Pi*siginv*Pj*matrix ]+		flag=0
trace [ iPi*siginv*Pj*matrix ]+		anti-symmetric, trace zero
trace [ Pi*siginv*iPj*matrix ]+		anti-symmetric, trace zero
trace [ iPi*siginv*iPj*matrix ]+	flag=2
trace [ Pi*siginv*Pj*imatrix ]+		anti-symmetric, trace zero
trace [ iPi*siginv*Pj*imatrix ]+	term not computed directly, use Aji 
trace [ Pi*siginv*iPj*imatrix ]+	flag=1
trace [ iPi*siginv*iPj*matrix ]+	imaginary, does not contribute
trace [ Pi*isiginv*Pj*matrix ]+		anti-symmetric, trace zero
trace [ iPi*isiginv*Pj*matrix ]+	term not computed directly, use Aji 
trace [ Pi*isiginv*iPj*matrix ]+	flag=1
trace [ iPi*isiginv*iPj*matrix ]+	imaginary, does not contribute
trace [ Pi*isiginv*Pj*imatrix ]+	flag=0
trace [ iPi*isiginv*Pj*imatrix ]+	imaginary, does not contribute
trace [ Pi*isiginv*iPj*imatrix ]+	imaginary, does not contribute
trace [ iPi*isiginv*iPj*imatrix ]	flag=2

The Aij matrix is tNi by tNi, the Ni by Ni portion is due to the real
Pi, the others parts involve the imaginary iPi.

*/
doBi(siginv,Pi,Bi,qp)	/* make Bi */
register double Bi[],siginv[];
register short *Pi;
register int qp;
/* Bi[i]=trace siginv*S*siginv*Pi */
{
	register int i,fNi,ftNi;
	register double Bt;
	double spdmtrace();

	fNi=Ni;
	ftNi=tNi;
	if(rflag)
	{
	dmmult(S,siginv,sigtemp,qp);/* make S*siginv */
	dmmult(siginv,sigtemp,sig1temp,qp);/* make siginv*S*siginv */
	if(anderson && (initflag==2 || initflag==3 || it_cnt>1) && rflag)
		{
		dmsub(sig1temp,siginv,sigtemp,qp);/* make siginv*S*siginv-siginv */
		/* sig1temp=siginv*S*siginv */
		/* sigtemp=siginv*S*siginv - siginv */
		}
	for(i=0;i<fNi;i++)
		{
		make(Pi,i);
		if(anderson && (initflag==2 || initflag==3 || it_cnt>1) && rflag)
			Bt=spdmtrace(sigtemp,Pi,qp);/*trace(siginv*S*siginv-siginv)Pi */
		else
			Bt=spdmtrace(sig1temp,Pi,qp);/*trace(siginv*S*siginv*Pi) */
		Bi[i]=Bt;
		}
	return;
	}

	/* complex case