gsm610.c

GSM6.10 算法优化及详细注释 本代码较原GSM6.10代码更容易理解 仅供学习

//-----------------------------------------------------------------------------
//	Copyright(c) 2005-2006 Xiamen Yealink Network Technology Co,.Ltd
//	project:sip-phone 
//	By:陈剑辉		2005/1
// 	Discription: 语音编码，采用gsm6.10编码
//-----------------------------------------------------------------------------

#include "gsm610.h"
#include "operate.h"
/************************************************************************************************/
/*************************Offset compensation,Pre-emphasis(开始)*********************************/

/* 
 * 输入：s_in[0..159]未处理最原始的数据
 * 输出：s_out[0..159]经过偏移与缩减后的数据
 * s_in为未压缩语音数据;s_in,s_out长度为160;缩减预处理原始数据
 * 运算过程sof(k) = so(k) - so(k-1) + 32735*2^(-15)*sof(k-1)      s(k) = sof(k) - 28180*2^(-15)*sof(k-1)
 */
void Gsm_Preprocess (struct gsm610_state * g_gsm, word * s_in, word * s_out )
{

	word			pre_so = g_gsm->pre_so;
	longword		l_pre_sof = g_gsm->l_pre_sof;
	word			pre_sof = (word) (g_gsm->pre_sof);
	word			temp;
	longword		l_temp;
	word			msp;		/* msp为l_pre_sof*2^-15 */
	word			lsp;		/* lsp存储l_pre_sof低15位 */
	word			so;
	register int	k = 160;
	longword		ltmp;			/* for	 ADD */
	ulongword		utmp;			/* for L_ADD */

	
	while (k--)
	{
		so = SASR( *s_in, 3 ) << 2;
		assert (so >= -0x4000); /* downscaled by	 */
		assert (so <=  0x3FFC); /* previous routine. */

		temp = so-pre_so;       /* so(k)-so(k-1) */
		assert (temp != MIN_WORD);
		
		l_temp = temp;
		l_temp <<= 15; /* l_temp=(so(k)-so(k-1))*2^15 */
		
		msp = (word) (SASR( l_pre_sof, 15 )); /* lsp为l_z2低15位,msp为l_z2 * 2^(-15) l_z2为Sof(k-1) */
		lsp = (word) (l_pre_sof-((longword)msp<<15)); /* gsm_L_sub(L_z2,(msp<<15)); */

		l_temp += GSM_MULT_R( lsp,32735);
		l_temp = GSM_L_ADD( (longword)msp * 32735, l_temp);
		l_pre_sof = l_temp;
		l_temp = GSM_L_ADD( l_temp, 16384 );     /* 四舍五入 */
		temp   = (word) (GSM_MULT_R( pre_sof, -28180 )); /* Sof(k-1)*2^-15*(-28180) */
		pre_sof	  = (word) (SASR( l_temp, 15 ));  /* mp=Sof(k)四舍五入后的Sof(k) */
		*s_out = (word) (GSM_ADD( pre_sof, temp ));
		
		pre_so = so;
		s_in++;
		s_out++;
	}

	g_gsm->pre_so	=  pre_so;      /* 上一个缩减的数据SO(k-1) */
	g_gsm->l_pre_sof = l_pre_sof;   /* Sof(k-1)*2^15 */
	g_gsm->pre_sof	= pre_sof;      /* Sof(k-1) */
}

/*
 * 输入：s_out[0..159]未还原的数据
 * 输出：s_out[0..159]扩展后的数据，也是完成一帧160个解码后的采样数据
 * 运算过程s_out(k) = s(k) + 28180*2^(-15)*s_out(k-1)
 */
static void Postprocessing (struct gsm610_state * g_gsm, register word * s_out)
{
	register int		k;
	register word		pre_s_out = g_gsm->pre_s_out;
	register longword	ltmp;
	register word		temp;
	
	for (k = 0; k<160; k++)
	{
		temp = (word) GSM_MULT_R( pre_s_out, 28180);
		pre_s_out = (word)GSM_ADD(*s_out, temp);
		*s_out = (word)(GSM_ADD(pre_s_out, pre_s_out) & 0xFFF8);
		g_gsm->pre_s_out = pre_s_out;
		s_out ++;
	}
}

