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