1 /* 2 * psymodel.c 3 * 4 * Copyright (c) 1999 Mark Taylor 5 * 6 * This library is free software; you can redistribute it and/or 7 * modify it under the terms of the GNU Library General Public 8 * License as published by the Free Software Foundation; either 9 * version 2 of the License, or (at your option) any later version. 10 * 11 * This library is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 14 * Library General Public License for more details. 15 * 16 * You should have received a copy of the GNU Library General Public 17 * License along with this library; if not, write to the 18 * Free Software Foundation, Inc., 59 Temple Place - Suite 330, 19 * Boston, MA 02111-1307, USA. 20 */ 21 22 /* $Id: psymodel.c,v 1.75 2001/03/25 21:37:00 shibatch Exp $ */ 23 24 25 /* 26 PSYCHO ACOUSTICS 27 28 29 This routine computes the psycho acoustics, delayed by 30 one granule. 31 32 Input: buffer of PCM data (1024 samples). 33 34 This window should be centered over the 576 sample granule window. 35 The routine will compute the psycho acoustics for 36 this granule, but return the psycho acoustics computed 37 for the *previous* granule. This is because the block 38 type of the previous granule can only be determined 39 after we have computed the psycho acoustics for the following 40 granule. 41 42 Output: maskings and energies for each scalefactor band. 43 block type, PE, and some correlation measures. 44 The PE is used by CBR modes to determine if extra bits 45 from the bit reservoir should be used. The correlation 46 measures are used to determine mid/side or regular stereo. 47 48 49 Notation: 50 51 barks: a non-linear frequency scale. Mapping from frequency to 52 barks is given by freq2bark() 53 54 scalefactor bands: The spectrum (frequencies) are broken into 55 SBMAX "scalefactor bands". Thes bands 56 are determined by the MPEG ISO spec. In 57 the noise shaping/quantization code, we allocate 58 bits among the partition bands to achieve the 59 best possible quality 60 61 partition bands: The spectrum is also broken into about 62 64 "partition bands". Each partition 63 band is about .34 barks wide. There are about 2-5 64 partition bands for each scalefactor band. 65 66 LAME computes all psycho acoustic information for each partition 67 band. Then at the end of the computations, this information 68 is mapped to scalefactor bands. The energy in each scalefactor 69 band is taken as the sum of the energy in all partition bands 70 which overlap the scalefactor band. The maskings can be computed 71 in the same way (and thus represent the average masking in that band) 72 or by taking the minmum value multiplied by the number of 73 partition bands used (which represents a minimum masking in that band). 74 75 76 The general outline is as follows: 77 78 79 1. compute the energy in each partition band 80 2. compute the tonality in each partition band 81 3. compute the strength of each partion band "masker" 82 4. compute the masking (via the spreading function applied to each masker) 83 5. Modifications for mid/side masking. 84 85 Each partition band is considiered a "masker". The strength 86 of the i'th masker in band j is given by: 87 88 s3(bark(i)-bark(j))*strength(i) 89 90 The strength of the masker is a function of the energy and tonality. 91 The more tonal, the less masking. LAME uses a simple linear formula 92 (controlled by NMT and TMN) which says the strength is given by the 93 energy divided by a linear function of the tonality. 94 95 96 s3() is the "spreading function". It is given by a formula 97 determined via listening tests. 98 99 The total masking in the j'th partition band is the sum over 100 all maskings i. It is thus given by the convolution of 101 the strength with s3(), the "spreading function." 102 103 masking(j) = sum_over_i s3(i-j)*strength(i) = s3 o strength 104 105 where "o" = convolution operator. s3 is given by a formula determined 106 via listening tests. It is normalized so that s3 o 1 = 1. 107 108 Note: instead of a simple convolution, LAME also has the 109 option of using "additive masking" 110 111 The most critical part is step 2, computing the tonality of each 112 partition band. LAME has two tonality estimators. The first 113 is based on the ISO spec, and measures how predictiable the 114 signal is over time. The more predictable, the more tonal. 115 The second measure is based on looking at the spectrum of 116 a single granule. The more peaky the spectrum, the more 117 tonal. By most indications, the latter approach is better. 118 119 Finally, in step 5, the maskings for the mid and side 120 channel are possibly increased. Under certain circumstances, 121 noise in the mid & side channels is assumed to also 122 be masked by strong maskers in the L or R channels. 123 124 125 Other data computed by the psy-model: 126 127 ms_ratio side-channel / mid-channel masking ratio (for previous granule) 128 ms_ratio_next side-channel / mid-channel masking ratio for this granule 129 130 percep_entropy[2] L and R values (prev granule) of PE - A measure of how 131 much pre-echo is in the previous granule 132 percep_entropy_MS[2] mid and side channel values (prev granule) of percep_entropy 133 energy[4] L,R,M,S energy in each channel, prev granule 134 blocktype_d[2] block type to use for previous granule 135 136 137 */ 138 139 140 141 142 #ifdef HAVE_CONFIG_H 143 # include <config.h> 144 #endif 145 146 #include "util.h" 147 #include "encoder.h" 148 #include "psymodel.h" 149 #include "l3side.h" 150 #include <assert.h> 151 #include "tables.h" 152 #include "fft.h" 153 154 #ifdef WITH_DMALLOC 155 #include <dmalloc.h> 156 #endif 157 158 #define NSFIRLEN 21 159 #define rpelev 2 160 #define rpelev2 16 161 162 /* size of each partition band, in barks: */ 163 #define DELBARK .34 164 165 166 #if 1 167 /* AAC values, results in more masking over MP3 values */ 168 # define TMN 18 169 # define NMT 6 170 #else 171 /* MP3 values */ 172 # define TMN 29 173 # define NMT 6 174 #endif 175 176 #define NBPSY_l (SBMAX_l) 177 #define NBPSY_s (SBMAX_s) 178 179 180 #ifdef M_LN10 181 #define LN_TO_LOG10 (M_LN10/10) 182 #else 183 #define LN_TO_LOG10 0.2302585093 184 #endif 185 186 187 188 189 190 /* 191 L3psycho_anal. Compute psycho acoustics. 192 193 Data returned to the calling program must be delayed by one 194 granule. 195 196 This is done in two places. 197 If we do not need to know the blocktype, the copying 198 can be done here at the top of the program: we copy the data for 199 the last granule (computed during the last call) before it is 200 overwritten with the new data. It looks like this: 201 202 0. static psymodel_data 203 1. calling_program_data = psymodel_data 204 2. compute psymodel_data 205 206 For data which needs to know the blocktype, the copying must be 207 done at the end of this loop, and the old values must be saved: 208 209 0. static psymodel_data_old 210 1. compute psymodel_data 211 2. compute possible block type of this granule 212 3. compute final block type of previous granule based on #2. 213 4. calling_program_data = psymodel_data_old 214 5. psymodel_data_old = psymodel_data 215 */ 216 int L3psycho_anal( lame_global_flags * gfp, 217 const sample_t *buffer[2], int gr_out, 218 FLOAT8 *ms_ratio, 219 FLOAT8 *ms_ratio_next, 220 III_psy_ratio masking_ratio[2][2], 221 III_psy_ratio masking_MS_ratio[2][2], 222 FLOAT8 percep_entropy[2],FLOAT8 percep_MS_entropy[2], 223 FLOAT8 energy[4], 224 int blocktype_d[2]) 225 { 226 /* to get a good cache performance, one has to think about 227 * the sequence, in which the variables are used. 228 * (Note: these static variables have been moved to the gfc-> struct, 229 * and their order in memory is layed out in util.h) 230 */ 231 lame_internal_flags *gfc=gfp->internal_flags; 232 233 234 /* fft and energy calculation */ 235 FLOAT (*wsamp_l)[BLKSIZE]; 236 FLOAT (*wsamp_s)[3][BLKSIZE_s]; 237 238 /* convolution */ 239 FLOAT8 eb[CBANDS]; 240 FLOAT8 cb[CBANDS]; 241 FLOAT8 thr[CBANDS]; 242 243 /* ratios */ 244 FLOAT8 ms_ratio_l=0,ms_ratio_s=0; 245 246 /* block type */ 247 int blocktype[2],uselongblock[2]; 248 249 /* usual variables like loop indices, etc.. */ 250 int numchn, chn; 251 int b, i, j, k; 252 int sb,sblock; 253 254 255 if(gfc->psymodel_init==0) { 256 psymodel_init(gfp); 257 init_fft(gfc); 258 gfc->psymodel_init=1; 259 } 260 261 262 263 264 265 numchn = gfc->channels_out; 266 /* chn=2 and 3 = Mid and Side channels */ 267 if (gfp->mode == JOINT_STEREO) numchn=4; 268 269 for (chn=0; chn<numchn; chn++) { 270 for (i=0; i<numchn; ++i) { 271 energy[chn]=gfc->tot_ener[chn]; 272 } 273 } 274 275 for (chn=0; chn<numchn; chn++) { 276 277 /* there is a one granule delay. Copy maskings computed last call 278 * into masking_ratio to return to calling program. 279 */ 280 if (chn < 2) { 281 /* LR maskings */ 282 percep_entropy [chn] = gfc -> pe [chn]; 283 masking_ratio [gr_out] [chn] .en = gfc -> en [chn]; 284 masking_ratio [gr_out] [chn] .thm = gfc -> thm [chn]; 285 } else { 286 /* MS maskings */ 287 percep_MS_entropy [chn-2] = gfc -> pe [chn]; 288 masking_MS_ratio [gr_out] [chn-2].en = gfc -> en [chn]; 289 masking_MS_ratio [gr_out] [chn-2].thm = gfc -> thm [chn]; 290 } 291 292 293 294 /********************************************************************** 295 * compute FFTs 296 **********************************************************************/ 297 wsamp_s = gfc->wsamp_S+(chn & 1); 298 wsamp_l = gfc->wsamp_L+(chn & 1); 299 if (chn<2) { 300 fft_long ( gfc, *wsamp_l, chn, buffer); 301 fft_short( gfc, *wsamp_s, chn, buffer); 302 } 303 /* FFT data for mid and side channel is derived from L & R */ 304 if (chn == 2) 305 { 306 for (j = BLKSIZE-1; j >=0 ; --j) 307 { 308 FLOAT l = gfc->wsamp_L[0][j]; 309 FLOAT r = gfc->wsamp_L[1][j]; 310 gfc->wsamp_L[0][j] = (l+r)*(FLOAT)(SQRT2*0.5); 311 gfc->wsamp_L[1][j] = (l-r)*(FLOAT)(SQRT2*0.5); 312 } 313 for (b = 2; b >= 0; --b) 314 { 315 for (j = BLKSIZE_s-1; j >= 0 ; --j) 316 { 317 FLOAT l = gfc->wsamp_S[0][b][j]; 318 FLOAT r = gfc->wsamp_S[1][b][j]; 319 gfc->wsamp_S[0][b][j] = (l+r)*(FLOAT)(SQRT2*0.5); 320 gfc->wsamp_S[1][b][j] = (l-r)*(FLOAT)(SQRT2*0.5); 321 } 322 } 323 } 324 325 326 /********************************************************************** 327 * compute energies 328 **********************************************************************/ 329 330 331 332 gfc->energy[0] = (*wsamp_l)[0]; 333 gfc->energy[0] *= gfc->energy[0]; 334 335 gfc->tot_ener[chn] = gfc->energy[0]; /* sum total energy at nearly no extra cost */ 336 337 for (j=BLKSIZE/2-1; j >= 0; --j) 338 { 339 FLOAT re = (*wsamp_l)[BLKSIZE/2-j]; 340 FLOAT im = (*wsamp_l)[BLKSIZE/2+j]; 341 gfc->energy[BLKSIZE/2-j] = (re * re + im * im) * (FLOAT)0.5; 342 343 if (BLKSIZE/2-j > 10) 344 gfc->tot_ener[chn] += gfc->energy[BLKSIZE/2-j]; 345 } 346 for (b = 2; b >= 0; --b) 347 { 348 gfc->energy_s[b][0] = (*wsamp_s)[b][0]; 349 gfc->energy_s[b][0] *= gfc->energy_s [b][0]; 350 for (j=BLKSIZE_s/2-1; j >= 0; --j) 351 { 352 FLOAT re = (*wsamp_s)[b][BLKSIZE_s/2-j]; 353 FLOAT im = (*wsamp_s)[b][BLKSIZE_s/2+j]; 354 gfc->energy_s[b][BLKSIZE_s/2-j] = (re * re + im * im) * (FLOAT)0.5; 355 } 356 } 357 358 359 if (gfp->analysis) { 360 for (j=0; j<HBLKSIZE ; j++) { 361 gfc->pinfo->energy[gr_out][chn][j]=gfc->energy_save[chn][j]; 362 gfc->energy_save[chn][j]=gfc->energy[j]; 363 } 364 } 365 366 /********************************************************************** 367 * compute unpredicatability of first six spectral lines * 368 **********************************************************************/ 369 for ( j = 0; j < gfc->cw_lower_index; j++ ) 370 { /* calculate unpredictability measure cw */ 371 FLOAT an, a1, a2; 372 FLOAT bn, b1, b2; 373 FLOAT rn, r1, r2; 374 FLOAT numre, numim, den; 375 376 a2 = gfc-> ax_sav[chn][1][j]; 377 b2 = gfc-> bx_sav[chn][1][j]; 378 r2 = gfc-> rx_sav[chn][1][j]; 379 a1 = gfc-> ax_sav[chn][1][j] = gfc-> ax_sav[chn][0][j]; 380 b1 = gfc-> bx_sav[chn][1][j] = gfc-> bx_sav[chn][0][j]; 381 r1 = gfc-> rx_sav[chn][1][j] = gfc-> rx_sav[chn][0][j]; 382 an = gfc-> ax_sav[chn][0][j] = (*wsamp_l)[j]; 383 bn = gfc-> bx_sav[chn][0][j] = j==0 ? (*wsamp_l)[0] : (*wsamp_l)[BLKSIZE-j]; 384 rn = gfc-> rx_sav[chn][0][j] = sqrt(gfc->energy[j]); 385 386 { /* square (x1,y1) */ 387 if( r1 != 0 ) { 388 numre = (a1*b1); 389 numim = (a1*a1-b1*b1)*(FLOAT)0.5; 390 den = r1*r1; 391 } else { 392 numre = 1; 393 numim = 0; 394 den = 1; 395 } 396 } 397 398 { /* multiply by (x2,-y2) */ 399 if( r2 != 0 ) { 400 FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5; 401 FLOAT tmp1 = -a2*numre+tmp2; 402 numre = -b2*numim+tmp2; 403 numim = tmp1; 404 den *= r2; 405 } else { 406 /* do nothing */ 407 } 408 } 409 410 { /* r-prime factor */ 411 FLOAT tmp = (2*r1-r2)/den; 412 numre *= tmp; 413 numim *= tmp; 414 } 415 den=rn+fabs(2*r1-r2); 416 if( den != 0 ) { 417 numre = (an+bn)*(FLOAT)0.5-numre; 418 numim = (an-bn)*(FLOAT)0.5-numim; 419 den = sqrt(numre*numre+numim*numim)/den; 420 } 421 gfc->cw[j] = den; 422 } 423 424 425 426 /********************************************************************** 427 * compute unpredicatibility of next 200 spectral lines * 428 **********************************************************************/ 429 for ( j = gfc->cw_lower_index; j < gfc->cw_upper_index; j += 4 ) 430 {/* calculate unpredictability measure cw */ 431 FLOAT rn, r1, r2; 432 FLOAT numre, numim, den; 433 434 k = (j+2) / 4; 435 436 { /* square (x1,y1) */ 437 r1 = gfc->energy_s[0][k]; 438 if( r1 != 0 ) { 439 FLOAT a1 = (*wsamp_s)[0][k]; 440 FLOAT b1 = (*wsamp_s)[0][BLKSIZE_s-k]; /* k is never 0 */ 441 numre = (a1*b1); 442 numim = (a1*a1-b1*b1)*(FLOAT)0.5; 443 den = r1; 444 r1 = sqrt(r1); 445 } else { 446 numre = 1; 447 numim = 0; 448 den = 1; 449 } 450 } 451 452 453 { /* multiply by (x2,-y2) */ 454 r2 = gfc->energy_s[2][k]; 455 if( r2 != 0 ) { 456 FLOAT a2 = (*wsamp_s)[2][k]; 457 FLOAT b2 = (*wsamp_s)[2][BLKSIZE_s-k]; 458 459 460 FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5; 461 FLOAT tmp1 = -a2*numre+tmp2; 462 numre = -b2*numim+tmp2; 463 numim = tmp1; 464 465 r2 = sqrt(r2); 466 den *= r2; 467 } else { 468 /* do nothing */ 469 } 470 } 471 472 { /* r-prime factor */ 473 FLOAT tmp = (2*r1-r2)/den; 474 numre *= tmp; 475 numim *= tmp; 476 } 477 478 rn = sqrt(gfc->energy_s[1][k]); 479 den=rn+fabs(2*r1-r2); 480 if( den != 0 ) { 481 FLOAT an = (*wsamp_s)[1][k]; 482 FLOAT bn = (*wsamp_s)[1][BLKSIZE_s-k]; 483 numre = (an+bn)*(FLOAT)0.5-numre; 484 numim = (an-bn)*(FLOAT)0.5-numim; 485 den = sqrt(numre*numre+numim*numim)/den; 486 } 487 488 gfc->cw[j+1] = gfc->cw[j+2] = gfc->cw[j+3] = gfc->cw[j] = den; 489 } 490 491 #if 0 492 for ( j = 14; j < HBLKSIZE-4; j += 4 ) 493 {/* calculate energy from short ffts */ 494 FLOAT8 tot,ave; 495 k = (j+2) / 4; 496 for (tot=0, sblock=0; sblock < 3; sblock++) 497 tot+=gfc->energy_s[sblock][k]; 498 ave = gfc->energy[j+1]+ gfc->energy[j+2]+ gfc->energy[j+3]+ gfc->energy[j]; 499 ave /= 4.; 500 gfc->energy[j+1] = gfc->energy[j+2] = gfc->energy[j+3] = gfc->energy[j]=tot; 501 } 502 #endif 503 504 /********************************************************************** 505 * Calculate the energy and the unpredictability in the threshold * 506 * calculation partitions * 507 **********************************************************************/ 508 509 b = 0; 510 for (j = 0; j < gfc->cw_upper_index && gfc->numlines_l[b] && b<gfc->npart_l_orig; ) 511 { 512 FLOAT8 ebb, cbb; 513 514 ebb = gfc->energy[j]; 515 cbb = gfc->energy[j] * gfc->cw[j]; 516 j++; 517 518 for (i = gfc->numlines_l[b] - 1; i > 0; i--) 519 { 520 ebb += gfc->energy[j]; 521 cbb += gfc->energy[j] * gfc->cw[j]; 522 j++; 523 } 524 eb[b] = ebb; 525 cb[b] = cbb; 526 b++; 527 } 528 529 for (; b < gfc->npart_l_orig; b++ ) 530 { 531 FLOAT8 ebb = gfc->energy[j++]; 532 assert(gfc->numlines_l[b]); 533 for (i = gfc->numlines_l[b] - 1; i > 0; i--) 534 { 535 ebb += gfc->energy[j++]; 536 } 537 eb[b] = ebb; 538 cb[b] = ebb * 0.4; 539 } 540 541 /********************************************************************** 542 * convolve the partitioned energy and unpredictability * 543 * with the spreading function, s3_l[b][k] * 544 ******************************************************************** */ 545 gfc->pe[chn] = 0; /* calculate percetual entropy */ 546 for ( b = 0;b < gfc->npart_l; b++ ) 547 { 548 FLOAT8 tbb,ecb,ctb; 549 550 ecb = 0; 551 ctb = 0; 552 553 for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ ) 554 { 555 ecb += gfc->s3_l[b][k] * eb[k]; /* sprdngf for Layer III */ 556 ctb += gfc->s3_l[b][k] * cb[k]; 557 } 558 559 /* calculate the tonality of each threshold calculation partition 560 * calculate the SNR in each threshold calculation partition 561 * tonality = -0.299 - .43*log(ctb/ecb); 562 * tonality = 0: use NMT (lots of masking) 563 * tonality = 1: use TMN (little masking) 564 */ 565 566 /* ISO values */ 567 #define CONV1 (-.299) 568 #define CONV2 (-.43) 569 570 tbb = ecb; 571 if (tbb != 0) 572 { 573 tbb = ctb / tbb; 574 if (tbb <= exp((1-CONV1)/CONV2)) 575 { /* tonality near 1 */ 576 tbb = exp(-LN_TO_LOG10 * TMN); 577 } 578 else if (tbb >= exp((0-CONV1)/CONV2)) 579 {/* tonality near 0 */ 580 tbb = exp(-LN_TO_LOG10 * NMT); 581 } 582 else 583 { 584 /* convert to tonality index */ 585 /* tonality small: tbb=1 */ 586 /* tonality large: tbb=-.299 */ 587 tbb = CONV1 + CONV2*log(tbb); 588 tbb = exp(-LN_TO_LOG10 * ( TMN*tbb + (1-tbb)*NMT) ); 589 } 590 } 591 592 /* at this point, tbb represents the amount the spreading function 593 * will be reduced. The smaller the value, the less masking. 594 * minval[] = 1 (0db) says just use tbb. 595 * minval[]= .01 (-20db) says reduce spreading function by 596 * at least 20db. 597 */ 598 599 tbb = Min(gfc->minval[b], tbb); 600 ecb *= tbb; 601 602 /* long block pre-echo control. */ 603 /* rpelev=2.0, rpelev2=16.0 */ 604 /* note: all surges in PE are because of this pre-echo formula 605 * for thr[b]. If it this is not used, PE is always around 600 606 */ 607 /* dont use long block pre-echo control if previous granule was 608 * a short block. This is to avoid the situation: 609 * frame0: quiet (very low masking) 610 * frame1: surge (triggers short blocks) 611 * frame2: regular frame. looks like pre-echo when compared to 612 * frame0, but all pre-echo was in frame1. 613 */ 614 /* chn=0,1 L and R channels 615 chn=2,3 S and M channels. 616 */ 617 618 if (vbr_mtrh == gfp->VBR) { 619 thr[b] = Min (rpelev*gfc->nb_1[chn][b], rpelev2*gfc->nb_2[chn][b]); 620 thr[b] = Max (thr[b], gfc->ATH->adjust * gfc->ATH->cb[b]); 621 thr[b] = Min (thr[b], ecb); 622 } 623 else if (gfc->blocktype_old[chn>1 ? chn-2 : chn] == SHORT_TYPE ) 624 thr[b] = ecb; /* Min(ecb, rpelev*gfc->nb_1[chn][b]); */ 625 else 626 thr[b] = Min(ecb, Min(rpelev*gfc->nb_1[chn][b],rpelev2*gfc->nb_2[chn][b]) ); 627 628 gfc->nb_2[chn][b] = gfc->nb_1[chn][b]; 629 gfc->nb_1[chn][b] = ecb; 630 631 { 632 FLOAT8 thrpe; 633 thrpe = Max(thr[b],gfc->ATH->cb[b]); 634 /* 635 printf("%i thr=%e ATH=%e \n",b,thr[b],gfc->ATH->cb[b]); 636 */ 637 if (thrpe < eb[b]) 638 gfc->pe[chn] -= gfc->numlines_l[b] * log(thrpe / eb[b]); 639 } 640 } 641 642 /*************************************************************** 643 * determine the block type (window type) based on L & R channels 644 * 645 ***************************************************************/ 646 { /* compute PE for all 4 channels */ 647 FLOAT mn,mx,ma=0,mb=0,mc=0; 648 649 for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++) 650 { 651 ma += gfc->energy_s[0][j]; 652 mb += gfc->energy_s[1][j]; 653 mc += gfc->energy_s[2][j]; 654 } 655 mn = Min(ma,mb); 656 mn = Min(mn,mc); 657 mx = Max(ma,mb); 658 mx = Max(mx,mc); 659 660 /* bit allocation is based on pe. */ 661 if (mx>mn) { 662 FLOAT8 tmp = 400*log(mx/(1e-12+mn)); 663 if (tmp>gfc->pe[chn]) gfc->pe[chn]=tmp; 664 } 665 666 /* block type is based just on L or R channel */ 667 if (chn<2) { 668 uselongblock[chn] = 1; 669 670 /* tuned for t1.wav. doesnt effect most other samples */ 671 if (gfc->pe[chn] > 3000) 672 uselongblock[chn]=0; 673 674 if ( mx > 30*mn ) 675 {/* big surge of energy - always use short blocks */ 676 uselongblock[chn] = 0; 677 } 678 else if ((mx > 10*mn) && (gfc->pe[chn] > 1000)) 679 {/* medium surge, medium pe - use short blocks */ 680 uselongblock[chn] = 0; 681 } 682 683 /* disable short blocks */ 684 if (gfp->no_short_blocks) 685 uselongblock[chn]=1; 686 } 687 } 688 689 if (gfp->analysis) { 690 FLOAT mn,mx,ma=0,mb=0,mc=0; 691 692 for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++) 693 { 694 ma += gfc->energy_s[0][j]; 695 mb += gfc->energy_s[1][j]; 696 mc += gfc->energy_s[2][j]; 697 } 698 mn = Min(ma,mb); 699 mn = Min(mn,mc); 700 mx = Max(ma,mb); 701 mx = Max(mx,mc); 702 703 gfc->pinfo->ers[gr_out][chn]=gfc->ers_save[chn]; 704 gfc->ers_save[chn]=(mx/(1e-12+mn)); 705 gfc->pinfo->pe[gr_out][chn]=gfc->pe_save[chn]; 706 gfc->pe_save[chn]=gfc->pe[chn]; 707 } 708 709 710 /*************************************************************** 711 * compute masking thresholds for both short and long blocks 712 ***************************************************************/ 713 /* longblock threshold calculation (part 2) */ 714 for ( sb = 0; sb < NBPSY_l; sb++ ) 715 { 716 FLOAT8 enn = gfc->w1_l[sb] * eb[gfc->bu_l[sb]] + gfc->w2_l[sb] * eb[gfc->bo_l[sb]]; 717 FLOAT8 thmm = gfc->w1_l[sb] *thr[gfc->bu_l[sb]] + gfc->w2_l[sb] * thr[gfc->bo_l[sb]]; 718 719 for ( b = gfc->bu_l[sb]+1; b < gfc->bo_l[sb]; b++ ) 720 { 721 enn += eb[b]; 722 thmm += thr[b]; 723 } 724 725 gfc->en [chn].l[sb] = enn; 726 gfc->thm[chn].l[sb] = thmm; 727 } 728 729 730 /* threshold calculation for short blocks */ 731 for ( sblock = 0; sblock < 3; sblock++ ) 732 { 733 j = 0; 734 for ( b = 0; b < gfc->npart_s_orig; b++ ) 735 { 736 FLOAT ecb = gfc->energy_s[sblock][j++]; 737 for (i = 1 ; i<gfc->numlines_s[b]; ++i) 738 { 739 ecb += gfc->energy_s[sblock][j++]; 740 } 741 eb[b] = ecb; 742 } 743 744 for ( b = 0; b < gfc->npart_s; b++ ) 745 { 746 FLOAT8 ecb = 0; 747 for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) { 748 ecb += gfc->s3_s[b][k] * eb[k]; 749 } 750 ecb *= gfc->SNR_s[b]; 751 thr[b] = Max (1e-6, ecb); 752 } 753 754 for ( sb = 0; sb < NBPSY_s; sb++ ){ 755 FLOAT8 enn = gfc->w1_s[sb] * eb[gfc->bu_s[sb]] + gfc->w2_s[sb] * eb[gfc->bo_s[sb]]; 756 FLOAT8 thmm = gfc->w1_s[sb] *thr[gfc->bu_s[sb]] + gfc->w2_s[sb] * thr[gfc->bo_s[sb]]; 757 758 for ( b = gfc->bu_s[sb]+1; b < gfc->bo_s[sb]; b++ ) { 759 enn += eb[b]; 760 thmm += thr[b]; 761 } 762 763 gfc->en [chn].s[sb][sblock] = enn; 764 gfc->thm[chn].s[sb][sblock] = thmm; 765 } 766 } 767 } /* end loop over chn */ 768 769 770 /* compute M/S thresholds from Johnston & Ferreira 1992 ICASSP paper */ 771 if (gfp->mode == JOINT_STEREO) { 772 FLOAT8 rside,rmid,mld; 773 int chmid=2,chside=3; 774 775 for ( sb = 0; sb < NBPSY_l; sb++ ) { 776 /* use this fix if L & R masking differs by 2db or less */ 777 /* if db = 10*log10(x2/x1) < 2 */ 778 /* if (x2 < 1.58*x1) { */ 779 if (gfc->thm[0].l[sb] <= 1.58*gfc->thm[1].l[sb] 780 && gfc->thm[1].l[sb] <= 1.58*gfc->thm[0].l[sb]) { 781 782 mld = gfc->mld_l[sb]*gfc->en[chside].l[sb]; 783 rmid = Max(gfc->thm[chmid].l[sb], Min(gfc->thm[chside].l[sb],mld)); 784 785 mld = gfc->mld_l[sb]*gfc->en[chmid].l[sb]; 786 rside = Max(gfc->thm[chside].l[sb],Min(gfc->thm[chmid].l[sb],mld)); 787 788 gfc->thm[chmid].l[sb]=rmid; 789 gfc->thm[chside].l[sb]=rside; 790 } 791 } 792 for ( sb = 0; sb < NBPSY_s; sb++ ) { 793 for ( sblock = 0; sblock < 3; sblock++ ) { 794 if (gfc->thm[0].s[sb][sblock] <= 1.58*gfc->thm[1].s[sb][sblock] 795 && gfc->thm[1].s[sb][sblock] <= 1.58*gfc->thm[0].s[sb][sblock]) { 796 797 mld = gfc->mld_s[sb]*gfc->en[chside].s[sb][sblock]; 798 rmid = Max(gfc->thm[chmid].s[sb][sblock],Min(gfc->thm[chside].s[sb][sblock],mld)); 799 800 mld = gfc->mld_s[sb]*gfc->en[chmid].s[sb][sblock]; 801 rside = Max(gfc->thm[chside].s[sb][sblock],Min(gfc->thm[chmid].s[sb][sblock],mld)); 802 803 gfc->thm[chmid].s[sb][sblock]=rmid; 804 gfc->thm[chside].s[sb][sblock]=rside; 805 } 806 } 807 } 808 } 809 810 if (gfp->mode == JOINT_STEREO) { 811 /* determin ms_ratio from masking thresholds*/ 812 /* use ms_stereo (ms_ratio < .35) if average thresh. diff < 5 db */ 813 FLOAT8 db,x1,x2,sidetot=0,tot=0; 814 for (sb= NBPSY_l/4 ; sb< NBPSY_l; sb ++ ) { 815 x1 = Min(gfc->thm[0].l[sb],gfc->thm[1].l[sb]); 816 x2 = Max(gfc->thm[0].l[sb],gfc->thm[1].l[sb]); 817 /* thresholds difference in db */ 818 if (x2 >= 1000*x1) db=3; 819 else db = log10(x2/x1); 820 /* DEBUGF(gfc,"db = %f %e %e \n",db,gfc->thm[0].l[sb],gfc->thm[1].l[sb]);*/ 821 sidetot += db; 822 tot++; 823 } 824 ms_ratio_l= (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */ 825 ms_ratio_l = Min(ms_ratio_l,0.5); 826 827 sidetot=0; tot=0; 828 for ( sblock = 0; sblock < 3; sblock++ ) 829 for ( sb = NBPSY_s/4; sb < NBPSY_s; sb++ ) { 830 x1 = Min(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]); 831 x2 = Max(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]); 832 /* thresholds difference in db */ 833 if (x2 >= 1000*x1) db=3; 834 else db = log10(x2/x1); 835 sidetot += db; 836 tot++; 837 } 838 ms_ratio_s = (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */ 839 ms_ratio_s = Min(ms_ratio_s,.5); 840 } 841 842 /*************************************************************** 843 * determine final block type 844 ***************************************************************/ 845 846 for (chn=0; chn<gfc->channels_out; chn++) { 847 blocktype[chn] = NORM_TYPE; 848 } 849 850 851 if (gfc->channels_out==2) { 852 if (!gfp->allow_diff_short || gfp->mode==JOINT_STEREO) { 853 /* force both channels to use the same block type */ 854 /* this is necessary if the frame is to be encoded in ms_stereo. */ 855 /* But even without ms_stereo, FhG does this */ 856 int bothlong= (uselongblock[0] && uselongblock[1]); 857 if (!bothlong) { 858 uselongblock[0]=0; 859 uselongblock[1]=0; 860 } 861 } 862 } 863 864 865 866 /* update the blocktype of the previous granule, since it depends on what 867 * happend in this granule */ 868 for (chn=0; chn<gfc->channels_out; chn++) { 869 if ( uselongblock[chn]) 870 { /* no attack : use long blocks */ 871 assert( gfc->blocktype_old[chn] != START_TYPE ); 872 switch( gfc->blocktype_old[chn] ) 873 { 874 case NORM_TYPE: 875 case STOP_TYPE: 876 blocktype[chn] = NORM_TYPE; 877 break; 878 case SHORT_TYPE: 879 blocktype[chn] = STOP_TYPE; 880 break; 881 } 882 } else { 883 /* attack : use short blocks */ 884 blocktype[chn] = SHORT_TYPE; 885 if ( gfc->blocktype_old[chn] == NORM_TYPE ) { 886 gfc->blocktype_old[chn] = START_TYPE; 887 } 888 if ( gfc->blocktype_old[chn] == STOP_TYPE ) { 889 gfc->blocktype_old[chn] = SHORT_TYPE ; 890 } 891 } 892 893 blocktype_d[chn] = gfc->blocktype_old[chn]; /* value returned to calling program */ 894 gfc->blocktype_old[chn] = blocktype[chn]; /* save for next call to l3psy_anal */ 895 } 896 897 if (blocktype_d[0]==2) 898 *ms_ratio = gfc->ms_ratio_s_old; 899 else 900 *ms_ratio = gfc->ms_ratio_l_old; 901 902 gfc->ms_ratio_s_old = ms_ratio_s; 903 gfc->ms_ratio_l_old = ms_ratio_l; 904 905 /* we dont know the block type of this frame yet - assume long */ 906 *ms_ratio_next = ms_ratio_l; 907 908 return 0; 909 } 910 911 /* addition of simultaneous masking Naoki Shibata 2000/7 */ 912 inline static FLOAT8 mask_add(FLOAT8 m1,FLOAT8 m2,int k,int b, lame_internal_flags * const gfc) 913 { 914 static const FLOAT8 table1[] = { 915 3.3246 *3.3246 ,3.23837*3.23837,3.15437*3.15437,3.00412*3.00412,2.86103*2.86103,2.65407*2.65407,2.46209*2.46209,2.284 *2.284 , 916 2.11879*2.11879,1.96552*1.96552,1.82335*1.82335,1.69146*1.69146,1.56911*1.56911,1.46658*1.46658,1.37074*1.37074,1.31036*1.31036, 917 1.25264*1.25264,1.20648*1.20648,1.16203*1.16203,1.12765*1.12765,1.09428*1.09428,1.0659 *1.0659 ,1.03826*1.03826,1.01895*1.01895, 918 1 919 }; 920 921 static const FLOAT8 table2[] = { 922 1.33352*1.33352,1.35879*1.35879,1.38454*1.38454,1.39497*1.39497,1.40548*1.40548,1.3537 *1.3537 ,1.30382*1.30382,1.22321*1.22321, 923 1.14758*1.14758 924 }; 925 926 static const FLOAT8 table3[] = { 927 2.35364*2.35364,2.29259*2.29259,2.23313*2.23313,2.12675*2.12675,2.02545*2.02545,1.87894*1.87894,1.74303*1.74303,1.61695*1.61695, 928 1.49999*1.49999,1.39148*1.39148,1.29083*1.29083,1.19746*1.19746,1.11084*1.11084,1.03826*1.03826 929 }; 930 931 932 int i; 933 FLOAT8 m; 934 935 if (m1 == 0) return m2; 936 937 if (b < 0) b = -b; 938 939 i = 10*log10(m2 / m1)/10*16; 940 m = 10*log10((m1+m2)/gfc->ATH->cb[k]); 941 942 if (i < 0) i = -i; 943 944 if (b <= 3) { /* approximately, 1 bark = 3 partitions */ 945 if (i > 8) return m1+m2; 946 return (m1+m2)*table2[i]; 947 } 948 949 if (m<15) { 950 if (m > 0) { 951 FLOAT8 f=1.0,r; 952 if (i > 24) return m1+m2; 953 if (i > 13) f = 1; else f = table3[i]; 954 r = (m-0)/15; 955 return (m1+m2)*(table1[i]*r+f*(1-r)); 956 } 957 if (i > 13) return m1+m2; 958 return (m1+m2)*table3[i]; 959 } 960 961 if (i > 24) return m1+m2; 962 return (m1+m2)*table1[i]; 963 } 964 965 int L3psycho_anal_ns( lame_global_flags * gfp, 966 const sample_t *buffer[2], int gr_out, 967 FLOAT8 *ms_ratio, 968 FLOAT8 *ms_ratio_next, 969 III_psy_ratio masking_ratio[2][2], 970 III_psy_ratio masking_MS_ratio[2][2], 971 FLOAT8 percep_entropy[2],FLOAT8 percep_MS_entropy[2], 972 FLOAT8 energy[4], 973 int blocktype_d[2]) 974 { 975 /* to get a good cache performance, one has to think about 976 * the sequence, in which the variables are used. 977 * (Note: these static variables have been moved to the gfc-> struct, 978 * and their order in memory is layed out in util.h) 979 */ 980 lame_internal_flags *gfc=gfp->internal_flags; 981 982 /* fft and energy calculation */ 983 FLOAT (*wsamp_l)[BLKSIZE]; 984 FLOAT (*wsamp_s)[3][BLKSIZE_s]; 985 986 /* convolution */ 987 FLOAT8 eb[CBANDS],eb2[CBANDS]; 988 FLOAT8 thr[CBANDS]; 989 990 FLOAT8 max[CBANDS],avg[CBANDS],tonality2[CBANDS]; 991 992 /* ratios */ 993 FLOAT8 ms_ratio_l=0,ms_ratio_s=0; 994 995 /* block type */ 996 int blocktype[2],uselongblock[2]; 997 998 /* usual variables like loop indices, etc.. */ 999 int numchn, chn; 1000 int b, i, j, k; 1001 int sb,sblock; 1002 1003 /* variables used for --nspsytune */ 1004 int ns_attacks[4]; 1005 FLOAT ns_hpfsmpl[4][576+576/3+NSFIRLEN]; 1006 FLOAT pe_l[4],pe_s[4]; 1007 FLOAT pcfact; 1008 1009 1010 if(gfc->psymodel_init==0) { 1011 psymodel_init(gfp); 1012 init_fft(gfc); 1013 gfc->psymodel_init=1; 1014 } 1015 1016 1017 numchn = gfc->channels_out; 1018 /* chn=2 and 3 = Mid and Side channels */ 1019 if (gfp->mode == JOINT_STEREO) numchn=4; 1020 1021 if (gfp->VBR==vbr_off) pcfact = gfc->ResvMax == 0 ? 0 : ((FLOAT)gfc->ResvSize)/gfc->ResvMax*0.5; 1022 else if (gfp->VBR == vbr_rh || gfp->VBR == vbr_mtrh || gfp->VBR == vbr_mt) { 1023 static const FLOAT8 pcQns[10]={1.0,1.0,1.0,0.8,0.6,0.5,0.4,0.3,0.2,0.1}; 1024 pcfact = pcQns[gfp->VBR_q]; 1025 } else pcfact = 1; 1026 1027 /********************************************************************** 1028 * Apply HPF of fs/4 to the input signal. 1029 * This is used for attack detection / handling. 1030 **********************************************************************/ 1031 1032 { 1033 static const FLOAT fircoef[] = { 1034 -8.65163e-18,-0.00851586,-6.74764e-18, 0.0209036, 1035 -3.36639e-17,-0.0438162 ,-1.54175e-17, 0.0931738, 1036 -5.52212e-17,-0.313819 , 0.5 ,-0.313819, 1037 -5.52212e-17, 0.0931738 ,-1.54175e-17,-0.0438162, 1038 -3.36639e-17, 0.0209036 ,-6.74764e-18,-0.00851586, 1039 -8.65163e-18, 1040 }; 1041 1042 for(chn=0;chn<gfc->channels_out;chn++) 1043 { 1044 FLOAT firbuf[576+576/3+NSFIRLEN]; 1045 1046 /* apply high pass filter of fs/4 */ 1047 1048 for(i=-NSFIRLEN;i<576+576/3;i++) 1049 firbuf[i+NSFIRLEN] = buffer[chn][576-350+(i)]; 1050 1051 for(i=0;i<576+576/3-NSFIRLEN;i++) 1052 { 1053 FLOAT sum = 0; 1054 for(j=0;j<NSFIRLEN;j++) 1055 sum += fircoef[j] * firbuf[i+j]; 1056 ns_hpfsmpl[chn][i] = sum; 1057 } 1058 for(;i<576+576/3;i++) 1059 ns_hpfsmpl[chn][i] = 0; 1060 } 1061 1062 if (gfp->mode == JOINT_STEREO) { 1063 for(i=0;i<576+576/3;i++) 1064 { 1065 ns_hpfsmpl[2][i] = ns_hpfsmpl[0][i]+ns_hpfsmpl[1][i]; 1066 ns_hpfsmpl[3][i] = ns_hpfsmpl[0][i]-ns_hpfsmpl[1][i]; 1067 } 1068 } 1069 } 1070 1071 1072 1073 /* there is a one granule delay. Copy maskings computed last call 1074 * into masking_ratio to return to calling program. 1075 */ 1076 for (chn=0; chn<numchn; chn++) { 1077 for (i=0; i<numchn; ++i) { 1078 energy[chn]=gfc->tot_ener[chn]; 1079 } 1080 } 1081 1082 1083 for (chn=0; chn<numchn; chn++) { 1084 /* there is a one granule delay. Copy maskings computed last call 1085 * into masking_ratio to return to calling program. 1086 */ 1087 pe_l[chn] = gfc->nsPsy.pe_l[chn]; 1088 pe_s[chn] = gfc->nsPsy.pe_s[chn]; 1089 1090 if (chn < 2) { 1091 /* LR maskings */ 1092 //percep_entropy [chn] = gfc -> pe [chn]; 1093 masking_ratio [gr_out] [chn] .en = gfc -> en [chn]; 1094 masking_ratio [gr_out] [chn] .thm = gfc -> thm [chn]; 1095 } else { 1096 /* MS maskings */ 1097 //percep_MS_entropy [chn-2] = gfc -> pe [chn]; 1098 masking_MS_ratio [gr_out] [chn-2].en = gfc -> en [chn]; 1099 masking_MS_ratio [gr_out] [chn-2].thm = gfc -> thm [chn]; 1100 } 1101 } 1102 1103 for (chn=0; chn<numchn; chn++) { 1104 1105 /********************************************************************** 1106 * compute FFTs 1107 **********************************************************************/ 1108 1109 wsamp_s = gfc->wsamp_S+(chn & 1); 1110 wsamp_l = gfc->wsamp_L+(chn & 1); 1111 1112 if (chn<2) { 1113 fft_long ( gfc, *wsamp_l, chn, buffer); 1114 fft_short( gfc, *wsamp_s, chn, buffer); 1115 } 1116 1117 /* FFT data for mid and side channel is derived from L & R */ 1118 1119 if (chn == 2) 1120 { 1121 for (j = BLKSIZE-1; j >=0 ; --j) 1122 { 1123 FLOAT l = gfc->wsamp_L[0][j]; 1124 FLOAT r = gfc->wsamp_L[1][j]; 1125 gfc->wsamp_L[0][j] = (l+r)*(FLOAT)(SQRT2*0.5); 1126 gfc->wsamp_L[1][j] = (l-r)*(FLOAT)(SQRT2*0.5); 1127 } 1128 for (b = 2; b >= 0; --b) 1129 { 1130 for (j = BLKSIZE_s-1; j >= 0 ; --j) 1131 { 1132 FLOAT l = gfc->wsamp_S[0][b][j]; 1133 FLOAT r = gfc->wsamp_S[1][b][j]; 1134 gfc->wsamp_S[0][b][j] = (l+r)*(FLOAT)(SQRT2*0.5); 1135 gfc->wsamp_S[1][b][j] = (l-r)*(FLOAT)(SQRT2*0.5); 1136 } 1137 } 1138 } 1139 1140 1141 /********************************************************************** 1142 * compute energies for each spectral line 1143 **********************************************************************/ 1144 1145 /* long block */ 1146 1147 gfc->energy[0] = (*wsamp_l)[0]; 1148 gfc->energy[0] *= gfc->energy[0]; 1149 1150 gfc->tot_ener[chn] = gfc->energy[0]; /* sum total energy at nearly no extra cost */ 1151 1152 for (j=BLKSIZE/2-1; j >= 0; --j) 1153 { 1154 FLOAT re = (*wsamp_l)[BLKSIZE/2-j]; 1155 FLOAT im = (*wsamp_l)[BLKSIZE/2+j]; 1156 gfc->energy[BLKSIZE/2-j] = (re * re + im * im) * (FLOAT)0.5; 1157 1158 if (BLKSIZE/2-j > 10) 1159 gfc->tot_ener[chn] += gfc->energy[BLKSIZE/2-j]; 1160 } 1161 1162 1163 /* short block */ 1164 1165 for (b = 2; b >= 0; --b) 1166 { 1167 gfc->energy_s[b][0] = (*wsamp_s)[b][0]; 1168 gfc->energy_s[b][0] *= gfc->energy_s [b][0]; 1169 for (j=BLKSIZE_s/2-1; j >= 0; --j) 1170 { 1171 FLOAT re = (*wsamp_s)[b][BLKSIZE_s/2-j]; 1172 FLOAT im = (*wsamp_s)[b][BLKSIZE_s/2+j]; 1173 gfc->energy_s[b][BLKSIZE_s/2-j] = (re * re + im * im) * (FLOAT)0.5; 1174 } 1175 } 1176 1177 1178 /* output data for analysis */ 1179 1180 if (gfp->analysis) { 1181 for (j=0; j<HBLKSIZE ; j++) { 1182 gfc->pinfo->energy[gr_out][chn][j]=gfc->energy_save[chn][j]; 1183 gfc->energy_save[chn][j]=gfc->energy[j]; 1184 } 1185 } 1186 1187 1188 /********************************************************************** 1189 * Calculate the energy and the tonality of each partition. 1190 **********************************************************************/ 1191 1192 for (b=0, j=0; b<gfc->npart_l_orig; b++) 1193 { 1194 FLOAT8 ebb,m,a; 1195 1196 ebb = gfc->energy[j]; 1197 m = a = gfc->energy[j]; 1198 j++; 1199 1200 for (i = gfc->numlines_l[b] - 1; i > 0; i--) 1201 { 1202 FLOAT8 el = gfc->energy[j]; 1203 ebb += gfc->energy[j]; 1204 a += el; 1205 m = m < el ? el : m; 1206 j++; 1207 } 1208 eb[b] = ebb; 1209 max[b] = m; 1210 avg[b] = a / gfc->numlines_l[b]; 1211 } 1212 1213 j = 0; 1214 for (b=0; b < gfc->npart_l_orig; b++ ) 1215 { 1216 int c1,c2; 1217 FLOAT8 m,a; 1218 tonality2[b] = 0; 1219 c1 = c2 = 0; 1220 m = a = 0; 1221 for(k=b-1;k<=b+1;k++) 1222 { 1223 if (k >= 0 && k < gfc->npart_l_orig) { 1224 c1++; 1225 c2 += gfc->numlines_l[k]; 1226 a += avg[k]; 1227 m = m < max[k] ? max[k] : m; 1228 } 1229 } 1230 1231 a /= c1; 1232 tonality2[b] = a == 0 ? 0 : (m / a - 1)/(c2-1); 1233 } 1234 1235 for (b=0; b < gfc->npart_l_orig; b++ ) 1236 { 1237 #if 0 1238 static FLOAT8 tab[20] = 1239 { 0, 1, 2, 2, 2, 2, 2, 6,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3}; 1240 1241 static int init = 1; 1242 if (init) { 1243 int j; 1244 for(j=0;j<20;j++) { 1245 tab[j] = pow(10.0,-tab[j]/10.0); 1246 } 1247 init = 0; 1248 } 1249 #else 1250 static FLOAT8 tab[20] = { 1251 1,0.79433,0.63096,0.63096,0.63096,0.63096,0.63096,0.25119,0.11749,0.11749, 1252 0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749 1253 }; 1254 #endif 1255 1256 int t = 20*tonality2[b]; 1257 if (t > 19) t = 19; 1258 eb2[b] = eb[b] * tab[t]; 1259 } 1260 1261 1262 /********************************************************************** 1263 * convolve the partitioned energy and unpredictability 1264 * with the spreading function, s3_l[b][k] 1265 ******************************************************************** */ 1266 1267 for ( b = 0;b < gfc->npart_l; b++ ) 1268 { 1269 FLOAT8 ecb; 1270 1271 /**** convolve the partitioned energy with the spreading function ****/ 1272 1273 ecb = 0; 1274 1275 #if 1 1276 for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ ) 1277 { 1278 ecb = mask_add(ecb,gfc->s3_l[b][k] * eb2[k],k,k-b,gfc); 1279 } 1280 1281 ecb *= 0.158489319246111; // pow(10,-0.8) 1282 #endif 1283 1284 #if 0 1285 for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ ) 1286 { 1287 ecb += gfc->s3_l[k][b] * eb2[k]; 1288 } 1289 1290 ecb *= 0.223872113856834; // pow(10,-0.65); 1291 #endif 1292 1293 /**** long block pre-echo control ****/ 1294 1295 /* dont use long block pre-echo control if previous granule was 1296 * a short block. This is to avoid the situation: 1297 * frame0: quiet (very low masking) 1298 * frame1: surge (triggers short blocks) 1299 * frame2: regular frame. looks like pre-echo when compared to 1300 * frame0, but all pre-echo was in frame1. 1301 */ 1302 1303 /* chn=0,1 L and R channels 1304 chn=2,3 S and M channels. 1305 */ 1306 1307 #define NS_INTERP(x,y,r) (pow((x),(r))*pow((y),1-(r))) 1308 1309 if (gfc->blocktype_old[chn>1 ? chn-2 : chn] == SHORT_TYPE ) 1310 thr[b] = ecb; /* Min(ecb, rpelev*gfc->nb_1[chn][b]); */ 1311 else 1312 thr[b] = NS_INTERP(Min(ecb, Min(rpelev*gfc->nb_1[chn][b],rpelev2*gfc->nb_2[chn][b])),ecb,pcfact); 1313 1314 gfc->nb_2[chn][b] = gfc->nb_1[chn][b]; 1315 gfc->nb_1[chn][b] = ecb; 1316 } 1317 1318 1319 /*************************************************************** 1320 * determine the block type (window type) 1321 ***************************************************************/ 1322 1323 { 1324 static int count=0; 1325 FLOAT en_subshort[12]; 1326 FLOAT attack_intensity[12]; 1327 int ns_uselongblock = 1; 1328 1329 /* calculate energies of each sub-shortblocks */ 1330 1331 k = 0; 1332 for(i=0;i<12;i++) 1333 { 1334 en_subshort[i] = 0; 1335 for(j=0;j<576/9;j++) 1336 { 1337 en_subshort[i] += ns_hpfsmpl[chn][k] * ns_hpfsmpl[chn][k]; 1338 k++; 1339 } 1340 1341 if (en_subshort[i] < 100) en_subshort[i] = 100; 1342 } 1343 1344 /* compare energies between sub-shortblocks */ 1345 1346 #define NSATTACKTHRE 150 1347 #define NSATTACKTHRE_S 300 1348 1349 for(i=0;i<2;i++) { 1350 attack_intensity[i] = en_subshort[i] / gfc->nsPsy.last_en_subshort[chn][7+i]; 1351 } 1352 1353 for(;i<12;i++) { 1354 attack_intensity[i] = en_subshort[i] / en_subshort[i-2]; 1355 } 1356 1357 ns_attacks[0] = ns_attacks[1] = ns_attacks[2] = ns_attacks[3] = 0; 1358 1359 for(i=0;i<12;i++) 1360 { 1361 if (!ns_attacks[i/3] && attack_intensity[i] > (chn == 3 ? NSATTACKTHRE_S : NSATTACKTHRE)) ns_attacks[i/3] = (i % 3)+1; 1362 } 1363 1364 if (ns_attacks[0] && gfc->nsPsy.last_attacks[chn][2]) ns_attacks[0] = 0; 1365 if (ns_attacks[1] && ns_attacks[0]) ns_attacks[1] = 0; 1366 if (ns_attacks[2] && ns_attacks[1]) ns_attacks[2] = 0; 1367 if (ns_attacks[3] && ns_attacks[2]) ns_attacks[3] = 0; 1368 1369 if (gfc->nsPsy.last_attacks[chn][2] == 3 || 1370 ns_attacks[0] || ns_attacks[1] || ns_attacks[2] || ns_attacks[3]) ns_uselongblock = 0; 1371 1372 if (chn < 4) count++; 1373 1374 for(i=0;i<9;i++) 1375 { 1376 gfc->nsPsy.last_en_subshort[chn][i] = en_subshort[i]; 1377 gfc->nsPsy.last_attack_intensity[chn][i] = attack_intensity[i]; 1378 } 1379 1380 if (gfp->no_short_blocks) { 1381 uselongblock[chn] = 1; 1382 } else { 1383 if (chn < 2) { 1384 uselongblock[chn] = ns_uselongblock; 1385 } else { 1386 if (ns_uselongblock == 0) uselongblock[0] = uselongblock[1] = 0; 1387 } 1388 } 1389 } 1390 1391 if (gfp->analysis) { 1392 FLOAT mn,mx,ma=0,mb=0,mc=0; 1393 1394 for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++) 1395 { 1396 ma += gfc->energy_s[0][j]; 1397 mb += gfc->energy_s[1][j]; 1398 mc += gfc->energy_s[2][j]; 1399 } 1400 mn = Min(ma,mb); 1401 mn = Min(mn,mc); 1402 mx = Max(ma,mb); 1403 mx = Max(mx,mc); 1404 1405 gfc->pinfo->ers[gr_out][chn]=gfc->ers_save[chn]; 1406 gfc->ers_save[chn]=(mx/(1e-12+mn)); 1407 } 1408 1409 1410 /*************************************************************** 1411 * compute masking thresholds for long blocks 1412 ***************************************************************/ 1413 1414 for ( sb = 0; sb < NBPSY_l; sb++ ) 1415 { 1416 FLOAT8 enn = gfc->w1_l[sb] * eb[gfc->bu_l[sb]] + gfc->w2_l[sb] * eb[gfc->bo_l[sb]]; 1417 FLOAT8 thmm = gfc->w1_l[sb] *thr[gfc->bu_l[sb]] + gfc->w2_l[sb] * thr[gfc->bo_l[sb]]; 1418 1419 for ( b = gfc->bu_l[sb]+1; b < gfc->bo_l[sb]; b++ ) 1420 { 1421 enn += eb[b]; 1422 thmm += thr[b]; 1423 } 1424 1425 gfc->en [chn].l[sb] = enn; 1426 gfc->thm[chn].l[sb] = thmm; 1427 } 1428 1429 1430 /*************************************************************** 1431 * compute masking thresholds for short blocks 1432 ***************************************************************/ 1433 1434 for ( sblock = 0; sblock < 3; sblock++ ) 1435 { 1436 j = 0; 1437 for ( b = 0; b < gfc->npart_s_orig; b++ ) 1438 { 1439 FLOAT ecb = gfc->energy_s[sblock][j++]; 1440 for (i = 1 ; i<gfc->numlines_s[b]; ++i) 1441 { 1442 ecb += gfc->energy_s[sblock][j++]; 1443 } 1444 eb[b] = ecb; 1445 } 1446 1447 for ( b = 0; b < gfc->npart_s; b++ ) 1448 { 1449 FLOAT8 ecb = 0; 1450 for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) 1451 { 1452 ecb += gfc->s3_s[b][k] * eb[k]; 1453 } 1454 thr[b] = Max (1e-6, ecb); 1455 } 1456 1457 for ( sb = 0; sb < NBPSY_s; sb++ ) 1458 { 1459 FLOAT8 enn = gfc->w1_s[sb] * eb[gfc->bu_s[sb]] + gfc->w2_s[sb] * eb[gfc->bo_s[sb]]; 1460 FLOAT8 thmm = gfc->w1_s[sb] *thr[gfc->bu_s[sb]] + gfc->w2_s[sb] * thr[gfc->bo_s[sb]]; 1461 1462 for ( b = gfc->bu_s[sb]+1; b < gfc->bo_s[sb]; b++ ) 1463 { 1464 enn += eb[b]; 1465 thmm += thr[b]; 1466 } 1467 1468 /**** short block pre-echo control ****/ 1469 1470 #define NS_PREECHO_ATT0 0.8 1471 #define NS_PREECHO_ATT1 0.6 1472 #define NS_PREECHO_ATT2 0.3 1473 1474 thmm *= NS_PREECHO_ATT0; 1475 1476 if (ns_attacks[sblock] >= 2) { 1477 if (sblock != 0) { 1478 double p = NS_INTERP(gfc->thm[chn].s[sb][sblock-1],thmm,NS_PREECHO_ATT1*pcfact); 1479 thmm = Min(thmm,p); 1480 } else { 1481 double p = NS_INTERP(gfc->nsPsy.last_thm[chn][sb][2],thmm,NS_PREECHO_ATT1*pcfact); 1482 thmm = Min(thmm,p); 1483 } 1484 } else if (ns_attacks[sblock+1] == 1) { 1485 if (sblock != 0) { 1486 double p = NS_INTERP(gfc->thm[chn].s[sb][sblock-1],thmm,NS_PREECHO_ATT1*pcfact); 1487 thmm = Min(thmm,p); 1488 } else { 1489 double p = NS_INTERP(gfc->nsPsy.last_thm[chn][sb][2],thmm,NS_PREECHO_ATT1*pcfact); 1490 thmm = Min(thmm,p); 1491 } 1492 } 1493 1494 if (ns_attacks[sblock] == 1) { 1495 double p = sblock == 0 ? gfc->nsPsy.last_thm[chn][sb][2] : gfc->thm[chn].s[sb][sblock-1]; 1496 p = NS_INTERP(p,thmm,NS_PREECHO_ATT2*pcfact); 1497 thmm = Min(thmm,p); 1498 } else if ((sblock != 0 && ns_attacks[sblock-1] == 3) || 1499 (sblock == 0 && gfc->nsPsy.last_attacks[chn][2] == 3)) { 1500 double p = sblock <= 1 ? gfc->nsPsy.last_thm[chn][sb][sblock+1] : gfc->thm[chn].s[sb][0]; 1501 p = NS_INTERP(p,thmm,NS_PREECHO_ATT2*pcfact); 1502 thmm = Min(thmm,p); 1503 } 1504 1505 gfc->en [chn].s[sb][sblock] = enn; 1506 gfc->thm[chn].s[sb][sblock] = thmm; 1507 } 1508 } 1509 1510 1511 /*************************************************************** 1512 * save some values for analysis of the next granule 1513 ***************************************************************/ 1514 1515 for ( sblock = 0; sblock < 3; sblock++ ) 1516 { 1517 for ( sb = 0; sb < NBPSY_s; sb++ ) 1518 { 1519 gfc->nsPsy.last_thm[chn][sb][sblock] = gfc->thm[chn].s[sb][sblock]; 1520 } 1521 } 1522 1523 for(i=0;i<3;i++) 1524 gfc->nsPsy.last_attacks[chn][i] = ns_attacks[i]; 1525 1526 } /* end loop over chn */ 1527 1528 1529 1530 /*************************************************************** 1531 * compute M/S thresholds 1532 ***************************************************************/ 1533 1534 /* from Johnston & Ferreira 1992 ICASSP paper */ 1535 1536 if ( numchn==4 /* mid/side and r/l */) { 1537 FLOAT8 rside,rmid,mld; 1538 int chmid=2,chside=3; 1539 1540 for ( sb = 0; sb < NBPSY_l; sb++ ) { 1541 /* use this fix if L & R masking differs by 2db or less */ 1542 /* if db = 10*log10(x2/x1) < 2 */ 1543 /* if (x2 < 1.58*x1) { */ 1544 if (gfc->thm[0].l[sb] <= 1.58*gfc->thm[1].l[sb] 1545 && gfc->thm[1].l[sb] <= 1.58*gfc->thm[0].l[sb]) { 1546 1547 mld = gfc->mld_l[sb]*gfc->en[chside].l[sb]; 1548 rmid = Max(gfc->thm[chmid].l[sb], Min(gfc->thm[chside].l[sb],mld)); 1549 1550 mld = gfc->mld_l[sb]*gfc->en[chmid].l[sb]; 1551 rside = Max(gfc->thm[chside].l[sb],Min(gfc->thm[chmid].l[sb],mld)); 1552 1553 gfc->thm[chmid].l[sb]=rmid; 1554 gfc->thm[chside].l[sb]=rside; 1555 } 1556 } 1557 for ( sb = 0; sb < NBPSY_s; sb++ ) { 1558 for ( sblock = 0; sblock < 3; sblock++ ) { 1559 if (gfc->thm[0].s[sb][sblock] <= 1.58*gfc->thm[1].s[sb][sblock] 1560 && gfc->thm[1].s[sb][sblock] <= 1.58*gfc->thm[0].s[sb][sblock]) { 1561 1562 mld = gfc->mld_s[sb]*gfc->en[chside].s[sb][sblock]; 1563 rmid = Max(gfc->thm[chmid].s[sb][sblock],Min(gfc->thm[chside].s[sb][sblock],mld)); 1564 1565 mld = gfc->mld_s[sb]*gfc->en[chmid].s[sb][sblock]; 1566 rside = Max(gfc->thm[chside].s[sb][sblock],Min(gfc->thm[chmid].s[sb][sblock],mld)); 1567 1568 gfc->thm[chmid].s[sb][sblock]=rmid; 1569 gfc->thm[chside].s[sb][sblock]=rside; 1570 } 1571 } 1572 } 1573 } 1574 1575 1576 /* Naoki Shibata 2000 */ 1577 1578 #define NS_MSFIX 3.5 1579 1580 if (numchn == 4) { 1581 FLOAT msfix = NS_MSFIX; 1582 if (gfc->nsPsy.safejoint) msfix = 1; 1583 1584 for ( sb = 0; sb < NBPSY_l; sb++ ) 1585 { 1586 FLOAT8 thmL,thmR,thmM,thmS,ath; 1587 ath = (gfc->ATH->cb[(gfc->bu_l[sb] + gfc->bo_l[sb])/2])*pow(10,-gfp->ATHlower/10.0); 1588 thmL = Max(gfc->thm[0].l[sb],ath); 1589 thmR = Max(gfc->thm[1].l[sb],ath); 1590 thmM = Max(gfc->thm[2].l[sb],ath); 1591 thmS = Max(gfc->thm[3].l[sb],ath); 1592 1593 if (thmL*msfix < (thmM+thmS)/2) { 1594 FLOAT8 f = thmL*msfix / ((thmM+thmS)/2); 1595 thmM *= f; 1596 thmS *= f; 1597 } 1598 if (thmR*msfix < (thmM+thmS)/2) { 1599 FLOAT8 f = thmR*msfix / ((thmM+thmS)/2); 1600 thmM *= f; 1601 thmS *= f; 1602 } 1603 1604 gfc->thm[2].l[sb] = Min(thmM,gfc->thm[2].l[sb]); 1605 gfc->thm[3].l[sb] = Min(thmS,gfc->thm[3].l[sb]); 1606 } 1607 1608 for ( sb = 0; sb < NBPSY_s; sb++ ) { 1609 for ( sblock = 0; sblock < 3; sblock++ ) { 1610 FLOAT8 thmL,thmR,thmM,thmS,ath; 1611 ath = (gfc->ATH->cb[(gfc->bu_s[sb] + gfc->bo_s[sb])/2])*pow(10,-gfp->ATHlower/10.0); 1612 thmL = Max(gfc->thm[0].s[sb][sblock],ath); 1613 thmR = Max(gfc->thm[1].s[sb][sblock],ath); 1614 thmM = Max(gfc->thm[2].s[sb][sblock],ath); 1615 thmS = Max(gfc->thm[3].s[sb][sblock],ath); 1616 1617 if (thmL*msfix < (thmM+thmS)/2) { 1618 FLOAT8 f = thmL*msfix / ((thmM+thmS)/2); 1619 thmM *= f; 1620 thmS *= f; 1621 } 1622 if (thmR*msfix < (thmM+thmS)/2) { 1623 FLOAT8 f = thmR*msfix / ((thmM+thmS)/2); 1624 thmM *= f; 1625 thmS *= f; 1626 } 1627 1628 gfc->thm[2].s[sb][sblock] = Min(gfc->thm[2].s[sb][sblock],thmM); 1629 gfc->thm[3].s[sb][sblock] = Min(gfc->thm[3].s[sb][sblock],thmS); 1630 } 1631 } 1632 } 1633 1634 ms_ratio_l = 0; 1635 ms_ratio_s = 0; 1636 1637 1638 /*************************************************************** 1639 * compute estimation of the amount of bit used in the granule 1640 ***************************************************************/ 1641 1642 for(chn=0;chn<numchn;chn++) 1643 { 1644 { 1645 static FLOAT8 regcoef[] = { 1646 1124.23,10.0583,10.7484,7.29006,16.2714,6.2345,4.09743,3.05468,3.33196,2.54688, 1647 3.68168,5.83109,2.93817,-8.03277,-10.8458,8.48777,9.13182,2.05212,8.6674,50.3937,73.267,97.5664,0 1648 }; 1649 1650 FLOAT8 msum = regcoef[0]/4; 1651 int sb; 1652 1653 for ( sb = 0; sb < NBPSY_l; sb++ ) 1654 { 1655 FLOAT8 t; 1656 1657 if (gfc->thm[chn].l[sb]*gfc->masking_lower != 0 && 1658 gfc->en[chn].l[sb]/(gfc->thm[chn].l[sb]*gfc->masking_lower) > 1) 1659 t = log(gfc->en[chn].l[sb]/(gfc->thm[chn].l[sb]*gfc->masking_lower)); 1660 else 1661 t = 0; 1662 msum += regcoef[sb+1] * t; 1663 } 1664 1665 gfc->nsPsy.pe_l[chn] = msum; 1666 } 1667 1668 { 1669 static FLOAT8 regcoef[] = { 1670 1236.28,0,0,0,0.434542,25.0738,0,0,0,19.5442,19.7486,60,100,0 1671 }; 1672 1673 FLOAT8 msum = regcoef[0]/4; 1674 int sb,sblock; 1675 1676 for(sblock=0;sblock<3;sblock++) 1677 { 1678 for ( sb = 0; sb < NBPSY_s; sb++ ) 1679 { 1680 FLOAT8 t; 1681 1682 if (gfc->thm[chn].s[sb][sblock] * gfc->masking_lower != 0 && 1683 gfc->en[chn].s[sb][sblock] / (gfc->thm[chn].s[sb][sblock] * gfc->masking_lower) > 1) 1684 t = log(gfc->en[chn].s[sb][sblock] / (gfc->thm[chn].s[sb][sblock] * gfc->masking_lower)); 1685 else 1686 t = 0; 1687 msum += regcoef[sb+1] * t; 1688 } 1689 } 1690 1691 gfc->nsPsy.pe_s[chn] = msum; 1692 } 1693 1694 //gfc->pe[chn] -= 150; 1695 } 1696 1697 /*************************************************************** 1698 * determine final block type 1699 ***************************************************************/ 1700 1701 for (chn=0; chn<gfc->channels_out; chn++) { 1702 blocktype[chn] = NORM_TYPE; 1703 } 1704 1705 if (gfc->channels_out==2) { 1706 if (!gfp->allow_diff_short || gfp->mode==JOINT_STEREO) { 1707 /* force both channels to use the same block type */ 1708 /* this is necessary if the frame is to be encoded in ms_stereo. */ 1709 /* But even without ms_stereo, FhG does this */ 1710 int bothlong= (uselongblock[0] && uselongblock[1]); 1711 if (!bothlong) { 1712 uselongblock[0]=0; 1713 uselongblock[1]=0; 1714 } 1715 } 1716 } 1717 1718 /* update the blocktype of the previous granule, since it depends on what 1719 * happend in this granule */ 1720 for (chn=0; chn<gfc->channels_out; chn++) { 1721 if ( uselongblock[chn]) 1722 { /* no attack : use long blocks */ 1723 assert( gfc->blocktype_old[chn] != START_TYPE ); 1724 switch( gfc->blocktype_old[chn] ) 1725 { 1726 case NORM_TYPE: 1727 case STOP_TYPE: 1728 blocktype[chn] = NORM_TYPE; 1729 break; 1730 case SHORT_TYPE: 1731 blocktype[chn] = STOP_TYPE; 1732 break; 1733 } 1734 } else { 1735 /* attack : use short blocks */ 1736 blocktype[chn] = SHORT_TYPE; 1737 if ( gfc->blocktype_old[chn] == NORM_TYPE ) { 1738 gfc->blocktype_old[chn] = START_TYPE; 1739 } 1740 if ( gfc->blocktype_old[chn] == STOP_TYPE ) { 1741 gfc->blocktype_old[chn] = SHORT_TYPE ; 1742 } 1743 } 1744 1745 blocktype_d[chn] = gfc->blocktype_old[chn]; /* value returned to calling program */ 1746 gfc->blocktype_old[chn] = blocktype[chn]; /* save for next call to l3psy_anal */ 1747 } 1748 1749 /********************************************************************* 1750 * compute the value of PE to return (one granule delay) 1751 *********************************************************************/ 1752 for(chn=0;chn<numchn;chn++) 1753 { 1754 if (chn < 2) { 1755 if (blocktype_d[chn] == SHORT_TYPE) { 1756 percep_entropy[chn] = pe_s[chn]; 1757 } else { 1758 percep_entropy[chn] = pe_l[chn]; 1759 } 1760 if (gfp->analysis) gfc->pinfo->pe[gr_out][chn] = percep_entropy[chn]; 1761 } else { 1762 if (blocktype_d[0] == SHORT_TYPE) { 1763 percep_MS_entropy[chn-2] = pe_s[chn]; 1764 } else { 1765 percep_MS_entropy[chn-2] = pe_l[chn]; 1766 } 1767 if (gfp->analysis) gfc->pinfo->pe[gr_out][chn] = percep_MS_entropy[chn-2]; 1768 } 1769 } 1770 1771 1772 return 0; 1773 } 1774 1775 1776 1777 1778 1779 /* 1780 * The spreading function. Values returned in units of energy 1781 */ 1782 FLOAT8 s3_func(FLOAT8 bark) { 1783 1784 FLOAT8 tempx,x,tempy,temp; 1785 tempx = bark; 1786 if (tempx>=0) tempx *= 3; 1787 else tempx *=1.5; 1788 1789 if(tempx>=0.5 && tempx<=2.5) 1790 { 1791 temp = tempx - 0.5; 1792 x = 8.0 * (temp*temp - 2.0 * temp); 1793 } 1794 else x = 0.0; 1795 tempx += 0.474; 1796 tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx); 1797 1798 if (tempy <= -60.0) return 0.0; 1799 1800 tempx = exp( (x + tempy)*LN_TO_LOG10 ); 1801 1802 /* Normalization. The spreading function should be normalized so that: 1803 +inf 1804 / 1805 | s3 [ bark ] d(bark) = 1 1806 / 1807 -inf 1808 */ 1809 tempx /= .6609193; 1810 return tempx; 1811 1812 } 1813 1814 1815 1816 1817 1818 1819 1820 1821 int L3para_read(lame_global_flags * gfp, FLOAT8 sfreq, int *numlines_l,int *numlines_s, 1822 FLOAT8 *minval, 1823 FLOAT8 s3_l[CBANDS][CBANDS], FLOAT8 s3_s[CBANDS][CBANDS], 1824 FLOAT8 *SNR, 1825 int *bu_l, int *bo_l, FLOAT8 *w1_l, FLOAT8 *w2_l, 1826 int *bu_s, int *bo_s, FLOAT8 *w1_s, FLOAT8 *w2_s, 1827 int *npart_l_orig,int *npart_l,int *npart_s_orig,int *npart_s) 1828 { 1829 lame_internal_flags *gfc=gfp->internal_flags; 1830 1831 1832 FLOAT8 freq_tp; 1833 FLOAT8 bval_l[CBANDS], bval_s[CBANDS]; 1834 FLOAT8 bval_l_width[CBANDS], bval_s_width[CBANDS]; 1835 int cbmax=0, cbmax_tp; 1836 int sbmax ; 1837 int i,j; 1838 int freq_scale=1; 1839 int partition[HBLKSIZE]; 1840 int loop, k2; 1841 1842 1843 1844 /* compute numlines, the number of spectral lines in each partition band */ 1845 /* each partition band should be about DELBARK wide. */ 1846 j=0; 1847 for(i=0;i<CBANDS;i++) 1848 { 1849 FLOAT8 ji, bark1,bark2; 1850 int k,j2; 1851 1852 j2 = j; 1853 j2 = Min(j2,BLKSIZE/2); 1854 1855 do { 1856 /* find smallest j2 >= j so that (bark - bark_l[i-1]) < DELBARK */ 1857 ji = j; 1858 bark1 = freq2bark(sfreq*ji/BLKSIZE); 1859 1860 ++j2; 1861 ji = j2; 1862 bark2 = freq2bark(sfreq*ji/BLKSIZE); 1863 } while ((bark2 - bark1) < DELBARK && j2<=BLKSIZE/2); 1864 1865 for (k=j; k<j2; ++k) 1866 partition[k]=i; 1867 numlines_l[i]=(j2-j); 1868 j = j2; 1869 if (j > BLKSIZE/2) break; 1870 } 1871 *npart_l_orig = i+1; 1872 assert(*npart_l_orig <= CBANDS); 1873 1874 /* compute which partition bands are in which scalefactor bands */ 1875 { int i1,i2,sfb,start,end; 1876 FLOAT8 freq1,freq2; 1877 for ( sfb = 0; sfb < SBMAX_l; sfb++ ) { 1878 start = gfc->scalefac_band.l[ sfb ]; 1879 end = gfc->scalefac_band.l[ sfb+1 ]; 1880 freq1 = sfreq*(start-.5)/(2*576); 1881 freq2 = sfreq*(end-1+.5)/(2*576); 1882 1883 i1 = floor(.5 + BLKSIZE*freq1/sfreq); 1884 if (i1<0) i1=0; 1885 i2 = floor(.5 + BLKSIZE*freq2/sfreq); 1886 if (i2>BLKSIZE/2) i2=BLKSIZE/2; 1887 1888 // DEBUGF(gfc,"longblock: old: (%i,%i) new: (%i,%i) %i %i \n",bu_l[sfb],bo_l[sfb],partition[i1],partition[i2],i1,i2); 1889 1890 w1_l[sfb]=.5; 1891 w2_l[sfb]=.5; 1892 bu_l[sfb]=partition[i1]; 1893 bo_l[sfb]=partition[i2]; 1894 1895 } 1896 } 1897 1898 1899 /* compute bark value and ATH of each critical band */ 1900 j = 0; 1901 for ( i = 0; i < *npart_l_orig; i++ ) { 1902 int k; 1903 FLOAT8 bark1,bark2; 1904 /* FLOAT8 mval,freq; */ 1905 1906 // Calculating the medium bark scaled frequency of the spectral lines 1907 // from j ... j + numlines[i]-1 (=spectral lines in parition band i) 1908 1909 k = numlines_l[i] - 1; 1910 bark1 = freq2bark(sfreq*(j+0)/BLKSIZE); 1911 bark2 = freq2bark(sfreq*(j+k)/BLKSIZE); 1912 bval_l[i] = .5*(bark1+bark2); 1913 1914 bark1 = freq2bark(sfreq*(j+0-.5)/BLKSIZE); 1915 bark2 = freq2bark(sfreq*(j+k+.5)/BLKSIZE); 1916 bval_l_width[i] = bark2-bark1; 1917 1918 gfc->ATH->cb [i] = 1.e37; // preinit for minimum search 1919 for (k=0; k < numlines_l[i]; k++, j++) { 1920 FLOAT8 freq = sfreq*j/(1000.0*BLKSIZE); 1921 FLOAT8 level; 1922 assert( freq <= 24 ); // or only '<' 1923 // freq = Min(.1,freq); // ATH below 100 Hz constant, not further climbing 1924 level = ATHformula (freq*1000, gfp) - 20; // scale to FFT units; returned value is in dB 1925 level = pow ( 10., 0.1*level ); // convert from dB -> energy 1926 level *= numlines_l [i]; 1927 if ( level < gfc->ATH->cb [i] ) 1928 gfc->ATH->cb [i] = level; 1929 } 1930 1931 1932 } 1933 1934 /* MINVAL. For low freq, the strength of the masking is limited by minval 1935 * this is an ISO MPEG1 thing, dont know if it is really needed */ 1936 for(i=0;i<*npart_l_orig;i++){ 1937 double x = (-20+bval_l[i]*20.0/10.0); 1938 if (bval_l[i]>10) x = 0; 1939 minval[i]=pow(10.0,x/10); 1940 } 1941 1942 1943 1944 1945 1946 1947 1948 /************************************************************************/ 1949 /* SHORT BLOCKS */ 1950 /************************************************************************/ 1951 1952 /* compute numlines */ 1953 j=0; 1954 for(i=0;i<CBANDS;i++) 1955 { 1956 FLOAT8 ji, bark1,bark2; 1957 int k,j2; 1958 1959 j2 = j; 1960 j2 = Min(j2,BLKSIZE_s/2); 1961 1962 do { 1963 /* find smallest j2 >= j so that (bark - bark_s[i-1]) < DELBARK */ 1964 ji = j; 1965 bark1 = freq2bark(sfreq*ji/BLKSIZE_s); 1966 1967 ++j2; 1968 ji = j2; 1969 bark2 = freq2bark(sfreq*ji/BLKSIZE_s); 1970 1971 } while ((bark2 - bark1) < DELBARK && j2<=BLKSIZE_s/2); 1972 1973 for (k=j; k<j2; ++k) 1974 partition[k]=i; 1975 numlines_s[i]=(j2-j); 1976 j = j2; 1977 if (j > BLKSIZE_s/2) break; 1978 } 1979 *npart_s_orig = i+1; 1980 assert(*npart_s_orig <= CBANDS); 1981 1982 /* compute which partition bands are in which scalefactor bands */ 1983 { int i1,i2,sfb,start,end; 1984 FLOAT8 freq1,freq2; 1985 for ( sfb = 0; sfb < SBMAX_s; sfb++ ) { 1986 start = gfc->scalefac_band.s[ sfb ]; 1987 end = gfc->scalefac_band.s[ sfb+1 ]; 1988 freq1 = sfreq*(start-.5)/(2*192); 1989 freq2 = sfreq*(end-1+.5)/(2*192); 1990 1991 i1 = floor(.5 + BLKSIZE_s*freq1/sfreq); 1992 if (i1<0) i1=0; 1993 i2 = floor(.5 + BLKSIZE_s*freq2/sfreq); 1994 if (i2>BLKSIZE_s/2) i2=BLKSIZE_s/2; 1995 1996 //DEBUGF(gfc,"shortblock: old: (%i,%i) new: (%i,%i) %i %i \n",bu_s[sfb],bo_s[sfb], partition[i1],partition[i2],i1,i2); 1997 1998 w1_s[sfb]=.5; 1999 w2_s[sfb]=.5; 2000 bu_s[sfb]=partition[i1]; 2001 bo_s[sfb]=partition[i2]; 2002 2003 } 2004 } 2005 2006 2007 2008 2009 2010 /* compute bark values of each critical band */ 2011 j = 0; 2012 for(i=0;i<*npart_s_orig;i++) 2013 { 2014 int k; 2015 FLOAT8 bark1,bark2,snr; 2016 k = numlines_s[i] - 1; 2017 2018 bark1 = freq2bark (sfreq*(j+0)/BLKSIZE_s); 2019 bark2 = freq2bark (sfreq*(j+k)/BLKSIZE_s); 2020 bval_s[i] = .5*(bark1+bark2); 2021 2022 bark1 = freq2bark (sfreq*(j+0-.5)/BLKSIZE_s); 2023 bark2 = freq2bark (sfreq*(j+k+.5)/BLKSIZE_s); 2024 bval_s_width[i] = bark2-bark1; 2025 j += k+1; 2026 2027 /* SNR formula */ 2028 if (bval_s[i]<13) 2029 snr=-8.25; 2030 else 2031 snr = -4.5 * (bval_s[i]-13)/(24.0-13.0) + 2032 -8.25*(bval_s[i]-24)/(13.0-24.0); 2033 2034 SNR[i]=pow(10.0,snr/10.0); 2035 } 2036 2037 2038 2039 2040 2041 2042 /************************************************************************ 2043 * Now compute the spreading function, s[j][i], the value of the spread-* 2044 * ing function, centered at band j, for band i, store for later use * 2045 ************************************************************************/ 2046 /* i.e.: sum over j to spread into signal barkval=i 2047 NOTE: i and j are used opposite as in the ISO docs */ 2048 for(i=0;i<*npart_l_orig;i++) { 2049 for(j=0;j<*npart_l_orig;j++) { 2050 s3_l[i][j]=s3_func(bval_l[i]-bval_l[j])*bval_l_width[j]; 2051 } 2052 } 2053 for(i=0;i<*npart_s_orig;i++) { 2054 for(j=0;j<*npart_s_orig;j++) { 2055 s3_s[i][j]=s3_func(bval_s[i]-bval_s[j])*bval_s_width[j]; 2056 } 2057 } 2058 2059 2060 2061 2062 /* compute: */ 2063 /* npart_l_orig = number of partition bands before convolution */ 2064 /* npart_l = number of partition bands after convolution */ 2065 2066 *npart_l=bo_l[NBPSY_l-1]+1; 2067 *npart_s=bo_s[NBPSY_s-1]+1; 2068 2069 assert(*npart_l <= *npart_l_orig); 2070 assert(*npart_s <= *npart_s_orig); 2071 2072 2073 /* setup stereo demasking thresholds */ 2074 /* formula reverse enginerred from plot in paper */ 2075 for ( i = 0; i < NBPSY_s; i++ ) { 2076 FLOAT8 arg,mld; 2077 arg = freq2bark(sfreq*gfc->scalefac_band.s[i]/(2*192)); 2078 arg = (Min(arg, 15.5)/15.5); 2079 2080 mld = 1.25*(1-cos(PI*arg))-2.5; 2081 gfc->mld_s[i] = pow(10.0,mld); 2082 } 2083 for ( i = 0; i < NBPSY_l; i++ ) { 2084 FLOAT8 arg,mld; 2085 arg = freq2bark(sfreq*gfc->scalefac_band.l[i]/(2*576)); 2086 arg = (Min(arg, 15.5)/15.5); 2087 2088 mld = 1.25*(1-cos(PI*arg))-2.5; 2089 gfc->mld_l[i] = pow(10.0,mld); 2090 } 2091 2092 #define temporalmask_sustain_sec 0.01 2093 2094 /* setup temporal masking */ 2095 gfc->decay = exp(-1.0*LOG10/(temporalmask_sustain_sec*sfreq/192.0)); 2096 2097 return 0; 2098 } 2099 2100 2101 2102 2103 2104 2105 2106 2107 2108 2109 2110 2111 2112 2113 2114 2115 2116 2117 2118 int psymodel_init(lame_global_flags *gfp) 2119 { 2120 lame_internal_flags *gfc=gfp->internal_flags; 2121 int i,j,b,sb,k,samplerate; 2122 FLOAT cwlimit; 2123 2124 samplerate = gfp->out_samplerate; 2125 gfc->ms_ener_ratio_old=.25; 2126 gfc->blocktype_old[0]=STOP_TYPE; 2127 gfc->blocktype_old[1]=STOP_TYPE; 2128 gfc->blocktype_old[0]=SHORT_TYPE; 2129 gfc->blocktype_old[1]=SHORT_TYPE; 2130 2131 for (i=0; i<4; ++i) { 2132 for (j=0; j<CBANDS; ++j) { 2133 gfc->nb_1[i][j]=1e20; 2134 gfc->nb_2[i][j]=1e20; 2135 } 2136 for ( sb = 0; sb < NBPSY_l; sb++ ) { 2137 gfc->en[i].l[sb] = 1e20; 2138 gfc->thm[i].l[sb] = 1e20; 2139 } 2140 for (j=0; j<3; ++j) { 2141 for ( sb = 0; sb < NBPSY_s; sb++ ) { 2142 gfc->en[i].s[sb][j] = 1e20; 2143 gfc->thm[i].s[sb][j] = 1e20; 2144 } 2145 } 2146 } 2147 for (i=0; i<4; ++i) { 2148 for (j=0; j<3; ++j) { 2149 for ( sb = 0; sb < NBPSY_s; sb++ ) { 2150 gfc->nsPsy.last_thm[i][sb][j] = 1e20; 2151 } 2152 } 2153 } 2154 for(i=0;i<4;i++) { 2155 for(j=0;j<9;j++) 2156 gfc->nsPsy.last_en_subshort[i][j] = 100; 2157 for(j=0;j<3;j++) 2158 gfc->nsPsy.last_attacks[i][j] = 0; 2159 gfc->nsPsy.pe_l[i] = gfc->nsPsy.pe_s[i] = 0; 2160 } 2161 2162 2163 2164 2165 /* gfp->cwlimit = sfreq*j/1024.0; */ 2166 gfc->cw_lower_index=6; 2167 if (gfp->cwlimit>0) 2168 cwlimit=gfp->cwlimit; 2169 else 2170 cwlimit=(FLOAT)8871.7; 2171 gfc->cw_upper_index = cwlimit*1024.0/gfp->out_samplerate; 2172 gfc->cw_upper_index=Min(HBLKSIZE-4,gfc->cw_upper_index); /* j+3 < HBLKSIZE-1 */ 2173 gfc->cw_upper_index=Max(6,gfc->cw_upper_index); 2174 2175 for ( j = 0; j < HBLKSIZE; j++ ) 2176 gfc->cw[j] = 0.4f; 2177 2178 2179 2180 i=L3para_read( gfp,(FLOAT8) samplerate,gfc->numlines_l,gfc->numlines_s, 2181 gfc->minval,gfc->s3_l,gfc->s3_s,gfc->SNR_s,gfc->bu_l, 2182 gfc->bo_l,gfc->w1_l,gfc->w2_l, gfc->bu_s,gfc->bo_s, 2183 gfc->w1_s,gfc->w2_s,&gfc->npart_l_orig,&gfc->npart_l, 2184 &gfc->npart_s_orig,&gfc->npart_s ); 2185 if (i!=0) return -1; 2186 2187 2188 2189 /* npart_l_orig = number of partition bands before convolution */ 2190 /* npart_l = number of partition bands after convolution */ 2191 2192 for (i=0; i<gfc->npart_l; i++) { 2193 for (j = 0; j < gfc->npart_l_orig; j++) { 2194 if (gfc->s3_l[i][j] != 0.0) 2195 break; 2196 } 2197 gfc->s3ind[i][0] = j; 2198 2199 for (j = gfc->npart_l_orig - 1; j > 0; j--) { 2200 if (gfc->s3_l[i][j] != 0.0) 2201 break; 2202 } 2203 gfc->s3ind[i][1] = j; 2204 } 2205 2206 2207 for (i=0; i<gfc->npart_s; i++) { 2208 for (j = 0; j < gfc->npart_s_orig; j++) { 2209 if (gfc->s3_s[i][j] != 0.0) 2210 break; 2211 } 2212 gfc->s3ind_s[i][0] = j; 2213 2214 for (j = gfc->npart_s_orig - 1; j > 0; j--) { 2215 if (gfc->s3_s[i][j] != 0.0) 2216 break; 2217 } 2218 gfc->s3ind_s[i][1] = j; 2219 } 2220 2221 2222 /* 2223 #include "debugscalefac.c" 2224 */ 2225 2226 2227 2228 if (gfc->nsPsy.use) { 2229 /* long block spreading function normalization */ 2230 for ( b = 0;b < gfc->npart_l; b++ ) { 2231 for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ ) { 2232 // spreading function has been properly normalized by 2233 // multiplying by DELBARK/.6609193 = .515. 2234 // It looks like Naoki was 2235 // way ahead of me and added this factor here! 2236 // it is no longer needed. 2237 //gfc->s3_l[b][k] *= 0.5; 2238 } 2239 } 2240 /* short block spreading function normalization */ 2241 // no longer needs to be normalized, but nspsytune wants 2242 // SNR_s applied here istead of later to save CPU cycles 2243 for ( b = 0;b < gfc->npart_s; b++ ) { 2244 FLOAT8 norm=0; 2245 for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) { 2246 norm += gfc->s3_s[b][k]; 2247 } 2248 for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) { 2249 gfc->s3_s[b][k] *= gfc->SNR_s[b] /* / norm */; 2250 } 2251 } 2252 } 2253 2254 2255 2256 if (gfc->nsPsy.use) { 2257 #if 1 2258 /* spread only from npart_l bands. Normally, we use the spreading 2259 * function to convolve from npart_l_orig down to npart_l bands 2260 */ 2261 for(b=0;b<gfc->npart_l;b++) 2262 if (gfc->s3ind[b][1] > gfc->npart_l-1) gfc->s3ind[b][1] = gfc->npart_l-1; 2263 #endif 2264 } 2265 2266 return 0; 2267 } 2268 2269 2270 2271 2272 2273 2274