Statistics
| Branch: | Revision:

ffmpeg / libavcodec / dct.c @ 2912e87a

History | View | Annotate | Download (5.22 KB)

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 Libav.
8
 *
9
 * Libav 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
 * Libav 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 Libav; 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
#include "x86/fft.h"
34

    
35
#define DCT32_FLOAT
36
#include "dct32.c"
37

    
38
/* sin((M_PI * x / (2*n)) */
39
#define SIN(s,n,x) (s->costab[(n) - (x)])
40

    
41
/* cos((M_PI * x / (2*n)) */
42
#define COS(s,n,x) (s->costab[x])
43

    
44
static void ff_dst_calc_I_c(DCTContext *ctx, FFTSample *data)
45
{
46
    int n = 1 << ctx->nbits;
47
    int i;
48

    
49
    data[0] = 0;
50
    for(i = 1; i < n/2; i++) {
51
        float tmp1 = data[i    ];
52
        float tmp2 = data[n - i];
53
        float s = SIN(ctx, n, 2*i);
54

    
55
        s *= tmp1 + tmp2;
56
        tmp1 = (tmp1 - tmp2) * 0.5f;
57
        data[i    ] = s + tmp1;
58
        data[n - i] = s - tmp1;
59
    }
60

    
61
    data[n/2] *= 2;
62
    ff_rdft_calc(&ctx->rdft, data);
63

    
64
    data[0] *= 0.5f;
65

    
66
    for(i = 1; i < n-2; i += 2) {
67
        data[i + 1] += data[i - 1];
68
        data[i    ] = -data[i + 2];
69
    }
70

    
71
    data[n-1] = 0;
72
}
73

    
74
static void ff_dct_calc_I_c(DCTContext *ctx, FFTSample *data)
75
{
76
    int n = 1 << ctx->nbits;
77
    int i;
78
    float next = -0.5f * (data[0] - data[n]);
79

    
80
    for(i = 0; i < n/2; i++) {
81
        float tmp1 = data[i    ];
82
        float tmp2 = data[n - i];
83
        float s = SIN(ctx, n, 2*i);
84
        float c = COS(ctx, n, 2*i);
85

    
86
        c *= tmp1 - tmp2;
87
        s *= tmp1 - tmp2;
88

    
89
        next += c;
90

    
91
        tmp1 = (tmp1 + tmp2) * 0.5f;
92
        data[i    ] = tmp1 - s;
93
        data[n - i] = tmp1 + s;
94
    }
95

    
96
    ff_rdft_calc(&ctx->rdft, data);
97
    data[n] = data[1];
98
    data[1] = next;
99

    
100
    for(i = 3; i <= n; i += 2)
101
        data[i] = data[i - 2] - data[i];
102
}
103

    
104
static void ff_dct_calc_III_c(DCTContext *ctx, FFTSample *data)
105
{
106
    int n = 1 << ctx->nbits;
107
    int i;
108

    
109
    float next = data[n - 1];
110
    float inv_n = 1.0f / n;
111

    
112
    for (i = n - 2; i >= 2; i -= 2) {
113
        float val1 = data[i    ];
114
        float val2 = data[i - 1] - data[i + 1];
115
        float c = COS(ctx, n, i);
116
        float s = SIN(ctx, n, i);
117

    
118
        data[i    ] = c * val1 + s * val2;
119
        data[i + 1] = s * val1 - c * val2;
120
    }
121

    
122
    data[1] = 2 * next;
123

    
124
    ff_rdft_calc(&ctx->rdft, data);
125

    
126
    for (i = 0; i < n / 2; i++) {
127
        float tmp1 = data[i        ] * inv_n;
128
        float tmp2 = data[n - i - 1] * inv_n;
129
        float csc = ctx->csc2[i] * (tmp1 - tmp2);
130

    
131
        tmp1 += tmp2;
132
        data[i        ] = tmp1 + csc;
133
        data[n - i - 1] = tmp1 - csc;
134
    }
135
}
136

    
137
static void ff_dct_calc_II_c(DCTContext *ctx, FFTSample *data)
138
{
139
    int n = 1 << ctx->nbits;
140
    int i;
141
    float next;
142

    
143
    for (i=0; i < n/2; i++) {
144
        float tmp1 = data[i        ];
145
        float tmp2 = data[n - i - 1];
146
        float s = SIN(ctx, n, 2*i + 1);
147

    
148
        s *= tmp1 - tmp2;
149
        tmp1 = (tmp1 + tmp2) * 0.5f;
150

    
151
        data[i    ] = tmp1 + s;
152
        data[n-i-1] = tmp1 - s;
153
    }
154

    
155
    ff_rdft_calc(&ctx->rdft, data);
156

    
157
    next = data[1] * 0.5;
158
    data[1] *= -1;
159

    
160
    for (i = n - 2; i >= 0; i -= 2) {
161
        float inr = data[i    ];
162
        float ini = data[i + 1];
163
        float c = COS(ctx, n, i);
164
        float s = SIN(ctx, n, i);
165

    
166
        data[i  ] = c * inr + s * ini;
167

    
168
        data[i+1] = next;
169

    
170
        next +=     s * inr - c * ini;
171
    }
172
}
173

    
174
static void dct32_func(DCTContext *ctx, FFTSample *data)
175
{
176
    ctx->dct32(data, data);
177
}
178

    
179
void ff_dct_calc(DCTContext *s, FFTSample *data)
180
{
181
    s->dct_calc(s, data);
182
}
183

    
184
av_cold int ff_dct_init(DCTContext *s, int nbits, enum DCTTransformType inverse)
185
{
186
    int n = 1 << nbits;
187
    int i;
188

    
189
    s->nbits    = nbits;
190
    s->inverse  = inverse;
191

    
192
    ff_init_ff_cos_tabs(nbits+2);
193

    
194
    s->costab = ff_cos_tabs[nbits+2];
195

    
196
    s->csc2 = av_malloc(n/2 * sizeof(FFTSample));
197

    
198
    if (ff_rdft_init(&s->rdft, nbits, inverse == DCT_III) < 0) {
199
        av_free(s->csc2);
200
        return -1;
201
    }
202

    
203
    for (i = 0; i < n/2; i++)
204
        s->csc2[i] = 0.5 / sin((M_PI / (2*n) * (2*i + 1)));
205

    
206
    switch(inverse) {
207
    case DCT_I  : s->dct_calc = ff_dct_calc_I_c; break;
208
    case DCT_II : s->dct_calc = ff_dct_calc_II_c ; break;
209
    case DCT_III: s->dct_calc = ff_dct_calc_III_c; break;
210
    case DST_I  : s->dct_calc = ff_dst_calc_I_c; break;
211
    }
212

    
213
    if (inverse == DCT_II && nbits == 5)
214
        s->dct_calc = dct32_func;
215

    
216
    s->dct32 = dct32;
217
    if (HAVE_MMX)     ff_dct_init_mmx(s);
218

    
219
    return 0;
220
}
221

    
222
av_cold void ff_dct_end(DCTContext *s)
223
{
224
    ff_rdft_end(&s->rdft);
225
    av_free(s->csc2);
226
}