Mercurial > libavcodec.hg

changeset 10944:0985f1f7ab72 libavcodec

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Floating point discrete cosine transform
author vitor
date Wed, 20 Jan 2010 00:39:47 +0000
parents 3a723e8dcd26
children 0d93bbc17950
files Makefile dct.c dsputil.h fft-test.c
diffstat 4 files changed, 225 insertions(+), 1 deletions(-) [+]
line wrap: on
line diff
--- a/Makefile	Wed Jan 20 00:34:10 2010 +0000
+++ b/Makefile	Wed Jan 20 00:39:47 2010 +0000
@@ -27,6 +27,7 @@
 # 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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dct.c	Wed Jan 20 00:39:47 2010 +0000
@@ -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);
+}
--- a/dsputil.h	Wed Jan 20 00:34:10 2010 +0000
+++ b/dsputil.h	Wed Jan 20 00:39:47 2010 +0000
@@ -899,6 +899,26 @@
 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)\
--- a/fft-test.c	Wed Jan 20 00:34:10 2010 +0000
+++ b/fft-test.c	Wed Jan 20 00:39:47 2010 +0000
@@ -128,6 +128,39 @@
     }
 }
 
+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 @@
            "-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 @@
     TRANSFORM_FFT,
     TRANSFORM_MDCT,
     TRANSFORM_RDFT,
+    TRANSFORM_DCT,
 };
 
 int main(int argc, char **argv)
@@ -190,6 +225,7 @@
     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 @@
 
     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 @@
         case 'r':
             transform = TRANSFORM_RDFT;
             break;
+        case 'd':
+            transform = TRANSFORM_DCT;
+            break;
         case 'n':
             fft_nbits = atoi(optarg);
             break;
@@ -257,6 +296,13 @@
         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 @@
             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 @@
                     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 @@
     case TRANSFORM_RDFT:
         ff_rdft_end(r);
         break;
+    case TRANSFORM_DCT:
+        ff_dct_end(d);
+        break;
     }
     return 0;
 }