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 103344 : 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 103344 : 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 103344 : E_LPC_f_lsp_a_conversion( w, oldA, M );
93 :
94 103344 : 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 103344 : IF( EQ_16( L_frame, L_FRAME ) )
104 : {
105 5226 : powerspect_fx( grid50_fx, N50, R, S, G, DOWNCONV );
106 5226 : 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 98118 : powerspect_fx( grid40_fx, N40, R, S, G, UPCONV );
118 :
119 1079298 : FOR( i = N40; i < N50; i++ )
120 : {
121 981180 : G[i] = G[N40 - 1];
122 981180 : move32();
123 : }
124 :
125 98118 : 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 103344 : flag = E_LPC_lev_dur( rh, rl, A, epsP, M, oldA );
134 103344 : E_LPC_a_lsp_conversion( A, w, stable_LSP_fx, M );
135 :
136 :
137 103344 : 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 103344 : 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 103344 : IF( mode == DOWNCONV )
199 : {
200 5226 : iuni = ( GRID50_POINTS - 1 ) / 5 - 1;
201 5226 : move16();
202 5226 : imid = ( GRID50_POINTS - 1 ) / 2;
203 5226 : 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 98118 : iuni = 0;
214 98118 : move16();
215 98118 : imid = ( GRID40_POINTS - 1 ) / 2;
216 98118 : move16();
217 :
218 98118 : G[N - 1] = S[0];
219 98118 : move32();
220 :
221 883062 : FOR( j = 1; j <= NC; j++ )
222 : {
223 784944 : G[N - 1] = L_sub( S[j], G[N - 1] );
224 784944 : move32();
225 : }
226 98118 : G[N - 1] = b_inv_sq( G[N - 1], 19 );
227 98118 : move32();
228 : }
229 :
230 : /*---------------------------------------------------------------------*
231 : * Power spectrum x = cos(0) = 1.
232 : *---------------------------------------------------------------------*/
233 :
234 103344 : G[0] = R[0];
235 103344 : move32();
236 930096 : FOR( j = 1; j <= NC; j++ )
237 : {
238 826752 : G[0] = L_add( R[j], G[0] );
239 826752 : move32();
240 : }
241 :
242 103344 : G[0] = b_inv_sq( L_max( G[0], 1 ), 19 );
243 103344 : move32();
244 :
245 : /*---------------------------------------------------------------------*
246 : * Power spectrum at x = cos(pi/2) = 0.
247 : *---------------------------------------------------------------------*/
248 103344 : G[imid] = inv_pow( R[NC], S[NC], 0 );
249 103344 : 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 150378 : FOR( i = 1; i <= iuni; i++ )
257 : {
258 47034 : Mpy_32_16_ss( R[0], x[i - 1], &mh, &ml );
259 47034 : r0 = L_add( R[1], mh );
260 47034 : Mpy_32_16_ss( S[0], x[i - 1], &mh, &ml );
261 47034 : s0 = L_add( S[1], mh );
262 :
263 :
264 376272 : FOR( j = 2; j <= NC; j++ )
265 : {
266 329238 : Mpy_32_16_ss( r0, x[i - 1], &mh, &ml );
267 329238 : r0 = L_add( R[j], mh );
268 :
269 329238 : Mpy_32_16_ss( s0, x[i - 1], &mh, &ml );
270 329238 : s0 = L_add( S[j], mh );
271 : }
272 :
273 47034 : G[i] = inv_pow( r0, s0, x[i - 1] );
274 47034 : 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 2045976 : FOR( ; i < imid; i++ )
284 : {
285 1942632 : x2 = mult_r( x[i - 1], x[i - 1] );
286 :
287 1942632 : Mpy_32_16_ss( R[0], x2, &mh, &ml );
288 1942632 : re = L_add( R[2], mh );
289 1942632 : Mpy_32_16_ss( R[1], x2, &mh, &ml );
290 1942632 : ro = L_add( R[3], mh );
291 :
292 1942632 : Mpy_32_16_ss( S[0], x2, &mh, &ml );
293 1942632 : se = L_add( S[2], mh );
294 1942632 : Mpy_32_16_ss( S[1], x2, &mh, &ml );
295 1942632 : so = L_add( S[3], mh );
296 :
297 5827896 : FOR( j = 4; j < NC; j += 2 )
298 : {
299 3885264 : Mpy_32_16_ss( re, x2, &mh, &ml );
300 3885264 : re = L_add( R[j], mh );
301 3885264 : Mpy_32_16_ss( ro, x2, &mh, &ml );
302 3885264 : ro = L_add( R[j + 1], mh );
303 3885264 : Mpy_32_16_ss( se, x2, &mh, &ml );
304 3885264 : se = L_add( S[j], mh );
305 3885264 : Mpy_32_16_ss( so, x2, &mh, &ml );
306 3885264 : so = L_add( S[j + 1], mh );
307 : }
308 :
309 1942632 : Mpy_32_16_ss( re, x2, &mh, &ml );
310 1942632 : L_tmp = L_add( R[j], mh );
311 1942632 : Mpy_32_16_ss( ro, x[i - 1], &mh, &ml );
312 1942632 : re = L_add( L_tmp, mh );
313 1942632 : ro = L_sub( L_tmp, mh );
314 :
315 1942632 : Mpy_32_16_ss( se, x2, &mh, &ml );
316 1942632 : L_tmp = L_add( S[j], mh );
317 1942632 : Mpy_32_16_ss( so, x[i - 1], &mh, &ml );
318 1942632 : se = L_add( L_tmp, mh );
319 1942632 : so = L_sub( L_tmp, mh );
320 :
321 1942632 : G[i] = inv_pow( re, se, x[i - 1] );
322 1942632 : move32();
323 1942632 : G[N - i - 1] = inv_pow( so, ro, x[i - 1] );
324 1942632 : move32();
325 : }
326 :
327 103344 : return;
328 : }
329 :
330 201462 : 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 :
339 201462 : exp_den = norm_l( in32 );
340 201462 : m_den = extract_h( L_shl( in32, exp_den ) );
341 201462 : exp_den = add( sub( 30, exp_den ), sub( 16, exp_in ) );
342 :
343 201462 : m_den = mult_r( m_den, m_den );
344 201462 : exp_den = shl( exp_den, 1 );
345 :
346 201462 : div_out = div_s( 8192, m_den );
347 201462 : Ltmp = L_shl_sat( div_out, add( sub( 30 - 13, exp_den ), 15 ) ); /*Q15*/
348 :
349 201462 : return Ltmp;
350 : }
351 :
352 4035642 : static Word32 inv_pow(
353 : const Word32 re,
354 : const Word32 se,
355 : const Word16 x )
356 : {
357 : Word16 exp1, exp2;
358 : Word16 tmp;
359 : Word32 L_tmp;
360 : Word32 mh;
361 : UWord16 ml;
362 : Word32 r0, s0;
363 :
364 4035642 : IF( re == 0 )
365 : {
366 52 : exp1 = 30;
367 52 : move16();
368 52 : r0 = L_deposit_l( 0 );
369 : }
370 : ELSE
371 : {
372 4035590 : exp1 = norm_l( re );
373 4035590 : tmp = extract_h( L_shl( re, exp1 ) );
374 4035590 : L_tmp = L_shr( L_mult_sat( tmp, tmp ), 1 );
375 4035590 : Mpy_32_16_ss( L_tmp, x, &mh, &ml );
376 4035590 : r0 = L_add( L_tmp, mh );
377 : }
378 :
379 4035642 : IF( se == 0 )
380 : {
381 12 : exp2 = 30;
382 12 : move16();
383 12 : s0 = L_deposit_l( 0 );
384 : }
385 : ELSE
386 : {
387 4035630 : exp2 = norm_l( se );
388 4035630 : tmp = extract_h( L_shl( se, exp2 ) );
389 4035630 : L_tmp = L_shr( L_mult_sat( tmp, tmp ), 1 );
390 4035630 : Mpy_32_16_ss( L_tmp, x, &mh, &ml );
391 4035630 : s0 = L_sub( L_tmp, mh );
392 : }
393 :
394 4035642 : IF( exp1 > exp2 )
395 : {
396 2699730 : exp1 = shl( sub( exp1, exp2 ), 1 );
397 2699730 : r0 = L_shr( r0, exp1 );
398 :
399 2699730 : exp2 = add( add( exp2, exp2 ), 8 );
400 : }
401 : ELSE
402 : {
403 1335912 : exp2 = shl( sub( exp2, exp1 ), 1 );
404 1335912 : s0 = L_shr( s0, exp2 );
405 :
406 1335912 : exp2 = add( add( exp1, exp1 ), 8 );
407 : }
408 :
409 4035642 : r0 = L_add( r0, s0 );
410 4035642 : exp1 = norm_l( r0 );
411 4035642 : L_tmp = L_shl( r0, exp1 );
412 4035642 : tmp = extract_h( L_tmp );
413 4035642 : IF( tmp == 0 )
414 : {
415 0 : return MAX_32;
416 : }
417 4035642 : tmp = div_s( (Word16) ( ( 1 << 14 ) - 1 ), tmp );
418 4035642 : exp1 = add( exp1, exp2 );
419 4035642 : L_tmp = L_shr_sat( tmp, sub( 31, exp1 ) ); /* result in Q15 */
420 :
421 4035642 : return ( L_tmp );
422 : }
423 :
424 :
425 : /*---------------------------------------------------------------------*
426 : * spectautocorr()
427 : *
428 : * Computes the autocorrelation r[j] for j = 0, 1, ..., M from
429 : * the power spectrum P(w) by using rectangle rule to approximate
430 : * the integral
431 : *
432 : * 1 pi
433 : * r[j] = --- I P(w) cos(j*w) dw.
434 : * 2*pi -pi
435 : *
436 : * It is sufficient to evaluate the integrand only from w = 0 to
437 : * w = pi due to the symmetry P(-w) = P(w). We can further
438 : * employ the relation
439 : *
440 : * cos(j*(pi - w)) = (-1)^j cos(j*w)
441 : *
442 : * to use symmetries relative to w = pi/2.
443 : *
444 : * When applying the rectangle rule, it is useful to separate w = 0,
445 : * w = pi/2, and w = pi. By using a frequency grid of N points, we
446 : * can express the rectangle rule as
447 : *
448 : * r[j] = G[0] + 2*a*G[(N-1)/2] + b*G[N-1]
449 : *
450 : * M
451 : * + 2 sum (G[i] - G[N-i-1]) cos(j*x[i])
452 : * i=1
453 : *
454 : * where G[i] is the power spectrum at the grid point cos(i*pi/N)
455 : * and M = (N-1)/2 - 1 is the number of the grid points in the
456 : * interval(0, pi/2).
457 : *
458 : * The coefficients
459 : *
460 : * b = (-1)^j
461 : * a = (1 + (-1)^(j+1))(-1)^floor(j/2)
462 : *
463 : * follow from the properties of cosine. The computation further
464 : * uses the recursion
465 : *
466 : * cos(j*w) = 2*cos(w)*cos((j-1)*w) - cos((j-2)*w)
467 : *
468 : * Note that the autocorrelation can be scaled for convenience,
469 : * because this scaling has no impact on the LP coefficients to be
470 : * calculated from the autocorrelation. The expression of r[j] thus
471 : * omits the division by N.
472 : *
473 : * See the powerspect function on the definition of the grid.
474 : *
475 : * References
476 : * J. Makhoul, "Spectral linear prediction: properties and
477 : * applications," IEEE Trans. on Acoustics, Speech and Signal
478 : * Processing, Vol. 23, No. 3, pp.283-296, June 1975
479 : *---------------------------------------------------------------------*/
480 :
481 103344 : static void spectautocorr_fx(
482 : const Word16 x[], /* i: Grid points x[0:m-1] */
483 : const Word16 N, /* i: Number of grid points */
484 : const Word32 G[], /* i: Power spectrum G[0:N-1] */
485 : Word16 rh[], /* o: Autocorrelation r[0:M] */
486 : Word16 rl[] /* o: Autocorrelation r[0:M] */
487 : )
488 : {
489 : Word16 c[M + 1]; /* c[j] = cos(j*w) */
490 : Word32 gp, gn;
491 : Word16 i, j;
492 : Word16 imid;
493 : Word32 mh;
494 : UWord16 ml;
495 : Word32 r[M + 1];
496 : Word16 exp0;
497 :
498 : /*---------------------------------------------------------------------*
499 : * The mid point of the cosine table x of m entries assuming an odd m.
500 : * Only the entries x[0] = cos(pi/m), x[1] = cos(2*pi/m), ...,
501 : * x[imid-1] = cos((imid-1)*pi/m) need to be stored due to trivial
502 : * cos(0), cos(pi/2), cos(pi), and symmetry relative to pi/2.
503 : * Here m = 51.
504 : *---------------------------------------------------------------------*/
505 :
506 103344 : imid = shr( sub( N, 1 ), 1 );
507 103344 : move16();
508 :
509 : /*---------------------------------------------------------------------*
510 : * Autocorrelation r[j] at zero lag j = 0 for the upper half of the
511 : * unit circle, but excluding the points x = cos(0) and x = cos(pi).
512 : *---------------------------------------------------------------------*/
513 :
514 103344 : r[0] = G[1];
515 103344 : move32();
516 5011596 : FOR( i = 2; i < N - 1; i++ )
517 : {
518 4908252 : r[0] = L_add_sat( r[0], G[i] );
519 4908252 : move32();
520 : }
521 :
522 : /*---------------------------------------------------------------------*
523 : * Initialize the autocorrelation r[j] at lags greater than zero
524 : * by adding the midpoint x = cos(pi/2) = 0.
525 : *---------------------------------------------------------------------*/
526 :
527 103344 : r[1] = L_deposit_l( 0 );
528 103344 : move32();
529 :
530 103344 : r[2] = -G[imid];
531 103344 : move32();
532 :
533 826752 : FOR( i = 3; i < M; i += 2 )
534 : {
535 723408 : r[i] = L_deposit_l( 0 );
536 723408 : r[i + 1] = L_negate( r[i - 1] );
537 723408 : move32();
538 723408 : move32();
539 : }
540 :
541 : /*---------------------------------------------------------------------*
542 : * Autocorrelation r[j] at lags j = 1, 2, ..., M. The computation
543 : * employes the relation cos(j*(pi - w)) = (-1)^j cos(j*w) and
544 : * cos(j*w) = 2*cos(w)*cos((j-1)*w) - cos((j-2)*w) for obtaining
545 : * the cosine c[j] = cos(j*w).
546 : *---------------------------------------------------------------------*/
547 :
548 103344 : c[0] = (Word16) 32767;
549 103344 : move16(); /* 1.0 in Q15 */
550 2557470 : FOR( i = 1; i < imid; i++ )
551 : {
552 2454126 : gp = L_add_sat( G[i], G[N - i - 1] );
553 2454126 : gn = L_sub( G[i], G[N - i - 1] );
554 :
555 : /*r[1] = L_mac(r[1], x[i-1], gn);*/
556 2454126 : Mpy_32_16_ss( gn, x[i - 1], &mh, &ml );
557 2454126 : r[1] = L_add_sat( r[1], mh );
558 2454126 : move32();
559 2454126 : c[1] = x[i - 1];
560 2454126 : move16();
561 :
562 19633008 : FOR( j = 2; j < M; j += 2 )
563 : {
564 17178882 : c[j] = mult_r( c[j - 1], x[i - 1] );
565 17178882 : move16();
566 17178882 : c[j] = add_sat( c[j], sub_sat( c[j], c[j - 2] ) );
567 17178882 : move16();
568 :
569 : /*r[j] = L_mac(r[j], c[j], gp);*/
570 17178882 : Mpy_32_16_ss( gp, c[j], &mh, &ml );
571 17178882 : r[j] = L_add_sat( r[j], mh );
572 17178882 : move32();
573 :
574 17178882 : c[j + 1] = mult_r( c[j], x[i - 1] );
575 17178882 : move16();
576 17178882 : c[j + 1] = add_sat( c[j + 1], sub_sat( c[j + 1], c[j - 1] ) );
577 17178882 : move16();
578 :
579 : /*r[j+1] = L_mac(r[j+1], c[j+1], gn);*/
580 17178882 : Mpy_32_16_ss( gn, c[j + 1], &mh, &ml );
581 17178882 : r[j + 1] = L_add_sat( r[j + 1], mh );
582 17178882 : move32();
583 : }
584 2454126 : c[j] = mult_r( c[j - 1], x[i - 1] );
585 2454126 : move16();
586 2454126 : c[j] = add_sat( c[j], sub_sat( c[j], c[j - 2] ) );
587 2454126 : move16();
588 :
589 2454126 : Mpy_32_16_ss( gp, c[j], &mh, &ml );
590 2454126 : r[j] = L_add_sat( r[j], mh );
591 2454126 : move32();
592 : }
593 :
594 : /*---------------------------------------------------------------------*
595 : * Add the endpoints x = cos(0) = 1 and x = cos(pi) = -1 as
596 : * well as the lower half of the unit circle.
597 : *---------------------------------------------------------------------*/
598 :
599 103344 : gp = L_shr( L_add_sat( G[0], G[N - 1] ), 1 );
600 103344 : gn = L_shr( L_sub( G[0], G[N - 1] ), 1 );
601 :
602 103344 : r[0] = L_add_sat( r[0], gp );
603 103344 : move32();
604 103344 : exp0 = norm_l( r[0] );
605 103344 : L_Extract( L_shl( r[0], exp0 ), &rh[0], &rl[0] );
606 :
607 930096 : FOR( j = 1; j < M; j += 2 )
608 : {
609 826752 : L_Extract( L_shl( L_add_sat( r[j], gn ), exp0 ), &rh[j], &rl[j] );
610 826752 : L_Extract( L_shl( L_add_sat( r[j + 1], gp ), exp0 ), &rh[j + 1], &rl[j + 1] );
611 : }
612 :
613 103344 : return;
614 : }
615 :
616 : /*---------------------------------------------------------------------*
617 : * zeros2poly()
618 : *
619 : * Computes the coefficients of the polynomials
620 : *
621 : * R(x) = prod (x - x[i]),
622 : * i = 0,2,4,...
623 : *
624 : * S(x) = prod (x - x[i]),
625 : * i = 1,3,5,...
626 : *
627 : * when their zeros x[i] are given for i = 0, 1, ..., n-1. The
628 : * routine assumes n = 1 or even n greater than or equal to 4.
629 : *
630 : * The polynomial coefficients are returned in R[0:n/2-1] and
631 : * S[0:n/2-1]. The leading coefficients are in R[0] and S[0].
632 : *---------------------------------------------------------------------*/
633 103344 : static void zeros2poly_fx(
634 : Word16 x[], /* i: Q15 Zeros of R(x) and S(x) */
635 : Word32 R[], /* o: Q22 Coefficients of R(x) */
636 : Word32 S[] /* o: Q22 Coefficients of S(x) */
637 : )
638 : {
639 : Word16 xr, xs;
640 : Word16 i, j;
641 : Word32 mh;
642 : UWord16 ml;
643 :
644 103344 : R[0] = ( 1 << 27 ) - 1;
645 103344 : move32();
646 103344 : S[0] = ( 1 << 27 ) - 1;
647 103344 : move32();
648 103344 : R[1] = L_msu( 0, x[0], 1 << 11 );
649 103344 : move32();
650 103344 : S[1] = L_msu( 0, x[1], 1 << 11 );
651 103344 : move32();
652 :
653 826752 : FOR( i = 2; i <= NC; i++ )
654 : {
655 723408 : xr = negate( x[2 * i - 2] );
656 723408 : xs = negate( x[2 * i - 1] );
657 :
658 723408 : Mpy_32_16_ss( R[i - 1], xr, &R[i], &ml );
659 723408 : Mpy_32_16_ss( S[i - 1], xs, &S[i], &ml );
660 3617040 : FOR( j = i - 1; j > 0; j-- )
661 : {
662 2893632 : Mpy_32_16_ss( R[j - 1], xr, &mh, &ml );
663 2893632 : R[j] = L_add( R[j], mh );
664 2893632 : move32();
665 2893632 : Mpy_32_16_ss( S[j - 1], xs, &mh, &ml );
666 2893632 : S[j] = L_add( S[j], mh );
667 2893632 : move32();
668 : }
669 : }
670 :
671 103344 : return;
672 : }
673 :
674 : /*---------------------------------------------------------------------*
675 : * polydecomp()
676 : *
677 : * Computes the coefficients of the symmetric and antisymmetric
678 : * polynomials P(z) and Q(z) that define the line spectrum pair
679 : * decomposition of a given polynomial A(z) of order n. For even n,
680 : *
681 : * P(z) = [A(z) + z^(n+1) A(1/z)]/(1/z + 1),
682 : * Q(z) = [A(z) - z^(n+1) A(1/z)]/(1/z - 1),
683 : *
684 : * These polynomials are then expressed in their direct form,
685 : * respectively, R(x) and S(x), on the real axis x = cos w using
686 : * explicit Chebyshev polynomials of the first kind.
687 : *
688 : * The coefficients of the polynomials R(x) and S(x) are returned
689 : * in R[0:n/2] and S[0:n/2] for the given linear prediction
690 : * coefficients A[0:n/2]. Note that R(x) and S(x) are formed in
691 : * place such that P(z) is stored in the same array than R(x),
692 : * and Q(z) is stored in the same array than S(x).
693 : *
694 : * The routines assumes n = 16.
695 : *---------------------------------------------------------------------*/
696 :
697 0 : static void polydecomp_fx(
698 : Word16 A[], /* i: Q12 linear prediction coefficients */
699 : Word32 R[], /* o: Q20 coefficients of R(x) */
700 : Word32 S[] /* o: Q20 coefficients of S(x) */
701 : )
702 : {
703 : Word16 scale;
704 : Word16 i;
705 : Word32 Ltmp1, Ltmp2;
706 :
707 0 : scale = shl( ( 1 << 5 ), norm_s( A[0] ) );
708 :
709 0 : R[0] = ( 1 << 20 ) - 1;
710 0 : move32(); /* Q20 */
711 0 : S[0] = ( 1 << 20 ) - 1;
712 0 : move32();
713 :
714 0 : FOR( i = 0; i < NC; i++ )
715 : {
716 0 : Ltmp1 = L_mult( A[i + 1], scale );
717 :
718 0 : Ltmp2 = L_msu( Ltmp1, A[M - i], scale );
719 0 : Ltmp1 = L_mac( Ltmp1, A[M - i], scale );
720 :
721 0 : R[i + 1] = L_sub( Ltmp1, R[i] );
722 0 : move32();
723 0 : S[i + 1] = L_add( Ltmp2, S[i] );
724 0 : move32();
725 : }
726 :
727 0 : cheb2poly_fx( R );
728 0 : cheb2poly_fx( S );
729 :
730 0 : return;
731 : }
732 :
733 : /*---------------------------------------------------------------------*
734 : * cheb2poly_fx()
735 : *
736 : * Computes the coefficients of the explicit Chebyshev polynomial
737 : * P(x) = P[0]*x^n + P[1]*x^(n-1) + ... + P[n] given the coefficients
738 : * of the series
739 : *
740 : * C(x) = C[0]*T_n(x) + C[1]*T_n-1(x) + ... + C[n]*T_0(x)
741 : *
742 : * where T_n(x) is the nth Chebyshev polynomial of the first kind.
743 : * This implementation assumes C[0] = 1. Only value n = 8 is
744 : * supported.
745 : *
746 : * The conversion from C(x) to P(x) is done in place such that the
747 : * coefficients of C(x) are given in P[0:8] and those of P(x) are
748 : * returned in the same array.
749 : *---------------------------------------------------------------------*/
750 :
751 0 : static void cheb2poly_fx(
752 : Word32 L_P[] /* i/o Q20: The coefficients of C(x) and P(x) */
753 : )
754 : {
755 : Word32 L_C[NC + 1], L_tmp;
756 : Word16 i;
757 :
758 0 : FOR( i = 1; i <= NC; i++ )
759 : {
760 0 : L_C[i] = L_P[i];
761 0 : move32();
762 : }
763 :
764 0 : L_P[0] = ( 1 << 27 ) - 1;
765 0 : move32();
766 :
767 0 : L_P[1] = L_shl( L_C[1], 6 );
768 0 : move32(); /* 64.0*C[1] */
769 0 : L_P[8] = L_shl( L_C[1], 3 );
770 0 : move32(); /* 8.0*C[1] */
771 :
772 0 : L_P[5] = L_sub( L_P[1], L_P[8] );
773 0 : move32(); /* 56.0*C[1] */
774 :
775 0 : L_P[2] = L_shl( L_C[3], 2 );
776 0 : move32();
777 0 : L_tmp = L_add( L_C[3], L_P[2] ); /* 5.0*C[3] */
778 0 : L_P[7] = L_sub( L_tmp, L_sub( L_P[8], L_C[1] ) );
779 0 : move32(); /* -7.0*C[1] */
780 :
781 0 : L_P[8] = L_shl( L_C[3], 4 );
782 0 : move32(); /* 16.0*C[3] */
783 0 : L_P[3] = L_sub( L_P[8], L_shl( L_P[5], 1 ) );
784 0 : move32(); /*-112.0*C[1] */
785 :
786 0 : L_P[5] = L_sub( L_P[5], L_add( L_P[8], L_P[2] ) );
787 0 : move32(); /* -20.0*C[3] */
788 :
789 0 : L_P[2] = L_shl( L_C[5], 2 );
790 0 : move32();
791 0 : L_P[5] = L_add( L_P[5], L_P[2] );
792 0 : move32(); /* 4.0*C[5] */
793 :
794 0 : L_tmp = L_sub( L_P[7], L_sub( L_P[2], L_C[5] ) ); /* -3.0*C[5] */
795 0 : L_P[7] = L_add( L_tmp, L_C[7] );
796 0 : move32(); /* C[7] */
797 :
798 0 : L_P[6] = L_shl( L_C[2], 4 );
799 0 : move32();
800 0 : L_tmp = L_sub( ( 160 << 20 ), L_P[6] );
801 0 : L_P[4] = L_sub( L_tmp, L_shl( L_C[2], 5 ) );
802 0 : move32(); /* -48.0*C[2] */
803 :
804 0 : L_tmp = L_add( L_P[6], L_shl( L_C[2], 1 ) ); /* 18.0*C[2] */
805 0 : L_P[6] = L_sub( L_tmp, ( 32 << 20 ) );
806 0 : move32();
807 :
808 0 : L_P[8] = L_shl( L_C[4], 3 );
809 0 : move32();
810 0 : L_P[4] = L_add( L_P[4], L_P[8] );
811 0 : move32(); /* 8.0*C[4] */
812 :
813 0 : L_tmp = L_sub( L_P[6], L_P[8] ); /* -8.0*C[4] */
814 0 : L_P[6] = L_add( L_tmp, L_shl( L_C[6], 1 ) );
815 0 : move32(); /* 2.0*C[6] */
816 :
817 0 : L_tmp = L_shl( L_C[2], 5 ); /* 32.0*C[2] */
818 0 : L_P[2] = L_sub( L_tmp, ( 256 << 20 ) );
819 0 : move32();
820 :
821 0 : L_tmp = L_add( ( 1 << 21 ) + 1, L_C[8] );
822 0 : L_tmp = L_shr( L_tmp, 1 ); /* 1+0.5*C[8] */
823 0 : L_tmp = L_sub( L_tmp, L_C[2] );
824 0 : L_tmp = L_add( L_tmp, L_C[4] );
825 0 : L_P[8] = L_sub( L_tmp, L_C[6] );
826 0 : move32();
827 :
828 0 : return;
829 : }
|
