l3loop.c

mp3 source code decoder & encoder

/****************************************************************/
/* Optimized by 8hz, May 1998                                   */
/* Added precalc tables and rewritten some heavy calculations   */
/****************************************************************/

/* MP3STEGO                                                     */
/* Some lines added or modified for the information hidding     */
/*                                                              */
/* Fabien A.P. Petitcolas, August 1998                          */
/*                                                              */
/* $Header: /MP3Stego/MP3Stego Encoder/l3loop.c 6     30/11/00 15:57 Fabienpe $                                                   */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include "../../stegolib/stego.h"

#include "types.h"
#include "error.h"
#include "tables.h"
#include "layer3.h"
#include "l3loop.h"
#include "huffman.h"
#include "l3bitstream.h"
#include "reservoir.h"
#include "assert.h"

int bin_search_StepSize(int desired_rate, double start, int ix[576],
           double xrs[576], gr_info * cod_info);

int count_bits();

float worst_xfsf_to_xmin_ratio();

#define PRECALC_SIZE 1024 /* WAS 256 !!! */
static double pow43[PRECALC_SIZE];
#include "sqrttab.h"

#ifndef HAVE_NINT
int nint(double in)
{
    if (in<0) return (int)(in - 0.5);
    else     return (int)(in + 0.5);
}

#endif


/*
  The following table is used to implement the scalefactor partitioning 
  for MPEG2 as described in section 2.4.3.2 of the IS. 
  The indexing corresponds to the way the tables are presented in the IS:
  [table_number][row_in_table][column of nr_of_sfb]
*/
static unsigned nr_of_sfb_block[6][3][4] =
{
  { {6, 5, 5, 5}, {9, 9, 9, 9}, {6, 9, 9, 9} },
  { {6, 5, 7, 3}, {9, 9, 12, 6}, {6, 9, 12, 6} },
  { {11, 10, 0, 0}, {18, 18, 0, 0}, {15,18,0,0} },
  { {7, 7, 7, 0}, {12, 12, 12, 0}, {6, 15, 12, 0} },
  { {6, 6, 6, 3}, {12, 9, 9, 6}, {6, 12, 9, 6} },
  { {8, 8, 5, 0}, {15,12,9,0}, {6,18,9,0} }
};

/* Table B.6: layer3 preemphasis */
int  pretab[21] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                    1, 1, 1, 1, 2, 2, 3, 3, 3, 2 };

/* This is the scfsi_band table from 2.4.2.7 of the IS */
int scfsi_band_long[5] = { 0, 6, 11, 16, 21 };

struct
{
    unsigned region0_count;
    unsigned region1_count;
} subdv_table[ 23 ] =
{
{0, 0}, /* 0 bands */
{0, 0}, /* 1 bands */
{0, 0}, /* 2 bands */
{0, 0}, /* 3 bands */
{0, 0}, /* 4 bands */
{0, 1}, /* 5 bands */
{1, 1}, /* 6 bands */
{1, 1}, /* 7 bands */
{1, 2}, /* 8 bands */
{2, 2}, /* 9 bands */
{2, 3}, /* 10 bands */
{2, 3}, /* 11 bands */
{3, 4}, /* 12 bands */
{3, 4}, /* 13 bands */
{3, 4}, /* 14 bands */
{4, 5}, /* 15 bands */
{4, 5}, /* 16 bands */
{4, 6}, /* 17 bands */
{5, 6}, /* 18 bands */
{5, 6}, /* 19 bands */
{5, 7}, /* 20 bands */
{6, 7}, /* 21 bands */
{6, 7}, /* 22 bands */
};

int *scalefac_band_long  = &sfBandIndex[3].l[0];
int *scalefac_band_short = &sfBandIndex[3].s[0];

int quantanf_init(double xr[576]);
void calc_scfsi(double  xr[576], L3_side_info_t *l3_side, L3_psy_xmin_t *l3_xmin, int ch, int gr);
int part2_length(L3_scalefac_t *scalefac, int gr, int ch, L3_side_info_t *si);
int scale_bitcount(L3_scalefac_t *scalefac, gr_info *cod_info, int gr, int ch);
void calc_noise(double xr[576], int ix[576], gr_info *cod_info, double xfsf[4][CBLIMIT]);