/*************************Offset compensation,Pre-emphasis(结束)*********************************/
/************************************************************************************************/
/**************************LPC自相关分析(开始)***************************************************/

/*
 * 根据160个采样数据s_in,计算L_ACF[0..8]
 */
static void Autocorrelation (
		word	 * s_in,		/* [0..159] 	IN/OUT	*/
		longword * L_ACF)		/* [0..8]		OUT 	*/
{
	register int	k, i;
	word			temp;
	word			smax;
	word			scalauto;
	word			*sp;
	word			sl;
	
	smax = 0;
	for (k = 0; k < 160; k++) {
			temp = GSM_ABS( s_in[k] );
			if (temp > smax) smax = temp;
	}

	/*	Computation of the scaling factor.
	 */
	if (smax == 0) 
		scalauto = 0;
	else 
	{
			assert(smax > 0);
			scalauto = 4 - gsm_norm( (longword)smax << 16 );/* sub(4,..) */
	}

	/*	Scaling of the array s[0...159]
	 */

	if (scalauto > 0)
	{
		for (k=0; k<160; k++)
			s_in[k] = GSM_MULT_R( s_in[k], 16384 >> (scalauto-1)); /* 右移scalauto位保证计算结果不溢出 */
	}

	sp = s_in;
	sl = *sp;
	
	for ( k=0; k<9; k++)
	{
		L_ACF[k]=0;
		for ( i=k; i<160; i++)
		{
			sl=*(sp+i);
			L_ACF[k] += (longword)((longword) sl * sp[ i-k ]);
		}
		L_ACF[k] <<= 1;
	}

	if (scalauto > 0) 
	{
		assert(scalauto <= 4); 
		for (k = 0; k<160; k++)
		{
			*s_in <<= scalauto;
			s_in++;
		}
	}
}

/*
 * 输入L_ACF[0..9]根据schur recursion计算输出r[0..7]
 */
static void Reflection_coefficients (longword * L_ACF, register word * r)
{
	register int	i, m, n;
	register word	temp;
	register longword ltmp;
	word			ACF[9]; /* 0..8 */
	word			P[	9]; /* 0..8 */
	word			K[	9]; /* 2..8 */

	/*	Schur recursion with 16 bits arithmetic.
	 */

	if (L_ACF[0] == 0) /* everything is the same. */
	{	
		for (i = 0; i<8; i++)
			*r++ = 0 ;
		return;
	}

	assert( L_ACF[0] != 0 );
	temp = gsm_norm( L_ACF[0] );
	
	assert(temp >= 0 && temp < 32);
	for (i = 0; i <= 8; i++) 
		ACF[i] = ((word) (SASR( L_ACF[i] << temp, 16 )));

	/* Initialize array P[..] and K[..] for the recursion. */
	for (i = 1; i <= 7; i++) 
		K[ i ] = ACF[ i ];
	for (i = 0; i <= 8; i++) 
		P[ i ] = ACF[ i ];
	
	/* Compute reflection coefficients */
	for (n = 1; n <= 8; n++, r++) 
	{
		temp = P[1];
		temp = GSM_ABS(temp);
		if (P[0] < temp) 
		{
			for (i = n; i <= 8; i++) 
				*r++ = 0;
			return;
		}

		*r = gsm_div( temp, P[0] );

		assert(*r >= 0);
		if (P[1] > 0) 
			*r = -*r; 		/* r[n] = sub(0, r[n]) */
		assert (*r != MIN_WORD);
		if (n == 8)
			return; 

		/*	Schur recursion
		 */
		temp = GSM_MULT_R( P[1], *r );
		P[0] = (word) GSM_ADD( P[0], temp );

		for (m = 1; m <= 8 - n; m++) 
		{
			temp = GSM_MULT_R( K[ m ], *r );
			P[m] = (word) GSM_ADD( P[ m+1 ], temp );

			temp = GSM_MULT_R( P[ m+1 ], *r );
			K[m] = (word) GSM_ADD( K[ m ],  temp );
		}
	}
}


