1/* 2 * (I)DCT Transforms 3 * Copyright (c) 2009 Peter Ross <pross@xvid.org> 4 * Copyright (c) 2010 Alex Converse <alex.converse@gmail.com> 5 * Copyright (c) 2010 Vitor Sessak 6 * 7 * This file is part of Libav. 8 * 9 * Libav is free software; you can redistribute it and/or 10 * modify it under the terms of the GNU Lesser General Public 11 * License as published by the Free Software Foundation; either 12 * version 2.1 of the License, or (at your option) any later version. 13 * 14 * Libav is distributed in the hope that it will be useful, 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 * Lesser General Public License for more details. 18 * 19 * You should have received a copy of the GNU Lesser General Public 20 * License along with Libav; if not, write to the Free Software 21 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 22 */ 23 24/** 25 * @file 26 * (Inverse) Discrete Cosine Transforms. These are also known as the 27 * type II and type III DCTs respectively. 28 */ 29 30#include <math.h> 31 32#include "libavutil/mathematics.h" 33#include "dct.h" 34#include "dct32.h" 35 36/* sin((M_PI * x / (2 * n)) */ 37#define SIN(s, n, x) (s->costab[(n) - (x)]) 38 39/* cos((M_PI * x / (2 * n)) */ 40#define COS(s, n, x) (s->costab[x]) 41 42static void ff_dst_calc_I_c(DCTContext *ctx, FFTSample *data) 43{ 44 int n = 1 << ctx->nbits; 45 int i; 46 47 data[0] = 0; 48 for (i = 1; i < n / 2; i++) { 49 float tmp1 = data[i ]; 50 float tmp2 = data[n - i]; 51 float s = SIN(ctx, n, 2 * i); 52 53 s *= tmp1 + tmp2; 54 tmp1 = (tmp1 - tmp2) * 0.5f; 55 data[i] = s + tmp1; 56 data[n - i] = s - tmp1; 57 } 58 59 data[n / 2] *= 2; 60 ctx->rdft.rdft_calc(&ctx->rdft, data); 61 62 data[0] *= 0.5f; 63 64 for (i = 1; i < n - 2; i += 2) { 65 data[i + 1] += data[i - 1]; 66 data[i] = -data[i + 2]; 67 } 68 69 data[n - 1] = 0; 70} 71 72static void ff_dct_calc_I_c(DCTContext *ctx, FFTSample *data) 73{ 74 int n = 1 << ctx->nbits; 75 int i; 76 float next = -0.5f * (data[0] - data[n]); 77 78 for (i = 0; i < n / 2; i++) { 79 float tmp1 = data[i]; 80 float tmp2 = data[n - i]; 81 float s = SIN(ctx, n, 2 * i); 82 float c = COS(ctx, n, 2 * i); 83 84 c *= tmp1 - tmp2; 85 s *= tmp1 - tmp2; 86 87 next += c; 88 89 tmp1 = (tmp1 + tmp2) * 0.5f; 90 data[i] = tmp1 - s; 91 data[n - i] = tmp1 + s; 92 } 93 94 ctx->rdft.rdft_calc(&ctx->rdft, data); 95 data[n] = data[1]; 96 data[1] = next; 97 98 for (i = 3; i <= n; i += 2) 99 data[i] = data[i - 2] - data[i]; 100} 101 102static void ff_dct_calc_III_c(DCTContext *ctx, FFTSample *data) 103{ 104 int n = 1 << ctx->nbits; 105 int i; 106 107 float next = data[n - 1]; 108 float inv_n = 1.0f / n; 109 110 for (i = n - 2; i >= 2; i -= 2) { 111 float val1 = data[i]; 112 float val2 = data[i - 1] - data[i + 1]; 113 float c = COS(ctx, n, i); 114 float s = SIN(ctx, n, i); 115 116 data[i] = c * val1 + s * val2; 117 data[i + 1] = s * val1 - c * val2; 118 } 119 120 data[1] = 2 * next; 121 122 ctx->rdft.rdft_calc(&ctx->rdft, data); 123 124 for (i = 0; i < n / 2; i++) { 125 float tmp1 = data[i] * inv_n; 126 float tmp2 = data[n - i - 1] * inv_n; 127 float csc = ctx->csc2[i] * (tmp1 - tmp2); 128 129 tmp1 += tmp2; 130 data[i] = tmp1 + csc; 131 data[n - i - 1] = tmp1 - csc; 132 } 133} 134 135static void ff_dct_calc_II_c(DCTContext *ctx, FFTSample *data) 136{ 137 int n = 1 << ctx->nbits; 138 int i; 139 float next; 140 141 for (i = 0; i < n / 2; i++) { 142 float tmp1 = data[i]; 143 float tmp2 = data[n - i - 1]; 144 float s = SIN(ctx, n, 2 * i + 1); 145 146 s *= tmp1 - tmp2; 147 tmp1 = (tmp1 + tmp2) * 0.5f; 148 149 data[i] = tmp1 + s; 150 data[n-i-1] = tmp1 - s; 151 } 152 153 ctx->rdft.rdft_calc(&ctx->rdft, data); 154 155 next = data[1] * 0.5; 156 data[1] *= -1; 157 158 for (i = n - 2; i >= 0; i -= 2) { 159 float inr = data[i ]; 160 float ini = data[i + 1]; 161 float c = COS(ctx, n, i); 162 float s = SIN(ctx, n, i); 163 164 data[i] = c * inr + s * ini; 165 data[i + 1] = next; 166 167 next += s * inr - c * ini; 168 } 169} 170 171static void dct32_func(DCTContext *ctx, FFTSample *data) 172{ 173 ctx->dct32(data, data); 174} 175 176av_cold int ff_dct_init(DCTContext *s, int nbits, enum DCTTransformType inverse) 177{ 178 int n = 1 << nbits; 179 int i; 180 181 memset(s, 0, sizeof(*s)); 182 183 s->nbits = nbits; 184 s->inverse = inverse; 185 186 if (inverse == DCT_II && nbits == 5) { 187 s->dct_calc = dct32_func; 188 } else { 189 ff_init_ff_cos_tabs(nbits + 2); 190 191 s->costab = ff_cos_tabs[nbits + 2]; 192 s->csc2 = av_malloc(n / 2 * sizeof(FFTSample)); 193 194 if (ff_rdft_init(&s->rdft, nbits, inverse == DCT_III) < 0) { 195 av_free(s->csc2); 196 return -1; 197 } 198 199 for (i = 0; i < n / 2; i++) 200 s->csc2[i] = 0.5 / sin((M_PI / (2 * n) * (2 * i + 1))); 201 202 switch (inverse) { 203 case DCT_I : s->dct_calc = ff_dct_calc_I_c; break; 204 case DCT_II : s->dct_calc = ff_dct_calc_II_c; break; 205 case DCT_III: s->dct_calc = ff_dct_calc_III_c; break; 206 case DST_I : s->dct_calc = ff_dst_calc_I_c; break; 207 } 208 } 209 210 s->dct32 = ff_dct32_float; 211 if (HAVE_MMX) 212 ff_dct_init_mmx(s); 213 214 return 0; 215} 216 217av_cold void ff_dct_end(DCTContext *s) 218{ 219 ff_rdft_end(&s->rdft); 220 av_free(s->csc2); 221} 222