cnormal.c

speech signal process tools

cnormal(R,iR,a,ia,p,n)
	/* rtflag=1 singular, rtflag=2 non-pos-def */
double R[],iR[],a[],ia[],p[];
register int n;
{
	register int i,j,ss,s,iM;
	int rtflag;
	register double ci,ici,po,pf,Rt,iRt,at,iat,temp;

	rtflag=0;

	for(i=0,iM=0;i<n;i++,iM+=n)
		{
		a[iM]=1;
		ia[iM]=0;
		for(j=i+1;j<n;j++)
			a[iM+j]=ia[iM+j]=0;
		}
	p[0]=R[0];
	for(i=1,s=n;i<n;i++,s+=n)
		{
		ss=s-n;
		ci=R[i];
		ici=iR[i];
		for(j=1;j<i;j++)
			{
			Rt=R[i-j];
			iRt=iR[i-j];
			at=a[ss+j];
			iat=ia[ss+j];
			ci+= at*Rt-iat*iRt;
			ici+= at*iRt+iat*Rt;
			}
		if((po=p[i-1])==0)return(1);
		ci= -ci/po;
		ici= -ici/po;
		a[s+i]=ci;
		ia[s+i]=ici;
		for(j=1;j<i;j++)
			{
			temp=a[ss+j]+ci*a[ss+i-j]+ici*ia[ss+i-j];
			a[s+j]=temp;
			temp=ia[ss+j]-ci*ia[ss+i-j]+ici*a[ss+i-j];
			ia[s+j]=temp;
			}
		pf=(1-(ci*ci+ici*ici))*po;
		p[i]=pf;

		if(pf<0)rtflag++;

		if(pf==0)return(1);

		}
	if(rtflag>0)return(2);
	return(0);
}