lib/bladeenc/l3psy.c

   1 /*
   2                         (c) Copyright 1998-2001 - Tord Jansson
   3                         ======================================
   4
   5                 This file is part of the BladeEnc MP3 Encoder, based on
   6                 ISO's reference code for MPEG Layer 3 compression, and might
   7                 contain smaller or larger sections that are directly taken
   8                 from ISO's reference code.
   9
  10                 All changes to the ISO reference code herein are either
  11                 copyrighted by Tord Jansson (tord.jansson@swipnet.se)
  12                 or sublicensed to Tord Jansson by a third party.
  13
  14         BladeEnc is free software; you can redistribute this file
  15         and/or modify it under the terms of the GNU Lesser General Public
  16         License as published by the Free Software Foundation; either
  17         version 2.1 of the License, or (at your option) any later version.
  18
  19
  20
  21         ------------    Changes    ------------
  22
  23         2000-12-05  Andre Piotrowski
  24
  25         -       speed up: implemented prepacking of fft-data
  26
  27         2000-12-11  ap
  28
  29         -       speed up:  faster psycho_anal()
  30         -       optional bug fix: integrated better norm calclulation and block selecting
  31
  32         2000-12-12  ap
  33
  34         -       use SHORT_FFT_MIN_IDX to allow switching of "ORG_SHORT_CW_LIMIT" in "l3psy.h"
  35
  36         2001-01-12  ap
  37
  38         -       use some explicit type casting to avoid compiler warnings
  39 */
  40
  41 #define         RING_BUFFER                             1
  42
  43
  44
  45
  46
  47 #include        <stdlib.h>
  48
  49 #include        "common.h"
  50 #include        "tables.h"
  51 #include        "encoder.h"
  52 #include        "l3side.h"
  53 #include        "l3psy.h"
  54
  55
  56
  57
  58
  59 /*____ Global Static Variables ______________________________________________*/
  60
  61 /* The static variables "r", "phi_sav", "new_", "old" and "oldest" have    */
  62 /* to be remembered for the unpredictability measure.  For "r" and        */
  63 /* "phi_sav", the first index from the left is the channel select and     */
  64 /* the second index is the "age" of the data.                             */
  65
  66
  67 static  int                             new_, old, oldest;
  68 static  int                             flush, sync_flush, syncsize;
  69
  70 #if RING_BUFFER==1
  71 static  int                             savebuf_start_idx[2];
  72 #endif
  73
  74
  75
  76 #if NEW_L3PARM_TABLES
  77
  78 static  double                  *minval, *qthr_l;
  79 static  double                  *qthr_s, *SNR_s;
  80 static  int                             *cbw_l, *bu_l, *bo_l;
  81 static  int                             *cbw_s, *bu_s, *bo_s;
  82 static  double                  *w1_l, *w2_l;
  83 static  double                  *w1_s, *w2_s;
  84
  85 #if ORG_NUMLINES_NORM
  86
  87 static  int                             cbmax_l = CBANDS, cbmax_s = CBANDS_s;
  88 static  int                             numlines_l  [CBANDS];
  89
  90 static  int                             partition_l [HBLKSIZE];
  91 static  int                             partition_s [HBLKSIZE_s];
  92 static  double                  s3_l        [CBANDS][CBANDS];
  93 static  double                  *norm_l, *norm_s;
  94
  95 #else
  96
  97 static  int                             cbmax_l, cbmax_s;
  98 static  int                             *numlines_l;
  99 static  int                             *numlines_s;
 100
 101                                                 /* the non-zero entries of norm_l[i] * s3_l[i][j] */
 102 static  FLOAT                   normed_s3_l [900];   /* a bit more space than needed [799|855|735] */
 103 static  int                             lo_s3_l     [CBANDS];
 104 static  int                             hi_s3_l         [CBANDS];
 105
 106 static  FLOAT                   normed_s3_s [500];   /* a bit more space than needed [445|395|378] */
 107 static  int                             lo_s3_s     [CBANDS_s];
 108 static  int                             hi_s3_s         [CBANDS_s];
 109
 110 #endif          /* ORG_NUMLINES_NORM */
 111
 112 #else
 113
 114 static  double                  minval[CBANDS], qthr_l[CBANDS], norm_l[CBANDS];
 115 static  double                  qthr_s[CBANDS_s], norm_s[CBANDS_s], SNR_s[CBANDS_s];
 116 static  int                             cbw_l[SBMAX_l],bu_l[SBMAX_l],bo_l[SBMAX_l];
 117 static  int                             cbw_s[SBMAX_s],bu_s[SBMAX_s],bo_s[SBMAX_s];
 118 static  double                  w1_l[SBMAX_l], w2_l[SBMAX_l];
 119 static  double                  w1_s[SBMAX_s], w2_s[SBMAX_s];
 120
 121 static  int                             numlines_l  [CBANDS];
 122
 123 static  int                             partition_l [HBLKSIZE];
 124 static  int                             partition_s [HBLKSIZE_s];
 125 static  double                  s3_l        [CBANDS][CBANDS];
 126
 127 #endif          /* NEW_L3PARM_TABLES */
 128
 129
 130
 131 /* Scale Factor Bands */
 132 static  int                             blocktype_old[2];
 133
 134
 135
 136 static  double                  nb_1        [2][CBANDS];
 137 static  double                  nb_2        [2][CBANDS];
 138
 139 static  double                  cw          [HBLKSIZE];
 140
 141 static  FLOAT                   window      [BLKSIZE];
 142 static  FLOAT                   r                       [2][2][6];
 143 static  FLOAT                   phi_sav         [2][2][6];
 144
 145 static  FLOAT                   window_s    [BLKSIZE_s];
 146
 147 static  double                  ratio       [2][SBMAX_l];
 148 static  double                  ratio_s     [2][SBMAX_s][3];
 149
 150
 151
 152
 153
 154 #if NEW_L3PARM_TABLES
 155
 156 static  void                    L3para_read (int sfreq);
 157
 158 #if !ORG_NUMLINES_NORM
 159 static  void                    calc_normed_spreading
 160 (
 161         int                                             cbmax,                  /* number of lines and rows           */
 162         const double                    bval[],                 /* input values to compute the matrix */
 163         FLOAT                                   s3_ptr[],               /* the resulting non-zero entries     */
 164         int                                             lo_s3[],
 165         int                                             hi_s3[],
 166         const double                    norm[]
 167 );
 168 #endif
 169
 170 #else
 171
 172 static  void                    L3para_read
 173 (
 174         int                                             sfreq,
 175         int                                             numlines_l[CBANDS],
 176         int                                             partition_l[HBLKSIZE],
 177         double                                  minval[CBANDS],
 178         double                                  qthr_l[CBANDS],
 179         double                                  norm_l[CBANDS],
 180         double                                  s3_l[CBANDS][CBANDS],
 181         int                                             partition_s[HBLKSIZE_s],
 182         double                                  qthr_s[CBANDS_s],
 183         double                                  norm_s[CBANDS_s],
 184         double                                  SNR_s[CBANDS_s],
 185         int                                             cbw_l[SBMAX_l],
 186         int                                             bu_l[SBMAX_l],
 187         int                                             bo_l[SBMAX_l],
 188         double                                  w1_l[SBMAX_l],
 189         double                                  w2_l[SBMAX_l],
 190         int                                             cbw_s[SBMAX_s],
 191         int                                             bu_s[SBMAX_s],
 192         int                                             bo_s[SBMAX_s],
 193         double                                  w1_s[SBMAX_s],
 194         double                                  w2_s[SBMAX_s]
 195 );
 196
 197 #endif
 198
 199
 200
 201
 202
 203 /*____ psycho_anal_init() ___________________________________________________*/
 204
 205 void                                    psycho_anal_init (double sfreq)
 206 {
 207         unsigned int                    ch, sfb, b, i, j;
 208
 209
 210         /* reset the r, phi_sav "ring buffer" indices */
 211
 212         old = 1 - (new_ = oldest = 0);
 213
 214
 215         /* clear the ratio arrays */
 216
 217         for (ch = 0;  ch < 2;  ch++)
 218         {
 219                 for (sfb = 0;  sfb < SBMAX_l;  sfb++)
 220                         ratio[ch][sfb] = 0.0;
 221
 222                 for (sfb = 0;  sfb < SBMAX_s;  sfb++)
 223                         for (b = 0;  b < 3;  b++)
 224                                 ratio_s[ch][sfb][b] = 0.0;
 225         }
 226
 227
 228         /* clear preecho arrays */
 229
 230         for (ch = 0;  ch < 2;  ch++)
 231         {
 232                 for (i = 0;  i < CBANDS;  i++)
 233                 {
 234                         nb_1[ch][i] = 0;
 235                         nb_2[ch][i] = 0;
 236                 }
 237         }
 238
 239
 240         /* clear blocktype information */
 241
 242         for (ch = 0;  ch < 2;  ch++)
 243                 blocktype_old[ch] = NORM_TYPE;
 244
 245
 246         sync_flush =  768;
 247         flush      =  576;
 248         syncsize   = 1344;   /* sync_flush + flush */
 249
 250
 251 #if RING_BUFFER==1
 252         for (ch = 0;  ch < 2;  ch++)
 253                 savebuf_start_idx[ch] = 0;
 254 #endif
 255
 256
 257         /* calculate HANN window coefficients */
 258
 259     for (i = 0;  i < BLKSIZE;  i++)
 260                 window[i] = (FLOAT) (0.5 * (1 - cos (2.0 * PI * (i - 0.5) / BLKSIZE)));
 261
 262     for (i = 0;  i < BLKSIZE_s;  i++)
 263                 window_s[i] = (FLOAT) (0.5 * (1 - cos (2.0 * PI * (i - 0.5) / BLKSIZE_s)));
 264
 265
 266         /* reset states used in unpredictability measure */
 267
 268         for (ch = 0;  ch < 2;  ch++)
 269         {
 270                 for (i = 0;  i < 2;  i++)
 271                 {
 272                         for (j = 0;  j < 6;  j++)
 273                         {
 274                                       r[ch][i][j] = 0.0;
 275                                 phi_sav[ch][i][j] = 0.0;
 276                         }
 277                 }
 278         }
 279
 280
 281 #if NEW_L3PARM_TABLES
 282         L3para_read ((int) sfreq);
 283 #else
 284         L3para_read
 285         (
 286                 (int) sfreq,
 287                 numlines_l, partition_l, minval, qthr_l, norm_l, s3_l,
 288                 partition_s, qthr_s, norm_s, SNR_s,
 289                 cbw_l, bu_l, bo_l, w1_l, w2_l,
 290                 cbw_s, bu_s, bo_s, w1_s, w2_s
 291         );
 292 #endif
 293
 294
 295         /* Set unpredicatiblility of remaining spectral lines to 0.4 */
 296
 297         for (j = 206;  j < HBLKSIZE;  j++)
 298                 cw[j] = 0.4;
 299 }
 300
 301
 302
 303
 304
 305 /*____ psycho_anal_exit() ___________________________________________________*/
 306
 307 void psycho_anal_exit( void )
 308 {
 309         /* nothing to do */
 310 }
 311
 312
 313
 314
 315
 316 /*____ psycho_anal() ________________________________________________________*/
 317
 318 void                                    psycho_anal
 319 (
 320 #if ORG_BUFFERS
 321         short int                               *buffer,
 322         short int                               savebuf[2048],
 323 #else
 324         FLOAT                                   *buffer,
 325         int                                             buffer_idx,
 326 #endif
 327         int                                             ch,
 328         int                                             lay,
 329 /*      float                                   snr32[32], */
 330         double                                  ratio_d[SBMAX_l],
 331         double                                  ratio_ds[SBMAX_s][3],
 332         double                                  *pe,
 333         gr_info                                 *cod_info
 334 )
 335 {
 336         int                                             blocktype;
 337         unsigned int                    sfb, b, j, k;
 338         double                                  r_prime, phi_prime; /* not FLOAT */
 339         double                                  temp1, temp2, temp3;
 340
 341 #if !ORG_NUMLINES_NORM && NEW_L3PARM_TABLES
 342         FLOAT                                   *s3_ptr;
 343 #endif
 344
 345         int                                             sblock;
 346
 347         double                                  thr         [CBANDS];
 348         double                                  eb          [CBANDS];
 349         FLOAT                                   cb          [CBANDS];
 350         FLOAT                                   wsamp_r     [HBLKSIZE];
 351         FLOAT                                   wsamp_i     [HBLKSIZE];
 352
 353         FLOAT                                   energy      [HBLKSIZE];
 354         FLOAT                                   phi         [6];
 355         FLOAT                                   energy_s    [3][BLKSIZE_s];
 356         FLOAT                                   phi_s       [3][52];
 357
 358 #if ORG_BUFFERS
 359 #if RING_BUFFER==1
 360         int                                             beg, idx, fin;
 361 #endif
 362 #else
 363 #       define                                  savebuf         buffer
 364 #       define                                  beg                     buffer_idx
 365         int                                             idx, fin;
 366 #endif
 367
 368
 369         for (sfb = 0;  sfb < SBMAX_l;  sfb++)
 370                 ratio_d[sfb] = ratio[ch][sfb];
 371
 372         for (sfb = 0;  sfb < SBMAX_s;  sfb++)
 373                 for (b = 0;  b < 3;  b++)
 374                         ratio_ds[sfb][b] = ratio_s[ch][sfb][b];
 375
 376
 377         if (ch == 0)
 378                 old = 1 - (new_ = oldest = old);
 379
 380
 381 #if ORG_BUFFERS
 382         /**********************************************************************
 383         *  Delay signal by sync_flush=768 samples                             *
 384         **********************************************************************/
 385
 386 #       if RING_BUFFER==0
 387                 for (j = 0;  j < sync_flush;  j++)   /* for long window samples */
 388                         savebuf[j] = savebuf[j+flush];
 389
 390                 for (j = sync_flush;  j < syncsize;  j++)
 391                         savebuf[j] = *buffer++;
 392 #       else
 393                 beg = savebuf_start_idx[ch] = (savebuf_start_idx[ch] + flush) & 2047;
 394
 395                 idx = (beg + sync_flush) & 2047;
 396                 fin = (idx + flush) & 2047;
 397                 if (idx >= fin)
 398                 {
 399                         while (idx < 2048)
 400                                 savebuf[idx++] = *buffer++;
 401                         idx = 0;
 402                 }
 403                 while (idx < fin)
 404                         savebuf[idx++] = *buffer++;
 405 #       endif
 406 #endif
 407
 408
 409 /**********************************************************************
 410 *    compute unpredicatability of first six spectral lines            *
 411 **********************************************************************/
 412
 413 #if RING_BUFFER==0
 414         for (j = 0, k = 0, idx = 0;  j < BLKSIZE/2;  j++)
 415         {
 416                 wsamp_r[j] = window[k++] * savebuf[idx++];
 417                 wsamp_i[j] = window[k++] * savebuf[idx++];
 418         }
 419 #else
 420         j = 0;  k = 0;
 421         idx = beg;
 422         fin = (idx + BLKSIZE) & 2047;
 423         if (idx >= fin)
 424         {
 425                 while (idx < 2048)
 426                 {
 427                         wsamp_r[j] = window[k++] * savebuf[idx++];
 428                         wsamp_i[j] = window[k++] * savebuf[idx++];
 429                         j++;
 430                 }
 431                 idx = 0;
 432         }
 433         while (idx < fin)
 434         {
 435                 wsamp_r[j] = window[k++] * savebuf[idx++];
 436                 wsamp_i[j] = window[k++] * savebuf[idx++];
 437                 j++;
 438         }
 439 #endif
 440
 441         fft(wsamp_r, wsamp_i, energy, phi, BLKSIZE);   /* long FFT */
 442
 443         for (j = 0;  j < 6;  j++)
 444         {       /* calculate unpredictability measure cw */
 445                 double r1, phi1;
 446                                      r_prime = 2.0 *       r[ch][old][j] -       r[ch][oldest][j];
 447                                    phi_prime = 2.0 * phi_sav[ch][old][j] - phi_sav[ch][oldest][j];
 448                       r[ch][new_][j] = (FLOAT) (  r1 = sqrt((double) energy[j]));
 449                 phi_sav[ch][new_][j] = (FLOAT) (phi1 =                  phi[j] );
 450
 451                 temp3 = r1 + fabs(r_prime);
 452                 if (temp3 != 0.0)
 453                 {
 454                         temp1 = r1 * cos(phi1) - r_prime * cos(phi_prime);
 455                         temp2 = r1 * sin(phi1) - r_prime * sin(phi_prime);
 456                         cw[j] = sqrt(temp1*temp1 + temp2*temp2) / temp3;
 457                 }
 458                 else
 459                         cw[j] = 0;
 460         }
 461
 462
 463 /**********************************************************************
 464 *     compute unpredicatibility of next 200 spectral lines            *
 465 **********************************************************************/
 466
 467         for (b = 0;  b < 3;  b++)
 468         {
 469 #if RING_BUFFER==0
 470                 for (j = 0, k = 0, idx = 128*(2 + b);  j < BLKSIZE_s/2;  j++)
 471                 {       /* window data with HANN window */
 472                         wsamp_r[j] = window_s[k++] * savebuf[idx++];
 473                         wsamp_i[j] = window_s[k++] * savebuf[idx++];
 474                 }
 475 #else
 476                 j = 0;  k = 0;
 477                 idx = (beg + 128*(2 + b)) & 2047;
 478                 fin = (idx + BLKSIZE_s) & 2047;
 479                 if (idx >= fin)
 480                 {
 481                         while (idx < 2048)
 482                         {
 483                                 wsamp_r[j] = window_s[k++] * savebuf[idx++];
 484                                 wsamp_i[j] = window_s[k++] * savebuf[idx++];
 485                                 j++;
 486                         }
 487                         idx = 0;
 488                 }
 489                 while (idx < fin)
 490                 {
 491                         wsamp_r[j] = window_s[k++] * savebuf[idx++];
 492                         wsamp_i[j] = window_s[k++] * savebuf[idx++];
 493                         j++;
 494                 }
 495 #endif
 496
 497                 fft (wsamp_r, wsamp_i, energy_s[b], phi_s[b], BLKSIZE_s);   /* short FFT*/
 498         }
 499
 500         for (j = 6, k = SHORT_FFT_MIN_IDX;  j < 206;  j += 4, k++)
 501         {       /* calculate unpredictability measure cw */
 502                 double r1, phi1;
 503
 504                   r_prime = 2.0 * sqrt((double) energy_s[0][k]) - sqrt((double) energy_s[2][k]);
 505                 phi_prime = 2.0 *                  phi_s[0][k]  -                  phi_s[2][k];
 506                        r1 = sqrt((double) energy_s[1][k]);
 507                      phi1 =                  phi_s[1][k];
 508
 509                 temp3 = r1 + fabs(r_prime);
 510                 if (temp3 != 0.0)
 511                 {
 512                         temp1 = r1 * cos(phi1) - r_prime * cos(phi_prime);
 513                         temp2 = r1 * sin(phi1) - r_prime * sin(phi_prime);
 514                         cw[j] = sqrt(temp1*temp1 + temp2*temp2) / temp3;
 515                 }
 516                 else
 517                         cw[j] = 0.0;
 518
 519                 cw[j+1] = cw[j+2] = cw[j+3] = cw[j];
 520         }
 521
 522
 523 /**********************************************************************
 524 *    Calculate the energy and the unpredictability in the threshold   *
 525 *    calculation partitions                                           *
 526 **********************************************************************/
 527
 528 #if ORG_NUMLINES_NORM || !NEW_L3PARM_TABLES
 529
 530         for (b = 0;  b < cbmax_l;  b++)
 531         {
 532                 eb[b] = 0.0;
 533                 cb[b] = 0.0;
 534         }
 535         for (j = 0;  j < HBLKSIZE;  j++)
 536         {
 537                 int tp = partition_l[j];
 538                 if (tp >= 0)
 539                 {
 540                         eb[tp] += energy[j];
 541                         cb[tp] += cw[j] * energy[j];
 542                 }
 543         }
 544
 545 #else
 546
 547         j = 0;
 548         for (b = 0;  b < cbmax_l;  b++)
 549         {
 550                 eb[b] = 0.0;
 551                 cb[b] = 0.0;
 552
 553                 /*
 554                         Calculate the energy and the unpredictability in the threshold
 555                         calculation partitions
 556
 557                         cbmax_l holds the number of valid numlines_l entries
 558                 */
 559                 k = numlines_l[b];
 560                 do {
 561                         eb[b] += energy[j];
 562                         cb[b] += cw[j] * energy[j];
 563                 } while (j++, --k);
 564         }
 565
 566         s3_ptr = normed_s3_l;
 567
 568 #endif
 569
 570
 571         *pe = 0.0;
 572
 573         for (b = 0;  b < cbmax_l;  b++)
 574         {
 575                 FLOAT                                   nb;
 576                 FLOAT                                   ecb = 0.0;
 577                 double                                  ctb = 0.0;
 578                 double                                  SNR_l;
 579                 double                                  cbb, tbb;
 580
 581
 582                 /*
 583                         convolve the partitioned energy and unpredictability
 584                         with the spreading function, normed_s3_l[b][k]
 585                 */
 586 #if ORG_NUMLINES_NORM || !NEW_L3PARM_TABLES
 587                 for (k = 0;  k < cbmax_l;  k++)
 588                 {
 589                         ecb += s3_l[b][k] * eb[k];  /* sprdngf for Layer III */
 590                         ctb += s3_l[b][k] * cb[k];
 591                 }
 592 #else
 593                 for (k = lo_s3_l[b];  k < hi_s3_l[b];  k++)
 594                 {
 595                         ecb += *s3_ptr   * eb[k];  /* sprdngf for Layer III */
 596                         ctb += *s3_ptr++ * cb[k];
 597                 }
 598 #endif
 599
 600
 601                 /*
 602                         calculate the tonality of each threshold calculation partition
 603                         calculate the SNR in each threshhold calculation partition
 604                 */
 605                 if (ecb != 0.0)
 606                 {
 607                         cbb = ctb / ecb;
 608                         if (cbb < 0.01)
 609                                 cbb = 0.01;
 610                         tbb = -0.299 - 0.43 * log(cbb);   /* conv1=-0.299, conv2=-0.43 */
 611                         tbb = MIN(MAX (0.0, tbb), 1.0) ;  /* 0<=tbb<=1 */
 612                 }
 613                 else
 614                         tbb = 0.0;  /* cbb==0 => -0.299-0.43*cbb<0 => tbb=0*/
 615
 616                 /* TMN=29.0,NMT=6.0 for all calculation partitions */
 617                 SNR_l = MAX (minval[b], 23.0 * tbb + 6.0);   /* 29*tbb + 6*(1-tbb) */
 618
 619                 /* calculate the threshold for each partition */
 620 #if ORG_NUMLINES_NORM || !NEW_L3PARM_TABLES
 621             nb = ecb * norm_l[b] * exp(-SNR_l * LN_TO_LOG10);
 622 #else
 623             nb = ecb * exp(-SNR_l * LN_TO_LOG10);   /* our ecb is already normed */
 624 #endif
 625
 626                 /*
 627                         pre-echo control
 628                 */
 629                 thr[b] = MAX (qthr_l[b], MIN(nb, nb_2[ch][b]));
 630                 nb_2[ch][b] = MIN(2.0 * nb, 16.0 * nb_1[ch][b]);
 631             nb_1[ch][b] = nb;
 632
 633
 634                 /*
 635                         calculate percetual entropy
 636
 637                         thr[b] -> thr[b]+1.0 : for non sound portition
 638                 */
 639                 if (eb[b] > thr[b])
 640                         *pe += numlines_l[b] * log((eb[b]+1.0) / (thr[b]+1.0));
 641         }
 642
 643
 644 #define switch_pe  1800
 645
 646
 647         if (*pe < switch_pe)
 648         {
 649                 /* no attack : use long blocks */
 650
 651                 if (blocktype_old[ch] == SHORT_TYPE)
 652                         blocktype = STOP_TYPE;
 653                 else   /* NORM_TYPE, STOP_TYPE */
 654                         blocktype = NORM_TYPE;
 655
 656
 657                 /* threshold calculation (part 2) */
 658
 659                 for (sfb = 0;  sfb < SBMAX_l;  sfb++)
 660                 {
 661                         int             bu = bu_l[sfb];
 662                         int             bo = bo_l[sfb];
 663                         double  en = w1_l[sfb] * eb[bu] + w2_l[sfb] * eb[bo];
 664
 665                         for (b = bu+1;  b < bo;  b++)
 666                                 en += eb[b];
 667
 668                         if (en != 0.0)
 669                         {
 670                                 double  thm = w1_l[sfb] * thr[bu] + w2_l[sfb] * thr[bo];
 671
 672                                 for (b = bu+1;  b < bo;  b++)
 673                                         thm += thr[b];
 674
 675                                 ratio[ch][sfb] = thm / en;
 676                         }
 677                         else
 678                                 ratio[ch][sfb] = 0.0;
 679                 }
 680         }
 681         else
 682         {
 683                 /* attack : use short blocks */
 684                 blocktype = SHORT_TYPE;
 685 #if ORG_BLOCK_SELECT
 686                 if (blocktype_old[ch] == NORM_TYPE)
 687                         blocktype_old[ch] = START_TYPE;
 688                 else   /* SHORT_TYPE, STOP_TYPE */
 689                         blocktype_old[ch] = SHORT_TYPE;
 690 #else   /* ISO */
 691                 if (blocktype_old[ch] == SHORT_TYPE)
 692                         blocktype_old[ch] = SHORT_TYPE;
 693                 else   /* NORM_TYPE, STOP_TYPE */
 694                         blocktype_old[ch] = START_TYPE;
 695 #endif
 696
 697
 698                 /* threshold calculation for short blocks */
 699
 700                 for (sblock = 0;  sblock < 3;  sblock++)
 701                 {
 702 #if ORG_NUMLINES_NORM || !NEW_L3PARM_TABLES
 703
 704                         for (b = 0;  b < cbmax_s;  b++)
 705                                 eb[b] = 0.0;
 706
 707                         for (j = 0;  j < HBLKSIZE_s;  j++)
 708                                 eb[partition_s[j]] += energy_s[sblock][j];
 709
 710 #else
 711
 712                         j = 0;
 713                         for (b = 0;  b < cbmax_s;  b++)
 714                         {
 715                                 eb[b] = 0.0;
 716
 717                                 /*
 718                                         Calculate the energy and the unpredictability in the threshold
 719                                         calculation partitions
 720
 721                                         cbmax_s holds the number of valid numlines_s entries
 722                                 */
 723                                 k = numlines_s[b];
 724                                 do {
 725                                         eb[b] += energy_s[sblock][j];
 726                                 } while (j++, --k);
 727                         }
 728
 729                         s3_ptr = normed_s3_s;
 730 #endif
 731
 732                         for (b = 0;  b < cbmax_s;  b++)
 733                         {
 734                                 FLOAT                                   nb;
 735                                 FLOAT                                   ecb = 0.0;
 736
 737 #if ORG_NUMLINES_NORM || !NEW_L3PARM_TABLES
 738                                 for (k = 0;  k < cbmax_s;  k++)
 739                                         ecb += s3_l[b][k] * eb[k];
 740
 741                                 nb = ecb * norm_l[b] * exp((double) SNR_s[b] * LN_TO_LOG10);
 742 #else
 743                                 for (k = lo_s3_s[b];  k < hi_s3_s[b];  k++)
 744                                         ecb += *s3_ptr++ * eb[k];
 745
 746                                 nb = ecb * exp((double) SNR_s[b] * LN_TO_LOG10);   /* our ecb is already normed */
 747 #endif
 748                                 thr[b] = MAX(qthr_s[b], nb);
 749                         }
 750
 751                         for (sfb = 0;  sfb < SBMAX_s;  sfb++)
 752                         {
 753                                 int             bu = bu_s[sfb];
 754                                 int             bo = bo_s[sfb];
 755                                 double  en = w1_s[sfb] * eb[bu] + w2_s[sfb] * eb[bo];
 756
 757                                 for (b = bu+1;  b < bo;  b++)
 758                                         en += eb[b];
 759                                 if (en != 0.0)
 760                                 {
 761                                         double  thm = w1_s[sfb] * thr[bu] + w2_s[sfb] * thr[bo];
 762
 763                                         for (b = bu+1;  b < bo;  b++)
 764                                                 thm += thr[b];
 765
 766                                         ratio_s[ch][sfb][sblock] = thm / en;
 767                                 }
 768                                 else
 769                                         ratio_s[ch][sfb][sblock] = 0.0;
 770                         }
 771                 }
 772         }
 773
 774         cod_info->block_type = blocktype_old[ch];
 775         blocktype_old[ch] = blocktype;
 776
 777         if ( cod_info->block_type == NORM_TYPE )
 778             cod_info->window_switching_flag = 0;
 779         else
 780             cod_info->window_switching_flag = 1;
 781
 782         cod_info->mixed_block_flag = 0;
 783 }
 784
 785
 786
 787
 788
 789 /*____ L3para_read() __________________________________________________________*/
 790
 791 #if NEW_L3PARM_TABLES
 792
 793 static void                             L3para_read (int sfreq)
 794 {
 795         int                                             sfreq_idx;
 796         l3_parm_block                   *parm;
 797         double                                  *bval_l, *bval_s;
 798
 799 #if ORG_NUMLINES_NORM
 800         int                                             cbmax_l, cbmax_s;
 801         int                                             i, j, k;
 802 #else
 803         double                                  *norm_l, *norm_s;
 804 #endif
 805
 806
 807         /*
 808                 Set parameter block
 809         */
 810         switch (sfreq)
 811         {
 812                 case 32000:  sfreq_idx = 2;  break;
 813                 case 44100:  sfreq_idx = 0;  break;
 814                 case 48000:  sfreq_idx = 1;  break;
 815                 default   :  return;  /* Just to avoid compiler warnings */
 816         }
 817         parm = l3_parm + sfreq_idx;
 818
 819
 820         /*
 821                 Read long block data
 822         */
 823         cbmax_l    = parm->long_data.cbmax_l;
 824
 825 #if ORG_NUMLINES_NORM
 826         for (i = 0, j = 0;  i < cbmax_l;  i++)
 827         {
 828                 numlines_l[i] = parm->long_data.numlines_l[i];
 829
 830                 for (k = 0;  k < numlines_l[i];  k++)
 831                         partition_l[j++] = i;
 832         }
 833 #else
 834         numlines_l = parm->long_data.numlines_l;
 835 #endif
 836
 837         minval     = parm->long_data.minval;
 838         qthr_l     = parm->long_data.qthr_l;
 839         norm_l     = parm->long_data.norm_l;
 840         bval_l     = parm->long_data.bval_l;
 841
 842
 843         /*
 844                 Compute the normed spreading function norm_l[i] * s3_l[i][j]
 845         */
 846 #if ORG_NUMLINES_NORM
 847         for (i = 0;  i < cbmax_l;  i++)
 848         {
 849                 double x, temp, tempx, tempy;
 850
 851                 for (j = 0;  j < cbmax_l;  j++)
 852                 {
 853 /*                      tempx = (bval_l[i]-bval_l[j]) * 1.05; */
 854                         if (j >= i)
 855                                 tempx = (bval_l[i]-bval_l[j]) * 3.0;
 856                         else
 857                                 tempx = (bval_l[i]-bval_l[j]) * 1.5;
 858 /*                      if (j >= i)  tempx = (bval_l[j]-bval_l[i]) * 3.0;
 859                         else         tempx = (bval_l[j]-bval_l[i]) * 1.5; */
 860                         if (tempx > 0.5  &&  tempx < 2.5)
 861                         {
 862                                 temp = tempx - 0.5;
 863                                 x = 8.0 * temp * (temp-2.0);
 864                         }
 865                         else  x = 0.0;
 866                         tempx += 0.474;
 867                         tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
 868                         if (tempy <= -60.0)  s3_l[i][j] = 0.0;
 869                         else                 s3_l[i][j] = exp((x + tempy) * LN_TO_LOG10);
 870                 }
 871         }
 872 #else
 873         calc_normed_spreading (cbmax_l, bval_l, normed_s3_l, lo_s3_l, hi_s3_l, norm_l);
 874 #endif
 875
 876
 877         /*
 878                 Read short block data
 879         */
 880         cbmax_s    = parm->short_data.cbmax_s;
 881
 882 #if ORG_NUMLINES_NORM
 883         for (i = 0, j = 0;  i < cbmax_s;  i++)
 884         {
 885                 numlines_l[i] = parm->short_data.numlines_s[i];
 886
 887                 for (k = 0;  k < numlines_l[i];  k++)
 888                         partition_s[j++] = i;
 889         }
 890 #else
 891         numlines_s = parm->short_data.numlines_s;
 892 #endif
 893
 894         qthr_s     = parm->short_data.qthr_s;
 895         norm_s     = parm->short_data.norm_s;
 896         SNR_s      = parm->short_data.SNR_s;
 897         bval_s     = parm->short_data.bval_s;
 898
 899
 900 #if !ORG_NUMLINES_NORM
 901
 902         /*
 903                 Compute the normed spreading function norm_s[i] * s3_s[i][j]
 904         */
 905         calc_normed_spreading (cbmax_s, bval_s, normed_s3_s, lo_s3_s, hi_s3_s, norm_s);
 906
 907 #endif
 908
 909
 910         /*
 911                 Read long block data for converting threshold
 912                 calculation partitions to scale factor bands
 913         */
 914         cbw_l = parm->long_thres.cbw_l;
 915         bu_l  = parm->long_thres.bu_l;
 916         bo_l  = parm->long_thres.bo_l;
 917         w1_l  = parm->long_thres.w1_l;
 918         w2_l  = parm->long_thres.w2_l;
 919
 920
 921         /*
 922                 Read short block data for converting threshold
 923                 calculation partitions to scale factor bands
 924         */
 925         cbw_s = parm->short_thres.cbw_s;
 926         bu_s  = parm->short_thres.bu_s;
 927         bo_s  = parm->short_thres.bo_s;
 928         w1_s  = parm->short_thres.w1_s;
 929         w2_s  = parm->short_thres.w2_s;
 930 }
 931
 932 #else           /* NEW_L3PARM_TABLES */
 933
 934 static  void                    L3para_read
 935 (
 936         int                                             sfreq,
 937         int                                             numlines_l[CBANDS],
 938         int                                             partition_l[HBLKSIZE],
 939         double                                  minval[CBANDS],
 940         double                                  qthr_l[CBANDS],
 941         double                                  norm_l[CBANDS],
 942         double                                  s3_l[CBANDS][CBANDS],
 943         int                                             partition_s[HBLKSIZE_s],
 944         double                                  qthr_s[CBANDS_s],
 945         double                                  norm_s[CBANDS_s],
 946         double                                  SNR_s[CBANDS_s],
 947         int                                             cbw_l[SBMAX_l],
 948         int                                             bu_l[SBMAX_l],
 949         int                                             bo_l[SBMAX_l],
 950         double                                  w1_l[SBMAX_l],
 951         double                                  w2_l[SBMAX_l],
 952         int                                             cbw_s[SBMAX_s],
 953         int                                             bu_s[SBMAX_s],
 954         int                                             bo_s[SBMAX_s],
 955         double                                  w1_s[SBMAX_s],
 956         double                                  w2_s[SBMAX_s]
 957 )
 958 {
 959         static  double                  bval_l[CBANDS];
 960         int                                             cbmax_tp;
 961
 962         int                                             sbmax;
 963         int                                             i, j, k, k2;
 964
 965
 966         psyDataElem                             *rpa1;
 967         psyDataElem2                    *rpa2;
 968         psyDataElem3                    *rpa3;
 969
 970
 971 /* Read long block data */
 972
 973         switch (sfreq)
 974         {
 975                 case 32000:  rpa1 = psy_longBlock_32000_58;  cbmax_tp = 59;  break;
 976                 case 44100:  rpa1 = psy_longBlock_44100_62;  cbmax_tp = 63;  break;
 977                 case 48000:  rpa1 = psy_longBlock_48000_61;  cbmax_tp = 62;  break;
 978                 default   :  return;  /* Just to avoid compiler warnings */
 979         }
 980
 981         for (i = 0, k2 = 0;  i < cbmax_tp;  i++)
 982         {
 983                 numlines_l[i] = rpa1->lines;
 984                 minval[i]     = rpa1->minVal;
 985                 qthr_l[i]     = rpa1->qthr;
 986                 norm_l[i]     = rpa1->norm;
 987                 bval_l[i]     = rpa1->bVal;
 988                 rpa1++;
 989
 990                 for (k = 0;  k < numlines_l[i];  k++)
 991                         partition_l[k2++] = i;
 992         }
 993
 994
 995 /************************************************************************
 996  * Now compute the spreading function, s[j][i], the value of the spread-*
 997  * ing function, centered at band j, for band i, store for later use    *
 998  ************************************************************************/
 999
1000         for (i = 0;  i < cbmax_tp;  i++)
1001         {
1002                 double x, temp, tempx, tempy;
1003
1004                 for (j = 0;  j < cbmax_tp;  j++)
1005                 {
1006 /*                      tempx = (bval_l[i]-bval_l[j]) * 1.05; */
1007                         if (j >= i)
1008                                 tempx = (bval_l[i]-bval_l[j]) * 3.0;
1009                         else
1010                                 tempx = (bval_l[i]-bval_l[j]) * 1.5;
1011 /*                      if (j >= i)  tempx = (bval_l[j]-bval_l[i]) * 3.0;
1012                         else         tempx = (bval_l[j]-bval_l[i]) * 1.5; */
1013                         if (tempx > 0.5  &&  tempx < 2.5)
1014                         {
1015                                 temp = tempx - 0.5;
1016                                 x = 8.0 * temp * (temp-2.0);
1017                         }
1018                         else  x = 0.0;
1019                         tempx += 0.474;
1020                         tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
1021                         if (tempy <= -60.0)  s3_l[i][j] = 0.0;
1022                         else                 s3_l[i][j] = exp((x + tempy) * LN_TO_LOG10);
1023                 }
1024         }
1025
1026
1027 /* Read short block data */
1028
1029         switch (sfreq)
1030         {
1031                 case 32000:  rpa2 = psy_shortBlock_32000_41;  cbmax_tp = 42;  break;
1032                 case 44100:  rpa2 = psy_shortBlock_44100_38;  cbmax_tp = 39;  break;
1033                 case 48000:  rpa2 = psy_shortBlock_48000_37;  cbmax_tp = 38;  break;
1034                 default   :  return;  /* Just to avoid compiler warnings */
1035         }
1036
1037         for (i = 0, k2 = 0;  i < cbmax_tp;  i++)
1038         {
1039                 numlines_l[i] = rpa2->lines;
1040                 qthr_s[i]     = rpa2->qthr;
1041                 norm_s[i]     = rpa2->norm;
1042                 SNR_s[i]      = rpa2->snr;
1043                 rpa2++;
1044
1045                 for (k = 0;  k < numlines_l[i];  k++)
1046                         partition_s[k2++] = i;
1047         }
1048
1049
1050 /* Read long block data for converting threshold calculation
1051    partitions to scale factor bands */
1052
1053         switch (sfreq)
1054         {
1055                 case 32000:  rpa3 = psy_data3_32000_20;  break;
1056                 case 44100:  rpa3 = psy_data3_44100_20;  break;
1057                 case 48000:  rpa3 = psy_data3_48000_20;  break;
1058                 default   :  return;  /* Just to avoid compiler warnings */
1059         }
1060         sbmax = SBMAX_l;
1061
1062         for (i = 0;  i < sbmax;  i++)
1063         {
1064                 cbw_l[i] = rpa3->cbw;
1065                 bu_l[i] = rpa3->bu;
1066                 bo_l[i] = rpa3->bo;
1067                 w1_l[i] = rpa3->w1;
1068                 w2_l[i] = rpa3->w2;
1069                 rpa3++;
1070         }
1071
1072
1073 /* Read short block data for converting threshold calculation
1074    partitions to scale factor bands */
1075
1076         switch (sfreq)
1077         {
1078                 case 32000:  rpa3 = psy_data4_32000_11;  break;
1079                 case 44100:  rpa3 = psy_data4_44100_11;  break;
1080                 case 48000:  rpa3 = psy_data4_48000_11;  break;
1081                 default   :  return;  /* Just to avoid compiler warnings */
1082         }
1083         sbmax = SBMAX_s;
1084
1085         for (i = 0;  i < sbmax;  i++)
1086         {
1087                 cbw_s[i] = rpa3->cbw;
1088                 bu_s[i] = rpa3->bu;
1089                 bo_s[i] = rpa3->bo;
1090                 w1_s[i] = rpa3->w1;
1091                 w2_s[i] = rpa3->w2;
1092                 rpa3++;
1093         }
1094 }
1095
1096 #endif          /* NEW_L3PARM_TABLES */
1097
1098
1099
1100
1101
1102 #if !ORG_NUMLINES_NORM && NEW_L3PARM_TABLES
1103
1104 /*  ========================================================================================  */
1105 /*              calc_normed_spreading                                                         */
1106 /*  ========================================================================================  */
1107 /*
1108         Compute the normed spreading function,
1109         the normed value of the spreading function,
1110         centered at band j, for band i, store for later use
1111
1112         Since this is a band matrix, we store only the non-zero entries
1113         in linear order in the single dimension array normed_s3.
1114
1115         The array has to be accessed in linear order, too, starting with line 0,
1116         up to line cbmax-1. For line b, the current entries represent
1117
1118                 norm[b] * s3[b][lo_s3[b]]  ...  norm[b] * s3[b][hi_s3[b]-1]
1119
1120         Normally, we could easily compute the norm [building the reciprocal of the line sum].
1121         Alas, dist10 uses somewhat (strange and) different, that made our norm differring too
1122         much at the last few lines. Thus, we renounce and use the original values.
1123 */
1124
1125 static  void                    calc_normed_spreading
1126 (
1127         int                                             cbmax,                  /* number of lines and rows           */
1128         const double                    bval[],                 /* input values to compute the matrix */
1129         FLOAT                                   s3_ptr[],               /* the resulting non-zero entries     */
1130         int                                             lo_s3[],
1131         int                                             hi_s3[],
1132         const double                    norm[]
1133 )
1134 {
1135         double                                  arg, x, y;
1136         double                                  s3[CBANDS];
1137         int                                             i, j;
1138         int                                             non_zero_part;
1139
1140
1141
1142         for (i = 0;  i < cbmax;  i++)
1143         {
1144                 non_zero_part = FALSE;
1145                 hi_s3[i] = cbmax;   /* we preset this value for the case that the line ends with a non-zero entry */
1146
1147                 for (j = 0;  j < cbmax;  j++)
1148                 {
1149                         if (j >= i)
1150                                 arg = (bval[i] - bval[j]) * 3.0;
1151                         else
1152                                 arg = (bval[i] - bval[j]) * 1.5;
1153
1154                         if (arg > 0.5  &&  arg < 2.5)
1155                                 x = 8.0 * (arg - 0.5) * (arg - 2.5);
1156                         else
1157                                 x = 0.0;
1158
1159                         arg += 0.474;
1160
1161                         y = 15.811389 + 7.5 * arg - 17.5 * sqrt(1.0 + arg * arg);
1162
1163                         if (y <= -60.0)
1164                         {
1165                                 if (non_zero_part)   /* only zeroes will follow */
1166                                 {
1167                                         hi_s3[i] = j;
1168                                         break;   /* so cut the computing for this line */
1169                                 }
1170                         }
1171                         else
1172                         {
1173                                 s3[j] = exp((x + y) * LN_TO_LOG10);
1174
1175                                 if (! non_zero_part)
1176                                 {
1177                                         lo_s3[i] = j;
1178                                         non_zero_part = TRUE;   /* the first non-zero entry ends the non_zero_part */
1179                                 }
1180                         }
1181                 }
1182
1183                 for (j = lo_s3[i];  j < hi_s3[i];  j++)
1184                         *s3_ptr++ = s3[j] * norm[i];
1185         }
1186 }
1187
1188 #endif          /* ORG_NUMLINES_NORM */