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