ikdistn.c

来自「speech signal process tools」· C语言 代码 · 共 105 行

/*	

 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 SCCS
static 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 distortion
measure 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];
}