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 20 21------------------------------------------------------------------------------ 22 REVISION HISTORY 23 24 Description: Modified from original shareware code 25 26 Description: Modified to remove instances of pow() and sqrt(), and 27 optimized for inclusion in fixed-point version of decoder. 28 29 Description: Modified to include comments/optimizations from code review. 30 Also, declared appropriate variables as type "const" 31 32 Description: Adopted strategy of "one q-format per sfb" strategy, which 33 eliminated the array q-format from this function. The q-format the 34 random vector is stored in is now returned from the function. 35 36 Description: Completely redesigned the routine to allow a simplified 37 calculation of the adjusted noise, by eliminating the dependency 38 on the band_length. Added polynomial approximation for the 39 function 1/sqrt(power). Updated comments and pseudo-code 40 41 Description: Modified function description, pseudocode, etc. 42 43 Description: 44 Modified casting to ensure proper operations for different platforms 45 46 Description: 47 Eliminiated access to memory for noise seed. Now a local variable is 48 used. Also unrolled loops to speed up code. 49 50 Description: 51 Modified pointer decrement to a pointer increment, to ensure proper 52 compiler behavior 53 54 Description: 55------------------------------------------------------------------------------ 56 INPUT AND OUTPUT DEFINITIONS 57 58 Inputs: random_array[] = Array for storage of the power-scaled 59 random values of length "band_length" 60 Int32 61 62 band_length = Length of random_array[] 63 const Int 64 65 pSeed = seed for random number generator 66 Int32* 67 68 power_scale = scale factor for this particular band 69 const Int 70 71 72 Local Stores/Buffers/Pointers Needed: 73 None 74 75 Global Stores/Buffers/Pointers Needed: 76 None 77 78 Outputs: Function returns the q-format the random vector is stored in. 79 80 Pointers and Buffers Modified: 81 random_array[] = filled with random numbers scaled 82 to the correct power as defined by the input value power_scale. 83 84 Local Stores Modified: 85 None 86 87 Global Stores Modified: 88 None 89 90------------------------------------------------------------------------------ 91 FUNCTION DESCRIPTION 92 93 This function generates a vector of uniformly distributed random numbers for 94 the PNS block. The random numbers are each scaled by a scale_factor, 95 defined in Ref(2) as 96 97 2^(scale_factor/4) 98 ------------------ 99 sqrt(N*MEAN_NRG) 100 101 where N == band_length, and MEAN_NRG is defined as... 102 103 N-1 104 ___ 105 1 \ 106 --- > x(i)^2 107 N /__ 108 i=0 109 110 And x is the unscaled vector from the random number generator. 111 112 This function takes advantage of the fact that the portion of the 113 scale_factor that is divisible by 4 can be simply accounted for by varying 114 the q-format. 115 116 The scaling of the random numbers is thus broken into the 117 equivalent equation below. 118 119 2^(scale_factor%4) 2^(floor(scale_factor/4)) 120 ------------------ * 121 sqrt(N*MEAN_NRG) 122 123 124 2^(scale_factor%4) is stored in a simple 4-element table. 125 2^(floor(scale_factor/4) is accounted for by adjusting the q-format. 126 sqrt(N*MEAN_NRG) is calculated and implemented via a polynomial approximation. 127------------------------------------------------------------------------------ 128 REQUIREMENTS 129 130 This function shall produce uniformly distributed random 32-bit integers, 131 with signed random values of average energy equal to the results of the ISO 132 code's multiplying factor discussed in the FUNCTION DESCRIPTION section. 133 134 Please see Ref (2) for a detailed description of the requirements. 135------------------------------------------------------------------------------ 136 REFERENCES 137 138 (1) Numerical Recipes in C Second Edition 139 William H. Press Saul A. Teukolsky 140 William T. Vetterling Brian P. Flannery 141 Page 284 142 143 (2) ISO/IEC 14496-3:1999(E) 144 Part 3 145 Subpart 4.6.12 (Perceptual Noise Substitution) 146 147 (3) MPEG-2 NBC Audio Decoder 148 "This software module was originally developed by AT&T, Dolby 149 Laboratories, Fraunhofer Gesellschaft IIS in the course of development 150 of the MPEG-2 NBC/MPEG-4 Audio standard ISO/IEC 13818-7, 14496-1,2 and 151 3. This software module is an implementation of a part of one or more 152 MPEG-2 NBC/MPEG-4 Audio tools as specified by the MPEG-2 NBC/MPEG-4 153 Audio standard. ISO/IEC gives users of the MPEG-2 NBC/MPEG-4 Audio 154 standards free license to this software module or modifications thereof 155 for use in hardware or software products claiming conformance to the 156 MPEG-2 NBC/MPEG-4 Audio standards. Those intending to use this software 157 module in hardware or software products are advised that this use may 158 infringe existing patents. The original developer of this software 159 module and his/her company, the subsequent editors and their companies, 160 and ISO/IEC have no liability for use of this software module or 161 modifications thereof in an implementation. Copyright is not released 162 for non MPEG-2 NBC/MPEG-4 Audio conforming products.The original 163 developer retains full right to use the code for his/her own purpose, 164 assign or donate the code to a third party and to inhibit third party 165 from using the code for non MPEG-2 NBC/MPEG-4 Audio conforming products. 166 This copyright notice must be included in all copies or derivative 167 works." 168 Copyright(c)1996. 169 170------------------------------------------------------------------------------ 171 PSEUDO-CODE 172 173 power_adj = scale_mod_4[power_scale & 3]; 174 175 power = 0; 176 177 FOR (k=band_length; k > 0; k--) 178 179 *(pSeed) = *(pSeed) * 1664525L; 180 *(pSeed) = *(pSeed) + 1013904223L; 181 182 temp = (Int)(*(pSeed) >> 16); 183 184 power = power + ((temp*temp) >> 6); 185 186 *(pArray) = (Int32)temp; 187 188 pArray = pArray + 1; 189 190 ENDFOR 191 192 k = 0; 193 q_adjust = 30; 194 195 IF (power) 196 THEN 197 198 WHILE ( power > 32767) 199 200 power = power >> 1; 201 k = k + 1; 202 203 ENDWHILE 204 205 k = k - 13; 206 207 IF (k < 0) 208 THEN 209 k = -k; 210 IF ( k & 1 ) 211 THEN 212 power_adj = (power_adj*SQRT_OF_2)>>14; 213 ENDIF 214 q_adjust = q_adjust - ( k >> 1); 215 216 ELSE IF (k > 0) 217 THEN 218 IF ( k & 1 ) 219 THEN 220 power_adj = (power_adj*INV_SQRT_OF_2)>>14; 221 ENDIF 222 q_adjust = q_adjust + ( k >> 1); 223 ENDIF 224 225 pInvSqrtCoeff = inv_sqrt_coeff; 226 227 inv_sqrt_power = (*(pInvSqrtCoeff)* power) >>15; 228 229 pInvSqrtCoeff = pInvSqrtCoeff + 1; 230 231 inv_sqrt_power = inv_sqrt_power + *(pInvSqrtCoeff); 232 233 pInvSqrtCoeff = pInvSqrtCoeff + 1; 234 235 FOR ( k=INV_SQRT_POLY_ORDER - 1; k>0; k--) 236 237 inv_sqrt_power = ( inv_sqrt_power * power)>>15; 238 239 inv_sqrt_power = inv_sqrt_power + *(pInvSqrtCoeff); 240 241 pInvSqrtCoeff = pInvSqrtCoeff + 1; 242 243 ENDFOR 244 245 inv_sqrt_power = (inv_sqrt_power*power_adj)>>13; 246 247 FOR (k=band_length; k > 0; k--) 248 249 pArray = pArray - 1; 250 251 *(pArray) = *(pArray)*inv_sqrt_power; 252 253 ENDFOR 254 255 ENDIF 256 257 q_adjust = q_adjust - (power_scale >> 2); 258 259 return q_adjust; 260 261------------------------------------------------------------------------------ 262 RESOURCES USED 263 When the code is written for a specific target processor 264 the resources used should be documented below. 265 266 STACK USAGE: [stack count for this module] + [variable to represent 267 stack usage for each subroutine called] 268 269 where: [stack usage variable] = stack usage for [subroutine 270 name] (see [filename].ext) 271 272 DATA MEMORY USED: x words 273 274 PROGRAM MEMORY USED: x words 275 276 CLOCK CYCLES: [cycle count equation for this module] + [variable 277 used to represent cycle count for each subroutine 278 called] 279 280 where: [cycle count variable] = cycle count for [subroutine 281 name] (see [filename].ext) 282 283------------------------------------------------------------------------------ 284*/ 285 286 287/*---------------------------------------------------------------------------- 288; INCLUDES 289----------------------------------------------------------------------------*/ 290#include "pv_audio_type_defs.h" 291#include "gen_rand_vector.h" 292#include "window_block_fxp.h" 293 294/*---------------------------------------------------------------------------- 295; MACROS 296; Define module specific macros here 297----------------------------------------------------------------------------*/ 298 299/*---------------------------------------------------------------------------- 300; DEFINES 301; Include all pre-processor statements here. Include conditional 302; compile variables also. 303----------------------------------------------------------------------------*/ 304#define SQRT_OF_2 23170 /* sqrt(2) in Q14 */ 305#define INV_SQRT_OF_2 11585 /* 1/sqrt(2) in Q14 */ 306#define INV_SQRT_POLY_ORDER 4 307 308/*---------------------------------------------------------------------------- 309; LOCAL FUNCTION DEFINITIONS 310; Function Prototype declaration 311----------------------------------------------------------------------------*/ 312 313/*---------------------------------------------------------------------------- 314; LOCAL STORE/BUFFER/POINTER DEFINITIONS 315; Variable declaration - defined here and used outside this module 316----------------------------------------------------------------------------*/ 317 318 319 320/* 321 * 2^([0:3]/4) = 1.0000 1.1892 1.4142 1.6818 322 */ 323const UInt scale_mod_4[4] = { 16384, 19484, 23170, 27554}; 324 325/* 326 * polynomial approx. in Q12 (type Int) 327 */ 328 329const Int inv_sqrt_coeff[INV_SQRT_POLY_ORDER+1] = 330 { 4680, -17935, 27697, -22326, 11980}; 331 332 333/*---------------------------------------------------------------------------- 334; EXTERNAL FUNCTION REFERENCES 335; Declare functions defined elsewhere and referenced in this module 336----------------------------------------------------------------------------*/ 337 338/*---------------------------------------------------------------------------- 339; EXTERNAL GLOBAL STORE/BUFFER/POINTER REFERENCES 340; Declare variables used in this module but defined elsewhere 341----------------------------------------------------------------------------*/ 342 343/*---------------------------------------------------------------------------- 344; FUNCTION CODE 345----------------------------------------------------------------------------*/ 346 347Int gen_rand_vector( 348 Int32 random_array[], 349 const Int band_length, 350 Int32* pSeed, 351 const Int power_scale) 352{ 353 354 Int k; 355 UInt power_adj; 356 Int q_adjust = 30; 357 358 Int32 temp; 359 Int32 seed; 360 Int32 power; 361 362 Int32* pArray = &random_array[0]; 363 364 Int32 inv_sqrt_power; 365 const Int *pInvSqrtCoeff; 366 367 /* 368 * The out of the random number generator is scaled is such a way 369 * that is independent of the band length. 370 * The output is computed as: 371 * 372 * x(i) 373 * output = ------------------ * 2^(power_scale%4) 2^(floor(power_scale/4)) 374 * bl 375 * sqrt( SUM x(i)^2 ) 376 * 0 377 * 378 * bl == band length 379 */ 380 381 382 /* 383 * get 2^(power_scale%4) 384 */ 385 386 387 power = 0; 388 389 seed = *pSeed; 390 391 /* 392 * band_length is always an even number (check tables in pg.66 IS0 14496-3) 393 */ 394 if (band_length < 0 || band_length > LONG_WINDOW) 395 { 396 return q_adjust; /* avoid any processing on error condition */ 397 } 398 399 for (k = (band_length >> 1); k != 0; k--) 400 { 401 /*------------------------------------------------ 402 Numerical Recipes in C 403 Page 284 404 ------------------------------------------------*/ 405 seed *= 1664525L; 406 seed += 1013904223L; 407 408 temp = seed >> 16; 409 410 seed *= 1664525L; 411 seed += 1013904223L; 412 413 /* shift by 6 make room for band length accumulation */ 414 power += ((temp * temp) >> 6); 415 *pArray++ = temp; 416 417 temp = seed >> 16; 418 power += ((temp * temp) >> 6); 419 *pArray++ = temp; 420 421 } /* END for (k=half_band_length; k > 0; k--) */ 422 423 424 *pSeed = seed; 425 426 /* 427 * If the distribution is uniform, the power is expected to use between 428 * 28 and 27 bits, by shifting down by 13 bits the power will be a 429 * Q15 number. 430 * For different band lengths, the power uses between 20 and 29 bits 431 */ 432 433 434 k = 0; 435 436 if (power) 437 { 438 /* 439 * approximation requires power between 0.5 < power < 1 in Q15. 440 */ 441 442 while (power > 32767) 443 { 444 power >>= 1; 445 k++; 446 } 447 448 /* 449 * expected power bit usage == 27 bits 450 */ 451 452 k -= 13; 453 454 power_adj = scale_mod_4[power_scale & 3]; 455 456 if (k < 0) 457 { 458 k = -k; 459 if (k & 1) 460 { /* multiply by sqrt(2) */ 461 power_adj = (UInt)(((UInt32) power_adj * SQRT_OF_2) >> 14); 462 } 463 q_adjust -= (k >> 1); /* adjust Q instead of shifting up */ 464 } 465 else if (k > 0) 466 { 467 if (k & 1) 468 { /* multiply by 1/sqrt(2) */ 469 power_adj = (UInt)(((UInt32) power_adj * INV_SQRT_OF_2) >> 14); 470 } 471 q_adjust += (k >> 1); /* adjust Q instead of shifting down */ 472 } 473 474 /* 475 * Compute 1/sqrt(power), where 0.5 < power < 1.0 is approximated 476 * using a polynomial order INV_SQRT_POLY_ORDER 477 */ 478 479 pInvSqrtCoeff = inv_sqrt_coeff; 480 481 inv_sqrt_power = (*(pInvSqrtCoeff++) * power) >> 15; 482 inv_sqrt_power += *(pInvSqrtCoeff++); 483 inv_sqrt_power = (inv_sqrt_power * power) >> 15; 484 inv_sqrt_power += *(pInvSqrtCoeff++); 485 inv_sqrt_power = (inv_sqrt_power * power) >> 15; 486 inv_sqrt_power += *(pInvSqrtCoeff++); 487 inv_sqrt_power = (inv_sqrt_power * power) >> 15; 488 inv_sqrt_power += *(pInvSqrtCoeff); 489 490 inv_sqrt_power = (inv_sqrt_power * power_adj) >> 13; 491 492 pArray = &random_array[0]; 493 494 for (k = (band_length >> 1); k != 0; k--) 495 { 496 temp = *(pArray) * inv_sqrt_power; 497 *(pArray++) = temp; 498 temp = *(pArray) * inv_sqrt_power; 499 *(pArray++) = temp; 500 } /* END for (k=half_band_length; k > 0; k--) */ 501 502 } /* if(power) */ 503 504 /* 505 * Adjust Q with the value corresponding to 2^(floor(power_scale/4)) 506 */ 507 508 q_adjust -= (power_scale >> 2); 509 510 return (q_adjust); 511 512} /* gen_rand_vector */ 513