Redo the forward 4x4 dct

The new fdct lowers the round trip sum squared error for a
4x4 block ~0.12. or ~0.008/pixel. For reference, the old
matrix multiply version has average round trip error 1.46
for a 4x4 block.

Thanks to "derf" for his suggestions and references.

Change-Id: I5559d1e81d333b319404ab16b336b739f87afc79

1 files changed, 16 insertions, 125 deletions

diff --git a/vp8/encoder/dct.c b/vp8/encoder/dct.c
index 3075e5853..58e36109c 100644
--- a/vp8/encoder/dct.c
+++ b/vp8/encoder/dct.c
@@ -11,163 +11,54 @@
 
 #include <math.h>
 
-
-static const short dct_matrix2[4][4] =
-{
-    { 23170,  30274,  23170, 12540 },
-    { 23170,  12540, -23170, -30274 },
-    { 23170, -12540, -23170, 30274 },
-    { 23170, -30274,  23170, -12540 }
-};
-
-static const short dct_matrix1[4][4] =
-{
-    { 23170,  23170,  23170,  23170 },
-    { 30274,  12540, -12540, -30274 },
-    { 23170, -23170, -23170,  23170 },
-    { 12540, -30274,  30274, -12540 }
-};
-
-
-#define _1STSTAGESHIFT           14
-#define _1STSTAGEROUNDING        (1<<( _1STSTAGESHIFT-1))
-#define _2NDSTAGESHIFT           16
-#define _2NDSTAGEROUNDING        (1<<( _2NDSTAGESHIFT-1))
-
-// using matrix multiply
 void vp8_short_fdct4x4_c(short *input, short *output, int pitch)
 {
-    int i, j, k;
-    short temp[4][4];
-    int sumtemp;
-    pitch >>= 1;
-
-    for (i = 0; i < 4; i++)
-    {
-        for (j = 0; j < 4; j++)
-        {
-            sumtemp = 0;
-
-            for (k = 0; k < 4; k++)
-            {
-                sumtemp += input[i*pitch+k] * dct_matrix2[k][j];
-
-            }
-
-            temp[i][j] = (short)((sumtemp + _1STSTAGEROUNDING) >> _1STSTAGESHIFT);
-        }
-    }
-
-
-    for (i = 0; i < 4; i++)
-    {
-        for (j = 0; j < 4; j++)
-        {
-            sumtemp = 0;
-
-            for (k = 0; k < 4; k++)
-            {
-                sumtemp += dct_matrix1[i][ k] * temp[k][ j];
-            }
-
-            output[i*4+j] = (short)((sumtemp + _2NDSTAGEROUNDING) >> _2NDSTAGESHIFT);
-        }
-    }
-
-}
-
-
-void vp8_short_fdct8x4_c(short *input, short *output, int pitch)
-{
-    vp8_short_fdct4x4_c(input,   output,    pitch);
-    vp8_short_fdct4x4_c(input + 4, output + 16, pitch);
-}
-
-
-static const signed short x_c1 = 60547;
-static const signed short x_c2 = 46341;
-static const signed short x_c3 = 25080;
-
-void vp8_fast_fdct4x4_c(short *input, short *output, int pitch)
-{
     int i;
     int a1, b1, c1, d1;
-    int a2, b2, c2, d2;
     short *ip = input;
-
     short *op = output;
-    int temp1, temp2;
 
     for (i = 0; i < 4; i++)
     {
-        a1 = (ip[0] + ip[3]) * 2;
-        b1 = (ip[1] + ip[2]) * 2;
-        c1 = (ip[1] - ip[2]) * 2;
-        d1 = (ip[0] - ip[3]) * 2;
-
-        temp1 = a1 + b1;
-        temp2 = a1 - b1;
-
-        op[0] = ((temp1 * x_c2) >> 16) + temp1;
-        op[2] = ((temp2 * x_c2) >> 16) + temp2;
-
-        temp1 = (c1 * x_c3) >> 16;
-        temp2 = ((d1 * x_c1) >> 16) + d1;
+        a1 = ((ip[0] + ip[3])<<3);
+        b1 = ((ip[1] + ip[2])<<3);
+        c1 = ((ip[1] - ip[2])<<3);
+        d1 = ((ip[0] - ip[3])<<3);
 
-        op[1] = temp1 + temp2;
-
-        temp1 = (d1 * x_c3) >> 16;
-        temp2 = ((c1 * x_c1) >> 16) + c1;
+        op[0] = a1 + b1;
+        op[2] = a1 - b1;
 
-        op[3] = temp1 - temp2;
+        op[1] = (c1 * 2217 + d1 * 5352 +  14500)>>12;
+        op[3] = (d1 * 2217 - c1 * 5352 +   7500)>>12;
 
         ip += pitch / 2;
         op += 4;
-    }
 
+    }
     ip = output;
     op = output;
-
     for (i = 0; i < 4; i++)
     {
-
         a1 = ip[0] + ip[12];
         b1 = ip[4] + ip[8];
         c1 = ip[4] - ip[8];
         d1 = ip[0] - ip[12];
 
+        op[0]  = ( a1 + b1 + 7)>>4;
+        op[8]  = ( a1 - b1 + 7)>>4;
 
-        temp1 = a1 + b1;
-        temp2 = a1 - b1;
-
-        a2 = ((temp1 * x_c2) >> 16) + temp1;
-        c2 = ((temp2 * x_c2) >> 16) + temp2;
-
-        temp1 = (c1 * x_c3) >> 16;
-        temp2 = ((d1 * x_c1) >> 16) + d1;
-
-        b2 = temp1 + temp2;
-
-        temp1 = (d1 * x_c3) >> 16;
-        temp2 = ((c1 * x_c1) >> 16) + c1;
-
-        d2 = temp1 - temp2;
-
-
-        op[0]   = (a2 + 1) >> 1;
-        op[4]   = (b2 + 1) >> 1;
-        op[8]   = (c2 + 1) >> 1;
-        op[12]  = (d2 + 1) >> 1;
+        op[4]  =((c1 * 2217 + d1 * 5352 +  12000)>>16) + (d1!=0);
+        op[12] = (d1 * 2217 - c1 * 5352 +  51000)>>16;
 
         ip++;
         op++;
     }
 }
 
-void vp8_fast_fdct8x4_c(short *input, short *output, int pitch)
+void vp8_short_fdct8x4_c(short *input, short *output, int pitch)
 {
-    vp8_fast_fdct4x4_c(input,   output,    pitch);
-    vp8_fast_fdct4x4_c(input + 4, output + 16, pitch);
+    vp8_short_fdct4x4_c(input,   output,    pitch);
+    vp8_short_fdct4x4_c(input + 4, output + 16, pitch);
 }
 
 void vp8_short_walsh4x4_c(short *input, short *output, int pitch)