• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /* -----------------------------------------------------------------------------
2 Software License for The Fraunhofer FDK AAC Codec Library for Android
3 
4 © Copyright  1995 - 2019 Fraunhofer-Gesellschaft zur Förderung der angewandten
5 Forschung e.V. All rights reserved.
6 
7  1.    INTRODUCTION
8 The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software
9 that implements the MPEG Advanced Audio Coding ("AAC") encoding and decoding
10 scheme for digital audio. This FDK AAC Codec software is intended to be used on
11 a wide variety of Android devices.
12 
13 AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient
14 general perceptual audio codecs. AAC-ELD is considered the best-performing
15 full-bandwidth communications codec by independent studies and is widely
16 deployed. AAC has been standardized by ISO and IEC as part of the MPEG
17 specifications.
18 
19 Patent licenses for necessary patent claims for the FDK AAC Codec (including
20 those of Fraunhofer) may be obtained through Via Licensing
21 (www.vialicensing.com) or through the respective patent owners individually for
22 the purpose of encoding or decoding bit streams in products that are compliant
23 with the ISO/IEC MPEG audio standards. Please note that most manufacturers of
24 Android devices already license these patent claims through Via Licensing or
25 directly from the patent owners, and therefore FDK AAC Codec software may
26 already be covered under those patent licenses when it is used for those
27 licensed purposes only.
28 
29 Commercially-licensed AAC software libraries, including floating-point versions
30 with enhanced sound quality, are also available from Fraunhofer. Users are
31 encouraged to check the Fraunhofer website for additional applications
32 information and documentation.
33 
34 2.    COPYRIGHT LICENSE
35 
36 Redistribution and use in source and binary forms, with or without modification,
37 are permitted without payment of copyright license fees provided that you
38 satisfy the following conditions:
39 
40 You must retain the complete text of this software license in redistributions of
41 the FDK AAC Codec or your modifications thereto in source code form.
42 
43 You must retain the complete text of this software license in the documentation
44 and/or other materials provided with redistributions of the FDK AAC Codec or
45 your modifications thereto in binary form. You must make available free of
46 charge copies of the complete source code of the FDK AAC Codec and your
47 modifications thereto to recipients of copies in binary form.
48 
49 The name of Fraunhofer may not be used to endorse or promote products derived
50 from this library without prior written permission.
51 
52 You may not charge copyright license fees for anyone to use, copy or distribute
53 the FDK AAC Codec software or your modifications thereto.
54 
55 Your modified versions of the FDK AAC Codec must carry prominent notices stating
56 that you changed the software and the date of any change. For modified versions
57 of the FDK AAC Codec, the term "Fraunhofer FDK AAC Codec Library for Android"
58 must be replaced by the term "Third-Party Modified Version of the Fraunhofer FDK
59 AAC Codec Library for Android."
60 
61 3.    NO PATENT LICENSE
62 
63 NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without
64 limitation the patents of Fraunhofer, ARE GRANTED BY THIS SOFTWARE LICENSE.
65 Fraunhofer provides no warranty of patent non-infringement with respect to this
66 software.
67 
68 You may use this FDK AAC Codec software or modifications thereto only for
69 purposes that are authorized by appropriate patent licenses.
70 
71 4.    DISCLAIMER
72 
73 This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright
74 holders and contributors "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES,
75 including but not limited to the implied warranties of merchantability and
76 fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
77 CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary,
78 or consequential damages, including but not limited to procurement of substitute
79 goods or services; loss of use, data, or profits, or business interruption,
80 however caused and on any theory of liability, whether in contract, strict
81 liability, or tort (including negligence), arising in any way out of the use of
82 this software, even if advised of the possibility of such damage.
83 
84 5.    CONTACT INFORMATION
85 
86 Fraunhofer Institute for Integrated Circuits IIS
87 Attention: Audio and Multimedia Departments - FDK AAC LL
88 Am Wolfsmantel 33
89 91058 Erlangen, Germany
90 
91 www.iis.fraunhofer.de/amm
92 amm-info@iis.fraunhofer.de
93 ----------------------------------------------------------------------------- */
94 
95 /**************************** AAC decoder library ******************************
96 
97    Author(s):   Matthias Hildenbrand, Manuel Jander
98 
99    Description: USAC LPC/AVQ decode
100 
101 *******************************************************************************/
102 
103 #include "usacdec_lpc.h"
104 
105 #include "usacdec_rom.h"
106 #include "FDK_trigFcts.h"
107 
108 #define NQ_MAX 36
109 
110 /*
111  * Helper functions.
112  */
113 
114 /**
115  * \brief Read unary code.
116  * \param hBs bitstream handle as data source.
117  * \return decoded value.
118  */
get_vlclbf(HANDLE_FDK_BITSTREAM hBs)119 static int get_vlclbf(HANDLE_FDK_BITSTREAM hBs) {
120   int result = 0;
121 
122   while (FDKreadBits(hBs, 1) && result <= NQ_MAX) {
123     result++;
124   }
125   return result;
126 }
127 
128 /**
129  * \brief Read bit count limited unary code.
130  * \param hBs bitstream handle as data source
131  * \param n max amount of bits to be read.
132  * \return decoded value.
133  */
get_vlclbf_n(HANDLE_FDK_BITSTREAM hBs,int n)134 static int get_vlclbf_n(HANDLE_FDK_BITSTREAM hBs, int n) {
135   int result = 0;
136 
137   while (FDKreadBits(hBs, 1)) {
138     result++;
139     n--;
140     if (n <= 0) {
141       break;
142     }
143   }
144 
145   return result;
146 }
147 
148 /*
149  * Algebraic Vector Quantizer
150  */
151 
152 /* ZF_SCALE must be greater than (number of FIXP_ZF)/2
153    because the loss of precision caused by fPow2Div2 in RE8_PPV() */
154 //#define ZF_SCALE ((NQ_MAX-3)>>1)
155 #define ZF_SCALE ((DFRACT_BITS / 2))
156 #define FIXP_ZF FIXP_DBL
157 #define INT2ZF(x, s) (FIXP_ZF)((x) << (ZF_SCALE - (s)))
158 #define ZF2INT(x) (INT)((x) >> ZF_SCALE)
159 
160 /* 1.0 in ZF format format */
161 #define ONEZF ((FIXP_ZF)INT2ZF(1, 0))
162 
163 /* static */
nearest_neighbor_2D8(FIXP_ZF x[8],int y[8])164 void nearest_neighbor_2D8(FIXP_ZF x[8], int y[8]) {
165   FIXP_ZF s, em, e[8];
166   int i, j, sum;
167 
168   /* round x into 2Z^8 i.e. compute y=(y1,...,y8) such that yi = 2[xi/2]
169      where [.] is the nearest integer operator
170      in the mean time, compute sum = y1+...+y8
171   */
172   sum = 0;
173   for (i = 0; i < 8; i++) {
174     FIXP_ZF tmp;
175     /* round to ..., -2, 0, 2, ... ([-1..1[ --> 0) */
176     if (x[i] < (FIXP_ZF)0) {
177       tmp = ONEZF - x[i];
178       y[i] = -2 * ((ZF2INT(tmp)) >> 1);
179     } else {
180       tmp = ONEZF + x[i];
181       y[i] = 2 * ((ZF2INT(tmp)) >> 1);
182     }
183     sum += y[i];
184   }
185   /* check if y1+...+y8 is a multiple of 4
186      if not, y is not round xj in the wrong way where j is defined by
187         j = arg max_i | xi -yi|
188      (this is called the Wagner rule)
189   */
190   if (sum % 4) {
191     /* find j = arg max_i | xi -yi| */
192     em = (FIXP_SGL)0;
193     j = 0;
194     for (i = 0; i < 8; i++) {
195       /* compute ei = xi-yi */
196       e[i] = x[i] - INT2ZF(y[i], 0);
197     }
198     for (i = 0; i < 8; i++) {
199       /* compute |ei| = | xi-yi | */
200       if (e[i] < (FIXP_ZF)0) {
201         s = -e[i];
202       } else {
203         s = e[i];
204       }
205       /* check if |ei| is maximal, if so, set j=i */
206       if (em < s) {
207         em = s;
208         j = i;
209       }
210     }
211     /* round xj in the "wrong way" */
212     if (e[j] < (FIXP_ZF)0) {
213       y[j] -= 2;
214     } else {
215       y[j] += 2;
216     }
217   }
218 }
219 
220 /*--------------------------------------------------------------
221   RE8_PPV(x,y)
222   NEAREST NEIGHBOR SEARCH IN INFINITE LATTICE RE8
223   the algorithm is based on the definition of RE8 as
224       RE8 = (2D8) U (2D8+[1,1,1,1,1,1,1,1])
225   it applies the coset decoding of Sloane and Conway
226   (i) x: point in R^8 in 32-ZF_SCALE.ZF_SCALE format
227   (o) y: point in RE8 (8-dimensional integer vector)
228   --------------------------------------------------------------
229 */
230 /* static */
RE8_PPV(FIXP_ZF x[],SHORT y[],int r)231 void RE8_PPV(FIXP_ZF x[], SHORT y[], int r) {
232   int i, y0[8], y1[8];
233   FIXP_ZF x1[8], tmp;
234   INT64 e;
235 
236   /* find the nearest neighbor y0 of x in 2D8 */
237   nearest_neighbor_2D8(x, y0);
238   /* find the nearest neighbor y1 of x in 2D8+(1,...,1) (by coset decoding) */
239   for (i = 0; i < 8; i++) {
240     x1[i] = x[i] - ONEZF;
241   }
242   nearest_neighbor_2D8(x1, y1);
243   for (i = 0; i < 8; i++) {
244     y1[i] += 1;
245   }
246 
247   /* compute e0=||x-y0||^2 and e1=||x-y1||^2 */
248   e = 0;
249   for (i = 0; i < 8; i++) {
250     tmp = x[i] - INT2ZF(y0[i], 0);
251     e += (INT64)fPow2Div2(
252         tmp << r); /* shift left to ensure that no fract part bits get lost. */
253     tmp = x[i] - INT2ZF(y1[i], 0);
254     e -= (INT64)fPow2Div2(tmp << r);
255   }
256   /* select best candidate y0 or y1 to minimize distortion */
257   if (e < 0) {
258     for (i = 0; i < 8; i++) {
259       y[i] = y0[i];
260     }
261   } else {
262     for (i = 0; i < 8; i++) {
263       y[i] = y1[i];
264     }
265   }
266 }
267 
268 /* table look-up of unsigned value: find i where index >= table[i]
269    Note: range must be >= 2, index must be >= table[0] */
table_lookup(const USHORT * table,unsigned int index,int range)270 static int table_lookup(const USHORT *table, unsigned int index, int range) {
271   int i;
272 
273   for (i = 4; i < range; i += 4) {
274     if (index < table[i]) {
275       break;
276     }
277   }
278   if (i > range) {
279     i = range;
280   }
281 
282   if (index < table[i - 2]) {
283     i -= 2;
284   }
285   if (index < table[i - 1]) {
286     i--;
287   }
288   i--;
289 
290   return (i); /* index >= table[i] */
291 }
292 
293 /*--------------------------------------------------------------------------
294   re8_decode_rank_of_permutation(rank, xs, x)
295   DECODING OF THE RANK OF THE PERMUTATION OF xs
296   (i) rank: index (rank) of a permutation
297   (i) xs:   signed leader in RE8 (8-dimensional integer vector)
298   (o) x:    point in RE8 (8-dimensional integer vector)
299   --------------------------------------------------------------------------
300  */
re8_decode_rank_of_permutation(int rank,int * xs,SHORT x[8])301 static void re8_decode_rank_of_permutation(int rank, int *xs, SHORT x[8]) {
302   INT a[8], w[8], B, fac, fac_B, target;
303   int i, j;
304 
305   /* --- pre-processing based on the signed leader xs ---
306      - compute the alphabet a=[a[0] ... a[q-1]] of x (q elements)
307        such that a[0]!=...!=a[q-1]
308        it is assumed that xs is sorted in the form of a signed leader
309        which can be summarized in 2 requirements:
310           a) |xs[0]| >= |xs[1]| >= |xs[2]| >= ... >= |xs[7]|
311           b) if |xs[i]|=|xs[i-1]|, xs[i]>=xs[i+1]
312        where |.| indicates the absolute value operator
313      - compute q (the number of symbols in the alphabet)
314      - compute w[0..q-1] where w[j] counts the number of occurences of
315        the symbol a[j] in xs
316      - compute B = prod_j=0..q-1 (w[j]!) where .! is the factorial */
317   /* xs[i], xs[i-1] and ptr_w/a*/
318   j = 0;
319   w[j] = 1;
320   a[j] = xs[0];
321   B = 1;
322   for (i = 1; i < 8; i++) {
323     if (xs[i] != xs[i - 1]) {
324       j++;
325       w[j] = 1;
326       a[j] = xs[i];
327     } else {
328       w[j]++;
329       B *= w[j];
330     }
331   }
332 
333   /* --- actual rank decoding ---
334      the rank of x (where x is a permutation of xs) is based on
335      Schalkwijk's formula
336      it is given by rank=sum_{k=0..7} (A_k * fac_k/B_k)
337      the decoding of this rank is sequential and reconstructs x[0..7]
338      element by element from x[0] to x[7]
339      [the tricky part is the inference of A_k for each k...]
340    */
341 
342   if (w[0] == 8) {
343     for (i = 0; i < 8; i++) {
344       x[i] = a[0]; /* avoid fac of 40320 */
345     }
346   } else {
347     target = rank * B;
348     fac_B = 1;
349     /* decode x element by element */
350     for (i = 0; i < 8; i++) {
351       fac = fac_B * fdk_dec_tab_factorial[i]; /* fac = 1..5040 */
352       j = -1;
353       do {
354         target -= w[++j] * fac;
355       } while (target >= 0); /* max of 30 tests / SV */
356       x[i] = a[j];
357       /* update rank, denominator B (B_k) and counter w[j] */
358       target += w[j] * fac; /* target = fac_B*B*rank */
359       fac_B *= w[j];
360       w[j]--;
361     }
362   }
363 }
364 
365 /*--------------------------------------------------------------------------
366   re8_decode_base_index(n, I, y)
367   DECODING OF AN INDEX IN Qn (n=0,2,3 or 4)
368   (i) n: codebook number (*n is an integer defined in {0,2,3,4})
369   (i) I: index of c (pointer to unsigned 16-bit word)
370   (o) y: point in RE8 (8-dimensional integer vector)
371   note: the index I is defined as a 32-bit word, but only
372   16 bits are required (long can be replaced by unsigned integer)
373   --------------------------------------------------------------------------
374  */
re8_decode_base_index(int * n,UINT index,SHORT y[8])375 static void re8_decode_base_index(int *n, UINT index, SHORT y[8]) {
376   int i, im, t, sign_code, ka, ks, rank, leader[8];
377 
378   if (*n < 2) {
379     for (i = 0; i < 8; i++) {
380       y[i] = 0;
381     }
382   } else {
383     // index = (unsigned int)*I;
384     /* search for the identifier ka of the absolute leader (table-lookup)
385        Q2 is a subset of Q3 - the two cases are considered in the same branch
386      */
387     switch (*n) {
388       case 2:
389       case 3:
390         i = table_lookup(fdk_dec_I3, index, NB_LDQ3);
391         ka = fdk_dec_A3[i];
392         break;
393       case 4:
394         i = table_lookup(fdk_dec_I4, index, NB_LDQ4);
395         ka = fdk_dec_A4[i];
396         break;
397       default:
398         FDK_ASSERT(0);
399         return;
400     }
401     /* reconstruct the absolute leader */
402     for (i = 0; i < 8; i++) {
403       leader[i] = fdk_dec_Da[ka][i];
404     }
405     /* search for the identifier ks of the signed leader (table look-up)
406        (this search is focused based on the identifier ka of the absolute
407         leader)*/
408     t = fdk_dec_Ia[ka];
409     im = fdk_dec_Ns[ka];
410     ks = table_lookup(fdk_dec_Is + t, index, im);
411 
412     /* reconstruct the signed leader from its sign code */
413     sign_code = 2 * fdk_dec_Ds[t + ks];
414     for (i = 7; i >= 0; i--) {
415       leader[i] *= (1 - (sign_code & 2));
416       sign_code >>= 1;
417     }
418 
419     /* compute and decode the rank of the permutation */
420     rank = index - fdk_dec_Is[t + ks]; /* rank = index - cardinality offset */
421 
422     re8_decode_rank_of_permutation(rank, leader, y);
423   }
424   return;
425 }
426 
427 /* re8_y2k(y,m,k)
428    VORONOI INDEXING (INDEX DECODING) k -> y
429    (i) k: Voronoi index k[0..7]
430    (i) m: Voronoi modulo (m = 2^r = 1<<r, where r is integer >=2)
431    (i) r: Voronoi order  (m = 2^r = 1<<r, where r is integer >=2)
432    (o) y: 8-dimensional point y[0..7] in RE8
433  */
re8_k2y(int * k,int r,SHORT * y)434 static void re8_k2y(int *k, int r, SHORT *y) {
435   int i, tmp, sum;
436   SHORT v[8];
437   FIXP_ZF zf[8];
438 
439   FDK_ASSERT(r <= ZF_SCALE);
440 
441   /* compute y = k M and z=(y-a)/m, where
442      M = [4        ]
443          [2 2      ]
444          [|   \    ]
445          [2     2  ]
446          [1 1 _ 1 1]
447      a=(2,0,...,0)
448      m = 1<<r
449   */
450   for (i = 0; i < 8; i++) {
451     y[i] = k[7];
452   }
453   zf[7] = INT2ZF(y[7], r);
454   sum = 0;
455   for (i = 6; i >= 1; i--) {
456     tmp = 2 * k[i];
457     sum += tmp;
458     y[i] += tmp;
459     zf[i] = INT2ZF(y[i], r);
460   }
461   y[0] += (4 * k[0] + sum);
462   zf[0] = INT2ZF(y[0] - 2, r);
463   /* find nearest neighbor v of z in infinite RE8 */
464   RE8_PPV(zf, v, r);
465   /* compute y -= m v */
466   for (i = 0; i < 8; i++) {
467     y[i] -= (SHORT)(v[i] << r);
468   }
469 }
470 
471 /*--------------------------------------------------------------------------
472   RE8_dec(n, I, k, y)
473   MULTI-RATE INDEXING OF A POINT y in THE LATTICE RE8 (INDEX DECODING)
474   (i) n: codebook number (*n is an integer defined in {0,2,3,4,..,n_max}). n_max
475   = 36 (i) I: index of c (pointer to unsigned 16-bit word) (i) k: index of v
476   (8-dimensional vector of binary indices) = Voronoi index (o) y: point in RE8
477   (8-dimensional integer vector) note: the index I is defined as a 32-bit word,
478   but only 16 bits are required (long can be replaced by unsigned integer)
479 
480   return 0 on success, -1 on error.
481   --------------------------------------------------------------------------
482  */
RE8_dec(int n,int I,int * k,FIXP_DBL * y)483 static int RE8_dec(int n, int I, int *k, FIXP_DBL *y) {
484   SHORT v[8];
485   SHORT _y[8];
486   UINT r;
487   int i;
488 
489   /* Check bound of codebook qn */
490   if (n > NQ_MAX) {
491     return -1;
492   }
493 
494   /* decode the sub-indices I and kv[] according to the codebook number n:
495      if n=0,2,3,4, decode I (no Voronoi extension)
496      if n>4, Voronoi extension is used, decode I and kv[] */
497   if (n <= 4) {
498     re8_decode_base_index(&n, I, _y);
499     for (i = 0; i < 8; i++) {
500       y[i] = (LONG)_y[i];
501     }
502   } else {
503     /* compute the Voronoi modulo m = 2^r where r is extension order */
504     r = ((n - 3) >> 1);
505 
506     while (n > 4) {
507       n -= 2;
508     }
509     /* decode base codebook index I into c (c is an element of Q3 or Q4)
510        [here c is stored in y to save memory] */
511     re8_decode_base_index(&n, I, _y);
512     /* decode Voronoi index k[] into v */
513     re8_k2y(k, r, v);
514     /* reconstruct y as y = m c + v (with m=2^r, r integer >=1) */
515     for (i = 0; i < 8; i++) {
516       y[i] = (LONG)((_y[i] << r) + v[i]);
517     }
518   }
519   return 0;
520 }
521 
522 /**************************/
523 /* start LPC decode stuff */
524 /**************************/
525 //#define M         16
526 #define FREQ_MAX 6400.0f
527 #define FREQ_DIV 400.0f
528 #define LSF_GAP 50.0f
529 
530 /**
531  * \brief calculate inverse weighting factor and add non-weighted residual
532  *        LSF vector to first stage LSF approximation
533  * \param lsfq first stage LSF approximation values.
534  * \param xq weighted residual LSF vector
535  * \param nk_mode code book number coding mode.
536  */
lsf_weight_2st(FIXP_LPC * lsfq,FIXP_DBL * xq,int nk_mode)537 static void lsf_weight_2st(FIXP_LPC *lsfq, FIXP_DBL *xq, int nk_mode) {
538   FIXP_LPC d[M_LP_FILTER_ORDER + 1];
539   FIXP_SGL factor;
540   LONG w; /* inverse weight factor */
541   int i;
542 
543   /* compute lsf distance */
544   d[0] = lsfq[0];
545   d[M_LP_FILTER_ORDER] =
546       FL2FXCONST_LPC(FREQ_MAX / (1 << LSF_SCALE)) - lsfq[M_LP_FILTER_ORDER - 1];
547   for (i = 1; i < M_LP_FILTER_ORDER; i++) {
548     d[i] = lsfq[i] - lsfq[i - 1];
549   }
550 
551   switch (nk_mode) {
552     case 0:
553       factor = FL2FXCONST_SGL(2.0f * 60.0f / FREQ_DIV);
554       break; /* abs */
555     case 1:
556       factor = FL2FXCONST_SGL(2.0f * 65.0f / FREQ_DIV);
557       break; /* mid */
558     case 2:
559       factor = FL2FXCONST_SGL(2.0f * 64.0f / FREQ_DIV);
560       break; /* rel1 */
561     default:
562       factor = FL2FXCONST_SGL(2.0f * 63.0f / FREQ_DIV);
563       break; /* rel2 */
564   }
565   /* add non-weighted residual LSF vector to LSF1st */
566   for (i = 0; i < M_LP_FILTER_ORDER; i++) {
567     w = (LONG)fMultDiv2(factor, sqrtFixp(fMult(d[i], d[i + 1])));
568     lsfq[i] = fAddSaturate(lsfq[i],
569                            FX_DBL2FX_LPC((FIXP_DBL)((INT64)w * (LONG)xq[i])));
570   }
571 
572   return;
573 }
574 
575 /**
576  * \brief decode nqn amount of code book numbers. These values determine the
577  * amount of following bits for nqn AVQ RE8 vectors.
578  * \param nk_mode quantization mode.
579  * \param nqn amount code book number to read.
580  * \param qn pointer to output buffer to hold decoded code book numbers qn.
581  */
decode_qn(HANDLE_FDK_BITSTREAM hBs,int nk_mode,int nqn,int qn[])582 static void decode_qn(HANDLE_FDK_BITSTREAM hBs, int nk_mode, int nqn,
583                       int qn[]) {
584   int n;
585 
586   if (nk_mode == 1) { /* nk mode 1 */
587     /* Unary code for mid LPC1/LPC3 */
588     /* Q0=0, Q2=10, Q3=110, ... */
589     for (n = 0; n < nqn; n++) {
590       qn[n] = get_vlclbf(hBs);
591       if (qn[n] > 0) {
592         qn[n]++;
593       }
594     }
595   } else { /* nk_mode 0, 3 and 2 */
596     /* 2 bits to specify Q2,Q3,Q4,ext */
597     for (n = 0; n < nqn; n++) {
598       qn[n] = 2 + FDKreadBits(hBs, 2);
599     }
600     if (nk_mode == 2) {
601       /* Unary code for rel LPC1/LPC3 */
602       /* Q0 = 0, Q5=10, Q6=110, ... */
603       for (n = 0; n < nqn; n++) {
604         if (qn[n] > 4) {
605           qn[n] = get_vlclbf(hBs);
606           if (qn[n] > 0) qn[n] += 4;
607         }
608       }
609     } else { /* nk_mode == (0 and 3) */
610       /* Unary code for abs and rel LPC0/LPC2 */
611       /* Q5 = 0, Q6=10, Q0=110, Q7=1110, ... */
612       for (n = 0; n < nqn; n++) {
613         if (qn[n] > 4) {
614           qn[n] = get_vlclbf(hBs);
615           switch (qn[n]) {
616             case 0:
617               qn[n] = 5;
618               break;
619             case 1:
620               qn[n] = 6;
621               break;
622             case 2:
623               qn[n] = 0;
624               break;
625             default:
626               qn[n] += 4;
627               break;
628           }
629         }
630       }
631     }
632   }
633 }
634 
635 /**
636  * \brief reorder LSF coefficients to minimum distance.
637  * \param lsf pointer to buffer containing LSF coefficients and where reordered
638  * LSF coefficients will be stored into, scaled by LSF_SCALE.
639  * \param min_dist min distance scaled by LSF_SCALE
640  * \param n number of LSF/LSP coefficients.
641  */
reorder_lsf(FIXP_LPC * lsf,FIXP_LPC min_dist,int n)642 static void reorder_lsf(FIXP_LPC *lsf, FIXP_LPC min_dist, int n) {
643   FIXP_LPC lsf_min;
644   int i;
645 
646   lsf_min = min_dist;
647   for (i = 0; i < n; i++) {
648     if (lsf[i] < lsf_min) {
649       lsf[i] = lsf_min;
650     }
651     lsf_min = fAddSaturate(lsf[i], min_dist);
652   }
653 
654   /* reverse */
655   lsf_min = FL2FXCONST_LPC(FREQ_MAX / (1 << LSF_SCALE)) - min_dist;
656   for (i = n - 1; i >= 0; i--) {
657     if (lsf[i] > lsf_min) {
658       lsf[i] = lsf_min;
659     }
660 
661     lsf_min = lsf[i] - min_dist;
662   }
663 }
664 
665 /**
666  * \brief First stage approximation
667  * \param hBs bitstream handle as data source
668  * \param lsfq pointer to output buffer to hold LPC coefficients scaled by
669  * LSF_SCALE.
670  */
vlpc_1st_dec(HANDLE_FDK_BITSTREAM hBs,FIXP_LPC * lsfq)671 static void vlpc_1st_dec(
672     HANDLE_FDK_BITSTREAM hBs, /* input:  codebook index                  */
673     FIXP_LPC *lsfq            /* i/o:    i:prediction   o:quantized lsf  */
674 ) {
675   const FIXP_LPC *p_dico;
676   int i, index;
677 
678   index = FDKreadBits(hBs, 8);
679   p_dico = &fdk_dec_dico_lsf_abs_8b[index * M_LP_FILTER_ORDER];
680   for (i = 0; i < M_LP_FILTER_ORDER; i++) {
681     lsfq[i] = p_dico[i];
682   }
683 }
684 
685 /**
686  * \brief Do first stage approximation weighting and multiply with AVQ
687  * refinement.
688  * \param hBs bitstream handle data ssource.
689  * \param lsfq buffer holding 1st stage approx, 2nd stage approx is added to
690  * this values.
691  * \param nk_mode quantization mode.
692  * \return 0 on success, -1 on error.
693  */
vlpc_2st_dec(HANDLE_FDK_BITSTREAM hBs,FIXP_LPC * lsfq,int nk_mode)694 static int vlpc_2st_dec(
695     HANDLE_FDK_BITSTREAM hBs,
696     FIXP_LPC *lsfq, /* i/o:    i:1st stage   o:1st+2nd stage   */
697     int nk_mode     /* input:  0=abs, >0=rel                   */
698 ) {
699   int err;
700   FIXP_DBL xq[M_LP_FILTER_ORDER]; /* weighted residual LSF vector */
701 
702   /* Decode AVQ refinement */
703   { err = CLpc_DecodeAVQ(hBs, xq, nk_mode, 2, 8); }
704   if (err != 0) {
705     return -1;
706   }
707 
708   /* add non-weighted residual LSF vector to LSF1st */
709   lsf_weight_2st(lsfq, xq, nk_mode);
710 
711   /* reorder */
712   reorder_lsf(lsfq, FL2FXCONST_LPC(LSF_GAP / (1 << LSF_SCALE)),
713               M_LP_FILTER_ORDER);
714 
715   return 0;
716 }
717 
718 /*
719  * Externally visible functions
720  */
721 
CLpc_DecodeAVQ(HANDLE_FDK_BITSTREAM hBs,FIXP_DBL * pOutput,int nk_mode,int no_qn,int length)722 int CLpc_DecodeAVQ(HANDLE_FDK_BITSTREAM hBs, FIXP_DBL *pOutput, int nk_mode,
723                    int no_qn, int length) {
724   int i, l;
725 
726   for (i = 0; i < length; i += 8 * no_qn) {
727     int qn[2], nk, n, I;
728     int kv[8] = {0};
729 
730     decode_qn(hBs, nk_mode, no_qn, qn);
731 
732     for (l = 0; l < no_qn; l++) {
733       if (qn[l] == 0) {
734         FDKmemclear(&pOutput[i + l * 8], 8 * sizeof(FIXP_DBL));
735       }
736 
737       /* Voronoi extension order ( nk ) */
738       nk = 0;
739       n = qn[l];
740       if (qn[l] > 4) {
741         nk = (qn[l] - 3) >> 1;
742         n = qn[l] - nk * 2;
743       }
744 
745       /* Base codebook index, in reverse bit group order (!) */
746       I = FDKreadBits(hBs, 4 * n);
747 
748       if (nk > 0) {
749         int j;
750 
751         for (j = 0; j < 8; j++) {
752           kv[j] = FDKreadBits(hBs, nk);
753         }
754       }
755 
756       if (RE8_dec(qn[l], I, kv, &pOutput[i + l * 8]) != 0) {
757         return -1;
758       }
759     }
760   }
761   return 0;
762 }
763 
CLpc_Read(HANDLE_FDK_BITSTREAM hBs,FIXP_LPC lsp[][M_LP_FILTER_ORDER],FIXP_LPC lpc4_lsf[M_LP_FILTER_ORDER],FIXP_LPC lsf_adaptive_mean_cand[M_LP_FILTER_ORDER],FIXP_SGL pStability[],UCHAR * mod,int first_lpd_flag,int last_lpc_lost,int last_frame_ok)764 int CLpc_Read(HANDLE_FDK_BITSTREAM hBs, FIXP_LPC lsp[][M_LP_FILTER_ORDER],
765               FIXP_LPC lpc4_lsf[M_LP_FILTER_ORDER],
766               FIXP_LPC lsf_adaptive_mean_cand[M_LP_FILTER_ORDER],
767               FIXP_SGL pStability[], UCHAR *mod, int first_lpd_flag,
768               int last_lpc_lost, int last_frame_ok) {
769   int i, k, err;
770   int mode_lpc_bin = 0; /* mode_lpc bitstream representation */
771   int lpc_present[5] = {0, 0, 0, 0, 0};
772   int lpc0_available = 1;
773   int s = 0;
774   int l = 3;
775   const int nbDiv = NB_DIV;
776 
777   lpc_present[4 >> s] = 1; /* LPC4 */
778 
779   /* Decode LPC filters in the following order: LPC 4,0,2,1,3 */
780 
781   /*** Decode LPC4 ***/
782   vlpc_1st_dec(hBs, lsp[4 >> s]);
783   err = vlpc_2st_dec(hBs, lsp[4 >> s], 0); /* nk_mode = 0 */
784   if (err != 0) {
785     return err;
786   }
787 
788   /*** Decode LPC0 and LPC2 ***/
789   k = 0;
790   if (!first_lpd_flag) {
791     lpc_present[0] = 1;
792     lpc0_available = !last_lpc_lost;
793     /* old LPC4 is new LPC0 */
794     for (i = 0; i < M_LP_FILTER_ORDER; i++) {
795       lsp[0][i] = lpc4_lsf[i];
796     }
797     /* skip LPC0 and continue with LPC2 */
798     k = 2;
799   }
800 
801   for (; k < l; k += 2) {
802     int nk_mode = 0;
803 
804     if ((k == 2) && (mod[0] == 3)) {
805       break; /* skip LPC2 */
806     }
807 
808     lpc_present[k >> s] = 1;
809 
810     mode_lpc_bin = FDKreadBit(hBs);
811 
812     if (mode_lpc_bin == 0) {
813       /* LPC0/LPC2: Abs */
814       vlpc_1st_dec(hBs, lsp[k >> s]);
815     } else {
816       /* LPC0/LPC2: RelR */
817       for (i = 0; i < M_LP_FILTER_ORDER; i++) {
818         lsp[k >> s][i] = lsp[4 >> s][i];
819       }
820       nk_mode = 3;
821     }
822 
823     err = vlpc_2st_dec(hBs, lsp[k >> s], nk_mode);
824     if (err != 0) {
825       return err;
826     }
827   }
828 
829   /*** Decode LPC1 ***/
830   if (mod[0] < 2) { /* else: skip LPC1 */
831     lpc_present[1] = 1;
832     mode_lpc_bin = get_vlclbf_n(hBs, 2);
833 
834     switch (mode_lpc_bin) {
835       case 1:
836         /* LPC1: abs */
837         vlpc_1st_dec(hBs, lsp[1]);
838         err = vlpc_2st_dec(hBs, lsp[1], 0);
839         if (err != 0) {
840           return err;
841         }
842         break;
843       case 2:
844         /* LPC1: mid0 (no second stage AVQ quantizer in this case) */
845         if (lpc0_available) { /* LPC0/lsf[0] might be zero some times */
846           for (i = 0; i < M_LP_FILTER_ORDER; i++) {
847             lsp[1][i] = (lsp[0][i] >> 1) + (lsp[2][i] >> 1);
848           }
849         } else {
850           for (i = 0; i < M_LP_FILTER_ORDER; i++) {
851             lsp[1][i] = lsp[2][i];
852           }
853         }
854         break;
855       case 0:
856         /* LPC1: RelR */
857         for (i = 0; i < M_LP_FILTER_ORDER; i++) {
858           lsp[1][i] = lsp[2][i];
859         }
860         err = vlpc_2st_dec(hBs, lsp[1], 2 << s);
861         if (err != 0) {
862           return err;
863         }
864         break;
865     }
866   }
867 
868   /*** Decode LPC3 ***/
869   if ((mod[2] < 2)) { /* else: skip LPC3 */
870     int nk_mode = 0;
871     lpc_present[3] = 1;
872 
873     mode_lpc_bin = get_vlclbf_n(hBs, 3);
874 
875     switch (mode_lpc_bin) {
876       case 1:
877         /* LPC3: abs */
878         vlpc_1st_dec(hBs, lsp[3]);
879         break;
880       case 0:
881         /* LPC3: mid */
882         for (i = 0; i < M_LP_FILTER_ORDER; i++) {
883           lsp[3][i] = (lsp[2][i] >> 1) + (lsp[4][i] >> 1);
884         }
885         nk_mode = 1;
886         break;
887       case 2:
888         /* LPC3: relL */
889         for (i = 0; i < M_LP_FILTER_ORDER; i++) {
890           lsp[3][i] = lsp[2][i];
891         }
892         nk_mode = 2;
893         break;
894       case 3:
895         /* LPC3: relR */
896         for (i = 0; i < M_LP_FILTER_ORDER; i++) {
897           lsp[3][i] = lsp[4][i];
898         }
899         nk_mode = 2;
900         break;
901     }
902     err = vlpc_2st_dec(hBs, lsp[3], nk_mode);
903     if (err != 0) {
904       return err;
905     }
906   }
907 
908   if (!lpc0_available && !last_frame_ok) {
909     /* LPC(0) was lost. Use next available LPC(k) instead */
910     for (k = 1; k < (nbDiv + 1); k++) {
911       if (lpc_present[k]) {
912         for (i = 0; i < M_LP_FILTER_ORDER; i++) {
913 #define LSF_INIT_TILT (0.25f)
914           if (mod[0] > 0) {
915             lsp[0][i] = FX_DBL2FX_LPC(
916                 fMult(lsp[k][i], FL2FXCONST_SGL(1.0f - LSF_INIT_TILT)) +
917                 fMult(fdk_dec_lsf_init[i], FL2FXCONST_SGL(LSF_INIT_TILT)));
918           } else {
919             lsp[0][i] = lsp[k][i];
920           }
921         }
922         break;
923       }
924     }
925   }
926 
927   for (i = 0; i < M_LP_FILTER_ORDER; i++) {
928     lpc4_lsf[i] = lsp[4 >> s][i];
929   }
930 
931   {
932     FIXP_DBL divFac;
933     int last, numLpc = 0;
934 
935     i = nbDiv;
936     do {
937       numLpc += lpc_present[i--];
938     } while (i >= 0 && numLpc < 3);
939 
940     last = i;
941 
942     switch (numLpc) {
943       case 3:
944         divFac = FL2FXCONST_DBL(1.0f / 3.0f);
945         break;
946       case 2:
947         divFac = FL2FXCONST_DBL(1.0f / 2.0f);
948         break;
949       default:
950         divFac = FL2FXCONST_DBL(1.0f);
951         break;
952     }
953 
954     /* get the adaptive mean for the next (bad) frame */
955     for (k = 0; k < M_LP_FILTER_ORDER; k++) {
956       FIXP_DBL tmp = (FIXP_DBL)0;
957       for (i = nbDiv; i > last; i--) {
958         if (lpc_present[i]) {
959           tmp = fMultAdd(tmp >> 1, lsp[i][k], divFac);
960         }
961       }
962       lsf_adaptive_mean_cand[k] = FX_DBL2FX_LPC(tmp);
963     }
964   }
965 
966   /* calculate stability factor Theta. Needed for ACELP decoder and concealment
967    */
968   {
969     FIXP_LPC *lsf_prev, *lsf_curr;
970     k = 0;
971 
972     FDK_ASSERT(lpc_present[0] == 1 && lpc_present[4 >> s] == 1);
973     lsf_prev = lsp[0];
974     for (i = 1; i < (nbDiv + 1); i++) {
975       if (lpc_present[i]) {
976         FIXP_DBL tmp = (FIXP_DBL)0;
977         int j;
978         lsf_curr = lsp[i];
979 
980         /* sum = tmp * 2^(LSF_SCALE*2 + 4) */
981         for (j = 0; j < M_LP_FILTER_ORDER; j++) {
982           tmp += fPow2Div2((FIXP_SGL)(lsf_curr[j] - lsf_prev[j])) >> 3;
983         }
984 
985         /* tmp = (float)(FL2FXCONST_DBL(1.25f) - fMult(tmp,
986          * FL2FXCONST_DBL(1/400000.0f))); */
987         tmp = FL2FXCONST_DBL(1.25f / (1 << LSF_SCALE)) -
988               fMult(tmp, FL2FXCONST_DBL((1 << (LSF_SCALE + 4)) / 400000.0f));
989         if (tmp >= FL2FXCONST_DBL(1.0f / (1 << LSF_SCALE))) {
990           pStability[k] = FL2FXCONST_SGL(1.0f / 2.0f);
991         } else if (tmp < FL2FXCONST_DBL(0.0f)) {
992           pStability[k] = FL2FXCONST_SGL(0.0f);
993         } else {
994           pStability[k] = FX_DBL2FX_SGL(tmp << (LSF_SCALE - 1));
995         }
996 
997         lsf_prev = lsf_curr;
998         k = i;
999       } else {
1000         /* Mark stability value as undefined. */
1001         pStability[i] = (FIXP_SGL)-1;
1002       }
1003     }
1004   }
1005 
1006   /* convert into LSP domain */
1007   for (i = 0; i < (nbDiv + 1); i++) {
1008     if (lpc_present[i]) {
1009       for (k = 0; k < M_LP_FILTER_ORDER; k++) {
1010         lsp[i][k] = FX_DBL2FX_LPC(
1011             fixp_cos(fMult(lsp[i][k],
1012                            FL2FXCONST_SGL((1 << LSPARG_SCALE) * M_PI / 6400.0)),
1013                      LSF_SCALE - LSPARG_SCALE));
1014       }
1015     }
1016   }
1017 
1018   return 0;
1019 }
1020 
CLpc_Conceal(FIXP_LPC lsp[][M_LP_FILTER_ORDER],FIXP_LPC lpc4_lsf[M_LP_FILTER_ORDER],FIXP_LPC lsf_adaptive_mean[M_LP_FILTER_ORDER],const int first_lpd_flag)1021 void CLpc_Conceal(FIXP_LPC lsp[][M_LP_FILTER_ORDER],
1022                   FIXP_LPC lpc4_lsf[M_LP_FILTER_ORDER],
1023                   FIXP_LPC lsf_adaptive_mean[M_LP_FILTER_ORDER],
1024                   const int first_lpd_flag) {
1025   int i, j;
1026 
1027 #define BETA (FL2FXCONST_SGL(0.25f))
1028 #define ONE_BETA (FL2FXCONST_SGL(0.75f))
1029 #define BFI_FAC (FL2FXCONST_SGL(0.90f))
1030 #define ONE_BFI_FAC (FL2FXCONST_SGL(0.10f))
1031 
1032   /* Frame loss concealment (could be improved) */
1033 
1034   if (first_lpd_flag) {
1035     /* Reset past LSF values */
1036     for (i = 0; i < M_LP_FILTER_ORDER; i++) {
1037       lsp[0][i] = lpc4_lsf[i] = fdk_dec_lsf_init[i];
1038     }
1039   } else {
1040     /* old LPC4 is new LPC0 */
1041     for (i = 0; i < M_LP_FILTER_ORDER; i++) {
1042       lsp[0][i] = lpc4_lsf[i];
1043     }
1044   }
1045 
1046   /* LPC1 */
1047   for (i = 0; i < M_LP_FILTER_ORDER; i++) {
1048     FIXP_LPC lsf_mean = FX_DBL2FX_LPC(fMult(BETA, fdk_dec_lsf_init[i]) +
1049                                       fMult(ONE_BETA, lsf_adaptive_mean[i]));
1050 
1051     lsp[1][i] = FX_DBL2FX_LPC(fMult(BFI_FAC, lpc4_lsf[i]) +
1052                               fMult(ONE_BFI_FAC, lsf_mean));
1053   }
1054 
1055   /* LPC2 - LPC4 */
1056   for (j = 2; j <= 4; j++) {
1057     for (i = 0; i < M_LP_FILTER_ORDER; i++) {
1058       /* lsf_mean[i] =  FX_DBL2FX_LPC(fMult((FIXP_LPC)(BETA + j *
1059          FL2FXCONST_LPC(0.1f)), fdk_dec_lsf_init[i])
1060                                     + fMult((FIXP_LPC)(ONE_BETA - j *
1061          FL2FXCONST_LPC(0.1f)), lsf_adaptive_mean[i])); */
1062 
1063       FIXP_LPC lsf_mean = FX_DBL2FX_LPC(
1064           fMult((FIXP_SGL)(BETA + (FIXP_SGL)(j * (INT)FL2FXCONST_SGL(0.1f))),
1065                 (FIXP_SGL)fdk_dec_lsf_init[i]) +
1066           fMult(
1067               (FIXP_SGL)(ONE_BETA - (FIXP_SGL)(j * (INT)FL2FXCONST_SGL(0.1f))),
1068               lsf_adaptive_mean[i]));
1069 
1070       lsp[j][i] = FX_DBL2FX_LPC(fMult(BFI_FAC, lsp[j - 1][i]) +
1071                                 fMult(ONE_BFI_FAC, lsf_mean));
1072     }
1073   }
1074 
1075   /* Update past values for the future */
1076   for (i = 0; i < M_LP_FILTER_ORDER; i++) {
1077     lpc4_lsf[i] = lsp[4][i];
1078   }
1079 
1080   /* convert into LSP domain */
1081   for (j = 0; j < 5; j++) {
1082     for (i = 0; i < M_LP_FILTER_ORDER; i++) {
1083       lsp[j][i] = FX_DBL2FX_LPC(fixp_cos(
1084           fMult(lsp[j][i], FL2FXCONST_SGL((1 << LSPARG_SCALE) * M_PI / 6400.0)),
1085           LSF_SCALE - LSPARG_SCALE));
1086     }
1087   }
1088 }
1089 
E_LPC_a_weight(FIXP_LPC * wA,const FIXP_LPC * A,int m)1090 void E_LPC_a_weight(FIXP_LPC *wA, const FIXP_LPC *A, int m) {
1091   FIXP_DBL f;
1092   int i;
1093 
1094   f = FL2FXCONST_DBL(0.92f);
1095   for (i = 0; i < m; i++) {
1096     wA[i] = FX_DBL2FX_LPC(fMult(A[i], f));
1097     f = fMult(f, FL2FXCONST_DBL(0.92f));
1098   }
1099 }
1100 
CLpd_DecodeGain(FIXP_DBL * gain,INT * gain_e,int gain_code)1101 void CLpd_DecodeGain(FIXP_DBL *gain, INT *gain_e, int gain_code) {
1102   /* gain * 2^(gain_e) = 10^(gain_code/28) */
1103   *gain = fLdPow(
1104       FL2FXCONST_DBL(3.3219280948873623478703194294894 / 4.0), /* log2(10)*/
1105       2,
1106       fMultDiv2((FIXP_DBL)gain_code << (DFRACT_BITS - 1 - 7),
1107                 FL2FXCONST_DBL(2.0f / 28.0f)),
1108       7, gain_e);
1109 }
1110 
1111   /**
1112    * \brief *   Find the polynomial F1(z) or F2(z) from the LSPs.
1113    * This is performed by expanding the product polynomials:
1114    *
1115    * F1(z) =   product   ( 1 - 2 LSP_i z^-1 + z^-2 )
1116    *         i=0,2,4,6,8
1117    * F2(z) =   product   ( 1 - 2 LSP_i z^-1 + z^-2 )
1118    *         i=1,3,5,7,9
1119    *
1120    * where LSP_i are the LSPs in the cosine domain.
1121    * R.A.Salami    October 1990
1122    * \param lsp input, line spectral freq. (cosine domain)
1123    * \param f output, the coefficients of F1 or F2, scaled by 8 bits
1124    * \param n no of coefficients (m/2)
1125    * \param flag 1 : F1(z) ; 2 : F2(z)
1126    */
1127 
1128 #define SF_F 8
1129 
get_lsppol(FIXP_LPC lsp[],FIXP_DBL f[],int n,int flag)1130 static void get_lsppol(FIXP_LPC lsp[], FIXP_DBL f[], int n, int flag) {
1131   FIXP_DBL b;
1132   FIXP_LPC *plsp;
1133   int i, j;
1134 
1135   plsp = lsp + flag - 1;
1136   f[0] = FL2FXCONST_DBL(1.0f / (1 << SF_F));
1137   b = -FX_LPC2FX_DBL(*plsp);
1138   f[1] = b >> (SF_F - 1);
1139   for (i = 2; i <= n; i++) {
1140     plsp += 2;
1141     b = -FX_LPC2FX_DBL(*plsp);
1142     f[i] = SATURATE_LEFT_SHIFT((fMultDiv2(b, f[i - 1]) + (f[i - 2] >> 1)), 2,
1143                                DFRACT_BITS);
1144     for (j = i - 1; j > 1; j--) {
1145       f[j] = SATURATE_LEFT_SHIFT(
1146           ((f[j] >> 2) + fMultDiv2(b, f[j - 1]) + (f[j - 2] >> 2)), 2,
1147           DFRACT_BITS);
1148     }
1149     f[1] = f[1] + (b >> (SF_F - 1));
1150   }
1151   return;
1152 }
1153 
1154 #define NC M_LP_FILTER_ORDER / 2
1155 
1156 /**
1157  * \brief lsp input LSP vector
1158  * \brief a output LP filter coefficient vector scaled by SF_A_COEFFS.
1159  */
E_LPC_f_lsp_a_conversion(FIXP_LPC * lsp,FIXP_LPC * a,INT * a_exp)1160 void E_LPC_f_lsp_a_conversion(FIXP_LPC *lsp, FIXP_LPC *a, INT *a_exp) {
1161   FIXP_DBL f1[NC + 1], f2[NC + 1];
1162   int i, k;
1163 
1164   /*-----------------------------------------------------*
1165    *  Find the polynomials F1(z) and F2(z)               *
1166    *-----------------------------------------------------*/
1167 
1168   get_lsppol(lsp, f1, NC, 1);
1169   get_lsppol(lsp, f2, NC, 2);
1170 
1171   /*-----------------------------------------------------*
1172    *  Multiply F1(z) by (1+z^-1) and F2(z) by (1-z^-1)   *
1173    *-----------------------------------------------------*/
1174   scaleValues(f1, NC + 1, -2);
1175   scaleValues(f2, NC + 1, -2);
1176 
1177   for (i = NC; i > 0; i--) {
1178     f1[i] += f1[i - 1];
1179     f2[i] -= f2[i - 1];
1180   }
1181 
1182   FIXP_DBL aDBL[M_LP_FILTER_ORDER];
1183 
1184   for (i = 1, k = M_LP_FILTER_ORDER - 1; i <= NC; i++, k--) {
1185     aDBL[i - 1] = f1[i] + f2[i];
1186     aDBL[k] = f1[i] - f2[i];
1187   }
1188 
1189   int headroom_a = getScalefactor(aDBL, M_LP_FILTER_ORDER);
1190 
1191   for (i = 0; i < M_LP_FILTER_ORDER; i++) {
1192     a[i] = FX_DBL2FX_LPC(aDBL[i] << headroom_a);
1193   }
1194 
1195   *a_exp = SF_F + (2 - 1) - headroom_a;
1196 }
1197