2 (c) Copyright 1998-2001 - 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-06 Andre Piotrowski
25 - speed up: complete new, much faster functions 'create_ana_filter()', 'windowFilterSubband()'
26 - 'WIND_SB_CHANGE_LEVEL' 3 requires a new 'enwindow[]' (see 'tables.c')
30 - bug fix: initWindowFilterSubband() -- Dividing the enwindow-entries is allowed to be done once only!
34 - speed up: WIND_SB_CHANGE_LEVEL 4
38 - bug fix: fixed some typo relevant in case of ORG_BUFFERS == 1
48 /* ======================================================================================== */
50 /* ======================================================================================== */
56 short buffer[2][1152],
58 unsigned int samples_read,
66 for (j = 0; j < samples_read/2; j++)
68 buffer[0][j] = insamp[2*j];
69 buffer[1][j] = insamp[2*j+1];
74 for (j = 0; j < samples_read; j++)
76 buffer[0][j] = insamp[j];
98 FLOAT buffer[2][2048],
100 unsigned int samples_read,
106 /* flush the last two read granules */
107 *buffer_idx = (*buffer_idx + FLUSH) & 2047;
108 idx = (*buffer_idx + DELAY) & 2047;
109 fin = (idx + FLUSH) & 2047;
113 med = (idx + samples_read/2) & 2047;
118 buffer[0][idx] = *insamp++;
119 buffer[1][idx] = *insamp++;
126 buffer[0][idx] = *insamp++;
127 buffer[1][idx] = *insamp++;
133 med = (idx + samples_read) & 2047;
138 buffer[0][idx] = *insamp++;
146 buffer[0][idx] = *insamp++;
180 #define WIND_SB_CHANGE_LEVEL 3
182 #define WIND_SB_CHANGE_LEVEL 4
187 #if WIND_SB_CHANGE_LEVEL==4
189 typedef double filter_matrix[SBLIMIT/4][32];
191 static filter_matrix m;
193 /* ======================================================================================== */
194 /* create_ana_filter */
195 /* ======================================================================================== */
197 Calculates the analysis filterbank coefficients and rounds to the 9th decimal place
198 accuracy of the filterbank tables in the ISO document.
199 The coefficients are stored in #new_filter#
202 Originally was computed
204 filter[i][j] := cos (PI64 * ((2*i+1) * (16-j)))
206 for i = 0..SBLIMIT-1 and j = 0..63 .
209 But for j = 0..15 is 16-j = -(16-(32-j)) which implies
211 (1) filter[i][j] = filter[i][32-j] .
213 (We know that filter[i][16] = cos(0) = 1 , but that is not so interesting.)
215 Furthermore, for j = 33..48 we get
218 = cos (PI64 * (2*i+1) * (j-80))
219 = cos (PI * (2*i+1) + PI64 * (2*i+1) * (16-j))
220 (2) = cos (PI * (2*i+1)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (2*i+1)) = 0
221 = -cos (PI64 * (2*i+1) * (16-j))
224 especially is filter[i][48] = 0 .
227 On the other hand there is
230 = cos (PI64 * (2*(31-i)+1) * (16-j))
231 = cos (PI64 * (64-(2*i+1)) * (16-j))
232 = cos (PI * (16-j) - PI64 * (2*i+1) * (16-j))
233 = cos (PI * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (16-j)) = 0
234 = (-1)^^j * cos (PI64 * (2*i+1) * (16-j))
235 = (-1)^^j * filter[i][j] .
238 There is a third somewhat more complication "symmetry":
241 = cos (PI64 * (2*(15-i)+1) * (16-j))
242 = cos (PI64 * (32-(2*i+1)) * (16-j))
243 = cos (PI/2 * (16-j) - PI64 * (2*i+1) * (16-j))
244 = cos (PI/2 * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) + sin (PI/2 * (16-j)) * sin (PI64 * (2*i+1) * (16-j))
245 = cos (PI/2 * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) + sin (PI/2 * (16-j)) * cos (PI64 * (2*i+1) * (-1)^^i * (16-(j+32)))
247 = (-1)^^( j /2 ) * filter[i][j ] for even j
248 = (-1)^^((j+1)/2 + i) * filter[i][j+32] for odd j with j < 32
249 = (-1)^^((j-1)/2 + i) * filter[i][j-32] for odd j with j >= 32
252 For these reasons we create a filter of the eighth size only and set
254 new_filter[i][j] := filter[i][j] for i = 0..SBLIMIT/4-1 and j = 0..16
255 new_filter[i][j] := filter[i][j+16] for i = 0..SBLIMIT/4-1 and j = 17..31
258 static void create_ana_filter (double new_filter[SBLIMIT/4][32])
263 for (i = 0; i < SBLIMIT/4; i++)
265 for (j = 0; j <= 16; j++)
267 t = 1e9 * cos (PI64 * (double) ((2*i+1) * (16-j)));
272 new_filter[i][j] = t * 1e-9;
274 for (j = 17; j < 32; j++)
276 t = 1e9 * cos (PI64 * (double) ((2*i+1) * j)); /* (16-(j+16)) */
281 new_filter[i][j] = t * 1e-9;
286 /* ======================================================================================== */
287 /* windowFilterSubband */
288 /* ======================================================================================== */
290 Calculates the analysis filter bank coefficients
292 The windowed samples #z# is filtered by the digital filter matrix #m# to produce the
293 subband samples #s#. This done by first selectively picking out values from the windowed
294 samples, and then multiplying them by the filter matrix, producing 32 subband samples.
297 void windowFilterSubband
314 for (k = 0; k <= 16; k++)
316 t = buffer[(buffer_idx+511) & 2047] * *pEnw++;
317 t += buffer[(buffer_idx+447) & 2047] * *pEnw++;
318 t += buffer[(buffer_idx+383) & 2047] * *pEnw++;
319 t += buffer[(buffer_idx+319) & 2047] * *pEnw++;
320 t += buffer[(buffer_idx+255) & 2047] * *pEnw++;
321 t += buffer[(buffer_idx+191) & 2047] * *pEnw++;
322 t += buffer[(buffer_idx+127) & 2047] * *pEnw++;
323 t += buffer[(buffer_idx+ 63) & 2047] * *pEnw++;
329 for (k = 15; k >= 0; k--)
331 t = buffer[(buffer_idx+511) & 2047] * *pEnw++;
332 t += buffer[(buffer_idx+447) & 2047] * *pEnw++;
333 t += buffer[(buffer_idx+383) & 2047] * *pEnw++;
334 t += buffer[(buffer_idx+319) & 2047] * *pEnw++;
335 t += buffer[(buffer_idx+255) & 2047] * *pEnw++;
336 t += buffer[(buffer_idx+191) & 2047] * *pEnw++;
337 t += buffer[(buffer_idx+127) & 2047] * *pEnw++;
338 t += buffer[(buffer_idx+ 63) & 2047] * *pEnw++;
344 for (k = 17; k < 32; k++)
346 t = buffer[(buffer_idx+511) & 2047] * *pEnw++;
347 t += buffer[(buffer_idx+447) & 2047] * *pEnw++;
348 t += buffer[(buffer_idx+383) & 2047] * *pEnw++;
349 t += buffer[(buffer_idx+319) & 2047] * *pEnw++;
350 t += buffer[(buffer_idx+255) & 2047] * *pEnw++;
351 t += buffer[(buffer_idx+191) & 2047] * *pEnw++;
352 t += buffer[(buffer_idx+127) & 2047] * *pEnw++;
353 t += buffer[(buffer_idx+ 63) & 2047] * *pEnw++;
359 /* normally y[48] was computed like the other entries, */
360 /* but this entry is not needed because m[i][48] = 0. */
361 buffer_idx--; /* But we have to increase the pointers. */
364 for (k = 31; k > 16; k--)
366 t = buffer[(buffer_idx+511) & 2047] * *pEnw++;
367 t += buffer[(buffer_idx+447) & 2047] * *pEnw++;
368 t += buffer[(buffer_idx+383) & 2047] * *pEnw++;
369 t += buffer[(buffer_idx+319) & 2047] * *pEnw++;
370 t += buffer[(buffer_idx+255) & 2047] * *pEnw++;
371 t += buffer[(buffer_idx+191) & 2047] * *pEnw++;
372 t += buffer[(buffer_idx+127) & 2047] * *pEnw++;
373 t += buffer[(buffer_idx+ 63) & 2047] * *pEnw++;
379 for (i = 0; i < SBLIMIT/4; i++)
384 for (k = 0; k < 16; k += 4)
387 v += p[k+ 2] * z[k+2];
388 u += p[k+ 1] * z[k+1] + p[k+ 3] * z[k+3];
389 w -= p[k+17] * z[k+1] - p[k+19] * z[k+3];
391 for ( ; k < 32; k += 4)
394 v += p[k+ 2] * z[k+2];
395 u += p[k+ 1] * z[k+1] + p[k+ 3] * z[k+3];
396 w += p[k-15] * z[k+1] - p[k-13] * z[k+3];
400 s[SBLIMIT-1-i] = (t + v) - u;
403 s[SBLIMIT/2-1-i] = (t - v) - w;
404 s[SBLIMIT/2 +i] = (t - v) + w;
408 s[SBLIMIT/2-1-i] = (t - v) + w;
409 s[SBLIMIT/2 +i] = (t - v) - w;
414 #elif WIND_SB_CHANGE_LEVEL==3
416 #define imax (SBLIMIT/4)
417 typedef double filter_matrix[imax][32];
421 static double x[2][512];
422 static filter_matrix m;
424 /* ======================================================================================== */
425 /* create_ana_filter */
426 /* ======================================================================================== */
428 Calculates the analysis filterbank coefficients and rounds to the 9th decimal place
429 accuracy of the filterbank tables in the ISO document.
430 The coefficients are stored in #new_filter#
433 Originally was computed
435 filter[i][j] := cos (PI64 * ((2*i+1) * (16-j)))
437 for i = 0..SBLIMIT-1 and j = 0..63 .
440 But for j = 0..15 is 16-j = -(16-(32-j)) which implies
442 (1) filter[i][j] = filter[i][32-j] .
444 (We know that filter[i][16] = cos(0) = 1 , but that is not so interesting.)
446 Furthermore, for j = 33..48 we get
449 = cos (PI64 * (2*i+1) * (j-80))
450 = cos (PI * (2*i+1) + PI64 * (2*i+1) * (16-j))
451 (2) = cos (PI * (2*i+1)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (2*i+1)) = 0
452 = -cos (PI64 * (2*i+1) * (16-j))
455 especially is filter[i][48] = 0 .
458 On the other hand there is
461 = cos (PI64 * (2*(31-i)+1) * (16-j))
462 = cos (PI64 * (64-(2*i+1)) * (16-j))
463 = cos (PI * (16-j) - PI64 * (2*i+1) * (16-j))
464 = cos (PI * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (16-j)) = 0
465 = (-1)^^j * cos (PI64 * (2*i+1) * (16-j))
466 = (-1)^^j * filter[i][j] .
469 There is a third somewhat more complication "symmetry":
472 = cos (PI64 * (2*(15-i)+1) * (16-j))
473 = cos (PI64 * (32-(2*i+1)) * (16-j))
474 = cos (PI/2 * (16-j) - PI64 * (2*i+1) * (16-j))
475 = cos (PI/2 * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) + sin (PI/2 * (16-j)) * sin (PI64 * (2*i+1) * (16-j))
476 = cos (PI/2 * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) + sin (PI/2 * (16-j)) * cos (PI64 * (2*i+1) * (-1)^^i * (16-(j+32)))
478 = (-1)^^( j /2 ) * filter[i][j ] for even j
479 = (-1)^^((j+1)/2 + i) * filter[i][j+32] for odd j with j < 32
480 = (-1)^^((j-1)/2 + i) * filter[i][j-32] for odd j with j >= 32
483 For these reasons we create a filter of the eighth size only and set
485 new_filter[i][j] := filter[i][j] for i = 0..SBLIMIT/4-1 and j = 0..16
486 new_filter[i][j] := filter[i][j+16] for i = 0..SBLIMIT/4-1 and j = 17..31
489 static void create_ana_filter (double new_filter[imax][32])
494 for (i = 0; i < imax; i++)
496 for (j = 0; j <= 16; j++)
498 t = 1e9 * cos (PI64 * (double) ((2*i+1) * (16-j)));
503 new_filter[i][j] = t * 1e-9;
505 for (j = 17; j < 32; j++)
507 t = 1e9 * cos (PI64 * (double) ((2*i+1) * j)); /* (16-(j+16)) */
512 new_filter[i][j] = t * 1e-9;
517 /* ======================================================================================== */
518 /* windowFilterSubband */
519 /* ======================================================================================== */
521 Calculates the analysis filter bank coefficients
523 The windowed samples #z# is filtered by the digital filter matrix #m# to produce the
524 subband samples #s#. This done by first selectively picking out values from the windowed
525 samples, and then multiplying them by the filter matrix, producing 32 subband samples.
528 void windowFilterSubband
546 px = x[ch] + half[ch];
549 /* replace 32 oldest samples with 32 new samples */
551 for (i = 0; i < 32; i++)
552 px[a + (31-i)*8] = (double) pBuffer[i]/*/SCALE*/;
557 for (k = 0; k <= 16; k++)
559 t = px[(a+0) & 7] * *pEnw++;
560 t += px[(a+1) & 7] * *pEnw++;
561 t += px[(a+2) & 7] * *pEnw++;
562 t += px[(a+3) & 7] * *pEnw++;
563 t += px[(a+4) & 7] * *pEnw++;
564 t += px[(a+5) & 7] * *pEnw++;
565 t += px[(a+6) & 7] * *pEnw++;
566 t += px[(a+7) & 7] * *pEnw++;
572 for (k = 15; k > 0; k--)
574 t = px[(a+0) & 7] * *pEnw++;
575 t += px[(a+1) & 7] * *pEnw++;
576 t += px[(a+2) & 7] * *pEnw++;
577 t += px[(a+3) & 7] * *pEnw++;
578 t += px[(a+4) & 7] * *pEnw++;
579 t += px[(a+5) & 7] * *pEnw++;
580 t += px[(a+6) & 7] * *pEnw++;
581 t += px[(a+7) & 7] * *pEnw++;
587 half[ch] ^= 256; /* toggling 0 or 256 */
590 off[ch] = (a + 7) & 7;
600 t = px[(a+0) & 7] * *pEnw++;
601 t += px[(a+1) & 7] * *pEnw++;
602 t += px[(a+2) & 7] * *pEnw++;
603 t += px[(a+3) & 7] * *pEnw++;
604 t += px[(a+4) & 7] * *pEnw++;
605 t += px[(a+5) & 7] * *pEnw++;
606 t += px[(a+6) & 7] * *pEnw++;
607 t += px[(a+7) & 7] * *pEnw++;
613 for (k = 17; k < 32; k++)
615 t = px[(a+0) & 7] * *pEnw++;
616 t += px[(a+1) & 7] * *pEnw++;
617 t += px[(a+2) & 7] * *pEnw++;
618 t += px[(a+3) & 7] * *pEnw++;
619 t += px[(a+4) & 7] * *pEnw++;
620 t += px[(a+5) & 7] * *pEnw++;
621 t += px[(a+6) & 7] * *pEnw++;
622 t += px[(a+7) & 7] * *pEnw++;
628 /* normally y[48] was computed like the other entries, */
629 /* but this entry is not needed because m[i][48] = 0. */
630 px += 8; /* But we have to increase the pointers. */
633 for (k = 31; k > 16; k--)
635 t = px[(a+0) & 7] * *pEnw++;
636 t += px[(a+1) & 7] * *pEnw++;
637 t += px[(a+2) & 7] * *pEnw++;
638 t += px[(a+3) & 7] * *pEnw++;
639 t += px[(a+4) & 7] * *pEnw++;
640 t += px[(a+5) & 7] * *pEnw++;
641 t += px[(a+6) & 7] * *pEnw++;
642 t += px[(a+7) & 7] * *pEnw++;
648 for (i = 0; i < imax; i++)
653 for (k = 0; k < 16; k += 4)
656 v += pm[k+ 2] * z[k+2];
657 u += pm[k+ 1] * z[k+1] + pm[k+ 3] * z[k+3];
658 w -= pm[k+17] * z[k+1] - pm[k+19] * z[k+3];
660 for ( ; k < 32; k += 4)
663 v += pm[k+ 2] * z[k+2];
664 u += pm[k+ 1] * z[k+1] + pm[k+ 3] * z[k+3];
665 w += pm[k-15] * z[k+1] - pm[k-13] * z[k+3];
669 s[SBLIMIT-1-i] = (t + v) - u;
672 s[SBLIMIT/2-1-i] = (t - v) - w;
673 s[SBLIMIT/2 +i] = (t - v) + w;
677 s[SBLIMIT/2-1-i] = (t - v) + w;
678 s[SBLIMIT/2 +i] = (t - v) - w;
683 #elif WIND_SB_CHANGE_LEVEL==2
685 typedef double filter_matrix[SBLIMIT/2][32];
689 static double x[2][512];
690 static filter_matrix m;
692 /* ======================================================================================== */
693 /* create_ana_filter */
694 /* ======================================================================================== */
696 Calculates the analysis filterbank coefficients and rounds to the 9th decimal place
697 accuracy of the filterbank tables in the ISO document.
698 The coefficients are stored in #new_filter#
701 Originally was computed
703 filter[i][j] := cos (PI64 * ((2*i+1) * (16-j)))
705 for i = 0..SBLIMIT-1 and j = 0..63 .
708 But for j = 0..15 is 16-j = -(16-(32-j)) which implies
710 (1) filter[i][j] = filter[i][32-j] .
712 (We know that filter[i][16] = cos(0) = 1 , but that is not so interesting.)
714 Furthermore, for j = 33..48 we get
717 = cos (PI64 * (2*i+1) * (j-80))
718 = cos (PI * (2*i+1) + PI64 * (2*i+1) * (16-j))
719 (2) = cos (PI * (2*i+1)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (2*i+1)) = 0
720 = -cos (PI64 * (2*i+1) * (16-j))
723 especially is filter[i][48] = 0 .
726 On the other hand there is
729 = cos (PI64 * (2*(31-i)+1) * (16-j))
730 = cos (PI64 * (64-(2*i+1)) * (16-j))
731 = cos (PI * (16-j) - PI64 * (2*i+1) * (16-j))
732 = cos (PI * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (16-j)) = 0
733 = (-1)^^j * cos (PI64 * (2*i+1) * (16-j))
734 = (-1)^^j * filter[i][j] .
737 For these reasons we create a filter of the quarter size only and set
739 new_filter[i][j] := filter[i][j] for i = 0..SBLIMIT/2-1 and j = 0..16
740 new_filter[i][j] := filter[i][j+16] for i = 0..SBLIMIT/2-1 and j = 17..31
743 static void create_ana_filter (double new_filter[SBLIMIT/2][32])
748 for (i = 0; i < SBLIMIT/2; i++)
750 for (j = 0; j < =16; j++)
752 t = 1e9 * cos (PI64 * (double) ((2*i+1) * (16-j)));
757 new_filter[i][j] = t * 1e-9;
759 for (j = 17; j < 32; j++)
761 t = 1e9 * cos (PI64 * (double) ((2*i+1) * j)); /* (16-(j+16)) */
766 new_filter[i][j] = t * 1e-9;
771 /* ======================================================================================== */
772 /* windowFilterSubband */
773 /* ======================================================================================== */
775 Calculates the analysis filter bank coefficients
777 The windowed samples #z# is filtered by the digital filter matrix #m# to produce the
778 subband samples #s#. This done by first selectively picking out values from the windowed
779 samples, and then multiplying them by the filter matrix, producing 32 subband samples.
782 void windowFilterSubband
794 double z[32]; /* y[64]; */
800 px = x[ch] + half[ch];
803 /* replace 32 oldest samples with 32 new samples */
805 for (i = 0; i < 32; i++)
806 px[a + (31-i)*8] = (double) pBuffer[i] /*/SCALE*/;
810 for k = 0..15 we compute
812 z[k] := y[k] + y[32-k]
821 for (k = 0; k <= 16; k++) /* for (j = 0; j < =16; j++) */
823 t = px[(a+0) & 7] * pEnw[64*0];
824 t += px[(a+1) & 7] * pEnw[64*1];
825 t += px[(a+2) & 7] * pEnw[64*2];
826 t += px[(a+3) & 7] * pEnw[64*3];
827 t += px[(a+4) & 7] * pEnw[64*4];
828 t += px[(a+5) & 7] * pEnw[64*5];
829 t += px[(a+6) & 7] * pEnw[64*6];
830 t += px[(a+7) & 7] * pEnw[64*7];
831 z[k] = t; /* y[j] = t; */
837 for (k = 15; k > 0; k--) /* for (j = 17; j < 32; j++) */
839 t = px[(a+0) & 7] * pEnw[64*0];
840 t += px[(a+1) & 7] * pEnw[64*1];
841 t += px[(a+2) & 7] * pEnw[64*2];
842 t += px[(a+3) & 7] * pEnw[64*3];
843 t += px[(a+4) & 7] * pEnw[64*4];
844 t += px[(a+5) & 7] * pEnw[64*5];
845 t += px[(a+6) & 7] * pEnw[64*6];
846 t += px[(a+7) & 7] * pEnw[64*7];
847 z[k] += t; /* y[j] = t; */
853 half[ch] ^= 256; /* toggling 0 or 256 */
856 off[ch] = (a + 7) & 7; /*offset is modulo (HAN_SIZE-1)*/
864 /* k = 0; not needed, actual value of k is 0 */ /* j = 32; */
866 t = px[(a+0) & 7] * pEnw[64*0];
867 t += px[(a+1) & 7] * pEnw[64*1];
868 t += px[(a+2) & 7] * pEnw[64*2];
869 t += px[(a+3) & 7] * pEnw[64*3];
870 t += px[(a+4) & 7] * pEnw[64*4];
871 t += px[(a+5) & 7] * pEnw[64*5];
872 t += px[(a+6) & 7] * pEnw[64*6];
873 t += px[(a+7) & 7] * pEnw[64*7];
874 z[k] += t; /* y[j] = t; */
881 for k = 17..31 we compute
883 z[k] := y[k+16] - y[80-k]
886 for (k = 17; k < 32; k++) /* for (j = 33; j < 47; j++) */
888 t = px[(a+0) & 7] * pEnw[64*0];
889 t += px[(a+1) & 7] * pEnw[64*1];
890 t += px[(a+2) & 7] * pEnw[64*2];
891 t += px[(a+3) & 7] * pEnw[64*3];
892 t += px[(a+4) & 7] * pEnw[64*4];
893 t += px[(a+5) & 7] * pEnw[64*5];
894 t += px[(a+6) & 7] * pEnw[64*6];
895 t += px[(a+7) & 7] * pEnw[64*7];
896 z[k] = t; /* y[j] = t; */
902 /* normally y[48] was computed like the other entries, */
903 /* but this entry is not needed because m[i][48] = 0. */
904 px += 8; /* But we have to increase the pointers. */
907 for (k = 31; k > 16; k--) /* for (j = 49; j < 64; j++) */
909 t = px[(a+0) & 7] * pEnw[64*0];
910 t += px[(a+1) & 7] * pEnw[64*1];
911 t += px[(a+2) & 7] * pEnw[64*2];
912 t += px[(a+3) & 7] * pEnw[64*3];
913 t += px[(a+4) & 7] * pEnw[64*4];
914 t += px[(a+5) & 7] * pEnw[64*5];
915 t += px[(a+6) & 7] * pEnw[64*6];
916 t += px[(a+7) & 7] * pEnw[64*7];
917 z[k] -= t; /* y[j] = t; */
923 pm = m[0]; /* pm = m[0]; */
924 for (i = 0; i < SBLIMIT/2; i++) /* for (i = 0; i < SBLIMIT; i++) */
927 t = u = 0.0; /* t = 0.0; */
928 for (k = 0; k < 32; ) /* for (j = 0; j < 64; j++) */
930 t += *pm++ * z[k++]; /* t += *pm++ * y[j]; */
934 s[ i] = t + u; /* s[i] = t; */
935 s[SBLIMIT-1-i] = t - u;
939 #elif WIND_SB_CHANGE_LEVEL==1
941 typedef double filter_matrix[SBLIMIT][64];
945 static double x[2][512];
946 static filter_matrix m;
948 /* ======================================================================================== */
949 /* create_ana_filter */
950 /* ======================================================================================== */
952 Calculates the analysis filterbank coefficients and rounds to the 9th decimal place
953 accuracy of the filterbank tables in the ISO document.
954 The coefficients are stored in #filter#
957 static void create_ana_filter (filter_matrix filter)
962 for (i = 0; i < SBLIMIT; i++)
964 for (j = 0; j < 64; j++)
966 t = 1e9 * cos (PI64 * (double) ((2*i+1) * (16-j)));
971 filter[i][j] = t * 1e-9;
976 /* ======================================================================================== */
977 /* windowFilterSubband */
978 /* ======================================================================================== */
980 Calculates the analysis filter bank coefficients
982 The windowed samples #y# is filtered by the digital filter matrix #m# to produce the
983 subband samples #s#. This done by first selectively picking out values from the windowed
984 samples, and then multiplying them by the filter matrix, producing 32 subband samples.
987 void windowFilterSubband
1005 px = x[ch] + half[ch];
1008 /* replace 32 oldest samples with 32 new samples */
1010 for (i = 0; i < 32; i++)
1011 px[a + (31-i)*8] = (double) pBuffer[i] /*/SCALE*/;
1016 for (j = 0; j < 32; j++)
1018 t = px[(a+0) & 7] * pEnw[64*0];
1019 t += px[(a+1) & 7] * pEnw[64*1];
1020 t += px[(a+2) & 7] * pEnw[64*2];
1021 t += px[(a+3) & 7] * pEnw[64*3];
1022 t += px[(a+4) & 7] * pEnw[64*4];
1023 t += px[(a+5) & 7] * pEnw[64*5];
1024 t += px[(a+6) & 7] * pEnw[64*6];
1025 t += px[(a+7) & 7] * pEnw[64*7];
1032 half[ch] ^= 256; /* toggling 0 or 256 */
1035 off[ch] = (a + 7) & 7; /*offset is modulo (HAN_SIZE-1)*/
1043 for (j = 32; j < 64; j++)
1045 t = px[(a+0) & 7] * pEnw[64*0];
1046 t += px[(a+1) & 7] * pEnw[64*1];
1047 t += px[(a+2) & 7] * pEnw[64*2];
1048 t += px[(a+3) & 7] * pEnw[64*3];
1049 t += px[(a+4) & 7] * pEnw[64*4];
1050 t += px[(a+5) & 7] * pEnw[64*5];
1051 t += px[(a+6) & 7] * pEnw[64*6];
1052 t += px[(a+7) & 7] * pEnw[64*7];
1060 for (i = 0; i < SBLIMIT; i++)
1063 for (j = 0; j < 64; j++)
1076 /* ======================================================================================== */
1077 /* initWindowFilterSubband */
1078 /* ======================================================================================== */
1080 void initWindowFilterSubband (void)
1082 static int initialized = 0;
1088 create_ana_filter (m);
1090 for (i = 0; i < 512; i++)
1091 enwindow[i] /= SCALE;
1097 half[0] = half[1] = 0;
1098 off[0] = off[1] = 0;
1099 memset (x, 0, sizeof(x));