int bin_search_StepSize(int desired_rate, double start, int *ix, double xrs[576], gr_info * cod_info);
int count_bits(int *ix /*int[576]*/, gr_info *cod_info);
void gr_deco(gr_info *cod_info);
int count_bit(int ix[576], unsigned int start, unsigned int end, unsigned int table);
int bigv_bitcount(int ix[576], gr_info *gi);
int choose_table(int max);
int new_choose_table(int ix[576], unsigned int begin, unsigned int end);
void bigv_tab_select(int ix[576], gr_info *cod_info);
void subdivide(gr_info *cod_info);
int count1_bitcount(int ix[ 576 ], gr_info *cod_info);
void calc_runlen(int ix[576], gr_info *cod_info);

void calc_xmin(double xr[2][2][576], L3_psy_ratio_t *ratio, gr_info *cod_info, L3_psy_xmin_t *l3_xmin, int gr, int ch);
int loop_break(L3_scalefac_t *scalefac, gr_info *cod_info, int gr, int ch);
void preemphasis(double xr[576], double xfsf[4][CBLIMIT], L3_psy_xmin_t  *l3_xmin, int gr, int ch, L3_side_info_t *l3_side);
int amp_scalefac_bands(double xr[576], double xfsf[4][CBLIMIT], L3_psy_xmin_t    *l3_xmin, L3_side_info_t *l3_side, L3_scalefac_t  *scalefac, int gr, int ch, int iteration);
void quantize(double xr[576], int ix[576], gr_info *cod_info);
int ix_max(int ix[576], unsigned int begin, unsigned int end);
double xr_max(double xr[576], unsigned int begin, unsigned int end);


void L3_loop_initialise()
{
    long i;

    for(i=0;i<PRECALC_SIZE;i++)
        pow43[i] = pow((double)i, 4.0/3.0);
}


/***************************************************************************/ 
/* The code selects the best quantizerStepSize for a particular set        */
/* of scalefacs                                                            */
/***************************************************************************/ 
/* MP3STEGO-> */
static int inner_loop(double xr[2][2][576],  int l3_enc[2][2][576], 
                      int max_bits, gr_info *cod_info, int gr, int ch ,
                      int hiddenBit, int part2length)
/* ->MP3STEGO */
{
    int bits, c1bits, bvbits;
    double *xrs;  /*  double[576] *xr; */
    int     *ix;  /*  int[576]    *ix; */
	int    embedRule = 0;

    xrs = &xr[gr][ch][0];
    ix  = l3_enc[gr][ch];

    if (max_bits<0) ERROR("Ehhh !?!, negative compression !?!");
    cod_info->quantizerStepSize -= 1.0;;
    do
    {
        do
        {
            cod_info->quantizerStepSize += 1.0;
            quantize(xrs,ix,cod_info);
        }
        while(ix_max(ix,0,576)>(8191+14));               /* within table range? */

        calc_runlen(ix,cod_info);                        /* rzero,count1,big_values*/
        bits = c1bits = count1_bitcount(ix,cod_info);    /* count1_table selection*/
        subdivide(cod_info);                             /* bigvalues sfb division */
        bigv_tab_select(ix,cod_info);                    /* codebook selection*/
        bits += bvbits = bigv_bitcount(ix, cod_info);  /* bit count */

/* MP3STEGO-> */
		switch (hiddenBit)
		{
		case 2:
			embedRule = 0;
			break;
		case 0:
		case 1:
			embedRule = ((bits + part2length) % 2) != hiddenBit;
			break;
		default:
			ERROR("inner_loop: unexpected hidden bit.");
		}
/* ->MP3STEGO */
    }
/* MP3STEGO-> */
    while((bits>max_bits) | embedRule);
/* ->MP3STEGO */
    return bits;
}



