newmdct.c

MPEG-4编解码的实现(包括MPEG4视音频编解码)

    s0 = s0 - s1 + s5;
    s2 = s2 - s3 - s4;
    /* 13 */
    out[13] = s0 * cos_l0[0] + s2 * cos_l0[1];

    /* 4 */
    out[4] = s0 * cos_l0[2] + s2 * cos_l0[3];
}


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;

    if ( init == 0 ) {
        void mdct_init48(void);
	mdct_init48();
	init++;
    }

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

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


	    /* apply filters on the polyphase filterbank outputs */
	    /* bands <= gfp->highpass_band will be zeroed out below */
	    /* bands >= gfp->lowpass_band  will be zeroed out below */
	    if (gfp->filter_type==0) {
	      FLOAT8 amp,freq;
	      for (band=gfp->highpass_band+1;  band < gfp->lowpass_band ; band++) { 
		freq = band/31.0;
		if (gfp->lowpass1 < freq && freq < gfp->lowpass2) {
		  amp = cos((PI/2)*(gfp->lowpass1-freq)/(gfp->lowpass2-gfp->lowpass1));
		  for (k=0; k<18; k++) 
		    sb_sample[ch][1-gr][k][band]*=amp;
		}
		if (gfp->highpass1 < freq && freq < gfp->highpass2) {
		  amp = cos((PI/2)*(gfp->highpass2-freq)/(gfp->highpass2-gfp->highpass1));
		  for (k=0; k<18; k++) 
		    sb_sample[ch][1-gr][k][band]*=amp;
		}
	      }
	    }
	    


	    /*
	     * 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;
#ifdef ALLOW_MIXED
		if (gi->mixed_block_flag && band < 2)
		    type = 0;
#endif
		if (band >= gfp->lowpass_band || band <= gfp->highpass_band) {
		    memset((char *)mdct_enc,0,18*sizeof(FLOAT8));
		}else {
		  if (type == SHORT_TYPE) {
		    for (k = 2; k >= 0; --k) {
		      FLOAT8 w1 = win[SHORT_TYPE][k];
		      work[k] =
			sb_sample[ch][gr][k+6][band] * w1 -
			sb_sample[ch][gr][11-k][band];
		      work[k+3] =
			sb_sample[ch][gr][k+12][band] +
			sb_sample[ch][gr][17-k][band] * w1;
		      
		      work[k+6] =
			sb_sample[ch][gr][k+12][band] * w1 -
			sb_sample[ch][gr][17-k][band];
		      work[k+9] =
			sb_sample[ch][1-gr][k][band] +
			sb_sample[ch][1-gr][5-k][band] * w1;
		      
		      work[k+12] =
			sb_sample[ch][1-gr][k][band] * w1 -
			sb_sample[ch][1-gr][5-k][band];
		      work[k+15] =
			sb_sample[ch][1-gr][k+6][band] +
			sb_sample[ch][1-gr][11-k][band] * w1;
		    }
		    mdct_short(mdct_enc, work);
		  } else {
		    for (k = 8; k >= 0; --k) {
		      work[k] =
			win[type][k  ] * sb_sample[ch][gr][k   ][band]
			- win[type][k+9] * sb_sample[ch][gr][17-k][band];
		      
		      work[9+k] =
			win[type][k+18] * sb_sample[ch][1-gr][k   ][band]
			+ win[type][k+27] * sb_sample[ch][1-gr][17-k][band];
		    }
		    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;
	if (gfp->mode_gr == 1) {
	    memcpy(sb_sample[ch][0], sb_sample[ch][1], 576 * sizeof(FLOAT8));
	}
    }
}



void mdct_init48(void)
{
    int i, k, m;
    FLOAT8 sq;
    FLOAT8 max;

    /* 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];

    sq = 4.0 / NL;
    {
	FLOAT8 *cos_l0 = cos_l;
	static const int d3[] = {1,7,10,16};
	static const int d9[] = {4,13};

	int j = sizeof(all) / sizeof(int) - 1;
	do {
	    m = all[j];
	    for (k = 0; k < NL / 4; k++) {
		*cos_l0++ = sq *
		    cos((PI / (4 * NL)) * (2 * m + 1) * (4 * k + 2 + NL));
	    }
	    for (k = 0; k < NL / 4; k++) {
		*cos_l0++ = sq *
		    cos((PI / (4 * NL)) * (2 * m + 1) * (4 * k + 2 + NL * 3));
	    }
	} while (--j >= 0);

	j = sizeof(d3) / sizeof(int) - 1;
	do {
	    m = d3[j];
	    for (k = 0; k < 3; k++) {
		*cos_l0++ = sq *
		    cos((PI / (4 * NL)) * (2 * m + 1) * (4 * k + 2 + NL));
	    }
	    for (k = 6; k < 9; k++) {
		*cos_l0++ = sq *
		    cos((PI / (4 * NL)) * (2 * m + 1) * (4 * k + 2 + NL));
	    }
	} while (--j >= 0);

	j = sizeof(d9) / sizeof(int) - 1;
	do {
	    m = d9[j];
	    *cos_l0++ = sq *
		cos((PI / (4 * NL)) * (2 * m + 1) * (2 + NL));
	    *cos_l0++ = sq *
		cos((PI / (4 * NL)) * (2 * m + 1) * (4 * 2 + 2 + NL));
	} while (--j >= 0);
    }

    max = enwindow[256 - 8];
    {
	FLOAT8 *wp = enwindow;
	FLOAT8 *wr = enwindow;
	FLOAT8 mmax[32 - 1];

	{
	    FLOAT8 w = *wp++;
	    mmax[15] = w / max;

	    for (k = 0; k < 7; k++) {
		*wr++ = *wp++ / w;
	    }
	}

	for (i = 14; i >= 0; --i) {
	    FLOAT8 w = *wp++;
	    mmax[i] = mmax[30 - i] = w / max;

	    for (k = 0; k < 15; k++) {
		*wr++ = *wp++ / w;
	    }
	}

	{
	    wp++;
	    for (k = 0; k < 7; k++) {
		*wr++ = *wp++ / max;
	    }
	}

	wp = &mm[0][0];
	for (i = 15; i >= 0; --i) {
	    for (k = 1; k < 32; k++) {
		*wp++ = cos((2 * i + 1) * k * PI/64) * mmax[k - 1];
	    }
	}
    }

    /* 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 < 36; i++) {
	win[0][i] *= max / SCALE;
	win[1][i] *= max / SCALE;
	win[3][i] *= max / SCALE;
    }

    /* type 2(short)*/
    sq = 4.0 / NS;
    for (i = 0; i < NS / 4; i++) {
	FLOAT8 w2 = cos(PI/12 * (i + 0.5)) * max / SCALE * sq;
	win[SHORT_TYPE][i] = tan(PI/12 * (i + 0.5));

	for (m = 0; m < NS / 2; m++) {
	    cos_s[m][i] = w2 *
		cos((PI / (4 * NS)) * (2 * m + 1) * (4 * i + 2 + NS));
	    cos_s[m][i + NS / 4] = w2 *
		cos((PI / (4 * NS)) * (2 * m + 1) * (4 * i + 2 + NS * 3));
	}
    }
}