mpeg2enc/idctdata.c

   1 //
   2 // MMX32 iDCT algorithm  (IEEE-1180 compliant) :: idct_mmx32()
   3 //
   4 // MPEG2AVI
   5 // --------
   6 //  v0.16B33 initial release
   7 //
   8 // This was one of the harder pieces of work to code.
   9 // Intel's app-note focuses on the numerical issues of the algorithm, but
  10 // assumes the programmer is familiar with IDCT mathematics, leaving the
  11 // form of the complete function up to the programmer's imagination.
  12
  13 #include "config.h"
  14 #include "mpeg2enc.h"
  15 #include "global.h"
  16 #include "simd.h"
  17
  18 //  ALGORITHM OVERVIEW
  19 //  ------------------
  20 // I played around with the code for quite a few hours.  I came up
  21 // with *A* working IDCT algorithm, however I'm not sure whether my routine
  22 // is "the correct one."  But rest assured, my code passes all six IEEE
  23 // accuracy tests with plenty of margin.
  24 //
  25 //   My IDCT algorithm consists of 4 steps:
  26 //
  27 //   1) IDCT-row transformation (using the IDCT-row function) on all 8 rows
  28 //      This yields an intermediate 8x8 matrix.
  29 //
  30 //   2) intermediate matrix transpose (mandatory)
  31 //
  32 //   3) IDCT-row transformation (2nd time) on all 8 rows of the intermediate
  33 //      matrix.  The output is the final-result, in transposed form.
  34 //
  35 //   4) post-transformation matrix transpose
  36 //      (not necessary if the input-data is already transposed, this could
  37 //       be done during the MPEG "zig-zag" scan, but since my algorithm
  38 //       requires at least one transpose operation, why not re-use the
  39 //       transpose-code.)
  40 //
  41 //   Although the (1st) and (3rd) steps use the SAME row-transform operation,
  42 //   the (3rd) step uses different shift&round constants (explained later.)
  43 //
  44 //   Also note that the intermediate transpose (2) would not be neccessary,
  45 //   if the subsequent operation were a iDCT-column transformation.  Since
  46 //   we only have the iDCT-row transform, we transpose the intermediate
  47 //   matrix and use the iDCT-row transform a 2nd time.
  48 //
  49 //   I had to change some constants/variables for my method to work :
  50 //
  51 //      As given by Intel, the #defines for SHIFT_INV_COL and RND_INV_COL are
  52 //      wrong.  Not surprising since I'm not using a true column-transform
  53 //      operation, but the row-transform operation (as mentioned earlier.)
  54 //      round_inv_col[], which is given as "4 int16_t" values, should have the
  55 //      same dimensions as round_inv_row[].  The corrected variables are
  56 //      shown.
  57 //
  58 //      Intel's code defines a different table for each each row operation.
  59 //      The tables given are 0/4, 1/7, 2/6, and 5/3.  My code only uses row#0.
  60 //      Using the other rows messes up the overall transform.
  61 //
  62 //   IMPLEMENTATION DETAILs
  63 //   ----------------------
  64 //
  65 //   I divided the algorithm's work into two subroutines,
  66 //    1) idct_mmx32_rows() - transforms 8 rows, then transpose
  67 //    2) idct_mmx32_cols() - transforms 8 rows, then transpose
  68 //       yields final result ("drop-in" direct replacement for INT32 IDCT)
  69 //
  70 //   The 2nd function is a clone of the 1st, with changes made only to the
  71 //   shift&rounding instructions.
  72 //
  73 //      In the 1st function (rows), the shift & round instructions use
  74 //       SHIFT_INV_ROW & round_inv_row[] (renamed to r_inv_row[])
  75 //
  76 //      In the 2nd function (cols)-> r_inv_col[], and
  77 //       SHIFT_INV_COL & round_inv_col[] (renamed to r_inv_col[])
  78 //
  79 //   Each function contains an integrated transpose-operator, which comes
  80 //   AFTER the primary transformation operation.  In the future, I'll optimize
  81 //   the code to do more of the transpose-work "in-place".  Right now, I've
  82 //   left the code as two subroutines and a main calling function, so other
  83 //   people can read the code more easily.
  84 //
  85 //   liaor@umcc.ais.org  http://members.tripod.com/~liaor
  86 //
  87
  88 //;=============================================================================
  89 //;
  90 //;  AP-922   http://developer.intel.com/vtune/cbts/strmsimd
  91 //; These examples contain code fragments for first stage iDCT 8x8
  92 //; (for rows) and first stage DCT 8x8 (for columns)
  93 //;
  94 //;=============================================================================
  95
  96 #define BITS_INV_ACC    4       //; 4 or 5 for IEEE
  97         // 5 yields higher accuracy, but lessens dynamic range on the input matrix
  98 #define SHIFT_INV_ROW   (16 - BITS_INV_ACC)
  99 #define SHIFT_INV_COL   (1 + BITS_INV_ACC +14 )  // changed from Intel's val)
 100
 101
 102 #define RND_INV_ROW             (1 << (SHIFT_INV_ROW-1))
 103 #define RND_INV_COL             (1 << (SHIFT_INV_COL-1))
 104 #define RND_INV_CORR    (RND_INV_COL - 1)               //; correction -1.0 and round
 105
 106 /* TODO: This should *really* be aligned on 16-byte boundaries... */
 107
 108 const int idct_r_inv_row[2] = { RND_INV_ROW, RND_INV_ROW};
 109 const int idct_r_inv_col[2] = {RND_INV_COL, RND_INV_COL};
 110 const int idct_r_inv_corr[2] = {RND_INV_CORR, RND_INV_CORR };
 111
 112 /* Unused and thus redundant...
 113 const long long dct_one_corr = 0x0001000100010001;
 114 */
 115
 116 /*
 117 ;=============================================================================
 118 ;
 119 ; The first stage iDCT 8x8 - inverse DCTs of rows
 120 ;
 121 ;-----------------------------------------------------------------------------
 122 ; The 8-point inverse DCT direct algorithm
 123 ;-----------------------------------------------------------------------------
 124 ;
 125 ; static const int16_t w[32] = {
 126 ; FIX(cos_4_16), FIX(cos_2_16), FIX(cos_4_16), FIX(cos_6_16),
 127 ; FIX(cos_4_16), FIX(cos_6_16), -FIX(cos_4_16), -FIX(cos_2_16),
 128 ; FIX(cos_4_16), -FIX(cos_6_16), -FIX(cos_4_16), FIX(cos_2_16),
 129 ; FIX(cos_4_16), -FIX(cos_2_16), FIX(cos_4_16), -FIX(cos_6_16),
 130 ; FIX(cos_1_16), FIX(cos_3_16), FIX(cos_5_16), FIX(cos_7_16),
 131 ; FIX(cos_3_16), -FIX(cos_7_16), -FIX(cos_1_16), -FIX(cos_5_16),
 132 ; FIX(cos_5_16), -FIX(cos_1_16), FIX(cos_7_16), FIX(cos_3_16),
 133 ; FIX(cos_7_16), -FIX(cos_5_16), FIX(cos_3_16), -FIX(cos_1_16) };
 134 ;
 135 ; #define DCT_8_INV_ROW(x, y)
 136
 137 ;{
 138 ; int a0, a1, a2, a3, b0, b1, b2, b3;
 139 ;
 140 ; a0 =x[0]*w[0]+x[2]*w[1]+x[4]*w[2]+x[6]*w[3];
 141 ; a1 =x[0]*w[4]+x[2]*w[5]+x[4]*w[6]+x[6]*w[7];
 142 ; a2 = x[0] * w[ 8] + x[2] * w[ 9] + x[4] * w[10] + x[6] * w[11];
 143 ; a3 = x[0] * w[12] + x[2] * w[13] + x[4] * w[14] + x[6] * w[15];
 144 ; b0 = x[1] * w[16] + x[3] * w[17] + x[5] * w[18] + x[7] * w[19];
 145 ; b1 = x[1] * w[20] + x[3] * w[21] + x[5] * w[22] + x[7] * w[23];
 146 ; b2 = x[1] * w[24] + x[3] * w[25] + x[5] * w[26] + x[7] * w[27];
 147 ; b3 = x[1] * w[28] + x[3] * w[29] + x[5] * w[30] + x[7] * w[31];
 148 ;
 149 ; y[0] = SHIFT_ROUND ( a0 + b0 );
 150 ; y[1] = SHIFT_ROUND ( a1 + b1 );
 151 ; y[2] = SHIFT_ROUND ( a2 + b2 );
 152 ; y[3] = SHIFT_ROUND ( a3 + b3 );
 153 ; y[4] = SHIFT_ROUND ( a3 - b3 );
 154 ; y[5] = SHIFT_ROUND ( a2 - b2 );
 155 ; y[6] = SHIFT_ROUND ( a1 - b1 );
 156 ; y[7] = SHIFT_ROUND ( a0 - b0 );
 157 ;}
 158 ;
 159 ;-----------------------------------------------------------------------------
 160 ;
 161 ; In this implementation the outputs of the iDCT-1D are multiplied
 162 ; for rows 0,4 - by cos_4_16,
 163 ; for rows 1,7 - by cos_1_16,
 164 ; for rows 2,6 - by cos_2_16,
 165 ; for rows 3,5 - by cos_3_16
 166 ; and are shifted to the left for better accuracy
 167 ;
 168 ; For the constants used,
 169 ; FIX(float_const) = (int16_t) (float_const * (1<<15) + 0.5)
 170 ;
 171 ;=============================================================================
 172 */
 173
 174 /* CONCATENATED TABLE, rows 0,1,2,3,4,5,6,7 (in order )
 175
 176    In our implementation, however, we only use row0 !
 177 */
 178
 179 const int16_t idct_tab_01234567[] = {
 180         //row0, this row is required
 181         16384, 16384, 16384, -16384,    // ; movq-> w06 w04 w02 w00
 182         21407, 8867, 8867, -21407,              // w07 w05 w03 w01
 183         16384, -16384, 16384, 16384,    //; w14 w12 w10 w08
 184         -8867, 21407, -21407, -8867,    //; w15 w13 w11 w09
 185         22725, 12873, 19266, -22725,    //; w22 w20 w18 w16
 186         19266, 4520, -4520, -12873,             //; w23 w21 w19 w17
 187         12873, 4520, 4520, 19266,               //; w30 w28 w26 w24
 188         -22725, 19266, -12873, -22725,  //w31 w29 w27 w25
 189
 190         // the rest of these rows (1-7), aren't used !
 191
 192         //row1
 193         22725, 22725, 22725, -22725,    // ; movq-> w06 w04 w02 w00
 194         29692, 12299, 12299, -29692,    //      ; w07 w05 w03 w01
 195         22725, -22725, 22725, 22725,    //; w14 w12 w10 w08
 196         -12299, 29692, -29692, -12299,  //; w15 w13 w11 w09
 197         31521, 17855, 26722, -31521,    //; w22 w20 w18 w16
 198         26722, 6270, -6270, -17855,             //; w23 w21 w19 w17
 199         17855, 6270, 6270, 26722,               //; w30 w28 w26 w24
 200         -31521, 26722, -17855, -31521,  // w31 w29 w27 w25
 201
 202         //row2
 203         21407, 21407, 21407, -21407,    // ; movq-> w06 w04 w02 w00
 204         27969, 11585, 11585, -27969,    // ; w07 w05 w03 w01
 205         21407, -21407, 21407, 21407,    // ; w14 w12 w10 w08
 206         -11585, 27969, -27969, -11585,  //  ;w15 w13 w11 w09
 207         29692, 16819, 25172, -29692,    // ;w22 w20 w18 w16
 208         25172, 5906, -5906, -16819,     // ;w23 w21 w19 w17
 209         16819, 5906, 5906, 25172,               // ;w30 w28 w26 w24
 210         -29692, 25172, -16819, -29692,  //  ;w31 w29 w27 w25
 211
 212         //row3
 213         19266, 19266, 19266, -19266,    //; movq-> w06 w04 w02 w00
 214         25172, 10426, 10426, -25172,    //; w07 w05 w03 w01
 215         19266, -19266, 19266, 19266,    //; w14 w12 w10 w08
 216         -10426, 25172, -25172, -10426,  //; w15 w13 w11 w09
 217         26722, 15137, 22654, -26722,    //; w22 w20 w18 w16
 218         22654, 5315, -5315, -15137,             //; w23 w21 w19 w17
 219         15137, 5315, 5315, 22654,               //; w30 w28 w26 w24
 220         -26722, 22654, -15137, -26722,  //; w31 w29 w27 w25
 221
 222         //row4
 223         16384, 16384, 16384, -16384,    // ; movq-> w06 w04 w02 w00
 224         21407, 8867, 8867, -21407,              // w07 w05 w03 w01
 225         16384, -16384, 16384, 16384,    //; w14 w12 w10 w08
 226         -8867, 21407, -21407, -8867,    //; w15 w13 w11 w09
 227         22725, 12873, 19266, -22725,    //; w22 w20 w18 w16
 228         19266, 4520, -4520, -12873,             //; w23 w21 w19 w17
 229         12873, 4520, 4520, 19266,               //; w30 w28 w26 w24
 230         -22725, 19266, -12873, -22725,  //w31 w29 w27 w25
 231
 232         //row5
 233         19266, 19266, 19266, -19266,    //; movq-> w06 w04 w02 w00
 234         25172, 10426, 10426, -25172,    //; w07 w05 w03 w01
 235         19266, -19266, 19266, 19266,    //; w14 w12 w10 w08
 236         -10426, 25172, -25172, -10426,  //; w15 w13 w11 w09
 237         26722, 15137, 22654, -26722,    //; w22 w20 w18 w16
 238         22654, 5315, -5315, -15137,             //; w23 w21 w19 w17
 239         15137, 5315, 5315, 22654,               //; w30 w28 w26 w24
 240         -26722, 22654, -15137, -26722,  //; w31 w29 w27 w25
 241
 242         //row6
 243         21407, 21407, 21407, -21407,    // ; movq-> w06 w04 w02 w00
 244         27969, 11585, 11585, -27969,    // ; w07 w05 w03 w01
 245         21407, -21407, 21407, 21407,    // ; w14 w12 w10 w08
 246         -11585, 27969, -27969, -11585,  //  ;w15 w13 w11 w09
 247         29692, 16819, 25172, -29692,    // ;w22 w20 w18 w16
 248         25172, 5906, -5906, -16819,     // ;w23 w21 w19 w17
 249         16819, 5906, 5906, 25172,               // ;w30 w28 w26 w24
 250         -29692, 25172, -16819, -29692,  //  ;w31 w29 w27 w25
 251
 252         //row7
 253         22725, 22725, 22725, -22725,    // ; movq-> w06 w04 w02 w00
 254         29692, 12299, 12299, -29692,    //      ; w07 w05 w03 w01
 255         22725, -22725, 22725, 22725,    //; w14 w12 w10 w08
 256         -12299, 29692, -29692, -12299,  //; w15 w13 w11 w09
 257         31521, 17855, 26722, -31521,    //; w22 w20 w18 w16
 258         26722, 6270, -6270, -17855,             //; w23 w21 w19 w17
 259         17855, 6270, 6270, 26722,               //; w30 w28 w26 w24
 260         -31521, 26722, -17855, -31521}; // w31 w29 w27 w25