1/*
2 *  Copyright (c) 2012 The WebRTC project authors. All Rights Reserved.
3 *
4 *  Use of this source code is governed by a BSD-style license
5 *  that can be found in the LICENSE file in the root of the source
6 *  tree. An additional intellectual property rights grant can be found
7 *  in the file PATENTS.  All contributing project authors may
8 *  be found in the AUTHORS file in the root of the source tree.
9 */
10
11/*
12 * entropy_coding.c
13 *
14 * This header file defines all of the functions used to arithmetically
15 * encode the iSAC bistream
16 *
17 */
18
19
20#include "entropy_coding.h"
21#include "settings.h"
22#include "arith_routines.h"
23#include "signal_processing_library.h"
24#include "spectrum_ar_model_tables.h"
25#include "lpc_tables.h"
26#include "pitch_gain_tables.h"
27#include "pitch_lag_tables.h"
28#include "encode_lpc_swb.h"
29#include "lpc_shape_swb12_tables.h"
30#include "lpc_shape_swb16_tables.h"
31#include "lpc_gain_swb_tables.h"
32#include "os_specific_inline.h"
33
34#include <math.h>
35#include <string.h>
36
37static const uint16_t kLpcVecPerSegmentUb12 = 5;
38static const uint16_t kLpcVecPerSegmentUb16 = 4;
39
40/* CDF array for encoder bandwidth (12 vs 16 kHz) indicator. */
41static const uint16_t kOneBitEqualProbCdf[3] = {
42    0, 32768, 65535 };
43
44/* Pointer to cdf array for encoder bandwidth (12 vs 16 kHz) indicator. */
45static const uint16_t* kOneBitEqualProbCdf_ptr[1] = {
46    kOneBitEqualProbCdf };
47
48/*
49 * Initial cdf index for decoder of encoded bandwidth
50 * (12 vs 16 kHz) indicator.
51 */
52static const uint16_t kOneBitEqualProbInitIndex[1] = { 1 };
53
54
55static const int kIsSWB12 = 1;
56
57/* compute correlation from power spectrum */
58static void FindCorrelation(int32_t* PSpecQ12, int32_t* CorrQ7) {
59  int32_t summ[FRAMESAMPLES / 8];
60  int32_t diff[FRAMESAMPLES / 8];
61  const int16_t* CS_ptrQ9;
62  int32_t sum;
63  int k, n;
64
65  for (k = 0; k < FRAMESAMPLES / 8; k++) {
66    summ[k] = (PSpecQ12[k] + PSpecQ12[FRAMESAMPLES_QUARTER - 1 - k] + 16) >> 5;
67    diff[k] = (PSpecQ12[k] - PSpecQ12[FRAMESAMPLES_QUARTER - 1 - k] + 16) >> 5;
68  }
69
70  sum = 2;
71  for (n = 0; n < FRAMESAMPLES / 8; n++) {
72    sum += summ[n];
73  }
74  CorrQ7[0] = sum;
75
76  for (k = 0; k < AR_ORDER; k += 2) {
77    sum = 0;
78    CS_ptrQ9 = WebRtcIsac_kCos[k];
79    for (n = 0; n < FRAMESAMPLES / 8; n++)
80      sum += (CS_ptrQ9[n] * diff[n] + 256) >> 9;
81    CorrQ7[k + 1] = sum;
82  }
83
84  for (k = 1; k < AR_ORDER; k += 2) {
85    sum = 0;
86    CS_ptrQ9 = WebRtcIsac_kCos[k];
87    for (n = 0; n < FRAMESAMPLES / 8; n++)
88      sum += (CS_ptrQ9[n] * summ[n] + 256) >> 9;
89    CorrQ7[k + 1] = sum;
90  }
91}
92
93/* compute inverse AR power spectrum */
94/* Changed to the function used in iSAC FIX for compatibility reasons */
95static void FindInvArSpec(const int16_t* ARCoefQ12,
96                          const int32_t gainQ10,
97                          int32_t* CurveQ16) {
98  int32_t CorrQ11[AR_ORDER + 1];
99  int32_t sum, tmpGain;
100  int32_t diffQ16[FRAMESAMPLES / 8];
101  const int16_t* CS_ptrQ9;
102  int k, n;
103  int16_t round, shftVal = 0, sh;
104
105  sum = 0;
106  for (n = 0; n < AR_ORDER + 1; n++) {
107    sum += WEBRTC_SPL_MUL(ARCoefQ12[n], ARCoefQ12[n]);   /* Q24 */
108  }
109  sum = ((sum >> 6) * 65 + 32768) >> 16;  /* Q8 */
110  CorrQ11[0] = (sum * gainQ10 + 256) >> 9;
111
112  /* To avoid overflow, we shift down gainQ10 if it is large.
113   * We will not lose any precision */
114  if (gainQ10 > 400000) {
115    tmpGain = gainQ10 >> 3;
116    round = 32;
117    shftVal = 6;
118  } else {
119    tmpGain = gainQ10;
120    round = 256;
121    shftVal = 9;
122  }
123
124  for (k = 1; k < AR_ORDER + 1; k++) {
125    sum = 16384;
126    for (n = k; n < AR_ORDER + 1; n++)
127      sum += WEBRTC_SPL_MUL(ARCoefQ12[n - k], ARCoefQ12[n]); /* Q24 */
128    sum >>= 15;
129    CorrQ11[k] = (sum * tmpGain + round) >> shftVal;
130  }
131  sum = CorrQ11[0] << 7;
132  for (n = 0; n < FRAMESAMPLES / 8; n++) {
133    CurveQ16[n] = sum;
134  }
135  for (k = 1; k < AR_ORDER; k += 2) {
136    for (n = 0; n < FRAMESAMPLES / 8; n++) {
137      CurveQ16[n] += (WebRtcIsac_kCos[k][n] * CorrQ11[k + 1] + 2) >> 2;
138    }
139  }
140
141  CS_ptrQ9 = WebRtcIsac_kCos[0];
142
143  /* If CorrQ11[1] too large we avoid getting overflow in the
144   * calculation by shifting */
145  sh = WebRtcSpl_NormW32(CorrQ11[1]);
146  if (CorrQ11[1] == 0) { /* Use next correlation */
147    sh = WebRtcSpl_NormW32(CorrQ11[2]);
148  }
149  if (sh < 9) {
150    shftVal = 9 - sh;
151  } else {
152    shftVal = 0;
153  }
154  for (n = 0; n < FRAMESAMPLES / 8; n++) {
155    diffQ16[n] = (CS_ptrQ9[n] * (CorrQ11[1] >> shftVal) + 2) >> 2;
156  }
157  for (k = 2; k < AR_ORDER; k += 2) {
158    CS_ptrQ9 = WebRtcIsac_kCos[k];
159    for (n = 0; n < FRAMESAMPLES / 8; n++) {
160      diffQ16[n] += (CS_ptrQ9[n] * (CorrQ11[k + 1] >> shftVal) + 2) >> 2;
161    }
162  }
163
164  for (k = 0; k < FRAMESAMPLES / 8; k++) {
165    CurveQ16[FRAMESAMPLES_QUARTER - 1 - k] = CurveQ16[k] -
166        (diffQ16[k] << shftVal);
167    CurveQ16[k] += diffQ16[k] << shftVal;
168  }
169}
170
171/* Generate array of dither samples in Q7. */
172static void GenerateDitherQ7Lb(int16_t* bufQ7, uint32_t seed,
173                               int length, int16_t AvgPitchGain_Q12) {
174  int   k, shft;
175  int16_t dither1_Q7, dither2_Q7, dither_gain_Q14;
176
177  /* This threshold should be equal to that in decode_spec(). */
178  if (AvgPitchGain_Q12 < 614) {
179    for (k = 0; k < length - 2; k += 3) {
180      /* New random unsigned int. */
181      seed = (seed * 196314165) + 907633515;
182
183      /* Fixed-point dither sample between -64 and 64 (Q7). */
184      /* dither = seed * 128 / 4294967295 */
185      dither1_Q7 = (int16_t)(((int)seed + 16777216) >> 25);
186
187      /* New random unsigned int. */
188      seed = (seed * 196314165) + 907633515;
189
190      /* Fixed-point dither sample between -64 and 64. */
191      dither2_Q7 = (int16_t)(((int)seed + 16777216) >> 25);
192
193      shft = (seed >> 25) & 15;
194      if (shft < 5) {
195        bufQ7[k]   = dither1_Q7;
196        bufQ7[k + 1] = dither2_Q7;
197        bufQ7[k + 2] = 0;
198      } else if (shft < 10) {
199        bufQ7[k]   = dither1_Q7;
200        bufQ7[k + 1] = 0;
201        bufQ7[k + 2] = dither2_Q7;
202      } else {
203        bufQ7[k]   = 0;
204        bufQ7[k + 1] = dither1_Q7;
205        bufQ7[k + 2] = dither2_Q7;
206      }
207    }
208  } else {
209    dither_gain_Q14 = (int16_t)(22528 - 10 * AvgPitchGain_Q12);
210
211    /* Dither on half of the coefficients. */
212    for (k = 0; k < length - 1; k += 2) {
213      /* New random unsigned int */
214      seed = (seed * 196314165) + 907633515;
215
216      /* Fixed-point dither sample between -64 and 64. */
217      dither1_Q7 = (int16_t)(((int)seed + 16777216) >> 25);
218
219      /* Dither sample is placed in either even or odd index. */
220      shft = (seed >> 25) & 1;     /* Either 0 or 1 */
221
222      bufQ7[k + shft] = (((dither_gain_Q14 * dither1_Q7) + 8192) >> 14);
223      bufQ7[k + 1 - shft] = 0;
224    }
225  }
226}
227
228
229
230/******************************************************************************
231 * GenerateDitherQ7LbUB()
232 *
233 * generate array of dither samples in Q7 There are less zeros in dither
234 * vector compared to GenerateDitherQ7Lb.
235 *
236 * A uniform random number generator with the range of [-64 64] is employed
237 * but the generated dithers are scaled by 0.35, a heuristic scaling.
238 *
239 * Input:
240 *      -seed               : the initial seed for the random number generator.
241 *      -length             : the number of dither values to be generated.
242 *
243 * Output:
244 *      -bufQ7              : pointer to a buffer where dithers are written to.
245 */
246static void GenerateDitherQ7LbUB(
247    int16_t* bufQ7,
248    uint32_t seed,
249    int length) {
250  int k;
251  for (k = 0; k < length; k++) {
252    /* new random unsigned int */
253    seed = (seed * 196314165) + 907633515;
254
255    /* Fixed-point dither sample between -64 and 64 (Q7). */
256    /* bufQ7 = seed * 128 / 4294967295 */
257    bufQ7[k] = (int16_t)(((int)seed + 16777216) >> 25);
258
259    /* Scale by 0.35. */
260    bufQ7[k] = (int16_t)WEBRTC_SPL_MUL_16_16_RSFT(bufQ7[k], 2048, 13);
261  }
262}
263
264/*
265 * Function to decode the complex spectrum from the bit stream
266 * returns the total number of bytes in the stream.
267 */
268int WebRtcIsac_DecodeSpec(Bitstr* streamdata, int16_t AvgPitchGain_Q12,
269                          enum ISACBand band, double* fr, double* fi) {
270  int16_t  DitherQ7[FRAMESAMPLES];
271  int16_t  data[FRAMESAMPLES];
272  int32_t  invARSpec2_Q16[FRAMESAMPLES_QUARTER];
273  uint16_t invARSpecQ8[FRAMESAMPLES_QUARTER];
274  int16_t  ARCoefQ12[AR_ORDER + 1];
275  int16_t  RCQ15[AR_ORDER];
276  int16_t  gainQ10;
277  int32_t  gain2_Q10, res;
278  int32_t  in_sqrt;
279  int32_t  newRes;
280  int k, len, i;
281  int is_12khz = !kIsSWB12;
282  int num_dft_coeff = FRAMESAMPLES;
283  /* Create dither signal. */
284  if (band == kIsacLowerBand) {
285    GenerateDitherQ7Lb(DitherQ7, streamdata->W_upper, FRAMESAMPLES,
286                       AvgPitchGain_Q12);
287  } else {
288    GenerateDitherQ7LbUB(DitherQ7, streamdata->W_upper, FRAMESAMPLES);
289    if (band == kIsacUpperBand12) {
290      is_12khz = kIsSWB12;
291      num_dft_coeff = FRAMESAMPLES_HALF;
292    }
293  }
294
295  /* Decode model parameters. */
296  if (WebRtcIsac_DecodeRc(streamdata, RCQ15) < 0)
297    return -ISAC_RANGE_ERROR_DECODE_SPECTRUM;
298
299  WebRtcSpl_ReflCoefToLpc(RCQ15, AR_ORDER, ARCoefQ12);
300
301  if (WebRtcIsac_DecodeGain2(streamdata, &gain2_Q10) < 0)
302    return -ISAC_RANGE_ERROR_DECODE_SPECTRUM;
303
304  /* Compute inverse AR power spectrum. */
305  FindInvArSpec(ARCoefQ12, gain2_Q10, invARSpec2_Q16);
306
307  /* Convert to magnitude spectrum,
308   * by doing square-roots (modified from SPLIB). */
309  res = 1 << (WebRtcSpl_GetSizeInBits(invARSpec2_Q16[0]) >> 1);
310  for (k = 0; k < FRAMESAMPLES_QUARTER; k++) {
311    in_sqrt = invARSpec2_Q16[k];
312    i = 10;
313
314    /* Negative values make no sense for a real sqrt-function. */
315    if (in_sqrt < 0)
316      in_sqrt = -in_sqrt;
317
318    newRes = (in_sqrt / res + res) >> 1;
319    do {
320      res = newRes;
321      newRes = (in_sqrt / res + res) >> 1;
322    } while (newRes != res && i-- > 0);
323
324    invARSpecQ8[k] = (int16_t)newRes;
325  }
326
327  len = WebRtcIsac_DecLogisticMulti2(data, streamdata, invARSpecQ8, DitherQ7,
328                                     num_dft_coeff, is_12khz);
329  /* Arithmetic decoding of spectrum. */
330  if (len < 1) {
331    return -ISAC_RANGE_ERROR_DECODE_SPECTRUM;
332  }
333
334  switch (band) {
335    case kIsacLowerBand: {
336      /* Scale down spectral samples with low SNR. */
337      int32_t p1;
338      int32_t p2;
339      if (AvgPitchGain_Q12 <= 614) {
340        p1 = 30 << 10;
341        p2 = 32768 + (33 << 16);
342      } else {
343        p1 = 36 << 10;
344        p2 = 32768 + (40 << 16);
345      }
346      for (k = 0; k < FRAMESAMPLES; k += 4) {
347        gainQ10 = WebRtcSpl_DivW32W16ResW16(p1, (int16_t)(
348            (invARSpec2_Q16[k >> 2] + p2) >> 16));
349        *fr++ = (double)((data[ k ] * gainQ10 + 512) >> 10) / 128.0;
350        *fi++ = (double)((data[k + 1] * gainQ10 + 512) >> 10) / 128.0;
351        *fr++ = (double)((data[k + 2] * gainQ10 + 512) >> 10) / 128.0;
352        *fi++ = (double)((data[k + 3] * gainQ10 + 512) >> 10) / 128.0;
353      }
354      break;
355    }
356    case kIsacUpperBand12: {
357      for (k = 0, i = 0; k < FRAMESAMPLES_HALF; k += 4) {
358        fr[i] = (double)data[ k ] / 128.0;
359        fi[i] = (double)data[k + 1] / 128.0;
360        i++;
361        fr[i] = (double)data[k + 2] / 128.0;
362        fi[i] = (double)data[k + 3] / 128.0;
363        i++;
364      }
365      /* The second half of real and imaginary coefficients is zero. This is
366       * due to using the old FFT module which requires two signals as input
367       * while in 0-12 kHz mode we only have 8-12 kHz band, and the second
368       * signal is set to zero. */
369      memset(&fr[FRAMESAMPLES_QUARTER], 0, FRAMESAMPLES_QUARTER *
370             sizeof(double));
371      memset(&fi[FRAMESAMPLES_QUARTER], 0, FRAMESAMPLES_QUARTER *
372             sizeof(double));
373      break;
374    }
375    case kIsacUpperBand16: {
376      for (i = 0, k = 0; k < FRAMESAMPLES; k += 4, i++) {
377        fr[i] = (double)data[ k ] / 128.0;
378        fi[i] = (double)data[k + 1] / 128.0;
379        fr[(FRAMESAMPLES_HALF) - 1 - i] = (double)data[k + 2] / 128.0;
380        fi[(FRAMESAMPLES_HALF) - 1 - i] = (double)data[k + 3] / 128.0;
381      }
382      break;
383    }
384  }
385  return len;
386}
387
388
389int WebRtcIsac_EncodeSpec(const int16_t* fr, const int16_t* fi,
390                          int16_t AvgPitchGain_Q12, enum ISACBand band,
391                          Bitstr* streamdata) {
392  int16_t ditherQ7[FRAMESAMPLES];
393  int16_t dataQ7[FRAMESAMPLES];
394  int32_t PSpec[FRAMESAMPLES_QUARTER];
395  int32_t invARSpec2_Q16[FRAMESAMPLES_QUARTER];
396  uint16_t invARSpecQ8[FRAMESAMPLES_QUARTER];
397  int32_t CorrQ7[AR_ORDER + 1];
398  int32_t CorrQ7_norm[AR_ORDER + 1];
399  int16_t RCQ15[AR_ORDER];
400  int16_t ARCoefQ12[AR_ORDER + 1];
401  int32_t gain2_Q10;
402  int16_t val;
403  int32_t nrg, res;
404  uint32_t sum;
405  int32_t in_sqrt;
406  int32_t newRes;
407  int16_t err;
408  uint32_t nrg_u32;
409  int shift_var;
410  int k, n, j, i;
411  int is_12khz = !kIsSWB12;
412  int num_dft_coeff = FRAMESAMPLES;
413
414  /* Create dither signal. */
415  if (band == kIsacLowerBand) {
416    GenerateDitherQ7Lb(ditherQ7, streamdata->W_upper, FRAMESAMPLES,
417                       AvgPitchGain_Q12);
418  } else {
419    GenerateDitherQ7LbUB(ditherQ7, streamdata->W_upper, FRAMESAMPLES);
420    if (band == kIsacUpperBand12) {
421      is_12khz = kIsSWB12;
422      num_dft_coeff = FRAMESAMPLES_HALF;
423    }
424  }
425
426  /* add dither and quantize, and compute power spectrum */
427  switch (band) {
428    case kIsacLowerBand: {
429      for (k = 0; k < FRAMESAMPLES; k += 4) {
430        val = ((*fr++ + ditherQ7[k]   + 64) & 0xFF80) - ditherQ7[k];
431        dataQ7[k] = val;
432        sum = val * val;
433
434        val = ((*fi++ + ditherQ7[k + 1] + 64) & 0xFF80) - ditherQ7[k + 1];
435        dataQ7[k + 1] = val;
436        sum += val * val;
437
438        val = ((*fr++ + ditherQ7[k + 2] + 64) & 0xFF80) - ditherQ7[k + 2];
439        dataQ7[k + 2] = val;
440        sum += val * val;
441
442        val = ((*fi++ + ditherQ7[k + 3] + 64) & 0xFF80) - ditherQ7[k + 3];
443        dataQ7[k + 3] = val;
444        sum += val * val;
445
446        PSpec[k >> 2] = sum >> 2;
447      }
448      break;
449    }
450    case kIsacUpperBand12: {
451      for (k = 0, j = 0; k < FRAMESAMPLES_HALF; k += 4) {
452        val = ((*fr++ + ditherQ7[k]   + 64) & 0xFF80) - ditherQ7[k];
453        dataQ7[k] = val;
454        sum = val * val;
455
456        val = ((*fi++ + ditherQ7[k + 1] + 64) & 0xFF80) - ditherQ7[k + 1];
457        dataQ7[k + 1] = val;
458        sum += val * val;
459
460        PSpec[j++] = sum >> 1;
461
462        val = ((*fr++ + ditherQ7[k + 2] + 64) & 0xFF80) - ditherQ7[k + 2];
463        dataQ7[k + 2] = val;
464        sum = val * val;
465
466        val = ((*fi++ + ditherQ7[k + 3] + 64) & 0xFF80) - ditherQ7[k + 3];
467        dataQ7[k + 3] = val;
468        sum += val * val;
469
470        PSpec[j++] = sum >> 1;
471      }
472      break;
473    }
474    case kIsacUpperBand16: {
475      for (j = 0, k = 0; k < FRAMESAMPLES; k += 4, j++) {
476        val = ((fr[j] + ditherQ7[k]   + 64) & 0xFF80) - ditherQ7[k];
477        dataQ7[k] = val;
478        sum = val * val;
479
480        val = ((fi[j] + ditherQ7[k + 1] + 64) & 0xFF80) - ditherQ7[k + 1];
481        dataQ7[k + 1] = val;
482        sum += val * val;
483
484        val = ((fr[(FRAMESAMPLES_HALF) - 1 - j] + ditherQ7[k + 2] + 64) &
485            0xFF80) - ditherQ7[k + 2];
486        dataQ7[k + 2] = val;
487        sum += val * val;
488
489        val = ((fi[(FRAMESAMPLES_HALF) - 1 - j] + ditherQ7[k + 3] + 64) &
490            0xFF80) - ditherQ7[k + 3];
491        dataQ7[k + 3] = val;
492        sum += val * val;
493
494        PSpec[k >> 2] = sum >> 2;
495      }
496      break;
497    }
498  }
499
500  /* compute correlation from power spectrum */
501  FindCorrelation(PSpec, CorrQ7);
502
503  /* Find AR coefficients */
504  /* Aumber of bit shifts to 14-bit normalize CorrQ7[0]
505   * (leaving room for sign) */
506  shift_var = WebRtcSpl_NormW32(CorrQ7[0]) - 18;
507
508  if (shift_var > 0) {
509    for (k = 0; k < AR_ORDER + 1; k++) {
510      CorrQ7_norm[k] = CorrQ7[k] << shift_var;
511    }
512  } else {
513    for (k = 0; k < AR_ORDER + 1; k++) {
514      CorrQ7_norm[k] = CorrQ7[k] >> (-shift_var);
515    }
516  }
517
518  /* Find RC coefficients. */
519  WebRtcSpl_AutoCorrToReflCoef(CorrQ7_norm, AR_ORDER, RCQ15);
520
521  /* Quantize & code RC Coefficient. */
522  WebRtcIsac_EncodeRc(RCQ15, streamdata);
523
524  /* RC -> AR coefficients */
525  WebRtcSpl_ReflCoefToLpc(RCQ15, AR_ORDER, ARCoefQ12);
526
527  /* Compute ARCoef' * Corr * ARCoef in Q19. */
528  nrg = 0;
529  for (j = 0; j <= AR_ORDER; j++) {
530    for (n = 0; n <= j; n++) {
531      nrg += (ARCoefQ12[j] * ((CorrQ7_norm[j - n] * ARCoefQ12[n] + 256) >> 9) +
532          4) >> 3;
533    }
534    for (n = j + 1; n <= AR_ORDER; n++) {
535      nrg += (ARCoefQ12[j] * ((CorrQ7_norm[n - j] * ARCoefQ12[n] + 256) >> 9) +
536          4) >> 3;
537    }
538  }
539
540  nrg_u32 = (uint32_t)nrg;
541  if (shift_var > 0) {
542    nrg_u32 = nrg_u32 >> shift_var;
543  } else {
544    nrg_u32 = nrg_u32 << (-shift_var);
545  }
546  if (nrg_u32 > 0x7FFFFFFF) {
547    nrg = 0x7FFFFFFF;
548  }  else {
549    nrg = (int32_t)nrg_u32;
550  }
551  /* Also shifts 31 bits to the left! */
552  gain2_Q10 = WebRtcSpl_DivResultInQ31(FRAMESAMPLES_QUARTER, nrg);
553
554  /* Quantize & code gain2_Q10. */
555  if (WebRtcIsac_EncodeGain2(&gain2_Q10, streamdata)) {
556    return -1;
557  }
558
559  /* Compute inverse AR power spectrum. */
560  FindInvArSpec(ARCoefQ12, gain2_Q10, invARSpec2_Q16);
561  /* Convert to magnitude spectrum, by doing square-roots
562   * (modified from SPLIB). */
563  res = 1 << (WebRtcSpl_GetSizeInBits(invARSpec2_Q16[0]) >> 1);
564  for (k = 0; k < FRAMESAMPLES_QUARTER; k++) {
565    in_sqrt = invARSpec2_Q16[k];
566    i = 10;
567    /* Negative values make no sense for a real sqrt-function. */
568    if (in_sqrt < 0) {
569      in_sqrt = -in_sqrt;
570    }
571    newRes = (in_sqrt / res + res) >> 1;
572    do {
573      res = newRes;
574      newRes = (in_sqrt / res + res) >> 1;
575    } while (newRes != res && i-- > 0);
576
577    invARSpecQ8[k] = (int16_t)newRes;
578  }
579  /* arithmetic coding of spectrum */
580  err = WebRtcIsac_EncLogisticMulti2(streamdata, dataQ7, invARSpecQ8,
581                                     num_dft_coeff, is_12khz);
582  if (err < 0) {
583    return (err);
584  }
585  return 0;
586}
587
588
589/* step-up */
590void WebRtcIsac_Rc2Poly(double* RC, int N, double* a) {
591  int m, k;
592  double tmp[MAX_AR_MODEL_ORDER];
593
594  a[0] = 1.0;
595  tmp[0] = 1.0;
596  for (m = 1; m <= N; m++) {
597    /* copy */
598    memcpy(&tmp[1], &a[1], (m - 1) * sizeof(double));
599    a[m] = RC[m - 1];
600    for (k = 1; k < m; k++) {
601      a[k] += RC[m - 1] * tmp[m - k];
602    }
603  }
604  return;
605}
606
607/* step-down */
608void WebRtcIsac_Poly2Rc(double* a, int N, double* RC) {
609  int m, k;
610  double tmp[MAX_AR_MODEL_ORDER];
611  double tmp_inv;
612
613  RC[N - 1] = a[N];
614  for (m = N - 1; m > 0; m--) {
615    tmp_inv = 1.0 / (1.0 - RC[m] * RC[m]);
616    for (k = 1; k <= m; k++) {
617      tmp[k] = (a[k] - RC[m] * a[m - k + 1]) * tmp_inv;
618    }
619
620    memcpy(&a[1], &tmp[1], (m - 1) * sizeof(double));
621    RC[m - 1] = tmp[m];
622  }
623  return;
624}
625
626
627#define MAX_ORDER 100
628
629/* Matlab's LAR definition */
630void WebRtcIsac_Rc2Lar(const double* refc, double* lar, int order) {
631  int k;
632  for (k = 0; k < order; k++) {
633    lar[k] = log((1 + refc[k]) / (1 - refc[k]));
634  }
635}
636
637
638void WebRtcIsac_Lar2Rc(const double* lar, double* refc,  int order) {
639  int k;
640  double tmp;
641
642  for (k = 0; k < order; k++) {
643    tmp = exp(lar[k]);
644    refc[k] = (tmp - 1) / (tmp + 1);
645  }
646}
647
648void WebRtcIsac_Poly2Lar(double* lowband, int orderLo, double* hiband,
649                         int orderHi, int Nsub, double* lars) {
650  int k;
651  double rc[MAX_ORDER], *inpl, *inph, *outp;
652
653  inpl = lowband;
654  inph = hiband;
655  outp = lars;
656  for (k = 0; k < Nsub; k++) {
657    /* gains */
658    outp[0] = inpl[0];
659    outp[1] = inph[0];
660    outp += 2;
661
662    /* Low band */
663    inpl[0] = 1.0;
664    WebRtcIsac_Poly2Rc(inpl, orderLo, rc);
665    WebRtcIsac_Rc2Lar(rc, outp, orderLo);
666    outp += orderLo;
667
668    /* High band */
669    inph[0] = 1.0;
670    WebRtcIsac_Poly2Rc(inph, orderHi, rc);
671    WebRtcIsac_Rc2Lar(rc, outp, orderHi);
672    outp += orderHi;
673
674    inpl += orderLo + 1;
675    inph += orderHi + 1;
676  }
677}
678
679
680int16_t WebRtcIsac_Poly2LarUB(double* lpcVecs, int16_t bandwidth) {
681  double      poly[MAX_ORDER];
682  double      rc[MAX_ORDER];
683  double*     ptrIO;
684  int16_t vecCntr;
685  int16_t vecSize;
686  int16_t numVec;
687
688  vecSize = UB_LPC_ORDER;
689  switch (bandwidth) {
690    case isac12kHz: {
691      numVec  = UB_LPC_VEC_PER_FRAME;
692      break;
693    }
694    case isac16kHz: {
695      numVec  = UB16_LPC_VEC_PER_FRAME;
696      break;
697    }
698    default:
699      return -1;
700  }
701
702  ptrIO = lpcVecs;
703  poly[0] = 1.0;
704  for (vecCntr = 0; vecCntr < numVec; vecCntr++) {
705    memcpy(&poly[1], ptrIO, sizeof(double) * vecSize);
706    WebRtcIsac_Poly2Rc(poly, vecSize, rc);
707    WebRtcIsac_Rc2Lar(rc, ptrIO, vecSize);
708    ptrIO += vecSize;
709  }
710  return 0;
711}
712
713
714void WebRtcIsac_Lar2Poly(double* lars, double* lowband, int orderLo,
715                         double* hiband, int orderHi, int Nsub) {
716  int k, orderTot;
717  double rc[MAX_ORDER], *outpl, *outph, *inp;
718
719  orderTot = (orderLo + orderHi + 2);
720  outpl = lowband;
721  outph = hiband;
722  /* First two elements of 'inp' store gains*/
723  inp = lars;
724  for (k = 0; k < Nsub; k++) {
725    /* Low band */
726    WebRtcIsac_Lar2Rc(&inp[2], rc, orderLo);
727    WebRtcIsac_Rc2Poly(rc, orderLo, outpl);
728
729    /* High band */
730    WebRtcIsac_Lar2Rc(&inp[orderLo + 2], rc, orderHi);
731    WebRtcIsac_Rc2Poly(rc, orderHi, outph);
732
733    /* gains */
734    outpl[0] = inp[0];
735    outph[0] = inp[1];
736
737    outpl += orderLo + 1;
738    outph += orderHi + 1;
739    inp += orderTot;
740  }
741}
742
743/*
744 *  assumes 2 LAR vectors interpolates to 'numPolyVec' A-polynomials
745 *  Note: 'numPolyVecs' includes the first and the last point of the interval
746 */
747void WebRtcIsac_Lar2PolyInterpolUB(double* larVecs, double* percepFilterParams,
748                                   int numPolyVecs) {
749  int polyCntr, coeffCntr;
750  double larInterpol[UB_LPC_ORDER];
751  double rc[UB_LPC_ORDER];
752  double delta[UB_LPC_ORDER];
753
754  /* calculate the step-size for linear interpolation coefficients */
755  for (coeffCntr = 0; coeffCntr < UB_LPC_ORDER; coeffCntr++) {
756    delta[coeffCntr] = (larVecs[UB_LPC_ORDER + coeffCntr] -
757        larVecs[coeffCntr]) / (numPolyVecs - 1);
758  }
759
760  for (polyCntr = 0; polyCntr < numPolyVecs; polyCntr++) {
761    for (coeffCntr = 0; coeffCntr < UB_LPC_ORDER; coeffCntr++) {
762      larInterpol[coeffCntr] = larVecs[coeffCntr] +
763          delta[coeffCntr] * polyCntr;
764    }
765    WebRtcIsac_Lar2Rc(larInterpol, rc, UB_LPC_ORDER);
766
767    /* convert to A-polynomial, the following function returns A[0] = 1;
768     * which is written where gains had to be written. Then we write the
769     * gain (outside this function). This way we say a memcpy. */
770    WebRtcIsac_Rc2Poly(rc, UB_LPC_ORDER, percepFilterParams);
771    percepFilterParams += (UB_LPC_ORDER + 1);
772  }
773}
774
775int WebRtcIsac_DecodeLpc(Bitstr* streamdata, double* LPCCoef_lo,
776                         double* LPCCoef_hi) {
777  double lars[KLT_ORDER_GAIN + KLT_ORDER_SHAPE];
778  int err;
779
780  err = WebRtcIsac_DecodeLpcCoef(streamdata, lars);
781  if (err < 0) {
782    return -ISAC_RANGE_ERROR_DECODE_LPC;
783  }
784  WebRtcIsac_Lar2Poly(lars, LPCCoef_lo, ORDERLO, LPCCoef_hi, ORDERHI,
785                      SUBFRAMES);
786  return 0;
787}
788
789int16_t WebRtcIsac_DecodeInterpolLpcUb(Bitstr* streamdata,
790                                       double* percepFilterParams,
791                                       int16_t bandwidth) {
792  double lpcCoeff[UB_LPC_ORDER * UB16_LPC_VEC_PER_FRAME];
793  int err;
794  int interpolCntr;
795  int subframeCntr;
796  int16_t numSegments;
797  int16_t numVecPerSegment;
798  int16_t numGains;
799
800  double percepFilterGains[SUBFRAMES << 1];
801  double* ptrOutParam = percepFilterParams;
802
803  err = WebRtcIsac_DecodeLpcCoefUB(streamdata, lpcCoeff, percepFilterGains,
804                                   bandwidth);
805  if (err < 0) {
806    return -ISAC_RANGE_ERROR_DECODE_LPC;
807  }
808
809  switch (bandwidth) {
810    case isac12kHz: {
811      numGains = SUBFRAMES;
812      numSegments = UB_LPC_VEC_PER_FRAME - 1;
813      numVecPerSegment = kLpcVecPerSegmentUb12;
814      break;
815    }
816    case isac16kHz: {
817      numGains = SUBFRAMES << 1;
818      numSegments = UB16_LPC_VEC_PER_FRAME - 1;
819      numVecPerSegment = kLpcVecPerSegmentUb16;
820      break;
821    }
822    default:
823      return -1;
824  }
825
826  for (interpolCntr = 0; interpolCntr < numSegments; interpolCntr++) {
827    WebRtcIsac_Lar2PolyInterpolUB(&lpcCoeff[interpolCntr * UB_LPC_ORDER],
828                                  ptrOutParam, numVecPerSegment + 1);
829    ptrOutParam += (numVecPerSegment * (UB_LPC_ORDER + 1));
830  }
831
832  ptrOutParam = percepFilterParams;
833
834  if (bandwidth == isac16kHz) {
835    ptrOutParam += (1 + UB_LPC_ORDER);
836  }
837
838  for (subframeCntr = 0; subframeCntr < numGains; subframeCntr++) {
839    *ptrOutParam = percepFilterGains[subframeCntr];
840    ptrOutParam += (1 + UB_LPC_ORDER);
841  }
842  return 0;
843}
844
845
846/* decode & dequantize LPC Coef */
847int WebRtcIsac_DecodeLpcCoef(Bitstr* streamdata, double* LPCCoef) {
848  int j, k, n, pos, pos2, posg, poss, offsg, offss, offs2;
849  int index_g[KLT_ORDER_GAIN], index_s[KLT_ORDER_SHAPE];
850  double tmpcoeffs_g[KLT_ORDER_GAIN], tmpcoeffs_s[KLT_ORDER_SHAPE];
851  double tmpcoeffs2_g[KLT_ORDER_GAIN], tmpcoeffs2_s[KLT_ORDER_SHAPE];
852  double sum;
853  int err;
854  int model = 1;
855
856  /* entropy decoding of model number */
857  /* We are keeping this for backward compatibility of bit-streams. */
858  err = WebRtcIsac_DecHistOneStepMulti(&model, streamdata,
859                                       WebRtcIsac_kQKltModelCdfPtr,
860                                       WebRtcIsac_kQKltModelInitIndex, 1);
861  if (err < 0) {
862    return err;
863  }
864  /* Only accepted value of model is 0. It is kept in bit-stream for backward
865   * compatibility. */
866  if (model != 0) {
867    return -ISAC_DISALLOWED_LPC_MODEL;
868  }
869
870  /* entropy decoding of quantization indices */
871  err = WebRtcIsac_DecHistOneStepMulti(
872      index_s, streamdata, WebRtcIsac_kQKltCdfPtrShape,
873      WebRtcIsac_kQKltInitIndexShape, KLT_ORDER_SHAPE);
874  if (err < 0) {
875    return err;
876  }
877  err = WebRtcIsac_DecHistOneStepMulti(
878      index_g, streamdata, WebRtcIsac_kQKltCdfPtrGain,
879      WebRtcIsac_kQKltInitIndexGain, KLT_ORDER_GAIN);
880  if (err < 0) {
881    return err;
882  }
883
884  /* find quantization levels for coefficients */
885  for (k = 0; k < KLT_ORDER_SHAPE; k++) {
886    tmpcoeffs_s[k] =
887        WebRtcIsac_kQKltLevelsShape[WebRtcIsac_kQKltOffsetShape[k] +
888                                    index_s[k]];
889  }
890  for (k = 0; k < KLT_ORDER_GAIN; k++) {
891    tmpcoeffs_g[k] = WebRtcIsac_kQKltLevelsGain[WebRtcIsac_kQKltOffsetGain[k] +
892                                                index_g[k]];
893  }
894
895  /* Inverse KLT  */
896
897  /* Left transform, transpose matrix!  */
898  offsg = 0;
899  offss = 0;
900  posg = 0;
901  poss = 0;
902  for (j = 0; j < SUBFRAMES; j++) {
903    offs2 = 0;
904    for (k = 0; k < LPC_GAIN_ORDER; k++) {
905      sum = 0;
906      pos = offsg;
907      pos2 = offs2;
908      for (n = 0; n < LPC_GAIN_ORDER; n++) {
909        sum += tmpcoeffs_g[pos++] * WebRtcIsac_kKltT1Gain[pos2++];
910      }
911      tmpcoeffs2_g[posg++] = sum;
912      offs2 += LPC_GAIN_ORDER;
913    }
914    offs2 = 0;
915    for (k = 0; k < LPC_SHAPE_ORDER; k++) {
916      sum = 0;
917      pos = offss;
918      pos2 = offs2;
919      for (n = 0; n < LPC_SHAPE_ORDER; n++) {
920        sum += tmpcoeffs_s[pos++] * WebRtcIsac_kKltT1Shape[pos2++];
921      }
922      tmpcoeffs2_s[poss++] = sum;
923      offs2 += LPC_SHAPE_ORDER;
924    }
925    offsg += LPC_GAIN_ORDER;
926    offss += LPC_SHAPE_ORDER;
927  }
928
929  /* Right transform, transpose matrix */
930  offsg = 0;
931  offss = 0;
932  posg = 0;
933  poss = 0;
934  for (j = 0; j < SUBFRAMES; j++) {
935    posg = offsg;
936    for (k = 0; k < LPC_GAIN_ORDER; k++) {
937      sum = 0;
938      pos = k;
939      pos2 = j;
940      for (n = 0; n < SUBFRAMES; n++) {
941        sum += tmpcoeffs2_g[pos] * WebRtcIsac_kKltT2Gain[pos2];
942        pos += LPC_GAIN_ORDER;
943        pos2 += SUBFRAMES;
944
945      }
946      tmpcoeffs_g[posg++] = sum;
947    }
948    poss = offss;
949    for (k = 0; k < LPC_SHAPE_ORDER; k++) {
950      sum = 0;
951      pos = k;
952      pos2 = j;
953      for (n = 0; n < SUBFRAMES; n++) {
954        sum += tmpcoeffs2_s[pos] * WebRtcIsac_kKltT2Shape[pos2];
955        pos += LPC_SHAPE_ORDER;
956        pos2 += SUBFRAMES;
957      }
958      tmpcoeffs_s[poss++] = sum;
959    }
960    offsg += LPC_GAIN_ORDER;
961    offss += LPC_SHAPE_ORDER;
962  }
963
964  /* scaling, mean addition, and gain restoration */
965  posg = 0;
966  poss = 0;
967  pos = 0;
968  for (k = 0; k < SUBFRAMES; k++) {
969    /* log gains */
970    LPCCoef[pos] = tmpcoeffs_g[posg] / LPC_GAIN_SCALE;
971    LPCCoef[pos] += WebRtcIsac_kLpcMeansGain[posg];
972    LPCCoef[pos] = exp(LPCCoef[pos]);
973    pos++;
974    posg++;
975    LPCCoef[pos] = tmpcoeffs_g[posg] / LPC_GAIN_SCALE;
976    LPCCoef[pos] += WebRtcIsac_kLpcMeansGain[posg];
977    LPCCoef[pos] = exp(LPCCoef[pos]);
978    pos++;
979    posg++;
980
981    /* Low-band LAR coefficients. */
982    for (n = 0; n < LPC_LOBAND_ORDER; n++, pos++, poss++) {
983      LPCCoef[pos] = tmpcoeffs_s[poss] / LPC_LOBAND_SCALE;
984      LPCCoef[pos] += WebRtcIsac_kLpcMeansShape[poss];
985    }
986
987    /* High-band LAR coefficients. */
988    for (n = 0; n < LPC_HIBAND_ORDER; n++, pos++, poss++) {
989      LPCCoef[pos] = tmpcoeffs_s[poss] / LPC_HIBAND_SCALE;
990      LPCCoef[pos] += WebRtcIsac_kLpcMeansShape[poss];
991    }
992  }
993  return 0;
994}
995
996/* Encode LPC in LAR domain. */
997void WebRtcIsac_EncodeLar(double* LPCCoef, Bitstr* streamdata,
998                          IsacSaveEncoderData* encData) {
999  int j, k, n, pos, pos2, poss, offss, offs2;
1000  int index_s[KLT_ORDER_SHAPE];
1001  int index_ovr_s[KLT_ORDER_SHAPE];
1002  double tmpcoeffs_s[KLT_ORDER_SHAPE];
1003  double tmpcoeffs2_s[KLT_ORDER_SHAPE];
1004  double sum;
1005  const int kModel = 0;
1006
1007  /* Mean removal and scaling. */
1008  poss = 0;
1009  pos = 0;
1010  for (k = 0; k < SUBFRAMES; k++) {
1011    /* First two element are gains, move over them. */
1012    pos += 2;
1013
1014    /* Low-band LAR coefficients. */
1015    for (n = 0; n < LPC_LOBAND_ORDER; n++, poss++, pos++) {
1016      tmpcoeffs_s[poss] = LPCCoef[pos] - WebRtcIsac_kLpcMeansShape[poss];
1017      tmpcoeffs_s[poss] *= LPC_LOBAND_SCALE;
1018    }
1019
1020    /* High-band LAR coefficients. */
1021    for (n = 0; n < LPC_HIBAND_ORDER; n++, poss++, pos++) {
1022      tmpcoeffs_s[poss] = LPCCoef[pos] - WebRtcIsac_kLpcMeansShape[poss];
1023      tmpcoeffs_s[poss] *= LPC_HIBAND_SCALE;
1024    }
1025  }
1026
1027  /* KLT  */
1028
1029  /* Left transform. */
1030  offss = 0;
1031  for (j = 0; j < SUBFRAMES; j++) {
1032    poss = offss;
1033    for (k = 0; k < LPC_SHAPE_ORDER; k++) {
1034      sum = 0;
1035      pos = offss;
1036      pos2 = k;
1037      for (n = 0; n < LPC_SHAPE_ORDER; n++) {
1038        sum += tmpcoeffs_s[pos++] * WebRtcIsac_kKltT1Shape[pos2];
1039        pos2 += LPC_SHAPE_ORDER;
1040      }
1041      tmpcoeffs2_s[poss++] = sum;
1042    }
1043    offss += LPC_SHAPE_ORDER;
1044  }
1045
1046  /* Right transform. */
1047  offss = 0;
1048  offs2 = 0;
1049  for (j = 0; j < SUBFRAMES; j++) {
1050    poss = offss;
1051    for (k = 0; k < LPC_SHAPE_ORDER; k++) {
1052      sum = 0;
1053      pos = k;
1054      pos2 = offs2;
1055      for (n = 0; n < SUBFRAMES; n++) {
1056        sum += tmpcoeffs2_s[pos] * WebRtcIsac_kKltT2Shape[pos2++];
1057        pos += LPC_SHAPE_ORDER;
1058      }
1059      tmpcoeffs_s[poss++] = sum;
1060    }
1061    offs2 += SUBFRAMES;
1062    offss += LPC_SHAPE_ORDER;
1063  }
1064
1065  /* Quantize coefficients. */
1066  for (k = 0; k < KLT_ORDER_SHAPE; k++) {
1067    index_s[k] = (WebRtcIsac_lrint(tmpcoeffs_s[k] / KLT_STEPSIZE)) +
1068        WebRtcIsac_kQKltQuantMinShape[k];
1069    if (index_s[k] < 0) {
1070      index_s[k] = 0;
1071    } else if (index_s[k] > WebRtcIsac_kQKltMaxIndShape[k]) {
1072      index_s[k] = WebRtcIsac_kQKltMaxIndShape[k];
1073    }
1074    index_ovr_s[k] = WebRtcIsac_kQKltOffsetShape[k] + index_s[k];
1075  }
1076
1077
1078  /* Only one model remains in this version of the code, kModel = 0. We
1079   * are keeping for bit-streams to be backward compatible. */
1080  /* entropy coding of model number */
1081  WebRtcIsac_EncHistMulti(streamdata, &kModel, WebRtcIsac_kQKltModelCdfPtr, 1);
1082
1083  /* Save data for creation of multiple bit streams */
1084  /* Entropy coding of quantization indices - shape only. */
1085  WebRtcIsac_EncHistMulti(streamdata, index_s, WebRtcIsac_kQKltCdfPtrShape,
1086                          KLT_ORDER_SHAPE);
1087
1088  /* Save data for creation of multiple bit streams. */
1089  for (k = 0; k < KLT_ORDER_SHAPE; k++) {
1090    encData->LPCindex_s[KLT_ORDER_SHAPE * encData->startIdx + k] = index_s[k];
1091  }
1092
1093  /* Find quantization levels for shape coefficients. */
1094  for (k = 0; k < KLT_ORDER_SHAPE; k++) {
1095    tmpcoeffs_s[k] = WebRtcIsac_kQKltLevelsShape[index_ovr_s[k]];
1096  }
1097  /* Inverse KLT.  */
1098  /* Left transform, transpose matrix.! */
1099  offss = 0;
1100  poss = 0;
1101  for (j = 0; j < SUBFRAMES; j++) {
1102    offs2 = 0;
1103    for (k = 0; k < LPC_SHAPE_ORDER; k++) {
1104      sum = 0;
1105      pos = offss;
1106      pos2 = offs2;
1107      for (n = 0; n < LPC_SHAPE_ORDER; n++) {
1108        sum += tmpcoeffs_s[pos++] * WebRtcIsac_kKltT1Shape[pos2++];
1109      }
1110      tmpcoeffs2_s[poss++] = sum;
1111      offs2 += LPC_SHAPE_ORDER;
1112    }
1113    offss += LPC_SHAPE_ORDER;
1114  }
1115
1116  /* Right transform, Transpose matrix */
1117  offss = 0;
1118  poss = 0;
1119  for (j = 0; j < SUBFRAMES; j++) {
1120    poss = offss;
1121    for (k = 0; k < LPC_SHAPE_ORDER; k++) {
1122      sum = 0;
1123      pos = k;
1124      pos2 = j;
1125      for (n = 0; n < SUBFRAMES; n++) {
1126        sum += tmpcoeffs2_s[pos] * WebRtcIsac_kKltT2Shape[pos2];
1127        pos += LPC_SHAPE_ORDER;
1128        pos2 += SUBFRAMES;
1129      }
1130      tmpcoeffs_s[poss++] = sum;
1131    }
1132    offss += LPC_SHAPE_ORDER;
1133  }
1134
1135  /* Scaling, mean addition, and gain restoration. */
1136  poss = 0;
1137  pos = 0;
1138  for (k = 0; k < SUBFRAMES; k++) {
1139    /* Ignore gains. */
1140    pos += 2;
1141
1142    /* Low band LAR coefficients. */
1143    for (n = 0; n < LPC_LOBAND_ORDER; n++, pos++, poss++) {
1144      LPCCoef[pos] = tmpcoeffs_s[poss] / LPC_LOBAND_SCALE;
1145      LPCCoef[pos] += WebRtcIsac_kLpcMeansShape[poss];
1146    }
1147
1148    /* High band LAR coefficients. */
1149    for (n = 0; n < LPC_HIBAND_ORDER; n++, pos++, poss++) {
1150      LPCCoef[pos] = tmpcoeffs_s[poss] / LPC_HIBAND_SCALE;
1151      LPCCoef[pos] += WebRtcIsac_kLpcMeansShape[poss];
1152    }
1153  }
1154}
1155
1156
1157void WebRtcIsac_EncodeLpcLb(double* LPCCoef_lo, double* LPCCoef_hi,
1158                            Bitstr* streamdata, IsacSaveEncoderData* encData) {
1159  double lars[KLT_ORDER_GAIN + KLT_ORDER_SHAPE];
1160  int k;
1161
1162  WebRtcIsac_Poly2Lar(LPCCoef_lo, ORDERLO, LPCCoef_hi, ORDERHI, SUBFRAMES,
1163                      lars);
1164  WebRtcIsac_EncodeLar(lars, streamdata, encData);
1165  WebRtcIsac_Lar2Poly(lars, LPCCoef_lo, ORDERLO, LPCCoef_hi, ORDERHI,
1166                      SUBFRAMES);
1167  /* Save data for creation of multiple bit streams (and transcoding). */
1168  for (k = 0; k < (ORDERLO + 1)*SUBFRAMES; k++) {
1169    encData->LPCcoeffs_lo[(ORDERLO + 1)*SUBFRAMES * encData->startIdx + k] =
1170        LPCCoef_lo[k];
1171  }
1172  for (k = 0; k < (ORDERHI + 1)*SUBFRAMES; k++) {
1173    encData->LPCcoeffs_hi[(ORDERHI + 1)*SUBFRAMES * encData->startIdx + k] =
1174        LPCCoef_hi[k];
1175  }
1176}
1177
1178
1179int16_t WebRtcIsac_EncodeLpcUB(double* lpcVecs, Bitstr* streamdata,
1180                               double* interpolLPCCoeff,
1181                               int16_t bandwidth,
1182                                     ISACUBSaveEncDataStruct* encData) {
1183  double    U[UB_LPC_ORDER * UB16_LPC_VEC_PER_FRAME];
1184  int     idx[UB_LPC_ORDER * UB16_LPC_VEC_PER_FRAME];
1185  int interpolCntr;
1186
1187  WebRtcIsac_Poly2LarUB(lpcVecs, bandwidth);
1188  WebRtcIsac_RemoveLarMean(lpcVecs, bandwidth);
1189  WebRtcIsac_DecorrelateIntraVec(lpcVecs, U, bandwidth);
1190  WebRtcIsac_DecorrelateInterVec(U, lpcVecs, bandwidth);
1191  WebRtcIsac_QuantizeUncorrLar(lpcVecs, idx, bandwidth);
1192
1193  WebRtcIsac_CorrelateInterVec(lpcVecs, U, bandwidth);
1194  WebRtcIsac_CorrelateIntraVec(U, lpcVecs, bandwidth);
1195  WebRtcIsac_AddLarMean(lpcVecs, bandwidth);
1196
1197  switch (bandwidth) {
1198    case isac12kHz: {
1199      /* Store the indices to be used for multiple encoding. */
1200      memcpy(encData->indexLPCShape, idx, UB_LPC_ORDER *
1201             UB_LPC_VEC_PER_FRAME * sizeof(int));
1202      WebRtcIsac_EncHistMulti(streamdata, idx, WebRtcIsac_kLpcShapeCdfMatUb12,
1203                              UB_LPC_ORDER * UB_LPC_VEC_PER_FRAME);
1204      for (interpolCntr = 0; interpolCntr < UB_INTERPOL_SEGMENTS;
1205          interpolCntr++) {
1206        WebRtcIsac_Lar2PolyInterpolUB(lpcVecs, interpolLPCCoeff,
1207                                      kLpcVecPerSegmentUb12 + 1);
1208        lpcVecs += UB_LPC_ORDER;
1209        interpolLPCCoeff += (kLpcVecPerSegmentUb12 * (UB_LPC_ORDER + 1));
1210      }
1211      break;
1212    }
1213    case isac16kHz: {
1214      /* Store the indices to be used for multiple encoding. */
1215      memcpy(encData->indexLPCShape, idx, UB_LPC_ORDER *
1216             UB16_LPC_VEC_PER_FRAME * sizeof(int));
1217      WebRtcIsac_EncHistMulti(streamdata, idx, WebRtcIsac_kLpcShapeCdfMatUb16,
1218                              UB_LPC_ORDER * UB16_LPC_VEC_PER_FRAME);
1219      for (interpolCntr = 0; interpolCntr < UB16_INTERPOL_SEGMENTS;
1220          interpolCntr++) {
1221        WebRtcIsac_Lar2PolyInterpolUB(lpcVecs, interpolLPCCoeff,
1222                                      kLpcVecPerSegmentUb16 + 1);
1223        lpcVecs += UB_LPC_ORDER;
1224        interpolLPCCoeff += (kLpcVecPerSegmentUb16 * (UB_LPC_ORDER + 1));
1225      }
1226      break;
1227    }
1228    default:
1229      return -1;
1230  }
1231  return 0;
1232}
1233
1234void WebRtcIsac_EncodeLpcGainLb(double* LPCCoef_lo, double* LPCCoef_hi,
1235                                Bitstr* streamdata,
1236                                IsacSaveEncoderData* encData) {
1237  int j, k, n, pos, pos2, posg, offsg, offs2;
1238  int index_g[KLT_ORDER_GAIN];
1239  int index_ovr_g[KLT_ORDER_GAIN];
1240  double tmpcoeffs_g[KLT_ORDER_GAIN];
1241  double tmpcoeffs2_g[KLT_ORDER_GAIN];
1242  double sum;
1243  /* log gains, mean removal and scaling */
1244  posg = 0;
1245  for (k = 0; k < SUBFRAMES; k++) {
1246    tmpcoeffs_g[posg] = log(LPCCoef_lo[(LPC_LOBAND_ORDER + 1) * k]);
1247    tmpcoeffs_g[posg] -= WebRtcIsac_kLpcMeansGain[posg];
1248    tmpcoeffs_g[posg] *= LPC_GAIN_SCALE;
1249    posg++;
1250    tmpcoeffs_g[posg] = log(LPCCoef_hi[(LPC_HIBAND_ORDER + 1) * k]);
1251    tmpcoeffs_g[posg] -= WebRtcIsac_kLpcMeansGain[posg];
1252    tmpcoeffs_g[posg] *= LPC_GAIN_SCALE;
1253    posg++;
1254  }
1255
1256  /* KLT  */
1257
1258  /* Left transform. */
1259  offsg = 0;
1260  for (j = 0; j < SUBFRAMES; j++) {
1261    posg = offsg;
1262    for (k = 0; k < LPC_GAIN_ORDER; k++) {
1263      sum = 0;
1264      pos = offsg;
1265      pos2 = k;
1266      for (n = 0; n < LPC_GAIN_ORDER; n++) {
1267        sum += tmpcoeffs_g[pos++] * WebRtcIsac_kKltT1Gain[pos2];
1268        pos2 += LPC_GAIN_ORDER;
1269      }
1270      tmpcoeffs2_g[posg++] = sum;
1271    }
1272    offsg += LPC_GAIN_ORDER;
1273  }
1274
1275  /* Right transform. */
1276  offsg = 0;
1277  offs2 = 0;
1278  for (j = 0; j < SUBFRAMES; j++) {
1279    posg = offsg;
1280    for (k = 0; k < LPC_GAIN_ORDER; k++) {
1281      sum = 0;
1282      pos = k;
1283      pos2 = offs2;
1284      for (n = 0; n < SUBFRAMES; n++) {
1285        sum += tmpcoeffs2_g[pos] * WebRtcIsac_kKltT2Gain[pos2++];
1286        pos += LPC_GAIN_ORDER;
1287      }
1288      tmpcoeffs_g[posg++] = sum;
1289    }
1290    offs2 += SUBFRAMES;
1291    offsg += LPC_GAIN_ORDER;
1292  }
1293
1294  /* Quantize coefficients. */
1295  for (k = 0; k < KLT_ORDER_GAIN; k++) {
1296    /* Get index. */
1297    pos2 = WebRtcIsac_lrint(tmpcoeffs_g[k] / KLT_STEPSIZE);
1298    index_g[k] = (pos2) + WebRtcIsac_kQKltQuantMinGain[k];
1299    if (index_g[k] < 0) {
1300      index_g[k] = 0;
1301    } else if (index_g[k] > WebRtcIsac_kQKltMaxIndGain[k]) {
1302      index_g[k] = WebRtcIsac_kQKltMaxIndGain[k];
1303    }
1304    index_ovr_g[k] = WebRtcIsac_kQKltOffsetGain[k] + index_g[k];
1305
1306    /* Find quantization levels for coefficients. */
1307    tmpcoeffs_g[k] = WebRtcIsac_kQKltLevelsGain[index_ovr_g[k]];
1308
1309    /* Save data for creation of multiple bit streams. */
1310    encData->LPCindex_g[KLT_ORDER_GAIN * encData->startIdx + k] = index_g[k];
1311  }
1312
1313  /* Entropy coding of quantization indices - gain. */
1314  WebRtcIsac_EncHistMulti(streamdata, index_g, WebRtcIsac_kQKltCdfPtrGain,
1315                          KLT_ORDER_GAIN);
1316
1317  /* Find quantization levels for coefficients. */
1318  /* Left transform. */
1319  offsg = 0;
1320  posg = 0;
1321  for (j = 0; j < SUBFRAMES; j++) {
1322    offs2 = 0;
1323    for (k = 0; k < LPC_GAIN_ORDER; k++) {
1324      sum = 0;
1325      pos = offsg;
1326      pos2 = offs2;
1327      for (n = 0; n < LPC_GAIN_ORDER; n++)
1328        sum += tmpcoeffs_g[pos++] * WebRtcIsac_kKltT1Gain[pos2++];
1329      tmpcoeffs2_g[posg++] = sum;
1330      offs2 += LPC_GAIN_ORDER;
1331    }
1332    offsg += LPC_GAIN_ORDER;
1333  }
1334
1335  /* Right transform, transpose matrix. */
1336  offsg = 0;
1337  posg = 0;
1338  for (j = 0; j < SUBFRAMES; j++) {
1339    posg = offsg;
1340    for (k = 0; k < LPC_GAIN_ORDER; k++) {
1341      sum = 0;
1342      pos = k;
1343      pos2 = j;
1344      for (n = 0; n < SUBFRAMES; n++) {
1345        sum += tmpcoeffs2_g[pos] * WebRtcIsac_kKltT2Gain[pos2];
1346        pos += LPC_GAIN_ORDER;
1347        pos2 += SUBFRAMES;
1348      }
1349      tmpcoeffs_g[posg++] = sum;
1350    }
1351    offsg += LPC_GAIN_ORDER;
1352  }
1353
1354
1355  /* Scaling, mean addition, and gain restoration. */
1356  posg = 0;
1357  for (k = 0; k < SUBFRAMES; k++) {
1358    sum = tmpcoeffs_g[posg] / LPC_GAIN_SCALE;
1359    sum += WebRtcIsac_kLpcMeansGain[posg];
1360    LPCCoef_lo[k * (LPC_LOBAND_ORDER + 1)] = exp(sum);
1361    pos++;
1362    posg++;
1363    sum = tmpcoeffs_g[posg] / LPC_GAIN_SCALE;
1364    sum += WebRtcIsac_kLpcMeansGain[posg];
1365    LPCCoef_hi[k * (LPC_HIBAND_ORDER + 1)] = exp(sum);
1366    pos++;
1367    posg++;
1368  }
1369
1370}
1371
1372void WebRtcIsac_EncodeLpcGainUb(double* lpGains, Bitstr* streamdata,
1373                                int* lpcGainIndex) {
1374  double U[UB_LPC_GAIN_DIM];
1375  int idx[UB_LPC_GAIN_DIM];
1376  WebRtcIsac_ToLogDomainRemoveMean(lpGains);
1377  WebRtcIsac_DecorrelateLPGain(lpGains, U);
1378  WebRtcIsac_QuantizeLpcGain(U, idx);
1379  /* Store the index for re-encoding for FEC. */
1380  memcpy(lpcGainIndex, idx, UB_LPC_GAIN_DIM * sizeof(int));
1381  WebRtcIsac_CorrelateLpcGain(U, lpGains);
1382  WebRtcIsac_AddMeanToLinearDomain(lpGains);
1383  WebRtcIsac_EncHistMulti(streamdata, idx, WebRtcIsac_kLpcGainCdfMat,
1384                          UB_LPC_GAIN_DIM);
1385}
1386
1387
1388void WebRtcIsac_StoreLpcGainUb(double* lpGains, Bitstr* streamdata) {
1389  double U[UB_LPC_GAIN_DIM];
1390  int idx[UB_LPC_GAIN_DIM];
1391  WebRtcIsac_ToLogDomainRemoveMean(lpGains);
1392  WebRtcIsac_DecorrelateLPGain(lpGains, U);
1393  WebRtcIsac_QuantizeLpcGain(U, idx);
1394  WebRtcIsac_EncHistMulti(streamdata, idx, WebRtcIsac_kLpcGainCdfMat,
1395                          UB_LPC_GAIN_DIM);
1396}
1397
1398
1399
1400int16_t WebRtcIsac_DecodeLpcGainUb(double* lpGains, Bitstr* streamdata) {
1401  double U[UB_LPC_GAIN_DIM];
1402  int idx[UB_LPC_GAIN_DIM];
1403  int err;
1404  err = WebRtcIsac_DecHistOneStepMulti(idx, streamdata,
1405                                       WebRtcIsac_kLpcGainCdfMat,
1406                                       WebRtcIsac_kLpcGainEntropySearch,
1407                                       UB_LPC_GAIN_DIM);
1408  if (err < 0) {
1409    return -1;
1410  }
1411  WebRtcIsac_DequantizeLpcGain(idx, U);
1412  WebRtcIsac_CorrelateLpcGain(U, lpGains);
1413  WebRtcIsac_AddMeanToLinearDomain(lpGains);
1414  return 0;
1415}
1416
1417
1418
1419/* decode & dequantize RC */
1420int WebRtcIsac_DecodeRc(Bitstr* streamdata, int16_t* RCQ15) {
1421  int k, err;
1422  int index[AR_ORDER];
1423
1424  /* entropy decoding of quantization indices */
1425  err = WebRtcIsac_DecHistOneStepMulti(index, streamdata,
1426                                       WebRtcIsac_kQArRcCdfPtr,
1427                                       WebRtcIsac_kQArRcInitIndex, AR_ORDER);
1428  if (err < 0)
1429    return err;
1430
1431  /* find quantization levels for reflection coefficients */
1432  for (k = 0; k < AR_ORDER; k++) {
1433    RCQ15[k] = *(WebRtcIsac_kQArRcLevelsPtr[k] + index[k]);
1434  }
1435  return 0;
1436}
1437
1438
1439/* quantize & code RC */
1440void WebRtcIsac_EncodeRc(int16_t* RCQ15, Bitstr* streamdata) {
1441  int k;
1442  int index[AR_ORDER];
1443
1444  /* quantize reflection coefficients (add noise feedback?) */
1445  for (k = 0; k < AR_ORDER; k++) {
1446    index[k] = WebRtcIsac_kQArRcInitIndex[k];
1447    // The safe-guards in following while conditions are to suppress gcc 4.8.3
1448    // warnings, Issue 2888. Otherwise, first and last elements of
1449    // |WebRtcIsac_kQArBoundaryLevels| are such that the following search
1450    // *never* cause an out-of-boundary read.
1451    if (RCQ15[k] > WebRtcIsac_kQArBoundaryLevels[index[k]]) {
1452      while (index[k] + 1 < NUM_AR_RC_QUANT_BAUNDARY &&
1453        RCQ15[k] > WebRtcIsac_kQArBoundaryLevels[index[k] + 1]) {
1454        index[k]++;
1455      }
1456    } else {
1457      while (index[k] > 0 &&
1458        RCQ15[k] < WebRtcIsac_kQArBoundaryLevels[--index[k]]) ;
1459    }
1460    RCQ15[k] = *(WebRtcIsac_kQArRcLevelsPtr[k] + index[k]);
1461  }
1462
1463  /* entropy coding of quantization indices */
1464  WebRtcIsac_EncHistMulti(streamdata, index, WebRtcIsac_kQArRcCdfPtr, AR_ORDER);
1465}
1466
1467
1468/* decode & dequantize squared Gain */
1469int WebRtcIsac_DecodeGain2(Bitstr* streamdata, int32_t* gainQ10) {
1470  int index, err;
1471
1472  /* entropy decoding of quantization index */
1473  err = WebRtcIsac_DecHistOneStepMulti(&index, streamdata,
1474                                       WebRtcIsac_kQGainCdf_ptr,
1475                                       WebRtcIsac_kQGainInitIndex, 1);
1476  if (err < 0) {
1477    return err;
1478  }
1479  /* find quantization level */
1480  *gainQ10 = WebRtcIsac_kQGain2Levels[index];
1481  return 0;
1482}
1483
1484
1485/* quantize & code squared Gain */
1486int WebRtcIsac_EncodeGain2(int32_t* gainQ10, Bitstr* streamdata) {
1487  int index;
1488
1489  /* find quantization index */
1490  index = WebRtcIsac_kQGainInitIndex[0];
1491  if (*gainQ10 > WebRtcIsac_kQGain2BoundaryLevels[index]) {
1492    while (*gainQ10 > WebRtcIsac_kQGain2BoundaryLevels[index + 1]) {
1493      index++;
1494    }
1495  } else {
1496    while (*gainQ10 < WebRtcIsac_kQGain2BoundaryLevels[--index]) ;
1497  }
1498  /* De-quantize */
1499  *gainQ10 = WebRtcIsac_kQGain2Levels[index];
1500
1501  /* entropy coding of quantization index */
1502  WebRtcIsac_EncHistMulti(streamdata, &index, WebRtcIsac_kQGainCdf_ptr, 1);
1503  return 0;
1504}
1505
1506
1507/* code and decode Pitch Gains and Lags functions */
1508
1509/* decode & dequantize Pitch Gains */
1510int WebRtcIsac_DecodePitchGain(Bitstr* streamdata,
1511                               int16_t* PitchGains_Q12) {
1512  int index_comb, err;
1513  const uint16_t* WebRtcIsac_kQPitchGainCdf_ptr[1];
1514
1515  /* Entropy decoding of quantization indices */
1516  *WebRtcIsac_kQPitchGainCdf_ptr = WebRtcIsac_kQPitchGainCdf;
1517  err = WebRtcIsac_DecHistBisectMulti(&index_comb, streamdata,
1518                                      WebRtcIsac_kQPitchGainCdf_ptr,
1519                                      WebRtcIsac_kQCdfTableSizeGain, 1);
1520  /* Error check, Q_mean_Gain.. tables are of size 144 */
1521  if ((err < 0) || (index_comb < 0) || (index_comb >= 144)) {
1522    return -ISAC_RANGE_ERROR_DECODE_PITCH_GAIN;
1523  }
1524  /* De-quantize back to pitch gains by table look-up. */
1525  PitchGains_Q12[0] = WebRtcIsac_kQMeanGain1Q12[index_comb];
1526  PitchGains_Q12[1] = WebRtcIsac_kQMeanGain2Q12[index_comb];
1527  PitchGains_Q12[2] = WebRtcIsac_kQMeanGain3Q12[index_comb];
1528  PitchGains_Q12[3] = WebRtcIsac_kQMeanGain4Q12[index_comb];
1529  return 0;
1530}
1531
1532
1533/* Quantize & code Pitch Gains. */
1534void WebRtcIsac_EncodePitchGain(int16_t* PitchGains_Q12,
1535                                Bitstr* streamdata,
1536                                IsacSaveEncoderData* encData) {
1537  int k, j;
1538  double C;
1539  double S[PITCH_SUBFRAMES];
1540  int index[3];
1541  int index_comb;
1542  const uint16_t* WebRtcIsac_kQPitchGainCdf_ptr[1];
1543  double PitchGains[PITCH_SUBFRAMES] = {0, 0, 0, 0};
1544
1545  /* Take the asin. */
1546  for (k = 0; k < PITCH_SUBFRAMES; k++) {
1547    PitchGains[k] = ((float)PitchGains_Q12[k]) / 4096;
1548    S[k] = asin(PitchGains[k]);
1549  }
1550
1551  /* Find quantization index; only for the first three
1552   * transform coefficients. */
1553  for (k = 0; k < 3; k++) {
1554    /*  transform */
1555    C = 0.0;
1556    for (j = 0; j < PITCH_SUBFRAMES; j++) {
1557      C += WebRtcIsac_kTransform[k][j] * S[j];
1558    }
1559    /* Quantize */
1560    index[k] = WebRtcIsac_lrint(C / PITCH_GAIN_STEPSIZE);
1561
1562    /* Check that the index is not outside the boundaries of the table. */
1563    if (index[k] < WebRtcIsac_kIndexLowerLimitGain[k]) {
1564      index[k] = WebRtcIsac_kIndexLowerLimitGain[k];
1565    } else if (index[k] > WebRtcIsac_kIndexUpperLimitGain[k]) {
1566      index[k] = WebRtcIsac_kIndexUpperLimitGain[k];
1567    }
1568    index[k] -= WebRtcIsac_kIndexLowerLimitGain[k];
1569  }
1570
1571  /* Calculate unique overall index. */
1572  index_comb = WebRtcIsac_kIndexMultsGain[0] * index[0] +
1573      WebRtcIsac_kIndexMultsGain[1] * index[1] + index[2];
1574
1575  /* unquantize back to pitch gains by table look-up */
1576  PitchGains_Q12[0] = WebRtcIsac_kQMeanGain1Q12[index_comb];
1577  PitchGains_Q12[1] = WebRtcIsac_kQMeanGain2Q12[index_comb];
1578  PitchGains_Q12[2] = WebRtcIsac_kQMeanGain3Q12[index_comb];
1579  PitchGains_Q12[3] = WebRtcIsac_kQMeanGain4Q12[index_comb];
1580
1581  /* entropy coding of quantization pitch gains */
1582  *WebRtcIsac_kQPitchGainCdf_ptr = WebRtcIsac_kQPitchGainCdf;
1583  WebRtcIsac_EncHistMulti(streamdata, &index_comb,
1584                          WebRtcIsac_kQPitchGainCdf_ptr, 1);
1585  encData->pitchGain_index[encData->startIdx] = index_comb;
1586}
1587
1588
1589
1590/* Pitch LAG */
1591/* Decode & de-quantize Pitch Lags. */
1592int WebRtcIsac_DecodePitchLag(Bitstr* streamdata, int16_t* PitchGain_Q12,
1593                              double* PitchLags) {
1594  int k, err;
1595  double StepSize;
1596  double C;
1597  int index[PITCH_SUBFRAMES];
1598  double mean_gain;
1599  const double* mean_val2, *mean_val3, *mean_val4;
1600  const int16_t* lower_limit;
1601  const uint16_t* init_index;
1602  const uint16_t* cdf_size;
1603  const uint16_t** cdf;
1604  double PitchGain[4] = {0, 0, 0, 0};
1605
1606  /* compute mean pitch gain */
1607  mean_gain = 0.0;
1608  for (k = 0; k < 4; k++) {
1609    PitchGain[k] = ((float)PitchGain_Q12[k]) / 4096;
1610    mean_gain += PitchGain[k];
1611  }
1612  mean_gain /= 4.0;
1613
1614  /* voicing classification. */
1615  if (mean_gain < 0.2) {
1616    StepSize = WebRtcIsac_kQPitchLagStepsizeLo;
1617    cdf = WebRtcIsac_kQPitchLagCdfPtrLo;
1618    cdf_size = WebRtcIsac_kQPitchLagCdfSizeLo;
1619    mean_val2 = WebRtcIsac_kQMeanLag2Lo;
1620    mean_val3 = WebRtcIsac_kQMeanLag3Lo;
1621    mean_val4 = WebRtcIsac_kQMeanLag4Lo;
1622    lower_limit = WebRtcIsac_kQIndexLowerLimitLagLo;
1623    init_index = WebRtcIsac_kQInitIndexLagLo;
1624  } else if (mean_gain < 0.4) {
1625    StepSize = WebRtcIsac_kQPitchLagStepsizeMid;
1626    cdf = WebRtcIsac_kQPitchLagCdfPtrMid;
1627    cdf_size = WebRtcIsac_kQPitchLagCdfSizeMid;
1628    mean_val2 = WebRtcIsac_kQMeanLag2Mid;
1629    mean_val3 = WebRtcIsac_kQMeanLag3Mid;
1630    mean_val4 = WebRtcIsac_kQMeanLag4Mid;
1631    lower_limit = WebRtcIsac_kQIndexLowerLimitLagMid;
1632    init_index = WebRtcIsac_kQInitIndexLagMid;
1633  } else {
1634    StepSize = WebRtcIsac_kQPitchLagStepsizeHi;
1635    cdf = WebRtcIsac_kQPitchLagCdfPtrHi;
1636    cdf_size = WebRtcIsac_kQPitchLagCdfSizeHi;
1637    mean_val2 = WebRtcIsac_kQMeanLag2Hi;
1638    mean_val3 = WebRtcIsac_kQMeanLag3Hi;
1639    mean_val4 = WebRtcIsac_kQMeanLag4Hi;
1640    lower_limit = WebRtcIsac_kQindexLowerLimitLagHi;
1641    init_index = WebRtcIsac_kQInitIndexLagHi;
1642  }
1643
1644  /* Entropy decoding of quantization indices. */
1645  err = WebRtcIsac_DecHistBisectMulti(index, streamdata, cdf, cdf_size, 1);
1646  if ((err < 0) || (index[0] < 0)) {
1647    return -ISAC_RANGE_ERROR_DECODE_PITCH_LAG;
1648  }
1649  err = WebRtcIsac_DecHistOneStepMulti(index + 1, streamdata, cdf + 1,
1650                                       init_index, 3);
1651  if (err < 0) {
1652    return -ISAC_RANGE_ERROR_DECODE_PITCH_LAG;
1653  }
1654
1655  /* Unquantize back to transform coefficients and do the inverse transform:
1656   * S = T'*C. */
1657  C = (index[0] + lower_limit[0]) * StepSize;
1658  for (k = 0; k < PITCH_SUBFRAMES; k++) {
1659    PitchLags[k] = WebRtcIsac_kTransformTranspose[k][0] * C;
1660  }
1661  C = mean_val2[index[1]];
1662  for (k = 0; k < PITCH_SUBFRAMES; k++) {
1663    PitchLags[k] += WebRtcIsac_kTransformTranspose[k][1] * C;
1664  }
1665  C = mean_val3[index[2]];
1666  for (k = 0; k < PITCH_SUBFRAMES; k++) {
1667    PitchLags[k] += WebRtcIsac_kTransformTranspose[k][2] * C;
1668  }
1669  C = mean_val4[index[3]];
1670  for (k = 0; k < PITCH_SUBFRAMES; k++) {
1671    PitchLags[k] += WebRtcIsac_kTransformTranspose[k][3] * C;
1672  }
1673  return 0;
1674}
1675
1676
1677
1678/* Quantize & code pitch lags. */
1679void WebRtcIsac_EncodePitchLag(double* PitchLags, int16_t* PitchGain_Q12,
1680                               Bitstr* streamdata,
1681                               IsacSaveEncoderData* encData) {
1682  int k, j;
1683  double StepSize;
1684  double C;
1685  int index[PITCH_SUBFRAMES];
1686  double mean_gain;
1687  const double* mean_val2, *mean_val3, *mean_val4;
1688  const int16_t* lower_limit, *upper_limit;
1689  const uint16_t** cdf;
1690  double PitchGain[4] = {0, 0, 0, 0};
1691
1692  /* compute mean pitch gain */
1693  mean_gain = 0.0;
1694  for (k = 0; k < 4; k++) {
1695    PitchGain[k] = ((float)PitchGain_Q12[k]) / 4096;
1696    mean_gain += PitchGain[k];
1697  }
1698  mean_gain /= 4.0;
1699
1700  /* Save data for creation of multiple bit streams */
1701  encData->meanGain[encData->startIdx] = mean_gain;
1702
1703  /* Voicing classification. */
1704  if (mean_gain < 0.2) {
1705    StepSize = WebRtcIsac_kQPitchLagStepsizeLo;
1706    cdf = WebRtcIsac_kQPitchLagCdfPtrLo;
1707    mean_val2 = WebRtcIsac_kQMeanLag2Lo;
1708    mean_val3 = WebRtcIsac_kQMeanLag3Lo;
1709    mean_val4 = WebRtcIsac_kQMeanLag4Lo;
1710    lower_limit = WebRtcIsac_kQIndexLowerLimitLagLo;
1711    upper_limit = WebRtcIsac_kQIndexUpperLimitLagLo;
1712  } else if (mean_gain < 0.4) {
1713    StepSize = WebRtcIsac_kQPitchLagStepsizeMid;
1714    cdf = WebRtcIsac_kQPitchLagCdfPtrMid;
1715    mean_val2 = WebRtcIsac_kQMeanLag2Mid;
1716    mean_val3 = WebRtcIsac_kQMeanLag3Mid;
1717    mean_val4 = WebRtcIsac_kQMeanLag4Mid;
1718    lower_limit = WebRtcIsac_kQIndexLowerLimitLagMid;
1719    upper_limit = WebRtcIsac_kQIndexUpperLimitLagMid;
1720  } else {
1721    StepSize = WebRtcIsac_kQPitchLagStepsizeHi;
1722    cdf = WebRtcIsac_kQPitchLagCdfPtrHi;
1723    mean_val2 = WebRtcIsac_kQMeanLag2Hi;
1724    mean_val3 = WebRtcIsac_kQMeanLag3Hi;
1725    mean_val4 = WebRtcIsac_kQMeanLag4Hi;
1726    lower_limit = WebRtcIsac_kQindexLowerLimitLagHi;
1727    upper_limit = WebRtcIsac_kQindexUpperLimitLagHi;
1728  }
1729
1730  /* find quantization index */
1731  for (k = 0; k < 4; k++) {
1732    /*  transform */
1733    C = 0.0;
1734    for (j = 0; j < PITCH_SUBFRAMES; j++) {
1735      C += WebRtcIsac_kTransform[k][j] * PitchLags[j];
1736    }
1737    /* quantize */
1738    index[k] = WebRtcIsac_lrint(C / StepSize);
1739
1740    /* check that the index is not outside the boundaries of the table */
1741    if (index[k] < lower_limit[k]) {
1742      index[k] = lower_limit[k];
1743    } else if (index[k] > upper_limit[k]) index[k] = upper_limit[k]; {
1744      index[k] -= lower_limit[k];
1745    }
1746    /* Save data for creation of multiple bit streams */
1747    encData->pitchIndex[PITCH_SUBFRAMES * encData->startIdx + k] = index[k];
1748  }
1749
1750  /* Un-quantize back to transform coefficients and do the inverse transform:
1751   * S = T'*C */
1752  C = (index[0] + lower_limit[0]) * StepSize;
1753  for (k = 0; k < PITCH_SUBFRAMES; k++) {
1754    PitchLags[k] = WebRtcIsac_kTransformTranspose[k][0] * C;
1755  }
1756  C = mean_val2[index[1]];
1757  for (k = 0; k < PITCH_SUBFRAMES; k++) {
1758    PitchLags[k] += WebRtcIsac_kTransformTranspose[k][1] * C;
1759  }
1760  C = mean_val3[index[2]];
1761  for (k = 0; k < PITCH_SUBFRAMES; k++) {
1762    PitchLags[k] += WebRtcIsac_kTransformTranspose[k][2] * C;
1763  }
1764  C = mean_val4[index[3]];
1765  for (k = 0; k < PITCH_SUBFRAMES; k++) {
1766    PitchLags[k] += WebRtcIsac_kTransformTranspose[k][3] * C;
1767  }
1768  /* entropy coding of quantization pitch lags */
1769  WebRtcIsac_EncHistMulti(streamdata, index, cdf, PITCH_SUBFRAMES);
1770}
1771
1772
1773
1774/* Routines for in-band signaling of bandwidth estimation */
1775/* Histograms based on uniform distribution of indices */
1776/* Move global variables later! */
1777
1778
1779/* cdf array for frame length indicator */
1780const uint16_t WebRtcIsac_kFrameLengthCdf[4] = {
1781    0, 21845, 43690, 65535 };
1782
1783/* pointer to cdf array for frame length indicator */
1784const uint16_t* WebRtcIsac_kFrameLengthCdf_ptr[1] = {
1785    WebRtcIsac_kFrameLengthCdf };
1786
1787/* initial cdf index for decoder of frame length indicator */
1788const uint16_t WebRtcIsac_kFrameLengthInitIndex[1] = { 1 };
1789
1790
1791int WebRtcIsac_DecodeFrameLen(Bitstr* streamdata, int16_t* framesamples) {
1792  int frame_mode, err;
1793  err = 0;
1794  /* entropy decoding of frame length [1:30ms,2:60ms] */
1795  err = WebRtcIsac_DecHistOneStepMulti(&frame_mode, streamdata,
1796                                       WebRtcIsac_kFrameLengthCdf_ptr,
1797                                       WebRtcIsac_kFrameLengthInitIndex, 1);
1798  if (err < 0)
1799    return -ISAC_RANGE_ERROR_DECODE_FRAME_LENGTH;
1800
1801  switch (frame_mode) {
1802    case 1:
1803      *framesamples = 480; /* 30ms */
1804      break;
1805    case 2:
1806      *framesamples = 960; /* 60ms */
1807      break;
1808    default:
1809      err = -ISAC_DISALLOWED_FRAME_MODE_DECODER;
1810  }
1811  return err;
1812}
1813
1814int WebRtcIsac_EncodeFrameLen(int16_t framesamples, Bitstr* streamdata) {
1815  int frame_mode, status;
1816
1817  status = 0;
1818  frame_mode = 0;
1819  /* entropy coding of frame length [1:480 samples,2:960 samples] */
1820  switch (framesamples) {
1821    case 480:
1822      frame_mode = 1;
1823      break;
1824    case 960:
1825      frame_mode = 2;
1826      break;
1827    default:
1828      status = - ISAC_DISALLOWED_FRAME_MODE_ENCODER;
1829  }
1830
1831  if (status < 0)
1832    return status;
1833
1834  WebRtcIsac_EncHistMulti(streamdata, &frame_mode,
1835                          WebRtcIsac_kFrameLengthCdf_ptr, 1);
1836  return status;
1837}
1838
1839/* cdf array for estimated bandwidth */
1840static const uint16_t kBwCdf[25] = {
1841    0, 2731, 5461, 8192, 10923, 13653, 16384, 19114, 21845, 24576, 27306, 30037,
1842    32768, 35498, 38229, 40959, 43690, 46421, 49151, 51882, 54613, 57343, 60074,
1843    62804, 65535 };
1844
1845/* pointer to cdf array for estimated bandwidth */
1846static const uint16_t* kBwCdfPtr[1] = { kBwCdf };
1847
1848/* initial cdf index for decoder of estimated bandwidth*/
1849static const uint16_t kBwInitIndex[1] = { 7 };
1850
1851
1852int WebRtcIsac_DecodeSendBW(Bitstr* streamdata, int16_t* BWno) {
1853  int BWno32, err;
1854
1855  /* entropy decoding of sender's BW estimation [0..23] */
1856  err = WebRtcIsac_DecHistOneStepMulti(&BWno32, streamdata, kBwCdfPtr,
1857                                       kBwInitIndex, 1);
1858  if (err < 0) {
1859    return -ISAC_RANGE_ERROR_DECODE_BANDWIDTH;
1860  }
1861  *BWno = (int16_t)BWno32;
1862  return err;
1863}
1864
1865void WebRtcIsac_EncodeReceiveBw(int* BWno, Bitstr* streamdata) {
1866  /* entropy encoding of receiver's BW estimation [0..23] */
1867  WebRtcIsac_EncHistMulti(streamdata, BWno, kBwCdfPtr, 1);
1868}
1869
1870
1871/* estimate code length of LPC Coef */
1872void WebRtcIsac_TranscodeLPCCoef(double* LPCCoef_lo, double* LPCCoef_hi,
1873                                 int* index_g) {
1874  int j, k, n, pos, pos2, posg, offsg, offs2;
1875  int index_ovr_g[KLT_ORDER_GAIN];
1876  double tmpcoeffs_g[KLT_ORDER_GAIN];
1877  double tmpcoeffs2_g[KLT_ORDER_GAIN];
1878  double sum;
1879
1880  /* log gains, mean removal and scaling */
1881  posg = 0;
1882  for (k = 0; k < SUBFRAMES; k++) {
1883    tmpcoeffs_g[posg] = log(LPCCoef_lo[(LPC_LOBAND_ORDER + 1) * k]);
1884    tmpcoeffs_g[posg] -= WebRtcIsac_kLpcMeansGain[posg];
1885    tmpcoeffs_g[posg] *= LPC_GAIN_SCALE;
1886    posg++;
1887    tmpcoeffs_g[posg] = log(LPCCoef_hi[(LPC_HIBAND_ORDER + 1) * k]);
1888    tmpcoeffs_g[posg] -= WebRtcIsac_kLpcMeansGain[posg];
1889    tmpcoeffs_g[posg] *= LPC_GAIN_SCALE;
1890    posg++;
1891  }
1892
1893  /* KLT  */
1894
1895  /* Left transform. */
1896  offsg = 0;
1897  for (j = 0; j < SUBFRAMES; j++) {
1898    posg = offsg;
1899    for (k = 0; k < LPC_GAIN_ORDER; k++) {
1900      sum = 0;
1901      pos = offsg;
1902      pos2 = k;
1903      for (n = 0; n < LPC_GAIN_ORDER; n++) {
1904        sum += tmpcoeffs_g[pos++] * WebRtcIsac_kKltT1Gain[pos2];
1905        pos2 += LPC_GAIN_ORDER;
1906      }
1907      tmpcoeffs2_g[posg++] = sum;
1908    }
1909    offsg += LPC_GAIN_ORDER;
1910  }
1911
1912  /* Right transform. */
1913  offsg = 0;
1914  offs2 = 0;
1915  for (j = 0; j < SUBFRAMES; j++) {
1916    posg = offsg;
1917    for (k = 0; k < LPC_GAIN_ORDER; k++) {
1918      sum = 0;
1919      pos = k;
1920      pos2 = offs2;
1921      for (n = 0; n < SUBFRAMES; n++) {
1922        sum += tmpcoeffs2_g[pos] * WebRtcIsac_kKltT2Gain[pos2++];
1923        pos += LPC_GAIN_ORDER;
1924      }
1925      tmpcoeffs_g[posg++] = sum;
1926    }
1927    offs2 += SUBFRAMES;
1928    offsg += LPC_GAIN_ORDER;
1929  }
1930
1931
1932  /* quantize coefficients */
1933  for (k = 0; k < KLT_ORDER_GAIN; k++) {
1934    /* Get index. */
1935    pos2 = WebRtcIsac_lrint(tmpcoeffs_g[k] / KLT_STEPSIZE);
1936    index_g[k] = (pos2) + WebRtcIsac_kQKltQuantMinGain[k];
1937    if (index_g[k] < 0) {
1938      index_g[k] = 0;
1939    } else if (index_g[k] > WebRtcIsac_kQKltMaxIndGain[k]) {
1940      index_g[k] = WebRtcIsac_kQKltMaxIndGain[k];
1941    }
1942    index_ovr_g[k] = WebRtcIsac_kQKltOffsetGain[k] + index_g[k];
1943
1944    /* find quantization levels for coefficients */
1945    tmpcoeffs_g[k] = WebRtcIsac_kQKltLevelsGain[index_ovr_g[k]];
1946  }
1947}
1948
1949
1950/* Decode & de-quantize LPC Coefficients. */
1951int WebRtcIsac_DecodeLpcCoefUB(Bitstr* streamdata, double* lpcVecs,
1952                               double* percepFilterGains,
1953                               int16_t bandwidth) {
1954  int  index_s[KLT_ORDER_SHAPE];
1955
1956  double U[UB_LPC_ORDER * UB16_LPC_VEC_PER_FRAME];
1957  int err;
1958
1959  /* Entropy decoding of quantization indices. */
1960  switch (bandwidth) {
1961    case isac12kHz: {
1962      err = WebRtcIsac_DecHistOneStepMulti(
1963          index_s, streamdata, WebRtcIsac_kLpcShapeCdfMatUb12,
1964          WebRtcIsac_kLpcShapeEntropySearchUb12, UB_LPC_ORDER *
1965          UB_LPC_VEC_PER_FRAME);
1966      break;
1967    }
1968    case isac16kHz: {
1969      err = WebRtcIsac_DecHistOneStepMulti(
1970          index_s, streamdata, WebRtcIsac_kLpcShapeCdfMatUb16,
1971          WebRtcIsac_kLpcShapeEntropySearchUb16, UB_LPC_ORDER *
1972          UB16_LPC_VEC_PER_FRAME);
1973      break;
1974    }
1975    default:
1976      return -1;
1977  }
1978
1979  if (err < 0) {
1980    return err;
1981  }
1982
1983  WebRtcIsac_DequantizeLpcParam(index_s, lpcVecs, bandwidth);
1984  WebRtcIsac_CorrelateInterVec(lpcVecs, U, bandwidth);
1985  WebRtcIsac_CorrelateIntraVec(U, lpcVecs, bandwidth);
1986  WebRtcIsac_AddLarMean(lpcVecs, bandwidth);
1987  WebRtcIsac_DecodeLpcGainUb(percepFilterGains, streamdata);
1988
1989  if (bandwidth == isac16kHz) {
1990    /* Decode another set of Gains. */
1991    WebRtcIsac_DecodeLpcGainUb(&percepFilterGains[SUBFRAMES], streamdata);
1992  }
1993  return 0;
1994}
1995
1996int16_t WebRtcIsac_EncodeBandwidth(enum ISACBandwidth bandwidth,
1997                                   Bitstr* streamData) {
1998  int bandwidthMode;
1999  switch (bandwidth) {
2000    case isac12kHz: {
2001      bandwidthMode = 0;
2002      break;
2003    }
2004    case isac16kHz: {
2005      bandwidthMode = 1;
2006      break;
2007    }
2008    default:
2009      return -ISAC_DISALLOWED_ENCODER_BANDWIDTH;
2010  }
2011  WebRtcIsac_EncHistMulti(streamData, &bandwidthMode, kOneBitEqualProbCdf_ptr,
2012                          1);
2013  return 0;
2014}
2015
2016int16_t WebRtcIsac_DecodeBandwidth(Bitstr* streamData,
2017                                   enum ISACBandwidth* bandwidth) {
2018  int bandwidthMode;
2019  if (WebRtcIsac_DecHistOneStepMulti(&bandwidthMode, streamData,
2020                                     kOneBitEqualProbCdf_ptr,
2021                                     kOneBitEqualProbInitIndex, 1) < 0) {
2022    return -ISAC_RANGE_ERROR_DECODE_BANDWITH;
2023  }
2024  switch (bandwidthMode) {
2025    case 0: {
2026      *bandwidth = isac12kHz;
2027      break;
2028    }
2029    case 1: {
2030      *bandwidth = isac16kHz;
2031      break;
2032    }
2033    default:
2034      return -ISAC_DISALLOWED_BANDWIDTH_MODE_DECODER;
2035  }
2036  return 0;
2037}
2038
2039int16_t WebRtcIsac_EncodeJitterInfo(int32_t jitterIndex,
2040                                    Bitstr* streamData) {
2041  /* This is to avoid LINUX warning until we change 'int' to 'Word32'. */
2042  int intVar;
2043
2044  if ((jitterIndex < 0) || (jitterIndex > 1)) {
2045    return -1;
2046  }
2047  intVar = (int)(jitterIndex);
2048  /* Use the same CDF table as for bandwidth
2049   * both take two values with equal probability.*/
2050  WebRtcIsac_EncHistMulti(streamData, &intVar, kOneBitEqualProbCdf_ptr, 1);
2051  return 0;
2052}
2053
2054int16_t WebRtcIsac_DecodeJitterInfo(Bitstr* streamData,
2055                                    int32_t* jitterInfo) {
2056  int intVar;
2057  /* Use the same CDF table as for bandwidth
2058   * both take two values with equal probability. */
2059  if (WebRtcIsac_DecHistOneStepMulti(&intVar, streamData,
2060                                     kOneBitEqualProbCdf_ptr,
2061                                     kOneBitEqualProbInitIndex, 1) < 0) {
2062    return -ISAC_RANGE_ERROR_DECODE_BANDWITH;
2063  }
2064  *jitterInfo = (int16_t)(intVar);
2065  return 0;
2066}
2067