* renamed logf to msg
[swftools.git] / 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 */