--- /dev/null
+/*
+ (c) Copyright 1998-2001 - 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-06 Andre Piotrowski
+
+ - speed up: complete new, much faster functions 'create_ana_filter()', 'windowFilterSubband()'
+ - 'WIND_SB_CHANGE_LEVEL' 3 requires a new 'enwindow[]' (see 'tables.c')
+
+ 2000-11-22 ap
+
+ - bug fix: initWindowFilterSubband() -- Dividing the enwindow-entries is allowed to be done once only!
+
+ 2000-12-11 ap
+
+ - speed up: WIND_SB_CHANGE_LEVEL 4
+
+ 2001-01-12 ap
+
+ - bug fix: fixed some typo relevant in case of ORG_BUFFERS == 1
+*/
+
+#include "common.h"
+#include "encoder.h"
+
+
+
+
+
+/* ======================================================================================== */
+/* rebuffer_audio */
+/* ======================================================================================== */
+
+#if ORG_BUFFERS
+
+void rebuffer_audio
+(
+ short buffer[2][1152],
+ short *insamp,
+ unsigned int samples_read,
+ int stereo
+)
+{
+ unsigned int j;
+
+ if (stereo == 2)
+ {
+ for (j = 0; j < samples_read/2; j++)
+ {
+ buffer[0][j] = insamp[2*j];
+ buffer[1][j] = insamp[2*j+1];
+ }
+ }
+ else
+ {
+ for (j = 0; j < samples_read; j++)
+ {
+ buffer[0][j] = insamp[j];
+ buffer[1][j] = 0;
+ }
+ }
+
+ for( ; j < 1152; j++)
+ {
+ buffer[0][j] = 0;
+ buffer[1][j] = 0;
+ }
+
+ return;
+}
+
+#else
+
+#define FLUSH 1152
+#define DELAY 768
+
+void rebuffer_audio
+(
+ const short *insamp,
+ FLOAT buffer[2][2048],
+ int *buffer_idx,
+ unsigned int samples_read,
+ int stereo
+)
+{
+ int idx, med, fin;
+
+ /* flush the last two read granules */
+ *buffer_idx = (*buffer_idx + FLUSH) & 2047;
+ idx = (*buffer_idx + DELAY) & 2047;
+ fin = (idx + FLUSH) & 2047;
+
+ if (stereo == 2)
+ {
+ med = (idx + samples_read/2) & 2047;
+ if (idx >= med)
+ {
+ while (idx < 2048)
+ {
+ buffer[0][idx] = *insamp++;
+ buffer[1][idx] = *insamp++;
+ idx++;
+ }
+ idx = 0;
+ }
+ while (idx < med)
+ {
+ buffer[0][idx] = *insamp++;
+ buffer[1][idx] = *insamp++;
+ idx++;
+ }
+ }
+ else
+ {
+ med = (idx + samples_read) & 2047;
+ if (idx >= med)
+ {
+ while (idx < 2048)
+ {
+ buffer[0][idx] = *insamp++;
+ buffer[1][idx] = 0;
+ idx++;
+ }
+ idx = 0;
+ }
+ while (idx < med)
+ {
+ buffer[0][idx] = *insamp++;
+ buffer[1][idx] = 0;
+ idx++;
+ }
+ }
+
+ if (idx != fin)
+ {
+ if (idx > fin)
+ {
+ while (idx < 2048)
+ {
+ buffer[0][idx] = 0;
+ buffer[1][idx] = 0;
+ idx++;
+ }
+ idx = 0;
+ }
+ while (idx < fin)
+ {
+ buffer[0][idx] = 0;
+ buffer[1][idx] = 0;
+ idx++;
+ }
+ }
+}
+
+#endif
+
+
+
+
+
+#if ORG_BUFFERS
+#define WIND_SB_CHANGE_LEVEL 3
+#else
+#define WIND_SB_CHANGE_LEVEL 4
+#endif
+
+
+
+#if WIND_SB_CHANGE_LEVEL==4
+
+ typedef double filter_matrix[SBLIMIT/4][32];
+
+ static filter_matrix m;
+
+ /* ======================================================================================== */
+ /* create_ana_filter */
+ /* ======================================================================================== */
+ /*
+ Calculates the analysis filterbank coefficients and rounds to the 9th decimal place
+ accuracy of the filterbank tables in the ISO document.
+ The coefficients are stored in #new_filter#
+
+
+ Originally was computed
+
+ filter[i][j] := cos (PI64 * ((2*i+1) * (16-j)))
+
+ for i = 0..SBLIMIT-1 and j = 0..63 .
+
+
+ But for j = 0..15 is 16-j = -(16-(32-j)) which implies
+
+ (1) filter[i][j] = filter[i][32-j] .
+
+ (We know that filter[i][16] = cos(0) = 1 , but that is not so interesting.)
+
+ Furthermore, for j = 33..48 we get
+
+ filter[i][96-j]
+ = cos (PI64 * (2*i+1) * (j-80))
+ = cos (PI * (2*i+1) + PI64 * (2*i+1) * (16-j))
+ (2) = cos (PI * (2*i+1)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (2*i+1)) = 0
+ = -cos (PI64 * (2*i+1) * (16-j))
+ = -filter[i][j] ,
+
+ especially is filter[i][48] = 0 .
+
+
+ On the other hand there is
+
+ filter[31-i][j]
+ = cos (PI64 * (2*(31-i)+1) * (16-j))
+ = cos (PI64 * (64-(2*i+1)) * (16-j))
+ = cos (PI * (16-j) - PI64 * (2*i+1) * (16-j))
+ = cos (PI * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (16-j)) = 0
+ = (-1)^^j * cos (PI64 * (2*i+1) * (16-j))
+ = (-1)^^j * filter[i][j] .
+
+
+ There is a third somewhat more complication "symmetry":
+
+ filter[15-i][j]
+ = cos (PI64 * (2*(15-i)+1) * (16-j))
+ = cos (PI64 * (32-(2*i+1)) * (16-j))
+ = cos (PI/2 * (16-j) - PI64 * (2*i+1) * (16-j))
+ = cos (PI/2 * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) + sin (PI/2 * (16-j)) * sin (PI64 * (2*i+1) * (16-j))
+ = 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)))
+
+ = (-1)^^( j /2 ) * filter[i][j ] for even j
+ = (-1)^^((j+1)/2 + i) * filter[i][j+32] for odd j with j < 32
+ = (-1)^^((j-1)/2 + i) * filter[i][j-32] for odd j with j >= 32
+
+
+ For these reasons we create a filter of the eighth size only and set
+
+ new_filter[i][j] := filter[i][j] for i = 0..SBLIMIT/4-1 and j = 0..16
+ new_filter[i][j] := filter[i][j+16] for i = 0..SBLIMIT/4-1 and j = 17..31
+ */
+
+ static void create_ana_filter (double new_filter[SBLIMIT/4][32])
+ {
+ register int i, j;
+ double t;
+
+ for (i = 0; i < SBLIMIT/4; i++)
+ {
+ for (j = 0; j <= 16; j++)
+ {
+ t = 1e9 * cos (PI64 * (double) ((2*i+1) * (16-j)));
+ if (t >= 0)
+ modf (t + 0.5, &t);
+ else
+ modf (t - 0.5, &t);
+ new_filter[i][j] = t * 1e-9;
+ }
+ for (j = 17; j < 32; j++)
+ {
+ t = 1e9 * cos (PI64 * (double) ((2*i+1) * j)); /* (16-(j+16)) */
+ if (t >= 0)
+ modf (t + 0.5, &t);
+ else
+ modf (t - 0.5, &t);
+ new_filter[i][j] = t * 1e-9;
+ }
+ }
+ }
+
+ /* ======================================================================================== */
+ /* windowFilterSubband */
+ /* ======================================================================================== */
+ /*
+ Calculates the analysis filter bank coefficients
+
+ The windowed samples #z# is filtered by the digital filter matrix #m# to produce the
+ subband samples #s#. This done by first selectively picking out values from the windowed
+ samples, and then multiplying them by the filter matrix, producing 32 subband samples.
+ */
+
+ void windowFilterSubband
+ (
+ const FLOAT *buffer,
+ int buffer_idx,
+ double s[SBLIMIT]
+ )
+ {
+ int i, k;
+
+ double t, u, v, w;
+
+ double z[32];
+ double *p;
+ double *pEnw;
+
+ pEnw = enwindow;
+
+ for (k = 0; k <= 16; k++)
+ {
+ t = buffer[(buffer_idx+511) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+447) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+383) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+319) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+255) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+191) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+127) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+ 63) & 2047] * *pEnw++;
+ z[k] = t;
+
+ buffer_idx--;
+ }
+
+ for (k = 15; k >= 0; k--)
+ {
+ t = buffer[(buffer_idx+511) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+447) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+383) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+319) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+255) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+191) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+127) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+ 63) & 2047] * *pEnw++;
+ z[k] += t;
+
+ buffer_idx--;
+ }
+
+ for (k = 17; k < 32; k++)
+ {
+ t = buffer[(buffer_idx+511) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+447) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+383) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+319) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+255) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+191) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+127) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+ 63) & 2047] * *pEnw++;
+ z[k] = t;
+
+ buffer_idx--;
+ }
+
+ /* normally y[48] was computed like the other entries, */
+ /* but this entry is not needed because m[i][48] = 0. */
+ buffer_idx--; /* But we have to increase the pointers. */
+ pEnw += 8;
+
+ for (k = 31; k > 16; k--)
+ {
+ t = buffer[(buffer_idx+511) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+447) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+383) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+319) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+255) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+191) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+127) & 2047] * *pEnw++;
+ t += buffer[(buffer_idx+ 63) & 2047] * *pEnw++;
+ z[k] -= t;
+
+ buffer_idx--;
+ }
+
+ for (i = 0; i < SBLIMIT/4; i++)
+ {
+ p = m[i];
+
+ t = u = v = w = 0.0;
+ for (k = 0; k < 16; k += 4)
+ {
+ t += p[k ] * z[k ];
+ v += p[k+ 2] * z[k+2];
+ u += p[k+ 1] * z[k+1] + p[k+ 3] * z[k+3];
+ w -= p[k+17] * z[k+1] - p[k+19] * z[k+3];
+ }
+ for ( ; k < 32; k += 4)
+ {
+ t += p[k ] * z[k ];
+ v += p[k+ 2] * z[k+2];
+ u += p[k+ 1] * z[k+1] + p[k+ 3] * z[k+3];
+ w += p[k-15] * z[k+1] - p[k-13] * z[k+3];
+ }
+
+ s[ i] = (t + v) + u;
+ s[SBLIMIT-1-i] = (t + v) - u;
+ if (i & 1)
+ {
+ s[SBLIMIT/2-1-i] = (t - v) - w;
+ s[SBLIMIT/2 +i] = (t - v) + w;
+ }
+ else
+ {
+ s[SBLIMIT/2-1-i] = (t - v) + w;
+ s[SBLIMIT/2 +i] = (t - v) - w;
+ }
+ }
+ }
+
+#elif WIND_SB_CHANGE_LEVEL==3
+
+ #define imax (SBLIMIT/4)
+ typedef double filter_matrix[imax][32];
+
+ static int half[2];
+ static int off[2];
+ static double x[2][512];
+ static filter_matrix m;
+
+ /* ======================================================================================== */
+ /* create_ana_filter */
+ /* ======================================================================================== */
+ /*
+ Calculates the analysis filterbank coefficients and rounds to the 9th decimal place
+ accuracy of the filterbank tables in the ISO document.
+ The coefficients are stored in #new_filter#
+
+
+ Originally was computed
+
+ filter[i][j] := cos (PI64 * ((2*i+1) * (16-j)))
+
+ for i = 0..SBLIMIT-1 and j = 0..63 .
+
+
+ But for j = 0..15 is 16-j = -(16-(32-j)) which implies
+
+ (1) filter[i][j] = filter[i][32-j] .
+
+ (We know that filter[i][16] = cos(0) = 1 , but that is not so interesting.)
+
+ Furthermore, for j = 33..48 we get
+
+ filter[i][96-j]
+ = cos (PI64 * (2*i+1) * (j-80))
+ = cos (PI * (2*i+1) + PI64 * (2*i+1) * (16-j))
+ (2) = cos (PI * (2*i+1)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (2*i+1)) = 0
+ = -cos (PI64 * (2*i+1) * (16-j))
+ = -filter[i][j] ,
+
+ especially is filter[i][48] = 0 .
+
+
+ On the other hand there is
+
+ filter[31-i][j]
+ = cos (PI64 * (2*(31-i)+1) * (16-j))
+ = cos (PI64 * (64-(2*i+1)) * (16-j))
+ = cos (PI * (16-j) - PI64 * (2*i+1) * (16-j))
+ = cos (PI * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (16-j)) = 0
+ = (-1)^^j * cos (PI64 * (2*i+1) * (16-j))
+ = (-1)^^j * filter[i][j] .
+
+
+ There is a third somewhat more complication "symmetry":
+
+ filter[15-i][j]
+ = cos (PI64 * (2*(15-i)+1) * (16-j))
+ = cos (PI64 * (32-(2*i+1)) * (16-j))
+ = cos (PI/2 * (16-j) - PI64 * (2*i+1) * (16-j))
+ = cos (PI/2 * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) + sin (PI/2 * (16-j)) * sin (PI64 * (2*i+1) * (16-j))
+ = 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)))
+
+ = (-1)^^( j /2 ) * filter[i][j ] for even j
+ = (-1)^^((j+1)/2 + i) * filter[i][j+32] for odd j with j < 32
+ = (-1)^^((j-1)/2 + i) * filter[i][j-32] for odd j with j >= 32
+
+
+ For these reasons we create a filter of the eighth size only and set
+
+ new_filter[i][j] := filter[i][j] for i = 0..SBLIMIT/4-1 and j = 0..16
+ new_filter[i][j] := filter[i][j+16] for i = 0..SBLIMIT/4-1 and j = 17..31
+ */
+
+ static void create_ana_filter (double new_filter[imax][32])
+ {
+ register int i, j;
+ double t;
+
+ for (i = 0; i < imax; i++)
+ {
+ for (j = 0; j <= 16; j++)
+ {
+ t = 1e9 * cos (PI64 * (double) ((2*i+1) * (16-j)));
+ if (t >= 0)
+ modf (t + 0.5, &t);
+ else
+ modf (t - 0.5, &t);
+ new_filter[i][j] = t * 1e-9;
+ }
+ for (j = 17; j < 32; j++)
+ {
+ t = 1e9 * cos (PI64 * (double) ((2*i+1) * j)); /* (16-(j+16)) */
+ if (t >= 0)
+ modf (t + 0.5, &t);
+ else
+ modf (t - 0.5, &t);
+ new_filter[i][j] = t * 1e-9;
+ }
+ }
+ }
+
+ /* ======================================================================================== */
+ /* windowFilterSubband */
+ /* ======================================================================================== */
+ /*
+ Calculates the analysis filter bank coefficients
+
+ The windowed samples #z# is filtered by the digital filter matrix #m# to produce the
+ subband samples #s#. This done by first selectively picking out values from the windowed
+ samples, and then multiplying them by the filter matrix, producing 32 subband samples.
+ */
+
+ void windowFilterSubband
+ (
+ short *pBuffer,
+ int ch,
+ double s[SBLIMIT]
+ )
+ {
+ int i, k;
+ int a;
+
+ double t, u, v, w;
+
+ double z[32];
+ double *pm;
+ double *px;
+ double *pEnw;
+
+
+ px = x[ch] + half[ch];
+ a = off[ch];
+
+ /* replace 32 oldest samples with 32 new samples */
+
+ for (i = 0; i < 32; i++)
+ px[a + (31-i)*8] = (double) pBuffer[i]/*/SCALE*/;
+
+
+ pEnw = enwindow;
+
+ for (k = 0; k <= 16; k++)
+ {
+ t = px[(a+0) & 7] * *pEnw++;
+ t += px[(a+1) & 7] * *pEnw++;
+ t += px[(a+2) & 7] * *pEnw++;
+ t += px[(a+3) & 7] * *pEnw++;
+ t += px[(a+4) & 7] * *pEnw++;
+ t += px[(a+5) & 7] * *pEnw++;
+ t += px[(a+6) & 7] * *pEnw++;
+ t += px[(a+7) & 7] * *pEnw++;
+ z[k] = t;
+
+ px += 8;
+ }
+
+ for (k = 15; k > 0; k--)
+ {
+ t = px[(a+0) & 7] * *pEnw++;
+ t += px[(a+1) & 7] * *pEnw++;
+ t += px[(a+2) & 7] * *pEnw++;
+ t += px[(a+3) & 7] * *pEnw++;
+ t += px[(a+4) & 7] * *pEnw++;
+ t += px[(a+5) & 7] * *pEnw++;
+ t += px[(a+6) & 7] * *pEnw++;
+ t += px[(a+7) & 7] * *pEnw++;
+ z[k] += t;
+
+ px += 8;
+ }
+
+ half[ch] ^= 256; /* toggling 0 or 256 */
+ if (half[ch])
+ {
+ off[ch] = (a + 7) & 7;
+ }
+ else
+ {
+ px = x[ch];
+ a = (a + 1) & 7;
+ }
+
+ /* k = 0; */
+ {
+ t = px[(a+0) & 7] * *pEnw++;
+ t += px[(a+1) & 7] * *pEnw++;
+ t += px[(a+2) & 7] * *pEnw++;
+ t += px[(a+3) & 7] * *pEnw++;
+ t += px[(a+4) & 7] * *pEnw++;
+ t += px[(a+5) & 7] * *pEnw++;
+ t += px[(a+6) & 7] * *pEnw++;
+ t += px[(a+7) & 7] * *pEnw++;
+ z[k] += t;
+
+ px += 8;
+ }
+
+ for (k = 17; k < 32; k++)
+ {
+ t = px[(a+0) & 7] * *pEnw++;
+ t += px[(a+1) & 7] * *pEnw++;
+ t += px[(a+2) & 7] * *pEnw++;
+ t += px[(a+3) & 7] * *pEnw++;
+ t += px[(a+4) & 7] * *pEnw++;
+ t += px[(a+5) & 7] * *pEnw++;
+ t += px[(a+6) & 7] * *pEnw++;
+ t += px[(a+7) & 7] * *pEnw++;
+ z[k] = t;
+
+ px += 8;
+ }
+
+ /* normally y[48] was computed like the other entries, */
+ /* but this entry is not needed because m[i][48] = 0. */
+ px += 8; /* But we have to increase the pointers. */
+ pEnw += 8;
+
+ for (k = 31; k > 16; k--)
+ {
+ t = px[(a+0) & 7] * *pEnw++;
+ t += px[(a+1) & 7] * *pEnw++;
+ t += px[(a+2) & 7] * *pEnw++;
+ t += px[(a+3) & 7] * *pEnw++;
+ t += px[(a+4) & 7] * *pEnw++;
+ t += px[(a+5) & 7] * *pEnw++;
+ t += px[(a+6) & 7] * *pEnw++;
+ t += px[(a+7) & 7] * *pEnw++;
+ z[k] -= t;
+
+ px += 8;
+ }
+
+ for (i = 0; i < imax; i++)
+ {
+ pm = m[i];
+
+ t = u = v = w = 0.0;
+ for (k = 0; k < 16; k += 4)
+ {
+ t += pm[k ] * z[k ];
+ v += pm[k+ 2] * z[k+2];
+ u += pm[k+ 1] * z[k+1] + pm[k+ 3] * z[k+3];
+ w -= pm[k+17] * z[k+1] - pm[k+19] * z[k+3];
+ }
+ for ( ; k < 32; k += 4)
+ {
+ t += pm[k ] * z[k ];
+ v += pm[k+ 2] * z[k+2];
+ u += pm[k+ 1] * z[k+1] + pm[k+ 3] * z[k+3];
+ w += pm[k-15] * z[k+1] - pm[k-13] * z[k+3];
+ }
+
+ s[ i] = (t + v) + u;
+ s[SBLIMIT-1-i] = (t + v) - u;
+ if (i & 1)
+ {
+ s[SBLIMIT/2-1-i] = (t - v) - w;
+ s[SBLIMIT/2 +i] = (t - v) + w;
+ }
+ else
+ {
+ s[SBLIMIT/2-1-i] = (t - v) + w;
+ s[SBLIMIT/2 +i] = (t - v) - w;
+ }
+ }
+ }
+
+#elif WIND_SB_CHANGE_LEVEL==2
+
+ typedef double filter_matrix[SBLIMIT/2][32];
+
+ static int off[2];
+ static int half[2];
+ static double x[2][512];
+ static filter_matrix m;
+
+ /* ======================================================================================== */
+ /* create_ana_filter */
+ /* ======================================================================================== */
+ /*
+ Calculates the analysis filterbank coefficients and rounds to the 9th decimal place
+ accuracy of the filterbank tables in the ISO document.
+ The coefficients are stored in #new_filter#
+
+
+ Originally was computed
+
+ filter[i][j] := cos (PI64 * ((2*i+1) * (16-j)))
+
+ for i = 0..SBLIMIT-1 and j = 0..63 .
+
+
+ But for j = 0..15 is 16-j = -(16-(32-j)) which implies
+
+ (1) filter[i][j] = filter[i][32-j] .
+
+ (We know that filter[i][16] = cos(0) = 1 , but that is not so interesting.)
+
+ Furthermore, for j = 33..48 we get
+
+ filter[i][96-j]
+ = cos (PI64 * (2*i+1) * (j-80))
+ = cos (PI * (2*i+1) + PI64 * (2*i+1) * (16-j))
+ (2) = cos (PI * (2*i+1)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (2*i+1)) = 0
+ = -cos (PI64 * (2*i+1) * (16-j))
+ = -filter[i][j] ,
+
+ especially is filter[i][48] = 0 .
+
+
+ On the other hand there is
+
+ filter[31-i][j]
+ = cos (PI64 * (2*(31-i)+1) * (16-j))
+ = cos (PI64 * (64-(2*i+1)) * (16-j))
+ = cos (PI * (16-j) - PI64 * (2*i+1) * (16-j))
+ = cos (PI * (16-j)) * cos (PI64 * (2*i+1) * (16-j)) // sin (PI * (16-j)) = 0
+ = (-1)^^j * cos (PI64 * (2*i+1) * (16-j))
+ = (-1)^^j * filter[i][j] .
+
+
+ For these reasons we create a filter of the quarter size only and set
+
+ new_filter[i][j] := filter[i][j] for i = 0..SBLIMIT/2-1 and j = 0..16
+ new_filter[i][j] := filter[i][j+16] for i = 0..SBLIMIT/2-1 and j = 17..31
+ */
+
+ static void create_ana_filter (double new_filter[SBLIMIT/2][32])
+ {
+ register int i, j;
+ double t;
+
+ for (i = 0; i < SBLIMIT/2; i++)
+ {
+ for (j = 0; j < =16; j++)
+ {
+ t = 1e9 * cos (PI64 * (double) ((2*i+1) * (16-j)));
+ if (t >= 0)
+ modf (t + 0.5, &t);
+ else
+ modf (t - 0.5, &t);
+ new_filter[i][j] = t * 1e-9;
+ }
+ for (j = 17; j < 32; j++)
+ {
+ t = 1e9 * cos (PI64 * (double) ((2*i+1) * j)); /* (16-(j+16)) */
+ if (t >= 0)
+ modf (t + 0.5, &t);
+ else
+ modf (t - 0.5, &t);
+ new_filter[i][j] = t * 1e-9;
+ }
+ }
+ }
+
+ /* ======================================================================================== */
+ /* windowFilterSubband */
+ /* ======================================================================================== */
+ /*
+ Calculates the analysis filter bank coefficients
+
+ The windowed samples #z# is filtered by the digital filter matrix #m# to produce the
+ subband samples #s#. This done by first selectively picking out values from the windowed
+ samples, and then multiplying them by the filter matrix, producing 32 subband samples.
+ */
+
+ void windowFilterSubband
+ (
+ short *pBuffer,
+ int ch,
+ double s[SBLIMIT]
+ )
+ {
+ int i, k;
+ int a;
+
+ double t, u;
+
+ double z[32]; /* y[64]; */
+ double *pm;
+ double *px;
+ double *pEnw;
+
+
+ px = x[ch] + half[ch];
+ a = off[ch];
+
+ /* replace 32 oldest samples with 32 new samples */
+
+ for (i = 0; i < 32; i++)
+ px[a + (31-i)*8] = (double) pBuffer[i] /*/SCALE*/;
+
+
+ /*
+ for k = 0..15 we compute
+
+ z[k] := y[k] + y[32-k]
+
+ and we set
+
+ z[16| := y[16] .
+ */
+
+ pEnw = enwindow;
+
+ for (k = 0; k <= 16; k++) /* for (j = 0; j < =16; j++) */
+ {
+ t = px[(a+0) & 7] * pEnw[64*0];
+ t += px[(a+1) & 7] * pEnw[64*1];
+ t += px[(a+2) & 7] * pEnw[64*2];
+ t += px[(a+3) & 7] * pEnw[64*3];
+ t += px[(a+4) & 7] * pEnw[64*4];
+ t += px[(a+5) & 7] * pEnw[64*5];
+ t += px[(a+6) & 7] * pEnw[64*6];
+ t += px[(a+7) & 7] * pEnw[64*7];
+ z[k] = t; /* y[j] = t; */
+
+ px += 8;
+ pEnw++;
+ }
+
+ for (k = 15; k > 0; k--) /* for (j = 17; j < 32; j++) */
+ {
+ t = px[(a+0) & 7] * pEnw[64*0];
+ t += px[(a+1) & 7] * pEnw[64*1];
+ t += px[(a+2) & 7] * pEnw[64*2];
+ t += px[(a+3) & 7] * pEnw[64*3];
+ t += px[(a+4) & 7] * pEnw[64*4];
+ t += px[(a+5) & 7] * pEnw[64*5];
+ t += px[(a+6) & 7] * pEnw[64*6];
+ t += px[(a+7) & 7] * pEnw[64*7];
+ z[k] += t; /* y[j] = t; */
+
+ px += 8;
+ pEnw++;
+ }
+
+ half[ch] ^= 256; /* toggling 0 or 256 */
+ if (half[ch])
+ {
+ off[ch] = (a + 7) & 7; /*offset is modulo (HAN_SIZE-1)*/
+ }
+ else
+ {
+ px = x[ch];
+ a = (a + 1) & 7;
+ }
+
+ /* k = 0; not needed, actual value of k is 0 */ /* j = 32; */
+ {
+ t = px[(a+0) & 7] * pEnw[64*0];
+ t += px[(a+1) & 7] * pEnw[64*1];
+ t += px[(a+2) & 7] * pEnw[64*2];
+ t += px[(a+3) & 7] * pEnw[64*3];
+ t += px[(a+4) & 7] * pEnw[64*4];
+ t += px[(a+5) & 7] * pEnw[64*5];
+ t += px[(a+6) & 7] * pEnw[64*6];
+ t += px[(a+7) & 7] * pEnw[64*7];
+ z[k] += t; /* y[j] = t; */
+
+ px += 8;
+ pEnw++;
+ }
+
+ /*
+ for k = 17..31 we compute
+
+ z[k] := y[k+16] - y[80-k]
+ */
+
+ for (k = 17; k < 32; k++) /* for (j = 33; j < 47; j++) */
+ {
+ t = px[(a+0) & 7] * pEnw[64*0];
+ t += px[(a+1) & 7] * pEnw[64*1];
+ t += px[(a+2) & 7] * pEnw[64*2];
+ t += px[(a+3) & 7] * pEnw[64*3];
+ t += px[(a+4) & 7] * pEnw[64*4];
+ t += px[(a+5) & 7] * pEnw[64*5];
+ t += px[(a+6) & 7] * pEnw[64*6];
+ t += px[(a+7) & 7] * pEnw[64*7];
+ z[k] = t; /* y[j] = t; */
+
+ px += 8;
+ pEnw++;
+ }
+
+ /* normally y[48] was computed like the other entries, */
+ /* but this entry is not needed because m[i][48] = 0. */
+ px += 8; /* But we have to increase the pointers. */
+ pEnw++;
+
+ for (k = 31; k > 16; k--) /* for (j = 49; j < 64; j++) */
+ {
+ t = px[(a+0) & 7] * pEnw[64*0];
+ t += px[(a+1) & 7] * pEnw[64*1];
+ t += px[(a+2) & 7] * pEnw[64*2];
+ t += px[(a+3) & 7] * pEnw[64*3];
+ t += px[(a+4) & 7] * pEnw[64*4];
+ t += px[(a+5) & 7] * pEnw[64*5];
+ t += px[(a+6) & 7] * pEnw[64*6];
+ t += px[(a+7) & 7] * pEnw[64*7];
+ z[k] -= t; /* y[j] = t; */
+
+ px += 8;
+ pEnw++;
+ }
+
+ pm = m[0]; /* pm = m[0]; */
+ for (i = 0; i < SBLIMIT/2; i++) /* for (i = 0; i < SBLIMIT; i++) */
+ {
+
+ t = u = 0.0; /* t = 0.0; */
+ for (k = 0; k < 32; ) /* for (j = 0; j < 64; j++) */
+ {
+ t += *pm++ * z[k++]; /* t += *pm++ * y[j]; */
+ u += *pm++ * z[k++];
+ }
+
+ s[ i] = t + u; /* s[i] = t; */
+ s[SBLIMIT-1-i] = t - u;
+ }
+ }
+
+#elif WIND_SB_CHANGE_LEVEL==1
+
+ typedef double filter_matrix[SBLIMIT][64];
+
+ static int off[2];
+ static int half[2];
+ static double x[2][512];
+ static filter_matrix m;
+
+ /* ======================================================================================== */
+ /* create_ana_filter */
+ /* ======================================================================================== */
+ /*
+ Calculates the analysis filterbank coefficients and rounds to the 9th decimal place
+ accuracy of the filterbank tables in the ISO document.
+ The coefficients are stored in #filter#
+ */
+
+ static void create_ana_filter (filter_matrix filter)
+ {
+ register int i, j;
+ double t;
+
+ for (i = 0; i < SBLIMIT; i++)
+ {
+ for (j = 0; j < 64; j++)
+ {
+ t = 1e9 * cos (PI64 * (double) ((2*i+1) * (16-j)));
+ if (t >= 0)
+ modf (t + 0.5, &t);
+ else
+ modf (t - 0.5, &t);
+ filter[i][j] = t * 1e-9;
+ }
+ }
+ }
+
+ /* ======================================================================================== */
+ /* windowFilterSubband */
+ /* ======================================================================================== */
+ /*
+ Calculates the analysis filter bank coefficients
+
+ The windowed samples #y# is filtered by the digital filter matrix #m# to produce the
+ subband samples #s#. This done by first selectively picking out values from the windowed
+ samples, and then multiplying them by the filter matrix, producing 32 subband samples.
+ */
+
+ void windowFilterSubband
+ (
+ short *pBuffer,
+ int ch,
+ double s[SBLIMIT]
+ )
+ {
+ int i, j;
+ int a;
+
+ double t;
+
+ double y[64];
+ double *pm;
+ double *px;
+ double *pEnw;
+
+
+ px = x[ch] + half[ch];
+ a = off[ch];
+
+ /* replace 32 oldest samples with 32 new samples */
+
+ for (i = 0; i < 32; i++)
+ px[a + (31-i)*8] = (double) pBuffer[i] /*/SCALE*/;
+
+
+ pEnw = enwindow;
+
+ for (j = 0; j < 32; j++)
+ {
+ t = px[(a+0) & 7] * pEnw[64*0];
+ t += px[(a+1) & 7] * pEnw[64*1];
+ t += px[(a+2) & 7] * pEnw[64*2];
+ t += px[(a+3) & 7] * pEnw[64*3];
+ t += px[(a+4) & 7] * pEnw[64*4];
+ t += px[(a+5) & 7] * pEnw[64*5];
+ t += px[(a+6) & 7] * pEnw[64*6];
+ t += px[(a+7) & 7] * pEnw[64*7];
+ y[j] = t;
+
+ px += 8;
+ pEnw++;
+ }
+
+ half[ch] ^= 256; /* toggling 0 or 256 */
+ if (half[ch])
+ {
+ off[ch] = (a + 7) & 7; /*offset is modulo (HAN_SIZE-1)*/
+ }
+ else
+ {
+ px = x[ch];
+ a = (a + 1) & 7;
+ }
+
+ for (j = 32; j < 64; j++)
+ {
+ t = px[(a+0) & 7] * pEnw[64*0];
+ t += px[(a+1) & 7] * pEnw[64*1];
+ t += px[(a+2) & 7] * pEnw[64*2];
+ t += px[(a+3) & 7] * pEnw[64*3];
+ t += px[(a+4) & 7] * pEnw[64*4];
+ t += px[(a+5) & 7] * pEnw[64*5];
+ t += px[(a+6) & 7] * pEnw[64*6];
+ t += px[(a+7) & 7] * pEnw[64*7];
+ y[j] = t;
+
+ px += 8;
+ pEnw++;
+ }
+
+ pm = m[0];
+ for (i = 0; i < SBLIMIT; i++)
+ {
+ t = 0.0;
+ for (j = 0; j < 64; j++)
+ t += *pm++ * y[j];
+
+ s[i] = t;
+ }
+ }
+
+#endif
+
+
+
+
+
+/* ======================================================================================== */
+/* initWindowFilterSubband */
+/* ======================================================================================== */
+
+void initWindowFilterSubband (void)
+{
+ static int initialized = 0;
+
+ int i;
+
+ if (!initialized)
+ {
+ create_ana_filter (m);
+
+ for (i = 0; i < 512; i++)
+ enwindow[i] /= SCALE;
+
+ initialized = 1;
+ }
+
+#if ORG_BUFFERS
+ half[0] = half[1] = 0;
+ off[0] = off[1] = 0;
+ memset (x, 0, sizeof(x));
+#endif
+}
+
+
+
+
+
git.asbjorn.biz / swftools.git / blobdiff