genburg.c

speech signal process tools

/*
The generalized Burg technique

Programmed by Daniel L. Wenger
408-425-5790 Santa Cruz, California

See:
"Estimation of Structured Covariance Matrices", J.P.Burg, D.G.Luenberger,
D.L.Wenger; Proceedings of the IEEE, Vol. 70, No. 9 September 1982

"Estimation of Covariance Matrices Which Are Linear Combinations or
Whose Inverses Are Linear Combinations of Given Matrices", T.W.Anderson

This program uses the algorithms discussed in the above papers to
find the best estimate for the covariance matrix of single channel real or
complex, multidimensional real or two channel real problems.

Input is an Hermitian sample covariance matrix 
Output is an Hermitian Toeplitz (block Toeplitz) matrix

sigma_in is the address of the real part of the sample matrix
isigma_in is the address of the imaginary part of the sample matrix

sigma_out is the address of the real part of the Toeplitze matrix
isigma_out is the address of the imaginary part of the Hermitian matrix

sigma_out and isigma_out are also used to pass the initial guess
              to the routine if initflag=1

*qd=qdim and *pd=pdim are integers

in the single channel case, qdim is the filter order plus one and pdim is one.

in the multi-dimensional case, qdim is the dimension of the block Toeplitz 
	matrix and pdim is the dimension of the individual element matrices.

in the multi-channel case, qdim is the filter order plus one and pdim
	is the number of channels ( restricted to 2 )

in all cases, the dimension of the input matrix is qdim times pdim

*c is a flag   =0  for a real matrix
	       =1  for a complex matrix

*sflag is a flag   =1  for the single channel case
	       =2  for the multi-dimensional case
               =3  for the multi-channel case

*mo=moflag is a flag  =0  for no monitoring
                      =1  for monitoring

if *inflag=initflag=0, initial guess is the identity matrix
if *inflag=initflag=1, initial guess is in sigma_out and possibly isigma_out
if *inflag=initflag=2, initial guess is the first projection of sigma_in
if *inflag=initflag=3, initial guess is John's average + first projection

the final distance or measure is returned in pdist, *pdist is a double

*small_in=small is a double used as a limit in the interpolation search
	        a reasonable value is 1e-8 to 1e-12

*Maxit_in=Maxit is a cutoff for maximum number of iterations with no decrease of
            	measure, a reasonable value is 4

a return flag is returned in returnflag
a value of zero is returned upon no decrease in measure after Maxit attempts
a value of one is returned upon Rinit or Rnew singular or Rinit neg definite
a value of two is returned upon sigma_in or Rinit non-pos-definite
a value of three is returned upon Aij singular
a value of four is returned upon unsuccessful interpolation, no pos-def Rnew
a value of five is returned upon measure ratio convergence test
a value of six is returned upon insufficient storage allocation

*convfact=cvfact is a convergence factor,
                 (R_measure-Rnew_measure)/Rnew_measure<convfact

*/


/* 
**********************************************************************
This material by no means is up to the SPS coding practices and
standards. It has been kludged together from other sources. Beware.
**********************************************************************
*/


/*
| This material contains proprietary software of Entropic Processing, Inc.   
| Any reproduction, distribution, or publication without the the prior	   
| written permission of Entropic Processing, Inc. is strictly prohibited.
| Any public distribution of copies of this work authorized in writing by
| Entropic Processing, Inc. must bear the notice			
|								
|              "Copyright 1986 Entropic Processing, Inc."
*/
#ifdef SCCS
static char *sccs_id = "@(#)genburg.c	1.1 5/21/87";
#endif



#include<stdio.h>
int qdim=0;
int pdim=1;
int rflag=1;	/* real flag */
int cflag=0;	/* complex */
int scflag=1;	/* single channel flag */
int mdflag=0;	/* multi-dimensional flag */
int mcflag=0;	/* multi-channel flag */
static int moflag = 0;	/* monitor flag */
static double small=1e-8;	/* convergence criterian */
static int cnt=1;		/* iterations for which measure does not decrease */
static int Maxit=4;	/* maximum iteration limit */
int Ni=0;
int tNi=0;
int qp=0;
int qpxqp=0;
static int it_cnt=0;/* total iteration count */

short *Pii=0,*Pi=0,*Pj=0,*Ptaba=0,*Ptabb=0,*Ptabc=0;

static double *R=0,*iR=0;	/* solution with the lowest measure todate */
static double R_det=0;		/* natural log of determinant of R */
static double R_measure=0;	/* measure of R */

