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 FFmpeg. 8 * 9 * FFmpeg 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 * FFmpeg 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 FFmpeg; 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#include <string.h> 32 33#include "libavutil/mathematics.h" 34#include "dct.h" 35#include "dct32.h" 36 37/* sin((M_PI * x / (2 * n)) */ 38#define SIN(s, n, x) (s->costab[(n) - (x)]) 39 40/* cos((M_PI * x / (2 * n)) */ 41#define COS(s, n, x) (s->costab[x]) 42 43static void dst_calc_I_c(DCTContext *ctx, FFTSample *data) 44{ 45 int n = 1 << ctx->nbits; 46 int i; 47 48 data[0] = 0; 49 for (i = 1; i < n / 2; i++) { 50 float tmp1 = data[i ]; 51 float tmp2 = data[n - i]; 52 float s = SIN(ctx, n, 2 * i); 53 54 s *= tmp1 + tmp2; 55 tmp1 = (tmp1 - tmp2) * 0.5f; 56 data[i] = s + tmp1; 57 data[n - i] = s - tmp1; 58 } 59 60 data[n / 2] *= 2; 61 ctx->rdft.rdft_calc(&ctx->rdft, data); 62 63 data[0] *= 0.5f; 64 65 for (i = 1; i < n - 2; i += 2) { 66 data[i + 1] += data[i - 1]; 67 data[i] = -data[i + 2]; 68 } 69 70 data[n - 1] = 0; 71} 72 73static void dct_calc_I_c(DCTContext *ctx, FFTSample *data) 74{ 75 int n = 1 << ctx->nbits; 76 int i; 77 float next = -0.5f * (data[0] - data[n]); 78 79 for (i = 0; i < n / 2; i++) { 80 float tmp1 = data[i]; 81 float tmp2 = data[n - i]; 82 float s = SIN(ctx, n, 2 * i); 83 float c = COS(ctx, n, 2 * i); 84 85 c *= tmp1 - tmp2; 86 s *= tmp1 - tmp2; 87 88 next += c; 89 90 tmp1 = (tmp1 + tmp2) * 0.5f; 91 data[i] = tmp1 - s; 92 data[n - i] = tmp1 + s; 93 } 94 95 ctx->rdft.rdft_calc(&ctx->rdft, data); 96 data[n] = data[1]; 97 data[1] = next; 98 99 for (i = 3; i <= n; i += 2) 100 data[i] = data[i - 2] - data[i]; 101} 102 103static void dct_calc_III_c(DCTContext *ctx, FFTSample *data) 104{ 105 int n = 1 << ctx->nbits; 106 int i; 107 108 float next = data[n - 1]; 109 float inv_n = 1.0f / n; 110 111 for (i = n - 2; i >= 2; i -= 2) { 112 float val1 = data[i]; 113 float val2 = data[i - 1] - data[i + 1]; 114 float c = COS(ctx, n, i); 115 float s = SIN(ctx, n, i); 116 117 data[i] = c * val1 + s * val2; 118 data[i + 1] = s * val1 - c * val2; 119 } 120 121 data[1] = 2 * next; 122 123 ctx->rdft.rdft_calc(&ctx->rdft, data); 124 125 for (i = 0; i < n / 2; i++) { 126 float tmp1 = data[i] * inv_n; 127 float tmp2 = data[n - i - 1] * inv_n; 128 float csc = ctx->csc2[i] * (tmp1 - tmp2); 129 130 tmp1 += tmp2; 131 data[i] = tmp1 + csc; 132 data[n - i - 1] = tmp1 - csc; 133 } 134} 135 136static void dct_calc_II_c(DCTContext *ctx, FFTSample *data) 137{ 138 int n = 1 << ctx->nbits; 139 int i; 140 float next; 141 142 for (i = 0; i < n / 2; i++) { 143 float tmp1 = data[i]; 144 float tmp2 = data[n - i - 1]; 145 float s = SIN(ctx, n, 2 * i + 1); 146 147 s *= tmp1 - tmp2; 148 tmp1 = (tmp1 + tmp2) * 0.5f; 149 150 data[i] = tmp1 + s; 151 data[n-i-1] = tmp1 - s; 152 } 153 154 ctx->rdft.rdft_calc(&ctx->rdft, data); 155 156 next = data[1] * 0.5; 157 data[1] *= -1; 158 159 for (i = n - 2; i >= 0; i -= 2) { 160 float inr = data[i ]; 161 float ini = data[i + 1]; 162 float c = COS(ctx, n, i); 163 float s = SIN(ctx, n, i); 164 165 data[i] = c * inr + s * ini; 166 data[i + 1] = next; 167 168 next += s * inr - c * ini; 169 } 170} 171 172static void dct32_func(DCTContext *ctx, FFTSample *data) 173{ 174 ctx->dct32(data, data); 175} 176 177av_cold int ff_dct_init(DCTContext *s, int nbits, enum DCTTransformType inverse) 178{ 179 int n = 1 << nbits; 180 int i; 181 182 memset(s, 0, sizeof(*s)); 183 184 s->nbits = nbits; 185 s->inverse = inverse; 186 187 if (inverse == DCT_II && nbits == 5) { 188 s->dct_calc = dct32_func; 189 } else { 190 ff_init_ff_cos_tabs(nbits + 2); 191 192 s->costab = ff_cos_tabs[nbits + 2]; 193 s->csc2 = av_malloc_array(n / 2, sizeof(FFTSample)); 194 195 if (ff_rdft_init(&s->rdft, nbits, inverse == DCT_III) < 0) { 196 av_free(s->csc2); 197 return -1; 198 } 199 200 for (i = 0; i < n / 2; i++) 201 s->csc2[i] = 0.5 / sin((M_PI / (2 * n) * (2 * i + 1))); 202 203 switch (inverse) { 204 case DCT_I : s->dct_calc = dct_calc_I_c; break; 205 case DCT_II : s->dct_calc = dct_calc_II_c; break; 206 case DCT_III: s->dct_calc = dct_calc_III_c; break; 207 case DST_I : s->dct_calc = dst_calc_I_c; break; 208 } 209 } 210 211 s->dct32 = ff_dct32_float; 212 if (ARCH_X86) 213 ff_dct_init_x86(s); 214 215 return 0; 216} 217 218av_cold void ff_dct_end(DCTContext *s) 219{ 220 ff_rdft_end(&s->rdft); 221 av_free(s->csc2); 222} 223