Statistics
| Branch: | Revision:

ffmpeg / libavcodec / mdct.c @ 18f77016

History | View | Annotate | Download (4.36 KB)

1 bb6f5690 Fabrice Bellard
/*
2
 * MDCT/IMDCT transforms
3
 * Copyright (c) 2002 Fabrice Bellard.
4
 *
5
 * This library is free software; you can redistribute it and/or
6
 * modify it under the terms of the GNU Lesser General Public
7
 * License as published by the Free Software Foundation; either
8
 * version 2 of the License, or (at your option) any later version.
9
 *
10
 * This library is distributed in the hope that it will be useful,
11
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13
 * Lesser General Public License for more details.
14
 *
15
 * You should have received a copy of the GNU Lesser General Public
16
 * License along with this library; if not, write to the Free Software
17
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
18
 */
19
#include "dsputil.h"
20
21
/*
22
 * init MDCT or IMDCT computation
23
 */
24 82696bee Fabrice Bellard
int ff_mdct_init(MDCTContext *s, int nbits, int inverse)
25 bb6f5690 Fabrice Bellard
{
26
    int n, n4, i;
27
    float alpha;
28
29
    memset(s, 0, sizeof(*s));
30
    n = 1 << nbits;
31
    s->nbits = nbits;
32
    s->n = n;
33
    n4 = n >> 2;
34
    s->tcos = malloc(n4 * sizeof(FFTSample));
35
    if (!s->tcos)
36
        goto fail;
37
    s->tsin = malloc(n4 * sizeof(FFTSample));
38
    if (!s->tsin)
39
        goto fail;
40
41
    for(i=0;i<n4;i++) {
42
        alpha = 2 * M_PI * (i + 1.0 / 8.0) / n;
43
        s->tcos[i] = -cos(alpha);
44
        s->tsin[i] = -sin(alpha);
45
    }
46
    if (fft_init(&s->fft, s->nbits - 2, inverse) < 0)
47
        goto fail;
48
    return 0;
49
 fail:
50
    av_freep(&s->tcos);
51
    av_freep(&s->tsin);
52
    return -1;
53
}
54
55
/* complex multiplication: p = a * b */
56
#define CMUL(pre, pim, are, aim, bre, bim) \
57
{\
58
    float _are = (are);\
59
    float _aim = (aim);\
60
    float _bre = (bre);\
61
    float _bim = (bim);\
62
    (pre) = _are * _bre - _aim * _bim;\
63
    (pim) = _are * _bim + _aim * _bre;\
64
}
65
66
/**
67
 * Compute inverse MDCT of size N = 2^nbits
68
 * @param output N samples
69
 * @param input N/2 samples
70
 * @param tmp N/2 samples
71
 */
72 82696bee Fabrice Bellard
void ff_imdct_calc(MDCTContext *s, FFTSample *output, 
73
                   const FFTSample *input, FFTSample *tmp)
74 bb6f5690 Fabrice Bellard
{
75
    int k, n8, n4, n2, n, j;
76
    const uint16_t *revtab = s->fft.revtab;
77
    const FFTSample *tcos = s->tcos;
78
    const FFTSample *tsin = s->tsin;
79
    const FFTSample *in1, *in2;
80
    FFTComplex *z = (FFTComplex *)tmp;
81
82
    n = 1 << s->nbits;
83
    n2 = n >> 1;
84
    n4 = n >> 2;
85
    n8 = n >> 3;
86
87
    /* pre rotation */
88
    in1 = input;
89
    in2 = input + n2 - 1;
90
    for(k = 0; k < n4; k++) {
91
        j=revtab[k];
92
        CMUL(z[j].re, z[j].im, *in2, *in1, tcos[k], tsin[k]);
93
        in1 += 2;
94
        in2 -= 2;
95
    }
96
    fft_calc(&s->fft, z);
97
98
    /* post rotation + reordering */
99
    /* XXX: optimize */
100
    for(k = 0; k < n4; k++) {
101
        CMUL(z[k].re, z[k].im, z[k].re, z[k].im, tcos[k], tsin[k]);
102
    }
103
    for(k = 0; k < n8; k++) {
104
        output[2*k] = -z[n8 + k].im;
105
        output[n2-1-2*k] = z[n8 + k].im;
106
107
        output[2*k+1] = z[n8-1-k].re;
108
        output[n2-1-2*k-1] = -z[n8-1-k].re;
109
110
        output[n2 + 2*k]=-z[k+n8].re;
111
        output[n-1- 2*k]=-z[k+n8].re;
112
113
        output[n2 + 2*k+1]=z[n8-k-1].im;
114
        output[n-2 - 2 * k] = z[n8-k-1].im;
115
    }
116
}
117
118
/**
119
 * Compute MDCT of size N = 2^nbits
120
 * @param input N samples
121
 * @param out N/2 samples
122
 * @param tmp temporary storage of N/2 samples
123
 */
124 82696bee Fabrice Bellard
void ff_mdct_calc(MDCTContext *s, FFTSample *out, 
125
                  const FFTSample *input, FFTSample *tmp)
126 bb6f5690 Fabrice Bellard
{
127
    int i, j, n, n8, n4, n2, n3;
128
    FFTSample re, im, re1, im1;
129
    const uint16_t *revtab = s->fft.revtab;
130
    const FFTSample *tcos = s->tcos;
131
    const FFTSample *tsin = s->tsin;
132
    FFTComplex *x = (FFTComplex *)tmp;
133
134
    n = 1 << s->nbits;
135
    n2 = n >> 1;
136
    n4 = n >> 2;
137
    n8 = n >> 3;
138
    n3 = 3 * n4;
139
140
    /* pre rotation */
141
    for(i=0;i<n8;i++) {
142
        re = -input[2*i+3*n4] - input[n3-1-2*i];
143
        im = -input[n4+2*i] + input[n4-1-2*i];
144
        j = revtab[i];
145
        CMUL(x[j].re, x[j].im, re, im, -tcos[i], tsin[i]);
146
147
        re = input[2*i] - input[n2-1-2*i];
148
        im = -(input[n2+2*i] + input[n-1-2*i]);
149
        j = revtab[n8 + i];
150
        CMUL(x[j].re, x[j].im, re, im, -tcos[n8 + i], tsin[n8 + i]);
151
    }
152
153
    fft_calc(&s->fft, x);
154
  
155
    /* post rotation */
156
    for(i=0;i<n4;i++) {
157
        re = x[i].re;
158
        im = x[i].im;
159
        CMUL(re1, im1, re, im, -tsin[i], -tcos[i]);
160
        out[2*i] = im1;
161
        out[n2-1-2*i] = re1;
162
    }
163
}
164
165 82696bee Fabrice Bellard
void ff_mdct_end(MDCTContext *s)
166 bb6f5690 Fabrice Bellard
{
167
    av_freep(&s->tcos);
168
    av_freep(&s->tsin);
169
    fft_end(&s->fft);
170
}