/************************************************************************/
/*  Function: The outer iteration loop controls the masking conditions  */
/*  of all scalefactorbands. It computes the best scalefac and          */
/*  global gain. This module calls the inner iteration loop             */
/************************************************************************/
static int outer_loop(double xr[2][2][576],     /*  magnitudes of the spectral values */
                       int max_bits,
                       L3_psy_xmin_t *l3_xmin,  /* the allowed distortion of the scalefactor */
                       int l3_enc[2][2][576],   /* vector of quantized values ix(0..575) */
                       L3_scalefac_t *scalefac, /* scalefactors */
                       int gr, int ch,
                       L3_side_info_t *side_info,
					   int hiddenBit, int mean_bits)
{
    int status ;
    int scalesave_l[CBLIMIT], scalesave_s[CBLIMIT][3];
    int sfb, bits, huff_bits, save_preflag, save_compress;
    double xfsf[4][CBLIMIT];
    int i, over, iteration;

    double  *xrs; 
    int     *ix;  
    gr_info *cod_info = &side_info->gr[gr].ch[ch].tt;

    xrs = (double *) &(xr[gr][ch][0]); 
    ix  = (int *) &(l3_enc[gr][ch][0]);

    iteration = 0;
    do 
    {
		iteration += 1;
		cod_info->part2_length = part2_length(scalefac,gr,ch,side_info);
        huff_bits = max_bits - cod_info->part2_length;
		
		if (iteration == 1)
        {
			bin_search_StepSize(max_bits,cod_info->quantizerStepSize,
								ix,xrs,cod_info); /* speeds things up a bit */
		}
/* MP3STEGO-> */
		bits = inner_loop(xr, l3_enc, huff_bits, cod_info, gr, ch, 
            hiddenBit, cod_info->part2_length);
/* ->MP3STEGO */
		
        /* distortion calculation */
        calc_noise(&xr[gr][ch][0], &l3_enc[gr][ch][0], cod_info, xfsf);


        for (sfb = 0; sfb < CBLIMIT; sfb++) /* save scaling factors */
			scalesave_l[sfb] = scalefac->l[gr][ch][sfb];

        for (sfb = 0; sfb < SFB_SMAX; sfb++)
			for (i = 0; i < 3; i++)
				scalesave_s[sfb][i] = scalefac->s[gr][ch][sfb][i];
        
        save_preflag  = cod_info->preflag;
        save_compress = cod_info->scalefac_compress;

        preemphasis(&xr[gr][ch][0],xfsf,l3_xmin,gr,ch,side_info);

        over = amp_scalefac_bands(&xr[gr][ch][0], xfsf, l3_xmin,
                                   side_info, scalefac, gr, ch, iteration);
		

		if ((status=loop_break(scalefac,cod_info,gr,ch))==0)
		{
			/*if (fr_ps->header->version == 1)                         */
			status = scale_bitcount(scalefac, cod_info, gr, ch);
			/*else                                                       */
			/*status = scale_bitcount_lsf(scalefac, cod_info, gr, ch); */
		}

	}
	while ((status==0)&&(over>0));

    cod_info->preflag = save_preflag;
    cod_info->scalefac_compress = save_compress;

    for(sfb=0;sfb<21;sfb++) scalefac->l[gr][ch][sfb] = scalesave_l[sfb];    

    for(i=0;i<3;i++)
        for(sfb=0;sfb<12;sfb++)
            scalefac->s[gr][ch][sfb][i] = scalesave_s[sfb][i];    

    cod_info->part2_length   = part2_length(scalefac,gr,ch,side_info);
    cod_info->part2_3_length = cod_info->part2_length + bits;

/* MP3STEGO-> */
	if (hiddenBit == 0 || hiddenBit == 1)
		assert(hiddenBit == (int)(cod_info->part2_3_length % 2));
/* ->MP3SETGO */
    return cod_info->part2_3_length;
}