static double *R0=0,*iR0=0;	/* last solution of algorithm or searches */
static double R0_det=0;		/* natural log of determinant of R0 */
static double R0_measure=0;	/* measure of R0 */

static double *Rnew=0,*iRnew=0;	/* current solution of algorithm or searches */
static double Rnew_det=0;		/* natural log of determinant of Rnew */
static double Rnew_measure=0;	/* measure of Rnew */

static double *del=0,*idel=0;	/* difference between Rnew and R0 */

static double *siginv=0;
static double *isiginv=0;
static double *sigtemp=0,*isigtemp=0;
static double *sig1temp=0,*isig1temp=0;
static double *sig2temp=0,*isig2temp=0;
static double *Aij=0;
static double *Ainv=0,*M=0;
static double *Bi=0,*tBi=0,*ttBi=0,*tttBi=0;

static double *S=0;	/* pointer to sigma */
static double *iS=0;	/* pointer to isigma */
static double logdetS=0;

/* relating to unroll of doAij */
static int ii=0;/* count for table length */
static int iiflag=0;/* status of table making */
static short *Aijtab=0;/* unroll table */
static short *pAijtab=0;/* address in unroll table */
static int qdimsave=0;/* previous qdim */
static int pdimsave=0;/* previous pdim */

static long t0=0,tf=0;
long time();

static double cvfact=1e-6;
	/* convergence factor, |(d(new)-d(previous))/d(new)|<cvfact */

int initflag;
	/* flag to select initial guess, 0=identity,1=input matrix, 2=first
	projection of sigma_in, 3= John's average + first projecton */

int anderson=0;/* 0=burg,1=anderson, no anderson for complex case */

int cntAij=0;/* count of Aij construction */

