npp.c

MELPe 1200 bps, fixed point

		/*	Gain[i] = ksi_vq * exp (0.5 * eiv); */

		L_sum = L_mpy_ls(L_sum, 23637);                  /* 0.72135 to base 2 */
		shift = shr(extract_h(L_sum), 8);                        /* get shift */
		if (comp_data_shift(ksi_vq, shift, 32767, 0) > 0){
			Gain[i] = 32767;
			continue;
		} else {
			temp1 = extract_l(L_shr(L_sum, 9));
			temp1 = (Shortword) (temp1 & 0x7fff);
			temp1 = shr(mult(temp1, 9864), 3);       /* change to base 10 Q12 */
			temp1 = pow10_fxp(temp1, 14);

			L_sum = L_shl(L_deposit_h(ksi_vq), shift);
			L_sum = L_mpy_ls(L_sum, temp1);
			if (L_sub(L_sum, 1073676288L) > 0)
				Gain[i] = 32767;
			else
				Gain[i] = extract_h(L_shl(L_sum, 1));
		}
	}
}


/*****************************************************************************/
/*	   Subroutine ksi_min_adapt: computes the adaptive ksi_min				 */
/*****************************************************************************/
static Shortword	ksi_min_adapt(BOOLEAN n_flag, Shortword ksi_min,
								  Shortword sn_lt, Shortword sn_lt_shift)
{
	Shortword	ksi_min_new, temp, shift;
	Longword	L_sum;


	if (n_flag)             /* RM: adaptive modification of ksi limit (10/98) */
		ksi_min_new = ksi_min;
	else {
		if (sn_lt_shift > 0){
			L_sum = L_add(L_deposit_l(sn_lt), L_shr(X05_Q15, sn_lt_shift));
			shift = sn_lt_shift;
		} else {
			L_sum = L_add(L_shl(L_deposit_l(sn_lt), sn_lt_shift), X05_Q15);
			shift = 0;
		}

		/* L_sum is (sn_lt * 2^sn_lt_shift + 0.5) */
		if (L_sum > ONE_Q15){
			L_sum = L_shr(L_sum, 1);
			shift = add(shift, 1);
		}
		temp = extract_l(L_sum);

		/* We want to compute 2^{0.65*[log2(0.5+sn_lt)]-7.2134752} */
		/* equiv. 2^{0.65*log10(temp)/log10(2) + 0.65*shift - 7.2134752} */
		temp = log10_fxp(temp, 15);                 /* log(0.5 + sn_lt), Q12 */
		L_sum = L_shr(L_mult(temp, 8844), 9);       /* Q12*Q12->Q25 --> Q16 */

		L_sum = L_add(L_sum, L_mult(shift, 21299)); /* 21299 = 0.65 Q15 */
		L_sum = L_sub(L_sum, 472742L);              /* 472742 = 7.2134752 Q16 */
		shift = extract_h(L_sum);                   /* the pure shift */
		temp = (Shortword) (extract_l(L_shr(L_sum, 1)) & 0x7fff);
		temp = mult(temp, 9864);                    /* change to base 10 */
		temp = shr(temp, 3);                        /* change to Q12 */
		temp = pow10_fxp(temp, 14);
		L_sum = L_shl(L_mult(ksi_min, temp), 1);    /* Now Q31 */
		temp = extract_h(L_sum);

		/*	if (ksi_min_new > 0.25) */
		if (comp_data_shift(temp, shift, 8192, 0) > 0)
			ksi_min_new = 8192;
		else
			ksi_min_new = shl(temp, shift);
	}

	return(ksi_min_new);
}


