/*	 This material contains proprietary software of Entropic Processing, Inc. Any reproduction, distribution, or publication without the the prior	    written permission of Entropic Processing, Inc. is strictly prohibited. Any public distribution of copies of this work authorized in writing by Entropic Processing, Inc. must bear the notice			       	Copyright 1986, Entropic Proccessing, Inc	(C) 1985, Entropic Processing, Inc.         */#ifdef SCCSstatic char *sccsid = "@(#)ikdistn.c	1.2 4/28/86";#endif	#include <stdio.h>extern int debug_level;IKdistn(dPrc,dQrc,order,pdistn, pgain)/*This routine computes Itakura-Saito distortionmeasure between two all-pole spectra*/double dPrc[],dQrc[];float *pgain,*pdistn;int order;{	int i,j;	float distn,t,gain;	double Qr[1000];	float Pa[1000],Prc[1000],Qrc[1000];	double log();	/* double to single precision conversion to make my routines 	compatible with wilson algorithm output. */	for (i = 0; i < 1; i++)	{		Qrc[i] = dQrc[i];		Prc[i] = dPrc[i];	}	for (i = 1; i <=  order; i++)	{		Qrc[i] = -dQrc[i];		Prc[i] = -dPrc[i];	}	rctoc(Prc,order,Pa,&gain);	rctor(Qrc,order,Qr);	Pa[0] = -1.0;	t = 0; 	for (i = 0; i <= order; i++)		for (j = 0; j <= i; j++)			if (i == j)				t = t + Pa[i] * Pa[j] * Qr[i-j];			else	t = t + 2.0 * Pa[i] * Pa[j] * Qr[i-j];	distn = t / Prc[0] - log((double) Qrc[0]/Prc[0]) -1;	gain = Prc[0]/t;	if(debug_level>=2)		fprintf(stderr,"D(Q,P) = %f   opt gain = %f   ",distn,gain); 	distn = -log((double) gain) -log((double) Qrc[0]/Prc[0]);	if(debug_level>=2)		fprintf(stderr,"D(kQ,P) = %f\n",distn);	*pdistn = distn; *pgain = gain;}float Newton_Raphson (y)float y;{	float t,xn;	int i;	if (y < .001) return(1.0);	xn = 1.0 + 2.0*y;	for (i = 0; i < 5; i++)	{		t = xn * (y + log((double) xn))/(xn - 1.0);		xn = t;	}	return(xn);}rctor(rc,n,r)int n;float rc[];double r[];{	float state[1000],nrc[1000],a[1000],num[1000];	float fowrd_filter(),gain;	int i;	rctoc (rc,n,a,&gain);	r[0] =  rc[0]/gain;	for (i = 1; i <= n; i++)	{		state[i] = 0.0;		nrc[i] = -rc[i];	}	rctoc (nrc,n,num,&gain);	num [0] = r[0]/2.0;	for (i = 1; i <= n; i++) num[i] = -0.5 * r[0] * num[i];	for (i = 0; i <= n; i++)		r[i] =  fowrd_filter(a,state,n,num[i]);	r[0] = 2.0 * r[0];}