lctoepinv.c

speech signal process tools

lctoepinv(R,iR,Rinv,iRinv,pdet,n)
	/* natural log of determinant is returned */
	/* R a vector of length n, Rinv a matrix of rank n */
register double R[],iR[],Rinv[],iRinv[],*pdet;
register int n;
{
	register int l,i,iM,j,jM,ss;

	int N,rtflag,toggle;
	static double *a,*ia,*p,Rt,iRt;
	register double det,pt,dtemp;
	double log();

	N=n-1;

	if(allocdp(&a,n*n))
		{
		printf("ctoepinv not possible due to storage limitations\n");
		exit(1);
		}
	if(allocdp(&ia,n*n))
		{
		freedp(&a);
		printf("ctoepinv not pollsible due to storage limitations\n");
		exit(1);
		}
	if(allocdp(&p,n))
		{
		freedp(&a);
		freedp(&ia);
		printf("ctoepinv not possible due to storage limitations\n");
		exit(1);
		}
	
	rtflag=cnormal(R,iR,a,ia,p,n);
	if(rtflag==1)
		{
		freedp(&a);
		freedp(ia);
		freedp(&p);
		*pdet=0;
		return(1);
		}
	
	/* if(rtflag==2) printf("ctoep non-pos-definite\n"); */

	det=0;
	if(rtflag==0)
		{
		for(i=0;i<n;i++)
			det+=log(p[i]);
		}
	if(rtflag==2)
		{
		toggle=1;
		for(i=0;i<n;i++)
			{
			if((pt=p[i])<0)
				{
				toggle= -toggle;
				pt= -pt;
				}
			det+=log(pt);
			}
		if(toggle!=1)
			det=0;
		}
	*pdet=det;

	for(i=0,iM=0;i<n;i++,iM+=n)
		{
		for(j=i,jM=i*n;j<n;j++,jM+=n)
			{
			Rt=0;
			iRt=0;
			for(l=0,ss=N*n;l<=i;l++,ss-=n)
				toepmult(a[ss+i-l],-ia[ss+i-l],a[ss+j-l],ia[ss+j-l],&Rt,&iRt,p[N-l]);
			Rinv[jM+i]=Rinv[iM+j]=Rt;
			iRinv[iM+j]=iRt;
			if(i!=j)
				{
				dtemp= -iRt;
				iRinv[jM+i]=dtemp;
				}
			}
		}
	
	freedp(&a);
	freedp(&ia);
	freedp(&p);
	return(rtflag);
}