2 (c) Copyright 1998-2000 - Tord Jansson
3 ======================================
5 This file is part of the BladeEnc MP3 Encoder, based on
6 ISO's reference code for MPEG Layer 3 compression, and might
7 contain smaller or larger sections that are directly taken
8 from ISO's reference code.
10 All changes to the ISO reference code herein are either
11 copyrighted by Tord Jansson (tord.jansson@swipnet.se)
12 or sublicensed to Tord Jansson by a third party.
14 BladeEnc is free software; you can redistribute this file
15 and/or modify it under the terms of the GNU Lesser General Public
16 License as published by the Free Software Foundation; either
17 version 2.1 of the License, or (at your option) any later version.
21 ------------ Changes ------------
23 2000-11-04 Andre Piotrowski
26 - included "l3psy.h" to use the defined block type constants
27 - negligible speed up of mdct_sub()
32 - speed up of mdct_sub() and mdct() - MDCT_CHANGE_LEVEL 5
35 /**********************************************************************
36 * ISO MPEG Audio Subgroup Software Simulation Group (1996)
37 * ISO 13818-3 MPEG-2 Audio Encoder - Lower Sampling Frequency Extension
39 * $Id: mdct.c,v 1.1 2004/05/08 12:19:56 kramm Exp $
42 * Revision 1.1 2004/05/08 12:19:56 kramm
43 * Version 0.94.1 of the bladeenc mp3 encoder
45 * Revision 1.1 2002/01/10 17:30:01 kramm
46 * Version 0.94.1 of the bladeenc mp3 encoder
48 * Revision 1.1 1996/02/14 04:04:23 rowlands
51 * Received from Mike Coleman
52 **********************************************************************/
65 This is table B.9: coefficients for aliasing reduction
67 static double c[8] = { -0.6, -0.535, -0.33, -0.185, -0.095, -0.041, -0.0142, -0.0037 };
73 int gr_idx[3] = {0,1,2};
93 double (*mdct_freq)[2][576],
95 III_side_info_t *l3_side,
100 #if MDCT_CHANGE_LEVEL < 5
103 int gr, ch, band, k, j;
106 double (*mdct_enc_long)[2][32][18] = (double (*)[2][32][18]) mdct_freq;
110 /* prepare the aliasing reduction butterflies */
111 for (k = 0; k < 8; k++)
113 double sq = sqrt (1.0 + c[k]*c[k]);
123 for (gr = 0; gr < mode_gr; gr++)
125 int pre_gr = gr_idx[gr ];
126 int cur_gr = gr_idx[gr+1];
128 for (ch = 0; ch < stereo; ch++)
130 double (*mdct_enc)[18] = mdct_enc_long[gr][ch];
132 cod_info = &l3_side->gr[gr].ch[ch].tt;
133 block_type = cod_info->block_type;
135 /* Compensate for inversion in the analysis filter */
136 for (k = 1; k < 18; k++)
138 for (band = 1; band < 32; band++)
140 (*sb_sample)[ch][cur_gr][k][band] *= -1.0;
142 Perform imdct of 18 previous subband samples
143 + 18 current subband samples
145 for (band = 0; band < 32; band++)
147 #if MDCT_CHANGE_LEVEL < 5
148 for (k = 0; k < 18; k++)
150 mdct_in[k] = (*sb_sample)[ch][pre_gr][k][band];
151 mdct_in[k+18] = (*sb_sample)[ch][cur_gr][k][band];
153 if (cod_info->mixed_block_flag && (band < 2))
154 block_type = NORM_TYPE; /* AND WHEN WILL THE BLOCK_TYPE BE SWITCHED BACK? */
156 /* &mdct_enc[gr][ch][band][0] */
157 mdct (mdct_in, mdct_enc/*[gr][ch]*/[band], block_type);
159 mdct ((*sb_sample)[ch][pre_gr], (*sb_sample)[ch][cur_gr], band, mdct_enc[band], block_type);
164 Perform aliasing reduction butterfly
168 if (block_type != SHORT_TYPE)
169 for (band = 0; band < 31; band++)
170 for (k = 0; k < 8; k++)
172 #if 1 /* This is faster because the calculation can be done more sequential */
173 bu = (mdct_enc[band ][17-k] + mdct_enc[band+1][ k] * c[k]) * cs[k];
174 bd = (mdct_enc[band+1][ k] - mdct_enc[band ][17-k] * c[k]) * cs[k];
175 mdct_enc[band ][17-k] = bu;
176 mdct_enc[band+1][ k] = bd;
178 bu = mdct_enc[band ][17-k] * cs[k] + mdct_enc[band+1][ k] * ca[k];
179 bd = mdct_enc[band+1][ k] * cs[k] - mdct_enc[band ][17-k] * ca[k];
180 mdct_enc[band ][17-k] = bu;
181 mdct_enc[band+1][ k] = bd;
188 Save latest granule's subband samples to be used in
192 for (k = mode_gr; k > 0; k--)
193 gr_idx[k] = gr_idx[k-1];
201 /*-------------------------------------------------------------------*/
203 /* Function: Calculation of the MDCT */
204 /* In the case of long blocks ( block_type 0,1,3 ) there are */
205 /* 36 coefficents in the time domain and 18 in the frequency */
207 /* In the case of short blocks (block_type 2 ) there are 3 */
208 /* transformations with short length. This leads to 12 coefficents */
209 /* in the time and 6 in the frequency domain. In this case the */
210 /* results are stored side by side in the vector out[]. */
214 /*-------------------------------------------------------------------*/
220 #if MDCT_CHANGE_LEVEL == 5
231 static double cos_l[18][18];
232 static double cos_s[ 6][ 6];
234 static double winNorm [18];
235 static double winShort[ 6];
238 double t1, t2, t3, t4, t5, t6, u1, u2;
244 for (k = 0; k < N/2; k++) winNorm [k] = sin (PI/36 * (k+0.5));
245 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);
246 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);
249 for (k = 0; k < N/2; k++) winShort[k] = sin (PI/12 * (k+0.5));
250 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);
251 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);
256 for (m = 0; m < 18; m++) out[m] = 0.0;
260 for (k = 0; k < 9; k++)
262 t1 = winNorm [ k] * inA[k][band] - winNorm [17-k] * inA[17-k][band];
263 t2 = winNorm [17-k] * inB[k][band] + winNorm [ k] * inB[17-k][band];
264 for (m = 0; m < 18; m++)
265 out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
270 for (k = 0; k < 6; k++)
272 t1 = winNorm [ k] * inA[k][band] - winNorm [17-k] * inA[17-k][band];
274 for (m = 0; m < 18; m++)
275 out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
277 for (k = 6; k < 9; k++)
279 t1 = winNorm [ k] * inA[k][band] - winNorm [17-k] * inA[17-k][band];
280 t2 = winShort[11-k] * inB[k][band] + winShort[k- 6] * inB[17-k][band];
281 for (m = 0; m < 18; m++)
282 out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
287 for (k = 0; k < 6; k++)
289 t1 = -inA[17-k][band];
290 t2 = winNorm [17-k] * inB[k][band] + winNorm [ k] * inB[17-k][band];
291 for (m = 0; m < 18; m++)
292 out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
294 for (k = 6; k < 9; k++)
296 t1 = winShort[k- 6] * inA[k][band] - winShort[11-k] * inA[17-k][band];
297 t2 = winNorm [17-k] * inB[k][band] + winNorm [ k] * inB[17-k][band];
298 for (m = 0; m < 18; m++)
299 out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
304 for (k = 0; k < 3; k++)
306 u1 = winShort[k]; u2 = winShort[5-k];
307 t1 = u1 * inA[k+ 6][band] - u2 * inA[11-k][band];
308 t2 = u2 * inA[k+12][band] + u1 * inA[17-k][band];
309 t3 = u1 * inA[k+12][band] - u2 * inA[17-k][band];
310 t4 = u2 * inB[k ][band] + u1 * inB[ 5-k][band];
311 t5 = u1 * inB[k ][band] - u2 * inB[ 5-k][band];
312 t6 = u2 * inB[k+ 6][band] + u1 * inB[11-k][band];
313 for (m = 0; m < 6; m++)
315 u1 = cos_s[k][m]; u2 = cos_s[k+3][m];
316 out[3*m ] += t1 * u1 + t2 * u2;
317 out[3*m+1] += t3 * u1 + t4 * u2;
318 out[3*m+2] += t5 * u1 + t6 * u2;
324 #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]) */
333 static double cos_l[36][18];
334 static double cos_s[12][ 6];
336 static double winNorm [36];
337 static double winShort[12];
338 static double *winStart = winShort-18;
339 static double *winStop = winShort- 6;
348 for (k = 0; k < N; k++) winNorm [k] = sin (PI/36 * (k+0.5));
349 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);
352 for (k = 0; k < N; k++) winShort[k] = sin (PI/12 * (k+0.5));
353 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);
358 for (m = 0; m < 18; m++) out[m] = 0.0;
362 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];}
363 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];}
367 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];}
368 for (k = 18; k < 24; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
369 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];}
373 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];}
374 for (k = 12; k < 18; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
375 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];}
379 for (l = 0; l < 3; l++)
381 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];}
382 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];}
387 #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 */
396 static double cos_l[36][18];
397 static double cos_s[12][ 6];
399 static double winNorm [36];
400 static double winShort[12];
401 static double *winStart = winShort-18;
402 static double *winStop = winShort- 6;
411 for (k = 0; k < N; k++) winNorm [k] = sin (PI/36 * (k+0.5));
412 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);
415 for (k = 0; k < N; k++) winShort[k] = sin (PI/12 * (k+0.5));
416 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);
421 for (m = 0; m < 18; m++) out[m] = 0.0;
425 for (k = 0; k < 36; k++) {temp = winNorm [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
429 for (k = 0; k < 18; k++) {temp = winNorm [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
430 for (k = 18; k < 24; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
431 for (k = 24; k < 30; k++) {temp = winStart[k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
435 for (k = 6; k < 12; k++) {temp = winStop [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
436 for (k = 12; k < 18; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
437 for (k = 18; k < 36; k++) {temp = winNorm [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
441 for (l = 0; l < 3; l++)
442 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];}
446 #elif MDCT_CHANGE_LEVEL == 2 /* avoid calculating of some redundant values; take care for some special values */
455 static double cos_l[18][36];
456 static double cos_s[ 6][12];
458 static double winNorm [36];
459 static double winShort[12];
460 static double *winStart = winShort-18;
461 static double *winStop = winShort- 6;
469 for (k = 0; k < 36; k++) winNorm [k] = sin (PI/36 * (k+0.5));
470 for (k = 0; k < 12; k++) winShort[k] = sin (PI/12 * (k+0.5));
472 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);
473 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);
481 for (m = 0; m < 18; m++)
484 for (k = 0; k < 36; k++) sum += winNorm [k] * in[k] * cos_l[m][k];
490 for (m = 0; m < 18; m++)
493 for (k = 0; k < 18; k++) sum += winNorm [k] * in[k] * cos_l[m][k];
494 for (k = 18; k < 24; k++) sum += in[k] * cos_l[m][k];
495 for (k = 24; k < 30; k++) sum += winStart[k] * in[k] * cos_l[m][k];
501 for (m = 0; m < 18; m++)
504 for (k = 6; k < 12; k++) sum += winStop [k] * in[k] * cos_l[m][k];
505 for (k = 12; k < 18; k++) sum += in[k] * cos_l[m][k];
506 for (k = 18; k < 36; k++) sum += winNorm [k] * in[k] * cos_l[m][k];
512 for (l = 0; l < 3; l++)
514 for (m = 0; m < 6; m++)
517 for (k = 0; k < 12; k++) sum += winShort[k] * in[k+6+l*6] * cos_s[m][k];
524 #elif MDCT_CHANGE_LEVEL == 1 /* reformatted for better overview; use block type constants */
533 static double cos_l[18][36];
534 static double cos_s[ 6][12];
536 static double winNorm [36];
537 static double winShort[12];
538 static double winStart[36];
539 static double winStop [36];
547 /* type 0 -- NORM_TYPE */
548 for (k = 0; k < 36; k++) winNorm [k] = sin (PI/36 * (k+0.5));
549 /* type 1 -- START_TYPE */
550 for (k = 0; k < 18; k++) winStart[k] = sin (PI/36 * (k+0.5));
551 for (k = 18; k < 24; k++) winStart[k] = 1.0;
552 for (k = 24; k < 30; k++) winStart[k] = sin (PI/12 * (k+0.5 - 18));
553 for (k = 30; k < 36; k++) winStart[k] = 0.0;
554 /* type 3 -- STOP_TYPE */
555 for (k = 0; k < 6; k++) winStop [k] = 0.0;
556 for (k = 6; k < 12; k++) winStop [k] = sin (PI/12 * (k+0.5 - 6));
557 for (k = 12; k < 18; k++) winStop [k] = 1.0;
558 for (k = 18; k < 36; k++) winStop [k] = sin (PI/36 * (k+0.5));
559 /* type 2 -- SHORT_TYPE */
560 for (k = 0; k < 12; k++) winShort[k] = sin (PI/12 * (k+0.5));
562 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);
563 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);
570 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;
571 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;
572 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;
573 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;}}
577 #endif /* MDCT_CHANGE_LEVEL */