genburg(sigma_in,isigma_in,qd,pdist,sigma_out,isigma_out,monitor_flag,returnflag,R_out,iR_out)
double *sigma_in,*isigma_in,*pdist,*sigma_out,*isigma_out;
double *R_out,*iR_out;
int *qd,*returnflag,monitor_flag;
{
	/* assign input values to internal variables */

	cntAij=0;
	S=sigma_in;
	iS=isigma_in;
	qdim= *qd;
	moflag = monitor_flag;
	set_up();/* calculate constants */

	if(do_alloc())/* allocate storage */
		{
		fprintf(stderr,"Storage allocation not possible for genburg().\n");
		gb_free();
		*returnflag=6;
		return;
		}

	doPtab();/* make tables for making the bases matrices */

	ckS();/* get ln|S| */

	init(initflag,sigma_out,isigma_out,R,iR);/* select initial R vector and put in R[] and R0[] */

	if(moflag)
		{
		fprintf(stderr,"Rinit=\n");
		dvprintf(R,Ni); if(cflag)dvprintf(iR,Ni);
		fprintf(stderr,"ln|Rinit|=%g\n",R_det);
		fprintf(stderr,"d(S,Rinit)=%.15g\n",R_measure);
		}

	*returnflag=Iterate();/* main iteration loop */

	makesigma(R,iR,sigma_out,isigma_out,Pi);/* construct output matrix
	from final R vector */

	dvmove(R,R_out,Ni);
	if(cflag)
		dvmove(iR,iR_out,Ni);

	*pdist=R_measure;/* output final measure */

	gb_free();/* free storage */
	return;
}
static
init(iflag,sigma,isigma,R,iR)
/* select initial R vector */
/* check for positive definiteness, revert to a positive definite initial R */
/* setup R and R0 */
register int iflag;
register double *sigma,*isigma,R[],iR[];
{



	if(iflag==1)
		{
		exR(sigma,isigma,R,iR,qp,Ni); /* extract R from sigma */
		ckRinit();
		}

	else if(iflag==2)
		{
		doinit(S,iS,R,iR);/* use S projection as Rinit */
		ckRinit();
		}

	else if(iflag==3 && rflag)
		{
		dojohnavg(S,siginv,qp);/* average matrix according to John */
		doinit(siginv,isiginv,R,iR);/* use S projection as Rinit */
		ckRinit();
		}

	else /* set initial R to the identity matrix */
		{
		initflag=0;
		dvclear(R,Ni);
		if(cflag)
			dvclear(iR,Ni);
		R[0]=1;
		doidentity(siginv,isiginv,qp);/* set inverse and measure */
		}

	dvmove(R,R0,Ni);
	dvclear(Rnew,Ni);
	if(cflag)
		{
		dvmove(iR,iR0,Ni);
		dvclear(iRnew,Ni);
		}
	return;
}
/* init()
	Based upon the input flag 'initflag', sets up R[] with the initial
	guess, if initflag==0, the initial guess is the identity matrix,
	if initflag==1, the initial guess is extracted via exR() from
	the sigma given in the input, if initflag==2, the initial guess
	is the projection R[i]=trace[S,Pi]/Pii[i], obtained via doinit()
	where S is the input sample covariance matrix, Pi are the basis
	matrices, and Pii[i] is the metric of the basis matrices, namely,
	Pii[i]=trace[Pi,Pi], the off diagonal elements of the metric
	are zero since the Pi are chosen to be orthogonal.
	The initial R vector is checked for positive definitness and
	reversion to a positive definite initial R is done if necessary.
*/
doidentity(siginv,isiginv,qp)/* set inverse and measure for identity start */
register double *siginv,*isiginv;
register int qp;
{
	register int i;
	double trace();

	dmclear(siginv,qp);/* set inverse to identity matrix */
	for(i=0;i<qp;i++,siginv+=qp+1)
		*siginv=1;
	if(cflag)
		dmclear(isiginv,qp);
	R_measure=trace(S,qp)-qp;
	R_measure-=logdetS;/* measure for identity start */
	return;
}
doinit(S,iS,R,iR)/* make an R vector from the input sample covariance matrix */
register double S[],iS[],R[],iR[];
{
	register int fNi,i;
	register double Rt,temp;
	double spdmtrace();

	fNi=Ni;
	for(i=0;i<fNi;i++)
		{
		make(Pi,i);
		Rt=spdmtrace(S,Pi,qp);/* trace of M1*M2, M1 double, M2 short */
		temp=Rt/Pii[i];
		R[i]=temp;
		}
	if(cflag)
		{
		iR[0]=0;
		for(i=1;i<fNi;i++)
			{
			make(Pi,i);
			/* make the symmetric matrix Pi antisymmetric */
			antisymmetrize(Pi,qp);
			Rt=spdmtrace(iS,Pi,qp);
			temp= -Rt/Pii[i];
			iR[i]=temp;
			}
		}
	return;
}
/* doinit()
		Forms the R[] vector from the input sample covariance
		matrix, R[i]=trace[S,Pi]/Pii[i], where S is the sample
		covariance matrix, Pi are the basis matrices, and Pii[i]
		is trace[Pi,Pi]. This corresponds to the initial guess
		(almost) of Anderson.
*/
ckS()/* check S for positive definitness and get ln|S| */
{
	register int rtflag;

	if(rflag)
		{
		dmmove(S,sigtemp,qp);

		/* test for positive definitness with cholesky and get ln|S| */
		rtflag=lsyminv(sigtemp,siginv,sig1temp,qp,&logdetS);

		if(rtflag==1)
			{
			fprintf(stderr,"sigma sample singular (or possibly negative definite)\n");
			logdetS=0;
			}
		if(rtflag==2) 
			fprintf(stderr,"sigma sample non-positive-definite\n");
		}

	if(cflag)
		{
		fprintf(stderr,"no ln|S| available, set to zero\n");
		logdetS=0;
		/*
		rtflag=lhdmsym(siginv,isiginv,siginv,isiginv,&logdetS,qp);
		*/
		}

	if(moflag)
		fprintf(stderr,"ln|S|=%g\n",logdetS);

	return(rtflag);
}
/* ckS()
		Check to see if the sample covariance matrix S is
		positive definite, singular or negative definite,
		by computing its inverse via a routine that also
		returns the determinant.
		The determination of the properties of S is not
		definitive.
*/