/*****************************************************************************/
/* Subroutine smoothing_win: applies the Parzen window.  The window applies  */
/* an inverse trapezoid window and wtr_front[] supplies the coefficients for */
/* the two edges.                                                            */
/*****************************************************************************/
static void		smoothing_win(Shortword initial_noise[])
{
	register Shortword	i;
	const static Shortword	wtr_front[WTR_FRONT_LEN] = {               /* Q15 */
		32767, 32582, 32048, 31202, 30080, 28718, 27152, 25418,
		23552, 21590, 19568, 17522, 15488, 13502, 11600,  9818,
		 8192,  6750,  5488,  4394,  3456,  2662,  2000,  1458,
		 1024,   686,   432,   250,   128,    54,    16,     2
	};


	for (i = 1; i < WTR_FRONT_LEN; i++)
		initial_noise[i] = mult(initial_noise[i], wtr_front[i]);
	for (i = ENH_WINLEN - WTR_FRONT_LEN + 1; i < ENH_WINLEN; i++)
		initial_noise[i] = mult(initial_noise[i], wtr_front[ENH_WINLEN - i]);

	/* Clearing the central part of initial_noise[]. */
	v_zap(&(initial_noise[WTR_FRONT_LEN]), ENH_WINLEN - 2*WTR_FRONT_LEN + 1);
}


