1 /* ------------------------------------------------------------------
2 * Copyright (C) 1998-2009 PacketVideo
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
13 * express or implied.
14 * See the License for the specific language governing permissions
15 * and limitations under the License.
16 * -------------------------------------------------------------------
17 */
18 /****************************************************************************************
19 Portions of this file are derived from the following 3GPP standard:
20
21 3GPP TS 26.073
22 ANSI-C code for the Adaptive Multi-Rate (AMR) speech codec
23 Available from http://www.3gpp.org
24
25 (C) 2004, 3GPP Organizational Partners (ARIB, ATIS, CCSA, ETSI, TTA, TTC)
26 Permission to distribute, modify and use this file under the standard license
27 terms listed above has been obtained from the copyright holder.
28 ****************************************************************************************/
29 /*
30
31 Pathname: ./audio/gsm-amr/c/src/az_lsp.c
32 Funtions: Chebps
33 Chebps_Wrapper
34 Az_lsp
35
36 ------------------------------------------------------------------------------
37 REVISION HISTORY
38
39 Description: Finished first pass of optimization.
40
41 Description: Made changes based on review comments.
42
43 Description: Made input to Chebps_Wrapper consistent with that of Chebps.
44
45 Description: Replaced current Pseudo-code with the UMTS code version 3.2.0.
46 Updated coding template.
47
48 Description: Replaced basic_op.h and oper_32b.h with the header files of the
49 math functions used by the file.
50
51 Description: Made the following changes per comments from Phase 2/3 review:
52 1. Used "-" operator instead of calling sub function in the
53 az_lsp() code.
54 2. Copied detailed function description of az_lsp from the
55 header file.
56 3. Modified local variable definition to one per line.
57 4. Used NC in the definition of f1 and f2 arrays.
58 5. Added curly brackets in the IF statement.
59
60 Description: Changed function interface to pass in a pointer to overflow
61 flag into the function instead of using a global flag. Removed
62 inclusion of unneeded header files.
63
64 Description: For Chebps() and Az_lsp()
65 1. Eliminated unused include files.
66 2. Replaced array addressing by pointers
67 3. Eliminated math operations that unnecessary checked for
68 saturation.
69 4. Eliminated not needed variables
70 5. Eliminated if-else checks for saturation
71 6. Deleted unused function cheps_wraper
72
73 Description: Added casting to eliminate warnings
74
75
76 Who: Date:
77 Description:
78
79 ------------------------------------------------------------------------------
80 MODULE DESCRIPTION
81
82 These modules compute the LSPs from the LP coefficients.
83
84 ------------------------------------------------------------------------------
85 */
86
87
88 /*----------------------------------------------------------------------------
89 ; INCLUDES
90 ----------------------------------------------------------------------------*/
91 #include "az_lsp.h"
92 #include "cnst.h"
93 #include "basic_op.h"
94
95 /*----------------------------------------------------------------------------
96 ; MACROS
97 ; Define module specific macros here
98 ----------------------------------------------------------------------------*/
99
100
101 /*----------------------------------------------------------------------------
102 ; DEFINES
103 ; Include all pre-processor statements here. Include conditional
104 ; compile variables also.
105 ----------------------------------------------------------------------------*/
106 #define NC (M/2) /* M = LPC order, NC = M/2 */
107
108 /*----------------------------------------------------------------------------
109 ; LOCAL FUNCTION DEFINITIONS
110 ; Function Prototype declaration
111 ----------------------------------------------------------------------------*/
112
113
114 /*----------------------------------------------------------------------------
115 ; LOCAL VARIABLE DEFINITIONS
116 ; Variable declaration - defined here and used outside this module
117 ----------------------------------------------------------------------------*/
118
119 /*
120 ------------------------------------------------------------------------------
121 FUNCTION NAME: Chebps
122 ------------------------------------------------------------------------------
123 INPUT AND OUTPUT DEFINITIONS
124
125 Inputs:
126 x = input value (Word16)
127 f = polynomial (Word16)
128 n = polynomial order (Word16)
129
130 pOverflow = pointer to overflow (Flag)
131
132 Outputs:
133 pOverflow -> 1 if the operations in the function resulted in saturation.
134
135 Returns:
136 cheb = Chebyshev polynomial for the input value x.(Word16)
137
138 Global Variables Used:
139 None.
140
141 Local Variables Needed:
142 None.
143
144 ------------------------------------------------------------------------------
145 FUNCTION DESCRIPTION
146
147 This module evaluates the Chebyshev polynomial series.
148 - The polynomial order is n = m/2 = 5
149 - The polynomial F(z) (F1(z) or F2(z)) is given by
150 F(w) = 2 exp(-j5w) C(x)
151 where
152 C(x) = T_n(x) + f(1)T_n-1(x) + ... +f(n-1)T_1(x) + f(n)/2
153 and T_m(x) = cos(mw) is the mth order Chebyshev
154 polynomial ( x=cos(w) )
155 - C(x) for the input x is returned.
156
157 ------------------------------------------------------------------------------
158 REQUIREMENTS
159
160 None.
161
162 ------------------------------------------------------------------------------
163 REFERENCES
164
165 az_lsp.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001
166
167 ------------------------------------------------------------------------------
168 PSEUDO-CODE
169
170 static Word16 Chebps (Word16 x,
171 Word16 f[], // (n)
172 Word16 n)
173 {
174 Word16 i, cheb;
175 Word16 b0_h, b0_l, b1_h, b1_l, b2_h, b2_l;
176 Word32 t0;
177
178 // The reference ETSI code uses a global flag for Overflow. However, in the
179 // actual implementation a pointer to Overflow flag is passed in as a
180 // parameter to the function. This pointer is passed into all the basic math
181 // functions invoked
182
183 b2_h = 256; // b2 = 1.0
184 b2_l = 0;
185
186 t0 = L_mult (x, 512); // 2*x
187 t0 = L_mac (t0, f[1], 8192); // + f[1]
188 L_Extract (t0, &b1_h, &b1_l); // b1 = 2*x + f[1]
189
190 for (i = 2; i < n; i++)
191 {
192 t0 = Mpy_32_16 (b1_h, b1_l, x); // t0 = 2.0*x*b1
193 t0 = L_shl (t0, 1);
194 t0 = L_mac (t0, b2_h, (Word16) 0x8000); // t0 = 2.0*x*b1 - b2
195 t0 = L_msu (t0, b2_l, 1);
196 t0 = L_mac (t0, f[i], 8192); // t0 = 2.0*x*b1 - b2 + f[i]
197
198 L_Extract (t0, &b0_h, &b0_l); // b0 = 2.0*x*b1 - b2 + f[i]
199
200 b2_l = b1_l; // b2 = b1;
201 b2_h = b1_h;
202 b1_l = b0_l; // b1 = b0;
203 b1_h = b0_h;
204 }
205
206 t0 = Mpy_32_16 (b1_h, b1_l, x); // t0 = x*b1;
207 t0 = L_mac (t0, b2_h, (Word16) 0x8000); // t0 = x*b1 - b2
208 t0 = L_msu (t0, b2_l, 1);
209 t0 = L_mac (t0, f[i], 4096); // t0 = x*b1 - b2 + f[i]/2
210
211 t0 = L_shl (t0, 6);
212
213 cheb = extract_h (t0);
214
215 return (cheb);
216 }
217
218 ------------------------------------------------------------------------------
219 RESOURCES USED [optional]
220
221 When the code is written for a specific target processor the
222 the resources used should be documented below.
223
224 HEAP MEMORY USED: x bytes
225
226 STACK MEMORY USED: x bytes
227
228 CLOCK CYCLES: (cycle count equation for this function) + (variable
229 used to represent cycle count for each subroutine
230 called)
231 where: (cycle count variable) = cycle count for [subroutine
232 name]
233
234 ------------------------------------------------------------------------------
235 CAUTION [optional]
236 [State any special notes, constraints or cautions for users of this function]
237
238 ------------------------------------------------------------------------------
239 */
240 #ifdef __clang__
241 __attribute__((no_sanitize("integer")))
242 #endif
Chebps(Word16 x,Word16 f[],Word16 n,Flag * pOverflow)243 static Word16 Chebps(Word16 x,
244 Word16 f[], /* (n) */
245 Word16 n,
246 Flag *pOverflow)
247 {
248 Word16 i;
249 Word16 cheb;
250 Word16 b1_h;
251 Word16 b1_l;
252 Word32 t0;
253 Word32 L_temp;
254 Word16 *p_f = &f[1];
255
256 OSCL_UNUSED_ARG(pOverflow);
257
258 /* L_temp = 1.0 */
259
260 L_temp = 0x01000000L;
261
262 t0 = ((Word32) x << 10) + ((Word32) * (p_f++) << 14);
263
264 /* b1 = t0 = 2*x + f[1] */
265
266 b1_h = (Word16)(t0 >> 16);
267 b1_l = (Word16)((t0 >> 1) - (b1_h << 15));
268
269
270 for (i = 2; i < n; i++)
271 {
272 /* t0 = 2.0*x*b1 */
273 t0 = ((Word32) b1_h * x);
274 t0 += ((Word32) b1_l * x) >> 15;
275 t0 <<= 2;
276
277 /* t0 = 2.0*x*b1 - b2 */
278 t0 -= L_temp;
279
280 /* t0 = 2.0*x*b1 - b2 + f[i] */
281 t0 += (Word32) * (p_f++) << 14;
282
283 L_temp = ((Word32) b1_h << 16) + ((Word32) b1_l << 1);
284
285 /* b0 = 2.0*x*b1 - b2 + f[i]*/
286 b1_h = (Word16)(t0 >> 16);
287 b1_l = (Word16)((t0 >> 1) - (b1_h << 15));
288
289 }
290
291 /* t0 = x*b1; */
292 t0 = ((Word32) b1_h * x);
293 t0 += ((Word32) b1_l * x) >> 15;
294 t0 <<= 1;
295
296
297 /* t0 = x*b1 - b2 */
298 t0 -= L_temp;
299
300 /* t0 = x*b1 - b2 + f[i]/2 */
301 t0 += (Word32) * (p_f) << 13;
302
303
304 if ((UWord32)(t0 - 0xfe000000L) < (UWord32)0x03ffffffL)
305 {
306 cheb = (Word16)(t0 >> 10);
307 }
308 else
309 {
310 if (t0 > (Word32) 0x01ffffffL)
311 {
312 cheb = MAX_16;
313
314 }
315 else
316 {
317 cheb = MIN_16;
318 }
319 }
320
321 return (cheb);
322 }
323
324
325 /*
326 ------------------------------------------------------------------------------
327 FUNCTION NAME: Az_lsp
328 ------------------------------------------------------------------------------
329 INPUT AND OUTPUT DEFINITIONS FOR Az_lsp
330
331 Inputs:
332 a = predictor coefficients (Word16)
333 lsp = line spectral pairs (Word16)
334 old_lsp = old line spectral pairs (Word16)
335
336 pOverflow = pointer to overflow (Flag)
337
338 Outputs:
339 pOverflow -> 1 if the operations in the function resulted in saturation.
340
341 Returns:
342 None.
343
344 Global Variables Used:
345 None.
346
347 Local Variables Needed:
348 None.
349
350 ------------------------------------------------------------------------------
351 FUNCTION DESCRIPTION
352
353 This function computes the LSPs from the LP coefficients.
354
355 The sum and difference filters are computed and divided by 1+z^{-1} and
356 1-z^{-1}, respectively.
357
358 f1[i] = a[i] + a[11-i] - f1[i-1] ; i=1,...,5
359 f2[i] = a[i] - a[11-i] + f2[i-1] ; i=1,...,5
360
361 The roots of F1(z) and F2(z) are found using Chebyshev polynomial evaluation.
362 The polynomials are evaluated at 60 points regularly spaced in the
363 frequency domain. The sign change interval is subdivided 4 times to better
364 track the root. The LSPs are found in the cosine domain [1,-1].
365
366 If less than 10 roots are found, the LSPs from the past frame are used.
367
368 ------------------------------------------------------------------------------
369 REQUIREMENTS
370
371 None.
372
373 ------------------------------------------------------------------------------
374 REFERENCES
375
376 az_lsp.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001
377
378 ------------------------------------------------------------------------------
379 PSEUDO-CODE
380
381 void Az_lsp (
382 Word16 a[], // (i) : predictor coefficients (MP1)
383 Word16 lsp[], // (o) : line spectral pairs (M)
384 Word16 old_lsp[] // (i) : old lsp[] (in case not found 10 roots) (M)
385 )
386 {
387 Word16 i, j, nf, ip;
388 Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
389 Word16 x, y, sign, exp;
390 Word16 *coef;
391 Word16 f1[M / 2 + 1], f2[M / 2 + 1];
392 Word32 t0;
393
394 *-------------------------------------------------------------*
395 * find the sum and diff. pol. F1(z) and F2(z) *
396 * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) *
397 * *
398 * f1[0] = 1.0; *
399 * f2[0] = 1.0; *
400 * *
401 * for (i = 0; i< NC; i++) *
402 * { *
403 * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; *
404 * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; *
405 * } *
406 *-------------------------------------------------------------*
407
408 f1[0] = 1024; // f1[0] = 1.0
409 f2[0] = 1024; // f2[0] = 1.0
410
411 // The reference ETSI code uses a global flag for Overflow. However, in the
412 // actual implementation a pointer to Overflow flag is passed in as a
413 // parameter to the function. This pointer is passed into all the basic math
414 // functions invoked
415
416 for (i = 0; i < NC; i++)
417 {
418 t0 = L_mult (a[i + 1], 8192); // x = (a[i+1] + a[M-i]) >> 2
419 t0 = L_mac (t0, a[M - i], 8192);
420 x = extract_h (t0);
421 // f1[i+1] = a[i+1] + a[M-i] - f1[i]
422 f1[i + 1] = sub (x, f1[i]);
423
424 t0 = L_mult (a[i + 1], 8192); // x = (a[i+1] - a[M-i]) >> 2
425 t0 = L_msu (t0, a[M - i], 8192);
426 x = extract_h (t0);
427 // f2[i+1] = a[i+1] - a[M-i] + f2[i]
428 f2[i + 1] = add (x, f2[i]);
429 }
430
431 *-------------------------------------------------------------*
432 * find the LSPs using the Chebychev pol. evaluation *
433 *-------------------------------------------------------------*
434
435 nf = 0; // number of found frequencies
436 ip = 0; // indicator for f1 or f2
437
438 coef = f1;
439
440 xlow = grid[0];
441 ylow = Chebps (xlow, coef, NC);
442
443 j = 0;
444 // while ( (nf < M) && (j < grid_points) )
445 while ((sub (nf, M) < 0) && (sub (j, grid_points) < 0))
446 {
447 j++;
448 xhigh = xlow;
449 yhigh = ylow;
450 xlow = grid[j];
451 ylow = Chebps (xlow, coef, NC);
452
453 if (L_mult (ylow, yhigh) <= (Word32) 0L)
454 {
455
456 // divide 4 times the interval
457
458 for (i = 0; i < 4; i++)
459 {
460 // xmid = (xlow + xhigh)/2
461 xmid = add (shr (xlow, 1), shr (xhigh, 1));
462 ymid = Chebps (xmid, coef, NC);
463
464 if (L_mult (ylow, ymid) <= (Word32) 0L)
465 {
466 yhigh = ymid;
467 xhigh = xmid;
468 }
469 else
470 {
471 ylow = ymid;
472 xlow = xmid;
473 }
474 }
475
476 *-------------------------------------------------------------*
477 * Linear interpolation *
478 * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); *
479 *-------------------------------------------------------------*
480
481 x = sub (xhigh, xlow);
482 y = sub (yhigh, ylow);
483
484 if (y == 0)
485 {
486 xint = xlow;
487 }
488 else
489 {
490 sign = y;
491 y = abs_s (y);
492 exp = norm_s (y);
493 y = shl (y, exp);
494 y = div_s ((Word16) 16383, y);
495 t0 = L_mult (x, y);
496 t0 = L_shr (t0, sub (20, exp));
497 y = extract_l (t0); // y= (xhigh-xlow)/(yhigh-ylow)
498
499 if (sign < 0)
500 y = negate (y);
501
502 t0 = L_mult (ylow, y);
503 t0 = L_shr (t0, 11);
504 xint = sub (xlow, extract_l (t0)); // xint = xlow - ylow*y
505 }
506
507 lsp[nf] = xint;
508 xlow = xint;
509 nf++;
510
511 if (ip == 0)
512 {
513 ip = 1;
514 coef = f2;
515 }
516 else
517 {
518 ip = 0;
519 coef = f1;
520 }
521 ylow = Chebps (xlow, coef, NC);
522
523 }
524 }
525
526 // Check if M roots found
527
528 if (sub (nf, M) < 0)
529 {
530 for (i = 0; i < M; i++)
531 {
532 lsp[i] = old_lsp[i];
533 }
534
535 }
536 return;
537 }
538
539 ------------------------------------------------------------------------------
540 RESOURCES USED [optional]
541
542 When the code is written for a specific target processor the
543 the resources used should be documented below.
544
545 HEAP MEMORY USED: x bytes
546
547 STACK MEMORY USED: x bytes
548
549 CLOCK CYCLES: (cycle count equation for this function) + (variable
550 used to represent cycle count for each subroutine
551 called)
552 where: (cycle count variable) = cycle count for [subroutine
553 name]
554
555 ------------------------------------------------------------------------------
556 CAUTION [optional]
557 [State any special notes, constraints or cautions for users of this function]
558
559 ------------------------------------------------------------------------------
560 */
561
Az_lsp(Word16 a[],Word16 lsp[],Word16 old_lsp[],Flag * pOverflow)562 void Az_lsp(
563 Word16 a[], /* (i) : predictor coefficients (MP1) */
564 Word16 lsp[], /* (o) : line spectral pairs (M) */
565 Word16 old_lsp[], /* (i) : old lsp[] (in case not found 10 roots) (M) */
566 Flag *pOverflow /* (i/o): overflow flag */
567 )
568 {
569 Word16 i;
570 Word16 j;
571 Word16 nf;
572 Word16 ip;
573 Word16 xlow;
574 Word16 ylow;
575 Word16 xhigh;
576 Word16 yhigh;
577 Word16 xmid;
578 Word16 ymid;
579 Word16 xint;
580 Word16 x;
581 Word16 y;
582 Word16 sign;
583 Word16 exp;
584 Word16 *coef;
585 Word16 f1[NC + 1];
586 Word16 f2[NC + 1];
587 Word32 L_temp1;
588 Word32 L_temp2;
589 Word16 *p_f1 = f1;
590 Word16 *p_f2 = f2;
591
592 /*-------------------------------------------------------------*
593 * find the sum and diff. pol. F1(z) and F2(z) *
594 * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) *
595 * *
596 * f1[0] = 1.0; *
597 * f2[0] = 1.0; *
598 * *
599 * for (i = 0; i< NC; i++) *
600 * { *
601 * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; *
602 * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; *
603 * } *
604 *-------------------------------------------------------------*/
605
606 *p_f1 = 1024; /* f1[0] = 1.0 */
607 *p_f2 = 1024; /* f2[0] = 1.0 */
608
609 for (i = 0; i < NC; i++)
610 {
611 L_temp1 = (Word32) * (a + i + 1);
612 L_temp2 = (Word32) * (a + M - i);
613 /* x = (a[i+1] + a[M-i]) >> 2 */
614 x = (Word16)((L_temp1 + L_temp2) >> 2);
615 /* y = (a[i+1] - a[M-i]) >> 2 */
616 y = (Word16)((L_temp1 - L_temp2) >> 2);
617 /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */
618 x -= *(p_f1++);
619 *(p_f1) = x;
620 /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */
621 y += *(p_f2++);
622 *(p_f2) = y;
623 }
624
625 /*-------------------------------------------------------------*
626 * find the LSPs using the Chebychev pol. evaluation *
627 *-------------------------------------------------------------*/
628
629 nf = 0; /* number of found frequencies */
630 ip = 0; /* indicator for f1 or f2 */
631
632 coef = f1;
633
634 xlow = *(grid);
635 ylow = Chebps(xlow, coef, NC, pOverflow);
636
637 j = 0;
638
639 while ((nf < M) && (j < grid_points))
640 {
641 j++;
642 xhigh = xlow;
643 yhigh = ylow;
644 xlow = *(grid + j);
645 ylow = Chebps(xlow, coef, NC, pOverflow);
646
647 if (((Word32)ylow*yhigh) <= 0)
648 {
649 /* divide 4 times the interval */
650 for (i = 4; i != 0; i--)
651 {
652 /* xmid = (xlow + xhigh)/2 */
653 x = xlow >> 1;
654 y = xhigh >> 1;
655 xmid = x + y;
656
657 ymid = Chebps(xmid, coef, NC, pOverflow);
658
659 if (((Word32)ylow*ymid) <= 0)
660 {
661 yhigh = ymid;
662 xhigh = xmid;
663 }
664 else
665 {
666 ylow = ymid;
667 xlow = xmid;
668 }
669 }
670
671 /*-------------------------------------------------------------*
672 * Linear interpolation *
673 * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); *
674 *-------------------------------------------------------------*/
675
676 x = xhigh - xlow;
677 y = yhigh - ylow;
678
679 if (y == 0)
680 {
681 xint = xlow;
682 }
683 else
684 {
685 sign = y;
686 y = abs_s(y);
687 exp = norm_s(y);
688 y <<= exp;
689 y = div_s((Word16) 16383, y);
690
691 y = ((Word32)x * y) >> (19 - exp);
692
693 if (sign < 0)
694 {
695 y = -y;
696 }
697
698 /* xint = xlow - ylow*y */
699 xint = xlow - (((Word32) ylow * y) >> 10);
700 }
701
702 *(lsp + nf) = xint;
703 xlow = xint;
704 nf++;
705
706 if (ip == 0)
707 {
708 ip = 1;
709 coef = f2;
710 }
711 else
712 {
713 ip = 0;
714 coef = f1;
715 }
716
717 ylow = Chebps(xlow, coef, NC, pOverflow);
718
719 }
720 }
721
722 /* Check if M roots found */
723
724 if (nf < M)
725 {
726 for (i = NC; i != 0 ; i--)
727 {
728 *lsp++ = *old_lsp++;
729 *lsp++ = *old_lsp++;
730 }
731 }
732
733 }
734
Chebps_Wrapper(Word16 x,Word16 f[],Word16 n,Flag * pOverflow)735 Word16 Chebps_Wrapper(Word16 x,
736 Word16 f[], /* (n) */
737 Word16 n,
738 Flag *pOverflow)
739 {
740 return Chebps(x, f, n, pOverflow);
741 }
742
743