Iterate()/* main interation loop */
{
	register int rtflag;

	cnt=0;	/* count of number of times no increase in measure */
	it_cnt=0;	/* iteration count */
	while(cnt<Maxit)
		{
		it_cnt++;

		/* with initial unit matrix, compute Rnew quickly */
		if(it_cnt==1 && initflag==0)
			doinit(S,iS,Rnew,iRnew);/* use S projection as Rnew */

		else
			{
			t0=time(0);
			rtflag=do_algorithm();/* compute a new Rnew using R0 */
			tf=time(0);
/*
			fprintf(stderr,"talgo=%ld\n",tf-t0);
*/
			if(rtflag==3)
				{
				fprintf(stderr,"Rnew(1) not possible, Aij singular\n");
				fprintf(stderr,"Quiting.\n");
				break;
				}

			}

		rtflag=ckRnew();/* check Rnew for positive definitness */

		if(rtflag==0)
			continue;
		if(rtflag==1)
			break;/* Rnew singular */
		if(rtflag==4)
			break;/* 1/10 interpolation failed */
		if(rtflag==5)
			break;/* converged by cvtest */
		}
	if(moflag)
		{
		if(cnt==Maxit) fprintf(stderr,"\nMaxit reached, stop\n");
		fprintf(stderr,"after %d iterations\n\n",it_cnt);
		fprintf(stderr,"R final =\n");
		dvprintf(R,Ni); if(cflag) dvprintf(iR,Ni);
		fprintf(stderr,"ln|Rfinal|=%g\n",R_det);
		fprintf(stderr,"d(S,Rfinal)=%.15g\n",R_measure);
		if (mcflag)fprintf(stderr,"coherence=%g\n",R[2]*R[2]/(R[0]*R[1]));
		}
	return(rtflag);
}
/* Iterate()
		Main loop of the program. Checks to see if 'Maxit' iterations
		have occured during which no new solution has been found
		with a reduced measure.
		Via do_algorithm(), computes a new Rnew[] vector and
		checks to see if it is positive definite via ckRnew().
		ckRnew() also calls an interpolation routine interp()
		if the new solution is not positive definite.
*/
ckRinit()/* check Rinit or first projection of S for positive definitness */
{
	register int rtflag;
	double domeasure();

	/* get inverse and ln|R| to do initial measure */

	rtflag=inv_det(R,iR,siginv,isiginv,&R_det);

	if(rtflag)
		{
		if(rtflag==1) fprintf(stderr,"initial R singular or negative definite\n");
		if(rtflag==2) fprintf(stderr,"initial R non-positive-definite\n");
		if(initflag==1)/* sigma_out start */
			{
			if(moflag) fprintf(stderr,"Reverting to first projection for initial quess\n");
			initflag=2;/* first projection */
			init(initflag,siginv,isiginv,R,iR);
			}
		else if(initflag==2 || initflag==3)/* first projection start */
			{
			if(moflag) fprintf(stderr,"Reverting to identity matrix for initial quess\n");
			initflag=0;
			init(initflag,siginv,isiginv,R,iR);
			rtflag=inv_det(R,iR,siginv,isiginv,&R_det);
			}

		}
	R_measure=domeasure(siginv,isiginv,R_det);
	return;
}
/* ckRinit()
		Check the initial guess for R[] to see if it is positive
		definit and in the process, get the measure for the
		initial guess. If the initial guess is not positive
		definite and is given as input, revert to the 
		first projection of S as the initial guess. If that
		is not positive definite, revert to the identity matrix
		as the initial guess.
*/
ckRnew()	/* check Rnew for positive definitness */
{
	register int rtflag,rttflag;
	double domeasure();
	double Rtemp_measure;/* R_measure saved for interpolation then cvtest */

	rtflag=inv_det(Rnew,iRnew,siginv,isiginv,&Rnew_det);
		 /* 0 Rnew pos-def
		    1 Rnew singular, Rnew_det=0
		    2 Rnew non-pos-def, Rnew_det=0 */
	Rnew_measure=domeasure(siginv,isiginv,Rnew_det);

	if(rtflag==1)
		{
		if(moflag)
			{
			fprintf(stderr,"\nRnew singular\n");
			fprintf(stderr,"Rnew=\n");
			dvprintf(Rnew,Ni);if(cflag)dvprintf(iRnew,Ni);
			fprintf(stderr,"interpolation not possible\n");
			}
		return(rtflag);/* break */
		}
	if(rtflag==2)
		{
		if(moflag) fprintf(stderr,"\nRnew npd\n");
		Rtemp_measure=R_measure;/* save measure for cvtest */
		rttflag=interp();/* do interpolation by 1/10 */
		if(rttflag==1)
			{
			if(moflag) fprintf(stderr,"stop, complete interpolation\n");
			return(4);
			}
		/* interp found a positive definite Rnew */
		if(cvtest(R_measure,Rtemp_measure))return(5);/* conclude iteration */
		return(0);
		}
	if(rtflag==0)
		{
		if(moflag)
			{
			fprintf(stderr,"\nRnew(1)=\n");
			dvprintf(Rnew,Ni); if(cflag) dvprintf(iRnew,Ni);
			fprintf(stderr,"ln|Rnew|=%g\n",Rnew_det);
			fprintf(stderr,"d(S,Rnew)=%.15g\n",Rnew_measure);
			}
		if(Rnew_measure>=R_measure)
			{
			Rtemp_measure=R_measure;/* save measure for cvtest */
			if(moflag)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;

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