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