libmpg123/src/libmpg123/tabinit.c

   1 /*
   2         tabinit.c: initialize tables...
   3
   4         copyright ?-2008 by the mpg123 project - free software under the terms of the LGPL 2.1
   5         see COPYING and AUTHORS files in distribution or http://mpg123.org
   6         initially written by Michael Hipp
   7 */
   8
   9 #include "mpg123lib_intern.h"
  10 #include "debug.h"
  11
  12 // The precomputed cos tables.
  13 #include "costabs.h"
  14 const real *pnts[] = { cos64,cos32,cos16,cos8,cos4 };
  15
  16 static const long intwinbase[] = {
  17      0,    -1,    -1,    -1,    -1,    -1,    -1,    -2,    -2,    -2,
  18     -2,    -3,    -3,    -4,    -4,    -5,    -5,    -6,    -7,    -7,
  19     -8,    -9,   -10,   -11,   -13,   -14,   -16,   -17,   -19,   -21,
  20    -24,   -26,   -29,   -31,   -35,   -38,   -41,   -45,   -49,   -53,
  21    -58,   -63,   -68,   -73,   -79,   -85,   -91,   -97,  -104,  -111,
  22   -117,  -125,  -132,  -139,  -147,  -154,  -161,  -169,  -176,  -183,
  23   -190,  -196,  -202,  -208,  -213,  -218,  -222,  -225,  -227,  -228,
  24   -228,  -227,  -224,  -221,  -215,  -208,  -200,  -189,  -177,  -163,
  25   -146,  -127,  -106,   -83,   -57,   -29,     2,    36,    72,   111,
  26    153,   197,   244,   294,   347,   401,   459,   519,   581,   645,
  27    711,   779,   848,   919,   991,  1064,  1137,  1210,  1283,  1356,
  28   1428,  1498,  1567,  1634,  1698,  1759,  1817,  1870,  1919,  1962,
  29   2001,  2032,  2057,  2075,  2085,  2087,  2080,  2063,  2037,  2000,
  30   1952,  1893,  1822,  1739,  1644,  1535,  1414,  1280,  1131,   970,
  31    794,   605,   402,   185,   -45,  -288,  -545,  -814, -1095, -1388,
  32  -1692, -2006, -2330, -2663, -3004, -3351, -3705, -4063, -4425, -4788,
  33  -5153, -5517, -5879, -6237, -6589, -6935, -7271, -7597, -7910, -8209,
  34  -8491, -8755, -8998, -9219, -9416, -9585, -9727, -9838, -9916, -9959,
  35  -9966, -9935, -9863, -9750, -9592, -9389, -9139, -8840, -8492, -8092,
  36  -7640, -7134, -6574, -5959, -5288, -4561, -3776, -2935, -2037, -1082,
  37    -70,   998,  2122,  3300,  4533,  5818,  7154,  8540,  9975, 11455,
  38  12980, 14548, 16155, 17799, 19478, 21189, 22929, 24694, 26482, 28289,
  39  30112, 31947, 33791, 35640, 37489, 39336, 41176, 43006, 44821, 46617,
  40  48390, 50137, 51853, 53534, 55178, 56778, 58333, 59838, 61289, 62684,
  41  64019, 65290, 66494, 67629, 68692, 69679, 70590, 71420, 72169, 72835,
  42  73415, 73908, 74313, 74630, 74856, 74992, 75038 };
  43
  44 #ifdef OPT_MMXORSSE
  45 #if !defined(OPT_X86_64) && !defined(OPT_NEON) && !defined(OPT_NEON64) && !defined(OPT_AVX)
  46 void make_decode_tables_mmx_asm(long scaleval, float* decwin_mmx, float *decwins);
  47 void make_decode_tables_mmx(mpg123_handle *fr)
  48 {
  49         debug("MMX decode tables");
  50         /* Take care: The scale should be like before, when we didn't have float output all around. */
  51         make_decode_tables_mmx_asm((long)((fr->lastscale < 0 ? fr->p.outscale : fr->lastscale)*SHORT_SCALE), fr->decwin_mmx, fr->decwins);
  52         debug("MMX decode tables done");
  53 }
  54 #else
  55
  56 /* This mimics round() as defined in C99. We stay C89. */
  57 static int rounded(double f)
  58 {
  59         return (int)(f>0 ? floor(f+0.5) : ceil(f-0.5));
  60 }
  61
  62 /* x86-64 doesn't use asm version */
  63 void make_decode_tables_mmx(mpg123_handle *fr)
  64 {
  65         int i,j,val;
  66         int idx = 0;
  67         short *ptr = (short *)fr->decwins;
  68         /* Scale is always based on 1.0 . */
  69         double scaleval = -0.5*(fr->lastscale < 0 ? fr->p.outscale : fr->lastscale);
  70         debug1("MMX decode tables with scaleval %g", scaleval);
  71         for(i=0,j=0;i<256;i++,j++,idx+=32)
  72         {
  73                 if(idx < 512+16)
  74                 fr->decwin_mmx[idx+16] = fr->decwin_mmx[idx] = DOUBLE_TO_REAL((double) intwinbase[j] * scaleval);
  75
  76                 if(i % 32 == 31)
  77                 idx -= 1023;
  78                 if(i % 64 == 63)
  79                 scaleval = - scaleval;
  80         }
  81
  82         for( /* i=256 */ ;i<512;i++,j--,idx+=32)
  83         {
  84                 if(idx < 512+16)
  85                 fr->decwin_mmx[idx+16] = fr->decwin_mmx[idx] = DOUBLE_TO_REAL((double) intwinbase[j] * scaleval);
  86
  87                 if(i % 32 == 31)
  88                 idx -= 1023;
  89                 if(i % 64 == 63)
  90                 scaleval = - scaleval;
  91         }
  92
  93         for(i=0; i<512; i++) {
  94                 if(i&1) val = rounded(fr->decwin_mmx[i]*0.5);
  95                 else val = rounded(fr->decwin_mmx[i]*-0.5);
  96                 if(val > 32767) val = 32767;
  97                 else if(val < -32768) val = -32768;
  98                 ptr[i] = val;
  99         }
 100         for(i=512; i<512+32; i++) {
 101                 if(i&1) val = rounded(fr->decwin_mmx[i]*0.5);
 102                 else val = 0;
 103                 if(val > 32767) val = 32767;
 104                 else if(val < -32768) val = -32768;
 105                 ptr[i] = val;
 106         }
 107         for(i=0; i<512; i++) {
 108                 val = rounded(fr->decwin_mmx[511-i]*-0.5);
 109                 if(val > 32767) val = 32767;
 110                 else if(val < -32768) val = -32768;
 111                 ptr[512+32+i] = val;
 112         }
 113         debug("decode tables done");
 114 }
 115 #endif
 116 #endif
 117
 118 #ifdef REAL_IS_FIXED
 119 /* Need saturating multiplication that keeps table values in 32 bit range,
 120    with the option to swap sign at will (so -2**31 is out).
 121    This code is far from the decoder core and so assembly optimization might
 122    be overkill. */
 123 static int32_t sat_mul32(int32_t a, int32_t b)
 124 {
 125         int64_t prod = (int64_t)a * (int64_t)b;
 126         /* TODO: record the clipping? An extra flag? */
 127         if(prod >  2147483647L) return  2147483647L;
 128         if(prod < -2147483647L) return -2147483647L;
 129         return (int32_t)prod;
 130 }
 131 #endif
 132
 133 void make_decode_tables(mpg123_handle *fr)
 134 {
 135         int i,j;
 136         int idx = 0;
 137         double scaleval;
 138 #ifdef REAL_IS_FIXED
 139         real scaleval_long;
 140 #endif
 141         /* Scale is always based on 1.0 . */
 142         scaleval = -0.5*(fr->lastscale < 0 ? fr->p.outscale : fr->lastscale);
 143         debug1("decode tables with scaleval %g", scaleval);
 144 #ifdef REAL_IS_FIXED
 145         scaleval_long = DOUBLE_TO_REAL_15(scaleval);
 146         debug1("decode table with fixed scaleval %li", (long)scaleval_long);
 147         if(scaleval_long > 28618 || scaleval_long < -28618)
 148         {
 149                 /* TODO: Limit the scaleval itself or limit the multiplication afterwards?
 150                    The former basically disables significant amplification for fixed-point
 151                    decoders, but avoids (possibly subtle) distortion. */
 152                 /* This would limit the amplification instead:
 153                    scaleval_long = scaleval_long < 0 ? -28618 : 28618; */
 154                 if(NOQUIET) warning("Desired amplification may introduce distortion.");
 155         }
 156 #endif
 157         for(i=0,j=0;i<256;i++,j++,idx+=32)
 158         {
 159                 if(idx < 512+16)
 160 #ifdef REAL_IS_FIXED
 161                 fr->decwin[idx+16] = fr->decwin[idx] =
 162                         REAL_SCALE_WINDOW(sat_mul32(intwinbase[j],scaleval_long));
 163 #else
 164                 fr->decwin[idx+16] = fr->decwin[idx] = DOUBLE_TO_REAL((double) intwinbase[j] * scaleval);
 165 #endif
 166
 167                 if(i % 32 == 31)
 168                 idx -= 1023;
 169                 if(i % 64 == 63)
 170 #ifdef REAL_IS_FIXED
 171                 scaleval_long = - scaleval_long;
 172 #else
 173                 scaleval = - scaleval;
 174 #endif
 175         }
 176
 177         for( /* i=256 */ ;i<512;i++,j--,idx+=32)
 178         {
 179                 if(idx < 512+16)
 180 #ifdef REAL_IS_FIXED
 181                 fr->decwin[idx+16] = fr->decwin[idx] =
 182                         REAL_SCALE_WINDOW(sat_mul32(intwinbase[j],scaleval_long));
 183 #else
 184                 fr->decwin[idx+16] = fr->decwin[idx] = DOUBLE_TO_REAL((double) intwinbase[j] * scaleval);
 185 #endif
 186
 187                 if(i % 32 == 31)
 188                 idx -= 1023;
 189                 if(i % 64 == 63)
 190 #ifdef REAL_IS_FIXED
 191                 scaleval_long = - scaleval_long;
 192 #else
 193                 scaleval = - scaleval;
 194 #endif
 195         }
 196 #if defined(OPT_X86_64) || defined(OPT_ALTIVEC) || defined(OPT_SSE) || defined(OPT_SSE_VINTAGE) || defined(OPT_ARM) || defined(OPT_NEON) || defined(OPT_NEON64) || defined(OPT_AVX)
 197         if(  fr->cpu_opts.type == x86_64
 198           || fr->cpu_opts.type == altivec
 199           || fr->cpu_opts.type == sse
 200           || fr->cpu_opts.type == sse_vintage
 201           || fr->cpu_opts.type == arm
 202           || fr->cpu_opts.type == neon
 203           || fr->cpu_opts.type == neon64
 204           || fr->cpu_opts.type == avx )
 205         { /* for float SSE / AltiVec / ARM decoder */
 206                 for(i=512; i<512+32; i++)
 207                 {
 208                         fr->decwin[i] = (i&1) ? fr->decwin[i] : 0;
 209                 }
 210                 for(i=0; i<512; i++)
 211                 {
 212                         fr->decwin[512+32+i] = -fr->decwin[511-i];
 213                 }
 214 #if defined(OPT_NEON) || defined(OPT_NEON64)
 215                 if(fr->cpu_opts.type == neon || fr->cpu_opts.type == neon64)
 216                 {
 217                         for(i=0; i<512; i+=2)
 218                         {
 219                                 fr->decwin[i] = -fr->decwin[i];
 220                         }
 221                 }
 222 #endif
 223         }
 224 #endif
 225         debug("decode tables done");
 226 }
 227
 228 #ifndef NO_8BIT
 229 int make_conv16to8_table(mpg123_handle *fr)
 230 {
 231   int i;
 232         int mode = fr->af.dec_enc;
 233
 234   /*
 235    * ????: 8.0 is right but on SB cards '2.0' is a better value ???
 236    */
 237   const double mul = 8.0;
 238
 239   if(!fr->conv16to8_buf){
 240     fr->conv16to8_buf = (unsigned char *) malloc(8192);
 241     if(!fr->conv16to8_buf) {
 242       fr->err = MPG123_ERR_16TO8TABLE;
 243       if(NOQUIET) error("Can't allocate 16 to 8 converter table!");
 244       return -1;
 245     }
 246     fr->conv16to8 = fr->conv16to8_buf + 4096;
 247   }
 248
 249         switch(mode)
 250         {
 251         case MPG123_ENC_ULAW_8:
 252         {
 253                 double m=127.0 / log(256.0);
 254                 int c1;
 255
 256                 for(i=-4096;i<4096;i++)
 257                 {
 258                         /* dunno whether this is a valid transformation rule ?!?!? */
 259                         if(i < 0)
 260                         c1 = 127 - (int) (log( 1.0 - 255.0 * (double) i*mul / 32768.0 ) * m);
 261                         else
 262                         c1 = 255 - (int) (log( 1.0 + 255.0 * (double) i*mul / 32768.0 ) * m);
 263                         if(c1 < 0 || c1 > 255)
 264                         {
 265                                 if(NOQUIET) error2("Converror %d %d",i,c1);
 266                                 return -1;
 267                         }
 268                         if(c1 == 0)
 269                         c1 = 2;
 270                         fr->conv16to8[i] = (unsigned char) c1;
 271                 }
 272         }
 273         break;
 274         case MPG123_ENC_SIGNED_8:
 275                 for(i=-4096;i<4096;i++)
 276                 fr->conv16to8[i] = i>>5;
 277         break;
 278         case MPG123_ENC_UNSIGNED_8:
 279                 for(i=-4096;i<4096;i++)
 280                 fr->conv16to8[i] = (i>>5)+128;
 281         break;
 282         case MPG123_ENC_ALAW_8:
 283         {
 284                 /*
 285                         Let's believe Wikipedia (http://en.wikipedia.org/wiki/G.711) that this
 286                         is the correct table:
 287
 288                         s0000000wxyza...        n000wxyz  [0-31] -> [0-15]
 289                         s0000001wxyza...        n001wxyz  [32-63] -> [16-31]
 290                         s000001wxyzab...        n010wxyz  [64-127] -> [32-47]
 291                         s00001wxyzabc...        n011wxyz  [128-255] -> [48-63]
 292                         s0001wxyzabcd...        n100wxyz  [256-511] -> [64-79]
 293                         s001wxyzabcde...        n101wxyz  [512-1023] -> [80-95]
 294                         s01wxyzabcdef...        n110wxyz  [1024-2047] -> [96-111]
 295                         s1wxyzabcdefg...        n111wxyz  [2048-4095] -> [112-127]
 296
 297                         Let's extend to -4096, too.
 298                         Also, bytes are xored with 0x55 for transmission.
 299
 300                         Since it sounds OK, I assume it is fine;-)
 301                 */
 302                 for(i=0; i<64; ++i)
 303                 fr->conv16to8[i] = ((unsigned int)i)>>1;
 304                 for(i=64; i<128; ++i)
 305                 fr->conv16to8[i] = ((((unsigned int)i)>>2) & 0xf) | (2<<4);
 306                 for(i=128; i<256; ++i)
 307                 fr->conv16to8[i] = ((((unsigned int)i)>>3) & 0xf) | (3<<4);
 308                 for(i=256; i<512; ++i)
 309                 fr->conv16to8[i] = ((((unsigned int)i)>>4) & 0xf) | (4<<4);
 310                 for(i=512; i<1024; ++i)
 311                 fr->conv16to8[i] = ((((unsigned int)i)>>5) & 0xf) | (5<<4);
 312                 for(i=1024; i<2048; ++i)
 313                 fr->conv16to8[i] = ((((unsigned int)i)>>6) & 0xf) | (6<<4);
 314                 for(i=2048; i<4096; ++i)
 315                 fr->conv16to8[i] = ((((unsigned int)i)>>7) & 0xf) | (7<<4);
 316
 317                 for(i=-4095; i<0; ++i)
 318                 fr->conv16to8[i] = fr->conv16to8[-i] | 0x80;
 319
 320                 fr->conv16to8[-4096] = fr->conv16to8[-4095];
 321
 322                 for(i=-4096;i<4096;i++)
 323                 {
 324                         /* fr->conv16to8[i] = - i>>5; */
 325                         /* fprintf(stderr, "table %i %i\n", i<<AUSHIFT, fr->conv16to8[i]); */
 326                         fr->conv16to8[i] ^= 0x55;
 327                 }
 328         }
 329         break;
 330         default:
 331                 fr->err = MPG123_ERR_16TO8TABLE;
 332                 if(NOQUIET) error("Unknown 8 bit encoding choice.");
 333                 return -1;
 334         break;
 335         }
 336
 337         return 0;
 338 }
 339 #endif
 340