/* (c) Copyright 1998-2000 - Tord Jansson ====================================== This file is part of the BladeEnc MP3 Encoder, based on ISO's reference code for MPEG Layer 3 compression, and might contain smaller or larger sections that are directly taken from ISO's reference code. All changes to the ISO reference code herein are either copyrighted by Tord Jansson (tord.jansson@swipnet.se) or sublicensed to Tord Jansson by a third party. BladeEnc is free software; you can redistribute this file and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. ------------ Changes ------------ 2000-11-04 Andre Piotrowski - some cosmetics - included "l3psy.h" to use the defined block type constants - negligible speed up of mdct_sub() - speed up of mdct() 2000-12-10 ap - speed up of mdct_sub() and mdct() - MDCT_CHANGE_LEVEL 5 */ /********************************************************************** * ISO MPEG Audio Subgroup Software Simulation Group (1996) * ISO 13818-3 MPEG-2 Audio Encoder - Lower Sampling Frequency Extension * * $Id: mdct.c,v 1.1 1996/02/14 04:04:23 rowlands Exp $ * * $Log: mdct.c,v $ * Revision 1.1 1996/02/14 04:04:23 rowlands * Initial revision * * Received from Mike Coleman **********************************************************************/ #include "common.h" #include #include "l3side.h" #include "l3psy.h" #include "mdct.h" /* This is table B.9: coefficients for aliasing reduction */ static double c[8] = { -0.6, -0.535, -0.33, -0.185, -0.095, -0.041, -0.0142, -0.0037 }; int gr_idx[3] = {0,1,2}; #if 0 static double ca[8]; #endif static double cs[8]; int fInit_mdct_sub; void mdct_sub ( L3SBS (*sb_sample), double (*mdct_freq)[2][576], int stereo, III_side_info_t *l3_side, int mode_gr ) { gr_info *cod_info; #if MDCT_CHANGE_LEVEL < 5 double mdct_in[36]; #endif int gr, ch, band, k, j; double bu, bd; int block_type; double (*mdct_enc_long)[2][32][18] = (double (*)[2][32][18]) mdct_freq; if (!fInit_mdct_sub) { /* prepare the aliasing reduction butterflies */ for (k = 0; k < 8; k++) { double sq = sqrt (1.0 + c[k]*c[k]); #if 0 ca[k] = c[k] / sq; #endif cs[k] = 1.0 / sq; } fInit_mdct_sub++; } for (gr = 0; gr < mode_gr; gr++) { int pre_gr = gr_idx[gr ]; int cur_gr = gr_idx[gr+1]; for (ch = 0; ch < stereo; ch++) { double (*mdct_enc)[18] = mdct_enc_long[gr][ch]; cod_info = &l3_side->gr[gr].ch[ch].tt; block_type = cod_info->block_type; /* Compensate for inversion in the analysis filter */ for (k = 1; k < 18; k++) if (k & 1) for (band = 1; band < 32; band++) if (band & 1) (*sb_sample)[ch][cur_gr][k][band] *= -1.0; /* Perform imdct of 18 previous subband samples + 18 current subband samples */ for (band = 0; band < 32; band++) { #if MDCT_CHANGE_LEVEL < 5 for (k = 0; k < 18; k++) { mdct_in[k] = (*sb_sample)[ch][pre_gr][k][band]; mdct_in[k+18] = (*sb_sample)[ch][cur_gr][k][band]; } if (cod_info->mixed_block_flag && (band < 2)) block_type = NORM_TYPE; /* AND WHEN WILL THE BLOCK_TYPE BE SWITCHED BACK? */ /* &mdct_enc[gr][ch][band][0] */ mdct (mdct_in, mdct_enc/*[gr][ch]*/[band], block_type); #else mdct ((*sb_sample)[ch][pre_gr], (*sb_sample)[ch][cur_gr], band, mdct_enc[band], block_type); #endif } /* Perform aliasing reduction butterfly on long blocks */ if (block_type != SHORT_TYPE) for (band = 0; band < 31; band++) for (k = 0; k < 8; k++) { #if 1 /* This is faster because the calculation can be done more sequential */ bu = (mdct_enc[band ][17-k] + mdct_enc[band+1][ k] * c[k]) * cs[k]; bd = (mdct_enc[band+1][ k] - mdct_enc[band ][17-k] * c[k]) * cs[k]; mdct_enc[band ][17-k] = bu; mdct_enc[band+1][ k] = bd; #else bu = mdct_enc[band ][17-k] * cs[k] + mdct_enc[band+1][ k] * ca[k]; bd = mdct_enc[band+1][ k] * cs[k] - mdct_enc[band ][17-k] * ca[k]; mdct_enc[band ][17-k] = bu; mdct_enc[band+1][ k] = bd; #endif } } } /* Save latest granule's subband samples to be used in the next mdct call */ j = gr_idx[mode_gr]; for (k = mode_gr; k > 0; k--) gr_idx[k] = gr_idx[k-1]; gr_idx[0] = j; } /*-------------------------------------------------------------------*/ /* */ /* Function: Calculation of the MDCT */ /* In the case of long blocks ( block_type 0,1,3 ) there are */ /* 36 coefficents in the time domain and 18 in the frequency */ /* domain. */ /* In the case of short blocks (block_type 2 ) there are 3 */ /* 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 */ /* */ /*-------------------------------------------------------------------*/ int fInit_mdct; #if MDCT_CHANGE_LEVEL == 5 void mdct ( double inA[18][32], double inB[18][32], int band, double *out, int block_type ) { static double cos_l[18][18]; static double cos_s[ 6][ 6]; static double winNorm [18]; static double winShort[ 6]; int k, m, N; double t1, t2, t3, t4, t5, t6, u1, u2; if (!fInit_mdct) { N = 36; for (k = 0; k < N/2; k++) winNorm [k] = sin (PI/36 * (k+0.5)); for (k = 0; k < N/4; k++) for (m = 0; m < N/2; m++) cos_l[k][m] = cos((PI/(2*N)) * (2* k +1+N/2) * (2*m+1)) / (N/4); for ( ; k < N/2; k++) for (m = 0; m < N/2; m++) cos_l[k][m] = cos((PI/(2*N)) * (2*(k+N/4)+1+N/2) * (2*m+1)) / (N/4); N = 12; for (k = 0; k < N/2; k++) winShort[k] = sin (PI/12 * (k+0.5)); for (k = 0; k < N/4; k++) for (m = 0; m < N/2; m++) cos_s[k][m] = cos((PI/(2*N)) * (2*k +1+N/2) * (2*m+1)) / (N/4); for ( ; k < N/2; k++) for (m = 0; m < N/2; m++) cos_s[k][m] = cos((PI/(2*N)) * (2*(k+N/4)+1+N/2) * (2*m+1)) / (N/4); fInit_mdct++; } for (m = 0; m < 18; m++) out[m] = 0.0; switch (block_type) { case NORM_TYPE: for (k = 0; k < 9; k++) { t1 = winNorm [ k] * inA[k][band] - winNorm [17-k] * inA[17-k][band]; t2 = winNorm [17-k] * inB[k][band] + winNorm [ k] * inB[17-k][band]; for (m = 0; m < 18; m++) out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m]; } break; case START_TYPE: for (k = 0; k < 6; k++) { t1 = winNorm [ k] * inA[k][band] - winNorm [17-k] * inA[17-k][band]; t2 = inB[k][band]; for (m = 0; m < 18; m++) out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m]; } for (k = 6; k < 9; k++) { t1 = winNorm [ k] * inA[k][band] - winNorm [17-k] * inA[17-k][band]; t2 = winShort[11-k] * inB[k][band] + winShort[k- 6] * inB[17-k][band]; for (m = 0; m < 18; m++) out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m]; } break; case STOP_TYPE: for (k = 0; k < 6; k++) { t1 = -inA[17-k][band]; t2 = winNorm [17-k] * inB[k][band] + winNorm [ k] * inB[17-k][band]; for (m = 0; m < 18; m++) out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m]; } for (k = 6; k < 9; k++) { t1 = winShort[k- 6] * inA[k][band] - winShort[11-k] * inA[17-k][band]; t2 = winNorm [17-k] * inB[k][band] + winNorm [ k] * inB[17-k][band]; for (m = 0; m < 18; m++) out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m]; } break; case SHORT_TYPE: for (k = 0; k < 3; k++) { u1 = winShort[k]; u2 = winShort[5-k]; t1 = u1 * inA[k+ 6][band] - u2 * inA[11-k][band]; t2 = u2 * inA[k+12][band] + u1 * inA[17-k][band]; t3 = u1 * inA[k+12][band] - u2 * inA[17-k][band]; t4 = u2 * inB[k ][band] + u1 * inB[ 5-k][band]; t5 = u1 * inB[k ][band] - u2 * inB[ 5-k][band]; t6 = u2 * inB[k+ 6][band] + u1 * inB[11-k][band]; for (m = 0; m < 6; m++) { u1 = cos_s[k][m]; u2 = cos_s[k+3][m]; out[3*m ] += t1 * u1 + t2 * u2; out[3*m+1] += t3 * u1 + t4 * u2; out[3*m+2] += t5 * u1 + t6 * u2; } } } } #elif MDCT_CHANGE_LEVEL == 4 /* reduce number of multiplications to nearly the half once more! (cos_x[9+k][m] = -cos_x[9-k][m] and cos_x[35+k][m] = cos_x[35-k][m]) */ void mdct ( double *in, double *out, int block_type ) { static double cos_l[36][18]; static double cos_s[12][ 6]; static double winNorm [36]; static double winShort[12]; static double *winStart = winShort-18; static double *winStop = winShort- 6; int l, k, m, N; double temp; if (!fInit_mdct) { N = 36; for (k = 0; k < N; k++) winNorm [k] = sin (PI/36 * (k+0.5)); for (k = 0; k < N; k++) for (m = 0; m < N/2; m++) cos_l[k][m] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4); N = 12; for (k = 0; k < N; k++) winShort[k] = sin (PI/12 * (k+0.5)); for (k = 0; k < N; k++) for (m = 0; m < N/2; m++) cos_s[k][m] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4); fInit_mdct++; } for (m = 0; m < 18; m++) out[m] = 0.0; switch (block_type) { case NORM_TYPE: for (k = 0; k < 9; k++) {temp = (winNorm [k] * in[k] - winNorm [17-k] * in[17-k]); for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} for (k = 18; k < 27; k++) {temp = (winNorm [k] * in[k] + winNorm [53-k] * in[53-k]); for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} break; case START_TYPE: for (k = 0; k < 9; k++) {temp = (winNorm [k] * in[k] - winNorm [17-k] * in[17-k]); for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} for (k = 18; k < 24; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} for (k = 24; k < 27; k++) {temp = (winStart[k] * in[k] + winStart[53-k] * in[53-k]); for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} break; case STOP_TYPE: for (k = 6; k < 9; k++) {temp = (winStop [k] * in[k] - winStop [17-k] * in[17-k]); for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} for (k = 12; k < 18; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} for (k = 18; k < 27; k++) {temp = (winNorm [k] * in[k] + winNorm [53-k] * in[53-k]); for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} break; case SHORT_TYPE: for (l = 0; l < 3; l++) { for (k = 0; k < 3; k++) {temp = (winShort[k] * in[k+6+l*6] - winShort[ 5-k] * in[ 5-k+6+l*6]); for (m = 0; m < 6; m++) out[3*m+l] += temp * cos_s[k][m];} for (k = 6; k < 9; k++) {temp = (winShort[k] * in[k+6+l*6] + winShort[17-k] * in[17-k+6+l*6]); for (m = 0; m < 6; m++) out[3*m+l] += temp * cos_s[k][m];} } } } #elif MDCT_CHANGE_LEVEL == 3 /* reduce number of multiplications to nearly the half (win_x[k]*in[k] is constant in m-loop)! flip cos_x components for faster access */ void mdct ( double *in, double *out, int block_type ) { static double cos_l[36][18]; static double cos_s[12][ 6]; static double winNorm [36]; static double winShort[12]; static double *winStart = winShort-18; static double *winStop = winShort- 6; int l, k, m, N; double temp; if (!fInit_mdct) { N = 36; for (k = 0; k < N; k++) winNorm [k] = sin (PI/36 * (k+0.5)); for (k = 0; k < N; k++) for (m = 0; m < N/2; m++) cos_l[k][m] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4); N = 12; for (k = 0; k < N; k++) winShort[k] = sin (PI/12 * (k+0.5)); for (k = 0; k < N; k++) for (m = 0; m < N/2; m++) cos_s[k][m] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4); fInit_mdct++; } for (m = 0; m < 18; m++) out[m] = 0.0; switch (block_type) { case NORM_TYPE: for (k = 0; k < 36; k++) {temp = winNorm [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} break; case START_TYPE: for (k = 0; k < 18; k++) {temp = winNorm [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} for (k = 18; k < 24; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} for (k = 24; k < 30; k++) {temp = winStart[k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} break; case STOP_TYPE: for (k = 6; k < 12; k++) {temp = winStop [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} for (k = 12; k < 18; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} for (k = 18; k < 36; k++) {temp = winNorm [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];} break; case SHORT_TYPE: for (l = 0; l < 3; l++) for (k = 0; k < 12; k++) {temp = winShort[k] * in[k+6+l*6]; for (m = 0; m < 6; m++) out[3*m+l] += temp * cos_s[k][m];} } } #elif MDCT_CHANGE_LEVEL == 2 /* avoid calculating of some redundant values; take care for some special values */ void mdct ( double *in, double *out, int block_type ) { static double cos_l[18][36]; static double cos_s[ 6][12]; static double winNorm [36]; static double winShort[12]; static double *winStart = winShort-18; static double *winStop = winShort- 6; int l, k, m, N; double sum; if (!fInit_mdct) { for (k = 0; k < 36; k++) winNorm [k] = sin (PI/36 * (k+0.5)); for (k = 0; k < 12; k++) winShort[k] = sin (PI/12 * (k+0.5)); N = 36; for (m = 0; m < N/2; m++) for (k = 0; k < N; k++) cos_l[m][k] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4); N = 12; for (m = 0; m < N/2; m++) for (k = 0; k < N; k++) cos_s[m][k] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4); fInit_mdct++; } switch (block_type) { case NORM_TYPE: for (m = 0; m < 18; m++) { sum = 0.0; for (k = 0; k < 36; k++) sum += winNorm [k] * in[k] * cos_l[m][k]; out[m] = sum; } break; case START_TYPE: for (m = 0; m < 18; m++) { sum = 0.0; for (k = 0; k < 18; k++) sum += winNorm [k] * in[k] * cos_l[m][k]; for (k = 18; k < 24; k++) sum += in[k] * cos_l[m][k]; for (k = 24; k < 30; k++) sum += winStart[k] * in[k] * cos_l[m][k]; out[m] = sum; } break; case STOP_TYPE: for (m = 0; m < 18; m++) { sum = 0.0; for (k = 6; k < 12; k++) sum += winStop [k] * in[k] * cos_l[m][k]; for (k = 12; k < 18; k++) sum += in[k] * cos_l[m][k]; for (k = 18; k < 36; k++) sum += winNorm [k] * in[k] * cos_l[m][k]; out[m] = sum; } break; case SHORT_TYPE: for (l = 0; l < 3; l++) { for (m = 0; m < 6; m++) { sum = 0.0; for (k = 0; k < 12; k++) sum += winShort[k] * in[k+6+l*6] * cos_s[m][k]; out[3*m+l] = sum; } } } } #elif MDCT_CHANGE_LEVEL == 1 /* reformatted for better overview; use block type constants */ void mdct ( double *in, double *out, int block_type ) { static double cos_l[18][36]; static double cos_s[ 6][12]; static double winNorm [36]; static double winShort[12]; static double winStart[36]; static double winStop [36]; int l, k, m, N; double sum; if (!fInit_mdct) { /* type 0 -- NORM_TYPE */ for (k = 0; k < 36; k++) winNorm [k] = sin (PI/36 * (k+0.5)); /* type 1 -- START_TYPE */ for (k = 0; k < 18; k++) winStart[k] = sin (PI/36 * (k+0.5)); for (k = 18; k < 24; k++) winStart[k] = 1.0; for (k = 24; k < 30; k++) winStart[k] = sin (PI/12 * (k+0.5 - 18)); for (k = 30; k < 36; k++) winStart[k] = 0.0; /* type 3 -- STOP_TYPE */ for (k = 0; k < 6; k++) winStop [k] = 0.0; for (k = 6; k < 12; k++) winStop [k] = sin (PI/12 * (k+0.5 - 6)); for (k = 12; k < 18; k++) winStop [k] = 1.0; for (k = 18; k < 36; k++) winStop [k] = sin (PI/36 * (k+0.5)); /* type 2 -- SHORT_TYPE */ for (k = 0; k < 12; k++) winShort[k] = sin (PI/12 * (k+0.5)); N = 12; for (m = 0; m < N/2; m++) for (k = 0; k < N; k++) cos_s[m][k] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4); N = 36; for (m = 0; m < N/2; m++) for (k = 0; k < N; k++) cos_l[m][k] = cos((PI/(2*N)) * (2*k+1+N/2) * (2*m+1)) / (N/4); fInit_mdct++; } switch (block_type) { case NORM_TYPE: N = 36; for (m = 0; m < N/2; m++) {sum = 0.0; for (k = 0; k < N; k++) sum += winNorm [k] * in[k ] * cos_l[m][k]; out[ m ] = sum;} break; case START_TYPE: N = 36; for (m = 0; m < N/2; m++) {sum = 0.0; for (k = 0; k < N; k++) sum += winStart[k] * in[k ] * cos_l[m][k]; out[ m ] = sum;} break; case STOP_TYPE: N = 36; for (m = 0; m < N/2; m++) {sum = 0.0; for (k = 0; k < N; k++) sum += winStop [k] * in[k ] * cos_l[m][k]; out[ m ] = sum;} break; case SHORT_TYPE: N = 12; for (l = 0; l < 3; l++) {for (m = 0; m < N/2; m++) {sum = 0.0; for (k = 0; k < N; k++) sum += winShort[k] * in[k+6+l*6] * cos_s[m][k]; out[3*m+l] = sum;}} } } #endif /* MDCT_CHANGE_LEVEL */