Line data Source code
1 : /******************************************************************************************************
2 :
3 : (C) 2022-2025 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB,
4 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
5 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
6 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
7 : contributors to this repository. All Rights Reserved.
8 :
9 : This software is protected by copyright law and by international treaties.
10 : The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB,
11 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
12 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
13 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
14 : contributors to this repository retain full ownership rights in their respective contributions in
15 : the software. This notice grants no license of any kind, including but not limited to patent
16 : license, nor is any license granted by implication, estoppel or otherwise.
17 :
18 : Contributors are required to enter into the IVAS codec Public Collaboration agreement before making
19 : contributions.
20 :
21 : This software is provided "AS IS", without any express or implied warranties. The software is in the
22 : development stage. It is intended exclusively for experts who have experience with such software and
23 : solely for the purpose of inspection. All implied warranties of non-infringement, merchantability
24 : and fitness for a particular purpose are hereby disclaimed and excluded.
25 :
26 : Any dispute, controversy or claim arising under or in relation to providing this software shall be
27 : submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in
28 : accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and
29 : the United Nations Convention on Contracts on the International Sales of Goods.
30 :
31 : *******************************************************************************************************/
32 :
33 : /*====================================================================================
34 : EVS Codec 3GPP TS26.443 Nov 04, 2021. Version 12.14.0 / 13.10.0 / 14.6.0 / 15.4.0 / 16.3.0
35 : ====================================================================================*/
36 :
37 : #include <stdint.h>
38 : #include "options.h"
39 : #ifdef DEBUGGING
40 : #include "debug.h"
41 : #endif
42 : #include "cnst.h"
43 : #include "rom_com.h"
44 : #include "prot.h"
45 : #include "wmc_auto.h"
46 :
47 : /*-------------------------------------------------------------------*
48 : * Local constants
49 : *-------------------------------------------------------------------*/
50 :
51 : /* The conversion modes. */
52 : #define DOWNCONV 0
53 : #define UPCONV 1
54 : /* The cap of the inverse power spectrum. */
55 : #define MAXPOWERSPECT 1e-5f
56 : #define N50 GRID50_POINTS
57 : #define N40 GRID40_POINTS
58 :
59 :
60 : /*-------------------------------------------------------------------*
61 : * Local function prototypes
62 : *-------------------------------------------------------------------*/
63 :
64 : static void powerspect( const float x[], const int16_t N, float R[], float S[], float G[], const int16_t mode );
65 :
66 : static void spectautocorr( const float x[], const int16_t N, const float G[], float r[] );
67 :
68 : static void zeros2poly( float x[], float R[], float S[] );
69 :
70 : static void polydecomp( float A[], float P[], float Q[] );
71 :
72 : static void cheb2poly( float P[] );
73 :
74 :
75 : /*---------------------------------------------------------------------*
76 : * lsp_convert_poly()
77 : *
78 : * Converts the LP filter estimated at 16.0 kHz sampling rate down
79 : * 12.8 kHz frequency scale or alternatively from 12.8 kHz up to
80 : * 16.0 kHz. The former is called down conversation and latter up
81 : * conversion. The resulting LP filter is characterized with its
82 : * line spectrum pairs. The original Lp filter can be either in
83 : * its immittance, used for the AMR-WB IO mode, or line spectrum
84 : * pair representation.
85 : *
86 : * The conversion is based the autocorrelation computed from the
87 : * power spectrum of the LP filter that is truncated or extrapolated
88 : * to the desired frequency scale.
89 : *---------------------------------------------------------------------*/
90 :
91 110024 : int16_t lsp_convert_poly(
92 : float w[], /* i/o: LSP or ISP parameters */
93 : const int16_t L_frame, /* i : flag for up or down conversion */
94 : const int16_t Opt_AMRWB /* i : flag for the AMR-WB IO mode */
95 : )
96 : {
97 : float epsP[M + 1];
98 : float G[GRID50_POINTS], r[M + 1], A[M + 1], A_old[M + 1], R[NC + 1], S[NC + 1];
99 : int16_t i;
100 : int16_t flag;
101 :
102 : /*---------------------------------------------------------------------*
103 : * Because AMR-WB IO uses immittance spectrum frequency representation
104 : * instead of line spectrum frequency representation, the input
105 : * parameters do not give the zeros of the polynomials R(x) and S(x).
106 : * Hence R(x) and S(x) are formed via the polynomial A(z) of the linear
107 : * prediction filter.
108 : *---------------------------------------------------------------------*/
109 :
110 110024 : if ( Opt_AMRWB )
111 : {
112 0 : isp2a( w, A_old, M );
113 0 : polydecomp( A_old, R, S );
114 : }
115 :
116 : /*---------------------------------------------------------------------*
117 : * Form the polynomials R(x) and S(x) from their zeros that are the
118 : * line spectrum pairs of A(z). The polynomial coefficients can be
119 : * scaled for convenience, because scaling will not affect the
120 : * resulting LP coefficients. Scaling by 128 gives the correct offset
121 : * to the power spectrum for n = 16.
122 : *---------------------------------------------------------------------*/
123 :
124 : else
125 : {
126 110024 : zeros2poly( w, R, S );
127 1100240 : for ( i = 0; i <= NC; i++ )
128 : {
129 990216 : R[i] *= 128.0f;
130 990216 : S[i] *= 128.0f;
131 : }
132 110024 : lsp2a_stab( w, A_old, M );
133 : }
134 :
135 : /*---------------------------------------------------------------------*
136 : * Conversion from 16.0 kHz down to 12.8 kHz. The power spectrum
137 : * needs to be computed only up to 6.4 kHz, because the upper band
138 : * is omitted.
139 : *---------------------------------------------------------------------*/
140 :
141 110024 : if ( L_frame == L_FRAME )
142 : {
143 10234 : powerspect( grid50, N50, R, S, G, DOWNCONV );
144 10234 : spectautocorr( grid40, N40, G, r );
145 : }
146 :
147 : /*---------------------------------------------------------------------*
148 : * Conversion from 12.8 kHz up to 16.0 kHz.
149 : * Compute the power spectrum of the LP filter, extrapolate the
150 : * power spectrum from 6.4 kHz to 8.0 kHz, and compute auto-
151 : * correlation on this power spectrum.
152 : *---------------------------------------------------------------------*/
153 :
154 : else
155 : {
156 99790 : powerspect( grid40, N40, R, S, G, UPCONV );
157 1097690 : for ( i = N40; i < N50; i++ )
158 : {
159 997900 : G[i] = G[N40 - 1];
160 : }
161 :
162 99790 : spectautocorr( grid50, N50, G, r );
163 : }
164 :
165 : /*---------------------------------------------------------------------*
166 : * Compute the linear prediction coefficients from the autocorrelation
167 : * and convert to line spectrum pairs.
168 : *---------------------------------------------------------------------*/
169 :
170 110024 : flag = lev_dur( A, r, M, epsP );
171 110024 : a2lsp_stab( A, w, stable_LSP );
172 :
173 110024 : return ( flag );
174 : }
175 :
176 :
177 : /*---------------------------------------------------------------------*
178 : * powerspect()
179 : *
180 : * Computes the power spectrum G(w) = 1/|A(w)|^2 at N points on
181 : * the real axis x = cos w by utilizing the line spectrum frequency
182 : * decomposition
183 : *
184 : * A(z) = (P(z) + Q(z))/2,
185 : *
186 : * where assuming A(z) of an even degree n,
187 : *
188 : * P(z) = [A(z) + z^(n+1) A(1/z)]/(1/z + 1),
189 : * Q(z) = [A(z) - z^(n+1) A(1/z)]/(1/z - 1).
190 : *
191 : * The zeros of these polynomials give the line spectrum frequencies
192 : * of A(z). It can be shown that for an even n,
193 : *
194 : * |A(x)|^2 = 2 (1 + x) R(x)^2 + 2 (1 - x) S(x)^2,
195 : *
196 : * where x = cos w, and R(x) and S(x) are the direct polynomials
197 : * resulting from the Chebyshev series representation of P(z)
198 : * and Q(z).
199 : *
200 : * This routine assumes the grid X = 1, x[0], x[1], .., x[m-1],
201 : * -, ..., -x[1], -x[0], -1 such that x[i] = cos((i+1)*pi/N) for
202 : * evaluating the power spectrum. Only m = (N-1)/2 - 1 grid points
203 : * need to be stored, because cos(0) and cos(pi/2) are trivial,
204 : * and the points above pi/2 are obtained readily using the symmetry
205 : * of cosine.
206 : *
207 : * The power spectrum can be scaled as a*G[], where a is chosen
208 : * for convenience. This is because the scaling has no impact on
209 : * the LP coefficients to be determined based on the power spectrum.
210 : *---------------------------------------------------------------------*/
211 :
212 110024 : static void powerspect(
213 : const float x[], /* i : Grid points x[0:N-1] */
214 : const int16_t N, /* i : Number of grid points */
215 : float R[], /* i : Coefficients of R(x) in R[0:NC] */
216 : float S[], /* i : Coefficients of S(x) in S[0:NC] */
217 : float G[], /* o : Power spectrum G[0:N] */
218 : const int16_t mode /* i : Flag for up or down conversion */
219 : )
220 : {
221 : float re, ro, se, so;
222 : float s0, s1, r0, r1, x0, x1, x2;
223 : Word16 i, j;
224 : Word16 iuni, imid;
225 :
226 : /* init buffer */
227 4723348 : for ( i = 0; i < N; i++ )
228 : {
229 4613324 : G[i] = 0;
230 : }
231 : /*---------------------------------------------------------------------*
232 : * Down conversion yields iuni unique grid points that do not have
233 : * symmetric counterparts above x = cos(pi/2) = 0.
234 : * Set the mid point of the frequency grid.
235 : *---------------------------------------------------------------------*/
236 :
237 110024 : if ( mode == DOWNCONV )
238 : {
239 10234 : iuni = ( GRID50_POINTS - 1 ) / 5 - 1;
240 10234 : imid = ( GRID50_POINTS - 1 ) / 2;
241 : }
242 :
243 : /*---------------------------------------------------------------------*
244 : * Power spectrum at x = cos(pi) = -1 that is not needed in down
245 : * conversion. Set the mid point of the frequency grid.
246 : *---------------------------------------------------------------------*/
247 :
248 : else
249 : {
250 99790 : iuni = 0;
251 99790 : imid = ( GRID40_POINTS - 1 ) / 2;
252 :
253 99790 : s0 = S[0];
254 898110 : for ( j = 1; j <= NC; j++ )
255 : {
256 798320 : s0 = S[j] - s0;
257 : }
258 99790 : G[N - 1] = 4.0f * s0 * s0;
259 : }
260 :
261 : /*---------------------------------------------------------------------*
262 : * Power spectrum at x = cos(0) = 1.
263 : *---------------------------------------------------------------------*/
264 :
265 110024 : r0 = R[0];
266 990216 : for ( j = 1; j <= NC; j++ )
267 : {
268 880192 : r0 += R[j];
269 : }
270 110024 : r0 = max( r0, 0.000000953674316f );
271 110024 : G[0] = 4.0f * r0 * r0;
272 :
273 : /*---------------------------------------------------------------------*
274 : * Power spectrum at x = cos(pi/2) = 0.
275 : *---------------------------------------------------------------------*/
276 :
277 110024 : r0 = R[NC];
278 110024 : s0 = S[NC];
279 110024 : r0 = r0 * r0;
280 110024 : s0 = s0 * s0;
281 :
282 110024 : G[imid] = 2.0f * ( r0 + s0 );
283 :
284 : /*---------------------------------------------------------------------*
285 : * Power spectrum at unique points that do not have symmetric
286 : * counterparts at x > cos(pi/2) = 0.
287 : *---------------------------------------------------------------------*/
288 :
289 202130 : for ( i = 1; i <= iuni; i++ )
290 : {
291 92106 : x0 = x[i - 1];
292 92106 : r0 = R[0];
293 92106 : s0 = S[0];
294 :
295 828954 : for ( j = 1; j <= NC; j++ )
296 : {
297 736848 : r0 = x0 * r0 + R[j];
298 736848 : s0 = x0 * s0 + S[j];
299 : }
300 92106 : r0 = ( 1.0f + x0 ) * r0 * r0;
301 92106 : s0 = ( 1.0f - x0 ) * s0 * s0;
302 :
303 92106 : G[i] = 2.0f * ( r0 + s0 );
304 : }
305 :
306 : /*---------------------------------------------------------------------*
307 : * Power spectrum at points other than x = -1, 0, and 1 and unique
308 : * points is computed using the anti-symmetry of the grid relative
309 : * to the midpoint x = 0 in order to reduce looping.
310 : *---------------------------------------------------------------------*/
311 :
312 2159544 : for ( ; i < imid; i++ )
313 : {
314 2049520 : x0 = x[i - 1];
315 2049520 : x2 = x0 * x0;
316 :
317 2049520 : re = R[0];
318 2049520 : se = S[0];
319 2049520 : ro = R[1];
320 2049520 : so = S[1];
321 :
322 8198080 : for ( j = 2; j < NC; j += 2 )
323 : {
324 6148560 : re = x2 * re + R[j];
325 6148560 : ro = x2 * ro + R[j + 1];
326 6148560 : se = x2 * se + S[j];
327 6148560 : so = x2 * so + S[j + 1];
328 : }
329 :
330 2049520 : re = x2 * re + R[j];
331 2049520 : ro *= x0;
332 2049520 : se = x2 * se + S[j];
333 2049520 : so *= x0;
334 :
335 2049520 : r0 = re + ro;
336 2049520 : s0 = se + so;
337 2049520 : r1 = re - ro;
338 2049520 : s1 = se - so;
339 :
340 2049520 : x1 = 1.0f + x0;
341 2049520 : x2 = 1.0f - x0;
342 :
343 2049520 : r0 = x1 * r0 * r0;
344 2049520 : s0 = x2 * s0 * s0;
345 2049520 : r1 = x2 * r1 * r1;
346 2049520 : s1 = x1 * s1 * s1;
347 :
348 2049520 : G[i] = 2.0f * ( r0 + s0 );
349 2049520 : G[N - i - 1] = 2.0f * ( r1 + s1 );
350 : }
351 :
352 : /*---------------------------------------------------------------------*
353 : * Compute the power spectrum 1/|A(x)|^2 from |A(x)|^2 with logic
354 : * to prevent division by small number and upper bound the spectrum.
355 : * This upper bound is implicit in fixed point arithmetic, but is
356 : * needed explicitly in floating point implementations to avoid numeric
357 : * problems.
358 : *---------------------------------------------------------------------*/
359 :
360 4723348 : for ( i = 0; i < N; i++ )
361 : {
362 4613324 : if ( G[i] < MAXPOWERSPECT )
363 : {
364 102345 : G[i] = MAXPOWERSPECT;
365 : }
366 4613324 : G[i] = 1.0f / G[i];
367 : }
368 110024 : return;
369 : }
370 :
371 : /*---------------------------------------------------------------------*
372 : * spectautocorr()
373 : *
374 : * Computes the autocorrelation r[j] for j = 0, 1, ..., M from
375 : * the power spectrum P(w) by using the rectangle rule to
376 : * approximate the integral
377 : *
378 : * 1 pi
379 : * r[j] = --- I P(w) cos(j*w) dw.
380 : * 2*pi -pi
381 : *
382 : * It is sufficient to evaluate the integrand only from w = 0 to
383 : * w = pi due to the symmetry P(-w) = P(w). We can further
384 : * employ the relation
385 : *
386 : * cos(j*(pi - w)) = (-1)^j cos(j*w)
387 : *
388 : * to use symmetries relative to w = pi/2 for w in (0, pi/2).
389 : *
390 : * When applying the rectangle rule, it is convenient to evaluate
391 : * w = 0, w = pi/2, and w = pi separately. By using a frequency
392 : * grid of N points, we can express the rectangle rule as
393 : *
394 : * r[j] = G[0] + 2*a*G[(N-1)/2] + b*G[N-1]
395 : *
396 : * M
397 : * + 2 sum (G[i] - G[N-i-1]) cos(j*x[i])
398 : * i=1
399 : *
400 : * where G[i] is the power spectrum at x[i] = i*pi/(N-1) and
401 : * M = (N-1)/2 - 1 is the number of the grid points in the
402 : * interval(0, pi/2). The number of grid points N is assumed odd.
403 : *
404 : * The coefficients
405 : *
406 : * b = (-1)^j
407 : * a = (1 + (-1)^(j+1))(-1)^floor(j/2)
408 : *
409 : * follow from the properties of the cosine function. The
410 : * computation further uses the recursion
411 : *
412 : * cos(j*w) = 2*cos(w)*cos((j-1)*w) - cos((j-2)*w)
413 : *
414 : * Note that the autocorrelation can be scaled for convenience,
415 : * because this scaling has no impact on the LP coefficients to be
416 : * calculated from the autocorrelation. The expression of r[j] thus
417 : * omits the division by N.
418 : *
419 : * See the powerspect function on the definition of the grid.
420 : *
421 : * References
422 : * J. Makhoul, "Spectral linear prediction: properties and
423 : * applications," IEEE Trans. on Acoustics, Speech and Signal
424 : * Processing, Vol. 23, No. 3, pp.283-296, June 1975
425 : *---------------------------------------------------------------------*/
426 :
427 110024 : static void spectautocorr(
428 : const float x[], /* i : Grid points x[0:N-1] */
429 : const int16_t N, /* i : Number of grid points */
430 : const float G[], /* i : Power spectrum G[0:N-1] */
431 : float r[] /* o : Autocorrelation r[0:M] */
432 : )
433 : {
434 : float c[M + 1]; /* c[j] = cos(j*w) */
435 : float gp, gn, c2;
436 : Word16 i, j;
437 : Word16 imid;
438 :
439 : /*---------------------------------------------------------------------*
440 : * The midpoint of the cosine table x of N entries. Only the entries
441 : * x[0] = cos(d), x[1] = cos(2*d), ..., x[imid-1] = cos((imid-1)*d)
442 : * need to be stored due to trivial cos(0), cos(pi/2), cos(pi), and
443 : * symmetry relative to pi/2. Here d = pi/(N - 1).
444 : *---------------------------------------------------------------------*/
445 :
446 110024 : imid = ( N - 1 ) / 2;
447 :
448 : /*---------------------------------------------------------------------*
449 : * Autocorrelation r[j] at zero lag j = 0 for the upper half of the
450 : * unit circle, but excluding the points x = cos(0) and x = cos(pi).
451 : *---------------------------------------------------------------------*/
452 :
453 110024 : r[0] = G[1];
454 5288836 : for ( i = 2; i < N - 1; i++ )
455 : {
456 5178812 : r[0] += G[i];
457 : }
458 :
459 : /*---------------------------------------------------------------------*
460 : * Initialize the autocorrelation r[j] at lags greater than zero
461 : * by adding the midpoint x = cos(pi/2) = 0.
462 : *---------------------------------------------------------------------*/
463 :
464 110024 : r[1] = 0.0f;
465 110024 : r[2] = -G[imid];
466 :
467 880192 : for ( j = 3; j < M; j += 2 )
468 : {
469 770168 : r[j] = 0.0f;
470 770168 : r[j + 1] = -r[j - 1];
471 : }
472 :
473 : /*---------------------------------------------------------------------*
474 : * Autocorrelation r[j] at lags j = 1, 2, ..., M. The computation
475 : * employes the relation cos(j*(pi - w)) = (-1)^j cos(j*w) and
476 : * cos(j*w) = 2*cos(w)*cos((j-1)*w) - cos((j-2)*w) for obtaining
477 : * the cosine c[j] = cos(j*w).
478 : *---------------------------------------------------------------------*/
479 :
480 110024 : c[0] = 1.0f;
481 :
482 2699430 : for ( i = 1; i < imid; i++ )
483 : {
484 2589406 : gp = G[i] + G[N - i - 1];
485 2589406 : gn = G[i] - G[N - i - 1];
486 :
487 2589406 : c[1] = x[i - 1];
488 2589406 : c2 = 2.0f * c[1];
489 :
490 2589406 : r[1] += gn * c[1];
491 :
492 20715248 : for ( j = 2; j < M; j += 2 )
493 : {
494 18125842 : c[j] = c2 * c[j - 1] - c[j - 2];
495 18125842 : r[j] += gp * c[j];
496 18125842 : c[j + 1] = c2 * c[j] - c[j - 1];
497 18125842 : r[j + 1] += gn * c[j + 1];
498 : }
499 :
500 2589406 : c[j] = c2 * c[j - 1] - c[j - 2];
501 2589406 : r[j] += gp * c[j];
502 : }
503 :
504 : /*---------------------------------------------------------------------*
505 : * Add the endpoints x = cos(0) = 1 and x = cos(pi) = -1 as
506 : * well as the lower half of the unit circle.
507 : *---------------------------------------------------------------------*/
508 :
509 110024 : gp = G[0] + G[N - 1];
510 110024 : gn = G[0] - G[N - 1];
511 :
512 990216 : for ( j = 0; j < M; j += 2 )
513 : {
514 880192 : r[j] = 2.0f * r[j] + gp;
515 880192 : r[j + 1] = 2.0f * r[j + 1] + gn;
516 : }
517 110024 : r[j] = 2.0f * r[j] + gp;
518 :
519 110024 : return;
520 : }
521 :
522 : /*---------------------------------------------------------------------*
523 : * zeros2poly()
524 : *
525 : * Computes the coefficients of the polynomials
526 : *
527 : * R(x) = prod (x - x[i]),
528 : * i = 0,2,4,...
529 : *
530 : * S(x) = prod (x - x[i]),
531 : * i = 1,3,5,...
532 : *
533 : * given their zeros x[i] for i = 0, 1, ..., n-1. The routine
534 : * assumes n = 1 or even n greater than or equal to 4.
535 : *
536 : * The polynomial coefficients are returned in R[0:n/2-1] and
537 : * S[0:n/2-1]. The leading coefficients are in R[0] and S[0].
538 : *
539 : * In this implementation, n is set to a compile time constant.
540 : *---------------------------------------------------------------------*/
541 :
542 110024 : static void zeros2poly(
543 : float x[], /* i : Zeros of R(x) and S(x) */
544 : float R[], /* o : Coefficients of R(x) */
545 : float S[] /* o : Coefficients of S(x) */
546 : )
547 : {
548 : float xr, xs;
549 : Word16 i, j;
550 :
551 110024 : R[0] = 1.0f;
552 110024 : R[1] = -*x++;
553 110024 : S[0] = 1.0f;
554 110024 : S[1] = -*x++;
555 :
556 880192 : for ( i = 2; i <= NC; i++ )
557 : {
558 770168 : xr = -*x++;
559 770168 : xs = -*x++;
560 :
561 770168 : R[i] = xr * R[i - 1];
562 770168 : S[i] = xs * S[i - 1];
563 :
564 3850840 : for ( j = i - 1; j > 0; j-- )
565 : {
566 3080672 : R[j] += xr * R[j - 1];
567 3080672 : S[j] += xs * S[j - 1];
568 : }
569 : }
570 :
571 110024 : return;
572 : }
573 :
574 : /*---------------------------------------------------------------------*
575 : * polydecomp()
576 : *
577 : * Computes the coefficients of the symmetric and antisymmetric
578 : * polynomials P(z) and Q(z) that define the line spectrum pair
579 : * decomposition of a given polynomial A(z) of order n. For even n,
580 : *
581 : * P(z) = [A(z) + z^(n+1) A(1/z)]/(1/z + 1),
582 : * Q(z) = [A(z) - z^(n+1) A(1/z)]/(1/z - 1),
583 : *
584 : * These polynomials are then expressed in their direct form,
585 : * respectively, R(x) and S(x), on the real axis x = cos w using
586 : * explicit Chebyshev polynomials of the first kind.
587 : *
588 : * The coefficients of the polynomials R(x) and S(x) are returned
589 : * in R[0:n/2] and S[0:n/2] for the given linear prediction
590 : * coefficients A[0:n/2]. Note that R(x) and S(x) are formed in
591 : * place such that P(z) is stored in the same array than R(x),
592 : * and Q(z) is stored in the same array than S(x).
593 : *
594 : * The routines assumes n = 16.
595 : *---------------------------------------------------------------------*/
596 :
597 0 : static void polydecomp(
598 : float A[], /* i : linear prediction coefficients */
599 : float R[], /* o : coefficients of R(x) */
600 : float S[] /* o : coefficients of S(x) */
601 : )
602 : {
603 0 : float *P = &R[0], *Q = &S[0];
604 : Word16 i, j;
605 :
606 0 : P[0] = A[0];
607 0 : Q[0] = A[0];
608 0 : for ( i = 1, j = M; i <= NC; i++, j-- )
609 : {
610 0 : P[i] = A[i] + A[j] - P[i - 1];
611 0 : Q[i] = A[i] - A[j] + Q[i - 1];
612 : }
613 :
614 0 : cheb2poly( P );
615 0 : cheb2poly( Q );
616 :
617 0 : return;
618 : }
619 :
620 : /*---------------------------------------------------------------------*
621 : * cheb2poly()
622 : *
623 : * Computes the coefficients of the explicit Chebyshev polynomial
624 : * P(x) = P[0]*x^n + P[1]*x^(n-1) + ... + P[n] given the coefficients
625 : * of the series
626 : *
627 : * C(x) = C[0]*T_n(x) + C[1]*T_n-1(x) + ... + C[n]*T_0(x),
628 : *
629 : * where T_n(x) is the nth Chebyshev polynomial of the first kind.
630 : * This implementation assumes C[0] = 1. Only the value n = 8 is
631 : * supported.
632 : *
633 : * The conversion from C(x) to P(x) is done in place such that the
634 : * coefficients of C(x) are given in P[0:8] and those of P(x) are
635 : * returned in the same array.
636 : *---------------------------------------------------------------------*/
637 :
638 0 : static void cheb2poly(
639 : float P[] /* i/o: The coefficients of C(x) and P(x) */
640 : )
641 : {
642 : float c1, c2, c3, c4, c5, c6, c7, c8;
643 :
644 0 : c1 = P[1];
645 0 : c2 = P[2];
646 0 : c3 = P[3];
647 0 : c4 = P[4];
648 0 : c5 = P[5];
649 0 : c6 = P[6];
650 0 : c7 = P[7];
651 0 : c8 = P[8];
652 :
653 0 : P[0] = 128.0f;
654 0 : P[1] = 64.0f * c1;
655 0 : P[2] = 32.0f * c2 - 256.0f;
656 0 : P[3] = 16.0f * c3 - 112.0f * c1;
657 0 : P[4] = 160.0f - 48.0f * c2 + 8.0f * c4;
658 0 : P[5] = 56.0f * c1 - 20.0f * c3 + 4.0f * c5;
659 0 : P[6] = 18.0f * c2 - 32.0f - 8.0f * c4 + 2.0f * c6;
660 0 : P[7] = 5.0f * c3 - 7.0f * c1 - 3.0f * c5 + c7;
661 0 : P[8] = 1.0f - c2 + c4 - c6 + 0.5f * c8;
662 :
663 0 : return;
664 : }
|
