newmdct.c

MP3编码程序和资料

/*   transformations with short length. This leads to 12 coefficents */
/*   in the time and 6 in the frequency domain. In this case the     */
/*   results are stored side by side in the vector out[].            */
/*                                                                   */
/*   New layer3                                                      */
/*                                                                   */
/*-------------------------------------------------------------------*/

INLINE static void mdct_short(FLOAT8 *out, FLOAT8 *in)
{
    int l;
    for ( l = 0; l < 3; l++ ) {
	FLOAT tc0,tc1,tc2,ts0,ts1,ts2;

	ts0 = in[5];
	tc0 = in[3];
	tc1 = ts0 + tc0;
	tc2 = ts0 - tc0;

	ts0 = in[2];
	tc0 = in[0];
	ts1 = ts0 + tc0;
	ts2 =-ts0 + tc0;

	tc0 = in[4];
	ts0 = in[1];

	out[3*0] = tc1 + tc0;
	out[3*5] =-ts1 + ts0;

	tc2 = tc2 * 0.86602540378443870761;
	ts1 = ts1 * 0.5 + ts0;
	out[3*1] = tc2-ts1;
	out[3*2] = tc2+ts1;

	tc1 = tc1 * 0.5 - tc0;
	ts2 = ts2 * 0.86602540378443870761;
	out[3*3] = tc1+ts2;
	out[3*4] = tc1-ts2;

	in += 6; out++;
    }
}

INLINE static void mdct_long(FLOAT8 *out, FLOAT8 *in)
{
#define inc(x) in[17-(x)]
#define ins(x) in[8-(x)]

    const FLOAT8 c0=0.98480775301220802032, c1=0.64278760968653936292, c2=0.34202014332566882393;
    const FLOAT8 c3=0.93969262078590842791, c4=-0.17364817766693030343, c5=-0.76604444311897790243;
    FLOAT8 tc1 = inc(0)-inc(8),tc2 = (inc(1)-inc(7))*0.86602540378443870761, tc3 = inc(2)-inc(6), tc4 = inc(3)-inc(5);
    FLOAT8 tc5 = inc(0)+inc(8),tc6 = (inc(1)+inc(7))*0.5,                    tc7 = inc(2)+inc(6), tc8 = inc(3)+inc(5);
    FLOAT8 ts1 = ins(0)-ins(8),ts2 = (ins(1)-ins(7))*0.86602540378443870761, ts3 = ins(2)-ins(6), ts4 = ins(3)-ins(5);
    FLOAT8 ts5 = ins(0)+ins(8),ts6 = (ins(1)+ins(7))*0.5,                    ts7 = ins(2)+ins(6), ts8 = ins(3)+ins(5);
    FLOAT8 ct,st;

    ct = tc5+tc7+tc8+inc(1)+inc(4)+inc(7);
    out[0] = ct;

    ct = tc1*c0 + tc2 + tc3*c1 + tc4*c2;
    st = -ts5*c4 + ts6 - ts7*c5 + ts8*c3 + ins(4);
    out[1] = ct+st;
    out[2] = ct-st;

    ct =  tc5*c3 + tc6 + tc7*c4 + tc8*c5 - inc(4);
    st = ts1*c2 + ts2 + ts3*c0 + ts4*c1;
    out[3] = ct+st;
    out[4] = ct-st;

    ct = (tc1-tc3-tc4)*0.86602540378443870761;
    st = (ts5+ts7-ts8)*0.5+ins(1)-ins(4)+ins(7);
    out[5] = ct+st;
    out[6] = ct-st;

    ct = -tc5*c5 - tc6 - tc7*c3 - tc8*c4 + inc(4);
    st = ts1*c1 + ts2 - ts3*c2 - ts4*c0;
    out[7] = ct+st;
    out[8] = ct-st;

    ct = tc1*c1 - tc2 - tc3*c2 + tc4*c0;
    st = -ts5*c5 + ts6 - ts7*c3 + ts8*c4 + ins(4);
    out[ 9] = ct+st;
    out[10] = ct-st;

    ct = (tc5+tc7+tc8)*0.5-inc(1)-inc(4)-inc(7);
    st = (ts1-ts3+ts4)*0.86602540378443870761;
    out[11] = ct+st;
    out[12] = ct-st;

    ct = tc1*c2 - tc2 + tc3*c0 - tc4*c1;
    st =  ts5*c3 - ts6 + ts7*c4 - ts8*c5 - ins(4);
    out[13] = ct+st;
    out[14] = ct-st;

    ct = -tc5*c4 - tc6 - tc7*c5 - tc8*c3 + inc(4);
    st = ts1*c0 - ts2 + ts3*c1 - ts4*c2;
    out[15] = ct+st;
    out[16] = ct-st;

    st = ts5+ts7-ts8-ins(1)+ins(4)-ins(7);
    out[17] = st;
}

