Statistics
| Branch: | Revision:

ffmpeg / libavcodec / mdct.c @ a45fbda9

History | View | Annotate | Download (5.25 KB)

1
/*
2
 * MDCT/IMDCT transforms
3
 * Copyright (c) 2002 Fabrice Bellard
4
 *
5
 * This file is part of Libav.
6
 *
7
 * Libav is free software; you can redistribute it and/or
8
 * modify it under the terms of the GNU Lesser General Public
9
 * License as published by the Free Software Foundation; either
10
 * version 2.1 of the License, or (at your option) any later version.
11
 *
12
 * Libav is distributed in the hope that it will be useful,
13
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15
 * Lesser General Public License for more details.
16
 *
17
 * You should have received a copy of the GNU Lesser General Public
18
 * License along with Libav; if not, write to the Free Software
19
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20
 */
21

    
22
#include <stdlib.h>
23
#include <string.h>
24
#include "libavutil/common.h"
25
#include "libavutil/mathematics.h"
26
#include "fft.h"
27

    
28
/**
29
 * @file
30
 * MDCT/IMDCT transforms.
31
 */
32

    
33
#include "mdct_tablegen.h"
34

    
35
/**
36
 * init MDCT or IMDCT computation.
37
 */
38
av_cold int ff_mdct_init(FFTContext *s, int nbits, int inverse, double scale)
39
{
40
    int n, n4, i;
41
    double alpha, theta;
42
    int tstep;
43

    
44
    memset(s, 0, sizeof(*s));
45
    n = 1 << nbits;
46
    s->mdct_bits = nbits;
47
    s->mdct_size = n;
48
    n4 = n >> 2;
49
    s->mdct_permutation = FF_MDCT_PERM_NONE;
50

    
51
    if (ff_fft_init(s, s->mdct_bits - 2, inverse) < 0)
52
        goto fail;
53

    
54
    s->tcos = av_malloc(n/2 * sizeof(FFTSample));
55
    if (!s->tcos)
56
        goto fail;
57

    
58
    switch (s->mdct_permutation) {
59
    case FF_MDCT_PERM_NONE:
60
        s->tsin = s->tcos + n4;
61
        tstep = 1;
62
        break;
63
    case FF_MDCT_PERM_INTERLEAVE:
64
        s->tsin = s->tcos + 1;
65
        tstep = 2;
66
        break;
67
    default:
68
        goto fail;
69
    }
70

    
71
    theta = 1.0 / 8.0 + (scale < 0 ? n4 : 0);
72
    scale = sqrt(fabs(scale));
73
    for(i=0;i<n4;i++) {
74
        alpha = 2 * M_PI * (i + theta) / n;
75
        s->tcos[i*tstep] = -cos(alpha) * scale;
76
        s->tsin[i*tstep] = -sin(alpha) * scale;
77
    }
78
    return 0;
79
 fail:
80
    ff_mdct_end(s);
81
    return -1;
82
}
83

    
84
/* complex multiplication: p = a * b */
85
#define CMUL(pre, pim, are, aim, bre, bim) \
86
{\
87
    FFTSample _are = (are);\
88
    FFTSample _aim = (aim);\
89
    FFTSample _bre = (bre);\
90
    FFTSample _bim = (bim);\
91
    (pre) = _are * _bre - _aim * _bim;\
92
    (pim) = _are * _bim + _aim * _bre;\
93
}
94

    
95
/**
96
 * Compute the middle half of the inverse MDCT of size N = 2^nbits,
97
 * thus excluding the parts that can be derived by symmetry
98
 * @param output N/2 samples
99
 * @param input N/2 samples
100
 */