/***************************************************************************/
/*	Subroutine smoothed_periodogram: compute short time psd with optimal   */
/*						 smoothing										   */
/***************************************************************************/
static void		smoothed_periodogram(Shortword YY_av, Shortword yy_shift)
{
	register Shortword	i;
	Shortword	smoothed_av, alphacorr_new, alpha_N_min_1, alpha_num;
	Shortword	smav_shift, shift, temp, temp1, tmpns, tmpalpha;
	Shortword	noise__shift, temp_shift, tmpns_shift, max_shift;
	Longword	L_sum, L_max, L_tmp;


	/* ---- compute smoothed_av ---- */
	max_shift = SW_MIN;
	for (i = 0; i < ENH_VEC_LENF; i++){
		if (sm_shift[i] > max_shift)
			max_shift = sm_shift[i];
	}

	L_sum = L_shl(L_deposit_l(smoothedspect[0]),
				  sub(7, sub(max_shift, sm_shift[0])));
	L_sum = L_add(L_sum, L_shl(L_deposit_l(smoothedspect[ENH_VEC_LENF - 1]),
				  sub(7, sub(max_shift, sm_shift[ENH_VEC_LENF - 1]))));
	for (i = 1; i < ENH_VEC_LENF - 1; i++)
		L_sum = L_add(L_sum, L_shl(L_deposit_l(smoothedspect[i]),
					  sub(8, sub(max_shift, sm_shift[i]))));

	/* Now L_sum contains smoothed_av.  Here we do not multiply L_sum with    */
	/* win_len_inv because we do not do this on YY_av either.                 */

	if (L_sum == 0)
		L_sum = 1;
	temp = norm_l(L_sum);
	temp = sub(temp, 1);                     /* make sure smoothed_av < YY_av */
	smoothed_av = extract_h(L_shl(L_sum, temp));
	smav_shift = sub(add(max_shift, 1), temp);

	/* ---- alphacorr_new = smoothed_av / YY_av - 1.0 ---- */
	alphacorr_new = divide_s(smoothed_av, YY_av);
	shift = sub(smav_shift, yy_shift);
	if (shift <= 0){
		if (shift > -15)
			alphacorr_new = sub(shl(alphacorr_new, shift), ONE_Q15);
		else
			alphacorr_new = negate(ONE_Q15);
		shift = 0;
	} else {
		if (shift < 15)
			alphacorr_new = sub(alphacorr_new, shr(ONE_Q15, shift));
	}

	/* -- alphacorr_new = 1.0 / (1.0 + alphacorr_new * alphacorr_new) -- */
	alphacorr_new = mult(alphacorr_new, alphacorr_new);
	alphacorr_new = shr(alphacorr_new, 1);          /* avoid overflow when +1 */
	shift = shl(shift, 1);
	if (shift < 15)
		alphacorr_new = add(alphacorr_new, shr(ONE_Q14, shift));
	if (alphacorr_new == 0)
		alphacorr_new = ONE_Q15;
	else {
		if (shift < 15){
			temp = shr(ONE_Q14, shift);
			alphacorr_new = divide_s(temp, alphacorr_new);
		} else                                                   /* too small */
			alphacorr_new = 0;
	}

	/* -- alphacorr_new > 0.7 ? alphacorr_new : 0.7 -- */
	if (alphacorr_new < X07_Q15)
		alphacorr_new = X07_Q15;

	/* -- alphacorr = 0.7*alphacorr + 0.3*alphacorr_new -- */
	alphacorr = extract_h(L_add(L_mult(X07_Q15, alphacorr),
								L_mult(X03_Q15, alphacorr_new)));

	/* -- compute alpha_N_min_1 -- */
	/* -- alpha_N_min_1 = pow(SN_LT, PDECAY_NUM) */
	if (comp_data_shift(SN_LT, SN_LT_shift, 16384, 15) > 0)
		L_sum = 536870912L;                                            /* Q15 */
	else if (comp_data_shift(SN_LT, SN_LT_shift, 32, 15) < 0)
		L_sum = 1048576L;
	else
		L_sum = L_shl(L_deposit_l(SN_LT), SN_LT_shift);
	temp = L_log10_fxp(L_sum, 15);                              /* output Q11 */
	temp = shl(temp, 1);                                    /* convert to Q12 */
	temp = mult(temp, (Shortword) PDECAY_NUM);
	alpha_N_min_1 = pow10_fxp(temp, 15);

	/* -- alpha_N_min_1 = (0.3 < alpha_N_min_1 ? 0.3 : alpha_N_min_1) -- */
	/* -- alpha_N_min_1 = (alpha_N_min_1 < 0.05 ? 0.05 : alpha_N_min_1) -- */
	if (alpha_N_min_1 > X03_Q15)
		alpha_N_min_1 = X03_Q15;
	else if (alpha_N_min_1 < X005_Q15)
		alpha_N_min_1 = X005_Q15;

	/* -- alpha_num = ALPHA_N_MAX * alphacorr -- */
	alpha_num = mult(ALPHA_N_MAX, alphacorr);

	/* -- compute smoothed spectrum -- */
	L_max = 0;
	for (i = 0; i < ENH_VEC_LENF; i++){
		tmpns = noisespect[i];

		/*	temp = smoothedspect[i] - tmpns; */
		shift = sub(sm_shift[i], noise_shift[i]);
		if (shift > 0){
			L_tmp = L_sub(L_deposit_h(smoothedspect[i]),
						  L_shr(L_deposit_h(noisespect[i]), shift));
			shift = sm_shift[i];
		} else {
			L_tmp = L_sub(L_shr(L_deposit_h(smoothedspect[i]), abs_s(shift)),
						  L_deposit_h(tmpns));
			shift = noise_shift[i];
		}

		temp1 = norm_l(L_tmp);
		temp = extract_h(L_shl(L_tmp, temp1));
		shift = sub(shift, temp1);

		/*	tmpns = tmpns * tmpns; */
		L_sum = L_mult(tmpns, tmpns);                /* noise square, shift*2 */
		noise__shift = norm_l(L_sum);
		tmpns = extract_h(L_shl(L_sum, noise__shift));
		tmpns_shift = sub(shl(noise_shift[i], 1), noise__shift);
		noisespect2[i] = tmpns;
		noise2_shift[i] = tmpns_shift;

		if (temp == SW_MIN)
			L_sum = 0x7fffffff;
		else
			L_sum = L_mult(temp, temp);
		noise__shift = norm_l(L_sum);
		temp = extract_h(L_shl(L_sum, noise__shift));
		if (temp == 0)
			temp_shift = -20;
		else
			temp_shift = sub(shl(shift, 1), noise__shift);

		/* -- tmpalpha = alpha_num * tmpns / (tmpns + temp * temp) -- */
		temp1 = sub(temp_shift, tmpns_shift);
		if (temp1 > 0){
			temp = shr(temp, 1);
			tmpns = shr(tmpns, add(temp1, 1));
		} else {
			tmpns = shr(tmpns, 1);
			temp = shl(temp, sub(temp1, 1));
		}
		tmpalpha = divide_s(tmpns, add(temp, tmpns));
		tmpalpha = mult(tmpalpha, alpha_num);

		/* tmpalpha = (tmpalpha < alpha_N_min_1 ? alpha_N_min_1 : tmpalpha) */
		if (tmpalpha < alpha_N_min_1)
			tmpalpha = alpha_N_min_1;

		temp = sub(ONE_Q15, tmpalpha);                        /* 1 - tmpalpha */
		alpha_var[i] = tmpalpha;                                /* save alpha */

		/*	smoothedspect[i] = tmpalpha * smoothedspect[i] + */
		/*                     (1 - tmpalpha) * YY[i]; */
		smav_shift = sub(YY_shift[i], sm_shift[i]);
		if (smav_shift > 0){
			L_sum = L_shr(L_mult(tmpalpha, smoothedspect[i]), smav_shift);
			L_sum = L_add(L_sum, L_mult(temp, YY[i]));
			sm_shift[i] = YY_shift[i];
		} else {
			L_sum = L_mult(tmpalpha, smoothedspect[i]);
			L_sum = L_add(L_sum, L_shl(L_mult(temp, YY[i]), smav_shift));
		}
		if (L_sum < 1)
			L_sum = 1;
		shift = norm_l(L_sum);
		smoothedspect[i] = extract_h(L_shl(L_sum, shift));
		sm_shift[i] = sub(sm_shift[i], shift);
	}
}