static const int order[] = {
    0,  16,  8, 24,  4,  20,  12,  28,
    2,  18, 10, 26,  6,  22,  14,  30,
    1,  17,  9, 25,  5,  21,  13,  29,
    3,  19, 11, 27,  7,  23,  15,  31
};



void mdct_init48(lame_global_flags *gfp)
{
    int i, k;
    FLOAT8 sq;

    /* prepare the aliasing reduction butterflies */
    for (k = 0; k < 8; k++) {
	/*
	  This is table B.9: coefficients for aliasing reduction
	  */
	static const FLOAT8 c[8] = {
	    -0.6,-0.535,-0.33,-0.185,-0.095,-0.041,-0.0142, -0.0037
	};
	sq = 1.0 + c[k] * c[k];
	sq = sqrt(sq);
	ca[k] = c[k] / sq;
	cs[k] = 1.0 / sq;
    }

    /* type 0*/
    for (i = 0; i < 36; i++)
	win[0][i] = sin(PI/36 * (i + 0.5));
    /* type 1*/
    for (i = 0; i < 18; i++) 
	win[1][i] = win[0][i];
    for (; i < 24; i++)
	win[1][i] = 1.0;
    for (; i < 30; i++)
	win[1][i] = cos(PI/12 * (i + 0.5));
    for (; i < 36; i++)
	win[1][i] = 0.0;
    /* type 3*/
    for (i = 0; i < 36; i++)
	win[3][i] = win[1][35 - i];

    /* swap window data*/
    for (k = 0; k < 4; k++) {
	FLOAT8 a;

	a = win[0][17-k];
	win[0][17-k] = win[0][9+k];
	win[0][9+k] = a;

	a = win[0][35-k];
	win[0][35-k] = win[0][27+k];
	win[0][27+k] = a;

	a = win[1][17-k];
	win[1][17-k] = win[1][9+k];
	win[1][9+k] = a;

	a = win[1][35-k];
	win[1][35-k] = win[1][27+k];
	win[1][27+k] = a;

	a = win[3][17-k];
	win[3][17-k] = win[3][9+k];
	win[3][9+k] = a;

	a = win[3][35-k];
	win[3][35-k] = win[3][27+k];
	win[3][27+k] = a;
    }

    for (i = 0; i < NL; i++) {
	win[0][i] *= cos((NL/4+0.5+i%9)*PI/NL) / SCALE/(NL/4);
	win[1][i] *= cos((NL/4+0.5+i%9)*PI/NL) / SCALE/(NL/4);
	win[3][i] *= cos((NL/4+0.5+i%9)*PI/NL) / SCALE/(NL/4);
    }

    for (i = NL/2; i < NL; i++) {
	win[0][i] *= -1;
	win[1][i] *= -1;
	win[3][i] *= -1;
    }

    for (i = 0; i < NL/4; i++)
	tantab_l[i] = tan((NL/4+0.5+i)*PI/NL);

    /* type 2(short)*/
    for (i = 0; i < NS / 4; i++) {
	FLOAT8 w2 = cos(PI / NS * (i + 0.5)) * (4.0/NS) / SCALE;
	win[SHORT_TYPE][i] = tan(PI / NS * (i + 0.5));
	tritab_s[i*2  ] = cos((0.5+2-i)*PI/NS) * w2;
	tritab_s[i*2+1] = sin((0.5+2-i)*PI/NS) * w2;
    }
}

