libavcodec/dct.c
Go to the documentation of this file.
00001 /*
00002  * (I)DCT Transforms
00003  * Copyright (c) 2009 Peter Ross <pross@xvid.org>
00004  * Copyright (c) 2010 Alex Converse <alex.converse@gmail.com>
00005  * Copyright (c) 2010 Vitor Sessak
00006  *
00007  * This file is part of Libav.
00008  *
00009  * Libav is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  * Libav is distributed in the hope that it will be useful,
00015  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00016  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017  * Lesser General Public License for more details.
00018  *
00019  * You should have received a copy of the GNU Lesser General Public
00020  * License along with Libav; if not, write to the Free Software
00021  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00022  */
00023 
00030 #include <math.h>
00031 
00032 #include "libavutil/mathematics.h"
00033 #include "dct.h"
00034 #include "dct32.h"
00035 
00036 /* sin((M_PI * x / (2 * n)) */
00037 #define SIN(s, n, x) (s->costab[(n) - (x)])
00038 
00039 /* cos((M_PI * x / (2 * n)) */
00040 #define COS(s, n, x) (s->costab[x])
00041 
00042 static void ff_dst_calc_I_c(DCTContext *ctx, FFTSample *data)
00043 {
00044     int n = 1 << ctx->nbits;
00045     int i;
00046 
00047     data[0] = 0;
00048     for (i = 1; i < n / 2; i++) {
00049         float tmp1   = data[i    ];
00050         float tmp2   = data[n - i];
00051         float s      = SIN(ctx, n, 2 * i);
00052 
00053         s           *= tmp1 + tmp2;
00054         tmp1         = (tmp1 - tmp2) * 0.5f;
00055         data[i]      = s + tmp1;
00056         data[n - i]  = s - tmp1;
00057     }
00058 
00059     data[n / 2] *= 2;
00060     ctx->rdft.rdft_calc(&ctx->rdft, data);
00061 
00062     data[0] *= 0.5f;
00063 
00064     for (i = 1; i < n - 2; i += 2) {
00065         data[i + 1] +=  data[i - 1];
00066         data[i]      = -data[i + 2];
00067     }
00068 
00069     data[n - 1] = 0;
00070 }
00071 
00072 static void ff_dct_calc_I_c(DCTContext *ctx, FFTSample *data)
00073 {
00074     int n = 1 << ctx->nbits;
00075     int i;
00076     float next = -0.5f * (data[0] - data[n]);
00077 
00078     for (i = 0; i < n / 2; i++) {
00079         float tmp1 = data[i];
00080         float tmp2 = data[n - i];
00081         float s    = SIN(ctx, n, 2 * i);
00082         float c    = COS(ctx, n, 2 * i);
00083 
00084         c *= tmp1 - tmp2;
00085         s *= tmp1 - tmp2;
00086 
00087         next += c;
00088 
00089         tmp1        = (tmp1 + tmp2) * 0.5f;
00090         data[i]     = tmp1 - s;
00091         data[n - i] = tmp1 + s;
00092     }
00093 
00094     ctx->rdft.rdft_calc(&ctx->rdft, data);
00095     data[n] = data[1];
00096     data[1] = next;
00097 
00098     for (i = 3; i <= n; i += 2)
00099         data[i] = data[i - 2] - data[i];
00100 }
00101 
00102 static void ff_dct_calc_III_c(DCTContext *ctx, FFTSample *data)
00103 {
00104     int n = 1 << ctx->nbits;
00105     int i;
00106 
00107     float next  = data[n - 1];
00108     float inv_n = 1.0f / n;
00109 
00110     for (i = n - 2; i >= 2; i -= 2) {
00111         float val1 = data[i];
00112         float val2 = data[i - 1] - data[i + 1];
00113         float c    = COS(ctx, n, i);
00114         float s    = SIN(ctx, n, i);
00115 
00116         data[i]     = c * val1 + s * val2;
00117         data[i + 1] = s * val1 - c * val2;
00118     }
00119 
00120     data[1] = 2 * next;
00121 
00122     ctx->rdft.rdft_calc(&ctx->rdft, data);
00123 
00124     for (i = 0; i < n / 2; i++) {
00125         float tmp1 = data[i]         * inv_n;
00126         float tmp2 = data[n - i - 1] * inv_n;
00127         float csc  = ctx->csc2[i] * (tmp1 - tmp2);
00128 
00129         tmp1            += tmp2;
00130         data[i]          = tmp1 + csc;
00131         data[n - i - 1]  = tmp1 - csc;
00132     }
00133 }
00134 
00135 static void ff_dct_calc_II_c(DCTContext *ctx, FFTSample *data)
00136 {
00137     int n = 1 << ctx->nbits;
00138     int i;
00139     float next;
00140 
00141     for (i = 0; i < n / 2; i++) {
00142         float tmp1 = data[i];
00143         float tmp2 = data[n - i - 1];
00144         float s    = SIN(ctx, n, 2 * i + 1);
00145 
00146         s    *= tmp1 - tmp2;
00147         tmp1  = (tmp1 + tmp2) * 0.5f;
00148 
00149         data[i]     = tmp1 + s;
00150         data[n-i-1] = tmp1 - s;
00151     }
00152 
00153     ctx->rdft.rdft_calc(&ctx->rdft, data);
00154 
00155     next     = data[1] * 0.5;
00156     data[1] *= -1;
00157 
00158     for (i = n - 2; i >= 0; i -= 2) {
00159         float inr = data[i    ];
00160         float ini = data[i + 1];
00161         float c   = COS(ctx, n, i);
00162         float s   = SIN(ctx, n, i);
00163 
00164         data[i]     = c * inr + s * ini;
00165         data[i + 1] = next;
00166 
00167         next += s * inr - c * ini;
00168     }
00169 }
00170 
00171 static void dct32_func(DCTContext *ctx, FFTSample *data)
00172 {
00173     ctx->dct32(data, data);
00174 }
00175 
00176 av_cold int ff_dct_init(DCTContext *s, int nbits, enum DCTTransformType inverse)
00177 {
00178     int n = 1 << nbits;
00179     int i;
00180 
00181     memset(s, 0, sizeof(*s));
00182 
00183     s->nbits   = nbits;
00184     s->inverse = inverse;
00185 
00186     if (inverse == DCT_II && nbits == 5) {
00187         s->dct_calc = dct32_func;
00188     } else {
00189         ff_init_ff_cos_tabs(nbits + 2);
00190 
00191         s->costab = ff_cos_tabs[nbits + 2];
00192         s->csc2   = av_malloc(n / 2 * sizeof(FFTSample));
00193 
00194         if (ff_rdft_init(&s->rdft, nbits, inverse == DCT_III) < 0) {
00195             av_free(s->csc2);
00196             return -1;
00197         }
00198 
00199         for (i = 0; i < n / 2; i++)
00200             s->csc2[i] = 0.5 / sin((M_PI / (2 * n) * (2 * i + 1)));
00201 
00202         switch (inverse) {
00203         case DCT_I  : s->dct_calc = ff_dct_calc_I_c;   break;
00204         case DCT_II : s->dct_calc = ff_dct_calc_II_c;  break;
00205         case DCT_III: s->dct_calc = ff_dct_calc_III_c; break;
00206         case DST_I  : s->dct_calc = ff_dst_calc_I_c;   break;
00207         }
00208     }
00209 
00210     s->dct32 = ff_dct32_float;
00211     if (HAVE_MMX)
00212         ff_dct_init_mmx(s);
00213 
00214     return 0;
00215 }
00216 
00217 av_cold void ff_dct_end(DCTContext *s)
00218 {
00219     ff_rdft_end(&s->rdft);
00220     av_free(s->csc2);
00221 }