blob: cca51eeaf8acf47bb9bfefd5f33634be275d66a6
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 | |
43 | static 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 | |
73 | static 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 | |
103 | static 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 | |
136 | static 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 | |
172 | static void dct32_func(DCTContext *ctx, FFTSample *data) |
173 | { |
174 | ctx->dct32(data, data); |
175 | } |
176 | |
177 | av_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 | if (!s->csc2) |
195 | return AVERROR(ENOMEM); |
196 | |
197 | if (ff_rdft_init(&s->rdft, nbits, inverse == DCT_III) < 0) { |
198 | av_freep(&s->csc2); |
199 | return -1; |
200 | } |
201 | |
202 | for (i = 0; i < n / 2; i++) |
203 | s->csc2[i] = 0.5 / sin((M_PI / (2 * n) * (2 * i + 1))); |
204 | |
205 | switch (inverse) { |
206 | case DCT_I : s->dct_calc = dct_calc_I_c; break; |
207 | case DCT_II : s->dct_calc = dct_calc_II_c; break; |
208 | case DCT_III: s->dct_calc = dct_calc_III_c; break; |
209 | case DST_I : s->dct_calc = dst_calc_I_c; break; |
210 | } |
211 | } |
212 | |
213 | s->dct32 = ff_dct32_float; |
214 | if (ARCH_X86) |
215 | ff_dct_init_x86(s); |
216 | |
217 | return 0; |
218 | } |
219 | |
220 | av_cold void ff_dct_end(DCTContext *s) |
221 | { |
222 | ff_rdft_end(&s->rdft); |
223 | av_freep(&s->csc2); |
224 | } |
225 |