orthog.c

speech signal process tools

orthog(M,T,n)
/* routine returns in M an orthogonal matrix O
and in T an upper triangular matrix T such that M=O*T */
/* if M non-singular, T non-singular */
register double *M,*T;
register int n;
{
	register int i,j,k;
	register double *s,*c,*sT,*s1T,*sM;
	double sqrt(),dot();

	for(k=0,s1T=T,c=M;k<n;k++,s1T+=n+1,c++)
		{
		if(k>0)
			{
			for(j=0,sT=T+k,sM=M;j<k;j++,sT+=n,sM++)
				{
				*sT=dot(sM,n,c,n,n);
				sp1(*sT,sM,c,n);
				}
			}
		(*s1T)=dot(c,n,c,n,n);
		if((*s1T)<=0)
			return(1);

		*s1T=sqrt(*s1T);
		for(i=0,s=c;i<n;i++,s+=n)
			*s/= *s1T;
		}
	return(0);
}
sp1(p1,p2,p3,n)
register double p1,*p2,*p3;
register int n;
{
	register int i;

	for(i=0;i<n;i++,p2+=n,p3+=n)
		*p3 -= p1 * *p2;
	return;
}