/*****************************************************************************/
/* Subroutine normalized_variance: compute variance of smoothed periodogram, */
/* normalize it, and use it to compute a biased smoothed periodogram		 */
/*****************************************************************************/
static void		bias_compensation(Shortword biased_spect[],
								  Shortword bias_shift[],
								  Shortword biased_spect_sub[],
								  Shortword bias_sub_shift[])
{
	register Shortword	i;
	Shortword	tmp, tmp1, tmp2, tmp5;
	Shortword	beta_var, var_rel_av_sqrt;
	Shortword	av__shift, av2__shift, shift1, shift3, shift4;
	Longword	L_max1, L_max2, L_sum, L_var_sum, L_tmp3, L_tmp4;


	/* ---- compute variance of smoothed psd ---- */
	L_var_sum = 0;
	for (i = 0; i < ENH_VEC_LENF; i++){

		/*	tmp1 = alpha_var[i]*alpha_var[i]; */
		/*	beta_var = (tmp1 > 0.8 ? 0.8 : tmp1); */
		beta_var = mult(alpha_var[i], alpha_var[i]);
		if (beta_var > X08_Q15)
			beta_var = X08_Q15;

		/*	tmp2 = (1 - beta_var) * smoothedspect[i]; */
		L_sum = L_mult(sub(ONE_Q15, beta_var), smoothedspect[i]);

		/*	var_sp_av[i] = beta_var * var_sp_av[i] + tmp2; */
		av__shift = sub(sm_shift[i], av_shift[i]);
		if (av__shift > 0){
			L_max1 = L_add(L_shr(L_mult(beta_var, var_sp_av[i]), av__shift),
						   L_sum);
			av_shift[i] = sm_shift[i];
		} else
			L_max1 = L_add(L_mult(beta_var, var_sp_av[i]),
						   L_shl(L_sum, av__shift));
		if (L_max1 < 1)
			L_max1 = 1;
		shift1 = norm_l(L_max1);
		var_sp_av[i] = extract_h(L_shl(L_max1, shift1));
		av_shift[i] = sub(av_shift[i], shift1);

		/*	tmp4 = tmp2 * smoothedspect[i]; */
		/*	var_sp_2[i] = beta_var * var_sp_2[i] + tmp4; */
		av2__shift = sub(shl(sm_shift[i], 1), av2_shift[i]);
		if (av2__shift > 0){
			L_max2 = L_add(L_shr(L_mult(beta_var, var_sp_2[i]), av2__shift),
						   L_mpy_ls(L_sum, smoothedspect[i]));
			av2_shift[i] = shl(sm_shift[i], 1);
		} else
			L_max2 = L_add(L_mult(beta_var, var_sp_2[i]),
						   L_shl(L_mpy_ls(L_sum, smoothedspect[i]),
								 av2__shift));

		if (L_max2 < 1)
			L_max2 = 1;
		shift1 = norm_l(L_max2);
		var_sp_2[i] = extract_h(L_shl(L_max2, shift1));
		av2_shift[i] = sub(av2_shift[i], shift1);

		/*	tmp3 = var_sp_av[i] * var_sp_av[i]; */
		L_sum = L_mult(var_sp_av[i], var_sp_av[i]);

		/*	var_sp = var_sp_2[i] - tmp3; */
		shift3 = sub(av2_shift[i], shl(av_shift[i], 1));
		if (shift3 > 0){
			L_sum = L_sub(L_deposit_h(var_sp_2[i]), L_shr(L_sum, shift3));
			shift4 = av2_shift[i];
		} else {
			L_sum = L_sub(L_shl(L_deposit_h(var_sp_2[i]), shift3), L_sum);
			shift4 = shl(av_shift[i], 1);
		}
		shift1 = sub(norm_l(L_sum), 1);
		tmp1 = extract_h(L_shl(L_sum, shift1));
		tmp = sub(sub(shift4, shift1), noise2_shift[i]);

		/*	tmp1 = var_sp / noisespect2[i]; */
		var_rel[i] = divide_s(tmp1, noisespect2[i]);

		/*	var_rel[i] = (tmp1 > 0.5 ? 0.5 : tmp1); */
		if (comp_data_shift(var_rel[i], tmp, X05_Q15, 0) > 0)
			var_rel[i] = X05_Q15;
		else
			var_rel[i] = shl(var_rel[i], tmp);
		if (var_rel[i] < 0)
			var_rel[i] = 0;

		/*	var_sum += var_rel[i]; */
		L_var_sum = L_add(L_var_sum, L_deposit_l(var_rel[i]));
	}

	/*	var_rel_av = (2 * var_sum - var_rel[0] - var_rel[ENH_VEC_LENF - 1]);  */
	L_var_sum = L_shl(L_var_sum, 1);
	L_var_sum = L_sub(L_var_sum, L_deposit_l(var_rel[0]));
	L_var_sum = L_sub(L_var_sum, L_deposit_l(var_rel[ENH_VEC_LENF - 1]));
	var_rel_av = extract_l(L_shr(L_var_sum, 8));                       /* Q15 */

	/*	var_rel_av = (var_rel_av < 0 ? 0 : var_rel_av); */
	if (var_rel_av < 0)
		var_rel_av = 0;

	var_rel_av_sqrt = mult(X15_Q13, sqrt_Q15(var_rel_av));			   /* Q13 */
	var_rel_av_sqrt = add(ONE_Q13, var_rel_av_sqrt);                   /* Q13 */

	/*	tmp1 = var_rel_av_sqrt * Fvar; */
	/*	tmp2 = var_rel_av_sqrt * Fvar_sub; */
	tmp1 = extract_h(L_shl(L_mult(var_rel_av_sqrt, FVAR), 1));         /* Q10 */
	tmp2 = extract_h(L_shl(L_mult(var_rel_av_sqrt, FVAR_SUB), 1));     /* Q13 */

	L_max1 = 0;                                         /* for biased_spect[] */
	L_max2 = 0;                                     /* for biased_spect_sub[] */

	for (i = 0; i < ENH_VEC_LENF; i++){
		L_tmp3 = L_mult(var_rel_av_sqrt, smoothedspect[i]);            /* Q29 */
		tmp5 = var_rel[i];                                             /* Q15 */
		L_tmp4 = L_mult(tmp5, smoothedspect[i]);                       /* Q31 */

		/* quadratic approximation */
		tmp = add(MINV, shr(tmp5, 1));                                 /* Q14 */
		tmp = add(MINV2, shr(mult(tmp5, tmp), 1));                     /* Q13 */