Bug Summary

File:libs/broadvoice/src/floating/common/a2lsp.c
Location:line 191, column 21
Description:Value stored to 'xhigh' is never read

Annotated Source Code

1/*
2 * broadvoice - a library for the BroadVoice 16 and 32 codecs
3 *
4 * a2lsp.c -
5 *
6 * Adapted by Steve Underwood <steveu@coppice.org> from code which is
7 * Copyright 2000-2009 Broadcom Corporation
8 *
9 * All rights reserved.
10 *
11 * This program is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU Lesser General Public License version 2.1,
13 * as published by the Free Software Foundation.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU Lesser General Public License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with this program; if not, write to the Free Software
22 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
23 *
24 * $Id: a2lsp.c,v 1.1.1.1 2009/11/19 12:10:48 steveu Exp $
25 */
26
27/*! \file */
28
29#if defined(HAVE_CONFIG_H1)
30#include "config.h"
31#endif
32
33#include <stdio.h>
34#include <math.h>
35#include "typedef.h"
36#include "bvcommon.h"
37
38#define PI3.14159265358979 3.14159265358979
39
40#define NAB((8 >> 1) + 1) ((LPCO8 >> 1) + 1)
41
42static Float FNevChebP(Float x, Float *c, int nd2);
43
44#define NBIS4 4 /* number of bisections */
45
46/*----------------------------------------------------------------------------
47* a2lsp - Convert predictor coefficients to line spectral pairs
48*
49* Description:
50* The transfer function of the prediction error filter is formed from the
51* predictor coefficients. This polynomial is transformed into two reciprocal
52* polynomials having roots on the unit circle. The roots of these polynomials
53* interlace. It is these roots that determine the line spectral pairs.
54* The two reciprocal polynomials are expressed as series expansions in
55* Chebyshev polynomials with roots in the range -1 to +1. The inverse cosine
56* of the roots of the Chebyshev polynomial expansion gives the line spectral
57* pairs. If np line spectral pairs are not found, this routine
58* stops with an error message. This error occurs if the input coefficients
59* do not give a prediction error filter with minimum phase.
60*
61* Line spectral pairs and predictor coefficients are usually expressed
62* algebraically as vectors.
63* lsp[0] first (lowest frequency) line spectral pair
64* lsp[i] 1 <= i < np
65* pc[0]=1.0 predictor coefficient corresponding to lag 0
66* pc[i] 1 <= 1 <= np
67*
68* Parameters:
69* -> Float pc[]
70* Vector of predictor coefficients (Np+1 values). These are the
71* coefficients of the predictor filter, with pc[0] being the predictor
72* coefficient corresponding to lag 0, and pc[Np] corresponding to lag Np.
73* The predictor coeffs. must correspond to a minimum phase prediction
74* error filter.
75* <- Float lsp[]
76* Array of Np line spectral pairss (in ascending order). Each line
77* spectral pair lies in the range 0 to pi.
78* -> int Np
79* Number of coefficients (at most LPCO = 50)
80*----------------------------------------------------------------------------
81*/
82
83void a2lsp(Float pc[], /* (i) input the np+1 predictor coeff. */
84 Float lsp[], /* (o) line spectral pairs */
85 Float old_lsp[]) /* (i/o) old lsp[] (in case not found 10 roots) */
86{
87 Float fa[NAB((8 >> 1) + 1)], fb[NAB((8 >> 1) + 1)];
88 Float ta[NAB((8 >> 1) + 1)], tb[NAB((8 >> 1) + 1)];
89 Float *t;
90 Float xlow, xmid, xhigh;
91 Float ylow, ymid, yhigh;
92 Float xroot;
93 Float dx;
94 int i, j, nf, nd2, nab = ((LPCO8>>1) + 1), ngrd;
95
96 fb[0] = fa[0] = 1.0;
97 for (i = 1, j = LPCO8; i <= (LPCO8/2); i++, j--)
98 {
99 fa[i] = pc[i] + pc[j] - fa[i-1];
100 fb[i] = pc[i] - pc[j] + fb[i-1];
101 }
102
103 nd2 = LPCO8/2;
104
105 /*
106 * To look for roots on the unit circle, Ga(D) and Gb(D) are evaluated for
107 * D=exp(jw). Since Gz(D) and Gb(D) are symmetric, they can be expressed in
108 * terms of a series in cos(nw) for D on the unit circle. Since M is odd and
109 * D=exp(jw)
110 *
111 * M-1 n
112 * Ga(D) = SUM fa(n) D (symmetric, fa(n) = fa(M-1-n))
113 * n=0
114 * Mh-1
115 * = exp(j Mh w) [ f1(Mh) + 2 SUM fa(n) cos((Mh-n)w) ]
116 * n=0
117 * Mh
118 * = exp(j Mh w) SUM ta(n) cos(nw),
119 * n=0
120 *
121 * where Mh=(M-1)/2=Nc-1. The Nc=Mh+1 coefficients ta(n) are defined as
122 *
123 * ta(n) = fa(Nc-1), n=0,
124 * = 2 fa(Nc-1-n), n=1,...,Nc-1.
125 * The next step is to identify cos(nw) with the Chebyshev polynomial T(n,x).
126 * The Chebyshev polynomials satisfy the relationship T(n,cos(w)) = cos(nw).
127 * Omitting the exponential delay term, the series expansion in terms of
128 * Chebyshev polynomials is
129 *
130 * Nc-1
131 * Ta(x) = SUM ta(n) T(n,x)
132 * n=0
133 *
134 * The domain of Ta(x) is -1 < x < +1. For a given root of Ta(x), say x0,
135 * the corresponding position of the root of Fa(D) on the unit circle is
136 * exp(j arccos(x0)).
137 */
138 ta[0] = fa[nab-1];
139 tb[0] = fb[nab-1];
140 for (i = 1, j = nab - 2; i < nab; ++i, --j)
141 {
142 ta[i] = 2.0 * fa[j];
143 tb[i] = 2.0 * fb[j];
144 }
145
146 /*
147 * To find the roots, we sample the polynomials Ta(x) and Tb(x) looking for
148 * sign changes. An interval containing a root is successively bisected to
149 * narrow the interval and then linear interpolation is used to estimate the
150 * root. For a given root at x0, the line spectral pair is w0=acos(x0).
151 *
152 * Since the roots of the two polynomials interlace, the search for roots
153 * alternates between the polynomials Ta(x) and Tb(x). The sampling interval
154 * must be small enough to avoid having two cancelling sign changes in the
155 * same interval. The sampling (grid) points were trained from a large amount
156 * of LSP vectors derived with high accuracy and stored in a table.
157 */
158
159 nf = 0;
160 t = ta;
161 xroot = 2.0;
162 ngrd = 0;
163 xlow = grid[0];
164 ylow = FNevChebP(xlow, t, nd2);
165
166
167 /* Root search loop */
168 while (ngrd<(Ngrd60-1) && nf < LPCO8)
169 {
170
171 /* New trial point */
172 ngrd++;
173 xhigh = xlow;
174 yhigh = ylow;
175 xlow = grid[ngrd];
176 ylow = FNevChebP(xlow, t, nd2);
177
178 if (ylow * yhigh <= 0.0)
179 {
180
181 /* Bisections of the interval containing a sign change */
182 dx = xhigh - xlow;
183 for (i = 1; i <= NBIS4; ++i)
184 {
185 dx = 0.5 * dx;
186 xmid = xlow + dx;
187 ymid = FNevChebP(xmid, t, nd2);
188 if (ylow * ymid <= 0.0)
189 {
190 yhigh = ymid;
191 xhigh = xmid;
Value stored to 'xhigh' is never read
192 }
193 else
194 {
195 ylow = ymid;
196 xlow = xmid;
197 }
198 }
199
200 /*
201 * Linear interpolation in the subinterval with a sign change
202 * (take care if yhigh=ylow=0)
203 */
204 if (yhigh != ylow)
205 xmid = xlow + dx * ylow / (ylow - yhigh);
206 else
207 xmid = xlow + dx;
208
209 /* New root position */
210 lsp[nf] = acos(xmid)/PI3.14159265358979;
211 ++nf;
212
213 /* Start the search for the roots of the next polynomial at the estimated
214 * location of the root just found. We have to catch the case that the
215 * two polynomials have roots at the same place to avoid getting stuck at
216 * that root.
217 */
218 if (xmid >= xroot)
219 {
220 xmid = xlow - dx;
221 }
222 xroot = xmid;
223 if (t == ta)
224 t = tb;
225 else
226 t = ta;
227 xlow = xmid;
228 ylow = FNevChebP(xlow, t, nd2);
229 }
230 }
231
232 if (nf != LPCO8)
233 {
234 /* LPCO roots have not been found */
235 printf("\nWARNING: a2lsp failed to find all lsp nf=%d LPCO=%d\n", nf, LPCO8);
236 for (i = 0; i < LPCO8; i++)
237 lsp[i] = old_lsp[i];
238 }
239 else
240 {
241 /* Update LSP of previous frame with the new LSP */
242 for (i = 0; i < LPCO8; i++)
243 old_lsp[i] = lsp[i];
244 }
245}
246
247static Float FNevChebP(Float x, /* (i) value */
248 Float *c, /* (i) coefficient array */
249 int nd2)
250{
251 Float t;
252 Float b[NAB((8 >> 1) + 1)];
253 int i;
254
255 t = x*2;
256 b[0] = c[nd2];
257 b[1] = c[nd2 - 1] + t*b[0];
258 for (i = 2; i < nd2; i++)
259 b[i] = c[nd2 - i] - b[i - 2] + t * b[i - 1];
260 return (c[0] - b[nd2 - 2] + x * b[nd2 - 1]);
261}