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
Chebps(Word16 x,Word16 f[],Word16 n,Flag * pOverflow)241 static Word16 Chebps(Word16 x,
242 Word16 f[], /* (n) */
243 Word16 n,
244 Flag *pOverflow)
245 {
246 Word16 i;
247 Word16 cheb;
248 Word16 b1_h;
249 Word16 b1_l;
250 Word32 t0;
251 Word32 L_temp;
252 Word16 *p_f = &f[1];
253
254 OSCL_UNUSED_ARG(pOverflow);
255
256 /* L_temp = 1.0 */
257
258 L_temp = 0x01000000L;
259
260 t0 = ((Word32) x << 10) + ((Word32) * (p_f++) << 14);
261
262 /* b1 = t0 = 2*x + f[1] */
263
264 b1_h = (Word16)(t0 >> 16);
265 b1_l = (Word16)((t0 >> 1) - (b1_h << 15));
266
267
268 for (i = 2; i < n; i++)
269 {
270 /* t0 = 2.0*x*b1 */
271 t0 = ((Word32) b1_h * x);
272 t0 += ((Word32) b1_l * x) >> 15;
273 t0 <<= 2;
274
275 /* t0 = 2.0*x*b1 - b2 */
276 t0 -= L_temp;
277
278 /* t0 = 2.0*x*b1 - b2 + f[i] */
279 t0 += (Word32) * (p_f++) << 14;
280
281 L_temp = ((Word32) b1_h << 16) + ((Word32) b1_l << 1);
282
283 /* b0 = 2.0*x*b1 - b2 + f[i]*/
284 b1_h = (Word16)(t0 >> 16);
285 b1_l = (Word16)((t0 >> 1) - (b1_h << 15));
286
287 }
288
289 /* t0 = x*b1; */
290 t0 = ((Word32) b1_h * x);
291 t0 += ((Word32) b1_l * x) >> 15;
292 t0 <<= 1;
293
294
295 /* t0 = x*b1 - b2 */
296 t0 -= L_temp;
297
298 /* t0 = x*b1 - b2 + f[i]/2 */
299 t0 += (Word32) * (p_f) << 13;
300
301
302 if ((UWord32)(t0 - 0xfe000000L) < (UWord32)0x03ffffffL)
303 {
304 cheb = (Word16)(t0 >> 10);
305 }
306 else
307 {
308 if (t0 > (Word32) 0x01ffffffL)
309 {
310 cheb = MAX_16;
311
312 }
313 else
314 {
315 cheb = MIN_16;
316 }
317 }
318
319 return (cheb);
320 }
321
322
323 /*
324 ------------------------------------------------------------------------------
325 FUNCTION NAME: Az_lsp
326 ------------------------------------------------------------------------------
327 INPUT AND OUTPUT DEFINITIONS FOR Az_lsp
328
329 Inputs:
330 a = predictor coefficients (Word16)
331 lsp = line spectral pairs (Word16)
332 old_lsp = old line spectral pairs (Word16)
333
334 pOverflow = pointer to overflow (Flag)
335
336 Outputs:
337 pOverflow -> 1 if the operations in the function resulted in saturation.
338
339 Returns:
340 None.
341
342 Global Variables Used:
343 None.
344
345 Local Variables Needed:
346 None.
347
348 ------------------------------------------------------------------------------
349 FUNCTION DESCRIPTION
350
351 This function computes the LSPs from the LP coefficients.
352
353 The sum and difference filters are computed and divided by 1+z^{-1} and
354 1-z^{-1}, respectively.
355
356 f1[i] = a[i] + a[11-i] - f1[i-1] ; i=1,...,5
357 f2[i] = a[i] - a[11-i] + f2[i-1] ; i=1,...,5
358
359 The roots of F1(z) and F2(z) are found using Chebyshev polynomial evaluation.
360 The polynomials are evaluated at 60 points regularly spaced in the
361 frequency domain. The sign change interval is subdivided 4 times to better
362 track the root. The LSPs are found in the cosine domain [1,-1].
363
364 If less than 10 roots are found, the LSPs from the past frame are used.
365
366 ------------------------------------------------------------------------------
367 REQUIREMENTS
368
369 None.
370
371 ------------------------------------------------------------------------------
372 REFERENCES
373
374 az_lsp.c, UMTS GSM AMR speech codec, R99 - Version 3.2.0, March 2, 2001
375
376 ------------------------------------------------------------------------------
377 PSEUDO-CODE
378
379 void Az_lsp (
380 Word16 a[], // (i) : predictor coefficients (MP1)
381 Word16 lsp[], // (o) : line spectral pairs (M)
382 Word16 old_lsp[] // (i) : old lsp[] (in case not found 10 roots) (M)
383 )
384 {
385 Word16 i, j, nf, ip;
386 Word16 xlow, ylow, xhigh, yhigh, xmid, ymid, xint;
387 Word16 x, y, sign, exp;
388 Word16 *coef;
389 Word16 f1[M / 2 + 1], f2[M / 2 + 1];
390 Word32 t0;
391
392 *-------------------------------------------------------------*
393 * find the sum and diff. pol. F1(z) and F2(z) *
394 * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) *
395 * *
396 * f1[0] = 1.0; *
397 * f2[0] = 1.0; *
398 * *
399 * for (i = 0; i< NC; i++) *
400 * { *
401 * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; *
402 * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; *
403 * } *
404 *-------------------------------------------------------------*
405
406 f1[0] = 1024; // f1[0] = 1.0
407 f2[0] = 1024; // f2[0] = 1.0
408
409 // The reference ETSI code uses a global flag for Overflow. However, in the
410 // actual implementation a pointer to Overflow flag is passed in as a
411 // parameter to the function. This pointer is passed into all the basic math
412 // functions invoked
413
414 for (i = 0; i < NC; i++)
415 {
416 t0 = L_mult (a[i + 1], 8192); // x = (a[i+1] + a[M-i]) >> 2
417 t0 = L_mac (t0, a[M - i], 8192);
418 x = extract_h (t0);
419 // f1[i+1] = a[i+1] + a[M-i] - f1[i]
420 f1[i + 1] = sub (x, f1[i]);
421
422 t0 = L_mult (a[i + 1], 8192); // x = (a[i+1] - a[M-i]) >> 2
423 t0 = L_msu (t0, a[M - i], 8192);
424 x = extract_h (t0);
425 // f2[i+1] = a[i+1] - a[M-i] + f2[i]
426 f2[i + 1] = add (x, f2[i]);
427 }
428
429 *-------------------------------------------------------------*
430 * find the LSPs using the Chebychev pol. evaluation *
431 *-------------------------------------------------------------*
432
433 nf = 0; // number of found frequencies
434 ip = 0; // indicator for f1 or f2
435
436 coef = f1;
437
438 xlow = grid[0];
439 ylow = Chebps (xlow, coef, NC);
440
441 j = 0;
442 // while ( (nf < M) && (j < grid_points) )
443 while ((sub (nf, M) < 0) && (sub (j, grid_points) < 0))
444 {
445 j++;
446 xhigh = xlow;
447 yhigh = ylow;
448 xlow = grid[j];
449 ylow = Chebps (xlow, coef, NC);
450
451 if (L_mult (ylow, yhigh) <= (Word32) 0L)
452 {
453
454 // divide 4 times the interval
455
456 for (i = 0; i < 4; i++)
457 {
458 // xmid = (xlow + xhigh)/2
459 xmid = add (shr (xlow, 1), shr (xhigh, 1));
460 ymid = Chebps (xmid, coef, NC);
461
462 if (L_mult (ylow, ymid) <= (Word32) 0L)
463 {
464 yhigh = ymid;
465 xhigh = xmid;
466 }
467 else
468 {
469 ylow = ymid;
470 xlow = xmid;
471 }
472 }
473
474 *-------------------------------------------------------------*
475 * Linear interpolation *
476 * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); *
477 *-------------------------------------------------------------*
478
479 x = sub (xhigh, xlow);
480 y = sub (yhigh, ylow);
481
482 if (y == 0)
483 {
484 xint = xlow;
485 }
486 else
487 {
488 sign = y;
489 y = abs_s (y);
490 exp = norm_s (y);
491 y = shl (y, exp);
492 y = div_s ((Word16) 16383, y);
493 t0 = L_mult (x, y);
494 t0 = L_shr (t0, sub (20, exp));
495 y = extract_l (t0); // y= (xhigh-xlow)/(yhigh-ylow)
496
497 if (sign < 0)
498 y = negate (y);
499
500 t0 = L_mult (ylow, y);
501 t0 = L_shr (t0, 11);
502 xint = sub (xlow, extract_l (t0)); // xint = xlow - ylow*y
503 }
504
505 lsp[nf] = xint;
506 xlow = xint;
507 nf++;
508
509 if (ip == 0)
510 {
511 ip = 1;
512 coef = f2;
513 }
514 else
515 {
516 ip = 0;
517 coef = f1;
518 }
519 ylow = Chebps (xlow, coef, NC);
520
521 }
522 }
523
524 // Check if M roots found
525
526 if (sub (nf, M) < 0)
527 {
528 for (i = 0; i < M; i++)
529 {
530 lsp[i] = old_lsp[i];
531 }
532
533 }
534 return;
535 }
536
537 ------------------------------------------------------------------------------
538 RESOURCES USED [optional]
539
540 When the code is written for a specific target processor the
541 the resources used should be documented below.
542
543 HEAP MEMORY USED: x bytes
544
545 STACK MEMORY USED: x bytes
546
547 CLOCK CYCLES: (cycle count equation for this function) + (variable
548 used to represent cycle count for each subroutine
549 called)
550 where: (cycle count variable) = cycle count for [subroutine
551 name]
552
553 ------------------------------------------------------------------------------
554 CAUTION [optional]
555 [State any special notes, constraints or cautions for users of this function]
556
557 ------------------------------------------------------------------------------
558 */
559
Az_lsp(Word16 a[],Word16 lsp[],Word16 old_lsp[],Flag * pOverflow)560 void Az_lsp(
561 Word16 a[], /* (i) : predictor coefficients (MP1) */
562 Word16 lsp[], /* (o) : line spectral pairs (M) */
563 Word16 old_lsp[], /* (i) : old lsp[] (in case not found 10 roots) (M) */
564 Flag *pOverflow /* (i/o): overflow flag */
565 )
566 {
567 register Word16 i;
568 register Word16 j;
569 register Word16 nf;
570 register Word16 ip;
571 Word16 xlow;
572 Word16 ylow;
573 Word16 xhigh;
574 Word16 yhigh;
575 Word16 xmid;
576 Word16 ymid;
577 Word16 xint;
578 Word16 x;
579 Word16 y;
580 Word16 sign;
581 Word16 exp;
582 Word16 *coef;
583 Word16 f1[NC + 1];
584 Word16 f2[NC + 1];
585 Word32 L_temp1;
586 Word32 L_temp2;
587 Word16 *p_f1 = f1;
588 Word16 *p_f2 = f2;
589
590 /*-------------------------------------------------------------*
591 * find the sum and diff. pol. F1(z) and F2(z) *
592 * F1(z) <--- F1(z)/(1+z**-1) & F2(z) <--- F2(z)/(1-z**-1) *
593 * *
594 * f1[0] = 1.0; *
595 * f2[0] = 1.0; *
596 * *
597 * for (i = 0; i< NC; i++) *
598 * { *
599 * f1[i+1] = a[i+1] + a[M-i] - f1[i] ; *
600 * f2[i+1] = a[i+1] - a[M-i] + f2[i] ; *
601 * } *
602 *-------------------------------------------------------------*/
603
604 *p_f1 = 1024; /* f1[0] = 1.0 */
605 *p_f2 = 1024; /* f2[0] = 1.0 */
606
607 for (i = 0; i < NC; i++)
608 {
609 L_temp1 = (Word32) * (a + i + 1);
610 L_temp2 = (Word32) * (a + M - i);
611 /* x = (a[i+1] + a[M-i]) >> 2 */
612 x = (Word16)((L_temp1 + L_temp2) >> 2);
613 /* y = (a[i+1] - a[M-i]) >> 2 */
614 y = (Word16)((L_temp1 - L_temp2) >> 2);
615 /* f1[i+1] = a[i+1] + a[M-i] - f1[i] */
616 x -= *(p_f1++);
617 *(p_f1) = x;
618 /* f2[i+1] = a[i+1] - a[M-i] + f2[i] */
619 y += *(p_f2++);
620 *(p_f2) = y;
621 }
622
623 /*-------------------------------------------------------------*
624 * find the LSPs using the Chebychev pol. evaluation *
625 *-------------------------------------------------------------*/
626
627 nf = 0; /* number of found frequencies */
628 ip = 0; /* indicator for f1 or f2 */
629
630 coef = f1;
631
632 xlow = *(grid);
633 ylow = Chebps(xlow, coef, NC, pOverflow);
634
635 j = 0;
636
637 while ((nf < M) && (j < grid_points))
638 {
639 j++;
640 xhigh = xlow;
641 yhigh = ylow;
642 xlow = *(grid + j);
643 ylow = Chebps(xlow, coef, NC, pOverflow);
644
645 if (((Word32)ylow*yhigh) <= 0)
646 {
647 /* divide 4 times the interval */
648 for (i = 4; i != 0; i--)
649 {
650 /* xmid = (xlow + xhigh)/2 */
651 x = xlow >> 1;
652 y = xhigh >> 1;
653 xmid = x + y;
654
655 ymid = Chebps(xmid, coef, NC, pOverflow);
656
657 if (((Word32)ylow*ymid) <= 0)
658 {
659 yhigh = ymid;
660 xhigh = xmid;
661 }
662 else
663 {
664 ylow = ymid;
665 xlow = xmid;
666 }
667 }
668
669 /*-------------------------------------------------------------*
670 * Linear interpolation *
671 * xint = xlow - ylow*(xhigh-xlow)/(yhigh-ylow); *
672 *-------------------------------------------------------------*/
673
674 x = xhigh - xlow;
675 y = yhigh - ylow;
676
677 if (y == 0)
678 {
679 xint = xlow;
680 }
681 else
682 {
683 sign = y;
684 y = abs_s(y);
685 exp = norm_s(y);
686 y <<= exp;
687 y = div_s((Word16) 16383, y);
688
689 y = ((Word32)x * y) >> (19 - exp);
690
691 if (sign < 0)
692 {
693 y = -y;
694 }
695
696 /* xint = xlow - ylow*y */
697 xint = xlow - (((Word32) ylow * y) >> 10);
698 }
699
700 *(lsp + nf) = xint;
701 xlow = xint;
702 nf++;
703
704 if (ip == 0)
705 {
706 ip = 1;
707 coef = f2;
708 }
709 else
710 {
711 ip = 0;
712 coef = f1;
713 }
714
715 ylow = Chebps(xlow, coef, NC, pOverflow);
716
717 }
718 }
719
720 /* Check if M roots found */
721
722 if (nf < M)
723 {
724 for (i = NC; i != 0 ; i--)
725 {
726 *lsp++ = *old_lsp++;
727 *lsp++ = *old_lsp++;
728 }
729 }
730
731 }
732
Chebps_Wrapper(Word16 x,Word16 f[],Word16 n,Flag * pOverflow)733 Word16 Chebps_Wrapper(Word16 x,
734 Word16 f[], /* (n) */
735 Word16 n,
736 Flag *pOverflow)
737 {
738 return Chebps(x, f, n, pOverflow);
739 }
740
741