genburg.c

speech signal process tools

				}
			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;

	temp=(R_measure-Rnew_measure)/Rnew_measure;
	temp=(temp>=0)?temp:-temp;
	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);
				/*
				rtflag=lhdmsym(siginv,isiginv,siginv,isiginv,pdet,qp);
				return(rtflag);
				*/
				}
			}
	fprintf(stderr,"Error in inv_det\n");
	exit(1);
	/*NOTREACHED*/
}
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 /* if (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
	Bi is real

	Let M' be the complex transpose of M,
	if M=M', then M is hermitian, siginv and S and Pi are hermitian,