101
void ff_imdct_half_c(FFTContext *s, FFTSample *output, const FFTSample *input)
102
{
103
    int k, n8, n4, n2, n, j;
104
    const uint16_t *revtab = s->revtab;
105
    const FFTSample *tcos = s->tcos;
106
    const FFTSample *tsin = s->tsin;
107
    const FFTSample *in1, *in2;
108
    FFTComplex *z = (FFTComplex *)output;
109

    
110
    n = 1 << s->mdct_bits;
111
    n2 = n >> 1;
112
    n4 = n >> 2;
113
    n8 = n >> 3;
114

    
115
    /* pre rotation */
116
    in1 = input;
117
    in2 = input + n2 - 1;
118
    for(k = 0; k < n4; k++) {
119
        j=revtab[k];
120
        CMUL(z[j].re, z[j].im, *in2, *in1, tcos[k], tsin[k]);
121
        in1 += 2;
122
        in2 -= 2;
123
    }
124
    s->fft_calc(s, z);
125

    
126
    /* post rotation + reordering */
127
    for(k = 0; k < n8; k++) {
128
        FFTSample r0, i0, r1, i1;
129
        CMUL(r0, i1, z[n8-k-1].im, z[n8-k-1].re, tsin[n8-k-1], tcos[n8-k-1]);
130
        CMUL(r1, i0, z[n8+k  ].im, z[n8+k  ].re, tsin[n8+k  ], tcos[n8+k  ]);
131
        z[n8-k-1].re = r0;
132
        z[n8-k-1].im = i0;
133
        z[n8+k  ].re = r1;
134
        z[n8+k  ].im = i1;
135
    }
136
}
137

    
138
/**
139
 * Compute inverse MDCT of size N = 2^nbits
140
 * @param output N samples
141
 * @param input N/2 samples
142
 */
143
void ff_imdct_calc_c(FFTContext *s, FFTSample *output, const FFTSample *input)
144
{
145
    int k;
146
    int n = 1 << s->mdct_bits;
147
    int n2 = n >> 1;
148
    int n4 = n >> 2;
149

    
150
    ff_imdct_half_c(s, output+n4, input);
151

    
152
    for(k = 0; k < n4; k++) {
153
        output[k] = -output[n2-k-1];
154
        output[n-k-1] = output[n2+k];
155
    }
156
}
157

    
158
/**
159
 * Compute MDCT of size N = 2^nbits
160
 * @param input N samples
161
 * @param out N/2 samples
162
 */
163
void ff_mdct_calc_c(FFTContext *s, FFTSample *out, const FFTSample *input)
164
{
165
    int i, j, n, n8, n4, n2, n3;
166
    FFTSample re, im;
167
    const uint16_t *revtab = s->revtab;
168
    const FFTSample *tcos = s->tcos;
169
    const FFTSample *tsin = s->tsin;
170
    FFTComplex *x = (FFTComplex *)out;
171

    
172
    n = 1 << s->mdct_bits;
173
    n2 = n >> 1;
174
    n4 = n >> 2;
175
    n8 = n >> 3;
176
    n3 = 3 * n4;
177

    
178
    /* pre rotation */
179
    for(i=0;i<n8;i++) {
180
        re = -input[2*i+n3] - input[n3-1-2*i];
181
        im = -input[n4+2*i] + input[n4-1-2*i];
182
        j = revtab[i];
183
        CMUL(x[j].re, x[j].im, re, im, -tcos[i], tsin[i]);
184

    
185
        re = input[2*i] - input[n2-1-2*i];
186
        im = -(input[n2+2*i] + input[n-1-2*i]);
187
        j = revtab[n8 + i];
188
        CMUL(x[j].re, x[j].im, re, im, -tcos[n8 + i], tsin[n8 + i]);
189
    }
190

    
191
    s->fft_calc(s, x);
192

    
193
    /* post rotation */
194
    for(i=0;i<n8;i++) {
195
        FFTSample r0, i0, r1, i1;
196
        CMUL(i1, r0, x[n8-i-1].re, x[n8-i-1].im, -tsin[n8-i-1], -tcos[n8-i-1]);
197
        CMUL(i0, r1, x[n8+i  ].re, x[n8+i  ].im, -tsin[n8+i  ], -tcos[n8+i  ]);
198
        x[n8-i-1].re = r0;
199
        x[n8-i-1].im = i0;
200
        x[n8+i  ].re = r1;
201
        x[n8+i  ].im = i1;
202
    }
203
}
204

    
205
av_cold void ff_mdct_end(FFTContext *s)
206
{
207
    av_freep(&s->tcos);
208
    ff_fft_end(s);
209
}