/*
 *	LAR(i) = r(i) ; |r(i)| < 0.675
 *	LAR(i) = sign[r(i)]*[2|r(i)|-0.675] ; 0.675 <= |r(i)| < 0.950
 *	LAR(i) = sign[r(i)]*[8|r(i)|-6.375] ; 0.950 <= |r(i)| <= 1.000
 *  输出LAR(i)*1/2;
 *  The following scaling for r[..] and LAR[..] has been used:
 *	r[..]	= integer( real_r[..]*32768. ); -1 <= real_r < 1.
 *	LAR[..] = integer( real_LAR[..] * 16384 );
 *	with -1.625 <= real_LAR <= 1.625
 */
static void Transformation_to_Log_Area_Ratios (register word * r) /* 0..7    IN/OUT */


{
	register word	temp;
	register int	i;


	/*  Computation of the LAR[0..7] from the r[0..7]
	 */
	for (i = 1; i <= 8; i++, r++)
	{
		temp = *r;
		temp = GSM_ABS(temp);
		assert(temp >= 0);

		if (temp < 22118)
		{
			temp >>= 1;
		} 
		else if (temp < 31130)
		{
			assert( temp >= 11059 );
			temp -= 11059;
		} 
		else 
		{
			assert( temp >= 26112 );
			temp -= 26112;
			temp <<= 2;
		}

		*r = *r < 0 ? -temp : temp;
		assert( *r != MIN_WORD );
	}
}


/*
 * 量化LAR[0..7]输出LARc[0..7],LARc[i] = Nint{A[i]*LAR[i] + B(i)
 */
static void Quantization_and_coding (register word * LAR) 	/* [0..7]		IN/OUT	*/
{
	register word	temp;
	longword		ltmp;
	register int	i;
	/*	This procedure needs four tables; the following equations
	 *	give the optimum scaling for the constants:
	 *	
	 *	A[0..7] = integer( real_A[0..7] * 1024 )
	 *	B[0..7] = integer( real_B[0..7] *  512 )
	 *	MAC[0..7] = maximum of the LARc[0..7]
	 *	MIC[0..7] = minimum of the LARc[0..7]
	 */
    for (i = 0; i < 8; i++)
	{
		temp = (word) GSM_MULT( gsm_A[i],   *LAR ); 
		temp = (word) GSM_ADD(  temp,  gsm_B[i] );
		temp = (word) GSM_ADD(  temp, 256 );
		temp = SASR( temp,  9 );
		if (temp > gsm_MAC[i])
			*LAR = gsm_MAC[i] - gsm_MIC[i];
		else if (temp < gsm_MIC[i])
			*LAR = 0;
		else
			*LAR = temp - gsm_MIC[i];
		LAR++;
	}
}


/*
 * 输入160个经过预处理的采样值s_in,输出自相关系数LARc[0..7]
 */
void Gsm_LPC_Analysis (struct gsm610_state *g_gsm, word	* s_in, word  * LARc)
{
	longword		L_ACF[9];

	Autocorrelation (s_in ,	L_ACF );
	Reflection_coefficients (L_ACF, LARc);
	Transformation_to_Log_Area_Ratios (LARc);
	Quantization_and_coding (LARc);
}

/**************************LPC自相关分析(结束)***************************************************/
/************************************************************************************************/
/**************************Short Term Analysis(开始)*********************************************/

/*
 * 输入量化值LARc[0..7]计算输出LARpp[0..7]
 */
static void Decoding_of_the_coded_Log_Area_Ratios (word	* LARc, word * LARpp)
{
	register word	temp;
	register long	ltmp;	/* for GSM_ADD */
	int				i;
		
	for (i = 0; i < 8; i++)
	{
		temp = (word) (GSM_ADD( (long) *LARc, (long) gsm_MIC[i] ) << 10);
		temp = (word) GSM_SUB( (long) temp, (long) gsm_B[i] << 1 ); 
		temp = (word) GSM_MULT_R( (long) gsm_INVA[i], (long) temp );           \
		*LARpp = (word) GSM_ADD( (long) temp, (long) temp );
		LARc++;
		LARpp++;
	}
}

