diff options
Diffstat (limited to 'media/libstagefright/codecs/amrwbenc/src/levinson.c')
| -rw-r--r-- | media/libstagefright/codecs/amrwbenc/src/levinson.c | 250 |
1 files changed, 250 insertions, 0 deletions
diff --git a/media/libstagefright/codecs/amrwbenc/src/levinson.c b/media/libstagefright/codecs/amrwbenc/src/levinson.c new file mode 100644 index 0000000..8bc6f62 --- /dev/null +++ b/media/libstagefright/codecs/amrwbenc/src/levinson.c @@ -0,0 +1,250 @@ +/*
+ ** 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;
+}
+
+
+
|