ltoepinv.c

speech signal process tools

ltoepinv(R,Rinv,pdet,n)
	/* natural log of determinant is returned */
double R[],Rinv[],*pdet;
register int n;
{
	register int l,i,iM,j,jM,ss;
	int N;
	int rtflag,toggle;
	double det,Rt,pt,log();
	static double *a,*p;

	N=n-1;

	if(allocdp(&a,n*n))
		{
		printf("ltoepinv not possible due to storage limitations\n");
		exit(1);
		}
	if(allocdp(&p,n))
		{
		freedp(&a);
		printf("ltoepinv not possible due to storage limitations\n");
		exit(1);
		}

	rtflag=normal(R,a,p,n);
	/* rtflag=1 singular, rtflag=2 non-pos-def */

	if(rtflag==1)
		{
		freedp(&a);
		freedp(&p);
		*pdet=0;
		return(1);
		}

	/* if(rtflag==2) printf("non-minimum phase\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;
			for(l=0,ss=N*n;l<=i;l++,ss-=n)
				Rt+= a[ss+i-l]*a[ss+j-l]/p[N-l];
			Rinv[jM+i]=Rinv[iM+j]=Rt;
			}
		}

	freedp(&a);
	freedp(&p);
	return(rtflag);
}