void mdct_sub48(lame_global_flags *gfp,
    short *w0, short *w1,
    FLOAT8 mdct_freq[2][2][576],
    III_side_info_t *l3_side)
{
    int gr, k, ch;
    short *wk;
    static int init = 0;
    lame_internal_flags *gfc=gfp->internal_flags;

    FLOAT8 work[18];

    if ( gfc->mdct_sub48_init == 0 ) {
        gfc->mdct_sub48_init=1;
	mdct_init48(gfp);
	init++;
    }

    wk = w0 + 286;
    /* thinking cache performance, ch->gr loop is better than gr->ch loop */
    for (ch = 0; ch < gfc->stereo; ch++) {
	for (gr = 0; gr < gfc->mode_gr; gr++) {
	    int	band;
	    FLOAT8 *mdct_enc = &mdct_freq[gr][ch][0];
	    gr_info *gi = &(l3_side->gr[gr].ch[ch].tt);
	    FLOAT8 *samp = gfc->sb_sample[ch][1 - gr][0];

	    for (k = 0; k < 18 / 2; k++) {
		window_subband(wk, samp + 16);
		idct32(samp);
		window_subband(wk + 32, samp + 32+16);
		idct32(samp+32);
		samp += 64;
		wk += 64;
		/*
		 * Compensate for inversion in the analysis filter
		 */
		for (band = 16-32; band < 0; band++)
		    samp[band] *= -1;
	    }


	    /* apply filters on the polyphase filterbank outputs */
	    /* bands <= gfc->highpass_band will be zeroed out below */
	    /* bands >= gfc->lowpass_band  will be zeroed out below */
	    if (gfc->filter_type==0) {
              for (band=gfc->highpass_start_band;  band <= gfc->highpass_end_band; band++) { 
		  for (k=0; k<18; k++) 
		    gfc->sb_sample[ch][1-gr][k][order[band]]*=gfc->amp_highpass[band];
	      }
              for (band=gfc->lowpass_start_band;  band <= gfc->lowpass_end_band; band++) { 
		  for (k=0; k<18; k++) 
		    gfc->sb_sample[ch][1-gr][k][order[band]]*=gfc->amp_lowpass[band];
	      }
	    }
	    


	    /*
	     * Perform imdct of 18 previous subband samples
	     * + 18 current subband samples
	     */
	    for (band = 0; band < 32; band++, mdct_enc += 18) 
	    {
		int type = gi->block_type;
		int band_swapped;
		band_swapped = order[band];
#ifdef ALLOW_MIXED
		if (gi->mixed_block_flag && band < 2)
		    type = 0;
#endif
		if (band >= gfc->lowpass_band || band <= gfc->highpass_band) {
		    memset((char *)mdct_enc,0,18*sizeof(FLOAT8));
		}else {
		  if (type == SHORT_TYPE) {
		    for (k = 2; k >= 0; --k) {
			FLOAT8 win1 = win[SHORT_TYPE][k];
			FLOAT8 a, b;

			a = gfc->sb_sample[ch][gr][k+6][band_swapped] * win1 -
			    gfc->sb_sample[ch][gr][11-k][band_swapped];

			b = gfc->sb_sample[ch][gr][k+12][band_swapped] +
			    gfc->sb_sample[ch][gr][17-k][band_swapped] * win1;

			work[k+3] = -b*tritab_s[k*2  ] + a * tritab_s[k*2+1];
			work[k  ] =  b*tritab_s[k*2+1] + a * tritab_s[k*2  ];

			a = gfc->sb_sample[ch][gr][k+12][band_swapped] * win1 -
			    gfc->sb_sample[ch][gr][17-k][band_swapped];

			b = gfc->sb_sample[ch][1-gr][k][band_swapped] +
			    gfc->sb_sample[ch][1-gr][5-k][band_swapped] * win1;

			work[k+9] = -b*tritab_s[k*2  ] + a * tritab_s[k*2+1];
			work[k+6] =  b*tritab_s[k*2+1] + a * tritab_s[k*2  ];

			a = gfc->sb_sample[ch][1-gr][k][band_swapped] * win1 -
			    gfc->sb_sample[ch][1-gr][5-k][band_swapped];

			b = gfc->sb_sample[ch][1-gr][k+6][band_swapped] +
			    gfc->sb_sample[ch][1-gr][11-k][band_swapped] * win1;

			work[k+15] = -b*tritab_s[k*2  ] + a * tritab_s[k*2+1];
			work[k+12] =  b*tritab_s[k*2+1] + a * tritab_s[k*2  ];
		    }
		    mdct_short(mdct_enc, work);
		  } else {
		    for (k = -NL/4; k < 0; k++) {
			FLOAT8 a, b;
			a = win[type][k+27] * gfc->sb_sample[ch][1-gr][k+9][band_swapped]
			  + win[type][k+36] * gfc->sb_sample[ch][1-gr][8-k][band_swapped];
			b = win[type][k+ 9] * gfc->sb_sample[ch][gr][k+9][band_swapped]
			  - win[type][k+18] * gfc->sb_sample[ch][gr][8-k][band_swapped];
			work[k+ 9] = a - b*tantab_l[k+9];
			work[k+18] = a*tantab_l[k+9] + b;
		    }

		    mdct_long(mdct_enc, work);
		  }
		}
		
		
		/*
		  Perform aliasing reduction butterfly
		*/
		if (type != SHORT_TYPE) {
		  if (band == 0)
		    continue;
		  for (k = 7; k >= 0; --k) {
		    FLOAT8 bu,bd;
		    bu = mdct_enc[k] * ca[k] + mdct_enc[-1-k] * cs[k];
		    bd = mdct_enc[k] * cs[k] - mdct_enc[-1-k] * ca[k];
		    
		    mdct_enc[-1-k] = bu;
		    mdct_enc[k]    = bd;
		  }
		}
	      }
	}
	wk = w1 + 286;
	if (gfc->mode_gr == 1) {
	    memcpy(gfc->sb_sample[ch][0], gfc->sb_sample[ch][1], 576 * sizeof(FLOAT8));
	}
    }
}