toolame-02l/subband.c

   1 #include <stdio.h>
   2 #include <string.h>
   3 #include <math.h>
   4 #include "common.h"
   5 #include "encoder.h"
   6 #include "mem.h"
   7 #include "bitstream.h"
   8 #include "encode.h"
   9 #include "enwindow.h"
  10 #include "subband.h"
  11
  12
  13 #ifdef REFERENCECODE
  14 /************************************************************************
  15 *
  16 * window_subband()
  17 *
  18 * PURPOSE:  Overlapping window on PCM samples
  19 *
  20 * SEMANTICS:
  21 * 32 16-bit pcm samples are scaled to fractional 2's complement and
  22 * concatenated to the end of the window buffer #x#. The updated window
  23 * buffer #x# is then windowed by the analysis window #c# to produce the
  24 * windowed sample #z#
  25 *
  26 ************************************************************************/
  27
  28 void window_subband (short **buffer, double z[64], int k)
  29 {
  30   typedef double XX[2][HAN_SIZE];
  31   static XX *x;
  32   double *xk;
  33   int i;
  34   static int off[2] = { 0, 0 };
  35   static char init = 0;
  36   double t;
  37   double *ep0, *ep1, *ep2, *ep3, *ep4, *ep5, *ep6, *ep7;
  38   if (!init) {
  39     x = (XX *) mem_alloc (sizeof (XX), "x");
  40     memset (x, 0, 2 * HAN_SIZE * sizeof (double));
  41     init = 1;
  42   }
  43   xk = (*x)[k];
  44
  45   /* replace 32 oldest samples with 32 new samples */
  46   for (i = 0; i < 32; i++)
  47     xk[31 - i + off[k]] = (double) *(*buffer)++ / SCALE;
  48
  49   ep0 = &enwindow[0];
  50   ep1 = &enwindow[64];
  51   ep2 = &enwindow[128];
  52   ep3 = &enwindow[192];
  53   ep4 = &enwindow[256];
  54   ep5 = &enwindow[320];
  55   ep6 = &enwindow[384];
  56   ep7 = &enwindow[448];
  57
  58   /* shift samples into proper window positions */
  59   for (i = 0; i < 64; i++) {
  60     t = xk[(i + off[k]) & (512 - 1)] * *ep0++;
  61     t += xk[(i + 64 + off[k]) & (512 - 1)] * *ep1++;
  62     t += xk[(i + 128 + off[k]) & (512 - 1)] * *ep2++;
  63     t += xk[(i + 192 + off[k]) & (512 - 1)] * *ep3++;
  64     t += xk[(i + 256 + off[k]) & (512 - 1)] * *ep4++;
  65     t += xk[(i + 320 + off[k]) & (512 - 1)] * *ep5++;
  66     t += xk[(i + 384 + off[k]) & (512 - 1)] * *ep6++;
  67     t += xk[(i + 448 + off[k]) & (512 - 1)] * *ep7++;
  68     z[i] = t;
  69   }
  70
  71   off[k] += 480;                /*offset is modulo (HAN_SIZE-1) */
  72   off[k] &= HAN_SIZE - 1;
  73
  74 }
  75
  76
  77 /************************************************************************
  78 *
  79 * filter_subband()
  80 *
  81 * PURPOSE:  Calculates the analysis filter bank coefficients
  82 *
  83 * SEMANTICS:
  84 *      The windowed samples #z# is filtered by the digital filter matrix #m#
  85 * to produce the subband samples #s#. This done by first selectively
  86 * picking out values from the windowed samples, and then multiplying
  87 * them by the filter matrix, producing 32 subband samples.
  88 *
  89 ************************************************************************/
  90 void filter_subband (double z[HAN_SIZE], double s[SBLIMIT])
  91 {
  92   double yprime[32];
  93   register int i, j;
  94
  95   static double m[16][32];
  96   static int init = 0;
  97
  98   if (init == 0) {
  99     init++;
 100     create_dct_matrix (m);
 101   }
 102
 103   yprime[0] = z[16];
 104   for (i = 1; i <= 16; i++)
 105     yprime[i] = z[i + 16] + z[16 - i];
 106   for (i = 17; i <= 31; i++)
 107     yprime[i] = z[i + 16] - z[80 - i];
 108
 109   for (i = 15; i >= 0; i--) {
 110     register double s0 = 0.0, s1 = 0.0;
 111     register double *mp = m[i];
 112     register double *xinp = yprime;
 113     for (j = 0; j < 8; j++) {
 114       s0 += *mp++ * *xinp++;
 115       s1 += *mp++ * *xinp++;
 116       s0 += *mp++ * *xinp++;
 117       s1 += *mp++ * *xinp++;
 118     }
 119     s[i] = s0 + s1;
 120     s[31 - i] = s0 - s1;
 121   }
 122 }
 123 #endif //REFERENCECODE
 124
 125 void create_dct_matrix (double filter[16][32])
 126 {
 127   register int i, k;
 128
 129   for (i = 0; i < 16; i++)
 130     for (k = 0; k < 32; k++) {
 131       if ((filter[i][k] = 1e9 * cos ((double) ((2 * i + 1) * k * PI64))) >= 0)
 132         modf (filter[i][k] + 0.5, &filter[i][k]);
 133       else
 134         modf (filter[i][k] - 0.5, &filter[i][k]);
 135       filter[i][k] *= 1e-9;
 136     }
 137 }
 138
 139
 140 #ifdef NEWWS
 141 /***********************************************************************
 142  An implementation of a modified window subband as seen in Kumar & Zubair's
 143 "A high performance software implentation of mpeg audio encoder"
 144 I think from IEEE ASCAP 1996 proceedings
 145
 146 input: shift in 32*12 (384) new samples into a 864 point buffer.
 147 ch - which channel we're looking at.
 148
 149 This routine basically does 12 calls to window subband all in one go.
 150 Not yet called in code. here for testing only.
 151 ************************************************************************/
 152 void window_subband12 (short **buffer, int ch)
 153 {
 154   static double x[2][864];      /* 2 channels, 864 buffer for each */
 155   double *xk;
 156   double t[12];                 /* a temp buffer for summing values */
 157   double y[12][64];             /* 12 output arrays of 64 values */
 158   int i, j, k, m;
 159   static double c[512];         /* enwindow array */
 160   static int init = 0;
 161   double c0;
 162
 163   xk = x[ch];                   /* an easier way of referencing the array */
 164
 165   /* shift 384 new samples into the buffer */
 166   for (i = 863; i >= 384; i--)
 167     xk[i] = xk[i - 384];
 168   for (i = 383; i >= 0; i--)
 169     xk[i] = (double) *(*buffer)++ / SCALE;
 170
 171   for (j = 0; j < 64; j++) {
 172     for (k = 0; k < 12; k++)
 173       t[k] = 0;
 174     for (i = 0; i < 8; i++) {
 175       m = i * 64 + j;
 176       c0 = c[m];
 177       t[0] += c0 * xk[m + 352];
 178       t[1] += c0 * xk[m + 320];
 179       t[2] += c0 * xk[m + 288];
 180       t[3] += c0 * xk[m + 256];
 181       t[4] += c0 * xk[m + 224];
 182       t[5] += c0 * xk[m + 192];
 183       t[6] += c0 * xk[m + 160];
 184       t[7] += c0 * xk[m + 128];
 185       t[8] += c0 * xk[m + 96];
 186       t[9] += c0 * xk[m + 64];
 187       t[10] += c0 * xk[m + 32];
 188       t[11] += c0 * xk[m];
 189     }
 190     for (i = 0; i < 12; i++) {
 191       y[i][j] = t[i];
 192     }
 193   }
 194 }
 195 #endif /* NEWWS */
 196
 197
 198 //____________________________________________________________________________
 199 //____ WindowFilterSubband() _________________________________________
 200 //____ RS&A - Feb 2003 _______________________________________________________
 201 void WindowFilterSubband (short *pBuffer, int ch, double s[SBLIMIT])
 202 {
 203   register int i, j;
 204   int pa, pb, pc, pd, pe, pf, pg, ph;
 205   double t;
 206   double *dp, *dp2;
 207   double *pEnw;
 208   double y[64];
 209   double yprime[32];
 210
 211   static double x[2][512];
 212   static double m[16][32];
 213   static int init = 0;
 214   static int off[2];
 215   static int half[2];
 216
 217   if (init == 0) {
 218     init++;
 219     off[0] = 0;
 220     off[1] = 0;
 221     half[0] = 0;
 222     half[1] = 0;
 223     for (i = 0; i < 2; i++)
 224       for (j = 0; j < 512; j++)
 225         x[i][j] = 0;
 226     create_dct_matrix (m);
 227   }
 228
 229   dp = x[ch] + off[ch] + half[ch] * 256;
 230
 231   /* replace 32 oldest samples with 32 new samples */
 232   for (i = 0; i < 32; i++)
 233     dp[(31 - i) * 8] = (double) pBuffer[i] / SCALE;
 234
 235   // looks like "school example" but does faster ...
 236   dp = (x[ch] + half[ch] * 256);
 237   pa = off[ch];
 238   pb = (pa + 1) % 8;
 239   pc = (pa + 2) % 8;
 240   pd = (pa + 3) % 8;
 241   pe = (pa + 4) % 8;
 242   pf = (pa + 5) % 8;
 243   pg = (pa + 6) % 8;
 244   ph = (pa + 7) % 8;
 245
 246   for (i = 0; i < 32; i++) {
 247     dp2 = dp + i * 8;
 248     pEnw = enwindow + i;
 249     t = dp2[pa] * pEnw[0];
 250     t += dp2[pb] * pEnw[64];
 251     t += dp2[pc] * pEnw[128];
 252     t += dp2[pd] * pEnw[192];
 253     t += dp2[pe] * pEnw[256];
 254     t += dp2[pf] * pEnw[320];
 255     t += dp2[pg] * pEnw[384];
 256     t += dp2[ph] * pEnw[448];
 257     y[i] = t;
 258   }
 259
 260   yprime[0] = y[16];            // Michael Chen´s dct filter
 261
 262   dp = half[ch] ? x[ch] : (x[ch] + 256);
 263   pa = half[ch] ? (off[ch] + 1) & 7 : off[ch];
 264   pb = (pa + 1) % 8;
 265   pc = (pa + 2) % 8;
 266   pd = (pa + 3) % 8;
 267   pe = (pa + 4) % 8;
 268   pf = (pa + 5) % 8;
 269   pg = (pa + 6) % 8;
 270   ph = (pa + 7) % 8;
 271
 272   for (i = 0; i < 32; i++) {
 273     dp2 = dp + i * 8;
 274     pEnw = enwindow + i + 32;
 275     t = dp2[pa] * pEnw[0];
 276     t += dp2[pb] * pEnw[64];
 277     t += dp2[pc] * pEnw[128];
 278     t += dp2[pd] * pEnw[192];
 279     t += dp2[pe] * pEnw[256];
 280     t += dp2[pf] * pEnw[320];
 281     t += dp2[pg] * pEnw[384];
 282     t += dp2[ph] * pEnw[448];
 283     y[i + 32] = t;
 284     // 1st pass on Michael Chen´s dct filter
 285     if (i > 0 && i < 17)
 286       yprime[i] = y[i + 16] + y[16 - i];
 287   }
 288
 289   // 2nd pass on Michael Chen´s dct filter
 290   for (i = 17; i < 32; i++)
 291     yprime[i] = y[i + 16] - y[80 - i];
 292
 293   for (i = 15; i >= 0; i--) {
 294     register double s0 = 0.0, s1 = 0.0;
 295     register double *mp = m[i];
 296     register double *xinp = yprime;
 297     for (j = 0; j < 8; j++) {
 298       s0 += *mp++ * *xinp++;
 299       s1 += *mp++ * *xinp++;
 300       s0 += *mp++ * *xinp++;
 301       s1 += *mp++ * *xinp++;
 302     }
 303     s[i] = s0 + s1;
 304     s[31 - i] = s0 - s1;
 305   }
 306
 307   half[ch] = (half[ch] + 1) & 1;
 308   if (half[ch] == 1)
 309     off[ch] = (off[ch] + 7) & 7;
 310 }