@@ -99,6 +99,7 @@ Configuration options:

   --disable-fastdiv        disable table-based division

   --enable-small           optimize for size instead of speed

   --disable-aandct         disable AAN DCT code

+  --disable-dct            disable DCT code

   --disable-fft            disable FFT code

   --disable-golomb         disable Golomb code

   --disable-lpc            disable LPC code

@@ -862,6 +863,7 @@ CONFIG_LIST="

     avisynth

     beos_netserver

     bzlib

+    dct

     doc

     fastdiv

     ffmpeg

@@ -27,6 +27,7 @@ OBJS = allcodecs.o                                                      \

 # parts needed for many different codecs

 OBJS-$(CONFIG_AANDCT)                  += aandcttab.o

 OBJS-$(CONFIG_ENCODERS)                += faandct.o jfdctfst.o jfdctint.o

+OBJS-$(CONFIG_DCT)                     += dct.o

 FFT-OBJS-$(CONFIG_HARDCODED_TABLES)    += cos_tables.o

 OBJS-$(CONFIG_FFT)                     += fft.o $(FFT-OBJS-yes)

 OBJS-$(CONFIG_GOLOMB)                  += golomb.o

new file mode 100644

@@ -0,0 +1,139 @@

+/*

+ * (I)DCT Transforms

+ * Copyright (c) 2009 Peter Ross <pross@xvid.org>

+ * Copyright (c) 2010 Alex Converse <alex.converse@gmail.com>

+ * Copyright (c) 2010 Vitor Sessak

+ *

+ * This file is part of FFmpeg.

+ *

+ * FFmpeg is free software; you can redistribute it and/or

+ * modify it under the terms of the GNU Lesser General Public

+ * License as published by the Free Software Foundation; either

+ * version 2.1 of the License, or (at your option) any later version.

+ *

+ * FFmpeg is distributed in the hope that it will be useful,

+ * but WITHOUT ANY WARRANTY; without even the implied warranty of

+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

+ * Lesser General Public License for more details.

+ *

+ * You should have received a copy of the GNU Lesser General Public

+ * License along with FFmpeg; if not, write to the Free Software

+ * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

+ */

+

+/**

+ * @file libavcodec/dct.c

+ * (Inverse) Discrete Cosine Transforms. These are also known as the

+ * type II and type III DCTs respectively.

+ */

+

+#include <math.h>

+#include "dsputil.h"

+

+av_cold int ff_dct_init(DCTContext *s, int nbits, int inverse)

+{

+    int n = 1 << nbits;

+    int i;

+

+    s->nbits    = nbits;

+    s->inverse  = inverse;

+

+    ff_init_ff_cos_tabs(nbits+2);

+

+    s->costab = ff_cos_tabs[nbits+2];

+

+    s->csc2 = av_malloc(n/2 * sizeof(FFTSample));

+

+    if (ff_rdft_init(&s->rdft, nbits, inverse) < 0) {

+        av_free(s->csc2);

+        return -1;

+    }

+

+    for (i = 0; i < n/2; i++)

+        s->csc2[i] = 0.5 / sin((M_PI / (2*n) * (2*i + 1)));

+

+    return 0;

+}

+

+/* sin((M_PI * x / (2*n)) */

+#define SIN(s,n,x) (s->costab[(n) - (x)])

+

+/* cos((M_PI * x / (2*n)) */

+#define COS(s,n,x) (s->costab[x])

+

+static void ff_dct_calc_c(DCTContext *ctx, FFTSample *data)

+{

+    int n = 1 << ctx->nbits;

+    int i;

+

+    if (ctx->inverse) {

+        float next = data[n - 1];

+        float inv_n = 1.0f / n;

+

+        for (i = n - 2; i >= 2; i -= 2) {

+            float val1 = data[i    ];

+            float val2 = data[i - 1] - data[i + 1];

+            float c = COS(ctx, n, i);

+            float s = SIN(ctx, n, i);

+

+            data[i    ] = c * val1 + s * val2;

+            data[i + 1] = s * val1 - c * val2;

+        }

+

+        data[1] = 2 * next;

+

+        ff_rdft_calc(&ctx->rdft, data);

+

+        for (i = 0; i < n / 2; i++) {

+            float tmp1 = data[i        ] * inv_n;

+            float tmp2 = data[n - i - 1] * inv_n;

+            float csc = ctx->csc2[i] * (tmp1 - tmp2);

+

+            tmp1 += tmp2;

+            data[i        ] = tmp1 + csc;

+            data[n - i - 1] = tmp1 - csc;

+        }

+    } else {

+        float next;

+        for (i=0; i < n/2; i++) {

+            float tmp1 = data[i        ];

+            float tmp2 = data[n - i - 1];

+            float s = SIN(ctx, n, 2*i + 1);

+

+            s *= tmp1 - tmp2;

+            tmp1 = (tmp1 + tmp2) * 0.5f;

+

+            data[i    ] = tmp1 + s;

+            data[n-i-1] = tmp1 - s;

+        }

+

+        ff_rdft_calc(&ctx->rdft, data);

+

+        next = data[1] * 0.5;

+        data[1] *= -1;

+

+        for (i = n - 2; i >= 0; i -= 2) {

+            float inr = data[i    ];

+            float ini = data[i + 1];

+            float c = COS(ctx, n, i);

+            float s = SIN(ctx, n, i);

+

+            data[i  ] = c * inr + s * ini;

+

+            data[i+1] = next;

+

+            next +=     s * inr - c * ini;

+        }

+    }

+}

+

+void ff_dct_calc(DCTContext *s, FFTSample *data)

+{

+    ff_dct_calc_c(s, data);

+}

+

+av_cold void ff_dct_end(DCTContext *s)

+{

+    ff_rdft_end(&s->rdft);

+    av_free(s->csc2);

+}