/*
 * 结合上一帧LARpp_j_1[0..7]计算LARp[0..7],根据一帧当中不同数据位置计算如下：
 * 0 ...12 : 0.75*LAR' 'J-1(i) + 0.25*LAR' 'J(i)
 * 13...26 : 0.50*LAR' 'J-1(i) + 0.50*LAR' 'J(i)
 * 27...39 : 0.25*LAR' 'J-1(i) + 0.75*LAR' 'J(i)
 * 40..159 : LAR' 'J(i)
 */
void Coefficients(register word * LARpp_j_1,register word * LARpp_j,register word * LARp,int ntype)
{
	register int	i;
	register longword ltmp;
	for ( i = 0; i <8; i++)
	{
		switch (ntype)
		{
		case 0:
			*LARp = (word) GSM_ADD( SASR( *LARpp_j_1, 2 ), SASR( *LARpp_j, 2 ));
			*LARp = (word) GSM_ADD( *LARp,  SASR( *LARpp_j_1, 1));
			break;
		case 1:
			*LARp = (word) GSM_ADD( SASR( *LARpp_j_1, 1), SASR( *LARpp_j, 1 ));
			break;
		case 2:
			*LARp = (word) GSM_ADD( SASR( *LARpp_j_1, 2 ), SASR( *LARpp_j, 2 ));
			*LARp = (word) GSM_ADD( *LARp, SASR( *LARpp_j, 1 ));
			break;
		case 3:
			*LARp = *LARpp_j;
			break;
		}
		LARpp_j_1++;
		LARpp_j++;
		LARp++;
	}
}

/*
 * 输入估计值LARp[0..7]计算r'[0..7],
 * Transformation_to_Log_Area_Ratios的逆运算
 */
static void LARp_to_rp (register word * LARp)
{
	register int			i;
	register word			temp;
	register longword		ltmp;
	for (i=0 ; i < 8; i++)
	{
		temp = GSM_ABS( *LARp );
		if (temp < 11059)
			temp <<=1;
		else if (temp < 20070)
			temp += 11059;
		else
			temp = (word) GSM_ADD( temp >> 2, 26112 );
		if (*LARp < 0)
			temp = -temp;
		*LARp = temp;
		LARp++;
	}
}

/* 
 * 输入：rp[0..7],k_n 要计算的数据长度k_end-k_start
 * 输出: s[0..k_n] 根据自相关系数分析后的值
 */
static void Short_term_analysis_filtering (struct gsm610_state * g_gsm, register word * rp, int k_n, register word * s)
{
	register word			* u = g_gsm->u;
	register int			i;
	register word			di, zzz, ui, sav, rpi;
	register longword		ltmp;

	for (; k_n--; s++)
	{
		di = sav = *s;

		for (i = 0; i < 8; i++) 
		{			
			ui	  = u[i];
			rpi   = rp[i];
			u[i]  = sav;

			zzz   = (word) GSM_MULT_R(rpi, di);
			sav   = (word) GSM_ADD(   ui,	zzz);

			zzz   = (word) GSM_MULT_R(rpi, ui);
			di	  = (word) GSM_ADD(   di,	zzz );
		}

		*s = di;
	}
}

/* 
 * 输入：rp[0..7],k_n 要计算的数据长度k_end-k_start
 * 输出: s[0..k_n-1] 根据自相关系数合成后的值
 * Short_term_analysis_filtering()的逆过程
 */
static void Short_term_synthesis_filtering (struct gsm610_state * g_gsm, register word * rrp, register int k, register word * wt, register word * s_out)
{
	register word			* v = g_gsm->v;
	register int			i;
	register word			sri, tmp1, tmp2;
	register longword		ltmp;	/* for GSM_ADD	& GSM_SUB */

	while (k--) 
	{
		sri = *wt++;
		for (i = 8; i--;) 
		{

			/* sri = GSM_SUB( sri, gsm_mult_r( rrp[i], v[i] ) );
			 */
			tmp1 = rrp[i];
			tmp2 = v[i];
			tmp2 =	( tmp1 == MIN_WORD && tmp2 == MIN_WORD
					? MAX_WORD