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
7 * This file is part of FFmpeg.
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.
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.
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
26 * (Inverse) Discrete Cosine Transforms. These are also known as the
27 * type II and type III DCTs respectively.
33 #include "libavutil/mathematics.h"
37 /* sin((M_PI * x / (2 * n)) */
38 #define SIN(s, n, x) (s->costab[(n) - (x)])
40 /* cos((M_PI * x / (2 * n)) */
41 #define COS(s, n, x) (s->costab[x])
43 static void dst_calc_I_c(DCTContext
*ctx
, FFTSample
*data
)
45 int n
= 1 << ctx
->nbits
;
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
);
55 tmp1
= (tmp1
- tmp2
) * 0.5f
;
57 data
[n
- i
] = s
- tmp1
;
61 ctx
->rdft
.rdft_calc(&ctx
->rdft
, data
);
65 for (i
= 1; i
< n
- 2; i
+= 2) {
66 data
[i
+ 1] += data
[i
- 1];
67 data
[i
] = -data
[i
+ 2];
73 static void dct_calc_I_c(DCTContext
*ctx
, FFTSample
*data
)
75 int n
= 1 << ctx
->nbits
;
77 float next
= -0.5f
* (data
[0] - data
[n
]);
79 for (i
= 0; i
< n
/ 2; i
++) {
81 float tmp2
= data
[n
- i
];
82 float s
= SIN(ctx
, n
, 2 * i
);
83 float c
= COS(ctx
, n
, 2 * i
);
90 tmp1
= (tmp1
+ tmp2
) * 0.5f
;
92 data
[n
- i
] = tmp1
+ s
;
95 ctx
->rdft
.rdft_calc(&ctx
->rdft
, data
);
99 for (i
= 3; i
<= n
; i
+= 2)
100 data
[i
] = data
[i
- 2] - data
[i
];
103 static void dct_calc_III_c(DCTContext
*ctx
, FFTSample
*data
)
105 int n
= 1 << ctx
->nbits
;
108 float next
= data
[n
- 1];
109 float inv_n
= 1.0f
/ n
;
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
);
117 data
[i
] = c
* val1
+ s
* val2
;
118 data
[i
+ 1] = s
* val1
- c
* val2
;
123 ctx
->rdft
.rdft_calc(&ctx
->rdft
, data
);
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
);
131 data
[i
] = tmp1
+ csc
;
132 data
[n
- i
- 1] = tmp1
- csc
;
136 static void dct_calc_II_c(DCTContext
*ctx
, FFTSample
*data
)
138 int n
= 1 << ctx
->nbits
;
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);
148 tmp1
= (tmp1
+ tmp2
) * 0.5f
;
151 data
[n
-i
-1] = tmp1
- s
;
154 ctx
->rdft
.rdft_calc(&ctx
->rdft
, data
);
156 next
= data
[1] * 0.5;
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
);
165 data
[i
] = c
* inr
+ s
* ini
;
168 next
+= s
* inr
- c
* ini
;
172 static void dct32_func(DCTContext
*ctx
, FFTSample
*data
)
174 ctx
->dct32(data
, data
);
177 av_cold
int ff_dct_init(DCTContext
*s
, int nbits
, enum DCTTransformType inverse
)
182 memset(s
, 0, sizeof(*s
));
185 s
->inverse
= inverse
;
187 if (inverse
== DCT_II
&& nbits
== 5) {
188 s
->dct_calc
= dct32_func
;
190 ff_init_ff_cos_tabs(nbits
+ 2);
192 s
->costab
= ff_cos_tabs
[nbits
+ 2];
193 s
->csc2
= av_malloc_array(n
/ 2, sizeof(FFTSample
));
195 if (ff_rdft_init(&s
->rdft
, nbits
, inverse
== DCT_III
) < 0) {
200 for (i
= 0; i
< n
/ 2; i
++)
201 s
->csc2
[i
] = 0.5 / sin((M_PI
/ (2 * n
) * (2 * i
+ 1)));
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;
211 s
->dct32
= ff_dct32_float
;
218 av_cold
void ff_dct_end(DCTContext
*s
)
220 ff_rdft_end(&s
->rdft
);