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.173
22 ANSI-C code for the Adaptive Multi-Rate - Wideband (AMR-WB) speech codec
23 Available from http://www.3gpp.org
24
25 (C) 2007, 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 This file contains mathematic operations in fixed point.
32
33 mult_int16_r() : Same as mult_int16 with rounding
34 shr_rnd() : Same as shr(var1,var2) but with rounding
35 div_16by16() : fractional integer division
36 one_ov_sqrt() : Compute 1/sqrt(L_x)
37 one_ov_sqrt_norm() : Compute 1/sqrt(x)
38 power_of_2() : power of 2
39 Dot_product12() : Compute scalar product of <x[],y[]> using accumulator
40 Isqrt() : inverse square root (16 bits precision).
41 amrwb_log_2() : log2 (16 bits precision).
42
43 These operations are not standard double precision operations.
44 They are used where low complexity is important and the full 32 bits
45 precision is not necessary. For example, the function Div_32() has a
46 24 bits precision which is enough for our purposes.
47
48 In this file, the values use theses representations:
49
50 int32 L_32 : standard signed 32 bits format
51 int16 hi, lo : L_32 = hi<<16 + lo<<1 (DPF - Double Precision Format)
52 int32 frac, int16 exp : L_32 = frac << exp-31 (normalised format)
53 int16 int, frac : L_32 = int.frac (fractional format)
54 ----------------------------------------------------------------------------*/
55
56 #include "pv_amr_wb_type_defs.h"
57 #include "pvamrwbdecoder_basic_op.h"
58 #include "pvamrwb_math_op.h"
59
60
61 /*----------------------------------------------------------------------------
62
63 Function Name : mult_int16_r
64
65 Purpose :
66
67 Same as mult_int16 with rounding, i.e.:
68 mult_int16_r(var1,var2) = extract_l(L_shr(((var1 * var2) + 16384),15)) and
69 mult_int16_r(-32768,-32768) = 32767.
70
71 Complexity weight : 2
72
73 Inputs :
74
75 var1
76 16 bit short signed integer (int16) whose value falls in the
77 range : 0xffff 8000 <= var1 <= 0x0000 7fff.
78
79 var2
80 16 bit short signed integer (int16) whose value falls in the
81 range : 0xffff 8000 <= var1 <= 0x0000 7fff.
82
83 Outputs :
84
85 none
86
87 Return Value :
88
89 var_out
90 16 bit short signed integer (int16) whose value falls in the
91 range : 0xffff 8000 <= var_out <= 0x0000 7fff.
92 ----------------------------------------------------------------------------*/
93
mult_int16_r(int16 var1,int16 var2)94 int16 mult_int16_r(int16 var1, int16 var2)
95 {
96 int32 L_product_arr;
97
98 L_product_arr = (int32) var1 * (int32) var2; /* product */
99 L_product_arr += (int32) 0x00004000L; /* round */
100 L_product_arr >>= 15; /* shift */
101 if ((L_product_arr >> 15) != (L_product_arr >> 31))
102 {
103 L_product_arr = (L_product_arr >> 31) ^ MAX_16;
104 }
105
106 return ((int16)L_product_arr);
107 }
108
109
110
111 /*----------------------------------------------------------------------------
112
113 Function Name : shr_rnd
114
115 Purpose :
116
117 Same as shr(var1,var2) but with rounding. Saturate the result in case of|
118 underflows or overflows :
119 - If var2 is greater than zero :
120 if (sub(shl_int16(shr(var1,var2),1),shr(var1,sub(var2,1))))
121 is equal to zero
122 then
123 shr_rnd(var1,var2) = shr(var1,var2)
124 else
125 shr_rnd(var1,var2) = add_int16(shr(var1,var2),1)
126 - If var2 is less than or equal to zero :
127 shr_rnd(var1,var2) = shr(var1,var2).
128
129 Complexity weight : 2
130
131 Inputs :
132
133 var1
134 16 bit short signed integer (int16) whose value falls in the
135 range : 0xffff 8000 <= var1 <= 0x0000 7fff.
136
137 var2
138 16 bit short signed integer (int16) whose value falls in the
139 range : 0x0000 0000 <= var2 <= 0x0000 7fff.
140
141 Outputs :
142
143 none
144
145 Return Value :
146
147 var_out
148 16 bit short signed integer (int16) whose value falls in the
149 range : 0xffff 8000 <= var_out <= 0x0000 7fff.
150 ----------------------------------------------------------------------------*/
151
shr_rnd(int16 var1,int16 var2)152 int16 shr_rnd(int16 var1, int16 var2)
153 {
154 int16 var_out;
155
156 var_out = (int16)(var1 >> (var2 & 0xf));
157 if (var2)
158 {
159 if ((var1 & ((int16) 1 << (var2 - 1))) != 0)
160 {
161 var_out++;
162 }
163 }
164 return (var_out);
165 }
166
167
168 /*----------------------------------------------------------------------------
169
170 Function Name : div_16by16
171
172 Purpose :
173
174 Produces a result which is the fractional integer division of var1 by
175 var2; var1 and var2 must be positive and var2 must be greater or equal
176 to var1; the result is positive (leading bit equal to 0) and truncated
177 to 16 bits.
178 If var1 = var2 then div(var1,var2) = 32767.
179
180 Complexity weight : 18
181
182 Inputs :
183
184 var1
185 16 bit short signed integer (int16) whose value falls in the
186 range : 0x0000 0000 <= var1 <= var2 and var2 != 0.
187
188 var2
189 16 bit short signed integer (int16) whose value falls in the
190 range : var1 <= var2 <= 0x0000 7fff and var2 != 0.
191
192 Outputs :
193
194 none
195
196 Return Value :
197
198 var_out
199 16 bit short signed integer (int16) whose value falls in the
200 range : 0x0000 0000 <= var_out <= 0x0000 7fff.
201 It's a Q15 value (point between b15 and b14).
202 ----------------------------------------------------------------------------*/
203
div_16by16(int16 var1,int16 var2)204 int16 div_16by16(int16 var1, int16 var2)
205 {
206
207 int16 var_out = 0;
208 int16 iteration;
209 int32 L_num;
210 int32 L_denom;
211 int32 L_denom_by_2;
212 int32 L_denom_by_4;
213
214 if ((var1 > var2) || (var1 < 0))
215 {
216 return 0; // used to exit(0);
217 }
218 if (var1)
219 {
220 if (var1 != var2)
221 {
222
223 L_num = (int32) var1;
224 L_denom = (int32) var2;
225 L_denom_by_2 = (L_denom << 1);
226 L_denom_by_4 = (L_denom << 2);
227 for (iteration = 5; iteration > 0; iteration--)
228 {
229 var_out <<= 3;
230 L_num <<= 3;
231
232 if (L_num >= L_denom_by_4)
233 {
234 L_num -= L_denom_by_4;
235 var_out |= 4;
236 }
237
238 if (L_num >= L_denom_by_2)
239 {
240 L_num -= L_denom_by_2;
241 var_out |= 2;
242 }
243
244 if (L_num >= (L_denom))
245 {
246 L_num -= (L_denom);
247 var_out |= 1;
248 }
249
250 }
251 }
252 else
253 {
254 var_out = MAX_16;
255 }
256 }
257
258 return (var_out);
259
260 }
261
262
263
264 /*----------------------------------------------------------------------------
265
266 Function Name : one_ov_sqrt
267
268 Compute 1/sqrt(L_x).
269 if L_x is negative or zero, result is 1 (7fffffff).
270
271 Algorithm:
272
273 1- Normalization of L_x.
274 2- call Isqrt_n(L_x, exponant)
275 3- L_y = L_x << exponant
276 ----------------------------------------------------------------------------*/
one_ov_sqrt(int32 L_x)277 int32 one_ov_sqrt( /* (o) Q31 : output value (range: 0<=val<1) */
278 int32 L_x /* (i) Q0 : input value (range: 0<=val<=7fffffff) */
279 )
280 {
281 int16 exp;
282 int32 L_y;
283
284 exp = normalize_amr_wb(L_x);
285 L_x <<= exp; /* L_x is normalized */
286 exp = 31 - exp;
287
288 one_ov_sqrt_norm(&L_x, &exp);
289
290 L_y = shl_int32(L_x, exp); /* denormalization */
291
292 return (L_y);
293 }
294
295 /*----------------------------------------------------------------------------
296
297 Function Name : one_ov_sqrt_norm
298
299 Compute 1/sqrt(value).
300 if value is negative or zero, result is 1 (frac=7fffffff, exp=0).
301
302 Algorithm:
303
304 The function 1/sqrt(value) is approximated by a table and linear
305 interpolation.
306
307 1- If exponant is odd then shift fraction right once.
308 2- exponant = -((exponant-1)>>1)
309 3- i = bit25-b30 of fraction, 16 <= i <= 63 ->because of normalization.
310 4- a = bit10-b24
311 5- i -=16
312 6- fraction = table[i]<<16 - (table[i] - table[i+1]) * a * 2
313 ----------------------------------------------------------------------------*/
314 static const int16 table_isqrt[49] =
315 {
316 32767, 31790, 30894, 30070, 29309, 28602, 27945, 27330, 26755, 26214,
317 25705, 25225, 24770, 24339, 23930, 23541, 23170, 22817, 22479, 22155,
318 21845, 21548, 21263, 20988, 20724, 20470, 20225, 19988, 19760, 19539,
319 19326, 19119, 18919, 18725, 18536, 18354, 18176, 18004, 17837, 17674,
320 17515, 17361, 17211, 17064, 16921, 16782, 16646, 16514, 16384
321 };
322
one_ov_sqrt_norm(int32 * frac,int16 * exp)323 void one_ov_sqrt_norm(
324 int32 * frac, /* (i/o) Q31: normalized value (1.0 < frac <= 0.5) */
325 int16 * exp /* (i/o) : exponent (value = frac x 2^exponent) */
326 )
327 {
328 int16 i, a, tmp;
329
330
331 if (*frac <= (int32) 0)
332 {
333 *exp = 0;
334 *frac = 0x7fffffffL;
335 return;
336 }
337
338 if ((*exp & 1) == 1) /* If exponant odd -> shift right */
339 *frac >>= 1;
340
341 *exp = negate_int16((*exp - 1) >> 1);
342
343 *frac >>= 9;
344 i = extract_h(*frac); /* Extract b25-b31 */
345 *frac >>= 1;
346 a = (int16)(*frac); /* Extract b10-b24 */
347 a = (int16)(a & (int16) 0x7fff);
348
349 i -= 16;
350
351 *frac = L_deposit_h(table_isqrt[i]); /* table[i] << 16 */
352 tmp = table_isqrt[i] - table_isqrt[i + 1]; /* table[i] - table[i+1]) */
353
354 *frac = msu_16by16_from_int32(*frac, tmp, a); /* frac -= tmp*a*2 */
355
356 return;
357 }
358
359 /*----------------------------------------------------------------------------
360
361 Function Name : power_2()
362
363 L_x = pow(2.0, exponant.fraction) (exponant = interger part)
364 = pow(2.0, 0.fraction) << exponant
365
366 Algorithm:
367
368 The function power_2(L_x) is approximated by a table and linear
369 interpolation.
370
371 1- i = bit10-b15 of fraction, 0 <= i <= 31
372 2- a = bit0-b9 of fraction
373 3- L_x = table[i]<<16 - (table[i] - table[i+1]) * a * 2
374 4- L_x = L_x >> (30-exponant) (with rounding)
375 ----------------------------------------------------------------------------*/
376 const int16 table_pow2[33] =
377 {
378 16384, 16743, 17109, 17484, 17867, 18258, 18658, 19066, 19484, 19911,
379 20347, 20792, 21247, 21713, 22188, 22674, 23170, 23678, 24196, 24726,
380 25268, 25821, 26386, 26964, 27554, 28158, 28774, 29405, 30048, 30706,
381 31379, 32066, 32767
382 };
383
power_of_2(int16 exponant,int16 fraction)384 int32 power_of_2( /* (o) Q0 : result (range: 0<=val<=0x7fffffff) */
385 int16 exponant, /* (i) Q0 : Integer part. (range: 0<=val<=30) */
386 int16 fraction /* (i) Q15 : Fractionnal part. (range: 0.0<=val<1.0) */
387 )
388 {
389 int16 exp, i, a, tmp;
390 int32 L_x;
391
392 L_x = fraction << 5; /* L_x = fraction<<6 */
393 i = (fraction >> 10); /* Extract b10-b16 of fraction */
394 a = (int16)(L_x); /* Extract b0-b9 of fraction */
395 a = (int16)(a & (int16) 0x7fff);
396
397 L_x = ((int32)table_pow2[i]) << 15; /* table[i] << 16 */
398 tmp = table_pow2[i] - table_pow2[i + 1]; /* table[i] - table[i+1] */
399 L_x -= ((int32)tmp * a); /* L_x -= tmp*a*2 */
400
401 exp = 29 - exponant ;
402
403 if (exp)
404 {
405 L_x = ((L_x >> exp) + ((L_x >> (exp - 1)) & 1));
406 }
407
408 return (L_x);
409 }
410
411 /*----------------------------------------------------------------------------
412 *
413 * Function Name : Dot_product12()
414 *
415 * Compute scalar product of <x[],y[]> using accumulator.
416 *
417 * The result is normalized (in Q31) with exponent (0..30).
418 *
419 * Algorithm:
420 *
421 * dot_product = sum(x[i]*y[i]) i=0..N-1
422 ----------------------------------------------------------------------------*/
423
Dot_product12(int16 x[],int16 y[],int16 lg,int16 * exp)424 int32 Dot_product12( /* (o) Q31: normalized result (1 < val <= -1) */
425 int16 x[], /* (i) 12bits: x vector */
426 int16 y[], /* (i) 12bits: y vector */
427 int16 lg, /* (i) : vector length */
428 int16 * exp /* (o) : exponent of result (0..+30) */
429 )
430 {
431 int16 i, sft;
432 int32 L_sum;
433 int16 *pt_x = x;
434 int16 *pt_y = y;
435
436 L_sum = 1L;
437
438
439 for (i = lg >> 3; i != 0; i--)
440 {
441 L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++));
442 L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++));
443 L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++));
444 L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++));
445 L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++));
446 L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++));
447 L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++));
448 L_sum = mac_16by16_to_int32(L_sum, *(pt_x++), *(pt_y++));
449 }
450
451 /* Normalize acc in Q31 */
452
453 sft = normalize_amr_wb(L_sum);
454 L_sum <<= sft;
455
456 *exp = 30 - sft; /* exponent = 0..30 */
457
458 return (L_sum);
459 }
460
461 /* Table for Log2() */
462 const int16 Log2_norm_table[33] =
463 {
464 0, 1455, 2866, 4236, 5568, 6863, 8124, 9352, 10549, 11716,
465 12855, 13967, 15054, 16117, 17156, 18172, 19167, 20142, 21097, 22033,
466 22951, 23852, 24735, 25603, 26455, 27291, 28113, 28922, 29716, 30497,
467 31266, 32023, 32767
468 };
469
470 /*----------------------------------------------------------------------------
471 *
472 * FUNCTION: Lg2_normalized()
473 *
474 * PURPOSE: Computes log2(L_x, exp), where L_x is positive and
475 * normalized, and exp is the normalisation exponent
476 * If L_x is negative or zero, the result is 0.
477 *
478 * DESCRIPTION:
479 * The function Log2(L_x) is approximated by a table and linear
480 * interpolation. The following steps are used to compute Log2(L_x)
481 *
482 * 1- exponent = 30-norm_exponent
483 * 2- i = bit25-b31 of L_x; 32<=i<=63 (because of normalization).
484 * 3- a = bit10-b24
485 * 4- i -=32
486 * 5- fraction = table[i]<<16 - (table[i] - table[i+1]) * a * 2
487 *
488 ----------------------------------------------------------------------------*/
Lg2_normalized(int32 L_x,int16 exp,int16 * exponent,int16 * fraction)489 void Lg2_normalized(
490 int32 L_x, /* (i) : input value (normalized) */
491 int16 exp, /* (i) : norm_l (L_x) */
492 int16 *exponent, /* (o) : Integer part of Log2. (range: 0<=val<=30) */
493 int16 *fraction /* (o) : Fractional part of Log2. (range: 0<=val<1) */
494 )
495 {
496 int16 i, a, tmp;
497 int32 L_y;
498
499 if (L_x <= (int32) 0)
500 {
501 *exponent = 0;
502 *fraction = 0;;
503 return;
504 }
505
506 *exponent = 30 - exp;
507
508 L_x >>= 9;
509 i = extract_h(L_x); /* Extract b25-b31 */
510 L_x >>= 1;
511 a = (int16)(L_x); /* Extract b10-b24 of fraction */
512 a &= 0x7fff;
513
514 i -= 32;
515
516 L_y = L_deposit_h(Log2_norm_table[i]); /* table[i] << 16 */
517 tmp = Log2_norm_table[i] - Log2_norm_table[i + 1]; /* table[i] - table[i+1] */
518 L_y = msu_16by16_from_int32(L_y, tmp, a); /* L_y -= tmp*a*2 */
519
520 *fraction = extract_h(L_y);
521
522 return;
523 }
524
525
526
527 /*----------------------------------------------------------------------------
528 *
529 * FUNCTION: amrwb_log_2()
530 *
531 * PURPOSE: Computes log2(L_x), where L_x is positive.
532 * If L_x is negative or zero, the result is 0.
533 *
534 * DESCRIPTION:
535 * normalizes L_x and then calls Lg2_normalized().
536 *
537 ----------------------------------------------------------------------------*/
amrwb_log_2(int32 L_x,int16 * exponent,int16 * fraction)538 void amrwb_log_2(
539 int32 L_x, /* (i) : input value */
540 int16 *exponent, /* (o) : Integer part of Log2. (range: 0<=val<=30) */
541 int16 *fraction /* (o) : Fractional part of Log2. (range: 0<=val<1) */
542 )
543 {
544 int16 exp;
545
546 exp = normalize_amr_wb(L_x);
547 Lg2_normalized(shl_int32(L_x, exp), exp, exponent, fraction);
548 }
549
550
551 /*****************************************************************************
552 *
553 * These operations are not standard double precision operations. *
554 * They are used where single precision is not enough but the full 32 bits *
555 * precision is not necessary. For example, the function Div_32() has a *
556 * 24 bits precision which is enough for our purposes. *
557 * *
558 * The double precision numbers use a special representation: *
559 * *
560 * L_32 = hi<<16 + lo<<1 *
561 * *
562 * L_32 is a 32 bit integer. *
563 * hi and lo are 16 bit signed integers. *
564 * As the low part also contains the sign, this allows fast multiplication. *
565 * *
566 * 0x8000 0000 <= L_32 <= 0x7fff fffe. *
567 * *
568 * We will use DPF (Double Precision Format )in this file to specify *
569 * this special format. *
570 *****************************************************************************
571 */
572
573
574 /*----------------------------------------------------------------------------
575 *
576 * Function int32_to_dpf()
577 *
578 * Extract from a 32 bit integer two 16 bit DPF.
579 *
580 * Arguments:
581 *
582 * L_32 : 32 bit integer.
583 * 0x8000 0000 <= L_32 <= 0x7fff ffff.
584 * hi : b16 to b31 of L_32
585 * lo : (L_32 - hi<<16)>>1
586 *
587 ----------------------------------------------------------------------------*/
588
int32_to_dpf(int32 L_32,int16 * hi,int16 * lo)589 void int32_to_dpf(int32 L_32, int16 *hi, int16 *lo)
590 {
591 *hi = (int16)(L_32 >> 16);
592 *lo = (int16)((L_32 - (*hi << 16)) >> 1);
593 return;
594 }
595
596
597 /*----------------------------------------------------------------------------
598 * Function mpy_dpf_32()
599 *
600 * Multiply two 32 bit integers (DPF). The result is divided by 2**31
601 *
602 * L_32 = (hi1*hi2)<<1 + ( (hi1*lo2)>>15 + (lo1*hi2)>>15 )<<1
603 *
604 * This operation can also be viewed as the multiplication of two Q31
605 * number and the result is also in Q31.
606 *
607 * Arguments:
608 *
609 * hi1 hi part of first number
610 * lo1 lo part of first number
611 * hi2 hi part of second number
612 * lo2 lo part of second number
613 *
614 ----------------------------------------------------------------------------*/
615
mpy_dpf_32(int16 hi1,int16 lo1,int16 hi2,int16 lo2)616 int32 mpy_dpf_32(int16 hi1, int16 lo1, int16 hi2, int16 lo2)
617 {
618 int32 L_32;
619
620 L_32 = mul_16by16_to_int32(hi1, hi2);
621 L_32 = mac_16by16_to_int32(L_32, mult_int16(hi1, lo2), 1);
622 L_32 = mac_16by16_to_int32(L_32, mult_int16(lo1, hi2), 1);
623
624 return (L_32);
625 }
626
627
628