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