/* ** Copyright 2003-2010, VisualOn, Inc. ** ** Licensed under the Apache License, Version 2.0 (the "License"); ** you may not use this file except in compliance with the License. ** You may obtain a copy of the License at ** ** http://www.apache.org/licenses/LICENSE-2.0 ** ** Unless required by applicable law or agreed to in writing, software ** distributed under the License is distributed on an "AS IS" BASIS, ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. ** See the License for the specific language governing permissions and ** limitations under the License. */ /*********************************************************************** * File: levinson.c * * * * Description:LEVINSON-DURBIN algorithm in double precision * * * ************************************************************************/ /*---------------------------------------------------------------------------* * LEVINSON.C * *---------------------------------------------------------------------------* * * * LEVINSON-DURBIN algorithm in double precision * * * * * * Algorithm * * * * R[i] autocorrelations. * * A[i] filter coefficients. * * K reflection coefficients. * * Alpha prediction gain. * * * * Initialization: * * A[0] = 1 * * K = -R[1]/R[0] * * A[1] = K * * Alpha = R[0] * (1-K**2] * * * * Do for i = 2 to M * * * * S = SUM ( R[j]*A[i-j] ,j=1,i-1 ) + R[i] * * * * K = -S / Alpha * * * * An[j] = A[j] + K*A[i-j] for j=1 to i-1 * * where An[i] = new A[i] * * An[i]=K * * * * Alpha=Alpha * (1-K**2) * * * * END * * * * Remarks on the dynamics of the calculations. * * * * The numbers used are in double precision in the following format : * * A = AH <<16 + AL<<1. AH and AL are 16 bit signed integers. * * Since the LSB's also contain a sign bit, this format does not * * correspond to standard 32 bit integers. We use this format since * * it allows fast execution of multiplications and divisions. * * * * "DPF" will refer to this special format in the following text. * * See oper_32b.c * * * * The R[i] were normalized in routine AUTO (hence, R[i] < 1.0). * * The K[i] and Alpha are theoretically < 1.0. * * The A[i], for a sampling frequency of 8 kHz, are in practice * * always inferior to 16.0. * * * * These characteristics allow straigthforward fixed-point * * implementation. We choose to represent the parameters as * * follows : * * * * R[i] Q31 +- .99.. * * K[i] Q31 +- .99.. * * Alpha Normalized -> mantissa in Q31 plus exponent * * A[i] Q27 +- 15.999.. * * * * The additions are performed in 32 bit. For the summation used * * to calculate the K[i], we multiply numbers in Q31 by numbers * * in Q27, with the result of the multiplications in Q27, * * resulting in a dynamic of +- 16. This is sufficient to avoid * * overflow, since the final result of the summation is * * necessarily < 1.0 as both the K[i] and Alpha are * * theoretically < 1.0. * *___________________________________________________________________________*/ #include "typedef.h" #include "basic_op.h" #include "oper_32b.h" #include "acelp.h" #define M 16 #define NC (M/2) void Init_Levinson( Word16 * mem /* output :static memory (18 words) */ ) { Set_zero(mem, 18); /* old_A[0..M-1] = 0, old_rc[0..1] = 0 */ return; } void Levinson( Word16 Rh[], /* (i) : Rh[M+1] Vector of autocorrelations (msb) */ Word16 Rl[], /* (i) : Rl[M+1] Vector of autocorrelations (lsb) */ Word16 A[], /* (o) Q12 : A[M] LPC coefficients (m = 16) */ Word16 rc[], /* (o) Q15 : rc[M] Reflection coefficients. */ Word16 * mem /* (i/o) :static memory (18 words) */ ) { Word32 i, j; Word16 hi, lo; Word16 Kh, Kl; /* reflection coefficient; hi and lo */ Word16 alp_h, alp_l, alp_exp; /* Prediction gain; hi lo and exponent */ Word16 Ah[M + 1], Al[M + 1]; /* LPC coef. in double prec. */ Word16 Anh[M + 1], Anl[M + 1]; /* LPC coef.for next iteration in double prec. */ Word32 t0, t1, t2; /* temporary variable */ Word16 *old_A, *old_rc; /* Last A(z) for case of unstable filter */ old_A = mem; old_rc = mem + M; /* K = A[1] = -R[1] / R[0] */ t1 = ((Rh[1] << 16) + (Rl[1] << 1)); /* R[1] in Q31 */ t2 = L_abs(t1); /* abs R[1] */ t0 = Div_32(t2, Rh[0], Rl[0]); /* R[1]/R[0] in Q31 */ if (t1 > 0) t0 = -t0; /* -R[1]/R[0] */ Kh = t0 >> 16; Kl = (t0 & 0xffff)>>1; rc[0] = Kh; t0 = (t0 >> 4); /* A[1] in Q27 */ Ah[1] = t0 >> 16; Al[1] = (t0 & 0xffff)>>1; /* Alpha = R[0] * (1-K**2) */ t0 = Mpy_32(Kh, Kl, Kh, Kl); /* K*K in Q31 */ t0 = L_abs(t0); /* Some case <0 !! */ t0 = vo_L_sub((Word32) 0x7fffffffL, t0); /* 1 - K*K in Q31 */ hi = t0 >> 16; lo = (t0 & 0xffff)>>1; t0 = Mpy_32(Rh[0], Rl[0], hi, lo); /* Alpha in Q31 */ /* Normalize Alpha */ alp_exp = norm_l(t0); t0 = (t0 << alp_exp); alp_h = t0 >> 16; alp_l = (t0 & 0xffff)>>1; /*--------------------------------------* * ITERATIONS I=2 to M * *--------------------------------------*/ for (i = 2; i <= M; i++) { /* t0 = SUM ( R[j]*A[i-j] ,j=1,i-1 ) + R[i] */ t0 = 0; for (j = 1; j < i; j++) t0 = vo_L_add(t0, Mpy_32(Rh[j], Rl[j], Ah[i - j], Al[i - j])); t0 = t0 << 4; /* result in Q27 -> convert to Q31 */ /* No overflow possible */ t1 = ((Rh[i] << 16) + (Rl[i] << 1)); t0 = vo_L_add(t0, t1); /* add R[i] in Q31 */ /* K = -t0 / Alpha */ t1 = L_abs(t0); t2 = Div_32(t1, alp_h, alp_l); /* abs(t0)/Alpha */ if (t0 > 0) t2 = -t2; /* K =-t0/Alpha */ t2 = (t2 << alp_exp); /* denormalize; compare to Alpha */ Kh = t2 >> 16; Kl = (t2 & 0xffff)>>1; rc[i - 1] = Kh; /* Test for unstable filter. If unstable keep old A(z) */ if (abs_s(Kh) > 32750) { A[0] = 4096; /* Ai[0] not stored (always 1.0) */ for (j = 0; j < M; j++) { A[j + 1] = old_A[j]; } rc[0] = old_rc[0]; /* only two rc coefficients are needed */ rc[1] = old_rc[1]; return; } /*------------------------------------------* * Compute new LPC coeff. -> An[i] * * An[j]= A[j] + K*A[i-j] , j=1 to i-1 * * An[i]= K * *------------------------------------------*/ for (j = 1; j < i; j++) { t0 = Mpy_32(Kh, Kl, Ah[i - j], Al[i - j]); t0 = vo_L_add(t0, ((Ah[j] << 16) + (Al[j] << 1))); Anh[j] = t0 >> 16; Anl[j] = (t0 & 0xffff)>>1; } t2 = (t2 >> 4); /* t2 = K in Q31 ->convert to Q27 */ VO_L_Extract(t2, &Anh[i], &Anl[i]); /* An[i] in Q27 */ /* Alpha = Alpha * (1-K**2) */ t0 = Mpy_32(Kh, Kl, Kh, Kl); /* K*K in Q31 */ t0 = L_abs(t0); /* Some case <0 !! */ t0 = vo_L_sub((Word32) 0x7fffffffL, t0); /* 1 - K*K in Q31 */ hi = t0 >> 16; lo = (t0 & 0xffff)>>1; t0 = Mpy_32(alp_h, alp_l, hi, lo); /* Alpha in Q31 */ /* Normalize Alpha */ j = norm_l(t0); t0 = (t0 << j); alp_h = t0 >> 16; alp_l = (t0 & 0xffff)>>1; alp_exp += j; /* Add normalization to alp_exp */ /* A[j] = An[j] */ for (j = 1; j <= i; j++) { Ah[j] = Anh[j]; Al[j] = Anl[j]; } } /* Truncate A[i] in Q27 to Q12 with rounding */ A[0] = 4096; for (i = 1; i <= M; i++) { t0 = (Ah[i] << 16) + (Al[i] << 1); old_A[i - 1] = A[i] = vo_round((t0 << 1)); } old_rc[0] = rc[0]; old_rc[1] = rc[1]; return; }