cinelerra-5.1/libzmpeg3/audio/dct.C

   1 /*
   2  *
   3  *  This file is part of libmpeg3
   4  *
   5  *  libmpeg3 is free software; you can redistribute it and/or modify
   6  *  it under the terms of the GNU General Public License as published by
   7  *  the Free Software Foundation; either version 2, or (at your option)
   8  *  any later version.
   9  *
  10  *  libmpeg3 is distributed in the hope that it will be useful,
  11  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13  *  GNU General Public License for more details.
  14  *
  15  *  You should have received a copy of the GNU General Public License
  16  *  along with GNU Make; see the file COPYING.  If not, write to
  17  *  the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  18  *
  19  */
  20
  21 /*
  22  * Discrete Cosine Tansform (DCT) for subband synthesis
  23  * optimized for machines with no auto-increment.
  24  * The performance is highly compiler dependend. Maybe
  25  * the dct64.c version for 'normal' processor may be faster
  26  * even for Intel processors.
  27  */
  28
  29 #include "../libzmpeg3.h"
  30
  31 int zaudio_decoder_layer_t::
  32 dct64_1(float *out0, float *out1, float *b1, float *b2, float *samples)
  33 {
  34   float *costab = pnts[0];
  35
  36   b1[0x00] = samples[0x00] + samples[0x1F];
  37   b1[0x01] = samples[0x01] + samples[0x1E];
  38   b1[0x1F] = (samples[0x00] - samples[0x1F]) * costab[0x0];
  39   b1[0x1E] = (samples[0x01] - samples[0x1E]) * costab[0x1];
  40
  41   b1[0x02] = samples[0x02] + samples[0x1D];
  42   b1[0x03] = samples[0x03] + samples[0x1C];
  43   b1[0x1D] = (samples[0x02] - samples[0x1D]) * costab[0x2];
  44   b1[0x1C] = (samples[0x03] - samples[0x1C]) * costab[0x3];
  45
  46   b1[0x04] = samples[0x04] + samples[0x1B];
  47   b1[0x05] = samples[0x05] + samples[0x1A];
  48   b1[0x1B] = (samples[0x04] - samples[0x1B]) * costab[0x4];
  49   b1[0x1A] = (samples[0x05] - samples[0x1A]) * costab[0x5];
  50
  51   b1[0x06] = samples[0x06] + samples[0x19];
  52   b1[0x07] = samples[0x07] + samples[0x18];
  53   b1[0x19] = (samples[0x06] - samples[0x19]) * costab[0x6];
  54   b1[0x18] = (samples[0x07] - samples[0x18]) * costab[0x7];
  55
  56   b1[0x08] = samples[0x08] + samples[0x17];
  57   b1[0x09] = samples[0x09] + samples[0x16];
  58   b1[0x17] = (samples[0x08] - samples[0x17]) * costab[0x8];
  59   b1[0x16] = (samples[0x09] - samples[0x16]) * costab[0x9];
  60
  61   b1[0x0A] = samples[0x0A] + samples[0x15];
  62   b1[0x0B] = samples[0x0B] + samples[0x14];
  63   b1[0x15] = (samples[0x0A] - samples[0x15]) * costab[0xA];
  64   b1[0x14] = (samples[0x0B] - samples[0x14]) * costab[0xB];
  65
  66   b1[0x0C] = samples[0x0C] + samples[0x13];
  67   b1[0x0D] = samples[0x0D] + samples[0x12];
  68   b1[0x13] = (samples[0x0C] - samples[0x13]) * costab[0xC];
  69   b1[0x12] = (samples[0x0D] - samples[0x12]) * costab[0xD];
  70
  71   b1[0x0E] = samples[0x0E] + samples[0x11];
  72   b1[0x0F] = samples[0x0F] + samples[0x10];
  73   b1[0x11] = (samples[0x0E] - samples[0x11]) * costab[0xE];
  74   b1[0x10] = (samples[0x0F] - samples[0x10]) * costab[0xF];
  75
  76   costab = pnts[1];
  77
  78   b2[0x00] = b1[0x00] + b1[0x0F];
  79   b2[0x01] = b1[0x01] + b1[0x0E];
  80   b2[0x0F] = (b1[0x00] - b1[0x0F]) * costab[0];
  81   b2[0x0E] = (b1[0x01] - b1[0x0E]) * costab[1];
  82
  83   b2[0x02] = b1[0x02] + b1[0x0D];
  84   b2[0x03] = b1[0x03] + b1[0x0C];
  85   b2[0x0D] = (b1[0x02] - b1[0x0D]) * costab[2];
  86   b2[0x0C] = (b1[0x03] - b1[0x0C]) * costab[3];
  87
  88   b2[0x04] = b1[0x04] + b1[0x0B];
  89   b2[0x05] = b1[0x05] + b1[0x0A];
  90   b2[0x0B] = (b1[0x04] - b1[0x0B]) * costab[4];
  91   b2[0x0A] = (b1[0x05] - b1[0x0A]) * costab[5];
  92
  93   b2[0x06] = b1[0x06] + b1[0x09];
  94   b2[0x07] = b1[0x07] + b1[0x08];
  95   b2[0x09] = (b1[0x06] - b1[0x09]) * costab[6];
  96   b2[0x08] = (b1[0x07] - b1[0x08]) * costab[7];
  97
  98   /* */
  99
 100   b2[0x10] = b1[0x10] + b1[0x1F];
 101   b2[0x11] = b1[0x11] + b1[0x1E];
 102   b2[0x1F] = (b1[0x1F] - b1[0x10]) * costab[0];
 103   b2[0x1E] = (b1[0x1E] - b1[0x11]) * costab[1];
 104
 105   b2[0x12] = b1[0x12] + b1[0x1D];
 106   b2[0x13] = b1[0x13] + b1[0x1C];
 107   b2[0x1D] = (b1[0x1D] - b1[0x12]) * costab[2];
 108   b2[0x1C] = (b1[0x1C] - b1[0x13]) * costab[3];
 109
 110   b2[0x14] = b1[0x14] + b1[0x1B];
 111   b2[0x15] = b1[0x15] + b1[0x1A];
 112   b2[0x1B] = (b1[0x1B] - b1[0x14]) * costab[4];
 113   b2[0x1A] = (b1[0x1A] - b1[0x15]) * costab[5];
 114
 115   b2[0x16] = b1[0x16] + b1[0x19];
 116   b2[0x17] = b1[0x17] + b1[0x18];
 117   b2[0x19] = (b1[0x19] - b1[0x16]) * costab[6];
 118   b2[0x18] = (b1[0x18] - b1[0x17]) * costab[7];
 119
 120    costab = pnts[2];
 121
 122   b1[0x00] = b2[0x00] + b2[0x07];
 123   b1[0x07] = (b2[0x00] - b2[0x07]) * costab[0];
 124   b1[0x01] = b2[0x01] + b2[0x06];
 125   b1[0x06] = (b2[0x01] - b2[0x06]) * costab[1];
 126   b1[0x02] = b2[0x02] + b2[0x05];
 127   b1[0x05] = (b2[0x02] - b2[0x05]) * costab[2];
 128   b1[0x03] = b2[0x03] + b2[0x04];
 129   b1[0x04] = (b2[0x03] - b2[0x04]) * costab[3];
 130
 131   b1[0x08] = b2[0x08] + b2[0x0F];
 132   b1[0x0F] = (b2[0x0F] - b2[0x08]) * costab[0];
 133   b1[0x09] = b2[0x09] + b2[0x0E];
 134   b1[0x0E] = (b2[0x0E] - b2[0x09]) * costab[1];
 135   b1[0x0A] = b2[0x0A] + b2[0x0D];
 136   b1[0x0D] = (b2[0x0D] - b2[0x0A]) * costab[2];
 137   b1[0x0B] = b2[0x0B] + b2[0x0C];
 138   b1[0x0C] = (b2[0x0C] - b2[0x0B]) * costab[3];
 139
 140   b1[0x10] = b2[0x10] + b2[0x17];
 141   b1[0x17] = (b2[0x10] - b2[0x17]) * costab[0];
 142   b1[0x11] = b2[0x11] + b2[0x16];
 143   b1[0x16] = (b2[0x11] - b2[0x16]) * costab[1];
 144   b1[0x12] = b2[0x12] + b2[0x15];
 145   b1[0x15] = (b2[0x12] - b2[0x15]) * costab[2];
 146   b1[0x13] = b2[0x13] + b2[0x14];
 147   b1[0x14] = (b2[0x13] - b2[0x14]) * costab[3];
 148
 149   b1[0x18] = b2[0x18] + b2[0x1F];
 150   b1[0x1F] = (b2[0x1F] - b2[0x18]) * costab[0];
 151   b1[0x19] = b2[0x19] + b2[0x1E];
 152   b1[0x1E] = (b2[0x1E] - b2[0x19]) * costab[1];
 153   b1[0x1A] = b2[0x1A] + b2[0x1D];
 154   b1[0x1D] = (b2[0x1D] - b2[0x1A]) * costab[2];
 155   b1[0x1B] = b2[0x1B] + b2[0x1C];
 156   b1[0x1C] = (b2[0x1C] - b2[0x1B]) * costab[3];
 157
 158   {
 159     float const cos0 = pnts[3][0];
 160     float const cos1 = pnts[3][1];
 161
 162     b2[0x00] = b1[0x00] + b1[0x03];
 163     b2[0x03] = (b1[0x00] - b1[0x03]) * cos0;
 164     b2[0x01] = b1[0x01] + b1[0x02];
 165     b2[0x02] = (b1[0x01] - b1[0x02]) * cos1;
 166
 167     b2[0x04] = b1[0x04] + b1[0x07];
 168     b2[0x07] = (b1[0x07] - b1[0x04]) * cos0;
 169     b2[0x05] = b1[0x05] + b1[0x06];
 170     b2[0x06] = (b1[0x06] - b1[0x05]) * cos1;
 171
 172     b2[0x08] = b1[0x08] + b1[0x0B];
 173     b2[0x0B] = (b1[0x08] - b1[0x0B]) * cos0;
 174     b2[0x09] = b1[0x09] + b1[0x0A];
 175     b2[0x0A] = (b1[0x09] - b1[0x0A]) * cos1;
 176
 177     b2[0x0C] = b1[0x0C] + b1[0x0F];
 178     b2[0x0F] = (b1[0x0F] - b1[0x0C]) * cos0;
 179     b2[0x0D] = b1[0x0D] + b1[0x0E];
 180     b2[0x0E] = (b1[0x0E] - b1[0x0D]) * cos1;
 181
 182     b2[0x10] = b1[0x10] + b1[0x13];
 183     b2[0x13] = (b1[0x10] - b1[0x13]) * cos0;
 184     b2[0x11] = b1[0x11] + b1[0x12];
 185     b2[0x12] = (b1[0x11] - b1[0x12]) * cos1;
 186
 187     b2[0x14] = b1[0x14] + b1[0x17];
 188     b2[0x17] = (b1[0x17] - b1[0x14]) * cos0;
 189     b2[0x15] = b1[0x15] + b1[0x16];
 190     b2[0x16] = (b1[0x16] - b1[0x15]) * cos1;
 191
 192     b2[0x18] = b1[0x18] + b1[0x1B];
 193     b2[0x1B] = (b1[0x18] - b1[0x1B]) * cos0;
 194     b2[0x19] = b1[0x19] + b1[0x1A];
 195     b2[0x1A] = (b1[0x19] - b1[0x1A]) * cos1;
 196
 197     b2[0x1C] = b1[0x1C] + b1[0x1F];
 198     b2[0x1F] = (b1[0x1F] - b1[0x1C]) * cos0;
 199     b2[0x1D] = b1[0x1D] + b1[0x1E];
 200     b2[0x1E] = (b1[0x1E] - b1[0x1D]) * cos1;
 201   }
 202
 203   {
 204     float const cos0 = pnts[4][0];
 205
 206     b1[0x00] = b2[0x00] + b2[0x01];
 207     b1[0x01] = (b2[0x00] - b2[0x01]) * cos0;
 208     b1[0x02] = b2[0x02] + b2[0x03];
 209     b1[0x03] = (b2[0x03] - b2[0x02]) * cos0;
 210     b1[0x02] += b1[0x03];
 211
 212     b1[0x04] = b2[0x04] + b2[0x05];
 213     b1[0x05] = (b2[0x04] - b2[0x05]) * cos0;
 214     b1[0x06] = b2[0x06] + b2[0x07];
 215     b1[0x07] = (b2[0x07] - b2[0x06]) * cos0;
 216     b1[0x06] += b1[0x07];
 217     b1[0x04] += b1[0x06];
 218     b1[0x06] += b1[0x05];
 219     b1[0x05] += b1[0x07];
 220
 221     b1[0x08] = b2[0x08] + b2[0x09];
 222     b1[0x09] = (b2[0x08] - b2[0x09]) * cos0;
 223     b1[0x0A] = b2[0x0A] + b2[0x0B];
 224     b1[0x0B] = (b2[0x0B] - b2[0x0A]) * cos0;
 225     b1[0x0A] += b1[0x0B];
 226
 227     b1[0x0C] = b2[0x0C] + b2[0x0D];
 228     b1[0x0D] = (b2[0x0C] - b2[0x0D]) * cos0;
 229     b1[0x0E] = b2[0x0E] + b2[0x0F];
 230     b1[0x0F] = (b2[0x0F] - b2[0x0E]) * cos0;
 231     b1[0x0E] += b1[0x0F];
 232     b1[0x0C] += b1[0x0E];
 233     b1[0x0E] += b1[0x0D];
 234     b1[0x0D] += b1[0x0F];
 235
 236     b1[0x10] = b2[0x10] + b2[0x11];
 237     b1[0x11] = (b2[0x10] - b2[0x11]) * cos0;
 238     b1[0x12] = b2[0x12] + b2[0x13];
 239     b1[0x13] = (b2[0x13] - b2[0x12]) * cos0;
 240     b1[0x12] += b1[0x13];
 241
 242     b1[0x14] = b2[0x14] + b2[0x15];
 243     b1[0x15] = (b2[0x14] - b2[0x15]) * cos0;
 244     b1[0x16] = b2[0x16] + b2[0x17];
 245     b1[0x17] = (b2[0x17] - b2[0x16]) * cos0;
 246     b1[0x16] += b1[0x17];
 247     b1[0x14] += b1[0x16];
 248     b1[0x16] += b1[0x15];
 249     b1[0x15] += b1[0x17];
 250
 251     b1[0x18] = b2[0x18] + b2[0x19];
 252     b1[0x19] = (b2[0x18] - b2[0x19]) * cos0;
 253     b1[0x1A] = b2[0x1A] + b2[0x1B];
 254     b1[0x1B] = (b2[0x1B] - b2[0x1A]) * cos0;
 255     b1[0x1A] += b1[0x1B];
 256
 257     b1[0x1C] = b2[0x1C] + b2[0x1D];
 258     b1[0x1D] = (b2[0x1C] - b2[0x1D]) * cos0;
 259     b1[0x1E] = b2[0x1E] + b2[0x1F];
 260     b1[0x1F] = (b2[0x1F] - b2[0x1E]) * cos0;
 261     b1[0x1E] += b1[0x1F];
 262     b1[0x1C] += b1[0x1E];
 263     b1[0x1E] += b1[0x1D];
 264     b1[0x1D] += b1[0x1F];
 265   }
 266
 267   out0[0x10*16] = b1[0x00];
 268   out0[0x10*12] = b1[0x04];
 269   out0[0x10* 8] = b1[0x02];
 270   out0[0x10* 4] = b1[0x06];
 271   out0[0x10* 0] = b1[0x01];
 272   out1[0x10* 0] = b1[0x01];
 273   out1[0x10* 4] = b1[0x05];
 274   out1[0x10* 8] = b1[0x03];
 275   out1[0x10*12] = b1[0x07];
 276
 277   out0[0x10*14] = b1[0x08] + b1[0x0C];
 278   out0[0x10*10] = b1[0x0C] + b1[0x0a];
 279   out0[0x10* 6] = b1[0x0A] + b1[0x0E];
 280   out0[0x10* 2] = b1[0x0E] + b1[0x09];
 281   out1[0x10* 2] = b1[0x09] + b1[0x0D];
 282   out1[0x10* 6] = b1[0x0D] + b1[0x0B];
 283   out1[0x10*10] = b1[0x0B] + b1[0x0F];
 284   out1[0x10*14] = b1[0x0F];
 285
 286   {
 287     float tmp;
 288     tmp = b1[0x18] + b1[0x1C];
 289     out0[0x10*15] = tmp + b1[0x10];
 290     out0[0x10*13] = tmp + b1[0x14];
 291     tmp = b1[0x1C] + b1[0x1A];
 292     out0[0x10*11] = tmp + b1[0x14];
 293     out0[0x10* 9] = tmp + b1[0x12];
 294     tmp = b1[0x1A] + b1[0x1E];
 295     out0[0x10* 7] = tmp + b1[0x12];
 296     out0[0x10* 5] = tmp + b1[0x16];
 297     tmp = b1[0x1E] + b1[0x19];
 298     out0[0x10* 3] = tmp + b1[0x16];
 299     out0[0x10* 1] = tmp + b1[0x11];
 300     tmp = b1[0x19] + b1[0x1D];
 301     out1[0x10* 1] = tmp + b1[0x11];
 302     out1[0x10* 3] = tmp + b1[0x15];
 303     tmp = b1[0x1D] + b1[0x1B];
 304     out1[0x10* 5] = tmp + b1[0x15];
 305     out1[0x10* 7] = tmp + b1[0x13];
 306     tmp = b1[0x1B] + b1[0x1F];
 307     out1[0x10* 9] = tmp + b1[0x13];
 308     out1[0x10*11] = tmp + b1[0x17];
 309     out1[0x10*13] = b1[0x17] + b1[0x1F];
 310     out1[0x10*15] = b1[0x1F];
 311   }
 312   return 0;
 313 }
 314
 315 /*
 316  * the call via dct64 is a trick to force GCC to use
 317  * (new) registers for the b1,b2 pointer to the bufs[xx] field
 318  */
 319 int zaudio_decoder_layer_t::
 320 dct64(float *a, float *b, float *c)
 321 {
 322   float bufs[0x40];
 323   return dct64_1(a, b, bufs, bufs + 0x20, c);
 324 }
 325
 326 /*//////////////////////////////////////////////////////////////// */
 327 /* */
 328 /* 9 Point Inverse Discrete Cosine Transform */
 329 /* */
 330 /* This piece of code is Copyright 1997 Mikko Tommila and is freely usable */
 331 /* by anybody. The algorithm itself is of course in the public domain. */
 332 /* */
 333 /* Again derived heuristically from the 9-point WFTA. */
 334 /* */
 335 /* The algorithm is optimized (?) for speed, not for small rounding errors or */
 336 /* good readability. */
 337 /* */
 338 /* 36 additions, 11 multiplications */
 339 /* */
 340 /* Again this is very likely sub-optimal. */
 341 /* */
 342 /* The code is optimized to use a minimum number of temporary variables, */
 343 /* so it should compile quite well even on 8-register Intel x86 processors. */
 344 /* This makes the code quite obfuscated and very difficult to understand. */
 345 /* */
 346 /* References: */
 347 /* [1] S. Winograd: "On Computing the Discrete Fourier Transform", */
 348 /*     Mathematics of Computation, Volume 32, Number 141, January 1978, */
 349 /*     Pages 175-199 */
 350
 351
 352 /*------------------------------------------------------------------*/
 353 /*                                                                  */
 354 /*    Function: Calculation of the inverse MDCT                     */
 355 /*                                                                  */
 356 /*------------------------------------------------------------------*/
 357
 358 int zaudio_decoder_layer_t::
 359 dct36(float *inbuf, float *o1, float *o2, float *wintab, float *tsbuf)
 360 {
 361   float tmp[18];
 362
 363   {
 364     float *in = inbuf;
 365
 366     in[17]+=in[16]; in[16]+=in[15]; in[15]+=in[14];
 367     in[14]+=in[13]; in[13]+=in[12]; in[12]+=in[11];
 368     in[11]+=in[10]; in[10]+=in[9];  in[9] +=in[8];
 369     in[8] +=in[7];  in[7] +=in[6];  in[6] +=in[5];
 370     in[5] +=in[4];  in[4] +=in[3];  in[3] +=in[2];
 371     in[2] +=in[1];  in[1] +=in[0];
 372
 373     in[17]+=in[15]; in[15]+=in[13]; in[13]+=in[11]; in[11]+=in[9];
 374     in[9] +=in[7];  in[7] +=in[5];  in[5] +=in[3];  in[3] +=in[1];
 375
 376
 377     {
 378       float t3;
 379       {
 380         float t0 = COS6_2 * (in[8] + in[16] - in[4]);
 381         float t1 = COS6_2 * in[12];
 382         float t2 = (t3 = in[0]) - t1 - t1;
 383         tmp[1] = tmp[7] = t2 - t0;
 384         tmp[4] = t2 + t0 + t0;
 385         t2 = COS6_1 * (in[10] + in[14] - in[2]);
 386         t3 += t1;
 387         tmp[1] -= t2;
 388         tmp[7] += t2;
 389       }
 390       {
 391         float t0 = cos9[0] * (in[4] + in[8] );
 392         float t1 = cos9[1] * (in[8] - in[16]);
 393         float t2 = cos9[2] * (in[4] + in[16]);
 394         tmp[2] = tmp[6] = t3 - t0      - t2;
 395         tmp[0] = tmp[8] = t3 + t0 + t1;
 396         tmp[3] = tmp[5] = t3 - t1 + t2;
 397       }
 398     }
 399     {
 400       float t1 = cos18[0] * (in[2]  + in[10]);
 401       float t2 = cos18[1] * (in[10] - in[14]);
 402       float t3 = COS6_1   * in[6];
 403       {
 404         float t0 = t1 + t2 + t3;
 405         tmp[0] += t0;
 406         tmp[8] -= t0;
 407       }
 408
 409       t1 -= t3;
 410       t2 -= t3;
 411       t3 = cos18[2] * (in[2] + in[14]);
 412
 413       t1 += t3;
 414       tmp[3] += t1;
 415       tmp[5] -= t1;
 416
 417       t2 -= t3;
 418       tmp[2] += t2;
 419       tmp[6] -= t2;
 420     }
 421
 422     {
 423       float t1 = COS6_2 * in[13];
 424       float t2 = COS6_2 * (in[9] + in[17] - in[5]);
 425       float t3 = in[1] + t1;
 426       float t4 = in[1] - t1 - t1;
 427       float t5 = t4 - t2;
 428       float t0 = cos9[0] * (in[5] + in[9]);
 429       t1 = cos9[1] * (in[9] - in[17]);
 430       tmp[13] = (t4 + t2 + t2) * tfcos36[17-13];
 431       t2 = cos9[2] * (in[5] + in[17]);
 432       float t6 = t3 - t0 - t2;
 433       t0 += t3 + t1;
 434       t3 += t2 - t1;
 435       t2 = cos18[0] * (in[3]  + in[11]);
 436       t4 = cos18[1] * (in[11] - in[15]);
 437       float t7 = COS6_1 * in[7];
 438
 439       t1 = t2 + t4 + t7;
 440       tmp[17] = (t0 + t1) * tfcos36[17-17];
 441       tmp[9]  = (t0 - t1) * tfcos36[17-9];
 442       t1 = cos18[2] * (in[3] + in[15]);
 443       t2 += t1 - t7;
 444
 445       tmp[14] = (t3 + t2) * tfcos36[17-14];
 446       t0 = COS6_1 * (in[11] + in[15] - in[3]);
 447       tmp[12] = (t3 - t2) * tfcos36[17-12];
 448
 449       t4 -= t1 + t7;
 450
 451       tmp[16] = (t5 - t0) * tfcos36[17-16];
 452       tmp[10] = (t5 + t0) * tfcos36[17-10];
 453       tmp[15] = (t6 + t4) * tfcos36[17-15];
 454       tmp[11] = (t6 - t4) * tfcos36[17-11];
 455     }
 456 #define MACRO(v) { float tmpval; \
 457   tmpval = tmp[(v)] + tmp[17-(v)]; \
 458   out2[9+(v)] = tmpval * w[27+(v)]; \
 459   out2[8-(v)] = tmpval * w[26-(v)]; \
 460   tmpval = tmp[(v)] - tmp[17-(v)]; \
 461   ts[SBLIMIT*(8-(v))] = out1[8-(v)] + tmpval * w[8-(v)]; \
 462   ts[SBLIMIT*(9+(v))] = out1[9+(v)] + tmpval * w[9+(v)]; \
 463   }
 464     {
 465       float *out2 = o2;
 466       float *w = wintab;
 467       float *out1 = o1;
 468       float *ts = tsbuf;
 469
 470       MACRO(0);
 471       MACRO(1);
 472       MACRO(2);
 473       MACRO(3);
 474       MACRO(4);
 475       MACRO(5);
 476       MACRO(6);
 477       MACRO(7);
 478       MACRO(8);
 479     }
 480   }
 481   return 0;
 482 }
 483
 484 /*
 485  * new DCT12
 486  */
 487 int zaudio_decoder_layer_t::
 488 dct12(float *in,float *rawout1,float *rawout2,float *wi,float *ts)
 489 {
 490 #define DCT12_PART1 \
 491   in5 = in[5*3]; \
 492   in5 += (in4 = in[4*3]); \
 493   in4 += (in3 = in[3*3]); \
 494   in3 += (in2 = in[2*3]); \
 495   in2 += (in1 = in[1*3]); \
 496   in1 += (in0 = in[0*3]); \
 497   in5 += in3; in3 += in1; \
 498   in2 *= COS6_1; \
 499   in3 *= COS6_1; \
 500
 501 #define DCT12_PART2 \
 502   in0 += in4 * COS6_2; \
 503   in4 = in0 + in2;  \
 504   in0 -= in2;       \
 505   in1 += in5 * COS6_2; \
 506   in5 = (in1 + in3) * tfcos12[0]; \
 507   in1 = (in1 - in3) * tfcos12[2]; \
 508   in3 = in4 + in5;  \
 509   in4 -= in5;       \
 510   in2 = in0 + in1;  \
 511   in0 -= in1;
 512
 513   {
 514     float in0,in1,in2,in3,in4,in5;
 515     float *out1 = rawout1;
 516     ts[SBLIMIT*0] = out1[0]; ts[SBLIMIT*1] = out1[1]; ts[SBLIMIT*2] = out1[2];
 517     ts[SBLIMIT*3] = out1[3]; ts[SBLIMIT*4] = out1[4]; ts[SBLIMIT*5] = out1[5];
 518
 519     DCT12_PART1
 520     {
 521       float tmp0,tmp1 = (in0 - in4);
 522       {
 523         float tmp2 = (in1 - in5) * tfcos12[1];
 524         tmp0 = tmp1 + tmp2;
 525         tmp1 -= tmp2;
 526       }
 527       ts[(17-1)*SBLIMIT] = out1[17-1] + tmp0 * wi[11-1];
 528       ts[(12+1)*SBLIMIT] = out1[12+1] + tmp0 * wi[6+1];
 529       ts[(6 +1)*SBLIMIT] = out1[6 +1] + tmp1 * wi[1];
 530       ts[(11-1)*SBLIMIT] = out1[11-1] + tmp1 * wi[5-1];
 531     }
 532
 533     DCT12_PART2
 534
 535     ts[(17-0)*SBLIMIT] = out1[17-0] + in2 * wi[11-0];
 536     ts[(12+0)*SBLIMIT] = out1[12+0] + in2 * wi[6+0];
 537     ts[(12+2)*SBLIMIT] = out1[12+2] + in3 * wi[6+2];
 538     ts[(17-2)*SBLIMIT] = out1[17-2] + in3 * wi[11-2];
 539
 540     ts[(6+0)*SBLIMIT]  = out1[6+0] + in0 * wi[0];
 541     ts[(11-0)*SBLIMIT] = out1[11-0] + in0 * wi[5-0];
 542     ts[(6+2)*SBLIMIT]  = out1[6+2] + in4 * wi[2];
 543     ts[(11-2)*SBLIMIT] = out1[11-2] + in4 * wi[5-2];
 544   }
 545
 546   in++;
 547
 548   {
 549     float in0,in1,in2,in3,in4,in5;
 550     float *out2 = rawout2;
 551
 552     DCT12_PART1
 553
 554     {
 555       float tmp0,tmp1 = (in0 - in4);
 556       {
 557         float tmp2 = (in1 - in5) * tfcos12[1];
 558         tmp0 = tmp1 + tmp2;
 559         tmp1 -= tmp2;
 560       }
 561       out2[5-1] = tmp0 * wi[11-1];
 562       out2[0+1] = tmp0 * wi[6+1];
 563       ts[(12+1)*SBLIMIT] += tmp1 * wi[1];
 564       ts[(17-1)*SBLIMIT] += tmp1 * wi[5-1];
 565     }
 566
 567     DCT12_PART2
 568
 569     out2[5-0] = in2 * wi[11-0];
 570     out2[0+0] = in2 * wi[6+0];
 571     out2[0+2] = in3 * wi[6+2];
 572     out2[5-2] = in3 * wi[11-2];
 573
 574     ts[(12+0)*SBLIMIT] += in0 * wi[0];
 575     ts[(17-0)*SBLIMIT] += in0 * wi[5-0];
 576     ts[(12+2)*SBLIMIT] += in4 * wi[2];
 577     ts[(17-2)*SBLIMIT] += in4 * wi[5-2];
 578   }
 579
 580   in++;
 581
 582   {
 583     float in0,in1,in2,in3,in4,in5;
 584     float *out2 = rawout2;
 585     out2[12]=out2[13]=out2[14]=out2[15]=out2[16]=out2[17]=0.0;
 586
 587     DCT12_PART1
 588
 589     {
 590       float tmp0,tmp1 = (in0 - in4);
 591       {
 592         float tmp2 = (in1 - in5) * tfcos12[1];
 593         tmp0 = tmp1 + tmp2;
 594         tmp1 -= tmp2;
 595       }
 596       out2[11-1] = tmp0 * wi[11-1];
 597       out2[6 +1] = tmp0 * wi[6+1];
 598       out2[0+1] += tmp1 * wi[1];
 599       out2[5-1] += tmp1 * wi[5-1];
 600     }
 601
 602     DCT12_PART2
 603
 604     out2[11-0] = in2 * wi[11-0];
 605     out2[6 +0] = in2 * wi[6+0];
 606     out2[6 +2] = in3 * wi[6+2];
 607     out2[11-2] = in3 * wi[11-2];
 608
 609     out2[0+0] += in0 * wi[0];
 610     out2[5-0] += in0 * wi[5-0];
 611     out2[0+2] += in4 * wi[2];
 612     out2[5-2] += in4 * wi[5-2];
 613   }
 614
 615   return 0;
 616 }
 617