opencore/codecs_v2/audio/gsm_amr/amr_nb/enc/src/levinson.cpp

   1 /* ------------------------------------------------------------------
   2  * Copyright (C) 1998-2009 PacketVideo
   3  *
   4  * Licensed under the Apache License, Version 2.0 (the "License");
   5  * you may not use this file except in compliance with the License.
   6  * You may obtain a copy of the License at
   7  *
   8  *      http://www.apache.org/licenses/LICENSE-2.0
   9  *
  10  * Unless required by applicable law or agreed to in writing, software
  11  * distributed under the License is distributed on an "AS IS" BASIS,
  12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
  13  * express or implied.
  14  * See the License for the specific language governing permissions
  15  * and limitations under the License.
  16  * -------------------------------------------------------------------
  17  */
  18 /****************************************************************************************
  19 Portions of this file are derived from the following 3GPP standard:
  20
  21     3GPP TS 26.073
  22     ANSI-C code for the Adaptive Multi-Rate (AMR) speech codec
  23     Available from http://www.3gpp.org
  24
  25 (C) 2004, 3GPP Organizational Partners (ARIB, ATIS, CCSA, ETSI, TTA, TTC)
  26 Permission to distribute, modify and use this file under the standard license
  27 terms listed above has been obtained from the copyright holder.
  28 ****************************************************************************************/
  29 /*
  30 ------------------------------------------------------------------------------
  31
  32
  33
  34  Filename: levinson.cpp
  35  Funtions: Levinson_init
  36            Levinson_reset
  37            Levinson_exit
  38            Levinson
  39
  40 ------------------------------------------------------------------------------
  41  MODULE DESCRIPTION
  42
  43  This file contains the function the implements the Levinson-Durbin algorithm
  44  using double-precision arithmetic. This file also includes functions to
  45  initialize, allocate, and deallocate memory used by the Levinson function.
  46
  47 ------------------------------------------------------------------------------
  48 */
  49
  50 /*----------------------------------------------------------------------------
  51 ; INCLUDES
  52 ----------------------------------------------------------------------------*/
  53 #include "levinson.h"
  54 #include "basicop_malloc.h"
  55 #include "basic_op.h"
  56 #include "l_abs.h"
  57 #include "div_32.h"
  58 #include "cnst.h"
  59 #include "oscl_mem.h"
  60
  61 /*----------------------------------------------------------------------------
  62 ; MACROS
  63 ; Define module specific macros here
  64 ----------------------------------------------------------------------------*/
  65
  66 /*----------------------------------------------------------------------------
  67 ; DEFINES
  68 ; Include all pre-processor statements here. Include conditional
  69 ; compile variables also.
  70 ----------------------------------------------------------------------------*/
  71
  72 /*----------------------------------------------------------------------------
  73 ; LOCAL FUNCTION DEFINITIONS
  74 ; Function Prototype declaration
  75 ----------------------------------------------------------------------------*/
  76
  77 /*----------------------------------------------------------------------------
  78 ; LOCAL VARIABLE DEFINITIONS
  79 ; Variable declaration - defined here and used outside this module
  80 ----------------------------------------------------------------------------*/
  81
  82
  83 /*
  84 ------------------------------------------------------------------------------
  85  FUNCTION NAME: Levinson_init
  86 ------------------------------------------------------------------------------
  87  INPUT AND OUTPUT DEFINITIONS
  88
  89  Inputs:
  90     state = pointer to an array of pointers to structures of type
  91             LevinsonState
  92
  93  Outputs:
  94     pointer pointed to by state points to the newly allocated memory to
  95       be used by Levinson function
  96
  97  Returns:
  98     return_value = 0, if initialization was successful; -1, otherwise (int)
  99
 100  Global Variables Used:
 101     None
 102
 103  Local Variables Needed:
 104     None
 105
 106 ------------------------------------------------------------------------------
 107  FUNCTION DESCRIPTION
 108
 109  This function allocates and initializes the state memory used by the
 110  Levinson function.
 111
 112 ------------------------------------------------------------------------------
 113  REQUIREMENTS
 114
 115  None
 116
 117 ------------------------------------------------------------------------------
 118  REFERENCES
 119
 120  levinson.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001
 121
 122 ------------------------------------------------------------------------------
 123  PSEUDO-CODE
 124
 125 int Levinson_init (LevinsonState **state)
 126 {
 127   LevinsonState* s;
 128
 129   if (state == (LevinsonState **) NULL){
 130       //fprint(stderr, "Levinson_init: invalid parameter\n");
 131       return -1;
 132   }
 133   *state = NULL;
 134
 135   // allocate memory
 136   if ((s= (LevinsonState *) malloc(sizeof(LevinsonState))) == NULL){
 137       //fprint(stderr, "Levinson_init: can not malloc state structure\n");
 138       return -1;
 139   }
 140
 141   Levinson_reset(s);
 142   *state = s;
 143
 144   return 0;
 145 }
 146
 147 ------------------------------------------------------------------------------
 148  CAUTION [optional]
 149  [State any special notes, constraints or cautions for users of this function]
 150
 151 ------------------------------------------------------------------------------
 152 */
 153
 154 Word16 Levinson_init(LevinsonState **state)
 155 {
 156     LevinsonState* s;
 157
 158     if (state == (LevinsonState **) NULL)
 159     {
 160         /*  fprint(stderr, "Levinson_init: invalid parameter\n");  */
 161         return(-1);
 162     }
 163     *state = NULL;
 164
 165     /* allocate memory */
 166     if ((s = (LevinsonState *) oscl_malloc(sizeof(LevinsonState))) == NULL)
 167     {
 168         /*  fprint(stderr, "Levinson_init:
 169                             can not malloc state structure\n");  */
 170         return(-1);
 171     }
 172
 173     Levinson_reset(s);
 174     *state = s;
 175
 176     return(0);
 177 }
 178
 179 /****************************************************************************/
 180
 181
 182 /*
 183 ------------------------------------------------------------------------------
 184  FUNCTION NAME: Levinson_reset
 185 ------------------------------------------------------------------------------
 186  INPUT AND OUTPUT DEFINITIONS
 187
 188  Inputs:
 189     state = pointer to structures of type LevinsonState
 190
 191  Outputs:
 192     old_A field of structure pointed to by state is initialized to 4096
 193       (first location) and the rest to zeros
 194
 195  Returns:
 196     return_value = 0, if reset was successful; -1, otherwise (int)
 197
 198  Global Variables Used:
 199     None
 200
 201  Local Variables Needed:
 202     None
 203
 204 ------------------------------------------------------------------------------
 205  FUNCTION DESCRIPTION
 206
 207  This function initializes the state memory used by the Levinson function to
 208  zero.
 209
 210 ------------------------------------------------------------------------------
 211  REQUIREMENTS
 212
 213  None
 214
 215 ------------------------------------------------------------------------------
 216  REFERENCES
 217
 218  levinson.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001
 219
 220 ------------------------------------------------------------------------------
 221  PSEUDO-CODE
 222
 223 int Levinson_reset (LevinsonState *state)
 224 {
 225   Word16 i;
 226
 227   if (state == (LevinsonState *) NULL){
 228       fprint(stderr, "Levinson_reset: invalid parameter\n");
 229       return -1;
 230   }
 231
 232   state->old_A[0] = 4096;
 233   for(i = 1; i < M + 1; i++)
 234       state->old_A[i] = 0;
 235
 236   return 0;
 237 }
 238
 239 ------------------------------------------------------------------------------
 240  CAUTION [optional]
 241  [State any special notes, constraints or cautions for users of this function]
 242
 243 ------------------------------------------------------------------------------
 244 */
 245
 246 Word16 Levinson_reset(LevinsonState *state)
 247 {
 248     Word16 i;
 249
 250     if (state == (LevinsonState *) NULL)
 251     {
 252         /*  fprint(stderr, "Levinson_reset: invalid parameter\n");  */
 253         return(-1);
 254     }
 255
 256     state->old_A[0] = 4096;
 257     for (i = 1; i < M + 1; i++)
 258     {
 259         state->old_A[i] = 0;
 260     }
 261
 262     return(0);
 263 }
 264
 265 /****************************************************************************/
 266
 267 /*
 268 ------------------------------------------------------------------------------
 269  FUNCTION NAME: Levinson_exit
 270 ------------------------------------------------------------------------------
 271  INPUT AND OUTPUT DEFINITIONS
 272
 273  Inputs:
 274     state = pointer to an array of pointers to structures of type
 275             LevinsonState
 276
 277  Outputs:
 278     pointer pointed to by state is set to the NULL address
 279
 280  Returns:
 281     None
 282
 283  Global Variables Used:
 284     None
 285
 286  Local Variables Needed:
 287     None
 288
 289 ------------------------------------------------------------------------------
 290  FUNCTION DESCRIPTION
 291
 292  This function deallocates the state memory used by the Levinson function.
 293
 294 ------------------------------------------------------------------------------
 295  REQUIREMENTS
 296
 297  None
 298
 299 ------------------------------------------------------------------------------
 300  REFERENCES
 301
 302  levinson.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001
 303
 304 ------------------------------------------------------------------------------
 305  PSEUDO-CODE
 306
 307 void Levinson_exit (LevinsonState **state)
 308 {
 309   if (state == NULL || *state == NULL)
 310       return;
 311
 312   // deallocate memory
 313   free(*state);
 314   *state = NULL;
 315
 316   return;
 317 }
 318
 319 ------------------------------------------------------------------------------
 320  CAUTION [optional]
 321  [State any special notes, constraints or cautions for users of this function]
 322
 323 ------------------------------------------------------------------------------
 324 */
 325
 326 void Levinson_exit(LevinsonState **state)
 327 {
 328     if (state == NULL || *state == NULL)
 329     {
 330         return;
 331     }
 332
 333     /* deallocate memory */
 334     oscl_free(*state);
 335     *state = NULL;
 336
 337     return;
 338 }
 339
 340 /****************************************************************************/
 341
 342
 343 /*
 344 ------------------------------------------------------------------------------
 345  FUNCTION NAME: Levinson
 346 ------------------------------------------------------------------------------
 347  INPUT AND OUTPUT DEFINITIONS
 348
 349  Inputs:
 350     st = pointer to structures of type LevinsonState
 351     Rh = vector containing most significant byte of
 352          autocorrelation values (Word16)
 353     Rl = vector containing least significant byte of
 354          autocorrelation values (Word16)
 355     A = vector of LPC coefficients (10th order) (Word16)
 356     rc = vector containing first four reflection coefficients (Word16)
 357     pOverflow = pointer to overflow indicator (Flag)
 358
 359  Outputs:
 360     A contains the newly calculated LPC coefficients
 361     rc contains the newly calculated reflection coefficients
 362
 363  Returns:
 364     return_value = 0 (int)
 365
 366  Global Variables Used:
 367     None
 368
 369  Local Variables Needed:
 370     None
 371
 372 ------------------------------------------------------------------------------
 373  FUNCTION DESCRIPTION
 374
 375  This function implements the Levinson-Durbin algorithm using double-
 376  precision arithmetic. This is used to compute the Linear Predictive (LP)
 377  filter parameters from the speech autocorrelation values.
 378
 379  The algorithm implemented is as follows:
 380     A[0] = 1
 381     K    = -R[1]/R[0]
 382     A[1] = K
 383     Alpha = R[0] * (1-K**2]
 384
 385     FOR  i = 2 to M
 386
 387         S =  SUM ( R[j]*A[i-j] ,j=1,i-1 ) +  R[i]
 388         K = -S / Alpha
 389
 390         FOR j = 1 to  i-1
 391             An[j] = A[j] + K*A[i-j]  where   An[i] = new A[i]
 392         ENDFOR
 393
 394         An[i]=K
 395         Alpha=Alpha * (1-K**2)
 396
 397     END
 398
 399  where:
 400     R[i] = autocorrelations
 401     A[i] = filter coefficients
 402     K = reflection coefficient
 403     Alpha = prediction gain
 404
 405 ------------------------------------------------------------------------------
 406  REQUIREMENTS
 407
 408  None
 409
 410 ------------------------------------------------------------------------------
 411  REFERENCES
 412
 413  levinson.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001
 414
 415 ------------------------------------------------------------------------------
 416  PSEUDO-CODE
 417
 418 int Levinson (
 419     LevinsonState *st,
 420     Word16 Rh[],       // i : Rh[m+1] Vector of autocorrelations (msb)
 421     Word16 Rl[],       // i : Rl[m+1] Vector of autocorrelations (lsb)
 422     Word16 A[],        // o : A[m]    LPC coefficients  (m = 10)
 423     Word16 rc[]        // o : rc[4]   First 4 reflection coefficients
 424 )
 425 {
 426     Word16 i, j;
 427     Word16 hi, lo;
 428     Word16 Kh, Kl;                // reflexion coefficient; hi and lo
 429     Word16 alp_h, alp_l, alp_exp; // Prediction gain; hi lo and exponent
 430     Word16 Ah[M + 1], Al[M + 1];  // LPC coef. in double prec.
 431     Word16 Anh[M + 1], Anl[M + 1];// LPC coef.for next iteration in double
 432                                      prec.
 433     Word32 t0, t1, t2;            // temporary variable
 434
 435     // K = A[1] = -R[1] / R[0]
 436
 437     t1 = L_Comp (Rh[1], Rl[1]);
 438     t2 = L_abs (t1);                    // abs R[1]
 439     t0 = Div_32 (t2, Rh[0], Rl[0]);     // R[1]/R[0]
 440     if (t1 > 0)
 441        t0 = L_negate (t0);             // -R[1]/R[0]
 442     L_Extract (t0, &Kh, &Kl);           // K in DPF
 443
 444     rc[0] = pv_round (t0);
 445
 446     t0 = L_shr (t0, 4);                 // A[1] in
 447     L_Extract (t0, &Ah[1], &Al[1]);     // A[1] in DPF
 448
 449     //  Alpha = R[0] * (1-K**2)
 450
 451     t0 = Mpy_32 (Kh, Kl, Kh, Kl);       // K*K
 452     t0 = L_abs (t0);                    // Some case <0 !!
 453     t0 = L_sub ((Word32) 0x7fffffffL, t0); // 1 - K*K
 454     L_Extract (t0, &hi, &lo);           // DPF format
 455     t0 = Mpy_32 (Rh[0], Rl[0], hi, lo); // Alpha in
 456
 457     // Normalize Alpha
 458
 459     alp_exp = norm_l (t0);
 460     t0 = L_shl (t0, alp_exp);
 461     L_Extract (t0, &alp_h, &alp_l);     // DPF format
 462
 463      *--------------------------------------*
 464      * ITERATIONS  I=2 to M                 *
 465      *--------------------------------------*
 466
 467     for (i = 2; i <= M; i++)
 468     {
 469        // t0 = SUM ( R[j]*A[i-j] ,j=1,i-1 ) +  R[i]
 470
 471        t0 = 0;
 472        for (j = 1; j < i; j++)
 473        {
 474           t0 = L_add (t0, Mpy_32 (Rh[j], Rl[j], Ah[i - j], Al[i - j]));
 475        }
 476        t0 = L_shl (t0, 4);
 477
 478        t1 = L_Comp (Rh[i], Rl[i]);
 479        t0 = L_add (t0, t1);            // add R[i]
 480
 481        // K = -t0 / Alpha
 482
 483        t1 = L_abs (t0);
 484        t2 = Div_32 (t1, alp_h, alp_l); // abs(t0)/Alpha
 485        if (t0 > 0)
 486           t2 = L_negate (t2);         // K =-t0/Alpha
 487        t2 = L_shl (t2, alp_exp);       // denormalize; compare to Alpha
 488        L_Extract (t2, &Kh, &Kl);       // K in DPF
 489
 490        if (sub (i, 5) < 0)
 491        {
 492           rc[i - 1] = pv_round (t2);
 493        }
 494        // Test for unstable filter. If unstable keep old A(z)
 495
 496        if (sub (abs_s (Kh), 32750) > 0)
 497        {
 498           for (j = 0; j <= M; j++)
 499           {
 500              A[j] = st->old_A[j];
 501           }
 502
 503           for (j = 0; j < 4; j++)
 504           {
 505              rc[j] = 0;
 506           }
 507
 508           return 0;
 509        }
 510         *------------------------------------------*
 511         *  Compute new LPC coeff. -> An[i]         *
 512         *  An[j]= A[j] + K*A[i-j]     , j=1 to i-1 *
 513         *  An[i]= K                                *
 514         *------------------------------------------*
 515
 516        for (j = 1; j < i; j++)
 517        {
 518           t0 = Mpy_32 (Kh, Kl, Ah[i - j], Al[i - j]);
 519           t0 = L_add(t0, L_Comp(Ah[j], Al[j]));
 520           L_Extract (t0, &Anh[j], &Anl[j]);
 521        }
 522        t2 = L_shr (t2, 4);
 523        L_Extract (t2, &Anh[i], &Anl[i]);
 524
 525        //  Alpha = Alpha * (1-K**2)
 526
 527        t0 = Mpy_32 (Kh, Kl, Kh, Kl);           // K*K
 528        t0 = L_abs (t0);                        // Some case <0 !!
 529        t0 = L_sub ((Word32) 0x7fffffffL, t0);  // 1 - K*K
 530        L_Extract (t0, &hi, &lo);               // DPF format
 531        t0 = Mpy_32 (alp_h, alp_l, hi, lo);
 532
 533        // Normalize Alpha
 534
 535        j = norm_l (t0);
 536        t0 = L_shl (t0, j);
 537        L_Extract (t0, &alp_h, &alp_l);         // DPF format
 538        alp_exp = add (alp_exp, j);             // Add normalization to
 539                                                   alp_exp
 540
 541        // A[j] = An[j]
 542
 543        for (j = 1; j <= i; j++)
 544        {
 545           Ah[j] = Anh[j];
 546           Al[j] = Anl[j];
 547        }
 548     }
 549
 550     A[0] = 4096;
 551     for (i = 1; i <= M; i++)
 552     {
 553        t0 = L_Comp (Ah[i], Al[i]);
 554        st->old_A[i] = A[i] = pv_round (L_shl (t0, 1));
 555     }
 556
 557     return 0;
 558 }
 559
 560 ------------------------------------------------------------------------------
 561  CAUTION [optional]
 562  [State any special notes, constraints or cautions for users of this function]
 563
 564 ------------------------------------------------------------------------------
 565 */
 566
 567 Word16 Levinson(
 568     LevinsonState *st,
 569     Word16 Rh[],       /* i : Rh[m+1] Vector of autocorrelations (msb) */
 570     Word16 Rl[],       /* i : Rl[m+1] Vector of autocorrelations (lsb) */
 571     Word16 A[],        /* o : A[m]    LPC coefficients  (m = 10)       */
 572     Word16 rc[],       /* o : rc[4]   First 4 reflection coefficients  */
 573     Flag   *pOverflow
 574 )
 575 {
 576     register Word16 i;
 577     register Word16 j;
 578     Word16 hi;
 579     Word16 lo;
 580     Word16 Kh;                    /* reflexion coefficient; hi and lo   */
 581     Word16 Kl;
 582     Word16 alp_h;                 /* Prediction gain; hi lo and exponent*/
 583     Word16 alp_l;
 584     Word16 alp_exp;
 585     Word16 Ah[M + 1];             /* LPC coef. in double prec.          */
 586     Word16 Al[M + 1];
 587     Word16 Anh[M + 1];            /* LPC coef.for next iteration in     */
 588     Word16 Anl[M + 1];            /* double prec.                       */
 589     register Word32 t0;           /* temporary variable                 */
 590     register Word32 t1;           /* temporary variable                 */
 591     register Word32 t2;           /* temporary variable                 */
 592
 593     Word16 *p_Rh;
 594     Word16 *p_Rl;
 595     Word16 *p_Ah;
 596     Word16 *p_Al;
 597     Word16 *p_Anh;
 598     Word16 *p_Anl;
 599     Word16 *p_A;
 600
 601     /* K = A[1] = -R[1] / R[0] */
 602     t1 = ((Word32) * (Rh + 1)) << 16;
 603     t1 += *(Rl + 1) << 1;
 604
 605     t2 = L_abs(t1);         /* abs R[1] - required by Div_32 */
 606     t0 = Div_32(t2, *Rh, *Rl, pOverflow);  /* R[1]/R[0]        */
 607
 608     if (t1 > 0)
 609     {
 610         t0 = L_negate(t0);  /* -R[1]/R[0]       */
 611     }
 612
 613     /* K in DPF         */
 614     Kh = (Word16)(t0 >> 16);
 615     Kl = (Word16)((t0 >> 1) - ((Word32)(Kh) << 15));
 616
 617     *rc = pv_round(t0, pOverflow);
 618
 619     t0 = t0 >> 4;
 620
 621     /* A[1] in DPF      */
 622     *(Ah + 1) = (Word16)(t0 >> 16);
 623
 624     *(Al + 1) = (Word16)((t0 >> 1) - ((Word32)(*(Ah + 1)) << 15));
 625
 626     /*  Alpha = R[0] * (1-K**2) */
 627     t0 = Mpy_32(Kh, Kl, Kh, Kl, pOverflow);         /* K*K              */
 628     t0 = L_abs(t0);                                 /* Some case <0 !!  */
 629     t0 = 0x7fffffffL - t0;                          /* 1 - K*K          */
 630
 631     /* DPF format       */
 632     hi = (Word16)(t0 >> 16);
 633     lo = (Word16)((t0 >> 1) - ((Word32)(hi) << 15));
 634
 635     t0 = Mpy_32(*Rh, *Rl, hi, lo, pOverflow);      /* Alpha in         */
 636
 637     /* Normalize Alpha */
 638
 639     alp_exp = norm_l(t0);
 640     t0 = t0 << alp_exp;
 641
 642     /* DPF format       */
 643     alp_h = (Word16)(t0 >> 16);
 644     alp_l = (Word16)((t0 >> 1) - ((Word32)(alp_h) << 15));
 645
 646     /*--------------------------------------*
 647     * ITERATIONS  I=2 to M                 *
 648     *--------------------------------------*/
 649
 650     for (i = 2; i <= M; i++)
 651     {
 652         /* t0 = SUM ( R[j]*A[i-j] ,j=1,i-1 ) +  R[i] */
 653
 654         t0 = 0;
 655         p_Rh = &Rh[1];
 656         p_Rl = &Rl[1];
 657         p_Ah = &Ah[i-1];
 658         p_Al = &Al[i-1];
 659         for (j = 1; j < i; j++)
 660         {
 661             t0 += (((Word32) * (p_Rh)* *(p_Al--)) >> 15);
 662             t0 += (((Word32) * (p_Rl++)* *(p_Ah)) >> 15);
 663             t0 += ((Word32) * (p_Rh++)* *(p_Ah--));
 664         }
 665
 666         t0 = t0 << 5;
 667
 668         t1 = ((Word32) * (Rh + i) << 16) + ((Word32)(*(Rl + i)) << 1);
 669         t0 += t1;
 670
 671         /* K = -t0 / Alpha */
 672
 673         t1 = L_abs(t0);
 674         t2 = Div_32(t1, alp_h, alp_l, pOverflow);  /* abs(t0)/Alpha        */
 675
 676         if (t0 > 0)
 677         {
 678             t2 = L_negate(t2);          /* K =-t0/Alpha     */
 679         }
 680
 681         t2 = L_shl(t2, alp_exp, pOverflow);  /* denormalize; compare to Alpha */
 682         Kh = (Word16)(t2 >> 16);
 683         Kl = (Word16)((t2 >> 1) - ((Word32)(Kh) << 15));
 684
 685         if (i < 5)
 686         {
 687             *(rc + i - 1) = (Word16)((t2 + 0x00008000L) >> 16);
 688         }
 689         /* Test for unstable filter. If unstable keep old A(z) */
 690         if ((abs_s(Kh)) > 32750)
 691         {
 692             oscl_memcpy(A, &(st->old_A[0]), sizeof(Word16)*(M + 1));
 693             oscl_memset(rc, 0, sizeof(Word16)*4);
 694             return(0);
 695         }
 696         /*------------------------------------------*
 697         *  Compute new LPC coeff. -> An[i]         *
 698         *  An[j]= A[j] + K*A[i-j]     , j=1 to i-1 *
 699         *  An[i]= K                                *
 700         *------------------------------------------*/
 701         p_Ah = &Ah[i-1];
 702         p_Al = &Al[i-1];
 703         p_Anh = &Anh[1];
 704         p_Anl = &Anl[1];
 705         for (j = 1; j < i; j++)
 706         {
 707             t0  = (((Word32)Kh* *(p_Al--)) >> 15);
 708             t0 += (((Word32)Kl* *(p_Ah)) >> 15);
 709             t0 += ((Word32)Kh* *(p_Ah--));
 710
 711             t0 += (Ah[j] << 15) + Al[j];
 712
 713             *(p_Anh) = (Word16)(t0 >> 15);
 714             *(p_Anl++) = (Word16)(t0 - ((Word32)(*(p_Anh++)) << 15));
 715         }
 716
 717         *(p_Anh) = (Word16)(t2 >> 20);
 718         *(p_Anl) = (Word16)((t2 >> 5) - ((Word32)(*(Anh + i)) << 15));
 719
 720         /*  Alpha = Alpha * (1-K**2) */
 721
 722         t0 = Mpy_32(Kh, Kl, Kh, Kl, pOverflow);  /* K*K             */
 723         t0 = L_abs(t0);                          /* Some case <0 !! */
 724         t0 = 0x7fffffffL - t0;                   /* 1 - K*K          */
 725
 726         hi = (Word16)(t0 >> 16);
 727         lo = (Word16)((t0 >> 1) - ((Word32)(hi) << 15));
 728
 729         t0  = (((Word32)alp_h * lo) >> 15);
 730         t0 += (((Word32)alp_l * hi) >> 15);
 731         t0 += ((Word32)alp_h * hi);
 732
 733         t0 <<= 1;
 734         /* Normalize Alpha */
 735
 736         j = norm_l(t0);
 737         t0 <<= j;
 738         alp_h = (Word16)(t0 >> 16);
 739         alp_l = (Word16)((t0 >> 1) - ((Word32)(alp_h) << 15));
 740         alp_exp += j;             /* Add normalization to alp_exp */
 741
 742         /* A[j] = An[j] */
 743         oscl_memcpy(&Ah[1], &Anh[1], sizeof(Word16)*i);
 744         oscl_memcpy(&Al[1], &Anl[1], sizeof(Word16)*i);
 745     }
 746
 747     p_A = &A[0];
 748     *(p_A++) = 4096;
 749     p_Ah = &Ah[1];
 750     p_Al = &Al[1];
 751
 752     for (i = 1; i <= M; i++)
 753     {
 754         t0 = ((Word32) * (p_Ah++) << 15) + *(p_Al++);
 755         st->old_A[i] = *(p_A++) = (Word16)((t0 + 0x00002000) >> 14);
 756     }
 757
 758     return(0);
 759 }