genburg.c

speech signal process tools

	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,
	and siginv*S*siginv=M is hermitian, so
	trace[M*Pi]=trace[Pi*M]=trace[(M*Pi)']=trace[M*Pi]complex cojugate,Q.E.D

	Bi=trace[M*Pi], M=siginv*S*siginv=Ms+iMa
	Ms is real and symmetric, iMa is imaginary and anti-symmetric
	Bi=trace[Ms*Pi] i=0 to Ni-1
	Bi=trace[iMa*iPi] i=Ni to tNi-1
	iPi is the anti-symmetrized Pi, with an i added
	*/

	/* make S*siginv in sigtemp */
	cdmmult(S,iS,siginv,isiginv,sigtemp,isigtemp,sig1temp,isig1temp,qp);

	/* make siginv*S*siginv in sig1temp */
	cdmmult(siginv,isiginv,sigtemp,isigtemp,sig1temp,isig1temp,sig2temp,isig2temp,qp);

	for(i=0;i<fNi;i++)
		{
		make(Pi,i);
		Bt=spdmtrace(sig1temp,Pi,qp);
		Bi[i]=Bt;
		}
	for(i=fNi;i<ftNi;i++)
		{
		make(Pi,i-fNi+1);
		antisymmetrize(Pi,qp);
		Bt=spdmtrace(isig1temp,Pi,qp);
		Bi[i]= -Bt;/* -1 because of product of two i's */
		}
	return;
}
antisymmetrize(Pi,n)/* make the symmetric matric Pi antisymmetric, ones in upper triangle */
register short *Pi;
register int n;
{
	register int i,j,iM,jM;

	for(i=iM=0;i<n;i++,iM+=n)
		{
		for(j=i+1,jM=iM+n;j<n;j++,jM+=n)
			{
			if(*(Pi+iM+j))
				*(Pi+jM+i)= -1;
			}
		}
	return;
}
doRnew(Ainv,Bi)/* make Ainv*Bi */
register double Ainv[],Bi[];
{
	register int i,j,iM,fNi,ftNi;
	register double Rt,*s,iRt;
	double dot();

	fNi=Ni;
	ftNi=tNi;
	for(i=0,s=Ainv;i<fNi;i++,s+=ftNi)
		{
		Rt=dot(s,1,Bi,1,ftNi);
		if(anderson && (initflag==2 || initflag==3 || it_cnt>1) && rflag)
			Rt+=R0[i];
		Rnew[i]=Rt;
		}
	if(cflag)
		{
		iRnew[0]=0;
		for(i=fNi,iM=fNi*ftNi;i<ftNi;i++,iM+=ftNi)
			{
			iRt=0;
			for(j=0;j<tNi;j++)
				iRt+=Ainv[iM+j]*Bi[j];
			iRnew[i-Ni+1]=iRt;
			}
		}
	return;
}
makesc(Pi,s,qp)/* make single channel basis matrix */
register short Pi[];
register int s,qp;
{
	register int i,del,lim;
	register short *p,*pm,value;

	lim=qpxqp;
	value=0;
	for(i=0,p=Pi;i<lim;i++,p++)
		*p=value;
	lim=qp-s;
	value=1;
	del=qp+1;
	for(i=0,p=Pi+s,pm=Pi+s*qp;i<lim;i++,p+=del,pm+=del)
		{
		*p=value;
		*pm=value;
		}
	/*
	fprintf(stderr,"\n");
	for(i=0;i<qp;i++)
		{
		for(j=0;j<qp;j++)
			fprintf(stderr,"%d ",*(Pi+i*qp+j));
		fprintf(stderr,"\n");
		}
	*/
	return;
}
make(Pi,s)
register short Pi[];
register int s;
{
	int i,j,k,l,P,a,b,c,kp,lp;
	register int m,n;

	if(scflag)
		{
		makesc(Pi,s,qp);
		return;
		}
	if(mdflag)
		{
		a=Ptaba[s]; b=Ptabb[s];
		for(k=0,kp=0;k<qdim;k++,kp+=pdim)
			{
			for(l=k,lp=kp;l<qdim;l++,lp+=pdim)
				{
				for(i=0;i<pdim;i++)
					{
					j=0;
					if(l==k)j=i;
					for(;j<pdim;j++)
						{
						m=kp+i;
						n=lp+j;
						P=0;
						if( (a == (l-k)) && (b == (j-i))) P=1;
						Pi[m*qp+n]=Pi[n*qp+m]=P;
						}
					}
				}
			}
		return;
		}
	if(mcflag)
		{
		a=Ptaba[s];
		b=Ptabb[s];
		c=Ptabc[s];
		for(k=0,kp=0;k<qdim;k++,kp+=pdim)
			{
			for(l=k,lp=kp;l<qdim;l++,lp+=pdim)
				{
				for(i=0;i<pdim;i++)
					{
					j=0;
					if(l==k)j=i;
					for(;j<pdim;j++)
						{
						m=kp+i;
						n=lp+j;
						P=0;
						if((a == (l-k)) && (b == (j-i)))
							{
							if(b==0 && i==c)P=1;
							if(b!=0) P=1;
							}
						Pi[m*qp+n]=Pi[n*qp+m]=P;
						}
					}
				}
			}
		return;
		}
	fprintf(stderr,"trouble in make().\n");
	exit(1);
}
doPtab()
{
	int i,a,b,j,isum,k,l;

	if(scflag || mdflag)
		{
		i=0;
		a=0;
		for(b=0;b<pdim;b++)
			{
			Ptaba[i]=a;
			Ptabb[i]=b;
			i++;
			}
		for(a=1;a<qdim;a++)
			{
			for(b= -(pdim-1);b<pdim;b++)
				{
				Ptaba[i]=a;
				Ptabb[i]=b;
				i++;
				}
			}
		}
	if(mcflag)
		{
		i=0;
		a=0;
		b=0;
		for(j=0;j<pdim;j++)
			{
			Ptaba[i]=a;
			Ptabb[i]=b;
			Ptabc[i]=j;
			i++;
			}
		for(b=1;b<pdim;b++)
			{
			Ptaba[i]=a;
			Ptabb[i]=b;
			i++;
			}
		for(a=1;a<qdim;a++)
			{
			for(b= -(pdim-1);b<pdim;b++)
				{
				if(b==0)
					{
					for(j=0;j<pdim;j++)
						{
						Ptaba[i]=a;
						Ptabb[i]=b;
						Ptabc[i]=j;
						i++;
						}
					}
				else
					{
					Ptaba[i]=a;
					Ptabb[i]=b;
					i++;