/************************************************************************/
/*  iteration_loop()                                                    */
/************************************************************************/
void L3_iteration_loop(double           pe[][2], 
                       double           mdct_freq_org[2][2][576],
                       L3_psy_ratio_t   *ratio,
                       L3_side_info_t   *side_info, 
                       int              l3_enc[2][2][576],
                       int              mean_bits, 
                       L3_scalefac_t    *scalefactor) 
{
    L3_psy_xmin_t l3_xmin;
    gr_info *cod_info;
    int *main_data_begin;

    int max_bits;
    int ch, gr, sfb, i;
    static int firstcall = 1;

    double xr[2][2][576];

/* MP3STEGO-> */
	int hiddenBit;
/* ->MP3STEGO */

	main_data_begin = &side_info->main_data_begin;
    side_info->resvDrain = 0;

    if (firstcall)
    {
	*main_data_begin = 0;
        firstcall=0;
    }

    scalefac_band_long  = &sfBandIndex[config.mpeg.samplerate_index + (config.mpeg.type * 3)].l[0];
    scalefac_band_short = &sfBandIndex[config.mpeg.samplerate_index + (config.mpeg.type * 3)].s[0];

    for (gr = 0; gr < config.mpeg.mode_gr; gr++)
    {
        for (ch = 0; ch < config.wave.channels; ch++)
        {
            for (i = 0; i < 576; i++) 
                xr[gr][ch][i] = mdct_freq_org[gr][ch][i];
        }
    }

    ResvFrameBegin(side_info,mean_bits,config.mpeg.bits_per_frame);

    for(gr=0;gr<config.mpeg.mode_gr;gr++)
    {
        for(ch=0;ch<config.wave.channels;ch++)
        {
            cod_info = (gr_info *) &(side_info->gr[gr].ch[ch]);
            gr_deco(cod_info);
            calc_xmin(xr, ratio, cod_info, &l3_xmin, gr, ch);
	    
            if (config.mpeg.type==TYPE_MPEG_I) calc_scfsi(xr[gr][ch],side_info,&l3_xmin,ch,gr);
	    
/* calculation of number of available bit(per granule) */
		    max_bits = ResvMaxBits(side_info,&pe[gr][ch],mean_bits);
	    
/* reset of iteration variables */
            for(sfb=0;sfb<21;sfb++) scalefactor->l[gr][ch][sfb] = 0;
            for(sfb=0;sfb<13;sfb++)
                for (i=0;i<3;i++)
                    scalefactor->s[gr][ch][sfb][i] = 0;

			for (i = 0; i < 4; i++) cod_info->slen[i] = 0;
			cod_info->sfb_partition_table = &nr_of_sfb_block[0][0][0];
				cod_info->part2_3_length    = 0;
				cod_info->big_values        = 0;
				cod_info->count1            = 0;
				cod_info->scalefac_compress = 0;
				cod_info->table_select[0]   = 0;
				cod_info->table_select[1]   = 0;
				cod_info->table_select[2]   = 0;
				cod_info->subblock_gain[0]  = 0;
				cod_info->subblock_gain[1]  = 0;
				cod_info->subblock_gain[2]  = 0;
				cod_info->region0_count     = 0;
				cod_info->region1_count     = 0;
				cod_info->part2_length      = 0;
				cod_info->preflag           = 0;
				cod_info->scalefac_scale    = 0;
				cod_info->quantizerStepSize = 0.0;
				cod_info->count1table_select= 0;
			
/* MP3STEGO-> */
			hiddenBit = StegoGetNextBit();
#if defined(_DEBUG)
			fprintf(fLog, "Hidden bit: %d\n", hiddenBit);
#endif
/* ->MP3STEGO */

            /* all spectral values zero ? */
            if (fabs(xr_max(xr[gr][ch],0,576))!=0.0)
            {
                cod_info->quantizerStepSize = (double) quantanf_init(xr[gr][ch]);
/* MP3STEGO-> */
                cod_info->part2_3_length    = outer_loop(xr, max_bits, &l3_xmin, l3_enc,
                                                         scalefactor, gr, ch,side_info,
                                                         hiddenBit, mean_bits);
/* ->MP3STEGO */
            }
		
			ResvAdjust(cod_info, side_info, mean_bits);

			cod_info->global_gain = nint(cod_info->quantizerStepSize+210.0);

            if (cod_info->global_gain>=256) ERROR("gain>256, program error.");
        } /* for ch */
    } /* for gr */

    ResvFrameEnd(side_info,mean_bits);
}



/************************************************************************/
/*  quantanf_init                                                       */
/************************************************************************/
int quantanf_init(double xr[576])
/* Function: Calculate the first quantization step quantanf.       */
{
    int i, tp = 0;
    double system_const, minlimit;
    double sfm = 0.0, sum1 = 0.0, sum2 = 0.0;
    
    system_const = 8.0;
    minlimit     = -100.0;

    for (i = 0; i < 576; i++)
        if (xr[i] != 0)
	{
            double tpd = xr[i]*xr[i];
            sum1 += log(tpd);
            sum2 += tpd;
        }

    if (sum2 != 0.0)
    {
        sfm = exp(sum1/576.0)/(sum2/576.0);
        tp  = nint(system_const*log(sfm));
	if (tp<minlimit) tp = minlimit;
    }

    return(tp-70.0); /* SS 19-12-96. Starting value of
                        global_gain or quantizerStepSize 
                        has to be reduced for iteration_loop