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 <math.h>
43 : #include "cnst.h"
44 : #include "prot.h"
45 : #include "rom_com.h"
46 : #include "basop_proto_func.h"
47 : #include "wmc_auto.h"
48 :
49 : /*-------------------------------------------------------------------*
50 : * Local function prototypes
51 : *-------------------------------------------------------------------*/
52 :
53 : static float chebps2( const float x, const float *f, const int16_t n );
54 : static float LPC_chebyshev( float x, float f[], const int16_t n );
55 : static void get_isppol( const float *isp, float f[], const int16_t n );
56 :
57 : /*---------------------------------------------------------------------*
58 : * a2isp()
59 : *
60 : * Compute the ISPs from the LPC coefficients a[] using Chebyshev
61 : * polynomials. The found ISPs are in the cosine domain with values
62 : * in the range from 1 down to -1.
63 : * The table grid[] contains the points (in the cosine domain) at
64 : * which the polynomials are evaluated.
65 : *
66 : * The ISPs are the roots of the two polynomials F1(z) and F2(z)
67 : * defined as
68 : * F1(z) = [A(z) + z^-M A(z^-1)]
69 : * and F2(z) = [A(z) - z^-M A(z^-1)]/ (1-z^-2)
70 : *
71 : * For an even order M=2N, F1(z) has M/2 conjugate roots on the unit circle
72 : * and F2(z) has M/2-1 conjugate roots on the unit circle in addition to two
73 : * roots at 0 and pi.
74 : *
75 : * For a 16th order LP analysis (M=16), F1(z) and F2(z) can be written as
76 : *
77 : * F1(z) = (1 + a[16]) * PRODUCT (1 - 2 cos(w_i) z^-1 + z^-2 )
78 : * i=0,2,...,14
79 : *
80 : * F2(z) = (1 - a[16]) * PRODUCT (1 - 2 cos(w_i) z^-1 + z^-2 )
81 : * i=1,3,...,13
82 : *
83 : * The ISPs are frequencies w_i, i=0...M-2 plus the last predictor
84 : * coefficient a[M].
85 : *---------------------------------------------------------------------*/
86 :
87 151 : void a2isp(
88 : const float *a, /* i : LP filter coefficients */
89 : float *isp, /* o : Immittance spectral pairs */
90 : const float *old_isp /* i : ISP vector from past frame */
91 : )
92 : {
93 : int16_t j, i, nf, ip, order;
94 : float xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
95 : float f1[NC + 1], f2[NC];
96 : float *coef;
97 :
98 : /*-------------------------------------------------------------*
99 : * find the sum and diff polynomials F1(z) and F2(z)
100 : * F1(z) = [A(z) + z^M A(z^-1)]
101 : * F2(z) = [A(z) - z^M A(z^-1)]/(1-z^-2)
102 : *
103 : * for (i=0; i<NC; i++)
104 : * {
105 : * f1[i] = a[i] + a[M-i];
106 : * f2[i] = a[i] - a[M-i];
107 : * }
108 : * f1[NC] = 2.0*a[NC];
109 : *
110 : * for (i=2; i<NC; i++) Divide by (1-z^-2)
111 : * f2[i] += f2[i-2];
112 : *-------------------------------------------------------------*/
113 :
114 1359 : for ( i = 0; i < NC; i++ )
115 : {
116 1208 : f1[i] = a[i] + a[M - i];
117 1208 : f2[i] = a[i] - a[M - i];
118 : }
119 :
120 151 : f1[NC] = 2.0f * a[NC];
121 1057 : for ( i = 2; i < NC; i++ ) /* divide by (1 - z^-2) */
122 : {
123 906 : f2[i] += f2[i - 2];
124 : }
125 :
126 : /*-----------------------------------------------------------------*
127 : * Find the ISPs (roots of F1(z) and F2(z) ) using the
128 : * Chebyshev polynomial evaluation.
129 : * The roots of F1(z) and F2(z) are alternatively searched.
130 : * We start by finding the first root of F1(z) then we switch
131 : * to F2(z) then back to F1(z) and so on until all roots are found.
132 : *
133 : * - Evaluate Chebyshev pol. at grid points and check for sign change.
134 : * - If sign change track the root by subdividing the interval
135 : * 4 times and ckecking sign change.
136 : *-----------------------------------------------------------------*/
137 :
138 151 : nf = 0; /* number of found frequencies */
139 151 : ip = 0; /* flag to first polynomial */
140 :
141 151 : coef = f1; /* start with F1(z) */
142 151 : order = NC;
143 :
144 151 : xlow = grid100[0];
145 151 : ylow = chebps2( xlow, coef, order );
146 :
147 151 : j = 0;
148 :
149 16299 : while ( ( nf < M - 1 ) && ( j < GRID100_POINTS ) )
150 : {
151 16148 : j++;
152 16148 : xhigh = xlow;
153 16148 : yhigh = ylow;
154 16148 : xlow = grid100[j];
155 16148 : ylow = chebps2( xlow, coef, order );
156 16148 : if ( ylow * yhigh <= 0.0f ) /* if sign changes new root exists */
157 : {
158 2123 : j--;
159 :
160 : /* divide the interval of sign change by NO_ITER */
161 10615 : for ( i = 0; i < NO_ITER; i++ )
162 : {
163 8492 : xmid = (float) ( 0.5f * ( xlow + xhigh ) );
164 8492 : ymid = chebps2( xmid, coef, order );
165 :
166 8492 : if ( ylow * ymid <= 0.0f )
167 : {
168 4168 : yhigh = ymid;
169 4168 : xhigh = xmid;
170 : }
171 : else
172 : {
173 4324 : ylow = ymid;
174 4324 : xlow = xmid;
175 : }
176 : }
177 :
178 : /*--------------------------------------------------------*
179 : * Linear interpolation
180 : * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow)
181 : *--------------------------------------------------------*/
182 2123 : ymid = ( yhigh - ylow );
183 2123 : xint = xlow;
184 2123 : if ( ymid != 0 && ylow != 0 ) /* whenever ylow is 0, it doesn't make sense to compute the remaining part of the equation */
185 : {
186 2123 : xint -= ylow * ( xhigh - xlow ) / ( ymid );
187 : }
188 : #ifdef DEBUGGING
189 : else if ( ymid == 0 && ylow != 0 )
190 : {
191 : IVAS_ERROR( IVAS_ERR_INTERNAL, "issue in a2lsp_stab()" );
192 : }
193 : #endif
194 2123 : isp[nf] = xint; /* new root */
195 2123 : nf++;
196 :
197 2123 : ip = 1 - ip; /* flag to other polynomial */
198 2123 : coef = ip ? f2 : f1; /* pointer to other polynomial */
199 :
200 2123 : order = ip ? ( NC - 1 ) : NC; /* order of other polynomial */
201 :
202 2123 : xlow = xint;
203 2123 : ylow = chebps2( xlow, coef, order );
204 : }
205 : }
206 :
207 151 : isp[M - 1] = a[M];
208 :
209 : /*-----------------------------------------------------------------*
210 : * Check if m-1 roots found, if not use the ISPs from previous frame
211 : *-----------------------------------------------------------------*/
212 :
213 151 : if ( ( nf < M - 1 ) || ( a[M] > 1.0f ) || ( a[M] < -1.0f ) )
214 : {
215 646 : for ( i = 0; i < M; i++ )
216 : {
217 608 : isp[i] = old_isp[i];
218 : }
219 : }
220 :
221 151 : return;
222 : }
223 :
224 :
225 : /*-------------------------------------------------------------------*
226 : * isp2a()
227 : *
228 : * Convert ISPs to predictor coefficients a[]
229 : *-------------------------------------------------------------------*/
230 :
231 73560 : void isp2a(
232 : const float *isp, /* i : ISP vector (in the cosine domain) */
233 : float *a, /* o : LP filter coefficients */
234 : const int16_t m /* i : order of LP analysis */
235 : )
236 : {
237 : float f1[NC16k + 1], f2[NC16k];
238 : int16_t i, j;
239 : int16_t nc;
240 :
241 73560 : nc = m / 2;
242 :
243 : /*-----------------------------------------------------------------*
244 : * Find the polynomials F1(z) and F2(z) *
245 : *-----------------------------------------------------------------*/
246 :
247 73560 : get_isppol( &isp[0], f1, nc );
248 73560 : get_isppol( &isp[1], f2, (int16_t) ( nc - 1 ) );
249 :
250 : /*-----------------------------------------------------------------*
251 : * Multiply F2(z) by (1 - z^-2) *
252 : *-----------------------------------------------------------------*/
253 :
254 514920 : for ( i = nc - 1; i > 1; i-- )
255 : {
256 441360 : f2[i] -= f2[i - 2];
257 : }
258 :
259 : /*-----------------------------------------------------------------*
260 : * Scale F1(z) by (1+isp[m-1]) and F2(z) by (1-isp[m-1]) *
261 : *-----------------------------------------------------------------*/
262 :
263 662040 : for ( i = 0; i < nc; i++ )
264 : {
265 588480 : f1[i] *= ( 1.0f + isp[m - 1] );
266 588480 : f2[i] *= ( 1.0f - isp[m - 1] );
267 : }
268 :
269 : /*-----------------------------------------------------------------*
270 : * A(z) = (F1(z)+F2(z))/2 *
271 : * F1(z) is symmetric and F2(z) is asymmetric *
272 : *-----------------------------------------------------------------*/
273 :
274 73560 : a[0] = 1.0f;
275 588480 : for ( i = 1, j = m - 1; i < nc; i++, j-- )
276 : {
277 514920 : a[i] = (float) ( 0.5f * ( f1[i] + f2[i] ) );
278 514920 : a[j] = (float) ( 0.5f * ( f1[i] - f2[i] ) );
279 : }
280 73560 : a[nc] = (float) ( 0.5f * f1[nc] * ( 1.0f + isp[m - 1] ) );
281 73560 : a[m] = isp[m - 1];
282 :
283 73560 : return;
284 : }
285 :
286 : /*-------------------------------------------------------------------*
287 : * a2lsp()
288 : *
289 : * Convert predictor coefficients a[] to LSPs
290 : *-------------------------------------------------------------------*/
291 :
292 1213468 : int16_t a2lsp(
293 : float *freq, /* o : LSP vector */
294 : const float *a_in, /* i : predictor coefficients */
295 : const int16_t order /* i : order of LP analysis */
296 : )
297 : {
298 : int16_t i, iswitch, offset, STEPindex;
299 : int16_t lspnumber, root, notlast, order_by_2;
300 : float temp, temp2;
301 : float q[20], prev[2];
302 : float frequency, LastFreq, STEP;
303 : const float *a;
304 : float space_min;
305 :
306 1213468 : a = &( a_in[1] );
307 1213468 : order_by_2 = order / 2;
308 1213468 : LastFreq = 0;
309 :
310 : /* calculate q[z] and p[z] , they are all stored in q */
311 1213468 : offset = order_by_2;
312 1213468 : q[0] = (float) ( a[0] + a[order - 1] - 1.0 );
313 1213468 : q[offset] = (float) ( a[0] - a[order - 1] + 1.0 );
314 5753451 : for ( i = 1; i < order_by_2; i++ )
315 : {
316 4539983 : q[i] = a[i] + a[order - 1 - i] - q[i - 1];
317 4539983 : q[i + offset] = a[i] - a[order - 1 - i] + q[i - 1 + offset];
318 : }
319 :
320 1213468 : q[order_by_2 - 1] = q[order_by_2 - 1] / 2;
321 1213468 : q[order_by_2 - 1 + offset] = q[order_by_2 - 1 + offset] / 2;
322 :
323 1213468 : prev[0] = 9e9f;
324 1213468 : prev[1] = 9e9f;
325 1213468 : lspnumber = 0;
326 1213468 : notlast = 1;
327 1213468 : iswitch = 0;
328 1213468 : frequency = 0;
329 :
330 12720370 : while ( notlast )
331 : {
332 11506902 : root = 1;
333 11506902 : offset = iswitch * order_by_2;
334 11506902 : STEPindex = 0; /* Start with low resolution grid */
335 11506902 : STEP = STEPS[STEPindex];
336 161945301 : while ( root )
337 : {
338 150438399 : temp = (float) cos( frequency * 6.2832 );
339 150438399 : if ( order >= 4 )
340 : {
341 150438399 : temp2 = LPC_chebyshev( temp, q + offset, order_by_2 );
342 : }
343 : else
344 : {
345 0 : temp2 = temp + q[0 + offset];
346 : }
347 :
348 150438399 : if ( ( temp2 * prev[iswitch] ) <= 0.0 || frequency >= 0.5 )
349 : {
350 46027608 : if ( STEPindex == STEPSNUM - 1 )
351 : {
352 11506902 : if ( fabs( temp2 ) < fabs( prev[iswitch] ) )
353 : {
354 5751991 : freq[lspnumber] = frequency;
355 : }
356 : else
357 : {
358 5754911 : freq[lspnumber] = frequency - STEP;
359 : }
360 11506902 : if ( ( prev[iswitch] ) < 0.0 )
361 : {
362 4553560 : prev[iswitch] = 9e9f;
363 : }
364 : else
365 : {
366 6953342 : prev[iswitch] = -9e9f;
367 : }
368 11506902 : root = 0;
369 11506902 : frequency = LastFreq;
370 11506902 : STEPindex = 0;
371 : }
372 : else
373 : {
374 34520706 : if ( STEPindex == 0 )
375 : {
376 11506902 : LastFreq = frequency;
377 : }
378 34520706 : frequency -= STEPS[++STEPindex]; /* Go back one grid step */
379 34520706 : STEP = STEPS[STEPindex];
380 : }
381 : }
382 : else
383 : {
384 104410791 : prev[iswitch] = temp2;
385 104410791 : frequency += STEP;
386 : }
387 : } /* while(root) */
388 :
389 11506902 : lspnumber++;
390 11506902 : if ( lspnumber > order - 1 )
391 : {
392 1213468 : notlast = 0;
393 : }
394 11506902 : iswitch = 1 - iswitch;
395 : } /* while (notlast) */
396 :
397 : /* stability check */
398 1213468 : space_min = 1;
399 11506902 : for ( i = 1; i < order; i++ )
400 : {
401 10293434 : space_min = ( ( freq[i] - freq[i - 1] ) < space_min ) ? ( freq[i] - freq[i - 1] ) : space_min;
402 : }
403 :
404 1213468 : if ( space_min <= 0 )
405 : {
406 94 : return 0;
407 : }
408 :
409 1213374 : return 1;
410 : }
411 :
412 : /*-------------------------------------------------------------------*
413 : * lsp2a()
414 : *
415 : * Convert LSPs to predictor coefficients a[]
416 : *-------------------------------------------------------------------*/
417 :
418 6180204 : void lsp2a(
419 : float *pc_in, /* i/o: predictor coefficients */
420 : float *freq, /* i/o: LSP coefficients */
421 : const int16_t order /* i : order of LP analysis */
422 : )
423 : {
424 : float p[LPC_SHB_ORDER], q[LPC_SHB_ORDER];
425 : float a[LPC_SHB_ORDER], a1[LPC_SHB_ORDER], a2[LPC_SHB_ORDER];
426 : float b[LPC_SHB_ORDER], b1[LPC_SHB_ORDER], b2[LPC_SHB_ORDER];
427 :
428 : float xx;
429 : int16_t i, k;
430 : int16_t lspflag;
431 : float *pc;
432 : int16_t order_by_2;
433 :
434 6180204 : lspflag = 1;
435 6180204 : pc = &( pc_in[1] );
436 :
437 6180204 : order_by_2 = order / 2;
438 :
439 : /* check input for ill-conditioned cases */
440 6180204 : if ( ( freq[0] <= 0.0 ) || ( freq[order - 1] >= 0.5 ) )
441 : {
442 0 : lspflag = 0;
443 :
444 0 : if ( freq[0] <= 0 )
445 : {
446 0 : freq[0] = 0.022f;
447 : }
448 :
449 0 : if ( freq[order - 1] >= 0.5 )
450 : {
451 0 : freq[order - 1] = 0.499f;
452 : }
453 : }
454 :
455 6180204 : if ( !lspflag )
456 : {
457 0 : xx = ( freq[order - 1] - freq[0] ) * recip_order[order];
458 0 : for ( i = 1; i < order; i++ )
459 : {
460 0 : freq[i] = freq[i - 1] + xx;
461 : }
462 : }
463 :
464 42364847 : for ( i = 0; i <= order_by_2; i++ )
465 : {
466 36184643 : a[i] = 0.;
467 36184643 : a1[i] = 0.;
468 36184643 : a2[i] = 0.;
469 36184643 : b[i] = 0.;
470 36184643 : b1[i] = 0.;
471 36184643 : b2[i] = 0.;
472 : }
473 :
474 36184643 : for ( i = 0; i < order_by_2; i++ )
475 : {
476 30004439 : p[i] = (float) cos( 6.2832 * freq[2 * i] );
477 30004439 : q[i] = (float) cos( 6.2832 * freq[2 * i + 1] );
478 : }
479 :
480 6180204 : a[0] = 0.25f;
481 6180204 : b[0] = 0.25f;
482 :
483 36184643 : for ( i = 0; i < order_by_2; i++ )
484 : {
485 30004439 : a[i + 1] = a[i] - 2 * p[i] * a1[i] + a2[i];
486 30004439 : b[i + 1] = b[i] - 2 * q[i] * b1[i] + b2[i];
487 30004439 : a2[i] = a1[i];
488 30004439 : a1[i] = a[i];
489 30004439 : b2[i] = b1[i];
490 30004439 : b1[i] = b[i];
491 : }
492 6180204 : a[0] = 0.25f;
493 6180204 : b[0] = -0.25f;
494 :
495 36184643 : for ( i = 0; i < order_by_2; i++ )
496 : {
497 30004439 : a[i + 1] = a[i] - 2 * p[i] * a1[i] + a2[i];
498 30004439 : b[i + 1] = b[i] - 2 * q[i] * b1[i] + b2[i];
499 30004439 : a2[i] = a1[i];
500 30004439 : a1[i] = a[i];
501 30004439 : b2[i] = b1[i];
502 30004439 : b1[i] = b[i];
503 : }
504 :
505 6180204 : pc[0] = 2 * ( a[order_by_2] + b[order_by_2] );
506 :
507 60008878 : for ( k = 2; k <= order; k++ )
508 : {
509 53828674 : a[0] = 0.0f;
510 53828674 : b[0] = 0.0f;
511 :
512 319002193 : for ( i = 0; i < order_by_2; i++ )
513 : {
514 265173519 : a[i + 1] = a[i] - 2 * p[i] * a1[i] + a2[i];
515 265173519 : b[i + 1] = b[i] - 2 * q[i] * b1[i] + b2[i];
516 265173519 : a2[i] = a1[i];
517 265173519 : a1[i] = a[i];
518 265173519 : b2[i] = b1[i];
519 265173519 : b1[i] = b[i];
520 : }
521 53828674 : pc[k - 1] = 2 * ( a[order_by_2] + b[order_by_2] );
522 : }
523 :
524 6180204 : return;
525 : }
526 :
527 : /*-----------------------------------------------------------*
528 : * get_lsppol()
529 : *
530 : * Find the polynomial F1(z) or F2(z) from the LSPs.
531 : * This is performed by expanding the product polynomials:
532 : *
533 : * F1(z) = product ( 1 - 2 LSF_i z^-1 + z^-2 )
534 : * i=0,2,4,..,n-2
535 : * F2(z) = product ( 1 - 2 LSF_i z^-1 + z^-2 )
536 : * i=1,3,5,..,n-1
537 : *
538 : * where LSP_i are the LSPs in the cosine domain.
539 : *-----------------------------------------------------------*/
540 :
541 278741106 : static void get_lsppol(
542 : const float lsp[], /* i : line spectral freq. (cosine domain) */
543 : float f[], /* o : the coefficients of F1 or F2 */
544 : const int16_t n, /* i : no of coefficients (m/2) */
545 : int16_t flag /* i : 1 ---> F1(z); 2 ---> F2(z) */
546 : )
547 : {
548 : float b;
549 : const float *plsp;
550 : int16_t i, j;
551 :
552 278741106 : plsp = lsp + flag - 1;
553 :
554 278741106 : f[0] = 1.0f;
555 278741106 : b = -2.0f * *plsp;
556 278741106 : f[1] = b;
557 :
558 2229928848 : for ( i = 2; i <= n; i++ )
559 : {
560 1951187742 : plsp += 2;
561 1951187742 : b = -2.0f * *plsp;
562 1951187742 : f[i] = b * f[i - 1] + 2.0f * f[i - 2];
563 :
564 7804750968 : for ( j = i - 1; j > 1; j-- )
565 : {
566 5853563226 : f[j] += b * f[j - 1] + f[j - 2];
567 : }
568 :
569 1951187742 : f[1] += b;
570 : }
571 :
572 278741106 : return;
573 : }
574 :
575 : /*---------------------------------------------------------------------*
576 : * a2lsp_stab()
577 : *
578 : * Compute the LSPs from the LPC coefficients a[] using Chebyshev
579 : * polynomials. The found LSPs are in the cosine domain with values
580 : * in the range from 1 down to -1.
581 : * The table grid[] contains the points (in the cosine domain) at
582 : * which the polynomials are evaluated.
583 : *
584 : * The ISPs are the roots of the two polynomials F1(z) and F2(z)
585 : * defined as
586 : * F1(z) = [A(z) + z^-M A(z^-1)]/ (1+z^-1)
587 : * and F2(z) = [A(z) - z^-M A(z^-1)]/ (1-z^-1)
588 : *---------------------------------------------------------------------*/
589 :
590 39176438 : void a2lsp_stab(
591 : const float *a, /* i : LP filter coefficients */
592 : float *lsp, /* o : LSP vector */
593 : const float *old_lsp /* i : LSP vector from past frame */
594 : )
595 : {
596 : int16_t j, i, nf, ip;
597 : float xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
598 : float *pf1, *pf2;
599 : const float *pa1, *pa2;
600 : float f1[NC + 1], f2[NC + 1];
601 :
602 : /*-------------------------------------------------------------*
603 : * find the sum and diff polynomials F1(z) and F2(z) *
604 : * F1(z) = [A(z) + z^11 A(z^-1)]/(1+z^-1) *
605 : * F2(z) = [A(z) - z^11 A(z^-1)]/(1-z^-1) *
606 : *-------------------------------------------------------------*/
607 :
608 39176438 : pf1 = f1; /* Equivalent code using indices */
609 39176438 : pf2 = f2;
610 39176438 : *pf1++ = 1.0f; /* f1[0] = 1.0; */
611 39176438 : *pf2++ = 1.0f; /* f2[0] = 1.0; */
612 39176438 : pa1 = a + 1;
613 39176438 : pa2 = a + M;
614 :
615 352587942 : for ( i = 0; i <= NC - 1; i++ ) /* for (i=1, j=M; i<=NC; i++, j--) */
616 : {
617 313411504 : *pf1 = *pa1 + *pa2 - *( pf1 - 1 ); /* f1[i] = a[i]+a[j]-f1[i-1]; */
618 313411504 : *pf2 = *pa1++ - *pa2-- + *( pf2 - 1 ); /* f2[i] = a[i]-a[j]+f2[i-1]; */
619 313411504 : pf1++;
620 313411504 : pf2++;
621 : }
622 :
623 : /*---------------------------------------------------------------------*
624 : * Find the LSPs (roots of F1(z) and F2(z) ) using the *
625 : * Chebyshev polynomial evaluation. *
626 : * The roots of F1(z) and F2(z) are alternatively searched. *
627 : * We start by finding the first root of F1(z) then we switch *
628 : * to F2(z) then back to F1(z) and so on until all roots are found. *
629 : * *
630 : * - Evaluate Chebyshev pol. at grid points and check for sign change.*
631 : * - If sign change track the root by subdividing the interval *
632 : * 4 times and ckecking sign change. *
633 : *---------------------------------------------------------------------*/
634 :
635 39176438 : nf = 0; /* number of found frequencies */
636 39176438 : ip = 0; /* flag to first polynomial */
637 :
638 39176438 : pf1 = f1; /* start with F1(z) */
639 :
640 39176438 : xlow = grid100[0];
641 39176438 : ylow = chebps2( xlow, pf1, NC );
642 :
643 39176438 : j = 0;
644 4283945180 : while ( ( nf < M ) && ( j < GRID100_POINTS ) )
645 : {
646 4244768742 : j++;
647 4244768742 : xhigh = xlow;
648 4244768742 : yhigh = ylow;
649 4244768742 : xlow = grid100[j];
650 4244768742 : ylow = chebps2( xlow, pf1, NC );
651 :
652 4244768742 : if ( ylow * yhigh <= 0.0 ) /* if sign change new root exists */
653 : {
654 626816364 : j--;
655 : /* divide the interval of sign change by NO_ITER */
656 3134081820 : for ( i = 0; i < NO_ITER; i++ )
657 : {
658 2507265456 : xmid = 0.5f * ( xlow + xhigh );
659 2507265456 : ymid = chebps2( xmid, pf1, NC );
660 2507265456 : if ( ylow * ymid <= 0.0 )
661 : {
662 1251172553 : yhigh = ymid;
663 1251172553 : xhigh = xmid;
664 : }
665 : else
666 : {
667 1256092903 : ylow = ymid;
668 1256092903 : xlow = xmid;
669 : }
670 : }
671 :
672 : /* linear interpolation for evaluating the root */
673 626816364 : ymid = ( yhigh - ylow );
674 626816364 : xint = xlow;
675 626816364 : if ( ymid != 0 && ylow != 0 ) /* whenever ylow is 0, it doesn't make sense to compute the remaining part of the equation */
676 : {
677 626815646 : xint -= ylow * ( xhigh - xlow ) / ( ymid );
678 : }
679 : #ifdef DEBUGGING
680 : else if ( ymid == 0 && ylow != 0 )
681 : {
682 : IVAS_ERROR( IVAS_ERR_INTERNAL, "issue in a2lsp_stab()" );
683 : }
684 : #endif
685 626816364 : lsp[nf] = xint; /* new root */
686 626816364 : nf++;
687 626816364 : ip = 1 - ip; /* flag to other polynomial */
688 626816364 : pf1 = ip ? f2 : f1; /* pointer to other polynomial */
689 626816364 : xlow = xint;
690 626816364 : ylow = chebps2( xlow, pf1, NC );
691 : }
692 : }
693 :
694 : /* Check if M roots found */
695 : /* if not use the LSPs from previous frame */
696 39176438 : if ( nf < M )
697 : {
698 :
699 28237 : for ( i = 0; i < M; i++ )
700 : {
701 26576 : lsp[i] = old_lsp[i];
702 : }
703 : }
704 :
705 39176438 : return;
706 : }
707 :
708 :
709 : /*-------------------------------------------------------------------*
710 : * lsp2a_stab()
711 : *
712 : * Convert LSPs to predictor coefficients A[]
713 : *-------------------------------------------------------------------*/
714 :
715 139370553 : void lsp2a_stab(
716 : const float *lsp, /* i : LSF vector (in the cosine domain) */
717 : float *a, /* o : LP filter coefficients */
718 : const int16_t m /* i : order of LP analysis */
719 : )
720 : {
721 : float f1[NC + 1], f2[NC + 1];
722 : int16_t i, k, nc;
723 : float *pf1, *pf2, *pf1_1, *pf2_1, *pa1, *pa2;
724 :
725 139370553 : nc = m / 2;
726 :
727 : /*-----------------------------------------------------*
728 : * Find the polynomials F1(z) and F2(z) *
729 : *-----------------------------------------------------*/
730 :
731 139370553 : get_lsppol( lsp, f1, nc, 1 );
732 139370553 : get_lsppol( lsp, f2, nc, 2 );
733 :
734 : /*-----------------------------------------------------*
735 : * Multiply F1(z) by (1+z^-1) and F2(z) by (1-z^-1) *
736 : *-----------------------------------------------------*/
737 :
738 139370553 : pf1 = f1 + nc;
739 139370553 : pf1_1 = pf1 - 1;
740 139370553 : pf2 = f2 + nc; /* Version using indices */
741 139370553 : pf2_1 = pf2 - 1;
742 139370553 : k = nc - 1;
743 1254334977 : for ( i = 0; i <= k; i++ ) /* for (i = NC; i > 0; i--) */
744 : {
745 1114964424 : *pf1-- += *pf1_1--; /* f1[i] += f1[i-1]; */
746 1114964424 : *pf2-- -= *pf2_1--; /* f2[i] -= f2[i-1]; */
747 : }
748 :
749 : /*-----------------------------------------------------*
750 : * A(z) = (F1(z)+F2(z))/2 *
751 : * F1(z) is symmetric and F2(z) is antisymmetric *
752 : *-----------------------------------------------------*/
753 :
754 139370553 : pa1 = a;
755 139370553 : *pa1++ = 1.0; /* a[0] = 1.0; */
756 139370553 : pa2 = a + m;
757 139370553 : pf1 = f1 + 1;
758 139370553 : pf2 = f2 + 1;
759 1254334977 : for ( i = 0; i <= k; i++ ) /* for (i=1, j=M; i<=NC; i++, j--) */
760 : {
761 1114964424 : *pa1++ = 0.5f * ( *pf1 + *pf2 ); /* a[i] = 0.5*(f1[i] + f2[i]); */
762 1114964424 : *pa2-- = 0.5f * ( *pf1++ - *pf2++ ); /* a[j] = 0.5*(f1[i] - f2[i]); */
763 : }
764 :
765 139370553 : return;
766 : }
767 :
768 : /*---------------------------------------------------------------------------
769 : * reorder_lsf()
770 : *
771 : * To make sure that the LSFs are properly ordered and to keep a certain
772 : * minimum distance between consecutive LSFs.
773 : *--------------------------------------------------------------------------*/
774 :
775 23040389 : void reorder_lsf(
776 : float *lsf, /* i/o: LSF vector */
777 : const float min_dist, /* i : minimum required distance */
778 : const int16_t n, /* i : LPC order */
779 : const int32_t Fs /* i : sampling frequency */
780 : )
781 : {
782 : int16_t i;
783 : float lsf_min;
784 : float lsf_max;
785 :
786 : /*-----------------------------------------------------------------*
787 : * Verify the LSF ordering and minimum GAP
788 : *-----------------------------------------------------------------*/
789 :
790 23040389 : lsf_min = min_dist;
791 391686613 : for ( i = 0; i < n; i++ )
792 : {
793 368646224 : if ( lsf[i] < lsf_min )
794 : {
795 1519048 : lsf[i] = lsf_min;
796 : }
797 368646224 : lsf_min = lsf[i] + min_dist;
798 : }
799 :
800 : /*------------------------------------------------------------------------------------------*
801 : * Reverify the LSF ordering and minimum GAP in the reverse order (security)
802 : *------------------------------------------------------------------------------------------*/
803 :
804 23040389 : lsf_max = Fs / 2.0f - min_dist;
805 :
806 23040389 : if ( lsf[n - 1] > lsf_max ) /* If danger of unstable filter in case of resonance in HF */
807 : {
808 95846 : for ( i = n - 1; i >= 0; i-- ) /* Reverify the minimum LSF gap in the reverse sens */
809 : {
810 90208 : if ( lsf[i] > lsf_max )
811 : {
812 17932 : lsf[i] = lsf_max;
813 : }
814 :
815 90208 : lsf_max = lsf[i] - min_dist;
816 : }
817 : }
818 :
819 23040389 : return;
820 : }
821 :
822 : /*-------------------------------------------------------------------*
823 : * space_lsfs()
824 : *
825 : *-------------------------------------------------------------------*/
826 :
827 2896112 : void space_lsfs(
828 : float *lsfs, /* i/o: Line spectral frequencies */
829 : const int16_t order /* i : order of LP analysis */
830 : )
831 : {
832 : float delta;
833 2896112 : int16_t i, flag = 1;
834 :
835 6487807 : while ( flag == 1 )
836 : {
837 3591695 : flag = 0;
838 43100340 : for ( i = 0; i <= order; i++ )
839 : {
840 39508645 : delta = (float) ( i == 0 ? lsfs[0] : ( i == order ? 0.5f - lsfs[i - 1] : ( lsfs[i] - lsfs[i - 1] ) ) );
841 39508645 : if ( delta < SPC )
842 : {
843 1030414 : flag = 1;
844 1030414 : delta -= SPC_plus;
845 1030414 : if ( i == order )
846 : {
847 68355 : lsfs[i - 1] += delta;
848 : }
849 : else
850 : {
851 962059 : if ( i == 0 )
852 : {
853 0 : lsfs[i] -= delta;
854 : }
855 : else
856 : {
857 962059 : delta *= 0.5f;
858 962059 : lsfs[i - 1] += delta;
859 962059 : lsfs[i] -= delta;
860 : }
861 : }
862 : }
863 : }
864 : }
865 :
866 2896112 : return;
867 : }
868 :
869 : /*-------------------------------------------------------------------*
870 : * lsp_weights()
871 : *
872 : *-------------------------------------------------------------------*/
873 :
874 353146 : void lsp_weights(
875 : const float *lsps, /* i : Line spectral pairs */
876 : float *weight, /* o : weights */
877 : const int16_t order /* i : order of LP analysis */
878 : )
879 : {
880 : int16_t i;
881 : float delta1, delta2;
882 :
883 3557140 : for ( i = 0; i < order; i++ )
884 : {
885 3203994 : delta1 = (float) ( ( i == 0 ) ? lsps[i] : lsps[i] - lsps[i - 1] );
886 3203994 : delta2 = (float) ( ( i == order - 1 ) ? 0.5f - lsps[i] : lsps[i + 1] - lsps[i] );
887 :
888 3203994 : weight[i] = 250 * root_a_over_b( ALPHA_SQ, delta1 * delta2 );
889 : }
890 :
891 353146 : if ( order != LPC_SHB_ORDER_WB )
892 : {
893 291645 : weight[3] *= 1.1f;
894 291645 : weight[4] *= 1.1f;
895 : }
896 :
897 353146 : return;
898 : }
899 :
900 :
901 : /*-----------------------------------------------------------------------*
902 : * isf2isp()
903 : *
904 : * Transformation of ISF to ISP
905 : *
906 : * ISP are immitance spectral pairs in cosine domain (-1 to 1).
907 : * ISF are immitance spectral pairs in frequency domain (0 to fs/2).
908 : *-----------------------------------------------------------------------*/
909 :
910 16975 : void isf2isp(
911 : const float isf[], /* i : isf[m] normalized (range: 0<=val<=fs/2) */
912 : float isp[], /* o : isp[m] (range: -1<=val<1) */
913 : const int16_t m, /* i : LPC order */
914 : const int32_t Fs /* i : sampling frequency */
915 : )
916 : {
917 : int16_t i;
918 :
919 271600 : for ( i = 0; i < m - 1; i++ )
920 : {
921 254625 : isp[i] = (float) cos( isf[i] * EVS_PI / ( Fs / 2.f ) );
922 : }
923 16975 : isp[m - 1] = (float) cos( isf[m - 1] * EVS_PI / ( Fs / 2.f ) * 2.0 );
924 :
925 16975 : return;
926 : }
927 :
928 :
929 : /*-----------------------------------------------------------------------*
930 : * isp2isf()
931 : *
932 : * Transformation of ISP to ISF
933 : *
934 : * ISP are immitance spectral pair in cosine domain (-1 to 1).
935 : * ISF are immitance spectral pair in frequency domain (0 to fs/2).
936 : *-----------------------------------------------------------------------*/
937 :
938 65 : void isp2isf(
939 : const float isp[], /* i : isp[m] (range: -1<=val<1) */
940 : float isf[], /* o : isf[m] normalized (range: 0<=val<=fs/2) */
941 : const int16_t m, /* i : LPC order */
942 : const int32_t Fs /* i : sampling frequency */
943 : )
944 : {
945 : int16_t i;
946 :
947 : /* convert ISPs to frequency domain 0..6400 */
948 1040 : for ( i = 0; i < m - 1; i++ )
949 : {
950 975 : isf[i] = (float) ( acos( isp[i] ) * ( ( Fs / 2.f ) / EVS_PI ) );
951 : }
952 65 : isf[m - 1] = (float) ( acos( isp[m - 1] ) * ( ( Fs / 2.f ) / EVS_PI ) * 0.5f );
953 :
954 65 : return;
955 : }
956 :
957 : /*-------------------------------------------------------------------*
958 : * a2rc()
959 : *
960 : * Convert from LPC to reflection coeff
961 : *-------------------------------------------------------------------*/
962 : /*! r: stability flag */
963 3448649 : uint16_t a2rc(
964 : const float *a, /* i : LPC coefficients */
965 : float *refl, /* o : Reflection co-efficients */
966 : const int16_t lpcorder /* i : LPC order */
967 : )
968 : {
969 : float f[MAX_LP_FILTER_ORDER];
970 : int16_t m, j, n;
971 : float km, denom, x;
972 :
973 59169629 : for ( m = 0; m < lpcorder; m++ )
974 : {
975 55720980 : f[m] = -a[m];
976 : }
977 :
978 : /* Initialization */
979 58955490 : for ( m = lpcorder - 1; m >= 0; m-- )
980 : {
981 55575392 : km = f[m];
982 55575392 : if ( km <= -1.0 || km >= 1.0 )
983 : {
984 1439571 : for ( j = 0; j < lpcorder; j++ )
985 : {
986 1371020 : refl[j] = 0.f;
987 : }
988 68551 : return 0;
989 : }
990 :
991 55506841 : refl[m] = -km;
992 55506841 : denom = 1.0f / ( 1.0f - km * km );
993 :
994 253097526 : for ( j = 0; j < m / 2; j++ )
995 : {
996 197590685 : n = m - 1 - j;
997 197590685 : x = denom * f[j] + km * denom * f[n];
998 197590685 : f[n] = denom * f[n] + km * denom * f[j];
999 197590685 : f[j] = x;
1000 : }
1001 :
1002 55506841 : if ( m & 1 )
1003 : {
1004 27757352 : f[j] = denom * f[j] + km * denom * f[j];
1005 : }
1006 : }
1007 :
1008 3380098 : return 1;
1009 : }
1010 :
1011 :
1012 : /*----------------------------------------------------------------------------------*
1013 : * vq_dec_lvq()
1014 : *
1015 : * Multi-stage VQ decoder for LSF parameters, last stage LVQ
1016 : *----------------------------------------------------------------------------------*/
1017 :
1018 9766424 : int16_t vq_dec_lvq(
1019 : int16_t sf_flag, /* i : safety net flag */
1020 : float x[], /* o : Decoded vector */
1021 : int16_t indices[], /* i : Indices */
1022 : const int16_t stages, /* i : Number of stages */
1023 : const int16_t N, /* i : Vector dimension */
1024 : const int16_t mode, /* i : mode_lvq, or mode_lvq_p */
1025 : const int16_t no_bits /* i : no. bits for lattice */
1026 : )
1027 : {
1028 : float x_lvq[16];
1029 : int16_t i;
1030 : int16_t ber_flag;
1031 :
1032 : /* clear vector */
1033 9766424 : set_f( x, 0, N );
1034 :
1035 : /*-----------------------------------------------*
1036 : * add contribution of each stage
1037 : *-----------------------------------------------*/
1038 :
1039 9766424 : if ( sf_flag == 1 )
1040 : {
1041 2911009 : for ( i = 0; i < stages - 1; i++ )
1042 : {
1043 1803400 : v_add( x, &Quantizers[CB_lsf[mode] + i][indices[i] * N], x, N );
1044 : }
1045 :
1046 1107609 : ber_flag = deindex_lvq( &indices[stages - 1], x_lvq, mode, sf_flag, no_bits );
1047 : }
1048 : else
1049 : {
1050 17267501 : for ( i = 0; i < stages - 1; i++ )
1051 : {
1052 8608686 : v_add( x, &Quantizers_p[CB_p_lsf[mode] + i][indices[i] * N], x, N );
1053 : }
1054 :
1055 8658815 : ber_flag = deindex_lvq( &indices[stages - 1], x_lvq, mode, sf_flag, no_bits );
1056 : }
1057 :
1058 9766424 : v_add( x, x_lvq, x, N );
1059 :
1060 9766424 : return ber_flag;
1061 : }
1062 :
1063 :
1064 : /*----------------------------------------------------------------------------------*
1065 : * lsf_allocate()
1066 : *
1067 : * Calculate number of stages and levels for each stage based on the allowed bit allocation
1068 : *----------------------------------------------------------------------------------*/
1069 :
1070 13263659 : ivas_error lsf_allocate(
1071 : const int16_t nBits, /* i : Number of bits to use for quantization */
1072 : const int16_t framemode, /* i : LSF quantizer mode */
1073 : const int16_t framemode_p, /* i : ISF quantizer mode predmode (mode_lvq_p) */
1074 : int16_t *stages0, /* o : Number of stages for safety-net quantizer */
1075 : int16_t *stages1, /* o : Number of stages for predictive quantizer */
1076 : int16_t levels0[], /* o : Number of vectors for each stage for SFNET */
1077 : int16_t levels1[], /* o : Number of vectors for each stage for pred */
1078 : int16_t bits0[], /* o : Number of bits for each stage safety net */
1079 : int16_t bits1[] /* o : Number of bits for each stage pred */
1080 : )
1081 : {
1082 : int16_t i;
1083 : int16_t cumleft;
1084 : int16_t bits_lvq, n_stages, nbits0;
1085 : ivas_error error;
1086 :
1087 13263659 : error = IVAS_ERR_OK;
1088 :
1089 : /*---------------------------------------------------*
1090 : * Calculate bit allocation for safety-net quantizer
1091 : *---------------------------------------------------*/
1092 :
1093 13263659 : cumleft = BitsVQ[framemode];
1094 13263659 : bits_lvq = nBits - cumleft;
1095 13263659 : nbits0 = CBbits[framemode];
1096 13263659 : if ( nbits0 > -1 )
1097 : {
1098 4817288 : if ( nbits0 > 0 ) /* */
1099 : {
1100 4817288 : n_stages = 2;
1101 4817288 : levels0[0] = CBsizes[nbits0];
1102 4817288 : bits0[0] = nbits0;
1103 4817288 : bits0[1] = cumleft - nbits0;
1104 :
1105 4817288 : if ( bits0[1] == 0 )
1106 : {
1107 2082045 : n_stages--;
1108 : }
1109 : else
1110 : {
1111 2735243 : levels0[1] = CBsizes[cumleft - nbits0];
1112 : }
1113 : }
1114 : else /* no bits for VQ stage */
1115 : {
1116 0 : n_stages = 0;
1117 : }
1118 :
1119 4817288 : *stages0 = n_stages;
1120 4817288 : if ( bits_lvq > 0 )
1121 : {
1122 4817288 : bits0[n_stages] = bits_lvq;
1123 4817288 : levels0[n_stages] = bits_lvq; /* this is number of bits, not levels */
1124 4817288 : *stages0 = n_stages + 1;
1125 : }
1126 : }
1127 : else
1128 : {
1129 8446371 : *stages0 = 0;
1130 : }
1131 :
1132 : /*---------------------------------------------------*
1133 : * Calculate bit allocation for predictive quantizer
1134 : *---------------------------------------------------*/
1135 :
1136 13263659 : if ( framemode_p > -1 )
1137 : {
1138 12819595 : cumleft = BitsVQ_p[framemode_p];
1139 12819595 : bits_lvq = nBits - cumleft;
1140 12819595 : nbits0 = CBbits_p[framemode_p];
1141 :
1142 12819595 : if ( nbits0 > -1 )
1143 : {
1144 12819595 : if ( nbits0 > 0 )
1145 : {
1146 10737550 : if ( framemode_p == 7 )
1147 : {
1148 : /* for UNVOICED_WB only */
1149 53094 : n_stages = 3;
1150 212376 : for ( i = 0; i < n_stages; i++ )
1151 : {
1152 159282 : levels1[i] = CBsizes[nbits0];
1153 159282 : bits1[i] = nbits0;
1154 : }
1155 53094 : bits1[n_stages] = bits_lvq;
1156 53094 : levels1[n_stages] = bits_lvq;
1157 53094 : *stages1 = n_stages + 1;
1158 : }
1159 : else
1160 : {
1161 10684456 : n_stages = 1;
1162 10684456 : levels1[0] = CBsizes[nbits0];
1163 10684456 : bits1[0] = nbits0;
1164 10684456 : nbits0 = cumleft - nbits0;
1165 10684456 : if ( nbits0 > 0 )
1166 : {
1167 985979 : levels1[1] = CBsizes[nbits0];
1168 985979 : bits1[1] = nbits0;
1169 985979 : n_stages = 2;
1170 : }
1171 :
1172 10684456 : levels1[n_stages] = bits_lvq; /* this is number of bits, not levels */
1173 10684456 : bits1[n_stages] = bits_lvq;
1174 10684456 : *stages1 = n_stages + 1;
1175 : }
1176 : }
1177 : else
1178 : {
1179 2082045 : *stages1 = 1;
1180 2082045 : bits1[0] = bits_lvq;
1181 2082045 : levels1[0] = bits_lvq;
1182 : }
1183 : }
1184 : #ifdef DEBUGGING
1185 : else
1186 : {
1187 : return IVAS_ERROR( IVAS_ERR_INTERNAL_FATAL, "lsf_allocate(): invalid number of bits in used predictive mode\n" );
1188 : }
1189 : #endif
1190 : }
1191 :
1192 13263659 : return error;
1193 : }
1194 :
1195 : /*----------------------------------------------------------------------------------*
1196 : * find_pred_mode()
1197 : *
1198 : *
1199 : *----------------------------------------------------------------------------------*/
1200 :
1201 13257759 : ivas_error find_pred_mode(
1202 : int16_t *predmode,
1203 : const int16_t coder_type,
1204 : const int16_t bwidth,
1205 : const int32_t int_fs,
1206 : int16_t *p_mode_lvq,
1207 : int16_t *p_mode_lvq_p,
1208 : const int32_t core_brate )
1209 : {
1210 : int16_t idx;
1211 : ivas_error error;
1212 :
1213 13257759 : error = IVAS_ERR_OK;
1214 :
1215 13257759 : idx = bwidth;
1216 13257759 : if ( idx > 1 )
1217 : {
1218 4545415 : idx = 1;
1219 : }
1220 :
1221 13257759 : if ( int_fs == INT_FS_16k )
1222 : {
1223 8617953 : idx = 2;
1224 : }
1225 : else
1226 : {
1227 4639806 : if ( core_brate >= GENERIC_MA_LIMIT && coder_type == GENERIC && idx == 1 )
1228 : {
1229 1308924 : idx = 3;
1230 : }
1231 : }
1232 :
1233 13257759 : *predmode = predmode_tab[idx][coder_type];
1234 :
1235 13257759 : if ( idx <= 2 )
1236 : {
1237 11948835 : *p_mode_lvq = NO_CODING_MODES * idx + coder_type;
1238 11948835 : if ( *predmode > 0 )
1239 : {
1240 11504771 : *p_mode_lvq_p = *p_mode_lvq;
1241 : }
1242 : else
1243 : {
1244 444064 : *p_mode_lvq_p = -1;
1245 : }
1246 : }
1247 : else /* WB 12.8 with MA pred in GENERIC*/
1248 : {
1249 1308924 : *p_mode_lvq = NO_CODING_MODES + coder_type;
1250 1308924 : if ( coder_type == GENERIC )
1251 : {
1252 1308924 : *p_mode_lvq_p = 18;
1253 : }
1254 : else
1255 : {
1256 0 : if ( *predmode > 0 )
1257 : {
1258 0 : *p_mode_lvq_p = *p_mode_lvq;
1259 : }
1260 : else
1261 : {
1262 0 : *p_mode_lvq_p = -1;
1263 : }
1264 : }
1265 : }
1266 :
1267 : #ifdef DEBUGGING
1268 : if ( *predmode == -1 )
1269 : {
1270 : return IVAS_ERROR( IVAS_ERR_INTERNAL_FATAL, "\nfind_pred_mode(): incorrect coder_type specification: %d\n", coder_type );
1271 : }
1272 : #endif
1273 :
1274 13257759 : return error;
1275 : }
1276 :
1277 :
1278 : /*---------------------------------------------------------------------------
1279 : * reorder_isf()
1280 : *
1281 : * To make sure that the ISFs are properly ordered and to keep a certain
1282 : * minimum distance between consecutive ISFs.
1283 : *--------------------------------------------------------------------------*/
1284 :
1285 11486 : void reorder_isf(
1286 : float *isf, /* i/o: ISF vector */
1287 : const float min_dist, /* i : minimum required distance */
1288 : const int16_t n, /* i : LPC order */
1289 : const float Fs /* i : sampling frequency */
1290 : )
1291 : {
1292 : int16_t i;
1293 : float isf_min;
1294 : float isf_max;
1295 :
1296 : /*-----------------------------------------------------------------*
1297 : * Verify the ISF ordering and minimum GAP
1298 : *-----------------------------------------------------------------*/
1299 :
1300 11486 : isf_min = min_dist;
1301 183776 : for ( i = 0; i < n - 1; i++ )
1302 : {
1303 172290 : if ( isf[i] < isf_min )
1304 : {
1305 960 : isf[i] = isf_min;
1306 : }
1307 :
1308 172290 : isf_min = isf[i] + min_dist;
1309 : }
1310 :
1311 : /*------------------------------------------------------------------------------------------*
1312 : * Reverify the ISF ordering and minimum GAP in the reverse order (security)
1313 : *------------------------------------------------------------------------------------------*/
1314 :
1315 11486 : isf_max = Fs / 2.0f - min_dist;
1316 :
1317 11486 : if ( isf[n - 2] > isf_max ) /* If danger of unstable filter in case of resonance in HF */
1318 : {
1319 0 : for ( i = n - 2; i >= 0; i-- ) /* Reverify the minimum ISF gap in the reverse sens */
1320 : {
1321 0 : if ( isf[i] > isf_max )
1322 : {
1323 0 : isf[i] = isf_max;
1324 : }
1325 :
1326 0 : isf_max = isf[i] - min_dist;
1327 : }
1328 : }
1329 :
1330 11486 : return;
1331 : }
1332 :
1333 : /*----------------------------------------------------------------------------------*
1334 : * lsf_stab()
1335 : *
1336 : * Check LSF stability (distance between old LSFs and current LSFs)
1337 : *----------------------------------------------------------------------------------*/
1338 :
1339 : /*! r: LP filter stability */
1340 29716853 : float lsf_stab(
1341 : const float *lsf, /* i : LSF vector */
1342 : const float *lsfold, /* i : old LSF vector */
1343 : const int16_t Opt_AMR_WB, /* i : flag indicating AMR-WB IO mode*/
1344 : const int16_t L_frame /* i : frame length */
1345 : )
1346 : {
1347 : int16_t i, m;
1348 : float scale, stab_fac, tmp;
1349 :
1350 29716853 : tmp = 0.0f;
1351 29716853 : if ( Opt_AMR_WB )
1352 : {
1353 10336 : m = M - 1;
1354 : }
1355 : else
1356 : {
1357 29706517 : m = M;
1358 : }
1359 :
1360 505176165 : for ( i = 0; i < m; i++ )
1361 : {
1362 475459312 : tmp += ( lsf[i] - lsfold[i] ) * ( lsf[i] - lsfold[i] );
1363 : }
1364 :
1365 29716853 : scale = (float) L_FRAME / (float) L_frame;
1366 29716853 : scale *= scale;
1367 :
1368 29716853 : stab_fac = (float) ( 1.25f - ( scale * tmp / 400000.0f ) );
1369 :
1370 29716853 : if ( stab_fac > 1.0f )
1371 : {
1372 14257442 : stab_fac = 1.0f;
1373 : }
1374 :
1375 29716853 : if ( stab_fac < 0.0f )
1376 : {
1377 1896386 : stab_fac = 0.0f;
1378 : }
1379 29716853 : return stab_fac;
1380 : }
1381 :
1382 :
1383 : /*---------------------------------------------------------------------
1384 : * get_isppol()
1385 : *
1386 : * Find the polynomial F1(z) or F2(z) from the ISPs.
1387 : * This is performed by expanding the product polynomials:
1388 : *
1389 : * F1(z) = PRODUCT ( 1 - 2 isp_i z^-1 + z^-2 )
1390 : * i=0,2,...,14
1391 : * F2(z) = PRODUCT ( 1 - 2 isp_i z^-1 + z^-2 )
1392 : * i=1,3,...,13
1393 : *
1394 : * where isp_i are the ISPs in the cosine domain.
1395 : *---------------------------------------------------------------------*/
1396 :
1397 147120 : static void get_isppol(
1398 : const float *isp, /* i : Immitance spectral pairs (cosine domaine) */
1399 : float f[], /* o : the coefficients of F1 or F2 */
1400 : const int16_t n /* i : nb. of coefficients (m/2) */
1401 : )
1402 : {
1403 : float b;
1404 : int16_t i, j;
1405 :
1406 147120 : f[0] = 1.0f;
1407 147120 : b = (float) ( -2.0f * *isp );
1408 147120 : f[1] = b;
1409 1103400 : for ( i = 2; i <= n; i++ )
1410 : {
1411 956280 : isp += 2;
1412 956280 : b = (float) ( -2.0f * *isp );
1413 956280 : f[i] = (float) ( b * f[i - 1] + 2.0f * f[i - 2] );
1414 :
1415 3604440 : for ( j = i - 1; j > 1; j-- )
1416 : {
1417 2648160 : f[j] += b * f[j - 1] + f[j - 2];
1418 : }
1419 956280 : f[1] += b;
1420 : }
1421 :
1422 147120 : return;
1423 : }
1424 :
1425 :
1426 : /*---------------------------------------------------------------------
1427 : * chebps2()
1428 : *
1429 : * Evaluates the Chebyshev polynomial series
1430 : *
1431 : * The polynomial order is
1432 : * n = m/2 (m is the prediction order)
1433 : * The polynomial is given by
1434 : * C(x) = f(0)T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2
1435 : *---------------------------------------------------------------------*/
1436 :
1437 : /*! r: the value of the polynomial C(x) */
1438 7418053914 : static float chebps2(
1439 : const float x, /* i : value of evaluation; x=cos(freq) */
1440 : const float *f, /* i : coefficients of sum or diff polyn. */
1441 : const int16_t n /* i : order of polynomial */
1442 : )
1443 : {
1444 : float b1, b2, b0, x2;
1445 : int16_t i;
1446 :
1447 7418053914 : x2 = (float) ( 2.0f * x );
1448 7418053914 : b2 = f[0];
1449 :
1450 7418053914 : b1 = x2 * b2 + f[1];
1451 :
1452 51926363915 : for ( i = 2; i < n; i++ )
1453 : {
1454 44508310001 : b0 = x2 * b1 - b2 + f[i];
1455 44508310001 : b2 = b1;
1456 44508310001 : b1 = b0;
1457 : }
1458 :
1459 7418053914 : return (float) ( x * b1 - b2 + 0.5f * f[n] );
1460 : }
1461 :
1462 :
1463 : /*---------------------------------------------------------------------
1464 : * LPC_chebyshev()
1465 : *
1466 : * Evaluate a series expansion in Chebyshev polynomials
1467 : *
1468 : * The polynomial order is
1469 : * n = m/2 (m is the prediction order)
1470 : * The polynomial is given by
1471 : * C(x) = T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2
1472 : *---------------------------------------------------------------------*/
1473 :
1474 150438399 : static float LPC_chebyshev(
1475 : float x,
1476 : float f[],
1477 : const int16_t n )
1478 : {
1479 : float b1, b2, b0, x2, val;
1480 : int16_t i;
1481 :
1482 150438399 : x2 = 2.0f * x;
1483 150438399 : b2 = 1.0f;
1484 150438399 : b1 = x2 + f[0];
1485 :
1486 570715060 : for ( i = 2; i < n; i++ )
1487 : {
1488 420276661 : b0 = x2 * b1 - b2 + f[i - 1];
1489 :
1490 : /* was previously casting x2 into float to have this
1491 : equation evaluated in float to be same
1492 : same as EVRC-B only code which has 2.0 in equation
1493 : instead of a float x2. This was causing non-bit-exactness
1494 : in a very very very rare corner case.
1495 : Doesnt really matter, but just to be picky! */
1496 420276661 : b2 = b1;
1497 420276661 : b1 = b0;
1498 : }
1499 :
1500 150438399 : val = ( x * b1 - b2 + f[i - 1] );
1501 :
1502 150438399 : return val;
1503 : }
1504 :
1505 :
1506 : /*-------------------------------------------------------------------*
1507 : * lsp2isp()
1508 : *
1509 : * Convert LSPs to ISPs via predictor coefficients A[]
1510 : *-------------------------------------------------------------------*/
1511 :
1512 151 : void lsp2isp(
1513 : const float *lsp, /* i : LSP vector */
1514 : float *isp, /* o : ISP filter coefficients */
1515 : float *stable_isp, /* i/o: ISP filter coefficients */
1516 : const int16_t m /* i : order of LP analysis */
1517 : )
1518 : {
1519 : float a[M + 1];
1520 :
1521 : /* LSP --> A */
1522 151 : lsp2a_stab( lsp, a, m );
1523 :
1524 : /* A --> ISP */
1525 151 : a2isp( a, isp, stable_isp );
1526 :
1527 : /* Update to latest stable ISP */
1528 151 : mvr2r( isp, stable_isp, M );
1529 151 : }
1530 :
1531 : /*-------------------------------------------------------------------*
1532 : * isp2lsp()
1533 : *
1534 : * Convert ISPs to LSPs via predictor coefficients A[]
1535 : *-------------------------------------------------------------------*/
1536 :
1537 187 : void isp2lsp(
1538 : const float *isp, /* i : LSP vector */
1539 : float *lsp, /* o : ISP filter coefficients */
1540 : float *stable_lsp, /* i/o: stable LSP filter coefficients */
1541 : const int16_t m /* i : order of LP analysis */
1542 : )
1543 : {
1544 : float a[M + 1];
1545 :
1546 : /* ISP --> A */
1547 187 : isp2a( isp, a, m );
1548 :
1549 : /* A --> LSP */
1550 187 : a2lsp_stab( a, lsp, stable_lsp );
1551 :
1552 : /* Update to latest stable LSP */
1553 187 : mvr2r( lsp, stable_lsp, M );
1554 187 : }
1555 :
1556 :
1557 : /*-------------------------------------------------------------------*
1558 : * lsf2isf()
1559 : *
1560 : * Convert LSPs to ISPs
1561 : *-------------------------------------------------------------------*/
1562 :
1563 41 : void lsf2isf(
1564 : const float *lsf, /* i : LSF vector */
1565 : float *isf, /* o : ISF vector */
1566 : float *stable_isp, /* i/o: stable ISP filter coefficients */
1567 : const int16_t m, /* i : order of LP analysis */
1568 : const int32_t Fs )
1569 : {
1570 : float tmp_lsp[M];
1571 : float tmp_isp[M];
1572 :
1573 : /* LSF --> LSP */
1574 41 : lsf2lsp( lsf, tmp_lsp, m, Fs );
1575 :
1576 : /* LSP --> ISP */
1577 41 : lsp2isp( tmp_lsp, tmp_isp, stable_isp, m );
1578 :
1579 : /* ISP --> ISF */
1580 41 : isp2isf( tmp_isp, isf, m, Fs );
1581 :
1582 41 : return;
1583 : }
1584 :
1585 : /*-------------------------------------------------------------------*
1586 : * isf2lsf()
1587 : *
1588 : * Convert ISFs to LSFs
1589 : *-------------------------------------------------------------------*/
1590 :
1591 54 : void isf2lsf(
1592 : const float *isf, /* i : ISF vector */
1593 : float *lsf, /* o : LSF vector */
1594 : float *stable_lsp, /* i/o: stable LSP filter coefficients */
1595 : const int16_t m, /* i : order of LP analysis */
1596 : const int32_t Fs )
1597 : {
1598 : float tmp_isp[M];
1599 : float tmp_lsp[M];
1600 :
1601 : /* ISF --> ISP */
1602 54 : isf2isp( isf, tmp_isp, m, Fs );
1603 :
1604 : /* ISP --> LSP */
1605 54 : isp2lsp( tmp_isp, tmp_lsp, stable_lsp, m );
1606 :
1607 : /* LSP --> LSF */
1608 54 : lsp2lsf( tmp_lsp, lsf, m, Fs );
1609 :
1610 54 : return;
1611 : }
1612 :
1613 : /*-----------------------------------------------------------------------*
1614 : * lsp2lsf()
1615 : *
1616 : * Transformation of LSPs to LSFs
1617 : *
1618 : * LSP are line spectral pair in cosine domain (-1 to 1).
1619 : * LSF are line spectral frequencies (0 to fs/2).
1620 : *-----------------------------------------------------------------------*/
1621 :
1622 37821165 : void lsp2lsf(
1623 : const float lsp[], /* i : isp[m] (range: -1<=val<1) */
1624 : float lsf[], /* o : isf[m] normalized (range: 0<=val<=fs/2) */
1625 : const int16_t m, /* i : LPC order */
1626 : const int32_t Fs /* i : sampling frequency */
1627 : )
1628 : {
1629 : int16_t i;
1630 :
1631 : /* convert LSPs to LSFs */
1632 611111553 : for ( i = 0; i < m; i++ )
1633 : {
1634 573290388 : lsf[i] = (float) ( acos( lsp[i] ) * ( ( Fs / 2.f ) / EVS_PI ) );
1635 : }
1636 :
1637 37821165 : return;
1638 : }
1639 :
1640 : /*-----------------------------------------------------------------------*
1641 : * lsf2lsp()
1642 : *
1643 : * Transformation of LSFs to LSPs
1644 : *
1645 : * LSP are line spectral pairs in cosine domain (-1 to 1).
1646 : * LSF are line spectral frequencies (0 to fs/2).
1647 : *-----------------------------------------------------------------------*/
1648 :
1649 24603462 : void lsf2lsp(
1650 : const float lsf[], /* i : isf[m] normalized (range: 0<=val<=fs/2) */
1651 : float lsp[], /* o : isp[m] (range: -1<=val<1) */
1652 : const int16_t m, /* i : LPC order */
1653 : const int32_t Fs /* i : sampling frequency */
1654 : )
1655 : {
1656 : int16_t i;
1657 :
1658 : /* convert LSFs to LSPs */
1659 400025910 : for ( i = 0; i < m; i++ )
1660 : {
1661 375422448 : lsp[i] = (float) cos( lsf[i] * EVS_PI / ( Fs / 2.f ) );
1662 : }
1663 :
1664 24603462 : return;
1665 : }
1666 :
1667 :
1668 : /*---------------------------------------------------------------------------
1669 : * tcvq_Dec()
1670 : *
1671 : *
1672 : *--------------------------------------------------------------------------*/
1673 :
1674 13090 : static void tcvq_Dec(
1675 : int16_t *ind,
1676 : float *d_out,
1677 : const int16_t safety_net )
1678 : {
1679 : int16_t i;
1680 : int16_t index[9];
1681 : int16_t stage, state, branch[N_STAGE], codeword[N_STAGE];
1682 : int16_t fins, iwd;
1683 : float pred[N_DIM];
1684 : float D[N_STAGE_VQ][N_DIM];
1685 : const float( *TCVQ_CB_SUB1 )[128][2], ( *TCVQ_CB_SUB2 )[64][2], ( *TCVQ_CB_SUB3 )[32][2];
1686 : const float( *IntraCoeff )[2][2];
1687 :
1688 13090 : mvs2s( ind, index, 9 );
1689 :
1690 13090 : if ( safety_net )
1691 : {
1692 2154 : TCVQ_CB_SUB1 = SN_TCVQ_CB_SUB1;
1693 2154 : TCVQ_CB_SUB2 = SN_TCVQ_CB_SUB2;
1694 2154 : TCVQ_CB_SUB3 = SN_TCVQ_CB_SUB3;
1695 2154 : IntraCoeff = SN_IntraCoeff;
1696 : }
1697 : else
1698 : {
1699 10936 : TCVQ_CB_SUB1 = AR_TCVQ_CB_SUB1;
1700 10936 : TCVQ_CB_SUB2 = AR_TCVQ_CB_SUB2;
1701 10936 : TCVQ_CB_SUB3 = AR_TCVQ_CB_SUB3;
1702 10936 : IntraCoeff = AR_IntraCoeff;
1703 : }
1704 :
1705 : /* Decode Final state */
1706 13090 : fins = index[0] & 15;
1707 :
1708 : /* Decode Branch info */
1709 13090 : branch[0] = index[1] >> 4;
1710 13090 : branch[1] = index[2] >> 4;
1711 :
1712 39270 : for ( stage = 2; stage < N_STAGE_VQ - 4; stage++ )
1713 : {
1714 26180 : branch[stage] = index[stage + 1] >> 3;
1715 : }
1716 :
1717 13090 : branch[4] = fins & 0x1;
1718 13090 : branch[5] = ( fins >> 1 ) & 0x1;
1719 13090 : branch[6] = ( fins >> 2 ) & 0x1;
1720 13090 : branch[7] = ( fins >> 3 ) & 0x1;
1721 :
1722 : /* Decode Codeword info */
1723 39270 : for ( stage = 0; stage < 2; stage++ )
1724 : {
1725 26180 : codeword[stage] = ( index[stage + 1] & 15 ) << 3;
1726 : }
1727 :
1728 39270 : for ( stage = 2; stage < N_STAGE_VQ - 4; stage++ )
1729 : {
1730 26180 : codeword[stage] = ( index[stage + 1] & 7 ) << 3;
1731 : }
1732 :
1733 65450 : for ( stage = N_STAGE_VQ - 4; stage < N_STAGE_VQ; stage++ )
1734 : {
1735 52360 : codeword[stage] = ( index[stage + 1] & 3 ) << 3;
1736 : }
1737 :
1738 13090 : state = ( fins >> 2 ) << 2;
1739 :
1740 : /* stage #1 */
1741 13090 : iwd = NTRANS2[branch[0] + 2][state] + codeword[0];
1742 13090 : D[0][0] = TCVQ_CB_SUB1[0][iwd][0];
1743 13090 : D[0][1] = TCVQ_CB_SUB1[0][iwd][1];
1744 13090 : state = NTRANS2[branch[0]][state];
1745 :
1746 : /* stage #2 */
1747 13090 : pred[0] = IntraCoeff[0][0][0] * D[0][0] + IntraCoeff[0][0][1] * D[0][1];
1748 13090 : pred[1] = IntraCoeff[0][1][0] * D[0][0] + IntraCoeff[0][1][1] * D[0][1];
1749 :
1750 13090 : iwd = NTRANS2[branch[1] + 2][state] + codeword[1];
1751 13090 : D[1][0] = TCVQ_CB_SUB1[1][iwd][0] + pred[0];
1752 13090 : D[1][1] = TCVQ_CB_SUB1[1][iwd][1] + pred[1];
1753 13090 : state = NTRANS2[branch[1]][state];
1754 :
1755 : /* stage #3 - #4 */
1756 39270 : for ( stage = 2; stage < N_STAGE_VQ - 4; stage++ )
1757 : {
1758 26180 : pred[0] = IntraCoeff[stage - 1][0][0] * D[stage - 1][0] + IntraCoeff[stage - 1][0][1] * D[stage - 1][1];
1759 26180 : pred[1] = IntraCoeff[stage - 1][1][0] * D[stage - 1][0] + IntraCoeff[stage - 1][1][1] * D[stage - 1][1];
1760 :
1761 26180 : iwd = NTRANS2[branch[stage] + 2][state] + codeword[stage];
1762 26180 : D[stage][0] = TCVQ_CB_SUB2[stage - 2][iwd][0] + pred[0];
1763 26180 : D[stage][1] = TCVQ_CB_SUB2[stage - 2][iwd][1] + pred[1];
1764 26180 : state = NTRANS2[branch[stage]][state];
1765 : }
1766 :
1767 : /* stage #5 - #8 */
1768 65450 : for ( stage = N_STAGE_VQ - 4; stage < N_STAGE_VQ; stage++ )
1769 : {
1770 52360 : pred[0] = IntraCoeff[stage - 1][0][0] * D[stage - 1][0] + IntraCoeff[stage - 1][0][1] * D[stage - 1][1];
1771 52360 : pred[1] = IntraCoeff[stage - 1][1][0] * D[stage - 1][0] + IntraCoeff[stage - 1][1][1] * D[stage - 1][1];
1772 :
1773 52360 : iwd = NTRANS2[branch[stage] + 2][state] + codeword[stage];
1774 52360 : D[stage][0] = TCVQ_CB_SUB3[stage - 4][iwd][0] + pred[0];
1775 52360 : D[stage][1] = TCVQ_CB_SUB3[stage - 4][iwd][1] + pred[1];
1776 52360 : state = NTRANS2[branch[stage]][state];
1777 : }
1778 :
1779 117810 : for ( stage = 0; stage < N_STAGE_VQ; stage++ )
1780 : {
1781 314160 : for ( i = 0; i < N_DIM; i++ )
1782 : {
1783 209440 : d_out[( N_DIM * stage ) + i] = D[stage][i];
1784 : }
1785 : }
1786 13090 : return;
1787 : }
1788 :
1789 : /*---------------------------------------------------------------------------
1790 : * qlsf_ARSN_tcvq_Dec_16k()
1791 : *
1792 : * Predictive BC-TCQ encoder for LSF quantization
1793 : *--------------------------------------------------------------------------*/
1794 :
1795 13090 : int16_t qlsf_ARSN_tcvq_Dec_16k(
1796 : float *y, /* o : Quantized LSF vector */
1797 : int16_t *indice, /* i : Indices */
1798 : const int16_t nBits /* i : number of bits */
1799 : )
1800 : {
1801 : int16_t i;
1802 : float error_svq_q[M];
1803 : int16_t safety_net;
1804 :
1805 : /* Select Mode */
1806 13090 : safety_net = indice[0];
1807 :
1808 13090 : if ( safety_net == 1 )
1809 : {
1810 2154 : tcvq_Dec( &indice[1], y, safety_net );
1811 :
1812 2154 : if ( nBits > 30 )
1813 : {
1814 10170 : for ( i = 0; i < 8; i++ )
1815 : {
1816 9040 : error_svq_q[i] = AR_SVQ_CB1[indice[10]][i];
1817 9040 : error_svq_q[i + 8] = AR_SVQ_CB2[indice[11]][i];
1818 : }
1819 :
1820 19210 : for ( i = 0; i < M; i++ )
1821 : {
1822 18080 : y[i] = y[i] + error_svq_q[i] * scale_ARSN[i];
1823 : }
1824 : }
1825 : }
1826 : else
1827 : {
1828 10936 : tcvq_Dec( &indice[1], y, safety_net );
1829 :
1830 10936 : if ( nBits > 30 )
1831 : {
1832 54999 : for ( i = 0; i < 8; i++ )
1833 : {
1834 48888 : error_svq_q[i] = AR_SVQ_CB1[indice[10]][i];
1835 48888 : error_svq_q[i + 8] = AR_SVQ_CB2[indice[11]][i];
1836 : }
1837 :
1838 103887 : for ( i = 0; i < M; i++ )
1839 : {
1840 97776 : y[i] = y[i] + error_svq_q[i];
1841 : }
1842 : }
1843 : }
1844 :
1845 13090 : return safety_net;
1846 : }
1847 :
1848 :
1849 : /*-------------------------------------------------------------------*
1850 : * lsf_syn_mem_backup()
1851 : *
1852 : * back-up synthesis filter memory and LSF qunatizer memories (used in SC-VBR)
1853 : *-------------------------------------------------------------------*/
1854 :
1855 1658570 : void lsf_syn_mem_backup(
1856 : Encoder_State *st, /* i : state structure */
1857 : float *btilt_code, /* i : tilt code */
1858 : float *bgc_threshold, /* i : */
1859 : float *clip_var_bck, /* o : */
1860 : int16_t *next_force_sf_bck, /* o : */
1861 : float *lsp_new, /* i : LSP vector to quantize */
1862 : float *lsp_mid, /* i : mid-frame LSP vector */
1863 : float *clip_var, /* o : pitch clipping state var */
1864 : float *mem_AR, /* o : quantizer memory for AR model */
1865 : float *mem_MA, /* o : quantizer memory for AR model */
1866 : float *lsp_new_bck, /* o : LSP vector to quantize- backup */
1867 : float *lsp_mid_bck, /* o : mid-frame LSP vector - backup */
1868 : float *Bin_E, /* o : FFT Bin energy 128 *2 sets */
1869 : float *Bin_E_old, /* o : FFT Bin energy 128 sets */
1870 : float *mem_syn_bck, /* o : synthesis filter memory */
1871 : float *mem_w0_bck, /* o : memory of the weighting filter */
1872 : float *streaklimit,
1873 : int16_t *pstreaklen )
1874 : {
1875 : int16_t i;
1876 1658570 : LPD_state_HANDLE hLPDmem = st->hLPDmem;
1877 :
1878 1658570 : *clip_var = st->clip_var[0];
1879 :
1880 28195690 : for ( i = 0; i < M; i++ )
1881 : {
1882 26537120 : mem_AR[i] = st->mem_AR[i];
1883 26537120 : mem_MA[i] = st->mem_MA[i];
1884 26537120 : lsp_new_bck[i] = lsp_new[i];
1885 26537120 : lsp_mid_bck[i] = lsp_mid[i];
1886 : }
1887 :
1888 1658570 : *streaklimit = st->streaklimit;
1889 1658570 : *pstreaklen = st->pstreaklen;
1890 :
1891 426252490 : for ( i = 0; i < L_FFT; i++ )
1892 : {
1893 424593920 : Bin_E[i] = st->Bin_E[i];
1894 : }
1895 :
1896 213955530 : for ( i = 0; i < ( L_FFT / 2 ); i++ )
1897 : {
1898 212296960 : Bin_E_old[i] = st->Bin_E_old[i];
1899 : }
1900 :
1901 : /* back-up memories */
1902 28195690 : for ( i = 0; i < M; i++ )
1903 : {
1904 26537120 : mem_syn_bck[i] = hLPDmem->mem_syn[i];
1905 : }
1906 :
1907 1658570 : *mem_w0_bck = hLPDmem->mem_w0;
1908 :
1909 1658570 : *btilt_code = hLPDmem->tilt_code;
1910 1658570 : *bgc_threshold = hLPDmem->gc_threshold;
1911 1658570 : mvr2r( st->clip_var, clip_var_bck, 6 );
1912 1658570 : *next_force_sf_bck = st->next_force_safety_net;
1913 :
1914 1658570 : return;
1915 : }
1916 :
1917 :
1918 : /*-------------------------------------------------------------------*
1919 : * lsf_syn_mem_restore()
1920 : *
1921 : * restore synthesis filter memory and LSF quantizer memories
1922 : *-------------------------------------------------------------------*/
1923 :
1924 357 : void lsf_syn_mem_restore(
1925 : Encoder_State *st, /* o : state structure */
1926 : float btilt_code, /* i : */
1927 : float gc_threshold, /* i : */
1928 : float *clip_var_bck, /* i : */
1929 : int16_t next_force_sf_bck, /* i : */
1930 : float *lsp_new, /* o : LSP vector to quantize */
1931 : float *lsp_mid, /* o : mid-frame LSP vector */
1932 : float clip_var, /* i : pitch clipping state var */
1933 : float *mem_AR, /* i : quantizer memory for AR model */
1934 : float *mem_MA, /* i : quantizer memory for AR model */
1935 : float *lsp_new_bck, /* i : LSP vector to quantize- backup */
1936 : float *lsp_mid_bck, /* i : mid-frame LSP vector - backup */
1937 : float *Bin_E, /* i : FFT Bin energy 128 *2 sets */
1938 : float *Bin_E_old, /* i : FFT Bin energy 128 sets */
1939 : float *mem_syn_bck, /* i : synthesis filter memory */
1940 : float mem_w0_bck, /* i : memory of the weighting filter */
1941 : const float streaklimit,
1942 : const int16_t pstreaklen )
1943 : {
1944 : int16_t i;
1945 357 : LPD_state_HANDLE hLPDmem = st->hLPDmem;
1946 :
1947 : /* restore lsf memories */
1948 357 : st->clip_var[0] = clip_var;
1949 :
1950 6069 : for ( i = 0; i < M; i++ )
1951 : {
1952 5712 : st->mem_AR[i] = mem_AR[i];
1953 5712 : st->mem_MA[i] = mem_MA[i];
1954 5712 : lsp_new[i] = lsp_new_bck[i];
1955 5712 : lsp_mid[i] = lsp_mid_bck[i];
1956 : }
1957 :
1958 357 : st->streaklimit = streaklimit;
1959 357 : st->pstreaklen = pstreaklen;
1960 :
1961 91749 : for ( i = 0; i < L_FFT; i++ )
1962 : {
1963 91392 : st->Bin_E[i] = Bin_E[i];
1964 : }
1965 :
1966 46053 : for ( i = 0; i < ( L_FFT / 2 ); i++ )
1967 : {
1968 45696 : st->Bin_E_old[i] = Bin_E_old[i];
1969 : }
1970 :
1971 : /* restoring memories */
1972 357 : hLPDmem->mem_w0 = mem_w0_bck;
1973 :
1974 6069 : for ( i = 0; i < M; i++ )
1975 : {
1976 5712 : hLPDmem->mem_syn[i] = mem_syn_bck[i];
1977 : }
1978 :
1979 357 : hLPDmem->tilt_code = btilt_code;
1980 357 : hLPDmem->gc_threshold = gc_threshold;
1981 357 : mvr2r( clip_var_bck, st->clip_var, 6 );
1982 357 : st->next_force_safety_net = next_force_sf_bck;
1983 :
1984 357 : return;
1985 : }
1986 :
1987 :
1988 : /*--------------------------------------------------------------------------*
1989 : * lsf_update_memory()
1990 : *
1991 : *
1992 : *--------------------------------------------------------------------------*/
1993 :
1994 635263 : void lsf_update_memory(
1995 : const int16_t narrowband, /* i : narrowband flag */
1996 : const float qlsf[], /* i : quantized lsf coefficients */
1997 : float old_mem_MA[], /* i : MA memory */
1998 : float mem_MA[] /* o : updated MA memory */
1999 : )
2000 : {
2001 : int16_t i;
2002 :
2003 10799471 : for ( i = 0; i < M; ++i )
2004 : {
2005 10164208 : mem_MA[i] = qlsf[i] - lsf_means[narrowband][i] - MU_MA * old_mem_MA[i];
2006 : }
2007 :
2008 635263 : return;
2009 : }
2010 :
2011 : /*--------------------------------------------------------------------------*
2012 : * tcxlpc_get_cdk()
2013 : *
2014 : *
2015 : *--------------------------------------------------------------------------*/
2016 :
2017 : /*! r: codebook index */
2018 1006433 : int16_t tcxlpc_get_cdk(
2019 : const int16_t coder_type /* i : GC/VC indicator */
2020 : )
2021 : {
2022 : int16_t cdk;
2023 :
2024 1006433 : cdk = ( coder_type == VOICED );
2025 :
2026 1006433 : return cdk;
2027 : }
2028 :
2029 :
2030 : /*--------------------------------------------------------------------------*
2031 : * dec_FDCNG_MSVQ_stage1()
2032 : *
2033 : *
2034 : *--------------------------------------------------------------------------*/
2035 :
2036 127198 : void dec_FDCNG_MSVQ_stage1(
2037 : int16_t j_full, /* i : index full range */
2038 : int16_t n, /* i : dimension to generate */
2039 : const float *invTrfMatrix, /* i : IDCT matrix for synthesis */
2040 : const DCTTYPE idcttype, /* i : specify which IDCT */
2041 : float *uq, /* o : synthesized stage1 vector */
2042 : Word16 *uq_ind /* o : synthesized stage1 vector in BASOP */
2043 : )
2044 : {
2045 : int16_t col, segm_ind, j;
2046 : float dct_vec[FDCNG_VQ_MAX_LEN];
2047 : float idct_vec[FDCNG_VQ_MAX_LEN];
2048 : const Word8 *cbpW8;
2049 : const Word16 *dct_col_shift_tab;
2050 :
2051 127198 : assert( n <= FDCNG_VQ_MAX_LEN );
2052 127198 : assert( n >= FDCNG_VQ_DCT_MINTRUNC );
2053 :
2054 127198 : segm_ind = 0;
2055 635990 : for ( col = 1; col <= FDCNG_VQ_DCT_NSEGM; col++ )
2056 : {
2057 508792 : if ( j_full >= cdk1_ivas_cum_entries_per_segment[col] )
2058 : {
2059 270494 : segm_ind++;
2060 : }
2061 : }
2062 :
2063 127198 : j = j_full - cdk1_ivas_cum_entries_per_segment[segm_ind]; /* j is the local segment index */
2064 :
2065 127198 : assert( j < cdk1_ivas_entries_per_segment[segm_ind] );
2066 :
2067 : /* Word8 column variable Qx storage*/
2068 127198 : cbpW8 = cdk_37bits_ivas_stage1_W8Qx_dct_sections[segm_ind]; /* Word8 storage fixed ptr_init */
2069 127198 : cbpW8 += j * cdk1_ivas_cols_per_segment[segm_ind]; /* adaptive ptr init */
2070 127198 : dct_col_shift_tab = stage1_dct_col_syn_shift[segm_ind];
2071 :
2072 2040866 : for ( col = 0; col < cdk1_ivas_cols_per_segment[segm_ind]; col++ )
2073 : {
2074 1913668 : dct_vec[col] = (float) shl( (Word16) cbpW8[col], dct_col_shift_tab[col] );
2075 : /* LOGIC( 1 ) , SHIFT( 1 );
2076 : in BASOP: s_and(for W8->W16), shl()
2077 : */
2078 : }
2079 127198 : dctT2_N_apply_matrix( (const float *) dct_vec, idct_vec, cdk1_ivas_cols_per_segment[segm_ind], n, invTrfMatrix, FDCNG_VQ_DCT_MAXTRUNC, idcttype );
2080 :
2081 : /*scale down to original fdcngvq domain and move to Q0 */
2082 127198 : v_multc( idct_vec, fdcng_dct_scaleF[1], idct_vec, n );
2083 : /* fdcng_dct_scaleF[1] --> 0.0625-->scale down from search Q4 domain to Q0 ,
2084 : not really relevant for BASOP loop */
2085 :
2086 : /*add common mid fdcng vector, in fdcng bands domain */
2087 127198 : v_add( idct_vec, cdk1r_tr_midQ_truncQ, uq, n );
2088 127198 : assert( uq_ind == NULL );
2089 :
2090 127198 : return;
2091 : }
2092 :
2093 :
2094 : /*--------------------------------------------------------------------------*
2095 : * msvq_dec()
2096 : *
2097 : *
2098 : *--------------------------------------------------------------------------*/
2099 :
2100 4948252 : void msvq_dec(
2101 : const float *const *cb, /* i : Codebook (indexed cb[*stages][levels][p]) */
2102 : const int16_t dims[], /* i : Dimension of each codebook stage (NULL: full dim.) */
2103 : const int16_t offs[], /* i : Starting dimension of each codebook stage (NULL: 0) */
2104 : const int16_t stages, /* i : Number of stages */
2105 : const int16_t N, /* i : Vector dimension */
2106 : const int16_t maxN, /* i : Codebook dimension */
2107 : const int16_t Idx[], /* i : Indices */
2108 : const int16_t applyIDCT_flag, /* i : applyIDCT flag */
2109 : const float *invTrfMatrix, /* i : matrix for IDCT synthesis */
2110 : float *uq, /* o : quantized vector */
2111 : Word16 *uq_ind /* o : quantized vector (fixed point) */
2112 : )
2113 : {
2114 : int16_t i, n, maxn, start;
2115 : Word16 j;
2116 4948252 : set_zero( uq, N );
2117 :
2118 4948252 : if ( uq_ind )
2119 : {
2120 18805485 : for ( i = 0; i < N; ++i )
2121 : {
2122 17699280 : uq_ind[i] = 0;
2123 17699280 : move16();
2124 : }
2125 : }
2126 20574415 : for ( i = 0; i < stages; i++ )
2127 : {
2128 15626163 : if ( dims )
2129 : {
2130 3119071 : n = dims[i];
2131 3119071 : maxn = n;
2132 : }
2133 : else
2134 : {
2135 12507092 : n = N;
2136 12507092 : maxn = maxN;
2137 : }
2138 15626163 : if ( offs )
2139 : {
2140 3119071 : start = offs[i];
2141 : }
2142 : else
2143 : {
2144 12507092 : start = 0;
2145 : }
2146 :
2147 15626163 : if ( i == 0 && applyIDCT_flag != 0 )
2148 : {
2149 30302 : assert( start == 0 );
2150 30302 : dec_FDCNG_MSVQ_stage1( Idx[0], N, invTrfMatrix, IDCT_T2_XX_24, uq, uq_ind ); /* IDCT_T2 N=24 used for all synthesis */
2151 : }
2152 : else
2153 : {
2154 15595861 : v_add( uq + start, cb[i] + Idx[i] * maxn, uq + start, n );
2155 : }
2156 :
2157 : #define WMC_TOOL_SKIP
2158 15626163 : IF( uq_ind != NULL )
2159 : {
2160 35524471 : FOR( j = 0; j < n; ++j )
2161 : {
2162 32405400 : move16();
2163 32405400 : uq_ind[start + j] = add( uq_ind[start + j], (Word16) ( cb[i][Idx[i] * maxn + j] * 2.0f * 1.28f ) );
2164 : }
2165 : }
2166 : #undef WMC_TOOL_SKIP
2167 : }
2168 :
2169 4948252 : return;
2170 : }
2171 :
2172 :
2173 : /*--------------------------------------------------------------------------*
2174 : * spec2isf()
2175 : *
2176 : *
2177 : *--------------------------------------------------------------------------*/
2178 :
2179 5296 : static void spec2isf(
2180 : float spec_r[], /* input spectrum real part (only left half + one zero)*/
2181 : float spec_i[], /* input spectrum imag part (only left half+right halt with zeros)*/
2182 : const int16_t speclen, /* length of spectrum (only left half)*/
2183 : float lsf[],
2184 : /* locations of LSFs (buffer must be sufficiently long) */ /*15Q16*/
2185 : const float old_lsf[] /* locations of LSFs (buffer must be sufficiently long) */ /*15Q16*/
2186 : )
2187 : {
2188 : /*spec_r[] needs a 0 in the end!*/
2189 : float s;
2190 : int16_t specix, lsfix, i;
2191 :
2192 5296 : specix = lsfix = 0;
2193 5296 : s = spec_r[specix++];
2194 :
2195 47664 : while ( ( specix < speclen ) && lsfix <= 15 )
2196 : {
2197 :
2198 : /*check for next zero crossing*/
2199 372098 : for ( ; s * spec_r[specix] >= 0; specix++ )
2200 : {
2201 : ;
2202 : }
2203 :
2204 42368 : lsf[lsfix++] = ( specix - 1 + spec_r[specix - 1] / ( spec_r[specix - 1] - spec_r[specix] ) ) * ( 12800 / 256 );
2205 :
2206 : /*check for the next zero crossing*/
2207 385230 : for ( ; s * spec_i[specix] >= 0; specix++ )
2208 : {
2209 : ;
2210 : }
2211 :
2212 42368 : lsf[lsfix++] = ( specix - 1 + spec_i[specix - 1] / ( spec_i[specix - 1] - spec_i[specix] ) ) * ( 12800 / 256 );
2213 :
2214 42368 : spec_r[speclen] = s;
2215 42368 : spec_i[speclen] = s;
2216 :
2217 42368 : s = -s;
2218 : }
2219 :
2220 5296 : if ( lsfix < 16 )
2221 : {
2222 0 : for ( i = 0; i < 16; i++ )
2223 : {
2224 0 : lsf[i] = old_lsf[i];
2225 : }
2226 : }
2227 :
2228 5296 : return;
2229 : }
2230 :
2231 : /*--------------------------------------------------------------------------*
2232 : * a2isf()
2233 : *
2234 : *
2235 : *--------------------------------------------------------------------------*/
2236 :
2237 : #define SCALE1_HALF 1018.59161376953f
2238 :
2239 : typedef struct
2240 : {
2241 : float re;
2242 : float im;
2243 : } Pfloat;
2244 :
2245 5296 : void a2isf(
2246 : float *a,
2247 : float *isf,
2248 : const float *old_isf,
2249 : const int16_t lpcOrder )
2250 : {
2251 : float RealFFT[128];
2252 : float ImagFFT[128];
2253 : float RealOut[130];
2254 : float ImagOut[130];
2255 : float *ptrReal;
2256 : float *ptrImag;
2257 : int16_t n, i, j;
2258 : const Pfloat *ptwiddle;
2259 : Pfloat *pwn17, pwn[128], *pwn15, tmpw15;
2260 5296 : int16_t N = 256;
2261 : float s[4];
2262 : float L_tmp, L_tmp1;
2263 : float lpc[19], fftTmpRe[16], fftTmpIm[16];
2264 : Pfloat twid[128];
2265 : float c;
2266 :
2267 5296 : set_zero( fftTmpRe, 16 );
2268 5296 : set_zero( fftTmpIm, 16 );
2269 :
2270 : /* half length FFT */
2271 : /*c = [sum(a) ((-1).^(1:length(a)))*a];*/
2272 :
2273 5296 : L_tmp = 0;
2274 95328 : for ( j = 0; j <= lpcOrder; j++ )
2275 : {
2276 90032 : L_tmp += a[j];
2277 : }
2278 :
2279 5296 : L_tmp1 = 0;
2280 47664 : for ( j = 0; j < lpcOrder / 2; j++ )
2281 : {
2282 42368 : L_tmp1 -= a[2 * j];
2283 42368 : L_tmp1 += a[2 * j + 1];
2284 : }
2285 5296 : L_tmp1 -= a[2 * j];
2286 :
2287 : /*s = [1 -2*(c(1)-c(2))/(c(1)+c(2)) 1]';*/
2288 5296 : s[0] = 1;
2289 5296 : if ( ( L_tmp + L_tmp1 ) != 0 )
2290 : {
2291 5296 : s[1] = -2 * ( ( L_tmp - L_tmp1 ) / ( L_tmp + L_tmp1 ) );
2292 : }
2293 : else
2294 : {
2295 0 : s[1] = 1;
2296 : }
2297 5296 : s[2] = 1;
2298 :
2299 5296 : lpc[0] = a[0] * s[0];
2300 5296 : L_tmp = a[1] * s[0];
2301 5296 : lpc[1] = L_tmp + ( a[1 - 1] * s[1] );
2302 :
2303 :
2304 84736 : for ( n = 2; n < 17; n++ )
2305 : {
2306 79440 : L_tmp = a[n] * s[0];
2307 79440 : L_tmp += ( a[n - 1] * s[1] );
2308 79440 : lpc[n] = L_tmp + ( a[n - 2] * s[2] );
2309 : }
2310 5296 : lpc[18] = a[16] * s[0];
2311 5296 : L_tmp = a[15] * s[0];
2312 5296 : lpc[17] = L_tmp + ( a[16] * s[1] );
2313 :
2314 :
2315 5296 : ptrReal = RealFFT;
2316 5296 : ptrImag = ImagFFT;
2317 :
2318 52960 : for ( j = 0; j < 9; j++ )
2319 : {
2320 47664 : fftTmpRe[j] = lpc[2 * j];
2321 47664 : fftTmpIm[j] = lpc[2 * j + 1];
2322 : }
2323 5296 : fftTmpRe[j] = lpc[2 * j];
2324 5296 : fftTmpIm[j] = 0;
2325 5296 : j++;
2326 :
2327 37072 : for ( ; j < 16; j++ )
2328 : {
2329 31776 : fftTmpRe[j] = 0;
2330 31776 : fftTmpIm[j] = 0;
2331 : }
2332 :
2333 5296 : DoRTFTn( fftTmpRe, fftTmpIm, 16 );
2334 :
2335 90032 : for ( j = 0; j < 16; j++ )
2336 : {
2337 84736 : ptrReal[j * 8] = fftTmpRe[j];
2338 84736 : ptrImag[j * 8] = fftTmpIm[j];
2339 : }
2340 :
2341 5296 : ptrReal++;
2342 5296 : ptrImag++;
2343 :
2344 42368 : for ( i = 1; i < 8; i++ )
2345 : {
2346 : /*X=x(i:8:M/8) .* exp(-j*2*pi*i*(0:M/8-1)/M);*/
2347 37072 : ptwiddle = (const Pfloat *) w_a[i - 1];
2348 :
2349 37072 : fftTmpRe[0] = lpc[0];
2350 37072 : fftTmpIm[0] = lpc[1];
2351 :
2352 333648 : for ( j = 1; j < 9; j++ )
2353 : {
2354 296576 : fftTmpRe[j] = ( lpc[2 * j] * ptwiddle->re ) - ( lpc[2 * j + 1] * ptwiddle->im );
2355 296576 : fftTmpIm[j] = ( lpc[2 * j + 1] * ptwiddle->re ) + ( lpc[2 * j] * ptwiddle->im );
2356 296576 : ptwiddle++;
2357 : }
2358 :
2359 37072 : fftTmpRe[j] = lpc[2 * j] * ptwiddle->re;
2360 37072 : fftTmpIm[j] = lpc[2 * j] * ptwiddle->im;
2361 37072 : ptwiddle++;
2362 37072 : j++;
2363 259504 : for ( ; j < 16; j++ )
2364 : {
2365 222432 : fftTmpRe[j] = 0;
2366 222432 : fftTmpIm[j] = 0;
2367 222432 : ptwiddle++;
2368 : }
2369 :
2370 37072 : DoRTFTn( fftTmpRe, fftTmpIm, 16 );
2371 :
2372 630224 : for ( j = 0; j < 16; j++ )
2373 : {
2374 593152 : ptrReal[j * 8] = fftTmpRe[j];
2375 593152 : ptrImag[j * 8] = fftTmpIm[j];
2376 : }
2377 :
2378 37072 : ptrReal++;
2379 37072 : ptrImag++;
2380 : }
2381 :
2382 5296 : c = EVS_PI / ( 2.0f * (float) 128 );
2383 :
2384 677888 : for ( i = 1; i < 128; i++ )
2385 : {
2386 672592 : twid[i - 1].im = (float) sin( c * ( 2.0f * (float) i ) );
2387 672592 : twid[i - 1].re = (float) cos( c * ( 2.0f * (float) i ) );
2388 : }
2389 5296 : ptwiddle = twid;
2390 :
2391 : /* pre-twiddle */
2392 677888 : for ( i = 1; i < 128; i++ )
2393 : {
2394 672592 : pwn[i - 1].im = (float) sin( c * ( 18.0f * (float) i ) );
2395 672592 : pwn[i - 1].re = (float) cos( c * ( 18.0f * (float) i ) );
2396 : }
2397 :
2398 5296 : pwn17 = pwn;
2399 5296 : pwn15 = &tmpw15;
2400 :
2401 : /*Y(1) = real(X(1)) + imag(X(1));*/
2402 5296 : RealOut[0] = ( RealFFT[0] + ImagFFT[0] );
2403 5296 : ImagOut[0] = 0;
2404 :
2405 : /*Y(N/2+1) = 0.5*(X(1) + conj(X(1))).*exp(pi*i*128*(18)/N) - i*0.5*(X(1) - conj(X(1))).*exp(pi*i*128*(16)/N);*/
2406 5296 : RealOut[128] = 0;
2407 5296 : ImagOut[128] = ( RealFFT[0] + RealFFT[0] ) - ( ImagFFT[0] + ImagFFT[0] );
2408 :
2409 : /*Y(2:N/2) = (0.5*(X(2:N/2) + conj(X(N/2:-1:2))) - i*0.5*(X(2:N/2) - conj(X(N/2:-1:2))).*exp(-pi*i*r*(2)/N)).*exp(pi*i*r*(18)/N);*/
2410 677888 : for ( i = 1; i < N / 2; i++ )
2411 : {
2412 672592 : float ReAr = ( RealFFT[i] + RealFFT[N / 2 - i] );
2413 672592 : float ReBr = ( RealFFT[N / 2 - i] - RealFFT[i] );
2414 672592 : float ImAr = ( ImagFFT[i] - ImagFFT[N / 2 - i] );
2415 672592 : float ImBr = -( ImagFFT[i] + ImagFFT[N / 2 - i] );
2416 :
2417 672592 : tmpw15.re = ( ptwiddle->re * pwn17->re ) + ( ptwiddle->im * pwn17->im );
2418 672592 : tmpw15.im = ( ptwiddle->re * pwn17->im ) - ( ptwiddle->im * pwn17->re );
2419 :
2420 : /*RealOut[i] = mac_r(L_msu(L_msu(L_mult(ReAr, pwn17->re),ImAr, pwn17->im), ReBr, pwn15->v.im), ImBr, pwn15->re); move16();*/
2421 672592 : RealOut[i] = ( ReAr * pwn17->re ) - ( ImAr * pwn17->im ) - ( ( ReBr * pwn15->im ) + ( ImBr * pwn15->re ) );
2422 672592 : ImagOut[i] = ( ReAr * pwn17->im ) + ( ImAr * pwn17->re ) + ( ReBr * pwn15->re ) - ( ImBr * pwn15->im );
2423 :
2424 672592 : ptwiddle++;
2425 672592 : pwn17++;
2426 : }
2427 :
2428 5296 : spec2isf( RealOut, ImagOut, 128, isf, old_isf );
2429 :
2430 5296 : isf[lpcOrder - 1] = (Float32) ( acos( a[lpcOrder] ) * SCALE1_HALF );
2431 :
2432 5296 : return;
2433 : }
2434 :
2435 :
2436 : /*-------------------------------------------------------------------*
2437 : * dctT2_N_apply_matrix()
2438 : *
2439 : * dct/idct truncated matrix appl. for DCT basis vector lengths of N
2440 : *-------------------------------------------------------------------*/
2441 :
2442 142138 : void dctT2_N_apply_matrix(
2443 : const float *input, /* i : input in fdcng or DCT(fdcng) domain */
2444 : float *output, /* o : output in DCT(fdcng) or fdcng ordomain */
2445 : const int16_t dct_dim, /* i : dct processing dim possibly truncated */
2446 : const int16_t fdcngvq_dim, /* i : fdcng domain length */
2447 : const float *matrix, /* i : IDCT matrix */
2448 : const int16_t matrix_row_dim, /* i : */
2449 : const DCTTYPE dcttype /* i : matrix operation type */
2450 : )
2451 : {
2452 : int16_t i, j, dim_in, dim_out;
2453 : int16_t mat_step_col, mat_step_row, mat_step_col_flag;
2454 : const float *pt_x, *pt_A;
2455 : float tmp_y[FDCNG_VQ_MAX_LEN];
2456 : float *pt_y;
2457 :
2458 : /* non-square DCT_N and IDCT_N matrix application,
2459 : using a stored format of an IDCT_Nx(FDCNG_VQ_DCT_MAXTRUNC) matrix */
2460 : /* efficiently parallelized in SIMD */
2461 :
2462 142138 : assert( dct_dim <= FDCNG_VQ_DCT_MAXTRUNC );
2463 142138 : assert( fdcngvq_dim <= FDCNG_VQ_MAX_LEN );
2464 :
2465 142138 : if ( ( dcttype & 1 ) == 0 ) /* even entries are DCTs */
2466 : {
2467 : /* DCT_typeII 24,21 -> XX in worst case */
2468 14940 : dim_in = fdcngvq_dim;
2469 14940 : dim_out = dct_dim;
2470 14940 : mat_step_col = matrix_row_dim; /* matrix maximum storage size dependent, width of first row in matrix */
2471 14940 : mat_step_row = 0;
2472 14940 : mat_step_col_flag = 1;
2473 14940 : assert( dcttype == DCT_T2_21_XX || dcttype == DCT_T2_24_XX );
2474 : }
2475 : else
2476 : {
2477 127198 : assert( ( dcttype & 1 ) != 0 ); /* idct */
2478 127198 : dim_in = dct_dim;
2479 127198 : dim_out = fdcngvq_dim;
2480 127198 : mat_step_col = 1;
2481 127198 : mat_step_row = matrix_row_dim;
2482 127198 : mat_step_col_flag = 0;
2483 127198 : assert( dcttype == IDCT_T2_XX_24 );
2484 : }
2485 :
2486 142138 : pt_y = tmp_y;
2487 3441909 : for ( i = 0; i < dim_out; i++ )
2488 : {
2489 3299771 : pt_x = input;
2490 3299771 : *pt_y = 0;
2491 :
2492 : /* +i(DCT) or +i*maxTrunc(IDCT) */
2493 : #define WMC_TOOL_SKIP
2494 3299771 : pt_A = &( matrix[i * ( mat_step_row + mat_step_col_flag )] ); /* ptr indexing */
2495 : PTR_INIT( 1 );
2496 : #undef WMC_TOOL_SKIP
2497 55186319 : for ( j = 0; j < dim_in; j++ )
2498 : {
2499 : #define WMC_TOOL_SKIP
2500 51886548 : *pt_y += ( *pt_x++ ) * ( *pt_A );
2501 51886548 : pt_A += mat_step_col; /* step +maxtrunc or +1 */ /* ptr indexing*/
2502 : MAC( 1 );
2503 : #undef WMC_TOOL_SKIP
2504 : }
2505 3299771 : pt_y++;
2506 : }
2507 :
2508 142138 : mvr2r( tmp_y, output, dim_out );
2509 :
2510 142138 : return;
2511 : }
2512 :
2513 :
2514 : /*-------------------------------------------------------------------*
2515 : * extend_dctN_input()
2516 : *
2517 : * (inputN, dctN) -> idct(N_ext) idct_N matrix application loop for
2518 : * extending, extrapolating a DCT basis vector length of N to N_ext
2519 : *-------------------------------------------------------------------*/
2520 :
2521 2828 : void extend_dctN_input(
2522 : const float *input, /* i : input in fdcng domain */
2523 : const float *dct_input, /* i : input in dctN(fdcng) domain */
2524 : const int16_t in_dim, /* i : in_dim == N */
2525 : float *ext_sig, /* o : extended output in fdcng domain */
2526 : const int16_t out_dim, /* i : output total dim */
2527 : float *matrix, /* i : idct synthesis matrix N rows, n_cols columns */
2528 : const int16_t n_cols, /* i : number of columns == DCT truncation length */
2529 : const DCTTYPE dcttype /* i : matrix operation type */
2530 : )
2531 : {
2532 : int16_t i, j, i_rev;
2533 2828 : const float( *ptr )[FDCNG_VQ_DCT_MAXTRUNC] = (void *) matrix;
2534 :
2535 : /* stored format is an IDCT_Nx(FDCNG_VQ_DCT_MAXTRUNC) matrix */
2536 2828 : assert( in_dim < FDCNG_VQ_MAX_LEN );
2537 2828 : assert( out_dim <= FDCNG_VQ_MAX_LEN );
2538 2828 : assert( out_dim > in_dim );
2539 2828 : assert( n_cols == FDCNG_VQ_DCT_MAXTRUNC ); /* for *ptr[MAX_TRUNC] adressing*/
2540 2828 : assert( ( dcttype & 1 ) != 0 ); /* idct tables always in use for this basis vector extension */
2541 :
2542 2828 : mvr2r( input, ext_sig, in_dim ); /* copy initial part, i.e. only last/tail parts are extended */
2543 2828 : set_f( &( ext_sig[in_dim] ), 0.0, out_dim - in_dim );
2544 :
2545 2828 : i_rev = in_dim; /*ptr init*/
2546 11312 : for ( i = in_dim; i < out_dim; i++ )
2547 : { /* for each extension sample */
2548 : /* i = 21 22 23;
2549 : i_rev = 20 19 18; for odd dctII reflect basis vector
2550 : */
2551 8484 : i_rev--;
2552 :
2553 161196 : for ( j = 0; j < n_cols; j++ ) /* for each available DCT coeff */
2554 : {
2555 : /* DCTcoeff * reflected basis vector */
2556 : #define WMC_TOOL_SKIP
2557 : /* pure ptr MAC operations */
2558 152712 : ext_sig[i] += dct_input[j] * ptr[i_rev][j]; /* sum up scaled and extended basis vector */
2559 : MAC( 1 );
2560 : #undef WMC_TOOL_SKIP
2561 : }
2562 : }
2563 :
2564 2828 : return;
2565 : }
2566 :
2567 :
2568 : /*-------------------------------------------------------------------*
2569 : * create_IDCT_N_Matrix()
2570 : *
2571 : * inititate idct24 FDCNG_VQ_DCT_MAXTRUNCx N matrix in
2572 : * RAM from a quantized compressed ROM format
2573 : *-------------------------------------------------------------------*/
2574 :
2575 28553 : void create_IDCT_N_Matrix(
2576 : float *inv_matrixFloatQ, /* i/o: RAM buffer */
2577 : const int16_t N, /* i : DCT length, number of time samples */
2578 : const int16_t n_cols, /* i : number of dct coeffs (as DCT may be truncated) */
2579 : const int16_t alloc_size /* i : RAM buffer size in elements */
2580 : )
2581 : {
2582 : int16_t c, c1, r, r_flip, W16_val;
2583 : int16_t len;
2584 : int16_t mat_cpy_size;
2585 : const Word16 *absval_ptr;
2586 : const Word8 *idx_ptr;
2587 : Word16 idx;
2588 28553 : float( *ptr )[FDCNG_VQ_DCT_MAXTRUNC] = (void *) inv_matrixFloatQ; /* fixed number of columns pointers, to simplifies adressing in ANSIC */
2589 :
2590 28553 : absval_ptr = unique_idctT2_24coeffsQ16;
2591 28553 : idx_ptr = idctT2_24_compressed_idx;
2592 28553 : len = FDCNG_VQ_MAX_LEN;
2593 :
2594 28553 : if ( N == FDCNG_VQ_MAX_LEN_WB )
2595 : {
2596 2219 : absval_ptr = unique_idctT2_21coeffsQ16;
2597 2219 : idx_ptr = idctT2_21_compressed_idx;
2598 2219 : len = N;
2599 : }
2600 :
2601 28553 : assert( alloc_size >= ( n_cols * len ) ); /* enough space for the full expanded IDCT matrix */
2602 28553 : assert( N <= len );
2603 :
2604 28553 : mat_cpy_size = ( n_cols ) * ( len >> 1 ); /* NB integer division of "len" */
2605 :
2606 28553 : if ( ( len & 1 ) != 0 )
2607 : { /* odd sized DCT with a non-reflected center row */
2608 2219 : mat_cpy_size += n_cols;
2609 : }
2610 :
2611 6156059 : for ( c = 0; c < mat_cpy_size; c++ )
2612 : {
2613 6127506 : idx = (Word16) ( idx_ptr[c] );
2614 6127506 : W16_val = absval_ptr[abs( idx )];
2615 :
2616 6127506 : if ( idx < 0 )
2617 : {
2618 2624363 : W16_val = -( W16_val );
2619 : }
2620 6127506 : inv_matrixFloatQ[c] = ( +1.52587890625e-05f ) * ( (float) W16_val ); /* 1.0/2.^16 scaling to a float-"Q0" , a scaling that is not done in BASOP */
2621 : }
2622 :
2623 : /* for even number of coeffs DCT24,
2624 : flip symmetry for odd, even is used to save 50% IDCT Table ROM */
2625 : /* for an odd DCT center is not flipped e.g for DCT21 */
2626 :
2627 28553 : assert( n_cols == FDCNG_VQ_DCT_MAXTRUNC );
2628 28553 : assert( ( n_cols & 1 ) == 0 );
2629 :
2630 285530 : for ( c = 0; c < ( n_cols ); c += 2 )
2631 : {
2632 256977 : c1 = c + 1;
2633 256977 : r_flip = len - 1;
2634 3300759 : for ( r = 0; r < ( len / 2 ); r++, r_flip-- )
2635 : {
2636 : #define WMC_TOOL_SKIP
2637 3043782 : ptr[r_flip][c] = ptr[r][c]; /* flipped */
2638 3043782 : ptr[r_flip][c1] = -( ptr[r][c1] ); /* flipped and sign swapped */
2639 : MOVE( 2 );
2640 : MULT( 1 ); /* for negate */
2641 : #undef WMC_TOOL_SKIP
2642 : }
2643 : }
2644 :
2645 28553 : return;
2646 : }
|
