1 /*
2  *  Copyright (c) 2011 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 
13  iLBC Speech Coder ANSI-C Source Code
14 
15  WebRtcIlbcfix_Poly2Lsp.c
16 
17 ******************************************************************/
18 
19 #include "defines.h"
20 #include "constants.h"
21 #include "chebyshev.h"
22 
23 /*----------------------------------------------------------------*
24  * conversion from lpc coefficients to lsp coefficients
25  * function is only for 10:th order LPC
26  *---------------------------------------------------------------*/
27 
WebRtcIlbcfix_Poly2Lsp(int16_t * a,int16_t * lsp,int16_t * old_lsp)28 void WebRtcIlbcfix_Poly2Lsp(
29     int16_t *a,  /* (o) A coefficients in Q12 */
30     int16_t *lsp, /* (i) LSP coefficients in Q15 */
31     int16_t *old_lsp /* (i) old LSP coefficients that are used if the new
32                               coefficients turn out to be unstable */
33                             ) {
34   int16_t f[2][6]; /* f[0][] represents f1 and f[1][] represents f2 */
35   int16_t *a_i_ptr, *a_10mi_ptr;
36   int16_t *f1ptr, *f2ptr;
37   int32_t tmpW32;
38   int16_t x, y, xlow, ylow, xmid, ymid, xhigh, yhigh, xint;
39   int16_t shifts, sign;
40   int i, j;
41   int foundFreqs;
42   int fi_select;
43 
44   /*
45      Calculate the two polynomials f1(z) and f2(z)
46      (the sum and the diff polynomial)
47      f1[0] = f2[0] = 1.0;
48      f1[i+1] = a[i+1] + a[10-i] - f1[i];
49      f2[i+1] = a[i+1] - a[10-i] - f1[i];
50   */
51 
52   a_i_ptr = a + 1;
53   a_10mi_ptr = a + 10;
54   f1ptr = f[0];
55   f2ptr = f[1];
56   (*f1ptr) = 1024; /* 1.0 in Q10 */
57   (*f2ptr) = 1024; /* 1.0 in Q10 */
58   for (i = 0; i < 5; i++) {
59     *(f1ptr + 1) =
60         (int16_t)((((int32_t)(*a_i_ptr) + *a_10mi_ptr) >> 2) - *f1ptr);
61     *(f2ptr + 1) =
62         (int16_t)((((int32_t)(*a_i_ptr) - *a_10mi_ptr) >> 2) + *f2ptr);
63     a_i_ptr++;
64     a_10mi_ptr--;
65     f1ptr++;
66     f2ptr++;
67   }
68 
69   /*
70     find the LSPs using the Chebychev pol. evaluation
71   */
72 
73   fi_select = 0; /* selector between f1 and f2, start with f1 */
74 
75   foundFreqs = 0;
76 
77   xlow = WebRtcIlbcfix_kCosGrid[0];
78   ylow = WebRtcIlbcfix_Chebyshev(xlow, f[fi_select]);
79 
80   /*
81      Iterate until all the 10 LSP's have been found or
82      all the grid points have been tried. If the 10 LSP's can
83      not be found, set the LSP vector to previous LSP
84   */
85 
86   for (j = 1; j < COS_GRID_POINTS && foundFreqs < 10; j++) {
87     xhigh = xlow;
88     yhigh = ylow;
89     xlow = WebRtcIlbcfix_kCosGrid[j];
90     ylow = WebRtcIlbcfix_Chebyshev(xlow, f[fi_select]);
91 
92     if (ylow * yhigh <= 0) {
93       /* Run 4 times to reduce the interval */
94       for (i = 0; i < 4; i++) {
95         /* xmid =(xlow + xhigh)/2 */
96         xmid = (xlow >> 1) + (xhigh >> 1);
97         ymid = WebRtcIlbcfix_Chebyshev(xmid, f[fi_select]);
98 
99         if (ylow * ymid <= 0) {
100           yhigh = ymid;
101           xhigh = xmid;
102         } else {
103           ylow = ymid;
104           xlow = xmid;
105         }
106       }
107 
108       /*
109         Calculater xint by linear interpolation:
110         xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow);
111       */
112 
113       x = xhigh - xlow;
114       y = yhigh - ylow;
115 
116       if (y == 0) {
117         xint = xlow;
118       } else {
119         sign = y;
120         y = WEBRTC_SPL_ABS_W16(y);
121         shifts = (int16_t)WebRtcSpl_NormW32(y)-16;
122         y <<= shifts;
123         y = (int16_t)WebRtcSpl_DivW32W16(536838144, y); /* 1/(yhigh-ylow) */
124 
125         tmpW32 = (x * y) >> (19 - shifts);
126 
127         /* y=(xhigh-xlow)/(yhigh-ylow) */
128         y = (int16_t)(tmpW32&0xFFFF);
129 
130         if (sign < 0) {
131           y = -y;
132         }
133         /* tmpW32 = ylow*(xhigh-xlow)/(yhigh-ylow) */
134         tmpW32 = (ylow * y) >> 10;
135         xint = xlow-(int16_t)(tmpW32&0xFFFF);
136       }
137 
138       /* Store the calculated lsp */
139       lsp[foundFreqs] = (int16_t)xint;
140       foundFreqs++;
141 
142       /* if needed, set xlow and ylow for next recursion */
143       if (foundFreqs<10) {
144         xlow = xint;
145         /* Swap between f1 and f2 (f[0][] and f[1][]) */
146         fi_select = ((fi_select+1)&0x1);
147 
148         ylow = WebRtcIlbcfix_Chebyshev(xlow, f[fi_select]);
149       }
150     }
151   }
152 
153   /* Check if M roots found, if not then use the old LSP */
154   if (foundFreqs < 10) {
155     WEBRTC_SPL_MEMCPY_W16(lsp, old_lsp, 10);
156   }
157   return;
158 }
159