summaryrefslogtreecommitdiffstats
path: root/libtoolame-dab/pds_subband.c
diff options
context:
space:
mode:
Diffstat (limited to 'libtoolame-dab/pds_subband.c')
-rw-r--r--libtoolame-dab/pds_subband.c361
1 files changed, 361 insertions, 0 deletions
diff --git a/libtoolame-dab/pds_subband.c b/libtoolame-dab/pds_subband.c
new file mode 100644
index 0000000..04938c4
--- /dev/null
+++ b/libtoolame-dab/pds_subband.c
@@ -0,0 +1,361 @@
+#include "common.h"
+#include "encoder.h"
+
+void IDCT32 (double *, double *, int);
+
+/***********************************************************************
+ An implementation of a modified window subband as seen in Kumar & Zubair's
+"A high performance software implentation of mpeg audio encoder"
+I think from IEEE ASCAP 1996 proceedings
+
+input: shift in 32*12 (384) new samples into a 864 point buffer.
+ch - which channel we're looking at.
+
+This routine basically does 12 calls to window subband all in one go.
+Not yet called in code. here for testing only.
+************************************************************************/
+#define NEWWS
+void window_subband12 (short **buffer, int ch)
+{
+ static double x[2][864]; /* 2 channels, 864 buffer for each */
+ double *xk;
+ double t[12]; /* a temp buffer for summing values */
+ double y[12][64]; /* 12 output arrays of 64 values */
+ int i, j, k, m;
+ static double c[512]; /* enwindow array */
+ static int init = 0;
+ double c0;
+
+ if (!init) {
+ read_ana_window (c);
+ printf ("done init\n");
+ init++;
+ }
+
+ xk = x[ch]; /* an easier way of referencing the array */
+
+ /* shift 384 new samples into the buffer */
+ for (i = 863; i >= 384; i--)
+ xk[i] = xk[i - 384];
+ for (i = 383; i >= 0; i--)
+ xk[i] = (double) *(*buffer)++ / SCALE;
+
+ for (j = 0; j < 64; j++) {
+ for (k = 0; k < 12; k++)
+ t[k] = 0;
+ for (i = 0; i < 8; i++) {
+ m = i * 64 + j;
+ c0 = c[m];
+ t[0] += c0 * xk[m + 352];
+ t[1] += c0 * xk[m + 320];
+ t[2] += c0 * xk[m + 288];
+ t[3] += c0 * xk[m + 256];
+ t[4] += c0 * xk[m + 224];
+ t[5] += c0 * xk[m + 192];
+ t[6] += c0 * xk[m + 160];
+ t[7] += c0 * xk[m + 128];
+ t[8] += c0 * xk[m + 96];
+ t[9] += c0 * xk[m + 64];
+ t[10] += c0 * xk[m + 32];
+ t[11] += c0 * xk[m];
+ }
+ for (i = 0; i < 12; i++) {
+ y[i][j] = t[i];
+ }
+ }
+#define DB1x
+#ifdef DB1
+ for (i = 0; i < 12; i++) {
+ printf ("--%i--\n", i);
+ for (j = 0; j < 64; j++) {
+ printf ("%f\t", y[i][j]);
+ if ((j + 1) % 4 == 0)
+ printf ("\n");
+ }
+ }
+ exit (0);
+#endif
+}
+
+
+/************************************************************************/
+/*
+/* read_ana_window()
+/*
+/* PURPOSE: Reads encoder window file "enwindow" into array #ana_win#
+/*
+/************************************************************************/
+
+void read_ana_window (ana_win)
+ double ana_win[HAN_SIZE];
+{
+ int i, j[4];
+ FILE *fp;
+ double f[4];
+ char t[150];
+
+ if (!(fp = OpenTableFile ("enwindow"))) {
+ printf ("Please check analysis window table 'enwindow'\n");
+ exit (1);
+ }
+ for (i = 0; i < 512; i += 4) {
+ fgets (t, 150, fp);
+ sscanf (t, "C[%d] = %lf C[%d] = %lf C[%d] = %lf C[%d] = %lf\n", j, f,
+ j + 1, f + 1, j + 2, f + 2, j + 3, f + 3);
+ if (i == j[0]) {
+ ana_win[i] = f[0];
+ ana_win[i + 1] = f[1];
+ ana_win[i + 2] = f[2];
+ ana_win[i + 3] = f[3];
+ } else {
+ printf ("Check index in analysis window table\n");
+ exit (1);
+ }
+ fgets (t, 150, fp);
+ }
+ fclose (fp);
+}
+
+/************************************************************************/
+/*
+/* window_subband()
+/*
+/* PURPOSE: Overlapping window on PCM samples
+/*
+/* SEMANTICS:
+/* 32 16-bit pcm samples are scaled to fractional 2's complement and
+/* concatenated to the end of the window buffer #x#. The updated window
+/* buffer #x# is then windowed by the analysis window #c# to produce the
+/* windowed sample #z#
+/*
+/************************************************************************/
+#ifdef COMBWS
+void window_subband (short **buffer, double s[SBLIMIT], int k, int sblimit)
+{
+ typedef double XX[2][HAN_SIZE];
+ static XX *x;
+ double *xk;
+ int i;
+ static int off[2] = { 0, 0 };
+ static char init = 0;
+ double t;
+ static double enwindow[512];
+ double *ep0, *ep1, *ep2, *ep3, *ep4, *ep5, *ep6, *ep7;
+ double z[64];
+ double yprime[32];
+ if (!init) {
+ read_ana_window (enwindow);
+ x = (XX *) mem_alloc (sizeof (XX), "x");
+ memset (x, 0, 2 * HAN_SIZE * sizeof (double));
+ init = 1;
+ }
+ xk = (*x)[k];
+
+ /* replace 32 oldest samples with 32 new samples */
+ for (i = 0; i < 32; i++)
+ xk[31 - i + off[k]] = (double) *(*buffer)++ / SCALE;
+
+ ep0 = &enwindow[0];
+ ep1 = &enwindow[64];
+ ep2 = &enwindow[128];
+ ep3 = &enwindow[192];
+ ep4 = &enwindow[256];
+ ep5 = &enwindow[320];
+ ep6 = &enwindow[384];
+ ep7 = &enwindow[448];
+
+ /* shift samples into proper window positions */
+ for (i = 0; i < 64; i++) {
+ t = xk[(i + off[k]) & (512 - 1)] * *ep0++;
+ t += xk[(i + 64 + off[k]) & (512 - 1)] * *ep1++;
+ t += xk[(i + 128 + off[k]) & (512 - 1)] * *ep2++;
+ t += xk[(i + 192 + off[k]) & (512 - 1)] * *ep3++;
+ t += xk[(i + 256 + off[k]) & (512 - 1)] * *ep4++;
+ t += xk[(i + 320 + off[k]) & (512 - 1)] * *ep5++;
+ t += xk[(i + 384 + off[k]) & (512 - 1)] * *ep6++;
+ t += xk[(i + 448 + off[k]) & (512 - 1)] * *ep7++;
+ z[i] = t;
+ }
+
+ off[k] += 480; /*offset is modulo (HAN_SIZE-1) */
+ off[k] &= HAN_SIZE - 1;
+
+ yprime[0] = z[16];
+ for (i = 1; i <= 16; i++)
+ yprime[i] = z[i + 16] + z[16 - i];
+ for (i = 17; i <= 31; i++)
+ yprime[i] = z[i + 16] - z[80 - i];
+ IDCT32 (yprime, s, sblimit);
+ /* filter_subband (z,s); */
+}
+#else
+void window_subband (short **buffer, double z[64], int k)
+{
+ typedef double XX[2][HAN_SIZE];
+ static XX *x;
+ double *xk;
+ int i;
+ static int off[2] = { 0, 0 };
+ static char init = 0;
+ double t;
+ static double enwindow[512];
+ double *ep0, *ep1, *ep2, *ep3, *ep4, *ep5, *ep6, *ep7;
+ if (!init) {
+ read_ana_window (enwindow);
+ x = (XX *) mem_alloc (sizeof (XX), "x");
+ memset (x, 0, 2 * HAN_SIZE * sizeof (double));
+ init = 1;
+ }
+ xk = (*x)[k];
+
+ /* replace 32 oldest samples with 32 new samples */
+ /* PDS old code: */
+ /* for (i=0;i<32;i++)
+ xk[31-i+off[k]] = (double) *(*buffer)++/SCALE;
+ */
+ {
+ register double *xk_t = xk + off[k];
+ for (i = 32; i--;)
+ xk_t[i] = (double) *(*buffer)++;
+ }
+
+ ep0 = &enwindow[0];
+ ep1 = &enwindow[64];
+ ep2 = &enwindow[128];
+ ep3 = &enwindow[192];
+ ep4 = &enwindow[256];
+ ep5 = &enwindow[320];
+ ep6 = &enwindow[384];
+ ep7 = &enwindow[448];
+
+ /* shift samples into proper window positions */
+ for (i = 0; i < 64; i++) {
+ t = xk[(i + off[k]) & (512 - 1)] * *ep0++;
+ t += xk[(i + 64 + off[k]) & (512 - 1)] * *ep1++;
+ t += xk[(i + 128 + off[k]) & (512 - 1)] * *ep2++;
+ t += xk[(i + 192 + off[k]) & (512 - 1)] * *ep3++;
+ t += xk[(i + 256 + off[k]) & (512 - 1)] * *ep4++;
+ t += xk[(i + 320 + off[k]) & (512 - 1)] * *ep5++;
+ t += xk[(i + 384 + off[k]) & (512 - 1)] * *ep6++;
+ t += xk[(i + 448 + off[k]) & (512 - 1)] * *ep7++;
+ z[i] = t;
+ }
+
+ off[k] += 480; /*offset is modulo (HAN_SIZE-1) */
+ off[k] &= HAN_SIZE - 1;
+
+}
+#endif
+/************************************************************************/
+/*
+/* create_ana_filter()
+/*
+/* PURPOSE: Calculates the analysis filter bank coefficients
+/*
+/* SEMANTICS:
+/* 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#
+/*
+/************************************************************************/
+
+void create_ana_filter (filter)
+ double filter[SBLIMIT][64];
+{
+ register int i, k;
+
+ for (i = 0; i < 32; i++)
+ for (k = 0; k < 64; k++) {
+ if ((filter[i][k] =
+ 1e9 * cos ((double) ((2 * i + 1) * (16 - k) * PI64))) >= 0)
+ modf (filter[i][k] + 0.5, &filter[i][k]);
+ else
+ modf (filter[i][k] - 0.5, &filter[i][k]);
+ filter[i][k] *= 1e-9;
+ }
+}
+
+/************************************************************************
+*
+* filter_subband()
+*
+* PURPOSE: Calculates the analysis filter bank coefficients
+*
+* SEMANTICS:
+* 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 create_dct_matrix (filter)
+ double filter[16][32];
+{
+ register int i, k;
+
+ for (i = 0; i < 16; i++)
+ for (k = 0; k < 32; k++) {
+ if ((filter[i][k] = 1e9 * cos ((double) ((2 * i + 1) * k * PI64))) >= 0)
+ modf (filter[i][k] + 0.5, &filter[i][k]);
+ else
+ modf (filter[i][k] - 0.5, &filter[i][k]);
+ filter[i][k] *= 1e-9;
+ filter[i][k] /= (double) SCALE; /* PDS */
+ }
+ /* PDS this code could/should be replaced; use simple cos */
+ /* and don't do additional rounding ??? See LAME(e.g.3.34) */
+}
+
+void IDCT32 (xin, xout, sblimit)
+ double *xin, *xout;
+ int sblimit;
+{
+ int i, j;
+ double s0, s1;
+ typedef double MM[16][32];
+ static MM *m = 0;
+ if (m == 0) {
+ m = (MM *) mem_alloc (sizeof (MM), "filter");
+ create_dct_matrix (*m);
+ }
+ /* Only compute subband filter info for frequency ranges which */
+ /* will be really needed/encoded later. Code is "general", but */
+ /* only produces speed up in low qual/bps coding situations. */
+ /* Added/adapted by PDS Oct 1999. */
+ for (i = ((sblimit > 16) ? SBLIMIT - sblimit : sblimit); i--;) {
+ s0 = 0.0;
+ for (j = 0; j < 32; j++) {
+ s0 += (*m)[i][j] * xin[j + 0];
+ }
+ xout[i] = s0;
+ }
+ for (i = SBLIMIT - sblimit; i < 16; i++) {
+ s0 = s1 = 0.0;
+ for (j = 0; j < 32; j += 2) {
+ s0 += (*m)[i][j] * xin[j];
+ s1 += (*m)[i][j + 1] * xin[j + 1];
+ }
+ xout[i] = s0 + s1;
+ xout[31 - i] = s0 - s1;
+ }
+ /* PDS TODO: use pointers instead of arrays ??? */
+ /* PDS TODO: is '--' j loop faster then original */
+ /* code: ' for( j=0; j<32; j+=2 ) { ... } ' ? */
+}
+
+void filter_subband (z, s, sblimit)
+ double z[HAN_SIZE], s[SBLIMIT];
+ int sblimit;
+{
+ double yprime[32];
+ int i, j;
+
+ {
+ yprime[0] = z[16];
+ for (i = 1; i <= 16; i++)
+ yprime[i] = z[i + 16] + z[16 - i];
+ for (i = 17; i <= 31; i++)
+ yprime[i] = z[i + 16] - z[80 - i];
+ IDCT32 (yprime, s, sblimit);
+ }
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-06 07:26:03 +0000