@@ -899,6 +899,26 @@ int ff_rdft_init(RDFTContext *s, int nbits, enum RDFTransformType trans);

 void ff_rdft_calc(RDFTContext *s, FFTSample *data);

 void ff_rdft_end(RDFTContext *s);

+/* Discrete Cosine Transform */

+

+typedef struct {

+    int nbits;

+    int inverse;

+    FFTSample *data;

+    RDFTContext rdft;

+    const float *costab;

+    FFTSample *csc2;

+} DCTContext;

+

+/**

+ * Sets up (Inverse)DCT.

+ * @param nbits           log2 of the length of the input array

+ * @param inverse         >0 forward transform, <0 inverse transform

+ */

+int  ff_dct_init(DCTContext *s, int nbits, int inverse);

+void ff_dct_calc(DCTContext *s, FFTSample *data);

+void ff_dct_end (DCTContext *s);

+

 #define WRAPPER8_16(name8, name16)\

 static int name16(void /*MpegEncContext*/ *s, uint8_t *dst, uint8_t *src, int stride, int h){\

     return name8(s, dst           , src           , stride, h)\

@@ -128,6 +128,39 @@ static void mdct_ref(float *output, float *input, int nbits)

     }

 }

+static void idct_ref(float *output, float *input, int nbits)

+{

+    int n = 1<<nbits;

+    int k, i;

+    double a, s;

+

+    /* do it by hand */

+    for (i = 0; i < n; i++) {

+        s = 0.5 * input[0];

+        for (k = 1; k < n; k++) {

+            a = M_PI*k*(i+0.5) / n;

+            s += input[k] * cos(a);

+        }

+        output[i] = 2 * s / n;

+    }

+}

+static void dct_ref(float *output, float *input, int nbits)

+{

+    int n = 1<<nbits;

+    int k, i;

+    double a, s;

+

+    /* do it by hand */

+    for (k = 0; k < n; k++) {

+        s = 0;

+        for (i = 0; i < n; i++) {

+            a = M_PI*k*(i+0.5) / n;

+            s += input[i] * cos(a);

+        }

+        output[k] = s;

+    }

+}

+

 static float frandom(AVLFG *prng)

 {

@@ -166,6 +199,7 @@ static void help(void)

            "-h     print this help\n"

            "-s     speed test\n"

            "-m     (I)MDCT test\n"

+           "-d     (I)DCT test\n"

            "-i     inverse transform test\n"

            "-n b   set the transform size to 2^b\n"

            "-f x   set scale factor for output data of (I)MDCT to x\n"

@@ -177,6 +211,7 @@ enum tf_transform {

     TRANSFORM_FFT,

     TRANSFORM_MDCT,

     TRANSFORM_RDFT,

+    TRANSFORM_DCT,

 };

 int main(int argc, char **argv)

@@ -190,6 +225,7 @@ int main(int argc, char **argv)

     FFTContext s1, *s = &s1;

     FFTContext m1, *m = &m1;

     RDFTContext r1, *r = &r1;

+    DCTContext d1, *d = &d1;

     int fft_nbits, fft_size, fft_size_2;

     double scale = 1.0;

     AVLFG prng;

@@ -197,7 +233,7 @@ int main(int argc, char **argv)

     fft_nbits = 9;

     for(;;) {

-        c = getopt(argc, argv, "hsimrn:f:");

+        c = getopt(argc, argv, "hsimrdn:f:");

         if (c == -1)

             break;

         switch(c) {

@@ -216,6 +252,9 @@ int main(int argc, char **argv)

         case 'r':

             transform = TRANSFORM_RDFT;

             break;

+        case 'd':

+            transform = TRANSFORM_DCT;

+            break;

         case 'n':

             fft_nbits = atoi(optarg);

             break;

@@ -257,6 +296,13 @@ int main(int argc, char **argv)

         ff_rdft_init(r, fft_nbits, do_inverse ? IRDFT : RDFT);

         fft_ref_init(fft_nbits, do_inverse);

         break;

+    case TRANSFORM_DCT:

+        if (do_inverse)

+            av_log(NULL, AV_LOG_INFO,"IDCT");

+        else

+            av_log(NULL, AV_LOG_INFO,"DCT");

+        ff_dct_init(d, fft_nbits, do_inverse);

+        break;

     }

     av_log(NULL, AV_LOG_INFO," %d test\n", fft_size);

@@ -321,6 +367,17 @@ int main(int argc, char **argv)

             tab_ref[0].im = tab_ref[fft_size_2].re;

             check_diff((float *)tab_ref, (float *)tab2, fft_size, 1.0);

         }

+        break;

+    case TRANSFORM_DCT:

+        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));

+        ff_dct_calc(d, tab);

+        if (do_inverse) {

+            idct_ref(tab_ref, tab1, fft_nbits);

+        } else {

+            dct_ref(tab_ref, tab1, fft_nbits);

+        }

+        check_diff((float *)tab_ref, (float *)tab, fft_size, 1.0);

+        break;

     }

     /* do a speed test */

@@ -351,6 +408,10 @@ int main(int argc, char **argv)

                     memcpy(tab2, tab1, fft_size * sizeof(FFTSample));

                     ff_rdft_calc(r, tab2);

                     break;

+                case TRANSFORM_DCT:

+                    memcpy(tab2, tab1, fft_size * sizeof(FFTSample));

+                    ff_dct_calc(d, tab2);

+                    break;

                 }

             }

             duration = gettime() - time_start;

@@ -374,6 +435,9 @@ int main(int argc, char **argv)

     case TRANSFORM_RDFT:

         ff_rdft_end(r);

         break;

+    case TRANSFORM_DCT:

+        ff_dct_end(d);

+        break;

     }

     return 0;

 }