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 2002/01/10 17:30:01 kramm Exp $
42 * Revision 1.1 2002/01/10 17:30:01 kramm
43 * Version 0.94.1 of the bladeenc mp3 encoder
45 * Revision 1.1 1996/02/14 04:04:23 rowlands
48 * Received from Mike Coleman
49 **********************************************************************/
62 This is table B.9: coefficients for aliasing reduction
64 static double c[8] = { -0.6, -0.535, -0.33, -0.185, -0.095, -0.041, -0.0142, -0.0037 };
70 int gr_idx[3] = {0,1,2};
90 double (*mdct_freq)[2][576],
92 III_side_info_t *l3_side,
97 #if MDCT_CHANGE_LEVEL < 5
100 int gr, ch, band, k, j;
103 double (*mdct_enc_long)[2][32][18] = (double (*)[2][32][18]) mdct_freq;
107 /* prepare the aliasing reduction butterflies */
108 for (k = 0; k < 8; k++)
110 double sq = sqrt (1.0 + c[k]*c[k]);
120 for (gr = 0; gr < mode_gr; gr++)
122 int pre_gr = gr_idx[gr ];
123 int cur_gr = gr_idx[gr+1];
125 for (ch = 0; ch < stereo; ch++)
127 double (*mdct_enc)[18] = mdct_enc_long[gr][ch];
129 cod_info = &l3_side->gr[gr].ch[ch].tt;
130 block_type = cod_info->block_type;
132 /* Compensate for inversion in the analysis filter */
133 for (k = 1; k < 18; k++)
135 for (band = 1; band < 32; band++)
137 (*sb_sample)[ch][cur_gr][k][band] *= -1.0;
139 Perform imdct of 18 previous subband samples
140 + 18 current subband samples
142 for (band = 0; band < 32; band++)
144 #if MDCT_CHANGE_LEVEL < 5
145 for (k = 0; k < 18; k++)
147 mdct_in[k] = (*sb_sample)[ch][pre_gr][k][band];
148 mdct_in[k+18] = (*sb_sample)[ch][cur_gr][k][band];
150 if (cod_info->mixed_block_flag && (band < 2))
151 block_type = NORM_TYPE; /* AND WHEN WILL THE BLOCK_TYPE BE SWITCHED BACK? */
153 /* &mdct_enc[gr][ch][band][0] */
154 mdct (mdct_in, mdct_enc/*[gr][ch]*/[band], block_type);
156 mdct ((*sb_sample)[ch][pre_gr], (*sb_sample)[ch][cur_gr], band, mdct_enc[band], block_type);
161 Perform aliasing reduction butterfly
165 if (block_type != SHORT_TYPE)
166 for (band = 0; band < 31; band++)
167 for (k = 0; k < 8; k++)
169 #if 1 /* This is faster because the calculation can be done more sequential */
170 bu = (mdct_enc[band ][17-k] + mdct_enc[band+1][ k] * c[k]) * cs[k];
171 bd = (mdct_enc[band+1][ k] - mdct_enc[band ][17-k] * c[k]) * cs[k];
172 mdct_enc[band ][17-k] = bu;
173 mdct_enc[band+1][ k] = bd;
175 bu = mdct_enc[band ][17-k] * cs[k] + mdct_enc[band+1][ k] * ca[k];
176 bd = mdct_enc[band+1][ k] * cs[k] - mdct_enc[band ][17-k] * ca[k];
177 mdct_enc[band ][17-k] = bu;
178 mdct_enc[band+1][ k] = bd;
185 Save latest granule's subband samples to be used in
189 for (k = mode_gr; k > 0; k--)
190 gr_idx[k] = gr_idx[k-1];
198 /*-------------------------------------------------------------------*/
200 /* Function: Calculation of the MDCT */
201 /* In the case of long blocks ( block_type 0,1,3 ) there are */
202 /* 36 coefficents in the time domain and 18 in the frequency */
204 /* In the case of short blocks (block_type 2 ) there are 3 */
205 /* transformations with short length. This leads to 12 coefficents */
206 /* in the time and 6 in the frequency domain. In this case the */
207 /* results are stored side by side in the vector out[]. */
211 /*-------------------------------------------------------------------*/
217 #if MDCT_CHANGE_LEVEL == 5
228 static double cos_l[18][18];
229 static double cos_s[ 6][ 6];
231 static double winNorm [18];
232 static double winShort[ 6];
235 double t1, t2, t3, t4, t5, t6, u1, u2;
241 for (k = 0; k < N/2; k++) winNorm [k] = sin (PI/36 * (k+0.5));
242 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);
243 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);
246 for (k = 0; k < N/2; k++) winShort[k] = sin (PI/12 * (k+0.5));
247 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);
248 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);
253 for (m = 0; m < 18; m++) out[m] = 0.0;
257 for (k = 0; k < 9; k++)
259 t1 = winNorm [ k] * inA[k][band] - winNorm [17-k] * inA[17-k][band];
260 t2 = winNorm [17-k] * inB[k][band] + winNorm [ k] * inB[17-k][band];
261 for (m = 0; m < 18; m++)
262 out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
267 for (k = 0; k < 6; k++)
269 t1 = winNorm [ k] * inA[k][band] - winNorm [17-k] * inA[17-k][band];
271 for (m = 0; m < 18; m++)
272 out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
274 for (k = 6; k < 9; k++)
276 t1 = winNorm [ k] * inA[k][band] - winNorm [17-k] * inA[17-k][band];
277 t2 = winShort[11-k] * inB[k][band] + winShort[k- 6] * inB[17-k][band];
278 for (m = 0; m < 18; m++)
279 out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
284 for (k = 0; k < 6; k++)
286 t1 = -inA[17-k][band];
287 t2 = winNorm [17-k] * inB[k][band] + winNorm [ k] * inB[17-k][band];
288 for (m = 0; m < 18; m++)
289 out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
291 for (k = 6; k < 9; k++)
293 t1 = winShort[k- 6] * inA[k][band] - winShort[11-k] * inA[17-k][band];
294 t2 = winNorm [17-k] * inB[k][band] + winNorm [ k] * inB[17-k][band];
295 for (m = 0; m < 18; m++)
296 out[m] += t1 * cos_l[k][m] + t2 * cos_l[k+9][m];
301 for (k = 0; k < 3; k++)
303 u1 = winShort[k]; u2 = winShort[5-k];
304 t1 = u1 * inA[k+ 6][band] - u2 * inA[11-k][band];
305 t2 = u2 * inA[k+12][band] + u1 * inA[17-k][band];
306 t3 = u1 * inA[k+12][band] - u2 * inA[17-k][band];
307 t4 = u2 * inB[k ][band] + u1 * inB[ 5-k][band];
308 t5 = u1 * inB[k ][band] - u2 * inB[ 5-k][band];
309 t6 = u2 * inB[k+ 6][band] + u1 * inB[11-k][band];
310 for (m = 0; m < 6; m++)
312 u1 = cos_s[k][m]; u2 = cos_s[k+3][m];
313 out[3*m ] += t1 * u1 + t2 * u2;
314 out[3*m+1] += t3 * u1 + t4 * u2;
315 out[3*m+2] += t5 * u1 + t6 * u2;
321 #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]) */
330 static double cos_l[36][18];
331 static double cos_s[12][ 6];
333 static double winNorm [36];
334 static double winShort[12];
335 static double *winStart = winShort-18;
336 static double *winStop = winShort- 6;
345 for (k = 0; k < N; k++) winNorm [k] = sin (PI/36 * (k+0.5));
346 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);
349 for (k = 0; k < N; k++) winShort[k] = sin (PI/12 * (k+0.5));
350 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);
355 for (m = 0; m < 18; m++) out[m] = 0.0;
359 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];}
360 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];}
364 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];}
365 for (k = 18; k < 24; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
366 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];}
370 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];}
371 for (k = 12; k < 18; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
372 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];}
376 for (l = 0; l < 3; l++)
378 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];}
379 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];}
384 #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 */
393 static double cos_l[36][18];
394 static double cos_s[12][ 6];
396 static double winNorm [36];
397 static double winShort[12];
398 static double *winStart = winShort-18;
399 static double *winStop = winShort- 6;
408 for (k = 0; k < N; k++) winNorm [k] = sin (PI/36 * (k+0.5));
409 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);
412 for (k = 0; k < N; k++) winShort[k] = sin (PI/12 * (k+0.5));
413 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);
418 for (m = 0; m < 18; m++) out[m] = 0.0;
422 for (k = 0; k < 36; k++) {temp = winNorm [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
426 for (k = 0; k < 18; k++) {temp = winNorm [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
427 for (k = 18; k < 24; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
428 for (k = 24; k < 30; k++) {temp = winStart[k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
432 for (k = 6; k < 12; k++) {temp = winStop [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
433 for (k = 12; k < 18; k++) {temp = in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
434 for (k = 18; k < 36; k++) {temp = winNorm [k] * in[k]; for (m = 0; m < 18; m++) out[m] += temp * cos_l[k][m];}
438 for (l = 0; l < 3; l++)
439 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];}
443 #elif MDCT_CHANGE_LEVEL == 2 /* avoid calculating of some redundant values; take care for some special values */
452 static double cos_l[18][36];
453 static double cos_s[ 6][12];
455 static double winNorm [36];
456 static double winShort[12];
457 static double *winStart = winShort-18;
458 static double *winStop = winShort- 6;
466 for (k = 0; k < 36; k++) winNorm [k] = sin (PI/36 * (k+0.5));
467 for (k = 0; k < 12; k++) winShort[k] = sin (PI/12 * (k+0.5));
469 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);
470 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);
478 for (m = 0; m < 18; m++)
481 for (k = 0; k < 36; k++) sum += winNorm [k] * in[k] * cos_l[m][k];
487 for (m = 0; m < 18; m++)
490 for (k = 0; k < 18; k++) sum += winNorm [k] * in[k] * cos_l[m][k];
491 for (k = 18; k < 24; k++) sum += in[k] * cos_l[m][k];
492 for (k = 24; k < 30; k++) sum += winStart[k] * in[k] * cos_l[m][k];
498 for (m = 0; m < 18; m++)
501 for (k = 6; k < 12; k++) sum += winStop [k] * in[k] * cos_l[m][k];
502 for (k = 12; k < 18; k++) sum += in[k] * cos_l[m][k];
503 for (k = 18; k < 36; k++) sum += winNorm [k] * in[k] * cos_l[m][k];
509 for (l = 0; l < 3; l++)
511 for (m = 0; m < 6; m++)
514 for (k = 0; k < 12; k++) sum += winShort[k] * in[k+6+l*6] * cos_s[m][k];
521 #elif MDCT_CHANGE_LEVEL == 1 /* reformatted for better overview; use block type constants */
530 static double cos_l[18][36];
531 static double cos_s[ 6][12];
533 static double winNorm [36];
534 static double winShort[12];
535 static double winStart[36];
536 static double winStop [36];
544 /* type 0 -- NORM_TYPE */
545 for (k = 0; k < 36; k++) winNorm [k] = sin (PI/36 * (k+0.5));
546 /* type 1 -- START_TYPE */
547 for (k = 0; k < 18; k++) winStart[k] = sin (PI/36 * (k+0.5));
548 for (k = 18; k < 24; k++) winStart[k] = 1.0;
549 for (k = 24; k < 30; k++) winStart[k] = sin (PI/12 * (k+0.5 - 18));
550 for (k = 30; k < 36; k++) winStart[k] = 0.0;
551 /* type 3 -- STOP_TYPE */
552 for (k = 0; k < 6; k++) winStop [k] = 0.0;
553 for (k = 6; k < 12; k++) winStop [k] = sin (PI/12 * (k+0.5 - 6));
554 for (k = 12; k < 18; k++) winStop [k] = 1.0;
555 for (k = 18; k < 36; k++) winStop [k] = sin (PI/36 * (k+0.5));
556 /* type 2 -- SHORT_TYPE */
557 for (k = 0; k < 12; k++) winShort[k] = sin (PI/12 * (k+0.5));
559 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);
560 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);
567 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;
568 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;
569 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;
570 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;}}
574 #endif /* MDCT_CHANGE_LEVEL */