home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / sound / misc / mpeg1_iis.lha / mpeg1_iis / decode.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-11  |  52.5 KB  |  1,525 lines
  1. /**********************************************************************
  2. Copyright (c) 1991 MPEG/audio software simulation group, All Rights Reserved
  3. decode.c
  4. **********************************************************************/
  5. /**********************************************************************
  6.  * MPEG/audio coding/decoding software, work in progress              *
  7.  *   NOT for public distribution until verified and approved by the   *
  8.  *   MPEG/audio committee.  For further information, please contact   *
  9.  *   Davis Pan, 508-493-2241, e-mail: pan@gauss.enet.dec.com          *
  10.  *                                                                    *
  11.  * VERSION 4.1                                                        *
  12.  *   changes made since last update:                                  *
  13.  *   date   programmers         comment                               *
  14.  * 2/25/91  Douglas Wong,       start of version 1.0 records          *
  15.  *          Davis Pan                                                 *
  16.  * 3/06/91  Douglas Wong        rename: setup.h to dedef.h            *
  17.  *                                      dfilter to defilter           *
  18.  *                                      dwindow to dewindow           *
  19.  *                              integrated "quantizer", "scalefactor" *
  20.  *                              combined window_samples routine into  *
  21.  *                              filter samples                        *
  22.  * 3/31/91  Bill Aspromonte     replaced read_filter by               *
  23.  *                              create_syn_filter and introduced a    *
  24.  *                              new Sub-Band Synthesis routine called *
  25.  *                              SubBandSynthesis()                    *
  26.  * 5/10/91  Vish (PRISM)        Ported to Macintosh and Unix.         *
  27.  *                              Changed "out_fifo()" so that last     *
  28.  *                              unfilled block is also written out.   *
  29.  *                              "create_syn_filter()" was modified so *
  30.  *                              that calculation precision is same as *
  31.  *                              in specification tables.              *
  32.  *                              Changed "decode_scale()" to reflect   *
  33.  *                              specifications.                       *
  34.  *                              Removed all routines used by          *
  35.  *                              "synchronize_buffer()".  This is now  *
  36.  *                              replaced by "seek_sync()".            *
  37.  *                              Incorporated Jean-Georges Fritsch's   *
  38.  *                              "bitstream.c" package.                *
  39.  *                              Deleted "reconstruct_sample()".       *
  40.  * 27jun91  dpwe (Aware)        Passed outFile and &sampFrames as     *
  41.  *                              args to out_fifo() - were global.     *
  42.  *                              Moved "alloc_*" reader to common.c.   *
  43.  *                              alloc, sblimit, stereo passed via new *
  44.  *                              'frame_params struct (were globals).  *
  45.  *                              Added JOINT STEREO decoding, lyrs I&II*
  46.  *                              Affects: decode_bitalloc,buffer_samps *
  47.  *                              Plus a few other cleanups.            *
  48.  * 6/10/91   Earle Jennings     conditional expansion added in        *
  49.  *                              II_dequantize_sample to handle range  *
  50.  *                              problems in MSDOS version             *
  51.  * 8/8/91    Jens Spille        Change for MS-C6.00                   *
  52.  *10/1/91    S.I. Sudharsanan,  Ported to IBM AIX platform.           *
  53.  *           Don H. Lee,                                              *
  54.  *           Peter W. Farrett                                         *
  55.  *10/3/91    Don H. Lee         implemented CRC-16 error protection   *
  56.  *                              newly introduced functions are        *
  57.  *                              buffer_CRC and recover_CRC_error.     *
  58.  * 2/11/92  W. Joseph Carter    Ported new code to Macintosh.  Most   *
  59.  *                              important fixes involved changing     *
  60.  *                              16-bit ints to long or unsigned in    *
  61.  *                              bit alloc routines for quant of 65535 *
  62.  *                              and passing proper function args.     *
  63.  *                              Removed "Other Joint Stereo" option   *
  64.  *                              and made bitrate be total channel     *
  65.  *                              bitrate, irrespective of the mode.    *
  66.  *                              Fixed many small bugs & reorganized.  *
  67.  * 7/27/92  Juan Pineda         Bug fix in SubBandSynthesis()         *
  68.  *--------------------------------------------------------------------*
  69.  * 6/14/92  Juan Pineda         Layer III decoding routines added.    *
  70.  *          Amit Gulati         Follows CD 3-11172 rev2.  Contains    *
  71.  *                              hacks deal with evolving available    *
  72.  *                              layerIII bitstreams.  Some (minor)    *
  73.  *                              modification of prior LI&II code.     *
  74.  * 10/25/92 Amit Gulati         Updated layerIII routines. Added code *
  75.  *                              for subblock_gain, switched block     *
  76.  *                              modes, stereo pre-processing.         *
  77.  *                              Corrected sign bits for huffman       *
  78.  *                              decoding of quadruples region and     *
  79.  *                              adjusted gain factor in III_dequant.  *
  80.  * 11/21/92 Amit Gulati         Several layerIII bugs fixed.          *
  81.  * 12/15/92 Amit Gulati         Corrected reordering (indexing)       *
  82.  *          Stan Searing        within IMDCT routine.                 *
  83.  *  8/24/93 Masahiro Iwadare    Included IS modification in Layer III.*
  84.  *                              Changed for 1 pass decoding.          *
  85.  *  9/07/93 Toshiyuki Ishino    Integrated Layer III with Ver 3.9.    *
  86.  *--------------------------------------------------------------------*
  87.  * 11/20/93 Masahiro Iwadare    Integrated Layer III with Ver 4.0.    *
  88.  *--------------------------------------------------------------------*
  89.  *  7/14/94 Juergen Koller      Bug fixes in Layer III code           *
  90.  **********************************************************************/
  91.  
  92. #include        "common.h"
  93. #include        "decoder.h"
  94. #include        "huffman.h"
  95.  
  96. /***************************************************************
  97. /*
  98. /* This module contains the core of the decoder ie all the
  99. /* computational routines. (Layer I and II only)
  100. /* Functions are common to both layer unless
  101. /* otherwise specified.
  102. /*
  103. /***************************************************************/
  104.  
  105. /*****************************************************************
  106. /*
  107. /* The following routines decode the system information
  108. /*
  109. /****************************************************************/
  110.  
  111. /************ Layer I, Layer II & Layer III ******************/
  112.  
  113. void decode_info(bs, fr_ps)
  114. Bit_stream_struc *bs;
  115. frame_params *fr_ps;
  116. {
  117.     layer *hdr = fr_ps->header;
  118.  
  119.     hdr->version = get1bit(bs);
  120.     hdr->lay = 4-getbits(bs,2);
  121.     hdr->error_protection = !get1bit(bs); /* error protect. TRUE/FALSE */
  122.     hdr->bitrate_index = getbits(bs,4);
  123.     hdr->sampling_frequency = getbits(bs,2);
  124.     hdr->padding = get1bit(bs);
  125.     hdr->extension = get1bit(bs);
  126.     hdr->mode = getbits(bs,2);
  127.     hdr->mode_ext = getbits(bs,2);
  128.     hdr->copyright = get1bit(bs);
  129.     hdr->original = get1bit(bs);
  130.     hdr->emphasis = getbits(bs,2);
  131. }
  132.  
  133. /*******************************************************************
  134. /*
  135. /* The bit allocation information is decoded. Layer I
  136. /* has 4 bit per subband whereas Layer II is Ws and bit rate
  137. /* dependent.
  138. /*
  139. /********************************************************************/
  140.  
  141. /**************************** Layer II *************/
  142.  
  143. void II_decode_bitalloc(bs, bit_alloc, fr_ps)
  144. Bit_stream_struc *bs;
  145. unsigned int bit_alloc[2][SBLIMIT];
  146. frame_params *fr_ps;
  147. {
  148.     int i,j;
  149.     int stereo = fr_ps->stereo;
  150.     int sblimit = fr_ps->sblimit;
  151.     int jsbound = fr_ps->jsbound;
  152.     al_table *alloc = fr_ps->alloc;
  153.  
  154.     for (i=0;i<jsbound;i++) for (j=0;j<stereo;j++)
  155.         bit_alloc[j][i] = (char) getbits(bs,(*alloc)[i][0].bits);
  156.  
  157.     for (i=jsbound;i<sblimit;i++) /* expand to 2 channels */
  158.         bit_alloc[0][i] = bit_alloc[1][i] =
  159.             (char) getbits(bs,(*alloc)[i][0].bits);
  160.  
  161.     for (i=sblimit;i<SBLIMIT;i++) for (j=0;j<stereo;j++)
  162.         bit_alloc[j][i] = 0;
  163. }
  164.  
  165. /**************************** Layer I *************/
  166.  
  167. void I_decode_bitalloc(bs, bit_alloc, fr_ps)
  168. Bit_stream_struc *bs;
  169. unsigned int bit_alloc[2][SBLIMIT];
  170. frame_params *fr_ps;
  171. {
  172.     int i,j;
  173.     int stereo  = fr_ps->stereo;
  174.     int sblimit = fr_ps->sblimit;
  175.     int jsbound = fr_ps->jsbound;
  176.     int b;
  177.  
  178.     for (i=0;i<jsbound;i++) for (j=0;j<stereo;j++)
  179.         bit_alloc[j][i] = getbits(bs,4);
  180.     for (i=jsbound;i<SBLIMIT;i++) {
  181.         b = getbits(bs,4);
  182.         for (j=0;j<stereo;j++)
  183.             bit_alloc[j][i] = b;
  184.     }
  185. }
  186.  
  187. /*****************************************************************
  188. /*
  189. /* The following two functions implement the layer I and II
  190. /* format of scale factor extraction. Layer I involves reading
  191. /* 6 bit per subband as scale factor. Layer II requires reading
  192. /* first the scfsi which in turn indicate the number of scale factors
  193. /* transmitted.
  194. /*    Layer I : I_decode_scale
  195. /*   Layer II : II_decode_scale
  196. /*
  197. /****************************************************************/
  198.  
  199. /************************** Layer I stuff ************************/
  200.  
  201. void I_decode_scale(bs, bit_alloc, scale_index, fr_ps)
  202. Bit_stream_struc *bs;
  203. unsigned int bit_alloc[2][SBLIMIT], scale_index[2][3][SBLIMIT];
  204. frame_params *fr_ps;
  205. {
  206.     int i,j;
  207.     int stereo = fr_ps->stereo;
  208.     int sblimit = fr_ps->sblimit;
  209.  
  210.     for (i=0;i<SBLIMIT;i++) for (j=0;j<stereo;j++)
  211.         if (!bit_alloc[j][i])
  212.             scale_index[j][0][i] = SCALE_RANGE-1;
  213.         else                    /* 6 bit per scale factor */
  214.             scale_index[j][0][i] =  getbits(bs,6);
  215.  
  216. }
  217.  
  218. /*************************** Layer II stuff ***************************/
  219.  
  220. void II_decode_scale(bs,scfsi, bit_alloc,scale_index, fr_ps)
  221. Bit_stream_struc *bs;
  222. unsigned int scfsi[2][SBLIMIT], bit_alloc[2][SBLIMIT],
  223.              scale_index[2][3][SBLIMIT];
  224. frame_params *fr_ps;
  225. {
  226.     int i,j;
  227.     int stereo = fr_ps->stereo;
  228.     int sblimit = fr_ps->sblimit;
  229.    
  230.     for (i=0;i<sblimit;i++) for (j=0;j<stereo;j++) /* 2 bit scfsi */
  231.         if (bit_alloc[j][i]) scfsi[j][i] = (char) getbits(bs,2);
  232.     for (i=sblimit;i<SBLIMIT;i++) for (j=0;j<stereo;j++)   
  233.         scfsi[j][i] = 0;
  234.  
  235.     for (i=0;i<sblimit;i++) for (j=0;j<stereo;j++) {
  236.         if (bit_alloc[j][i])   
  237.             switch (scfsi[j][i]) {
  238.                 /* all three scale factors transmitted */
  239.              case 0 : scale_index[j][0][i] = getbits(bs,6);
  240.                 scale_index[j][1][i] = getbits(bs,6);
  241.                 scale_index[j][2][i] = getbits(bs,6);
  242.                 break;
  243.                 /* scale factor 1 & 3 transmitted */
  244.              case 1 : scale_index[j][0][i] =
  245.                  scale_index[j][1][i] = getbits(bs,6);
  246.                 scale_index[j][2][i] = getbits(bs,6);
  247.                 break;
  248.                 /* scale factor 1 & 2 transmitted */
  249.              case 3 : scale_index[j][0][i] = getbits(bs,6);
  250.                 scale_index[j][1][i] =
  251.                     scale_index[j][2][i] =  getbits(bs,6);
  252.                 break;
  253.                 /* only one scale factor transmitted */
  254.              case 2 : scale_index[j][0][i] =
  255.                  scale_index[j][1][i] =
  256.                      scale_index[j][2][i] = getbits(bs,6);
  257.                 break;
  258.                 default : break;
  259.             }
  260.         else {
  261.             scale_index[j][0][i] = scale_index[j][1][i] =
  262.                 scale_index[j][2][i] = SCALE_RANGE-1;
  263.         }
  264.     }
  265.     for (i=sblimit;i<SBLIMIT;i++) for (j=0;j<stereo;j++) {
  266.         scale_index[j][0][i] = scale_index[j][1][i] =
  267.             scale_index[j][2][i] = SCALE_RANGE-1;
  268.     }
  269. }
  270.  
  271. /**************************************************************
  272. /*
  273. /*   The following two routines take care of reading the
  274. /* compressed sample from the bit stream for both layer 1 and
  275. /* layer 2. For layer 1, read the number of bits as indicated
  276. /* by the bit_alloc information. For layer 2, if grouping is
  277. /* indicated for a particular subband, then the sample size has
  278. /* to be read from the bits_group and the merged samples has
  279. /* to be decompose into the three distinct samples. Otherwise,
  280. /* it is the same for as layer one.
  281. /*
  282. /**************************************************************/
  283.  
  284. /******************************* Layer I stuff ******************/
  285.  
  286. void I_buffer_sample(bs, sample, bit_alloc, fr_ps)
  287. unsigned int FAR sample[2][3][SBLIMIT];
  288. unsigned int bit_alloc[2][SBLIMIT];
  289. Bit_stream_struc *bs;
  290. frame_params *fr_ps;
  291. {
  292.     int i,j,k;
  293.     int stereo = fr_ps->stereo;
  294.     int sblimit = fr_ps->sblimit;
  295.     int jsbound = fr_ps->jsbound;
  296.     unsigned int s;
  297.  
  298.     for (i=0;i<jsbound;i++) for (j=0;j<stereo;j++)
  299.         if ( (k = bit_alloc[j][i]) == 0)
  300.             sample[j][0][i] = 0;
  301.         else 
  302.             sample[j][0][i] = (unsigned int) getbits(bs,k+1);
  303.     for (i=jsbound;i<SBLIMIT;i++) {
  304.         if ( (k = bit_alloc[0][i]) == 0)
  305.             s = 0;
  306.         else 
  307.             s = (unsigned int)getbits(bs,k+1);
  308.         for (j=0;j<stereo;j++)
  309.             sample[j][0][i]    = s;
  310.     }
  311. }
  312.  
  313. /*************************** Layer II stuff ************************/
  314.  
  315. void II_buffer_sample(bs,sample,bit_alloc,fr_ps)
  316. unsigned int FAR sample[2][3][SBLIMIT];
  317. unsigned int bit_alloc[2][SBLIMIT];
  318. Bit_stream_struc *bs;
  319. frame_params *fr_ps;
  320. {
  321.     int i,j,k,m;
  322.     int stereo = fr_ps->stereo;
  323.     int sblimit = fr_ps->sblimit;
  324.     int jsbound = fr_ps->jsbound;
  325.     al_table *alloc = fr_ps->alloc;
  326.  
  327.     for (i=0;i<sblimit;i++) for (j=0;j<((i<jsbound)?stereo:1);j++) {
  328.         if (bit_alloc[j][i]) {
  329.             /* check for grouping in subband */
  330.             if ((*alloc)[i][bit_alloc[j][i]].group==3)
  331.                 for (m=0;m<3;m++) {
  332.                     k = (*alloc)[i][bit_alloc[j][i]].bits;
  333.                     sample[j][m][i] = (unsigned int) getbits(bs,k);
  334.                 }         
  335.             else {              /* bit_alloc = 3, 5, 9 */
  336.                 unsigned int nlevels, c=0;
  337.  
  338.                 nlevels = (*alloc)[i][bit_alloc[j][i]].steps;
  339.                 k=(*alloc)[i][bit_alloc[j][i]].bits;
  340.                 c = (unsigned int) getbits(bs, k);
  341.                 for (k=0;k<3;k++) {
  342.                     sample[j][k][i] = c % nlevels;
  343.                     c /= nlevels;
  344.                 }
  345.             }
  346.         }
  347.         else {                  /* for no sample transmitted */
  348.             for (k=0;k<3;k++) sample[j][k][i] = 0;
  349.         }
  350.         if(stereo == 2 && i>= jsbound) /* joint stereo : copy L to R */
  351.             for (k=0;k<3;k++) sample[1][k][i] = sample[0][k][i];
  352.     }
  353.     for (i=sblimit;i<SBLIMIT;i++) for (j=0;j<stereo;j++) for (k=0;k<3;k++)
  354.         sample[j][k][i] = 0;
  355. }      
  356.  
  357. /**************************************************************
  358. /*
  359. /*   Restore the compressed sample to a factional number.
  360. /*   first complement the MSB of the sample
  361. /*    for layer I :
  362. /*    Use s = (s' + 2^(-nb+1) ) * 2^nb / (2^nb-1)
  363. /*   for Layer II :
  364. /*   Use the formula s = s' * c + d
  365. /*
  366. /**************************************************************/
  367.  
  368. static double c[17] = { 1.33333333333, 1.60000000000, 1.14285714286,
  369.                         1.77777777777, 1.06666666666, 1.03225806452,
  370.                         1.01587301587, 1.00787401575, 1.00392156863,
  371.                         1.00195694716, 1.00097751711, 1.00048851979,
  372.                         1.00024420024, 1.00012208522, 1.00006103888,
  373.                         1.00003051851, 1.00001525902 };
  374.  
  375. static double d[17] = { 0.500000000, 0.500000000, 0.250000000, 0.500000000,
  376.                         0.125000000, 0.062500000, 0.031250000, 0.015625000,
  377.                         0.007812500, 0.003906250, 0.001953125, 0.0009765625,
  378.                         0.00048828125, 0.00024414063, 0.00012207031,
  379.                         0.00006103516, 0.00003051758 };
  380.  
  381. /************************** Layer II stuff ************************/
  382.  
  383. void II_dequantize_sample(sample, bit_alloc, fraction, fr_ps)
  384. unsigned int FAR sample[2][3][SBLIMIT];
  385. unsigned int bit_alloc[2][SBLIMIT];
  386. double FAR fraction[2][3][SBLIMIT];
  387. frame_params *fr_ps;
  388. {
  389.     int i, j, k, x;
  390.     int stereo = fr_ps->stereo;
  391.     int sblimit = fr_ps->sblimit;
  392.     al_table *alloc = fr_ps->alloc;
  393.  
  394.     for (i=0;i<sblimit;i++)  for (j=0;j<3;j++) for (k=0;k<stereo;k++)
  395.         if (bit_alloc[k][i]) {
  396.  
  397.             /* locate MSB in the sample */
  398.             x = 0;
  399. #ifndef MS_DOS
  400.             while ((1L<<x) < (*alloc)[i][bit_alloc[k][i]].steps) x++;
  401. #else
  402.             /* microsoft C thinks an int is a short */
  403.             while (( (unsigned long) (1L<<(long)x) <
  404.                     (unsigned long)( (*alloc)[i][bit_alloc[k][i]].steps)
  405.                     ) && ( x < 16) ) x++;
  406. #endif
  407.  
  408.             /* MSB inversion */
  409.             if (((sample[k][j][i] >> x-1) & 1) == 1)
  410.                 fraction[k][j][i] = 0.0;
  411.             else  fraction[k][j][i] = -1.0;
  412.  
  413.             /* Form a 2's complement sample */
  414.             fraction[k][j][i] += (double) (sample[k][j][i] & ((1<<x-1)-1)) /
  415.                 (double) (1L<<x-1);
  416.  
  417.             /* Dequantize the sample */
  418.             fraction[k][j][i] += d[(*alloc)[i][bit_alloc[k][i]].quant];
  419.             fraction[k][j][i] *= c[(*alloc)[i][bit_alloc[k][i]].quant];
  420.         }
  421.         else fraction[k][j][i] = 0.0;   
  422.    
  423.     for (i=sblimit;i<SBLIMIT;i++) for (j=0;j<3;j++) for(k=0;k<stereo;k++)
  424.         fraction[k][j][i] = 0.0;
  425. }
  426.  
  427. /***************************** Layer I stuff ***********************/
  428.  
  429. void I_dequantize_sample(sample, fraction, bit_alloc, fr_ps)
  430. unsigned int FAR sample[2][3][SBLIMIT];
  431. unsigned int bit_alloc[2][SBLIMIT];
  432. double FAR fraction[2][3][SBLIMIT];
  433. frame_params *fr_ps;
  434. {
  435.     int i, nb, k;
  436.     int stereo = fr_ps->stereo;
  437.     int sblimit = fr_ps->sblimit;
  438.  
  439.     for (i=0;i<SBLIMIT;i++)
  440.         for (k=0;k<stereo;k++)
  441.             if (bit_alloc[k][i]) {
  442.                 nb = bit_alloc[k][i] + 1;
  443.                 if (((sample[k][0][i] >> nb-1) & 1) == 1) fraction[k][0][i] = 0.0;
  444.                 else fraction[k][0][i] = -1.0;
  445.                 fraction[k][0][i] += (double) (sample[k][0][i] & ((1<<nb-1)-1)) /
  446.                     (double) (1L<<nb-1);
  447.  
  448.                 fraction[k][0][i] =
  449.                     (double) (fraction[k][0][i]+1.0/(double)(1L<<nb-1)) *
  450.                         (double) (1L<<nb) / (double) ((1L<<nb)-1);
  451.             }
  452.             else fraction[k][0][i] = 0.0;
  453. }
  454.  
  455. /************************************************************
  456. /*
  457. /*   Restore the original value of the sample ie multiply
  458. /*    the fraction value by its scalefactor.
  459. /*
  460. /************************************************************/
  461.  
  462. /************************* Layer II Stuff **********************/
  463.  
  464. void II_denormalize_sample(fraction, scale_index,fr_ps,x)
  465. double FAR fraction[2][3][SBLIMIT];
  466. unsigned int scale_index[2][3][SBLIMIT];
  467. frame_params *fr_ps;
  468. int x;
  469. {
  470.     int i,j,k;
  471.     int stereo = fr_ps->stereo;
  472.     int sblimit = fr_ps->sblimit;
  473.  
  474.     for (i=0;i<sblimit;i++) for (j=0;j<stereo;j++) {
  475.         fraction[j][0][i] *= multiple[scale_index[j][x][i]];
  476.         fraction[j][1][i] *= multiple[scale_index[j][x][i]];
  477.         fraction[j][2][i] *= multiple[scale_index[j][x][i]];
  478.     }
  479. }
  480.  
  481. /**************************** Layer I stuff ******************************/
  482.  
  483. void I_denormalize_sample(fraction,scale_index,fr_ps)
  484. double FAR fraction[2][3][SBLIMIT];
  485. unsigned int scale_index[2][3][SBLIMIT];
  486. frame_params *fr_ps;
  487. {
  488.     int i,j,k;
  489.     int stereo = fr_ps->stereo;
  490.     int sblimit = fr_ps->sblimit;
  491.  
  492.     for (i=0;i<SBLIMIT;i++) for (j=0;j<stereo;j++)
  493.         fraction[j][0][i] *= multiple[scale_index[j][0][i]];
  494. }
  495.  
  496. /*****************************************************************
  497. /*
  498. /* The following are the subband synthesis routines. They apply
  499. /* to both layer I and layer II stereo or mono. The user has to
  500. /* decide what parameters are to be passed to the routines.
  501. /*
  502. /***************************************************************/
  503.  
  504. /*************************************************************
  505. /*
  506. /*   Pass the subband sample through the synthesis window
  507. /*
  508. /**************************************************************/
  509.  
  510. /* create in synthesis filter */
  511.  
  512. void create_syn_filter(filter)
  513. double FAR filter[64][SBLIMIT];
  514. {
  515.     register int i,k;
  516.  
  517.     for (i=0; i<64; i++)
  518.         for (k=0; k<32; k++) {
  519.             if ((filter[i][k] = 1e9*cos((double)((PI64*i+PI4)*(2*k+1)))) >= 0)
  520.                 modf(filter[i][k]+0.5, &filter[i][k]);
  521.             else
  522.                 modf(filter[i][k]-0.5, &filter[i][k]);
  523.             filter[i][k] *= 1e-9;
  524.         }
  525. }
  526.  
  527. /***************************************************************
  528. /*
  529. /*   Window the restored sample
  530. /*
  531. /***************************************************************/
  532.  
  533. /* read in synthesis window */
  534.  
  535. void read_syn_window(window)
  536. double FAR window[HAN_SIZE];
  537. {
  538.     int i,j[4];
  539.     FILE *fp;
  540.     double f[4];
  541.     char t[150];
  542.  
  543.     if (!(fp = OpenTableFile("dewindow") )) {
  544.         printf("Please check synthesis window table 'dewindow'\n");
  545.         exit(1);
  546.     }
  547.     for (i=0;i<512;i+=4) {
  548.         fgets(t, 150, fp);
  549.         sscanf(t,"D[%d] = %lf D[%d] = %lf D[%d] = %lf D[%d] = %lf\n",
  550.                j, f,j+1,f+1,j+2,f+2,j+3,f+3);
  551.         if (i==j[0]) {
  552.             window[i] = f[0];
  553.             window[i+1] = f[1];
  554.             window[i+2] = f[2];
  555.             window[i+3] = f[3];
  556.         }
  557.         else {
  558.             printf("Check index in synthesis window table\n");
  559.             exit(1);
  560.         }
  561.         fgets(t,150,fp);
  562.     }
  563.     fclose(fp);
  564. }
  565.  
  566. int SubBandSynthesis (bandPtr, channel, samples)
  567. double *bandPtr;
  568. int channel;
  569. short *samples;
  570. {
  571.     register int i,j,k;
  572.     register double *bufOffsetPtr, sum;
  573.     static int init = 1;
  574.     typedef double NN[64][32];
  575.     static NN FAR *filter;
  576.     typedef double BB[2][2*HAN_SIZE];
  577.     static BB FAR *buf;
  578.     static int bufOffset[2] = {64,64};
  579.     static double FAR *window;
  580.     int clip = 0;               /* count & return how many samples clipped */
  581.  
  582.     if (init) {
  583.         buf = (BB FAR *) mem_alloc(sizeof(BB),"BB");
  584.         filter = (NN FAR *) mem_alloc(sizeof(NN), "NN");
  585.         create_syn_filter(*filter);
  586.         window = (double FAR *) mem_alloc(sizeof(double) * HAN_SIZE, "WIN");
  587.         read_syn_window(window);
  588.         init = 0;
  589.     }
  590. /*    if (channel == 0) */
  591.     bufOffset[channel] = (bufOffset[channel] - 64) & 0x3ff;
  592.     bufOffsetPtr = &((*buf)[channel][bufOffset[channel]]);
  593.  
  594.     for (i=0; i<64; i++) {
  595.         sum = 0;
  596.         for (k=0; k<32; k++)
  597.             sum += bandPtr[k] * (*filter)[i][k];
  598.         bufOffsetPtr[i] = sum;
  599.     }
  600.     /*  S(i,j) = D(j+32i) * U(j+32i+((i+1)>>1)*64)  */
  601.     /*  samples(i,j) = MWindow(j+32i) * bufPtr(j+32i+((i+1)>>1)*64)  */
  602.     for (j=0; j<32; j++) {
  603.         sum = 0;
  604.         for (i=0; i<16; i++) {
  605.             k = j + (i<<5);
  606.             sum += window[k] * (*buf) [channel] [( (k + ( ((i+1)>>1) <<6) ) +
  607.                                                   bufOffset[channel]) & 0x3ff];
  608.         }
  609.  
  610. /*   {long foo = (sum > 0) ? sum * SCALE + 0.5 : sum * SCALE - 0.5; */
  611.      {long foo = sum * SCALE;
  612.      if (foo >= (long) SCALE)      {samples[j] = SCALE-1; ++clip;}
  613.      else if (foo < (long) -SCALE) {samples[j] = -SCALE;  ++clip;}
  614.      else                           samples[j] = foo;
  615.  }
  616.     }
  617.     return(clip);
  618. }
  619.  
  620. void out_fifo(pcm_sample, num, fr_ps, done, outFile, psampFrames)
  621. short FAR pcm_sample[2][SSLIMIT][SBLIMIT];
  622. int num;
  623. frame_params *fr_ps;
  624. int done;
  625. FILE *outFile;
  626. unsigned long *psampFrames;
  627. {
  628.     int i,j,l;
  629.     int stereo = fr_ps->stereo;
  630.     int sblimit = fr_ps->sblimit;
  631.     static short int outsamp[1600];
  632.     static long k = 0;
  633.  
  634.     if (!done)
  635.         for (i=0;i<num;i++) for (j=0;j<SBLIMIT;j++) {
  636.             (*psampFrames)++;
  637.             for (l=0;l<stereo;l++) {
  638.                 if (!(k%1600) && k) {
  639.                     fwrite(outsamp,2,1600,outFile);
  640.                     k = 0;
  641.                 }
  642.                 outsamp[k++] = pcm_sample[l][i][j];
  643.             }
  644.         }
  645.     else {
  646.         fwrite(outsamp,2,(int)k,outFile);
  647.         k = 0;
  648.     }
  649. }
  650.  
  651. void  buffer_CRC(bs, old_crc)
  652. Bit_stream_struc  *bs;
  653. unsigned int  *old_crc;
  654. {
  655.     *old_crc = getbits(bs, 16);
  656. }
  657.  
  658. void  recover_CRC_error(pcm_sample, error_count, fr_ps, outFile, psampFrames)
  659. short FAR pcm_sample[2][SSLIMIT][SBLIMIT];
  660. int error_count;
  661. frame_params *fr_ps;
  662. FILE *outFile;
  663. unsigned long *psampFrames;
  664. {
  665.     int  stereo = fr_ps->stereo;
  666.     int  num, done, i;
  667.     int  samplesPerFrame, samplesPerSlot;
  668.     layer  *hdr = fr_ps->header;
  669.     long  offset;
  670.     short  *temp;
  671.  
  672.     num = 3;
  673.     if (hdr->lay == 1) num = 1;
  674.  
  675.     samplesPerSlot = SBLIMIT * num * stereo;
  676.     samplesPerFrame = samplesPerSlot * 32;
  677.  
  678.     if (error_count == 1) {     /* replicate previous error_free frame */
  679.         done = 1;
  680.         /* flush out fifo */
  681.         out_fifo(pcm_sample, num, fr_ps, done, outFile, psampFrames);
  682.         /* go back to the beginning of the previous frame */
  683.         offset = sizeof(short int) * samplesPerFrame;
  684.         fseek(outFile, -offset, SEEK_CUR);
  685.         done = 0;
  686.         for (i = 0; i < SCALE_BLOCK; i++) {
  687.             fread(pcm_sample, 2, samplesPerSlot, outFile);
  688.             out_fifo(pcm_sample, num, fr_ps, done, outFile, psampFrames);
  689.         }
  690.     }
  691.     else{                       /* mute the frame */
  692.         temp = (short*) pcm_sample;
  693.         done = 0;
  694.         for (i = 0; i < 2*3*SBLIMIT; i++)
  695.             *temp++ = MUTE;     /* MUTE value is in decoder.h */
  696.         for (i = 0; i < SCALE_BLOCK; i++)
  697.             out_fifo(pcm_sample, num, fr_ps, done, outFile, psampFrames);
  698.     }
  699. }
  700.  
  701. /************************* Layer III routines **********************/
  702.  
  703. void III_get_side_info(bs, si, fr_ps)
  704. Bit_stream_struc *bs;
  705. III_side_info_t *si;
  706. frame_params *fr_ps;
  707. {
  708.    int ch, gr, i;
  709.    int stereo = fr_ps->stereo;
  710.    
  711.    si->main_data_begin = getbits(bs, 9);
  712.    if (stereo == 1)
  713.       si->private_bits = getbits(bs,5);
  714.       else si->private_bits = getbits(bs,3);
  715.  
  716.    for (ch=0; ch<stereo; ch++)
  717.       for (i=0; i<4; i++)
  718.       si->ch[ch].scfsi[i] = get1bit(bs);
  719.  
  720.    for (gr=0; gr<2; gr++) {
  721.       for (ch=0; ch<stereo; ch++) {
  722.          si->ch[ch].gr[gr].part2_3_length = getbits(bs, 12);
  723.          si->ch[ch].gr[gr].big_values = getbits(bs, 9);
  724.          si->ch[ch].gr[gr].global_gain = getbits(bs, 8);
  725.          si->ch[ch].gr[gr].scalefac_compress = getbits(bs, 4);
  726.          si->ch[ch].gr[gr].window_switching_flag = get1bit(bs);
  727.          if (si->ch[ch].gr[gr].window_switching_flag) {
  728.             si->ch[ch].gr[gr].block_type = getbits(bs, 2);
  729.             si->ch[ch].gr[gr].mixed_block_flag = get1bit(bs);
  730.             for (i=0; i<2; i++)
  731.                si->ch[ch].gr[gr].table_select[i] = getbits(bs, 5);
  732.             for (i=0; i<3; i++)
  733.                si->ch[ch].gr[gr].subblock_gain[i] = getbits(bs, 3);
  734.                
  735.             /* Set region_count parameters since they are implicit in this case. */
  736.             
  737.             if (si->ch[ch].gr[gr].block_type == 0) {
  738.                printf("Side info bad: block_type == 0 in split block.\n");
  739.                exit(0);
  740.                }
  741.             else if (si->ch[ch].gr[gr].block_type == 2
  742.                      && si->ch[ch].gr[gr].mixed_block_flag == 0)
  743.                si->ch[ch].gr[gr].region0_count = 8; /* MI 9; */
  744.             else si->ch[ch].gr[gr].region0_count = 7; /* MI 8; */
  745.             si->ch[ch].gr[gr].region1_count = 20 -
  746.                       si->ch[ch].gr[gr].region0_count;
  747.             }
  748.          else {
  749.             for (i=0; i<3; i++)
  750.                si->ch[ch].gr[gr].table_select[i] = getbits(bs, 5);
  751.             si->ch[ch].gr[gr].region0_count = getbits(bs, 4);
  752.             si->ch[ch].gr[gr].region1_count = getbits(bs, 3);
  753.             si->ch[ch].gr[gr].block_type = 0;
  754.             }      
  755.          si->ch[ch].gr[gr].preflag = get1bit(bs);
  756.          si->ch[ch].gr[gr].scalefac_scale = get1bit(bs);
  757.          si->ch[ch].gr[gr].count1table_select = get1bit(bs);
  758.          }
  759.       }
  760. }
  761.  
  762. void III_put_side_info(bs, si, fr_ps)
  763. frame_params *fr_ps;
  764. Bit_stream_struc *bs;
  765. III_side_info_t *si;
  766. {
  767.    int ch, gr, i;
  768.    int stereo = fr_ps->stereo;
  769.  
  770.    putbits(bs, si->main_data_begin,9);
  771.    if (stereo == 1)
  772.       putbits(bs, si->private_bits, 5);
  773.       else putbits(bs, si->private_bits, 3);
  774.  
  775.    for (ch=0; ch<stereo; ch++)
  776.       for (i=0; i<4; i++)
  777.          put1bit(bs, si->ch[ch].scfsi[i]);
  778.  
  779.    for (gr=0; gr<2; gr++) {
  780.       for (ch=0; ch<stereo; ch++) {
  781.          putbits(bs, si->ch[ch].gr[gr].part2_3_length, 12);
  782.          putbits(bs, si->ch[ch].gr[gr].big_values, 9);
  783.          putbits(bs, si->ch[ch].gr[gr].global_gain, 8);
  784.          putbits(bs, si->ch[ch].gr[gr].scalefac_compress, 4);
  785.          put1bit(bs, si->ch[ch].gr[gr].window_switching_flag);
  786.          if (si->ch[ch].gr[gr].window_switching_flag) {
  787.             putbits(bs, si->ch[ch].gr[gr].block_type, 2);
  788.             put1bit(bs, si->ch[ch].gr[gr].mixed_block_flag);
  789.             for (i=0; i<2; i++)
  790.                putbits(bs, si->ch[ch].gr[gr].table_select[i], 5);
  791.             for (i=0; i<3; i++)
  792.                putbits(bs, si->ch[ch].gr[gr].subblock_gain[i], 3);
  793.             }
  794.          else {
  795.             for (i=0; i<3; i++)
  796.             putbits(bs, si->ch[ch].gr[gr].table_select[i], 5);
  797.             putbits(bs, si->ch[ch].gr[gr].region0_count, 4);
  798.             putbits(bs, si->ch[ch].gr[gr].region1_count, 3);
  799.             }      
  800.  
  801.          put1bit(bs, si->ch[ch].gr[gr].preflag);
  802.          put1bit(bs, si->ch[ch].gr[gr].scalefac_scale);
  803.          put1bit(bs, si->ch[ch].gr[gr].count1table_select);
  804.          }
  805.       }
  806. }
  807.  
  808. struct {
  809.    int l[5];
  810.    int s[3];} sfbtable = {{0, 6, 11, 16, 21},
  811.                           {0, 6, 12}};
  812.                          
  813. int slen[2][16] = {{0, 0, 0, 0, 3, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4},
  814.                    {0, 1, 2, 3, 0, 1, 2, 3, 1, 2, 3, 1, 2, 3, 2, 3}};
  815. struct  {
  816.    int l[23];
  817.    int s[14];} sfBandIndex[3] =   
  818.    {{{0,4,8,12,16,20,24,30,36,44,52,62,74,90,110,134,162,196,238,288,342,418,576},
  819.      {0,4,8,12,16,22,30,40,52,66,84,106,136,192}},
  820.     {{0,4,8,12,16,20,24,30,36,42,50,60,72,88,106,128,156,190,230,276,330,384,576},
  821.      {0,4,8,12,16,22,28,38,50,64,80,100,126,192}},
  822.     {{0,4,8,12,16,20,24,30,36,44,54,66,82,102,126,156,194,240,296,364,448,550,576},
  823.      {0,4,8,12,16,22,30,42,58,78,104,138,180,192}}};
  824.  
  825.  
  826. void III_get_scale_factors(scalefac, si, gr, ch, fr_ps)
  827. III_scalefac_t *scalefac;
  828. III_side_info_t *si;
  829. int gr, ch;
  830. frame_params *fr_ps;
  831. {
  832. int sfb, i, window;
  833. struct gr_info_s *gr_info = &(si->ch[ch].gr[gr]);
  834.  
  835.     if (gr_info->window_switching_flag && (gr_info->block_type == 2)) { 
  836.       if (gr_info->mixed_block_flag) { /* MIXED */ /* NEW - ag 11/25 */
  837.          for (sfb = 0; sfb < 8; sfb++)
  838.             (*scalefac)[ch].l[sfb] = hgetbits( 
  839.                  slen[0][gr_info->scalefac_compress]);
  840.          for (sfb = 3; sfb < 6; sfb++)
  841.             for (window=0; window<3; window++)
  842.                (*scalefac)[ch].s[window][sfb] = hgetbits(
  843.                  slen[0][gr_info->scalefac_compress]);
  844.          for (sfb = 6; sfb < 12; sfb++)
  845.             for (window=0; window<3; window++)
  846.                (*scalefac)[ch].s[window][sfb] = hgetbits(
  847.                  slen[1][gr_info->scalefac_compress]);
  848.          for (sfb=12,window=0; window<3; window++)
  849.             (*scalefac)[ch].s[window][sfb] = 0;
  850.       }
  851.       else {  /* SHORT*/
  852.          for (i=0; i<2; i++) 
  853.             for (sfb = sfbtable.s[i]; sfb < sfbtable.s[i+1]; sfb++)
  854.                for (window=0; window<3; window++)
  855.                   (*scalefac)[ch].s[window][sfb] = hgetbits( 
  856.                     slen[i][gr_info->scalefac_compress]);
  857.          for (sfb=12,window=0; window<3; window++)
  858.             (*scalefac)[ch].s[window][sfb] = 0;
  859.       }
  860.     }          
  861.     else {   /* LONG types 0,1,3 */
  862.         for (i=0; i<4; i++) {
  863.            if ((si->ch[ch].scfsi[i] == 0) || (gr == 0))
  864.               for (sfb = sfbtable.l[i]; sfb < sfbtable.l[i+1]; sfb++)
  865.                   (*scalefac)[ch].l[sfb] = hgetbits(
  866.                  slen[(i<2)?0:1][gr_info->scalefac_compress]);
  867.         }
  868.         (*scalefac)[ch].l[22] = 0; 
  869.     }
  870. }
  871.  
  872.  
  873. /* Already declared in huffman.c
  874. struct huffcodetab ht[HTN];
  875. */
  876. int huffman_initialized = FALSE;
  877.  
  878. void initialize_huffman() {
  879.    FILE *fi;
  880.   
  881.    if (huffman_initialized) return;
  882.    if (!(fi = OpenTableFile("huffdec") )) {
  883.       printf("Please check huffman table 'huffdec'\n");
  884.       exit(1);
  885.    }
  886.  
  887.    if (fi==NULL) {
  888.  
  889.       fprintf(stderr,"decoder table open error\n");
  890.  
  891.       exit(3);
  892.  
  893.       }
  894.  
  895.    if (read_decoder_table(fi) != HTN) {
  896.       fprintf(stderr,"decoder table read error\n");
  897.       exit(4);
  898.       }
  899. huffman_initialized = TRUE;
  900. }
  901.  
  902. III_hufman_decode(is, si, ch, gr, part2_start, fr_ps)
  903. long int is[SBLIMIT][SSLIMIT];
  904. III_side_info_t *si;
  905. int gr, ch, part2_start;
  906. frame_params *fr_ps;
  907. {
  908.    int i, x, y;
  909.    int v, w;
  910.    struct huffcodetab *h;
  911.    int region1Start;
  912.    int region2Start;
  913.    int bt = (*si).ch[ch].gr[gr].window_switching_flag && ((*si).ch[ch].gr[gr].block_type == 2);
  914.  
  915.    initialize_huffman();
  916.  
  917.    /* Find region boundary for short block case. */
  918.    
  919.    if ( ((*si).ch[ch].gr[gr].window_switching_flag) && 
  920.         ((*si).ch[ch].gr[gr].block_type == 2) ) { 
  921.    
  922.       /* Region2. */
  923.  
  924.       region1Start = 36;  /* sfb[9/3]*3=36 */
  925.       region2Start = 576; /* No Region2 for short block case. */
  926.    }
  927.  
  928.  
  929.    else {          /* Find region boundary for long block case. */
  930.  
  931.       region1Start = sfBandIndex[fr_ps->header->sampling_frequency]
  932.                            .l[(*si).ch[ch].gr[gr].region0_count + 1]; /* MI */
  933.       region2Start = sfBandIndex[fr_ps->header->sampling_frequency]
  934.                               .l[(*si).ch[ch].gr[gr].region0_count +
  935.                               (*si).ch[ch].gr[gr].region1_count + 2]; /* MI */
  936.       }
  937.  
  938.  
  939.    /* Read bigvalues area. */
  940.    for (i=0; i<(*si).ch[ch].gr[gr].big_values*2; i+=2) {
  941.       if      (i<region1Start) h = &ht[(*si).ch[ch].gr[gr].table_select[0]];
  942.       else if (i<region2Start) h = &ht[(*si).ch[ch].gr[gr].table_select[1]];
  943.            else                h = &ht[(*si).ch[ch].gr[gr].table_select[2]];
  944.       huffman_decoder(h, &x, &y, &v, &w);
  945.       is[i/SSLIMIT][i%SSLIMIT] = x;
  946.       is[(i+1)/SSLIMIT][(i+1)%SSLIMIT] = y;
  947.       }
  948.  
  949.    /* Read count1 area. */
  950.    h = &ht[(*si).ch[ch].gr[gr].count1table_select+32];
  951.    while ((hsstell() < part2_start + (*si).ch[ch].gr[gr].part2_3_length ) &&
  952.      ( i < SSLIMIT*SBLIMIT )) {
  953.       huffman_decoder(h, &x, &y, &v, &w);
  954.       is[i/SSLIMIT][i%SSLIMIT] = v;
  955.       is[(i+1)/SSLIMIT][(i+1)%SSLIMIT] = w;
  956.       is[(i+2)/SSLIMIT][(i+2)%SSLIMIT] = x;
  957.       is[(i+3)/SSLIMIT][(i+3)%SSLIMIT] = y;
  958.       i += 4;
  959.       }
  960.  
  961.    if (hsstell() > part2_start + (*si).ch[ch].gr[gr].part2_3_length)
  962.    {  i -=4;
  963.       rewindNbits(hsstell()-part2_start - (*si).ch[ch].gr[gr].part2_3_length);
  964.    }
  965.  
  966.    /* Dismiss stuffing Bits */
  967.    if ( hsstell() < part2_start + (*si).ch[ch].gr[gr].part2_3_length )
  968.       hgetbits( part2_start + (*si).ch[ch].gr[gr].part2_3_length - hsstell());
  969.  
  970.    /* Zero out rest. */
  971.    for (; i<SSLIMIT*SBLIMIT; i++)
  972.       is[i/SSLIMIT][i%SSLIMIT] = 0;
  973. }
  974.  
  975.  
  976. int pretab[22] = {0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,2,2,3,3,3,2,0};
  977.  
  978. void III_dequantize_sample(is,xr,scalefac,gr_info, ch,fr_ps)
  979. long int is[SBLIMIT][SSLIMIT];
  980. double xr[SBLIMIT][SSLIMIT];
  981. struct gr_info_s *gr_info;
  982. III_scalefac_t *scalefac;
  983. frame_params *fr_ps;
  984. int ch;
  985. {
  986.    int ss,sb,cb=0,sfreq=fr_ps->header->sampling_frequency;
  987.    int stereo = fr_ps->stereo;
  988.    int next_cb_boundary, cb_begin, cb_width, sign;
  989.  
  990.    /* choose correct scalefactor band per block type, initalize boundary */
  991.  
  992.    if (gr_info->window_switching_flag && (gr_info->block_type == 2) )
  993.       if (gr_info->mixed_block_flag) 
  994.          next_cb_boundary=sfBandIndex[sfreq].l[1];  /* LONG blocks: 0,1,3 */
  995.       else {
  996.          next_cb_boundary=sfBandIndex[sfreq].s[1]*3; /* pure SHORT block */
  997.     cb_width = sfBandIndex[sfreq].s[1];
  998.     cb_begin = 0;
  999.       }  
  1000.    else 
  1001.       next_cb_boundary=sfBandIndex[sfreq].l[1];  /* LONG blocks: 0,1,3 */
  1002.  
  1003.    /* apply formula per block type */
  1004.  
  1005.    for (sb=0 ; sb < SBLIMIT ; sb++)
  1006.       for (ss=0 ; ss < SSLIMIT ; ss++) {
  1007.  
  1008.          if ( (sb*18)+ss == next_cb_boundary)  { /* Adjust critical band boundary */
  1009.             if (gr_info->window_switching_flag && (gr_info->block_type == 2)) {
  1010.                if (gr_info->mixed_block_flag)  {
  1011.                   if (((sb*18)+ss) == sfBandIndex[sfreq].l[8])  {
  1012.                      next_cb_boundary=sfBandIndex[sfreq].s[4]*3; 
  1013.                      cb = 3;
  1014.                      cb_width = sfBandIndex[sfreq].s[cb+1] - 
  1015.                                 sfBandIndex[sfreq].s[cb];
  1016.                      cb_begin = sfBandIndex[sfreq].s[cb]*3;      
  1017.                   }
  1018.                   else if (((sb*18)+ss) < sfBandIndex[sfreq].l[8]) 
  1019.                      next_cb_boundary = sfBandIndex[sfreq].l[(++cb)+1];
  1020.                   else {
  1021.                      next_cb_boundary = sfBandIndex[sfreq].s[(++cb)+1]*3;
  1022.                      cb_width = sfBandIndex[sfreq].s[cb+1] - 
  1023.                                     sfBandIndex[sfreq].s[cb];
  1024.                      cb_begin = sfBandIndex[sfreq].s[cb]*3;      
  1025.                   }   
  1026.                }
  1027.                else  {
  1028.                   next_cb_boundary = sfBandIndex[sfreq].s[(++cb)+1]*3;
  1029.                   cb_width = sfBandIndex[sfreq].s[cb+1] - 
  1030.                                sfBandIndex[sfreq].s[cb];
  1031.                 cb_begin = sfBandIndex[sfreq].s[cb]*3;      
  1032.                } 
  1033.             }
  1034.             else /* long blocks */
  1035.                next_cb_boundary = sfBandIndex[sfreq].l[(++cb)+1];
  1036.          }
  1037.  
  1038.          /* Compute overall (global) scaling. */
  1039.  
  1040.          xr[sb][ss] = pow( 2.0 , (0.25 * (gr_info->global_gain - 210.0)));
  1041.  
  1042.          /* Do long/short dependent scaling operations. */
  1043.         
  1044.          if (gr_info->window_switching_flag && (
  1045.             ((gr_info->block_type == 2) && (gr_info->mixed_block_flag == 0)) ||
  1046.             ((gr_info->block_type == 2) && gr_info->mixed_block_flag && (sb >= 2)) )) {
  1047.  
  1048.             xr[sb][ss] *= pow(2.0, 0.25 * -8.0 * 
  1049.                     gr_info->subblock_gain[(((sb*18)+ss) - cb_begin)/cb_width]);
  1050.             xr[sb][ss] *= pow(2.0, 0.25 * -2.0 * (1.0+gr_info->scalefac_scale)
  1051.               * (*scalefac)[ch].s[(((sb*18)+ss) - cb_begin)/cb_width][cb]);
  1052.          }
  1053.          else {   /* LONG block types 0,1,3 & 1st 2 subbands of switched blocks */
  1054.             xr[sb][ss] *= pow(2.0, -0.5 * (1.0+gr_info->scalefac_scale)
  1055.                                         * ((*scalefac)[ch].l[cb]
  1056.                                         + gr_info->preflag * pretab[cb]));
  1057.          }
  1058.  
  1059.          /* Scale quantized value. */
  1060.         
  1061.          sign = (is[sb][ss]<0) ? 1 : 0; 
  1062.          xr[sb][ss] *= pow( (double) abs(is[sb][ss]), ((double)4.0/3.0) );
  1063.          if (sign) xr[sb][ss] = -xr[sb][ss];
  1064.       }
  1065. }
  1066.  
  1067. III_reorder (xr, ro, gr_info, fr_ps) 
  1068. double xr[SBLIMIT][SSLIMIT]; 
  1069. double ro[SBLIMIT][SSLIMIT]; 
  1070. struct gr_info_s *gr_info;
  1071. frame_params *fr_ps;
  1072. {
  1073.    int sfreq=fr_ps->header->sampling_frequency;
  1074.    int sfb, sfb_start, sfb_lines;
  1075.    int sb, ss, window, freq, src_line, des_line;
  1076.  
  1077.    for(sb=0;sb<SBLIMIT;sb++)
  1078.       for(ss=0;ss<SSLIMIT;ss++) 
  1079.          ro[sb][ss] = 0;
  1080.  
  1081.    if (gr_info->window_switching_flag && (gr_info->block_type == 2)) {
  1082.       if (gr_info->mixed_block_flag) {
  1083.          /* NO REORDER FOR LOW 2 SUBBANDS */
  1084.          for (sb=0 ; sb < 2 ; sb++)
  1085.             for (ss=0 ; ss < SSLIMIT ; ss++) {
  1086.                ro[sb][ss] = xr[sb][ss];
  1087.             }
  1088.          /* REORDERING FOR REST SWITCHED SHORT */
  1089.          for(sfb=3,sfb_start=sfBandIndex[sfreq].s[3],
  1090.             sfb_lines=sfBandIndex[sfreq].s[4] - sfb_start; 
  1091.             sfb < 13; sfb++,sfb_start=sfBandIndex[sfreq].s[sfb],
  1092.             (sfb_lines=sfBandIndex[sfreq].s[sfb+1] - sfb_start))
  1093.                for(window=0; window<3; window++)
  1094.                   for(freq=0;freq<sfb_lines;freq++) {
  1095.                      src_line = sfb_start*3 + window*sfb_lines + freq; 
  1096.                      des_line = (sfb_start*3) + window + (freq*3);
  1097.                      ro[des_line/SSLIMIT][des_line%SSLIMIT] = 
  1098.                                     xr[src_line/SSLIMIT][src_line%SSLIMIT];
  1099.                }
  1100.       } 
  1101.       else {  /* pure short */
  1102.          for(sfb=0,sfb_start=0,sfb_lines=sfBandIndex[sfreq].s[1]; 
  1103.             sfb < 13; sfb++,sfb_start=sfBandIndex[sfreq].s[sfb],
  1104.             (sfb_lines=sfBandIndex[sfreq].s[sfb+1] - sfb_start))
  1105.                for(window=0; window<3; window++)
  1106.                   for(freq=0;freq<sfb_lines;freq++) {
  1107.                      src_line = sfb_start*3 + window*sfb_lines + freq; 
  1108.                      des_line = (sfb_start*3) + window + (freq*3);
  1109.                      ro[des_line/SSLIMIT][des_line%SSLIMIT] = 
  1110.                                     xr[src_line/SSLIMIT][src_line%SSLIMIT];
  1111.                }
  1112.       }
  1113.    }
  1114.    else {   /*long blocks */
  1115.       for (sb=0 ; sb < SBLIMIT ; sb++)
  1116.          for (ss=0 ; ss < SSLIMIT ; ss++) 
  1117.             ro[sb][ss] = xr[sb][ss];
  1118.    }
  1119. }
  1120.  
  1121.  
  1122.  
  1123. void III_stereo(xr, lr, scalefac, gr_info, fr_ps)
  1124. double xr[2][SBLIMIT][SSLIMIT];
  1125. double lr[2][SBLIMIT][SSLIMIT];
  1126. III_scalefac_t *scalefac;
  1127. struct gr_info_s *gr_info;
  1128. frame_params *fr_ps;
  1129. {
  1130.    int sfreq = fr_ps->header->sampling_frequency;
  1131.    int stereo = fr_ps->stereo;
  1132.    int ms_stereo = (fr_ps->header->mode == MPG_MD_JOINT_STEREO) &&
  1133.                    (fr_ps->header->mode_ext & 0x2); 
  1134.    int i_stereo = (fr_ps->header->mode == MPG_MD_JOINT_STEREO) &&
  1135.                   (fr_ps->header->mode_ext & 0x1);
  1136.    int js_bound;  /* frequency line that marks the beggining of the zero part */  
  1137.    int sfb,next_sfb_boundary;
  1138.    int i,j,sb,ss,ch,is_pos[576]; 
  1139.    double is_ratio[576];
  1140.   
  1141.    /* intialization */
  1142.    for ( i=0; i<576; i++ )
  1143.       is_pos[i] = 7;
  1144.  
  1145.    if ((stereo == 2) && i_stereo )
  1146.    {  if (gr_info->window_switching_flag && (gr_info->block_type == 2))
  1147.       {  if( gr_info->mixed_block_flag )
  1148.          {  int max_sfb = 0;
  1149.  
  1150.             for ( j=0; j<3; j++ )
  1151.             {  int sfbcnt;
  1152.                sfbcnt = 2;
  1153.                for( sfb=12; sfb >=3; sfb-- )
  1154.                {  int lines;
  1155.                   lines = sfBandIndex[sfreq].s[sfb+1]-sfBandIndex[sfreq].s[sfb];
  1156.                   i = 3*sfBandIndex[sfreq].s[sfb] + (j+1) * lines - 1;
  1157.                   while ( lines > 0 )
  1158.                   {  if ( xr[1][i/SSLIMIT][i%SSLIMIT] != 0.0 )
  1159.                      {  sfbcnt = sfb;
  1160.                         sfb = -10;
  1161.                         lines = -10;
  1162.                      }
  1163.                      lines--;
  1164.                      i--;
  1165.                   }
  1166.                }
  1167.                sfb = sfbcnt + 1;
  1168.  
  1169.                if ( sfb > max_sfb )
  1170.                   max_sfb = sfb;
  1171.  
  1172.                while( sfb<12 )
  1173.                {  sb = sfBandIndex[sfreq].s[sfb+1]-sfBandIndex[sfreq].s[sfb];
  1174.                   i = 3*sfBandIndex[sfreq].s[sfb] + j * sb;
  1175.                   for ( ; sb > 0; sb--)
  1176.                   {  is_pos[i] = (*scalefac)[1].s[j][sfb];
  1177.                      if ( is_pos[i] != 7 )
  1178.                         is_ratio[i] = tan( is_pos[i] * (PI / 12));
  1179.                      i++;
  1180.                   }
  1181.                   sfb++;
  1182.                }
  1183.                sb = sfBandIndex[sfreq].s[11]-sfBandIndex[sfreq].s[10];
  1184.                sfb = 3*sfBandIndex[sfreq].s[10] + j * sb;
  1185.                sb = sfBandIndex[sfreq].s[12]-sfBandIndex[sfreq].s[11];
  1186.                i = 3*sfBandIndex[sfreq].s[11] + j * sb;
  1187.                for ( ; sb > 0; sb-- )
  1188.                {  is_pos[i] = is_pos[sfb];
  1189.                   is_ratio[i] = is_ratio[sfb];
  1190.                   i++;
  1191.                }
  1192.              }
  1193.              if ( max_sfb <= 3 )
  1194.              {  i = 2;
  1195.                 ss = 17;
  1196.                 sb = -1;
  1197.                 while ( i >= 0 )
  1198.                 {  if ( xr[1][i][ss] != 0.0 )
  1199.                    {  sb = i*18+ss;
  1200.                       i = -1;
  1201.                    } else
  1202.                    {  ss--;
  1203.                       if ( ss < 0 )
  1204.                       {  i--;
  1205.                          ss = 17;
  1206.                       }
  1207.                    }
  1208.                 }
  1209.                 i = 0;
  1210.                 while ( sfBandIndex[sfreq].l[i] <= sb )
  1211.                    i++;
  1212.                 sfb = i;
  1213.                 i = sfBandIndex[sfreq].l[i];
  1214.                 for ( ; sfb<8; sfb++ )
  1215.                 {  sb = sfBandIndex[sfreq].l[sfb+1]-sfBandIndex[sfreq].l[sfb];
  1216.                    for ( ; sb > 0; sb--)
  1217.                    {  is_pos[i] = (*scalefac)[1].l[sfb];
  1218.                       if ( is_pos[i] != 7 )
  1219.                          is_ratio[i] = tan( is_pos[i] * (PI / 12));
  1220.                       i++;
  1221.                    }
  1222.                 }
  1223.             }
  1224.          } else
  1225.          {  for ( j=0; j<3; j++ )
  1226.             {  int sfbcnt;
  1227.                sfbcnt = -1;
  1228.                for( sfb=12; sfb >=0; sfb-- )
  1229.                {  int lines;
  1230.                   lines = sfBandIndex[sfreq].s[sfb+1]-sfBandIndex[sfreq].s[sfb];
  1231.                   i = 3*sfBandIndex[sfreq].s[sfb] + (j+1) * lines - 1;
  1232.                   while ( lines > 0 )
  1233.                   {  if ( xr[1][i/SSLIMIT][i%SSLIMIT] != 0.0 )
  1234.                      {  sfbcnt = sfb;
  1235.                         sfb = -10;
  1236.                         lines = -10;
  1237.                      }
  1238.                      lines--;
  1239.                      i--;
  1240.                   }
  1241.                }
  1242.                sfb = sfbcnt + 1;
  1243.                while( sfb<12 )
  1244.                {  sb = sfBandIndex[sfreq].s[sfb+1]-sfBandIndex[sfreq].s[sfb];
  1245.                   i = 3*sfBandIndex[sfreq].s[sfb] + j * sb;
  1246.                   for ( ; sb > 0; sb--)
  1247.                   {  is_pos[i] = (*scalefac)[1].s[j][sfb];
  1248.                      if ( is_pos[i] != 7 )
  1249.                         is_ratio[i] = tan( is_pos[i] * (PI / 12));
  1250.                      i++;
  1251.                   }
  1252.                   sfb++;
  1253.                }
  1254.  
  1255.                sb = sfBandIndex[sfreq].s[11]-sfBandIndex[sfreq].s[10];
  1256.                sfb = 3*sfBandIndex[sfreq].s[10] + j * sb;
  1257.                sb = sfBandIndex[sfreq].s[12]-sfBandIndex[sfreq].s[11];
  1258.                i = 3*sfBandIndex[sfreq].s[11] + j * sb;
  1259.                for ( ; sb > 0; sb-- )
  1260.                {  is_pos[i] = is_pos[sfb];
  1261.                   is_ratio[i] = is_ratio[sfb];
  1262.                   i++;
  1263.                }
  1264.             }
  1265.          }
  1266.       } else
  1267.       {  i = 31;
  1268.          ss = 17;
  1269.          sb = 0;
  1270.          while ( i >= 0 )
  1271.          {  if ( xr[1][i][ss] != 0.0 )
  1272.             {  sb = i*18+ss;
  1273.                i = -1;
  1274.             } else
  1275.             {  ss--;
  1276.                if ( ss < 0 )
  1277.                {  i--;
  1278.                   ss = 17;
  1279.                }
  1280.             }
  1281.          }
  1282.          i = 0;
  1283.          while ( sfBandIndex[sfreq].l[i] <= sb )
  1284.             i++;
  1285.          sfb = i;
  1286.          i = sfBandIndex[sfreq].l[i];
  1287.          for ( ; sfb<21; sfb++ )
  1288.          {  sb = sfBandIndex[sfreq].l[sfb+1] - sfBandIndex[sfreq].l[sfb];
  1289.             for ( ; sb > 0; sb--)
  1290.             {  is_pos[i] = (*scalefac)[1].l[sfb];
  1291.                if ( is_pos[i] != 7 )
  1292.                   is_ratio[i] = tan( is_pos[i] * (PI / 12));
  1293.                i++;
  1294.             }
  1295.          }
  1296.          sfb = sfBandIndex[sfreq].l[20];
  1297.          for ( sb = 576 - sfBandIndex[sfreq].l[21]; sb > 0; sb-- )
  1298.          {  is_pos[i] = is_pos[sfb];
  1299.             is_ratio[i] = is_ratio[sfb];
  1300.             i++;
  1301.          }
  1302.       }
  1303.    }
  1304.  
  1305.    for(ch=0;ch<2;ch++)
  1306.       for(sb=0;sb<SBLIMIT;sb++)
  1307.          for(ss=0;ss<SSLIMIT;ss++) 
  1308.             lr[ch][sb][ss] = 0;
  1309.  
  1310.    if (stereo==2) 
  1311.       for(sb=0;sb<SBLIMIT;sb++)
  1312.          for(ss=0;ss<SSLIMIT;ss++) {
  1313.             i = (sb*18)+ss;
  1314.             if ( is_pos[i] == 7 ) {
  1315.                if ( ms_stereo ) {
  1316.                   lr[0][sb][ss] = (xr[0][sb][ss]+xr[1][sb][ss])/1.41421356;
  1317.                   lr[1][sb][ss] = (xr[0][sb][ss]-xr[1][sb][ss])/1.41421356;
  1318.                }
  1319.                else {
  1320.                   lr[0][sb][ss] = xr[0][sb][ss];
  1321.                   lr[1][sb][ss] = xr[1][sb][ss];
  1322.                }
  1323.             }
  1324.             else if (i_stereo ) {
  1325.                lr[0][sb][ss] = xr[0][sb][ss] * (is_ratio[i]/(1+is_ratio[i]));
  1326.                lr[1][sb][ss] = xr[0][sb][ss] * (1/(1+is_ratio[i])); 
  1327.             }
  1328.             else {
  1329.                printf("Error in streo processing\n");
  1330.             }
  1331.          }
  1332.    else  /* mono , bypass xr[0][][] to lr[0][][]*/
  1333.       for(sb=0;sb<SBLIMIT;sb++)
  1334.          for(ss=0;ss<SSLIMIT;ss++)
  1335.             lr[0][sb][ss] = xr[0][sb][ss];
  1336.  
  1337. }
  1338.  
  1339. double Ci[8]={-0.6,-0.535,-0.33,-0.185,-0.095,-0.041,-0.0142,-0.0037};
  1340.  
  1341. void III_antialias(xr, hybridIn, gr_info, fr_ps)
  1342. double xr[SBLIMIT][SSLIMIT];    
  1343. double hybridIn[SBLIMIT][SSLIMIT];
  1344. struct gr_info_s *gr_info;             
  1345. frame_params *fr_ps;            
  1346. {
  1347.    static int    init = 1;
  1348.    static double ca[8],cs[8];
  1349.    double        bu,bd;  /* upper and lower butterfly inputs */
  1350.    int           ss,sb,sblim;
  1351.  
  1352.    if (init) {
  1353.       int i;
  1354.       double    sq;
  1355.       for (i=0;i<8;i++) {
  1356.          sq=sqrt(1.0+Ci[i]*Ci[i]);
  1357.          cs[i] = 1.0/sq;
  1358.          ca[i] = Ci[i]/sq;
  1359.       }
  1360.       init = 0;
  1361.    }
  1362.    
  1363.    /* clear all inputs */  
  1364.       
  1365.     for(sb=0;sb<SBLIMIT;sb++)
  1366.        for(ss=0;ss<SSLIMIT;ss++)
  1367.           hybridIn[sb][ss] = xr[sb][ss];
  1368.  
  1369.    if  (gr_info->window_switching_flag && (gr_info->block_type == 2) &&
  1370.        !gr_info->mixed_block_flag ) return;
  1371.  
  1372.    if ( gr_info->window_switching_flag && gr_info->mixed_block_flag &&
  1373.      (gr_info->block_type == 2))
  1374.       sblim = 1;
  1375.    else
  1376.       sblim = SBLIMIT-1;
  1377.  
  1378.    /* 31 alias-reduction operations between each pair of sub-bands */
  1379.    /* with 8 butterflies between each pair                         */
  1380.  
  1381.    for(sb=0;sb<sblim;sb++)   
  1382.       for(ss=0;ss<8;ss++) {      
  1383.          bu = xr[sb][17-ss];
  1384.          bd = xr[sb+1][ss];
  1385.          hybridIn[sb][17-ss] = (bu * cs[ss]) - (bd * ca[ss]);
  1386.          hybridIn[sb+1][ss] = (bd * cs[ss]) + (bu * ca[ss]);
  1387.          }  
  1388. }
  1389.  
  1390. void inv_mdct(in, out, block_type)
  1391. double in[18];
  1392. double out[36];
  1393. int block_type;
  1394. {
  1395. /*------------------------------------------------------------------*/
  1396. /*                                                                  */
  1397. /*    Function: Calculation of the inverse MDCT                     */
  1398. /*    In the case of short blocks the 3 output vectors are already  */
  1399. /*    overlapped and added in this modul.                           */
  1400. /*                                                                  */
  1401. /*    New layer3                                                    */
  1402. /*                                                                  */
  1403. /*------------------------------------------------------------------*/
  1404.  
  1405. int     k,i,m,N,p;
  1406. double  tmp[12],sum;
  1407. static  double  win[4][36];
  1408. static  int init=0;
  1409. static  double COS[4*36];
  1410.  
  1411.     if(init==0){
  1412.  
  1413.     /* type 0 */
  1414.       for(i=0;i<36;i++)
  1415.          win[0][i] = sin( PI/36 *(i+0.5) );
  1416.  
  1417.     /* type 1*/
  1418.       for(i=0;i<18;i++)
  1419.          win[1][i] = sin( PI/36 *(i+0.5) );
  1420.       for(i=18;i<24;i++)
  1421.          win[1][i] = 1.0;
  1422.       for(i=24;i<30;i++)
  1423.          win[1][i] = sin( PI/12 *(i+0.5-18) );
  1424.       for(i=30;i<36;i++)
  1425.          win[1][i] = 0.0;
  1426.  
  1427.     /* type 3*/
  1428.       for(i=0;i<6;i++)
  1429.          win[3][i] = 0.0;
  1430.       for(i=6;i<12;i++)
  1431.          win[3][i] = sin( PI/12 *(i+0.5-6) );
  1432.       for(i=12;i<18;i++)
  1433.          win[3][i] =1.0;
  1434.       for(i=18;i<36;i++)
  1435.          win[3][i] = sin( PI/36*(i+0.5) );
  1436.  
  1437.     /* type 2*/
  1438.       for(i=0;i<12;i++)
  1439.          win[2][i] = sin( PI/12*(i+0.5) ) ;
  1440.       for(i=12;i<36;i++)
  1441.          win[2][i] = 0.0 ;
  1442.  
  1443.       for (i=0; i<4*36; i++)
  1444.          COS[i] = cos(PI/(2*36) * i);
  1445.  
  1446.       init++;
  1447.     }
  1448.  
  1449.     for(i=0;i<36;i++)
  1450.        out[i]=0;
  1451.  
  1452.     if(block_type == 2){
  1453.        N=12;
  1454.        for(i=0;i<3;i++){
  1455.           for(p= 0;p<N;p++){
  1456.              sum = 0.0;
  1457.              for(m=0;m<N/2;m++)
  1458.                 sum += in[i+3*m] * cos( PI/(2*N)*(2*p+1+N/2)*(2*m+1) );
  1459.              tmp[p] = sum * win[block_type][p] ;
  1460.           }
  1461.           for(p=0;p<N;p++)
  1462.              out[6*i+p+6] += tmp[p];
  1463.        }
  1464.     }
  1465.     else{
  1466.       N=36;
  1467.       for(p= 0;p<N;p++){
  1468.          sum = 0.0;
  1469.          for(m=0;m<N/2;m++)
  1470.            sum += in[m] * COS[((2*p+1+N/2)*(2*m+1))%(4*36)];
  1471.          out[p] = sum * win[block_type][p];
  1472.       }
  1473.     }
  1474. }
  1475.  
  1476. void III_hybrid(fsIn, tsOut ,sb, ch, gr_info, fr_ps)
  1477. double fsIn[SSLIMIT];   /* freq samples per subband in */
  1478. double tsOut[SSLIMIT];  /* time samples per subband out */
  1479. int sb, ch;
  1480. struct gr_info_s *gr_info;             
  1481. frame_params *fr_ps;            
  1482. {
  1483.    int ss;
  1484.    double rawout[36];
  1485.    static double prevblck[2][SBLIMIT][SSLIMIT];
  1486.    static int init = 1;
  1487.    int bt;
  1488.  
  1489.    if (init) {
  1490.       int i,j,k;
  1491.       
  1492.       for(i=0;i<2;i++)
  1493.          for(j=0;j<SBLIMIT;j++)
  1494.             for(k=0;k<SSLIMIT;k++)
  1495.                prevblck[i][j][k]=0.0;
  1496.       init = 0;
  1497.    }
  1498.  
  1499.    bt = (gr_info->window_switching_flag && gr_info->mixed_block_flag &&
  1500.           (sb < 2)) ? 0 : gr_info->block_type; 
  1501.  
  1502.    inv_mdct( fsIn, rawout, bt);
  1503.  
  1504.    /* overlap addition */
  1505.    for(ss=0; ss<SSLIMIT; ss++) {
  1506.       tsOut[ss] = rawout[ss] + prevblck[ch][sb][ss];
  1507.       prevblck[ch][sb][ss] = rawout[ss+18];
  1508.    }
  1509. }
  1510.  
  1511. /* Return the number of slots for main data of current frame, */
  1512.  
  1513. int main_data_slots(fr_ps)
  1514. frame_params fr_ps;
  1515. {int nSlots;
  1516.  
  1517.    nSlots = (144 * bitrate[2][fr_ps.header->bitrate_index])
  1518.         / s_freq[fr_ps.header->sampling_frequency];
  1519.   if (fr_ps.header->padding) nSlots++;
  1520.   nSlots -= 4;
  1521.   if (fr_ps.header->error_protection) nSlots -= 2;
  1522.   if (fr_ps.stereo == 1) nSlots -= 17; else nSlots -=32;
  1523.   return(nSlots);
  1524. }
  1525.