Include patched FDK-AAC in the repository

The initial idea was to get the DAB+ patch into upstream, but since
that follows the android source releases, there is no place for a custom
DAB+ patch there.

So instead of having to maintain a patched fdk-aac that has to have the
same .so version as the distribution package on which it is installed,
we prefer having a separate fdk-aac-dab library to avoid collision.

At that point, there's no reason to keep fdk-aac in a separate
repository, as odr-audioenc is the only tool that needs DAB+ encoding
support. Including it here simplifies installation, and makes it
consistent with toolame-dab, also shipped in this repository.

DAB+ decoding support (needed by ODR-SourceCompanion, dablin, etisnoop,
welle.io and others) can be done using upstream FDK-AAC.

1 files changed, 1194 insertions, 0 deletions

diff --git a/fdk-aac/libAACdec/src/usacdec_lpc.cpp b/fdk-aac/libAACdec/src/usacdec_lpc.cpp
new file mode 100644
index 0000000..271463f
--- /dev/null
+++ b/fdk-aac/libAACdec/src/usacdec_lpc.cpp
@@ -0,0 +1,1194 @@
+/* -----------------------------------------------------------------------------
+Software License for The Fraunhofer FDK AAC Codec Library for Android
+
+© Copyright  1995 - 2018 Fraunhofer-Gesellschaft zur Förderung der angewandten
+Forschung e.V. All rights reserved.
+
+ 1.    INTRODUCTION
+The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software
+that implements the MPEG Advanced Audio Coding ("AAC") encoding and decoding
+scheme for digital audio. This FDK AAC Codec software is intended to be used on
+a wide variety of Android devices.
+
+AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient
+general perceptual audio codecs. AAC-ELD is considered the best-performing
+full-bandwidth communications codec by independent studies and is widely
+deployed. AAC has been standardized by ISO and IEC as part of the MPEG
+specifications.
+
+Patent licenses for necessary patent claims for the FDK AAC Codec (including
+those of Fraunhofer) may be obtained through Via Licensing
+(www.vialicensing.com) or through the respective patent owners individually for
+the purpose of encoding or decoding bit streams in products that are compliant
+with the ISO/IEC MPEG audio standards. Please note that most manufacturers of
+Android devices already license these patent claims through Via Licensing or
+directly from the patent owners, and therefore FDK AAC Codec software may
+already be covered under those patent licenses when it is used for those
+licensed purposes only.
+
+Commercially-licensed AAC software libraries, including floating-point versions
+with enhanced sound quality, are also available from Fraunhofer. Users are
+encouraged to check the Fraunhofer website for additional applications
+information and documentation.
+
+2.    COPYRIGHT LICENSE
+
+Redistribution and use in source and binary forms, with or without modification,
+are permitted without payment of copyright license fees provided that you
+satisfy the following conditions:
+
+You must retain the complete text of this software license in redistributions of
+the FDK AAC Codec or your modifications thereto in source code form.
+
+You must retain the complete text of this software license in the documentation
+and/or other materials provided with redistributions of the FDK AAC Codec or
+your modifications thereto in binary form. You must make available free of
+charge copies of the complete source code of the FDK AAC Codec and your
+modifications thereto to recipients of copies in binary form.
+
+The name of Fraunhofer may not be used to endorse or promote products derived
+from this library without prior written permission.
+
+You may not charge copyright license fees for anyone to use, copy or distribute
+the FDK AAC Codec software or your modifications thereto.
+
+Your modified versions of the FDK AAC Codec must carry prominent notices stating
+that you changed the software and the date of any change. For modified versions
+of the FDK AAC Codec, the term "Fraunhofer FDK AAC Codec Library for Android"
+must be replaced by the term "Third-Party Modified Version of the Fraunhofer FDK
+AAC Codec Library for Android."
+
+3.    NO PATENT LICENSE
+
+NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without
+limitation the patents of Fraunhofer, ARE GRANTED BY THIS SOFTWARE LICENSE.
+Fraunhofer provides no warranty of patent non-infringement with respect to this
+software.
+
+You may use this FDK AAC Codec software or modifications thereto only for
+purposes that are authorized by appropriate patent licenses.
+
+4.    DISCLAIMER
+
+This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright
+holders and contributors "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES,
+including but not limited to the implied warranties of merchantability and
+fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
+CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary,
+or consequential damages, including but not limited to procurement of substitute
+goods or services; loss of use, data, or profits, or business interruption,
+however caused and on any theory of liability, whether in contract, strict
+liability, or tort (including negligence), arising in any way out of the use of
+this software, even if advised of the possibility of such damage.
+
+5.    CONTACT INFORMATION
+
+Fraunhofer Institute for Integrated Circuits IIS
+Attention: Audio and Multimedia Departments - FDK AAC LL
+Am Wolfsmantel 33
+91058 Erlangen, Germany
+
+www.iis.fraunhofer.de/amm
+amm-info@iis.fraunhofer.de
+----------------------------------------------------------------------------- */
+
+/**************************** AAC decoder library ******************************
+
+   Author(s):   Matthias Hildenbrand, Manuel Jander
+
+   Description: USAC LPC/AVQ decode
+
+*******************************************************************************/
+
+#include "usacdec_lpc.h"
+
+#include "usacdec_rom.h"
+#include "FDK_trigFcts.h"
+
+#define NQ_MAX 36
+
+/*
+ * Helper functions.
+ */
+
+/**
+ * \brief Read unary code.
+ * \param hBs bitstream handle as data source.
+ * \return decoded value.
+ */
+static int get_vlclbf(HANDLE_FDK_BITSTREAM hBs) {
+  int result = 0;
+
+  while (FDKreadBits(hBs, 1) && result <= NQ_MAX) {
+    result++;
+  }
+  return result;
+}
+
+/**
+ * \brief Read bit count limited unary code.
+ * \param hBs bitstream handle as data source
+ * \param n max amount of bits to be read.
+ * \return decoded value.
+ */
+static int get_vlclbf_n(HANDLE_FDK_BITSTREAM hBs, int n) {
+  int result = 0;
+
+  while (FDKreadBits(hBs, 1)) {
+    result++;
+    n--;
+    if (n <= 0) {
+      break;
+    }
+  }
+
+  return result;
+}
+
+/*
+ * Algebraic Vector Quantizer
+ */
+
+/* ZF_SCALE must be greater than (number of FIXP_ZF)/2
+   because the loss of precision caused by fPow2Div2 in RE8_PPV() */
+//#define ZF_SCALE ((NQ_MAX-3)>>1)
+#define ZF_SCALE ((DFRACT_BITS / 2))
+#define FIXP_ZF FIXP_DBL
+#define INT2ZF(x, s) (FIXP_ZF)((x) << (ZF_SCALE - (s)))
+#define ZF2INT(x) (INT)((x) >> ZF_SCALE)
+
+/* 1.0 in ZF format format */
+#define ONEZF ((FIXP_ZF)INT2ZF(1, 0))
+
+/* static */
+void nearest_neighbor_2D8(FIXP_ZF x[8], int y[8]) {
+  FIXP_ZF s, em, e[8];
+  int i, j, sum;
+
+  /* round x into 2Z^8 i.e. compute y=(y1,...,y8) such that yi = 2[xi/2]
+     where [.] is the nearest integer operator
+     in the mean time, compute sum = y1+...+y8
+  */
+  sum = 0;
+  for (i = 0; i < 8; i++) {
+    FIXP_ZF tmp;
+    /* round to ..., -2, 0, 2, ... ([-1..1[ --> 0) */
+    if (x[i] < (FIXP_ZF)0) {
+      tmp = ONEZF - x[i];
+      y[i] = -2 * ((ZF2INT(tmp)) >> 1);
+    } else {
+      tmp = ONEZF + x[i];
+      y[i] = 2 * ((ZF2INT(tmp)) >> 1);
+    }
+    sum += y[i];
+  }
+  /* check if y1+...+y8 is a multiple of 4
+     if not, y is not round xj in the wrong way where j is defined by
+        j = arg max_i | xi -yi|
+     (this is called the Wagner rule)
+  */
+  if (sum % 4) {
+    /* find j = arg max_i | xi -yi| */
+    em = (FIXP_SGL)0;
+    j = 0;
+    for (i = 0; i < 8; i++) {
+      /* compute ei = xi-yi */
+      e[i] = x[i] - INT2ZF(y[i], 0);
+    }
+    for (i = 0; i < 8; i++) {
+      /* compute |ei| = | xi-yi | */
+      if (e[i] < (FIXP_ZF)0) {
+        s = -e[i];
+      } else {
+        s = e[i];
+      }
+      /* check if |ei| is maximal, if so, set j=i */
+      if (em < s) {
+        em = s;
+        j = i;
+      }
+    }
+    /* round xj in the "wrong way" */
+    if (e[j] < (FIXP_ZF)0) {
+      y[j] -= 2;
+    } else {
+      y[j] += 2;
+    }
+  }
+}
+
+/*--------------------------------------------------------------
+  RE8_PPV(x,y)
+  NEAREST NEIGHBOR SEARCH IN INFINITE LATTICE RE8
+  the algorithm is based on the definition of RE8 as
+      RE8 = (2D8) U (2D8+[1,1,1,1,1,1,1,1])
+  it applies the coset decoding of Sloane and Conway
+  (i) x: point in R^8 in 32-ZF_SCALE.ZF_SCALE format
+  (o) y: point in RE8 (8-dimensional integer vector)
+  --------------------------------------------------------------
+*/
+/* static */
+void RE8_PPV(FIXP_ZF x[], SHORT y[], int r) {
+  int i, y0[8], y1[8];
+  FIXP_ZF x1[8], tmp;
+  FIXP_DBL e;
+
+  /* find the nearest neighbor y0 of x in 2D8 */
+  nearest_neighbor_2D8(x, y0);
+  /* find the nearest neighbor y1 of x in 2D8+(1,...,1) (by coset decoding) */
+  for (i = 0; i < 8; i++) {
+    x1[i] = x[i] - ONEZF;
+  }
+  nearest_neighbor_2D8(x1, y1);
+  for (i = 0; i < 8; i++) {
+    y1[i] += 1;
+  }
+
+  /* compute e0=||x-y0||^2 and e1=||x-y1||^2 */
+  e = (FIXP_DBL)0;
+  for (i = 0; i < 8; i++) {
+    tmp = x[i] - INT2ZF(y0[i], 0);
+    e += fPow2Div2(
+        tmp << r); /* shift left to ensure that no fract part bits get lost. */
+    tmp = x[i] - INT2ZF(y1[i], 0);
+    e -= fPow2Div2(tmp << r);
+  }
+  /* select best candidate y0 or y1 to minimize distortion */
+  if (e < (FIXP_DBL)0) {
+    for (i = 0; i < 8; i++) {
+      y[i] = y0[i];
+    }
+  } else {
+    for (i = 0; i < 8; i++) {
+      y[i] = y1[i];
+    }
+  }
+}
+
+/* table look-up of unsigned value: find i where index >= table[i]
+   Note: range must be >= 2, index must be >= table[0] */
+static int table_lookup(const USHORT *table, unsigned int index, int range) {
+  int i;
+
+  for (i = 4; i < range; i += 4) {
+    if (index < table[i]) {
+      break;
+    }
+  }
+  if (i > range) {
+    i = range;
+  }
+
+  if (index < table[i - 2]) {
+    i -= 2;
+  }
+  if (index < table[i - 1]) {
+    i--;
+  }
+  i--;
+
+  return (i); /* index >= table[i] */
+}
+
+/*--------------------------------------------------------------------------
+  re8_decode_rank_of_permutation(rank, xs, x)
+  DECODING OF THE RANK OF THE PERMUTATION OF xs
+  (i) rank: index (rank) of a permutation
+  (i) xs:   signed leader in RE8 (8-dimensional integer vector)
+  (o) x:    point in RE8 (8-dimensional integer vector)
+  --------------------------------------------------------------------------
+ */
+static void re8_decode_rank_of_permutation(int rank, int *xs, SHORT x[8]) {
+  INT a[8], w[8], B, fac, fac_B, target;
+  int i, j;
+
+  /* --- pre-processing based on the signed leader xs ---
+     - compute the alphabet a=[a[0] ... a[q-1]] of x (q elements)
+       such that a[0]!=...!=a[q-1]
+       it is assumed that xs is sorted in the form of a signed leader
+       which can be summarized in 2 requirements:
+          a) |xs[0]| >= |xs[1]| >= |xs[2]| >= ... >= |xs[7]|
+          b) if |xs[i]|=|xs[i-1]|, xs[i]>=xs[i+1]
+       where |.| indicates the absolute value operator
+     - compute q (the number of symbols in the alphabet)
+     - compute w[0..q-1] where w[j] counts the number of occurences of
+       the symbol a[j] in xs
+     - compute B = prod_j=0..q-1 (w[j]!) where .! is the factorial */
+  /* xs[i], xs[i-1] and ptr_w/a*/
+  j = 0;
+  w[j] = 1;
+  a[j] = xs[0];
+  B = 1;
+  for (i = 1; i < 8; i++) {
+    if (xs[i] != xs[i - 1]) {
+      j++;
+      w[j] = 1;
+      a[j] = xs[i];
+    } else {
+      w[j]++;
+      B *= w[j];
+    }
+  }
+
+  /* --- actual rank decoding ---
+     the rank of x (where x is a permutation of xs) is based on
+     Schalkwijk's formula
+     it is given by rank=sum_{k=0..7} (A_k * fac_k/B_k)
+     the decoding of this rank is sequential and reconstructs x[0..7]
+     element by element from x[0] to x[7]
+     [the tricky part is the inference of A_k for each k...]
+   */
+
+  if (w[0] == 8) {
+    for (i = 0; i < 8; i++) {
+      x[i] = a[0]; /* avoid fac of 40320 */
+    }
+  } else {
+    target = rank * B;
+    fac_B = 1;
+    /* decode x element by element */
+    for (i = 0; i < 8; i++) {
+      fac = fac_B * fdk_dec_tab_factorial[i]; /* fac = 1..5040 */
+      j = -1;
+      do {
+        target -= w[++j] * fac;
+      } while (target >= 0); /* max of 30 tests / SV */
+      x[i] = a[j];
+      /* update rank, denominator B (B_k) and counter w[j] */
+      target += w[j] * fac; /* target = fac_B*B*rank */
+      fac_B *= w[j];
+      w[j]--;
+    }
+  }
+}
+
+/*--------------------------------------------------------------------------
+  re8_decode_base_index(n, I, y)
+  DECODING OF AN INDEX IN Qn (n=0,2,3 or 4)
+  (i) n: codebook number (*n is an integer defined in {0,2,3,4})
+  (i) I: index of c (pointer to unsigned 16-bit word)
+  (o) y: point in RE8 (8-dimensional integer vector)
+  note: the index I is defined as a 32-bit word, but only
+  16 bits are required (long can be replaced by unsigned integer)
+  --------------------------------------------------------------------------
+ */
+static void re8_decode_base_index(int *n, UINT index, SHORT y[8]) {
+  int i, im, t, sign_code, ka, ks, rank, leader[8];
+
+  if (*n < 2) {
+    for (i = 0; i < 8; i++) {
+      y[i] = 0;
+    }
+  } else {
+    // index = (unsigned int)*I;
+    /* search for the identifier ka of the absolute leader (table-lookup)
+       Q2 is a subset of Q3 - the two cases are considered in the same branch
+     */
+    switch (*n) {
+      case 2:
+      case 3:
+        i = table_lookup(fdk_dec_I3, index, NB_LDQ3);
+        ka = fdk_dec_A3[i];
+        break;
+      case 4:
+        i = table_lookup(fdk_dec_I4, index, NB_LDQ4);
+        ka = fdk_dec_A4[i];
+        break;
+      default:
+        FDK_ASSERT(0);
+        return;
+    }
+    /* reconstruct the absolute leader */
+    for (i = 0; i < 8; i++) {
+      leader[i] = fdk_dec_Da[ka][i];
+    }
+    /* search for the identifier ks of the signed leader (table look-up)
+       (this search is focused based on the identifier ka of the absolute
+        leader)*/
+    t = fdk_dec_Ia[ka];
+    im = fdk_dec_Ns[ka];
+    ks = table_lookup(fdk_dec_Is + t, index, im);
+
+    /* reconstruct the signed leader from its sign code */
+    sign_code = 2 * fdk_dec_Ds[t + ks];
+    for (i = 7; i >= 0; i--) {
+      leader[i] *= (1 - (sign_code & 2));
+      sign_code >>= 1;
+    }
+
+    /* compute and decode the rank of the permutation */
+    rank = index - fdk_dec_Is[t + ks]; /* rank = index - cardinality offset */
+
+    re8_decode_rank_of_permutation(rank, leader, y);
+  }
+  return;
+}
+
+/* re8_y2k(y,m,k)
+   VORONOI INDEXING (INDEX DECODING) k -> y
+   (i) k: Voronoi index k[0..7]
+   (i) m: Voronoi modulo (m = 2^r = 1<<r, where r is integer >=2)
+   (i) r: Voronoi order  (m = 2^r = 1<<r, where r is integer >=2)
+   (o) y: 8-dimensional point y[0..7] in RE8
+ */
+static void re8_k2y(int *k, int r, SHORT *y) {
+  int i, tmp, sum;
+  SHORT v[8];
+  FIXP_ZF zf[8];
+
+  FDK_ASSERT(r <= ZF_SCALE);
+
+  /* compute y = k M and z=(y-a)/m, where
+     M = [4        ]
+         [2 2      ]
+         [|   \    ]
+         [2     2  ]
+         [1 1 _ 1 1]
+     a=(2,0,...,0)
+     m = 1<<r
+  */
+  for (i = 0; i < 8; i++) {
+    y[i] = k[7];
+  }
+  zf[7] = INT2ZF(y[7], r);
+  sum = 0;
+  for (i = 6; i >= 1; i--) {
+    tmp = 2 * k[i];
+    sum += tmp;
+    y[i] += tmp;
+    zf[i] = INT2ZF(y[i], r);
+  }
+  y[0] += (4 * k[0] + sum);
+  zf[0] = INT2ZF(y[0] - 2, r);
+  /* find nearest neighbor v of z in infinite RE8 */
+  RE8_PPV(zf, v, r);
+  /* compute y -= m v */
+  for (i = 0; i < 8; i++) {
+    y[i] -= (SHORT)(v[i] << r);
+  }
+}
+
+/*--------------------------------------------------------------------------
+  RE8_dec(n, I, k, y)
+  MULTI-RATE INDEXING OF A POINT y in THE LATTICE RE8 (INDEX DECODING)
+  (i) n: codebook number (*n is an integer defined in {0,2,3,4,..,n_max}). n_max
+  = 36 (i) I: index of c (pointer to unsigned 16-bit word) (i) k: index of v
+  (8-dimensional vector of binary indices) = Voronoi index (o) y: point in RE8
+  (8-dimensional integer vector) note: the index I is defined as a 32-bit word,
+  but only 16 bits are required (long can be replaced by unsigned integer)
+
+  return 0 on success, -1 on error.
+  --------------------------------------------------------------------------
+ */
+static int RE8_dec(int n, int I, int *k, FIXP_DBL *y) {
+  SHORT v[8];
+  SHORT _y[8];
+  UINT r;
+  int i;
+
+  /* Check bound of codebook qn */
+  if (n > NQ_MAX) {
+    return -1;
+  }
+
+  /* decode the sub-indices I and kv[] according to the codebook number n:
+     if n=0,2,3,4, decode I (no Voronoi extension)
+     if n>4, Voronoi extension is used, decode I and kv[] */
+  if (n <= 4) {
+    re8_decode_base_index(&n, I, _y);
+    for (i = 0; i < 8; i++) {
+      y[i] = (LONG)_y[i];
+    }
+  } else {
+    /* compute the Voronoi modulo m = 2^r where r is extension order */
+    r = ((n - 3) >> 1);
+
+    while (n > 4) {
+      n -= 2;
+    }
+    /* decode base codebook index I into c (c is an element of Q3 or Q4)
+       [here c is stored in y to save memory] */
+    re8_decode_base_index(&n, I, _y);
+    /* decode Voronoi index k[] into v */
+    re8_k2y(k, r, v);
+    /* reconstruct y as y = m c + v (with m=2^r, r integer >=1) */
+    for (i = 0; i < 8; i++) {
+      y[i] = (LONG)((_y[i] << r) + v[i]);
+    }
+  }
+  return 0;
+}
+
+/**************************/
+/* start LPC decode stuff */
+/**************************/
+//#define M         16
+#define FREQ_MAX 6400.0f
+#define FREQ_DIV 400.0f
+#define LSF_GAP 50.0f
+
+/**
+ * \brief calculate inverse weighting factor and add non-weighted residual
+ *        LSF vector to first stage LSF approximation
+ * \param lsfq first stage LSF approximation values.
+ * \param xq weighted residual LSF vector
+ * \param nk_mode code book number coding mode.
+ */
+static void lsf_weight_2st(FIXP_LPC *lsfq, FIXP_DBL *xq, int nk_mode) {
+  FIXP_LPC d[M_LP_FILTER_ORDER + 1];
+  FIXP_SGL factor;
+  LONG w; /* inverse weight factor */
+  int i;
+
+  /* compute lsf distance */
+  d[0] = lsfq[0];
+  d[M_LP_FILTER_ORDER] =
+      FL2FXCONST_LPC(FREQ_MAX / (1 << LSF_SCALE)) - lsfq[M_LP_FILTER_ORDER - 1];
+  for (i = 1; i < M_LP_FILTER_ORDER; i++) {
+    d[i] = lsfq[i] - lsfq[i - 1];
+  }
+
+  switch (nk_mode) {
+    case 0:
+      factor = FL2FXCONST_SGL(2.0f * 60.0f / FREQ_DIV);
+      break; /* abs */
+    case 1:
+      factor = FL2FXCONST_SGL(2.0f * 65.0f / FREQ_DIV);
+      break; /* mid */
+    case 2:
+      factor = FL2FXCONST_SGL(2.0f * 64.0f / FREQ_DIV);
+      break; /* rel1 */
+    default:
+      factor = FL2FXCONST_SGL(2.0f * 63.0f / FREQ_DIV);
+      break; /* rel2 */
+  }
+  /* add non-weighted residual LSF vector to LSF1st */
+  for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+    w = (LONG)fMultDiv2(factor, sqrtFixp(fMult(d[i], d[i + 1])));
+    lsfq[i] = fAddSaturate(lsfq[i], FX_DBL2FX_LPC((FIXP_DBL)(w * (LONG)xq[i])));
+  }
+
+  return;
+}
+
+/**
+ * \brief decode nqn amount of code book numbers. These values determine the
+ * amount of following bits for nqn AVQ RE8 vectors.
+ * \param nk_mode quantization mode.
+ * \param nqn amount code book number to read.
+ * \param qn pointer to output buffer to hold decoded code book numbers qn.
+ */
+static void decode_qn(HANDLE_FDK_BITSTREAM hBs, int nk_mode, int nqn,
+                      int qn[]) {
+  int n;
+
+  if (nk_mode == 1) { /* nk mode 1 */
+    /* Unary code for mid LPC1/LPC3 */
+    /* Q0=0, Q2=10, Q3=110, ... */
+    for (n = 0; n < nqn; n++) {
+      qn[n] = get_vlclbf(hBs);
+      if (qn[n] > 0) {
+        qn[n]++;
+      }
+    }
+  } else { /* nk_mode 0, 3 and 2 */
+    /* 2 bits to specify Q2,Q3,Q4,ext */
+    for (n = 0; n < nqn; n++) {
+      qn[n] = 2 + FDKreadBits(hBs, 2);
+    }
+    if (nk_mode == 2) {
+      /* Unary code for rel LPC1/LPC3 */
+      /* Q0 = 0, Q5=10, Q6=110, ... */
+      for (n = 0; n < nqn; n++) {
+        if (qn[n] > 4) {
+          qn[n] = get_vlclbf(hBs);
+          if (qn[n] > 0) qn[n] += 4;
+        }
+      }
+    } else { /* nk_mode == (0 and 3) */
+      /* Unary code for abs and rel LPC0/LPC2 */
+      /* Q5 = 0, Q6=10, Q0=110, Q7=1110, ... */
+      for (n = 0; n < nqn; n++) {
+        if (qn[n] > 4) {
+          qn[n] = get_vlclbf(hBs);
+          switch (qn[n]) {
+            case 0:
+              qn[n] = 5;
+              break;
+            case 1:
+              qn[n] = 6;
+              break;
+            case 2:
+              qn[n] = 0;
+              break;
+            default:
+              qn[n] += 4;
+              break;
+          }
+        }
+      }
+    }
+  }
+}
+
+/**
+ * \brief reorder LSF coefficients to minimum distance.
+ * \param lsf pointer to buffer containing LSF coefficients and where reordered
+ * LSF coefficients will be stored into, scaled by LSF_SCALE.
+ * \param min_dist min distance scaled by LSF_SCALE
+ * \param n number of LSF/LSP coefficients.
+ */
+static void reorder_lsf(FIXP_LPC *lsf, FIXP_LPC min_dist, int n) {
+  FIXP_LPC lsf_min;
+  int i;
+
+  lsf_min = min_dist;
+  for (i = 0; i < n; i++) {
+    if (lsf[i] < lsf_min) {
+      lsf[i] = lsf_min;
+    }
+    lsf_min = fAddSaturate(lsf[i], min_dist);
+  }
+
+  /* reverse */
+  lsf_min = FL2FXCONST_LPC(FREQ_MAX / (1 << LSF_SCALE)) - min_dist;
+  for (i = n - 1; i >= 0; i--) {
+    if (lsf[i] > lsf_min) {
+      lsf[i] = lsf_min;
+    }
+
+    lsf_min = lsf[i] - min_dist;
+  }
+}
+
+/**
+ * \brief First stage approximation
+ * \param hBs bitstream handle as data source
+ * \param lsfq pointer to output buffer to hold LPC coefficients scaled by
+ * LSF_SCALE.
+ */
+static void vlpc_1st_dec(
+    HANDLE_FDK_BITSTREAM hBs, /* input:  codebook index                  */
+    FIXP_LPC *lsfq            /* i/o:    i:prediction   o:quantized lsf  */
+) {
+  const FIXP_LPC *p_dico;
+  int i, index;
+
+  index = FDKreadBits(hBs, 8);
+  p_dico = &fdk_dec_dico_lsf_abs_8b[index * M_LP_FILTER_ORDER];
+  for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+    lsfq[i] = p_dico[i];
+  }
+}
+
+/**
+ * \brief Do first stage approximation weighting and multiply with AVQ
+ * refinement.
+ * \param hBs bitstream handle data ssource.
+ * \param lsfq buffer holding 1st stage approx, 2nd stage approx is added to
+ * this values.
+ * \param nk_mode quantization mode.
+ * \return 0 on success, -1 on error.
+ */
+static int vlpc_2st_dec(
+    HANDLE_FDK_BITSTREAM hBs,
+    FIXP_LPC *lsfq, /* i/o:    i:1st stage   o:1st+2nd stage   */
+    int nk_mode     /* input:  0=abs, >0=rel                   */
+) {
+  int err;
+  FIXP_DBL xq[M_LP_FILTER_ORDER]; /* weighted residual LSF vector */
+
+  /* Decode AVQ refinement */
+  { err = CLpc_DecodeAVQ(hBs, xq, nk_mode, 2, 8); }
+  if (err != 0) {
+    return -1;
+  }
+
+  /* add non-weighted residual LSF vector to LSF1st */
+  lsf_weight_2st(lsfq, xq, nk_mode);
+
+  /* reorder */
+  reorder_lsf(lsfq, FL2FXCONST_LPC(LSF_GAP / (1 << LSF_SCALE)),
+              M_LP_FILTER_ORDER);
+
+  return 0;
+}
+
+/*
+ * Externally visible functions
+ */
+
+int CLpc_DecodeAVQ(HANDLE_FDK_BITSTREAM hBs, FIXP_DBL *pOutput, int nk_mode,
+                   int no_qn, int length) {
+  int i, l;
+
+  for (i = 0; i < length; i += 8 * no_qn) {
+    int qn[2], nk, n, I;
+    int kv[8] = {0};
+
+    decode_qn(hBs, nk_mode, no_qn, qn);
+
+    for (l = 0; l < no_qn; l++) {
+      if (qn[l] == 0) {
+        FDKmemclear(&pOutput[i + l * 8], 8 * sizeof(FIXP_DBL));
+      }
+
+      /* Voronoi extension order ( nk ) */
+      nk = 0;
+      n = qn[l];
+      if (qn[l] > 4) {
+        nk = (qn[l] - 3) >> 1;
+        n = qn[l] - nk * 2;
+      }
+
+      /* Base codebook index, in reverse bit group order (!) */
+      I = FDKreadBits(hBs, 4 * n);
+
+      if (nk > 0) {
+        int j;
+
+        for (j = 0; j < 8; j++) {
+          kv[j] = FDKreadBits(hBs, nk);
+        }
+      }
+
+      if (RE8_dec(qn[l], I, kv, &pOutput[i + l * 8]) != 0) {
+        return -1;
+      }
+    }
+  }
+  return 0;
+}
+
+int 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) {
+  int i, k, err;
+  int mode_lpc_bin = 0; /* mode_lpc bitstream representation */
+  int lpc_present[5] = {0, 0, 0, 0, 0};
+  int lpc0_available = 1;
+  int s = 0;
+  int l = 3;
+  const int nbDiv = NB_DIV;
+
+  lpc_present[4 >> s] = 1; /* LPC4 */
+
+  /* Decode LPC filters in the following order: LPC 4,0,2,1,3 */
+
+  /*** Decode LPC4 ***/
+  vlpc_1st_dec(hBs, lsp[4 >> s]);
+  err = vlpc_2st_dec(hBs, lsp[4 >> s], 0); /* nk_mode = 0 */
+  if (err != 0) {
+    return err;
+  }
+
+  /*** Decode LPC0 and LPC2 ***/
+  k = 0;
+  if (!first_lpd_flag) {
+    lpc_present[0] = 1;
+    lpc0_available = !last_lpc_lost;
+    /* old LPC4 is new LPC0 */
+    for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+      lsp[0][i] = lpc4_lsf[i];
+    }
+    /* skip LPC0 and continue with LPC2 */
+    k = 2;
+  }
+
+  for (; k < l; k += 2) {
+    int nk_mode = 0;
+
+    if ((k == 2) && (mod[0] == 3)) {
+      break; /* skip LPC2 */
+    }
+
+    lpc_present[k >> s] = 1;
+
+    mode_lpc_bin = FDKreadBit(hBs);
+
+    if (mode_lpc_bin == 0) {
+      /* LPC0/LPC2: Abs */
+      vlpc_1st_dec(hBs, lsp[k >> s]);
+    } else {
+      /* LPC0/LPC2: RelR */
+      for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+        lsp[k >> s][i] = lsp[4 >> s][i];
+      }
+      nk_mode = 3;
+    }
+
+    err = vlpc_2st_dec(hBs, lsp[k >> s], nk_mode);
+    if (err != 0) {
+      return err;
+    }
+  }
+
+  /*** Decode LPC1 ***/
+  if (mod[0] < 2) { /* else: skip LPC1 */
+    lpc_present[1] = 1;
+    mode_lpc_bin = get_vlclbf_n(hBs, 2);
+
+    switch (mode_lpc_bin) {
+      case 1:
+        /* LPC1: abs */
+        vlpc_1st_dec(hBs, lsp[1]);
+        err = vlpc_2st_dec(hBs, lsp[1], 0);
+        if (err != 0) {
+          return err;
+        }
+        break;
+      case 2:
+        /* LPC1: mid0 (no second stage AVQ quantizer in this case) */
+        if (lpc0_available) { /* LPC0/lsf[0] might be zero some times */
+          for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+            lsp[1][i] = (lsp[0][i] >> 1) + (lsp[2][i] >> 1);
+          }
+        } else {
+          for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+            lsp[1][i] = lsp[2][i];
+          }
+        }
+        break;
+      case 0:
+        /* LPC1: RelR */
+        for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+          lsp[1][i] = lsp[2][i];
+        }
+        err = vlpc_2st_dec(hBs, lsp[1], 2 << s);
+        if (err != 0) {
+          return err;
+        }
+        break;
+    }
+  }
+
+  /*** Decode LPC3 ***/
+  if ((mod[2] < 2)) { /* else: skip LPC3 */
+    int nk_mode = 0;
+    lpc_present[3] = 1;
+
+    mode_lpc_bin = get_vlclbf_n(hBs, 3);
+
+    switch (mode_lpc_bin) {
+      case 1:
+        /* LPC3: abs */
+        vlpc_1st_dec(hBs, lsp[3]);
+        break;
+      case 0:
+        /* LPC3: mid */
+        for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+          lsp[3][i] = (lsp[2][i] >> 1) + (lsp[4][i] >> 1);
+        }
+        nk_mode = 1;
+        break;
+      case 2:
+        /* LPC3: relL */
+        for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+          lsp[3][i] = lsp[2][i];
+        }
+        nk_mode = 2;
+        break;
+      case 3:
+        /* LPC3: relR */
+        for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+          lsp[3][i] = lsp[4][i];
+        }
+        nk_mode = 2;
+        break;
+    }
+    err = vlpc_2st_dec(hBs, lsp[3], nk_mode);
+    if (err != 0) {
+      return err;
+    }
+  }
+
+  if (!lpc0_available && !last_frame_ok) {
+    /* LPC(0) was lost. Use next available LPC(k) instead */
+    for (k = 1; k < (nbDiv + 1); k++) {
+      if (lpc_present[k]) {
+        for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+#define LSF_INIT_TILT (0.25f)
+          if (mod[0] > 0) {
+            lsp[0][i] = FX_DBL2FX_LPC(
+                fMult(lsp[k][i], FL2FXCONST_SGL(1.0f - LSF_INIT_TILT)) +
+                fMult(fdk_dec_lsf_init[i], FL2FXCONST_SGL(LSF_INIT_TILT)));
+          } else {
+            lsp[0][i] = lsp[k][i];
+          }
+        }
+        break;
+      }
+    }
+  }
+
+  for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+    lpc4_lsf[i] = lsp[4 >> s][i];
+  }
+
+  {
+    FIXP_DBL divFac;
+    int last, numLpc = 0;
+
+    i = nbDiv;
+    do {
+      numLpc += lpc_present[i--];
+    } while (i >= 0 && numLpc < 3);
+
+    last = i;
+
+    switch (numLpc) {
+      case 3:
+        divFac = FL2FXCONST_DBL(1.0f / 3.0f);
+        break;
+      case 2:
+        divFac = FL2FXCONST_DBL(1.0f / 2.0f);
+        break;
+      default:
+        divFac = FL2FXCONST_DBL(1.0f);
+        break;
+    }
+
+    /* get the adaptive mean for the next (bad) frame */
+    for (k = 0; k < M_LP_FILTER_ORDER; k++) {
+      FIXP_DBL tmp = (FIXP_DBL)0;
+      for (i = nbDiv; i > last; i--) {
+        if (lpc_present[i]) {
+          tmp = fMultAdd(tmp >> 1, lsp[i][k], divFac);
+        }
+      }
+      lsf_adaptive_mean_cand[k] = FX_DBL2FX_LPC(tmp);
+    }
+  }
+
+  /* calculate stability factor Theta. Needed for ACELP decoder and concealment
+   */
+  {
+    FIXP_LPC *lsf_prev, *lsf_curr;
+    k = 0;
+
+    FDK_ASSERT(lpc_present[0] == 1 && lpc_present[4 >> s] == 1);
+    lsf_prev = lsp[0];
+    for (i = 1; i < (nbDiv + 1); i++) {
+      if (lpc_present[i]) {
+        FIXP_DBL tmp = (FIXP_DBL)0;
+        int j;
+        lsf_curr = lsp[i];
+
+        /* sum = tmp * 2^(LSF_SCALE*2 + 4) */
+        for (j = 0; j < M_LP_FILTER_ORDER; j++) {
+          tmp += fPow2Div2((FIXP_SGL)(lsf_curr[j] - lsf_prev[j])) >> 3;
+        }
+
+        /* tmp = (float)(FL2FXCONST_DBL(1.25f) - fMult(tmp,
+         * FL2FXCONST_DBL(1/400000.0f))); */
+        tmp = FL2FXCONST_DBL(1.25f / (1 << LSF_SCALE)) -
+              fMult(tmp, FL2FXCONST_DBL((1 << (LSF_SCALE + 4)) / 400000.0f));
+        if (tmp >= FL2FXCONST_DBL(1.0f / (1 << LSF_SCALE))) {
+          pStability[k] = FL2FXCONST_SGL(1.0f / 2.0f);
+        } else if (tmp < FL2FXCONST_DBL(0.0f)) {
+          pStability[k] = FL2FXCONST_SGL(0.0f);
+        } else {
+          pStability[k] = FX_DBL2FX_SGL(tmp << (LSF_SCALE - 1));
+        }
+
+        lsf_prev = lsf_curr;
+        k = i;
+      } else {
+        /* Mark stability value as undefined. */
+        pStability[i] = (FIXP_SGL)-1;
+      }
+    }
+  }
+
+  /* convert into LSP domain */
+  for (i = 0; i < (nbDiv + 1); i++) {
+    if (lpc_present[i]) {
+      for (k = 0; k < M_LP_FILTER_ORDER; k++) {
+        lsp[i][k] = FX_DBL2FX_LPC(
+            fixp_cos(fMult(lsp[i][k],
+                           FL2FXCONST_SGL((1 << LSPARG_SCALE) * M_PI / 6400.0)),
+                     LSF_SCALE - LSPARG_SCALE));
+      }
+    }
+  }
+
+  return 0;
+}
+
+void 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) {
+  int i, j;
+
+#define BETA (FL2FXCONST_SGL(0.25f))
+#define ONE_BETA (FL2FXCONST_SGL(0.75f))
+#define BFI_FAC (FL2FXCONST_SGL(0.90f))
+#define ONE_BFI_FAC (FL2FXCONST_SGL(0.10f))
+
+  /* Frame loss concealment (could be improved) */
+
+  if (first_lpd_flag) {
+    /* Reset past LSF values */
+    for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+      lsp[0][i] = lpc4_lsf[i] = fdk_dec_lsf_init[i];
+    }
+  } else {
+    /* old LPC4 is new LPC0 */
+    for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+      lsp[0][i] = lpc4_lsf[i];
+    }
+  }
+
+  /* LPC1 */
+  for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+    FIXP_LPC lsf_mean = FX_DBL2FX_LPC(fMult(BETA, fdk_dec_lsf_init[i]) +
+                                      fMult(ONE_BETA, lsf_adaptive_mean[i]));
+
+    lsp[1][i] = FX_DBL2FX_LPC(fMult(BFI_FAC, lpc4_lsf[i]) +
+                              fMult(ONE_BFI_FAC, lsf_mean));
+  }
+
+  /* LPC2 - LPC4 */
+  for (j = 2; j <= 4; j++) {
+    for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+      /* lsf_mean[i] =  FX_DBL2FX_LPC(fMult((FIXP_LPC)(BETA + j *
+         FL2FXCONST_LPC(0.1f)), fdk_dec_lsf_init[i])
+                                    + fMult((FIXP_LPC)(ONE_BETA - j *
+         FL2FXCONST_LPC(0.1f)), lsf_adaptive_mean[i])); */
+
+      FIXP_LPC lsf_mean = FX_DBL2FX_LPC(
+          fMult((FIXP_SGL)(BETA + (FIXP_SGL)(j * (INT)FL2FXCONST_SGL(0.1f))),
+                (FIXP_SGL)fdk_dec_lsf_init[i]) +
+          fMult(
+              (FIXP_SGL)(ONE_BETA - (FIXP_SGL)(j * (INT)FL2FXCONST_SGL(0.1f))),
+              lsf_adaptive_mean[i]));
+
+      lsp[j][i] = FX_DBL2FX_LPC(fMult(BFI_FAC, lsp[j - 1][i]) +
+                                fMult(ONE_BFI_FAC, lsf_mean));
+    }
+  }
+
+  /* Update past values for the future */
+  for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+    lpc4_lsf[i] = lsp[4][i];
+  }
+
+  /* convert into LSP domain */
+  for (j = 0; j < 5; j++) {
+    for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+      lsp[j][i] = FX_DBL2FX_LPC(fixp_cos(
+          fMult(lsp[j][i], FL2FXCONST_SGL((1 << LSPARG_SCALE) * M_PI / 6400.0)),
+          LSF_SCALE - LSPARG_SCALE));
+    }
+  }
+}
+
+void E_LPC_a_weight(FIXP_LPC *wA, const FIXP_LPC *A, int m) {
+  FIXP_DBL f;
+  int i;
+
+  f = FL2FXCONST_DBL(0.92f);
+  for (i = 0; i < m; i++) {
+    wA[i] = FX_DBL2FX_LPC(fMult(A[i], f));
+    f = fMult(f, FL2FXCONST_DBL(0.92f));
+  }
+}
+
+void CLpd_DecodeGain(FIXP_DBL *gain, INT *gain_e, int gain_code) {
+  /* gain * 2^(gain_e) = 10^(gain_code/28) */
+  *gain = fLdPow(
+      FL2FXCONST_DBL(3.3219280948873623478703194294894 / 4.0), /* log2(10)*/
+      2,
+      fMultDiv2((FIXP_DBL)gain_code << (DFRACT_BITS - 1 - 7),
+                FL2FXCONST_DBL(2.0f / 28.0f)),
+      7, gain_e);
+}
+
+  /**
+   * \brief *   Find the polynomial F1(z) or F2(z) from the LSPs.
+   * This is performed by expanding the product polynomials:
+   *
+   * F1(z) =   product   ( 1 - 2 LSP_i z^-1 + z^-2 )
+   *         i=0,2,4,6,8
+   * F2(z) =   product   ( 1 - 2 LSP_i z^-1 + z^-2 )
+   *         i=1,3,5,7,9
+   *
+   * where LSP_i are the LSPs in the cosine domain.
+   * R.A.Salami    October 1990
+   * \param lsp input, line spectral freq. (cosine domain)
+   * \param f output, the coefficients of F1 or F2, scaled by 8 bits
+   * \param n no of coefficients (m/2)
+   * \param flag 1 : F1(z) ; 2 : F2(z)
+   */
+
+#define SF_F 8
+
+static void get_lsppol(FIXP_LPC lsp[], FIXP_DBL f[], int n, int flag) {
+  FIXP_DBL b;
+  FIXP_LPC *plsp;
+  int i, j;
+
+  plsp = lsp + flag - 1;
+  f[0] = FL2FXCONST_DBL(1.0f / (1 << SF_F));
+  b = -FX_LPC2FX_DBL(*plsp);
+  f[1] = b >> (SF_F - 1);
+  for (i = 2; i <= n; i++) {
+    plsp += 2;
+    b = -FX_LPC2FX_DBL(*plsp);
+    f[i] = ((fMultDiv2(b, f[i - 1]) << 1) + (f[i - 2])) << 1;
+    for (j = i - 1; j > 1; j--) {
+      f[j] = f[j] + (fMultDiv2(b, f[j - 1]) << 2) + f[j - 2];
+    }
+    f[1] = f[1] + (b >> (SF_F - 1));
+  }
+  return;
+}
+
+#define NC M_LP_FILTER_ORDER / 2
+
+/**
+ * \brief lsp input LSP vector
+ * \brief a output LP filter coefficient vector scaled by SF_A_COEFFS.
+ */
+void E_LPC_f_lsp_a_conversion(FIXP_LPC *lsp, FIXP_LPC *a, INT *a_exp) {
+  FIXP_DBL f1[NC + 1], f2[NC + 1];
+  int i, k;
+
+  /*-----------------------------------------------------*
+   *  Find the polynomials F1(z) and F2(z)               *
+   *-----------------------------------------------------*/
+
+  get_lsppol(lsp, f1, NC, 1);
+  get_lsppol(lsp, f2, NC, 2);
+
+  /*-----------------------------------------------------*
+   *  Multiply F1(z) by (1+z^-1) and F2(z) by (1-z^-1)   *
+   *-----------------------------------------------------*/
+  for (i = NC; i > 0; i--) {
+    f1[i] += f1[i - 1];
+    f2[i] -= f2[i - 1];
+  }
+
+  FIXP_DBL aDBL[M_LP_FILTER_ORDER];
+
+  for (i = 1, k = M_LP_FILTER_ORDER - 1; i <= NC; i++, k--) {
+    FIXP_DBL tmp1, tmp2;
+
+    tmp1 = f1[i] >> 1;
+    tmp2 = f2[i] >> 1;
+
+    aDBL[i - 1] = (tmp1 + tmp2);
+    aDBL[k] = (tmp1 - tmp2);
+  }
+
+  int headroom_a = getScalefactor(aDBL, M_LP_FILTER_ORDER);
+
+  for (i = 0; i < M_LP_FILTER_ORDER; i++) {
+    a[i] = FX_DBL2FX_LPC(aDBL[i] << headroom_a);
+  }
+
+  *a_exp = 8 - headroom_a;
+}