toolame-02l/pds_subband.c

   1 #include "common.h"
   2 #include "encoder.h"
   3
   4 void IDCT32 (double *, double *, int);
   5
   6 /***********************************************************************
   7  An implementation of a modified window subband as seen in Kumar & Zubair's
   8 "A high performance software implentation of mpeg audio encoder"
   9 I think from IEEE ASCAP 1996 proceedings
  10
  11 input: shift in 32*12 (384) new samples into a 864 point buffer.
  12 ch - which channel we're looking at.
  13
  14 This routine basically does 12 calls to window subband all in one go.
  15 Not yet called in code. here for testing only.
  16 ************************************************************************/
  17 #define NEWWS
  18 void window_subband12 (short **buffer, int ch)
  19 {
  20   static double x[2][864];      /* 2 channels, 864 buffer for each */
  21   double *xk;
  22   double t[12];                 /* a temp buffer for summing values */
  23   double y[12][64];             /* 12 output arrays of 64 values */
  24   int i, j, k, m;
  25   static double c[512];         /* enwindow array */
  26   static int init = 0;
  27   double c0;
  28
  29   if (!init) {
  30     read_ana_window (c);
  31     printf ("done init\n");
  32     init++;
  33   }
  34
  35   xk = x[ch];                   /* an easier way of referencing the array */
  36
  37   /* shift 384 new samples into the buffer */
  38   for (i = 863; i >= 384; i--)
  39     xk[i] = xk[i - 384];
  40   for (i = 383; i >= 0; i--)
  41     xk[i] = (double) *(*buffer)++ / SCALE;
  42
  43   for (j = 0; j < 64; j++) {
  44     for (k = 0; k < 12; k++)
  45       t[k] = 0;
  46     for (i = 0; i < 8; i++) {
  47       m = i * 64 + j;
  48       c0 = c[m];
  49       t[0] += c0 * xk[m + 352];
  50       t[1] += c0 * xk[m + 320];
  51       t[2] += c0 * xk[m + 288];
  52       t[3] += c0 * xk[m + 256];
  53       t[4] += c0 * xk[m + 224];
  54       t[5] += c0 * xk[m + 192];
  55       t[6] += c0 * xk[m + 160];
  56       t[7] += c0 * xk[m + 128];
  57       t[8] += c0 * xk[m + 96];
  58       t[9] += c0 * xk[m + 64];
  59       t[10] += c0 * xk[m + 32];
  60       t[11] += c0 * xk[m];
  61     }
  62     for (i = 0; i < 12; i++) {
  63       y[i][j] = t[i];
  64     }
  65   }
  66 #define DB1x
  67 #ifdef DB1
  68   for (i = 0; i < 12; i++) {
  69     printf ("--%i--\n", i);
  70     for (j = 0; j < 64; j++) {
  71       printf ("%f\t", y[i][j]);
  72       if ((j + 1) % 4 == 0)
  73         printf ("\n");
  74     }
  75   }
  76   exit (0);
  77 #endif
  78 }
  79
  80
  81 /************************************************************************/
  82 /*
  83 /* read_ana_window()
  84 /*
  85 /* PURPOSE:  Reads encoder window file "enwindow" into array #ana_win#
  86 /*
  87 /************************************************************************/
  88
  89 void read_ana_window (ana_win)
  90      double ana_win[HAN_SIZE];
  91 {
  92   int i, j[4];
  93   FILE *fp;
  94   double f[4];
  95   char t[150];
  96
  97   if (!(fp = OpenTableFile ("enwindow"))) {
  98     printf ("Please check analysis window table 'enwindow'\n");
  99     exit (1);
 100   }
 101   for (i = 0; i < 512; i += 4) {
 102     fgets (t, 150, fp);
 103     sscanf (t, "C[%d] = %lf C[%d] = %lf C[%d] = %lf C[%d] = %lf\n", j, f,
 104             j + 1, f + 1, j + 2, f + 2, j + 3, f + 3);
 105     if (i == j[0]) {
 106       ana_win[i] = f[0];
 107       ana_win[i + 1] = f[1];
 108       ana_win[i + 2] = f[2];
 109       ana_win[i + 3] = f[3];
 110     } else {
 111       printf ("Check index in analysis window table\n");
 112       exit (1);
 113     }
 114     fgets (t, 150, fp);
 115   }
 116   fclose (fp);
 117 }
 118
 119 /************************************************************************/
 120 /*
 121 /* window_subband()
 122 /*
 123 /* PURPOSE:  Overlapping window on PCM samples
 124 /*
 125 /* SEMANTICS:
 126 /* 32 16-bit pcm samples are scaled to fractional 2's complement and
 127 /* concatenated to the end of the window buffer #x#. The updated window
 128 /* buffer #x# is then windowed by the analysis window #c# to produce the
 129 /* windowed sample #z#
 130 /*
 131 /************************************************************************/
 132 #ifdef COMBWS
 133 void window_subband (short **buffer, double s[SBLIMIT], int k, int sblimit)
 134 {
 135   typedef double XX[2][HAN_SIZE];
 136   static XX *x;
 137   double *xk;
 138   int i;
 139   static int off[2] = { 0, 0 };
 140   static char init = 0;
 141   double t;
 142   static double enwindow[512];
 143   double *ep0, *ep1, *ep2, *ep3, *ep4, *ep5, *ep6, *ep7;
 144   double z[64];
 145   double yprime[32];
 146   if (!init) {
 147     read_ana_window (enwindow);
 148     x = (XX *) mem_alloc (sizeof (XX), "x");
 149     memset (x, 0, 2 * HAN_SIZE * sizeof (double));
 150     init = 1;
 151   }
 152   xk = (*x)[k];
 153
 154   /* replace 32 oldest samples with 32 new samples */
 155   for (i = 0; i < 32; i++)
 156     xk[31 - i + off[k]] = (double) *(*buffer)++ / SCALE;
 157
 158   ep0 = &enwindow[0];
 159   ep1 = &enwindow[64];
 160   ep2 = &enwindow[128];
 161   ep3 = &enwindow[192];
 162   ep4 = &enwindow[256];
 163   ep5 = &enwindow[320];
 164   ep6 = &enwindow[384];
 165   ep7 = &enwindow[448];
 166
 167   /* shift samples into proper window positions */
 168   for (i = 0; i < 64; i++) {
 169     t = xk[(i + off[k]) & (512 - 1)] * *ep0++;
 170     t += xk[(i + 64 + off[k]) & (512 - 1)] * *ep1++;
 171     t += xk[(i + 128 + off[k]) & (512 - 1)] * *ep2++;
 172     t += xk[(i + 192 + off[k]) & (512 - 1)] * *ep3++;
 173     t += xk[(i + 256 + off[k]) & (512 - 1)] * *ep4++;
 174     t += xk[(i + 320 + off[k]) & (512 - 1)] * *ep5++;
 175     t += xk[(i + 384 + off[k]) & (512 - 1)] * *ep6++;
 176     t += xk[(i + 448 + off[k]) & (512 - 1)] * *ep7++;
 177     z[i] = t;
 178   }
 179
 180   off[k] += 480;                /*offset is modulo (HAN_SIZE-1) */
 181   off[k] &= HAN_SIZE - 1;
 182
 183   yprime[0] = z[16];
 184   for (i = 1; i <= 16; i++)
 185     yprime[i] = z[i + 16] + z[16 - i];
 186   for (i = 17; i <= 31; i++)
 187     yprime[i] = z[i + 16] - z[80 - i];
 188   IDCT32 (yprime, s, sblimit);
 189   /*    filter_subband (z,s);  */
 190 }
 191 #else
 192 void window_subband (short **buffer, double z[64], int k)
 193 {
 194   typedef double XX[2][HAN_SIZE];
 195   static XX *x;
 196   double *xk;
 197   int i;
 198   static int off[2] = { 0, 0 };
 199   static char init = 0;
 200   double t;
 201   static double enwindow[512];
 202   double *ep0, *ep1, *ep2, *ep3, *ep4, *ep5, *ep6, *ep7;
 203   if (!init) {
 204     read_ana_window (enwindow);
 205     x = (XX *) mem_alloc (sizeof (XX), "x");
 206     memset (x, 0, 2 * HAN_SIZE * sizeof (double));
 207     init = 1;
 208   }
 209   xk = (*x)[k];
 210
 211   /* replace 32 oldest samples with 32 new samples */
 212   /* PDS old code:   */
 213   /* for (i=0;i<32;i++)
 214      xk[31-i+off[k]] = (double) *(*buffer)++/SCALE;
 215    */
 216   {
 217     register double *xk_t = xk + off[k];
 218     for (i = 32; i--;)
 219       xk_t[i] = (double) *(*buffer)++;
 220   }
 221
 222   ep0 = &enwindow[0];
 223   ep1 = &enwindow[64];
 224   ep2 = &enwindow[128];
 225   ep3 = &enwindow[192];
 226   ep4 = &enwindow[256];
 227   ep5 = &enwindow[320];
 228   ep6 = &enwindow[384];
 229   ep7 = &enwindow[448];
 230
 231   /* shift samples into proper window positions */
 232   for (i = 0; i < 64; i++) {
 233     t = xk[(i + off[k]) & (512 - 1)] * *ep0++;
 234     t += xk[(i + 64 + off[k]) & (512 - 1)] * *ep1++;
 235     t += xk[(i + 128 + off[k]) & (512 - 1)] * *ep2++;
 236     t += xk[(i + 192 + off[k]) & (512 - 1)] * *ep3++;
 237     t += xk[(i + 256 + off[k]) & (512 - 1)] * *ep4++;
 238     t += xk[(i + 320 + off[k]) & (512 - 1)] * *ep5++;
 239     t += xk[(i + 384 + off[k]) & (512 - 1)] * *ep6++;
 240     t += xk[(i + 448 + off[k]) & (512 - 1)] * *ep7++;
 241     z[i] = t;
 242   }
 243
 244   off[k] += 480;                /*offset is modulo (HAN_SIZE-1) */
 245   off[k] &= HAN_SIZE - 1;
 246
 247 }
 248 #endif
 249 /************************************************************************/
 250 /*
 251 /* create_ana_filter()
 252 /*
 253 /* PURPOSE:  Calculates the analysis filter bank coefficients
 254 /*
 255 /* SEMANTICS:
 256 /* Calculates the analysis filterbank coefficients and rounds to the
 257 /* 9th decimal place accuracy of the filterbank tables in the ISO
 258 /* document.  The coefficients are stored in #filter#
 259 /*
 260 /************************************************************************/
 261
 262 void create_ana_filter (filter)
 263      double filter[SBLIMIT][64];
 264 {
 265   register int i, k;
 266
 267   for (i = 0; i < 32; i++)
 268     for (k = 0; k < 64; k++) {
 269       if ((filter[i][k] =
 270            1e9 * cos ((double) ((2 * i + 1) * (16 - k) * PI64))) >= 0)
 271         modf (filter[i][k] + 0.5, &filter[i][k]);
 272       else
 273         modf (filter[i][k] - 0.5, &filter[i][k]);
 274       filter[i][k] *= 1e-9;
 275     }
 276 }
 277
 278 /************************************************************************
 279 *
 280 * filter_subband()
 281 *
 282 * PURPOSE:  Calculates the analysis filter bank coefficients
 283 *
 284 * SEMANTICS:
 285 *      The windowed samples #z# is filtered by the digital filter matrix #m#
 286 * to produce the subband samples #s#. This done by first selectively
 287 * picking out values from the windowed samples, and then multiplying
 288 * them by the filter matrix, producing 32 subband samples.
 289 *
 290 ************************************************************************/
 291 void create_dct_matrix (filter)
 292      double filter[16][32];
 293 {
 294   register int i, k;
 295
 296   for (i = 0; i < 16; i++)
 297     for (k = 0; k < 32; k++) {
 298       if ((filter[i][k] = 1e9 * cos ((double) ((2 * i + 1) * k * PI64))) >= 0)
 299         modf (filter[i][k] + 0.5, &filter[i][k]);
 300       else
 301         modf (filter[i][k] - 0.5, &filter[i][k]);
 302       filter[i][k] *= 1e-9;
 303       filter[i][k] /= (double) SCALE;   /* PDS */
 304     }
 305   /* PDS this code could/should be replaced; use simple cos  */
 306   /* and don't do additional rounding ??? See LAME(e.g.3.34) */
 307 }
 308
 309 void IDCT32 (xin, xout, sblimit)
 310      double *xin, *xout;
 311      int sblimit;
 312 {
 313   int i, j;
 314   double s0, s1;
 315   typedef double MM[16][32];
 316   static MM *m = 0;
 317   if (m == 0) {
 318     m = (MM *) mem_alloc (sizeof (MM), "filter");
 319     create_dct_matrix (*m);
 320   }
 321   /* Only compute subband filter info for frequency ranges which */
 322   /* will be really needed/encoded later. Code is "general", but */
 323   /* only produces speed up in low qual/bps coding situations.   */
 324   /* Added/adapted by PDS Oct 1999.                              */
 325   for (i = ((sblimit > 16) ? SBLIMIT - sblimit : sblimit); i--;) {
 326     s0 = 0.0;
 327     for (j = 0; j < 32; j++) {
 328       s0 += (*m)[i][j] * xin[j + 0];
 329     }
 330     xout[i] = s0;
 331   }
 332   for (i = SBLIMIT - sblimit; i < 16; i++) {
 333     s0 = s1 = 0.0;
 334     for (j = 0; j < 32; j += 2) {
 335       s0 += (*m)[i][j] * xin[j];
 336       s1 += (*m)[i][j + 1] * xin[j + 1];
 337     }
 338     xout[i] = s0 + s1;
 339     xout[31 - i] = s0 - s1;
 340   }
 341   /* PDS TODO: use pointers instead of arrays ???  */
 342   /* PDS TODO: is '--' j loop faster then original */
 343   /* code:   ' for( j=0; j<32; j+=2 )  { ... } ' ? */
 344 }
 345
 346 void filter_subband (z, s, sblimit)
 347      double z[HAN_SIZE], s[SBLIMIT];
 348      int sblimit;
 349 {
 350   double yprime[32];
 351   int i, j;
 352
 353   {
 354     yprime[0] = z[16];
 355     for (i = 1; i <= 16; i++)
 356       yprime[i] = z[i + 16] + z[16 - i];
 357     for (i = 17; i <= 31; i++)
 358       yprime[i] = z[i + 16] - z[80 - i];
 359     IDCT32 (yprime, s, sblimit);
 360   }
 361 }