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