Line data Source code
1 : /*====================================================================================
2 : EVS Codec 3GPP TS26.452 Aug 12, 2021. Version 16.3.0
3 : ====================================================================================*/
4 :
5 : #include <stdint.h>
6 : #include "options.h"
7 : #include "cnst.h"
8 : #include "rom_com.h"
9 : #include "prot_fx.h"
10 : #include "wmc_auto.h"
11 :
12 :
13 : /*-------------------------------------------------------------------*
14 : * Local constants
15 : *-------------------------------------------------------------------*/
16 :
17 : /* The conversion modes. */
18 : #define DOWNCONV 0
19 : #define UPCONV 1
20 : /* The cap of the inverse power spectrum. */
21 : #define MAXPOWERSPECT 1e-5f
22 : #define N50 GRID50_POINTS
23 : #define N40 GRID40_POINTS
24 :
25 : /*-------------------------------------------------------------------*
26 : * Local function prototypes
27 : *-------------------------------------------------------------------*/
28 : static void powerspect_fx( const Word16 x[], Word16 N, Word32 R[], Word32 S[], Word32 G[], Word16 mode );
29 : static void spectautocorr_fx( const Word16 x[], const Word16 N, const Word32 G[], Word16 rh[], Word16 rl[] );
30 : static Word32 b_inv_sq( const Word32 in32, const Word16 exp_in );
31 : static Word32 inv_pow( const Word32 re, const Word32 se, const Word16 x );
32 : static void zeros2poly_fx( Word16 x[], Word32 R[], Word32 S[] );
33 : static void polydecomp_fx( Word16 A[], Word32 P[], Word32 Q[] );
34 : static void cheb2poly_fx( Word32 L_P[] );
35 :
36 : /*---------------------------------------------------------------------*
37 : * lsp_convert_poly()
38 : *
39 : * Converts the LP filter estimated at 16.0 kHz sampling rate down
40 : * 12.8 kHz frequency scale or alternatively from 12.8 kHz up to
41 : * 16.0 kHz. The former is called down conversation and latter up
42 : * conversion. The resulting LP filter is characterized with its
43 : * line spectrum pairs. The original Lp filter can be either in
44 : * its immittance, used for the AMR-WB IO mode, or line spectrum
45 : * pair representation.
46 : *
47 : * The conversion is based the autocorrelation computed from the
48 : * power spectrum of the LP filter that is truncated or extrapolated
49 : * to the desired frequency scale.
50 : *---------------------------------------------------------------------*/
51 :
52 102567 : Word16 lsp_convert_poly_fx(
53 : Word16 w[], /* i/o: LSP or ISP parameters */
54 : const Word16 L_frame, /* i : flag for up or down conversion */
55 : const Word16 Opt_AMRWB /* i : flag for the AMR-WB IO mode */
56 : )
57 : {
58 : Word16 flag;
59 :
60 : Word32 G[GRID50_POINTS];
61 : Word16 i;
62 : Word16 A[M + 1];
63 : Word32 R[NC + 1], S[NC + 1];
64 : Word32 epsP[M + 1];
65 : Word16 rh[M + 1], rl[M + 1];
66 : Word16 oldA[M + 3];
67 :
68 : /*---------------------------------------------------------------------*
69 : * Because AMR-WB IO mode uses immittance spectrum frequency representation
70 : * instead of line spectrum frequency representation, the input
71 : * parameters do not give the zeros of the polynomials R(x) and S(x).
72 : * Hence R(x) and S(x) are formed via the polynomial A(z) of the linear
73 : * prediction filter.
74 : *---------------------------------------------------------------------*/
75 :
76 102567 : IF( Opt_AMRWB )
77 : {
78 0 : E_LPC_f_isp_a_conversion( w, oldA, M );
79 :
80 0 : polydecomp_fx( oldA, R, S );
81 : }
82 :
83 : /*---------------------------------------------------------------------*
84 : * Form the polynomials R(x) and S(x) from their zeros that are the
85 : * line spectrum pairs of A(z). The polynomial coefficients can be
86 : * scaled for convenience, because scaling will not affect the
87 : * resulting LP coefficients. Scaling by 128 gives the correct offset
88 : * to the power spectrum for n = 16.
89 : *---------------------------------------------------------------------*/
90 :
91 : ELSE
92 : {
93 102567 : E_LPC_f_lsp_a_conversion( w, oldA, M );
94 :
95 102567 : zeros2poly_fx( w, R, S );
96 : }
97 :
98 : /*---------------------------------------------------------------------*
99 : * Conversion from 16.0 kHz down to 12.8 kHz. The power spectrum
100 : * needs to be computed only up to 6.4 kHz, because the upper band
101 : * is omitted.
102 : *---------------------------------------------------------------------*/
103 :
104 102567 : IF( EQ_16( L_frame, L_FRAME ) )
105 : {
106 5140 : powerspect_fx( grid50_fx, N50, R, S, G, DOWNCONV );
107 5140 : spectautocorr_fx( grid40_fx, N40, G, rh, rl );
108 : }
109 : /*---------------------------------------------------------------------*
110 : * Conversion from 12.8 kHz up to 16.0 kHz.
111 : * Compute the power spectrum of the LP filter, extrapolate the
112 : * power spectrum from 6.4 kHz to 8.0 kHz, and compute auto-
113 : * correlation on this power spectrum.
114 : *---------------------------------------------------------------------*/
115 :
116 : ELSE
117 : {
118 97427 : powerspect_fx( grid40_fx, N40, R, S, G, UPCONV );
119 :
120 1071697 : FOR( i = N40; i < N50; i++ )
121 : {
122 974270 : G[i] = G[N40 - 1];
123 974270 : move32();
124 : }
125 :
126 97427 : spectautocorr_fx( grid50_fx, N50, G, rh, rl );
127 : }
128 :
129 :
130 : /*---------------------------------------------------------------------*
131 : * Compute the linear prediction coefficients from the autocorrelation
132 : * and convert to line spectrum pairs.
133 : *---------------------------------------------------------------------*/
134 102567 : flag = E_LPC_lev_dur( rh, rl, A, epsP, M, oldA );
135 102567 : E_LPC_a_lsp_conversion( A, w, stable_LSP_fx, M );
136 :
137 :
138 102567 : return ( flag );
139 : }
140 :
141 : /*---------------------------------------------------------------------*
142 : * powerspect()
143 : *
144 : * Computes the power spectrum G(w) = 1/|A(w)|^2 at N points on
145 : * the real axis x = cos w by utilizing the line spectrum frequency
146 : * decomposition
147 : *
148 : * A(z) = (P(z) + Q(z))/2,
149 : *
150 : * where assuming A(z) of an even degree n,
151 : *
152 : * P(z) = [A(z) + z^(n+1) A(1/z)]/(1/z + 1),
153 : * Q(z) = [A(z) - z^(n+1) A(1/z)]/(1/z - 1).
154 : *
155 : * The zeros of these polynomials give the line spectrum frequencies
156 : * of A(z). It can be shown that for an even n,
157 : *
158 : * |A(x)|^2 = 2 (1 + x) R(x)^2 + 2 (1 - x) S(x)^2,
159 : *
160 : * where x = cos w, and R(x) and S(x) are the direct polynomials
161 : * resulting from the Chebyshev series representation of P(z)
162 : * and Q(z).
163 : *
164 : * This routine assumes the grid X = 1, x[0], x[1], .., x[m-1],
165 : * -, ..., -x[1], -x[0], -1 such that x[i] = cos((i+1)*pi/N) for
166 : * evaluating the power spectrum. Only m = (N-1)/2 - 1 grid points
167 : * need to be stored, because cos(0) and cos(pi/2) are trivial,
168 : * and the points above pi/2 are obtained readily using the symmetry
169 : * of cosine.
170 : *
171 : * The power spectrum can be scaled as a*G[], where a is chosen
172 : * for convenience. This is because the scaling has no impact on
173 : * the LP coefficients to be determined based on the power spectrum.
174 : *---------------------------------------------------------------------*/
175 :
176 102567 : void powerspect_fx(
177 : const Word16 x[], /* i: Q15 Grid points x[0:m-1] */
178 : Word16 N, /* i: Number of grid points */
179 : Word32 R[], /* i: Q20 Coefficients of R(x) in R[0:NC] */
180 : Word32 S[], /* i: Q20 Coefficients of S(x) in S[0:NC] */
181 : Word32 G[], /* o: Q15 Power spectrum G[0:N] */
182 : Word16 mode /* i: Flag for up or down conversion */
183 : )
184 : {
185 : Word32 s0, se, so, r0, re, ro;
186 : Word16 i, j;
187 : Word16 iuni, imid;
188 : Word32 L_tmp;
189 : Word16 x2;
190 : Word32 mh;
191 : UWord16 ml;
192 :
193 : /*---------------------------------------------------------------------*
194 : * Down conversion yields iuni unique grid points that do not have
195 : * symmetric counterparts above x = cos(pi/2) = 0.
196 : * Set the mid point of the frequency grid.
197 : *---------------------------------------------------------------------*/
198 :
199 102567 : IF( mode == DOWNCONV )
200 : {
201 5140 : iuni = ( GRID50_POINTS - 1 ) / 5 - 1;
202 5140 : move16();
203 5140 : imid = ( GRID50_POINTS - 1 ) / 2;
204 5140 : move16();
205 : }
206 :
207 : /*---------------------------------------------------------------------*
208 : * Power spectrum x = cos(pi) = -1 that is not needed in down
209 : * conversion. Set the mid point of the frequency grid.
210 : *---------------------------------------------------------------------*/
211 :
212 : ELSE
213 : {
214 97427 : iuni = 0;
215 97427 : move16();
216 97427 : imid = ( GRID40_POINTS - 1 ) / 2;
217 97427 : move16();
218 :
219 97427 : G[N - 1] = S[0];
220 97427 : move32();
221 :
222 876843 : FOR( j = 1; j <= NC; j++ )
223 : {
224 779416 : G[N - 1] = L_sub( S[j], G[N - 1] );
225 779416 : move32();
226 : }
227 97427 : G[N - 1] = b_inv_sq( G[N - 1], 19 );
228 97427 : move32();
229 : }
230 :
231 : /*---------------------------------------------------------------------*
232 : * Power spectrum x = cos(0) = 1.
233 : *---------------------------------------------------------------------*/
234 :
235 102567 : G[0] = R[0];
236 102567 : move32();
237 923103 : FOR( j = 1; j <= NC; j++ )
238 : {
239 820536 : G[0] = L_add( R[j], G[0] );
240 820536 : move32();
241 : }
242 :
243 102567 : G[0] = b_inv_sq( L_max( G[0], 1 ), 19 );
244 102567 : move32();
245 :
246 : /*---------------------------------------------------------------------*
247 : * Power spectrum at x = cos(pi/2) = 0.
248 : *---------------------------------------------------------------------*/
249 102567 : G[imid] = inv_pow( R[NC], S[NC], 0 );
250 102567 : move32();
251 :
252 : /*---------------------------------------------------------------------*
253 : * Power spectrum at unique points that do not have symmetric
254 : * counterparts at x > cos(pi/2) = 0.
255 : *---------------------------------------------------------------------*/
256 :
257 148827 : FOR( i = 1; i <= iuni; i++ )
258 : {
259 46260 : Mpy_32_16_ss( R[0], x[i - 1], &mh, &ml );
260 46260 : r0 = L_add( R[1], mh );
261 46260 : Mpy_32_16_ss( S[0], x[i - 1], &mh, &ml );
262 46260 : s0 = L_add( S[1], mh );
263 :
264 :
265 370080 : FOR( j = 2; j <= NC; j++ )
266 : {
267 323820 : Mpy_32_16_ss( r0, x[i - 1], &mh, &ml );
268 323820 : r0 = L_add( R[j], mh );
269 :
270 323820 : Mpy_32_16_ss( s0, x[i - 1], &mh, &ml );
271 323820 : s0 = L_add( S[j], mh );
272 : }
273 :
274 46260 : G[i] = inv_pow( r0, s0, x[i - 1] );
275 46260 : move32();
276 : }
277 :
278 : /*---------------------------------------------------------------------*
279 : * Power spectrum at points other than x = -1, 0, and 1 and unique
280 : * points is computed using the anti-symmetry of the grid relative
281 : * to the midpoint x = 0 in order to reduce looping.
282 : *---------------------------------------------------------------------*/
283 :
284 2030780 : FOR( ; i < imid; i++ )
285 : {
286 1928213 : x2 = mult_r( x[i - 1], x[i - 1] );
287 :
288 1928213 : Mpy_32_16_ss( R[0], x2, &mh, &ml );
289 1928213 : re = L_add( R[2], mh );
290 1928213 : Mpy_32_16_ss( R[1], x2, &mh, &ml );
291 1928213 : ro = L_add( R[3], mh );
292 :
293 1928213 : Mpy_32_16_ss( S[0], x2, &mh, &ml );
294 1928213 : se = L_add( S[2], mh );
295 1928213 : Mpy_32_16_ss( S[1], x2, &mh, &ml );
296 1928213 : so = L_add( S[3], mh );
297 :
298 5784639 : FOR( j = 4; j < NC; j += 2 )
299 : {
300 3856426 : Mpy_32_16_ss( re, x2, &mh, &ml );
301 3856426 : re = L_add( R[j], mh );
302 3856426 : Mpy_32_16_ss( ro, x2, &mh, &ml );
303 3856426 : ro = L_add( R[j + 1], mh );
304 3856426 : Mpy_32_16_ss( se, x2, &mh, &ml );
305 3856426 : se = L_add( S[j], mh );
306 3856426 : Mpy_32_16_ss( so, x2, &mh, &ml );
307 3856426 : so = L_add( S[j + 1], mh );
308 : }
309 :
310 1928213 : Mpy_32_16_ss( re, x2, &mh, &ml );
311 1928213 : L_tmp = L_add( R[j], mh );
312 1928213 : Mpy_32_16_ss( ro, x[i - 1], &mh, &ml );
313 1928213 : re = L_add( L_tmp, mh );
314 1928213 : ro = L_sub( L_tmp, mh );
315 :
316 1928213 : Mpy_32_16_ss( se, x2, &mh, &ml );
317 1928213 : L_tmp = L_add( S[j], mh );
318 1928213 : Mpy_32_16_ss( so, x[i - 1], &mh, &ml );
319 1928213 : se = L_add( L_tmp, mh );
320 1928213 : so = L_sub( L_tmp, mh );
321 :
322 1928213 : G[i] = inv_pow( re, se, x[i - 1] );
323 1928213 : move32();
324 1928213 : G[N - i - 1] = inv_pow( so, ro, x[i - 1] );
325 1928213 : move32();
326 : }
327 :
328 102567 : return;
329 : }
330 :
331 199994 : static Word32 b_inv_sq(
332 : const Word32 in32, /* i : Input not normalized to inverse */
333 : const Word16 exp_in /* i : input current exponent */
334 : )
335 : {
336 : Word16 m_den, exp_den;
337 : Word16 div_out;
338 : Word32 Ltmp;
339 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
340 199994 : Flag Overflow = 0;
341 : #endif
342 :
343 199994 : exp_den = norm_l( in32 );
344 199994 : m_den = extract_h( L_shl( in32, exp_den ) );
345 199994 : exp_den = add( sub( 30, exp_den ), sub( 16, exp_in ) );
346 :
347 199994 : m_den = mult_r( m_den, m_den );
348 199994 : exp_den = shl_o( exp_den, 1, &Overflow );
349 :
350 199994 : div_out = div_s( 8192, m_den );
351 199994 : Ltmp = L_shl_o( div_out, add( sub( 30 - 13, exp_den ), 15 ), &Overflow ); /*Q15*/
352 :
353 199994 : return Ltmp;
354 : }
355 :
356 4005253 : static Word32 inv_pow(
357 : const Word32 re,
358 : const Word32 se,
359 : const Word16 x )
360 : {
361 : Word16 exp1, exp2;
362 : Word16 tmp;
363 : Word32 L_tmp;
364 : Word32 mh;
365 : UWord16 ml;
366 : Word32 r0, s0;
367 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
368 4005253 : Flag Overflow = 0;
369 : #endif
370 :
371 4005253 : IF( re == 0 )
372 : {
373 51 : exp1 = 30;
374 51 : move16();
375 51 : r0 = L_deposit_l( 0 );
376 : }
377 : ELSE
378 : {
379 4005202 : exp1 = norm_l( re );
380 4005202 : tmp = extract_h( L_shl( re, exp1 ) );
381 4005202 : L_tmp = L_shr( L_mult_o( tmp, tmp, &Overflow ), 1 );
382 4005202 : Mpy_32_16_ss( L_tmp, x, &mh, &ml );
383 4005202 : r0 = L_add( L_tmp, mh );
384 : }
385 :
386 4005253 : IF( se == 0 )
387 : {
388 18 : exp2 = 30;
389 18 : move16();
390 18 : s0 = L_deposit_l( 0 );
391 : }
392 : ELSE
393 : {
394 4005235 : exp2 = norm_l( se );
395 4005235 : tmp = extract_h( L_shl( se, exp2 ) );
396 4005235 : L_tmp = L_shr( L_mult_o( tmp, tmp, &Overflow ), 1 );
397 4005235 : Mpy_32_16_ss( L_tmp, x, &mh, &ml );
398 4005235 : s0 = L_sub( L_tmp, mh );
399 : }
400 :
401 4005253 : IF( exp1 > exp2 )
402 : {
403 2680026 : exp1 = shl( sub( exp1, exp2 ), 1 );
404 2680026 : r0 = L_shr( r0, exp1 );
405 :
406 2680026 : exp2 = add( add( exp2, exp2 ), 8 );
407 : }
408 : ELSE
409 : {
410 1325227 : exp2 = shl( sub( exp2, exp1 ), 1 );
411 1325227 : s0 = L_shr( s0, exp2 );
412 :
413 1325227 : exp2 = add( add( exp1, exp1 ), 8 );
414 : }
415 :
416 4005253 : r0 = L_add( r0, s0 );
417 4005253 : exp1 = norm_l( r0 );
418 4005253 : L_tmp = L_shl( r0, exp1 );
419 4005253 : tmp = extract_h( L_tmp );
420 4005253 : IF( tmp == 0 )
421 : {
422 0 : return MAX_32;
423 : }
424 4005253 : tmp = div_s( (Word16) ( ( 1 << 14 ) - 1 ), tmp );
425 4005253 : exp1 = add( exp1, exp2 );
426 4005253 : L_tmp = L_shr_o( tmp, sub( 31, exp1 ), &Overflow ); /* result in Q15 */
427 :
428 4005253 : return ( L_tmp );
429 : }
430 :
431 :
432 : /*---------------------------------------------------------------------*
433 : * spectautocorr()
434 : *
435 : * Computes the autocorrelation r[j] for j = 0, 1, ..., M from
436 : * the power spectrum P(w) by using rectangle rule to approximate
437 : * the integral
438 : *
439 : * 1 pi
440 : * r[j] = --- I P(w) cos(j*w) dw.
441 : * 2*pi -pi
442 : *
443 : * It is sufficient to evaluate the integrand only from w = 0 to
444 : * w = pi due to the symmetry P(-w) = P(w). We can further
445 : * employ the relation
446 : *
447 : * cos(j*(pi - w)) = (-1)^j cos(j*w)
448 : *
449 : * to use symmetries relative to w = pi/2.
450 : *
451 : * When applying the rectangle rule, it is useful to separate w = 0,
452 : * w = pi/2, and w = pi. By using a frequency grid of N points, we
453 : * can express the rectangle rule as
454 : *
455 : * r[j] = G[0] + 2*a*G[(N-1)/2] + b*G[N-1]
456 : *
457 : * M
458 : * + 2 sum (G[i] - G[N-i-1]) cos(j*x[i])
459 : * i=1
460 : *
461 : * where G[i] is the power spectrum at the grid point cos(i*pi/N)
462 : * and M = (N-1)/2 - 1 is the number of the grid points in the
463 : * interval(0, pi/2).
464 : *
465 : * The coefficients
466 : *
467 : * b = (-1)^j
468 : * a = (1 + (-1)^(j+1))(-1)^floor(j/2)
469 : *
470 : * follow from the properties of cosine. The computation further
471 : * uses the recursion
472 : *
473 : * cos(j*w) = 2*cos(w)*cos((j-1)*w) - cos((j-2)*w)
474 : *
475 : * Note that the autocorrelation can be scaled for convenience,
476 : * because this scaling has no impact on the LP coefficients to be
477 : * calculated from the autocorrelation. The expression of r[j] thus
478 : * omits the division by N.
479 : *
480 : * See the powerspect function on the definition of the grid.
481 : *
482 : * References
483 : * J. Makhoul, "Spectral linear prediction: properties and
484 : * applications," IEEE Trans. on Acoustics, Speech and Signal
485 : * Processing, Vol. 23, No. 3, pp.283-296, June 1975
486 : *---------------------------------------------------------------------*/
487 :
488 102567 : static void spectautocorr_fx(
489 : const Word16 x[], /* i: Grid points x[0:m-1] */
490 : const Word16 N, /* i: Number of grid points */
491 : const Word32 G[], /* i: Power spectrum G[0:N-1] */
492 : Word16 rh[], /* o: Autocorrelation r[0:M] */
493 : Word16 rl[] /* o: Autocorrelation r[0:M] */
494 : )
495 : {
496 : Word16 c[M + 1]; /* c[j] = cos(j*w) */
497 : Word32 gp, gn;
498 : Word16 i, j;
499 : Word16 imid;
500 : Word32 mh;
501 : UWord16 ml;
502 : Word32 r[M + 1];
503 : Word16 exp0;
504 : #ifdef BASOP_NOGLOB_DECLARE_LOCAL
505 102567 : Flag Overflow = 0;
506 : #endif
507 :
508 : /*---------------------------------------------------------------------*
509 : * The mid point of the cosine table x of m entries assuming an odd m.
510 : * Only the entries x[0] = cos(pi/m), x[1] = cos(2*pi/m), ...,
511 : * x[imid-1] = cos((imid-1)*pi/m) need to be stored due to trivial
512 : * cos(0), cos(pi/2), cos(pi), and symmetry relative to pi/2.
513 : * Here m = 51.
514 : *---------------------------------------------------------------------*/
515 :
516 102567 : imid = shr( sub( N, 1 ), 1 );
517 102567 : move16();
518 :
519 : /*---------------------------------------------------------------------*
520 : * Autocorrelation r[j] at zero lag j = 0 for the upper half of the
521 : * unit circle, but excluding the points x = cos(0) and x = cos(pi).
522 : *---------------------------------------------------------------------*/
523 :
524 102567 : r[0] = G[1];
525 102567 : move32();
526 4974383 : FOR( i = 2; i < N - 1; i++ )
527 : {
528 4871816 : r[0] = L_add_o( r[0], G[i], &Overflow );
529 4871816 : move32();
530 : }
531 :
532 : /*---------------------------------------------------------------------*
533 : * Initialize the autocorrelation r[j] at lags greater than zero
534 : * by adding the midpoint x = cos(pi/2) = 0.
535 : *---------------------------------------------------------------------*/
536 :
537 102567 : r[1] = L_deposit_l( 0 );
538 102567 : move32();
539 :
540 102567 : r[2] = -G[imid];
541 102567 : move32();
542 :
543 820536 : FOR( i = 3; i < M; i += 2 )
544 : {
545 717969 : r[i] = L_deposit_l( 0 );
546 717969 : r[i + 1] = L_negate( r[i - 1] );
547 717969 : move32();
548 717969 : move32();
549 : }
550 :
551 : /*---------------------------------------------------------------------*
552 : * Autocorrelation r[j] at lags j = 1, 2, ..., M. The computation
553 : * employes the relation cos(j*(pi - w)) = (-1)^j cos(j*w) and
554 : * cos(j*w) = 2*cos(w)*cos((j-1)*w) - cos((j-2)*w) for obtaining
555 : * the cosine c[j] = cos(j*w).
556 : *---------------------------------------------------------------------*/
557 :
558 102567 : c[0] = (Word16) 32767;
559 102567 : move16(); /* 1.0 in Q15 */
560 2538475 : FOR( i = 1; i < imid; i++ )
561 : {
562 2435908 : gp = L_add_o( G[i], G[N - i - 1], &Overflow );
563 2435908 : gn = L_sub( G[i], G[N - i - 1] );
564 :
565 : /*r[1] = L_mac(r[1], x[i-1], gn);*/
566 2435908 : Mpy_32_16_ss( gn, x[i - 1], &mh, &ml );
567 2435908 : r[1] = L_add_o( r[1], mh, &Overflow );
568 2435908 : move32();
569 2435908 : c[1] = x[i - 1];
570 2435908 : move16();
571 :
572 19487264 : FOR( j = 2; j < M; j += 2 )
573 : {
574 17051356 : c[j] = mult_r( c[j - 1], x[i - 1] );
575 17051356 : move16();
576 17051356 : c[j] = add_o( c[j], sub_o( c[j], c[j - 2], &Overflow ), &Overflow );
577 17051356 : move16();
578 :
579 : /*r[j] = L_mac(r[j], c[j], gp);*/
580 17051356 : Mpy_32_16_ss( gp, c[j], &mh, &ml );
581 17051356 : r[j] = L_add_o( r[j], mh, &Overflow );
582 17051356 : move32();
583 :
584 17051356 : c[j + 1] = mult_r( c[j], x[i - 1] );
585 17051356 : move16();
586 17051356 : c[j + 1] = add_o( c[j + 1], sub_o( c[j + 1], c[j - 1], &Overflow ), &Overflow );
587 17051356 : move16();
588 :
589 : /*r[j+1] = L_mac(r[j+1], c[j+1], gn);*/
590 17051356 : Mpy_32_16_ss( gn, c[j + 1], &mh, &ml );
591 17051356 : r[j + 1] = L_add_o( r[j + 1], mh, &Overflow );
592 17051356 : move32();
593 : }
594 2435908 : c[j] = mult_r( c[j - 1], x[i - 1] );
595 2435908 : move16();
596 2435908 : c[j] = add_o( c[j], sub_o( c[j], c[j - 2], &Overflow ), &Overflow );
597 2435908 : move16();
598 :
599 2435908 : Mpy_32_16_ss( gp, c[j], &mh, &ml );
600 2435908 : r[j] = L_add_o( r[j], mh, &Overflow );
601 2435908 : move32();
602 : }
603 :
604 : /*---------------------------------------------------------------------*
605 : * Add the endpoints x = cos(0) = 1 and x = cos(pi) = -1 as
606 : * well as the lower half of the unit circle.
607 : *---------------------------------------------------------------------*/
608 102567 : gp = L_shr( L_add_o( G[0], G[N - 1], &Overflow ), 1 );
609 102567 : gn = L_shr( L_sub( G[0], G[N - 1] ), 1 );
610 :
611 102567 : r[0] = L_add_o( r[0], gp, &Overflow );
612 102567 : move32();
613 102567 : exp0 = norm_l( r[0] );
614 102567 : L_Extract( L_shl( r[0], exp0 ), &rh[0], &rl[0] );
615 :
616 923103 : FOR( j = 1; j < M; j += 2 )
617 : {
618 820536 : L_Extract( L_shl( L_add_o( r[j], gn, &Overflow ), exp0 ), &rh[j], &rl[j] );
619 820536 : L_Extract( L_shl( L_add_o( r[j + 1], gp, &Overflow ), exp0 ), &rh[j + 1], &rl[j + 1] );
620 : }
621 :
622 102567 : return;
623 : }
624 :
625 : /*---------------------------------------------------------------------*
626 : * zeros2poly()
627 : *
628 : * Computes the coefficients of the polynomials
629 : *
630 : * R(x) = prod (x - x[i]),
631 : * i = 0,2,4,...
632 : *
633 : * S(x) = prod (x - x[i]),
634 : * i = 1,3,5,...
635 : *
636 : * when their zeros x[i] are given for i = 0, 1, ..., n-1. The
637 : * routine assumes n = 1 or even n greater than or equal to 4.
638 : *
639 : * The polynomial coefficients are returned in R[0:n/2-1] and
640 : * S[0:n/2-1]. The leading coefficients are in R[0] and S[0].
641 : *---------------------------------------------------------------------*/
642 102567 : static void zeros2poly_fx(
643 : Word16 x[], /* i: Q15 Zeros of R(x) and S(x) */
644 : Word32 R[], /* o: Q22 Coefficients of R(x) */
645 : Word32 S[] /* o: Q22 Coefficients of S(x) */
646 : )
647 : {
648 : Word16 xr, xs;
649 : Word16 i, j;
650 : Word32 mh;
651 : UWord16 ml;
652 :
653 102567 : R[0] = ( 1 << 27 ) - 1;
654 102567 : move32();
655 102567 : S[0] = ( 1 << 27 ) - 1;
656 102567 : move32();
657 102567 : R[1] = L_msu( 0, x[0], 1 << 11 );
658 102567 : move32();
659 102567 : S[1] = L_msu( 0, x[1], 1 << 11 );
660 102567 : move32();
661 :
662 820536 : FOR( i = 2; i <= NC; i++ )
663 : {
664 717969 : xr = negate( x[2 * i - 2] );
665 717969 : xs = negate( x[2 * i - 1] );
666 :
667 717969 : Mpy_32_16_ss( R[i - 1], xr, &R[i], &ml );
668 717969 : Mpy_32_16_ss( S[i - 1], xs, &S[i], &ml );
669 3589845 : FOR( j = i - 1; j > 0; j-- )
670 : {
671 2871876 : Mpy_32_16_ss( R[j - 1], xr, &mh, &ml );
672 2871876 : R[j] = L_add( R[j], mh );
673 2871876 : move32();
674 2871876 : Mpy_32_16_ss( S[j - 1], xs, &mh, &ml );
675 2871876 : S[j] = L_add( S[j], mh );
676 2871876 : move32();
677 : }
678 : }
679 :
680 102567 : return;
681 : }
682 :
683 : /*---------------------------------------------------------------------*
684 : * polydecomp()
685 : *
686 : * Computes the coefficients of the symmetric and antisymmetric
687 : * polynomials P(z) and Q(z) that define the line spectrum pair
688 : * decomposition of a given polynomial A(z) of order n. For even n,
689 : *
690 : * P(z) = [A(z) + z^(n+1) A(1/z)]/(1/z + 1),
691 : * Q(z) = [A(z) - z^(n+1) A(1/z)]/(1/z - 1),
692 : *
693 : * These polynomials are then expressed in their direct form,
694 : * respectively, R(x) and S(x), on the real axis x = cos w using
695 : * explicit Chebyshev polynomials of the first kind.
696 : *
697 : * The coefficients of the polynomials R(x) and S(x) are returned
698 : * in R[0:n/2] and S[0:n/2] for the given linear prediction
699 : * coefficients A[0:n/2]. Note that R(x) and S(x) are formed in
700 : * place such that P(z) is stored in the same array than R(x),
701 : * and Q(z) is stored in the same array than S(x).
702 : *
703 : * The routines assumes n = 16.
704 : *---------------------------------------------------------------------*/
705 :
706 0 : static void polydecomp_fx(
707 : Word16 A[], /* i: Q12 linear prediction coefficients */
708 : Word32 R[], /* o: Q20 coefficients of R(x) */
709 : Word32 S[] /* o: Q20 coefficients of S(x) */
710 : )
711 : {
712 : Word16 scale;
713 : Word16 i;
714 : Word32 Ltmp1, Ltmp2;
715 :
716 0 : scale = shl( ( 1 << 5 ), norm_s( A[0] ) );
717 :
718 0 : R[0] = ( 1 << 20 ) - 1;
719 0 : move32(); /* Q20 */
720 0 : S[0] = ( 1 << 20 ) - 1;
721 0 : move32();
722 :
723 0 : FOR( i = 0; i < NC; i++ )
724 : {
725 0 : Ltmp1 = L_mult( A[i + 1], scale );
726 :
727 0 : Ltmp2 = L_msu( Ltmp1, A[M - i], scale );
728 0 : Ltmp1 = L_mac( Ltmp1, A[M - i], scale );
729 :
730 0 : R[i + 1] = L_sub( Ltmp1, R[i] );
731 0 : move32();
732 0 : S[i + 1] = L_add( Ltmp2, S[i] );
733 0 : move32();
734 : }
735 :
736 0 : cheb2poly_fx( R );
737 0 : cheb2poly_fx( S );
738 :
739 0 : return;
740 : }
741 :
742 : /*---------------------------------------------------------------------*
743 : * cheb2poly_fx()
744 : *
745 : * Computes the coefficients of the explicit Chebyshev polynomial
746 : * P(x) = P[0]*x^n + P[1]*x^(n-1) + ... + P[n] given the coefficients
747 : * of the series
748 : *
749 : * C(x) = C[0]*T_n(x) + C[1]*T_n-1(x) + ... + C[n]*T_0(x)
750 : *
751 : * where T_n(x) is the nth Chebyshev polynomial of the first kind.
752 : * This implementation assumes C[0] = 1. Only value n = 8 is
753 : * supported.
754 : *
755 : * The conversion from C(x) to P(x) is done in place such that the
756 : * coefficients of C(x) are given in P[0:8] and those of P(x) are
757 : * returned in the same array.
758 : *---------------------------------------------------------------------*/
759 :
760 0 : static void cheb2poly_fx(
761 : Word32 L_P[] /* i/o Q20: The coefficients of C(x) and P(x) */
762 : )
763 : {
764 : Word32 L_C[NC + 1], L_tmp;
765 : Word16 i;
766 :
767 0 : FOR( i = 1; i <= NC; i++ )
768 : {
769 0 : L_C[i] = L_P[i];
770 0 : move32();
771 : }
772 :
773 0 : L_P[0] = ( 1 << 27 ) - 1;
774 0 : move32();
775 :
776 0 : L_P[1] = L_shl( L_C[1], 6 );
777 0 : move32(); /* 64.0*C[1] */
778 0 : L_P[8] = L_shl( L_C[1], 3 );
779 0 : move32(); /* 8.0*C[1] */
780 :
781 0 : L_P[5] = L_sub( L_P[1], L_P[8] );
782 0 : move32(); /* 56.0*C[1] */
783 :
784 0 : L_P[2] = L_shl( L_C[3], 2 );
785 0 : move32();
786 0 : L_tmp = L_add( L_C[3], L_P[2] ); /* 5.0*C[3] */
787 0 : L_P[7] = L_sub( L_tmp, L_sub( L_P[8], L_C[1] ) );
788 0 : move32(); /* -7.0*C[1] */
789 :
790 0 : L_P[8] = L_shl( L_C[3], 4 );
791 0 : move32(); /* 16.0*C[3] */
792 0 : L_P[3] = L_sub( L_P[8], L_shl( L_P[5], 1 ) );
793 0 : move32(); /*-112.0*C[1] */
794 :
795 0 : L_P[5] = L_sub( L_P[5], L_add( L_P[8], L_P[2] ) );
796 0 : move32(); /* -20.0*C[3] */
797 :
798 0 : L_P[2] = L_shl( L_C[5], 2 );
799 0 : move32();
800 0 : L_P[5] = L_add( L_P[5], L_P[2] );
801 0 : move32(); /* 4.0*C[5] */
802 :
803 0 : L_tmp = L_sub( L_P[7], L_sub( L_P[2], L_C[5] ) ); /* -3.0*C[5] */
804 0 : L_P[7] = L_add( L_tmp, L_C[7] );
805 0 : move32(); /* C[7] */
806 :
807 0 : L_P[6] = L_shl( L_C[2], 4 );
808 0 : move32();
809 0 : L_tmp = L_sub( ( 160 << 20 ), L_P[6] );
810 0 : L_P[4] = L_sub( L_tmp, L_shl( L_C[2], 5 ) );
811 0 : move32(); /* -48.0*C[2] */
812 :
813 0 : L_tmp = L_add( L_P[6], L_shl( L_C[2], 1 ) ); /* 18.0*C[2] */
814 0 : L_P[6] = L_sub( L_tmp, ( 32 << 20 ) );
815 0 : move32();
816 :
817 0 : L_P[8] = L_shl( L_C[4], 3 );
818 0 : move32();
819 0 : L_P[4] = L_add( L_P[4], L_P[8] );
820 0 : move32(); /* 8.0*C[4] */
821 :
822 0 : L_tmp = L_sub( L_P[6], L_P[8] ); /* -8.0*C[4] */
823 0 : L_P[6] = L_add( L_tmp, L_shl( L_C[6], 1 ) );
824 0 : move32(); /* 2.0*C[6] */
825 :
826 0 : L_tmp = L_shl( L_C[2], 5 ); /* 32.0*C[2] */
827 0 : L_P[2] = L_sub( L_tmp, ( 256 << 20 ) );
828 0 : move32();
829 :
830 0 : L_tmp = L_add( ( 1 << 21 ) + 1, L_C[8] );
831 0 : L_tmp = L_shr( L_tmp, 1 ); /* 1+0.5*C[8] */
832 0 : L_tmp = L_sub( L_tmp, L_C[2] );
833 0 : L_tmp = L_add( L_tmp, L_C[4] );
834 0 : L_P[8] = L_sub( L_tmp, L_C[6] );
835 0 : move32();
836 :
837 0 : return;
838 : }
|
