Floating point discrete cosine transform

Originally committed as revision 21338 to svn://svn.ffmpeg.org/ffmpeg/trunk

4 files changed, 225 insertions, 1 deletions

diff --git a/libavcodec/Makefile b/libavcodec/Makefile
index 96a703b57f..3888589511 100644
--- a/libavcodec/Makefile
+++ b/libavcodec/Makefile
@@ -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
diff --git a/libavcodec/dct.c b/libavcodec/dct.c
new file mode 100644
index 0000000000..3776349df9
--- /dev/null
+++ b/libavcodec/dct.c
@@ -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);
+}
diff --git a/libavcodec/dsputil.h b/libavcodec/dsputil.h
index 1476ae7ba3..65088c11ca 100644
--- a/libavcodec/dsputil.h
+++ b/libavcodec/dsputil.h
@@ -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)\
diff --git a/libavcodec/fft-test.c b/libavcodec/fft-test.c
index ed1b370a40..f07ef032f8 100644
--- a/libavcodec/fft-test.c
+++ b/libavcodec/fft-test.c
@@ -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;
 }