Update licensing information: The FSF changed postal address.
[ffmpeg.git] / libavcodec / mdct.c
1 /*
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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 #include "dsputil.h"
20
21 /**
22  * @file mdct.c
23  * MDCT/IMDCT transforms.
24  */
25
26 /**
27  * init MDCT or IMDCT computation.
28  */
29 int ff_mdct_init(MDCTContext *s, int nbits, int inverse)
30 {
31     int n, n4, i;
32     float alpha;
33
34     memset(s, 0, sizeof(*s));
35     n = 1 << nbits;
36     s->nbits = nbits;
37     s->n = n;
38     n4 = n >> 2;
39     s->tcos = av_malloc(n4 * sizeof(FFTSample));
40     if (!s->tcos)
41         goto fail;
42     s->tsin = av_malloc(n4 * sizeof(FFTSample));
43     if (!s->tsin)
44         goto fail;
45
46     for(i=0;i<n4;i++) {
47         alpha = 2 * M_PI * (i + 1.0 / 8.0) / n;
48         s->tcos[i] = -cos(alpha);
49         s->tsin[i] = -sin(alpha);
50     }
51     if (ff_fft_init(&s->fft, s->nbits - 2, inverse) < 0)
52         goto fail;
53     return 0;
54  fail:
55     av_freep(&s->tcos);
56     av_freep(&s->tsin);
57     return -1;
58 }
59
60 /* complex multiplication: p = a * b */
61 #define CMUL(pre, pim, are, aim, bre, bim) \
62 {\
63     float _are = (are);\
64     float _aim = (aim);\
65     float _bre = (bre);\
66     float _bim = (bim);\
67     (pre) = _are * _bre - _aim * _bim;\
68     (pim) = _are * _bim + _aim * _bre;\
69 }
70
71 /**
72  * Compute inverse MDCT of size N = 2^nbits
73  * @param output N samples
74  * @param input N/2 samples
75  * @param tmp N/2 samples
76  */
77 void ff_imdct_calc(MDCTContext *s, FFTSample *output,
78                    const FFTSample *input, FFTSample *tmp)
79 {
80     int k, n8, n4, n2, n, j;
81     const uint16_t *revtab = s->fft.revtab;
82     const FFTSample *tcos = s->tcos;
83     const FFTSample *tsin = s->tsin;
84     const FFTSample *in1, *in2;
85     FFTComplex *z = (FFTComplex *)tmp;
86
87     n = 1 << s->nbits;
88     n2 = n >> 1;
89     n4 = n >> 2;
90     n8 = n >> 3;
91
92     /* pre rotation */
93     in1 = input;
94     in2 = input + n2 - 1;
95     for(k = 0; k < n4; k++) {
96         j=revtab[k];
97         CMUL(z[j].re, z[j].im, *in2, *in1, tcos[k], tsin[k]);
98         in1 += 2;
99         in2 -= 2;
100     }
101     ff_fft_calc(&s->fft, z);
102
103     /* post rotation + reordering */
104     /* XXX: optimize */
105     for(k = 0; k < n4; k++) {
106         CMUL(z[k].re, z[k].im, z[k].re, z[k].im, tcos[k], tsin[k]);
107     }
108     for(k = 0; k < n8; k++) {
109         output[2*k] = -z[n8 + k].im;
110         output[n2-1-2*k] = z[n8 + k].im;
111
112         output[2*k+1] = z[n8-1-k].re;
113         output[n2-1-2*k-1] = -z[n8-1-k].re;
114
115         output[n2 + 2*k]=-z[k+n8].re;
116         output[n-1- 2*k]=-z[k+n8].re;
117
118         output[n2 + 2*k+1]=z[n8-k-1].im;
119         output[n-2 - 2 * k] = z[n8-k-1].im;
120     }
121 }
122
123 /**
124  * Compute MDCT of size N = 2^nbits
125  * @param input N samples
126  * @param out N/2 samples
127  * @param tmp temporary storage of N/2 samples
128  */
129 void ff_mdct_calc(MDCTContext *s, FFTSample *out,
130                   const FFTSample *input, FFTSample *tmp)
131 {
132     int i, j, n, n8, n4, n2, n3;
133     FFTSample re, im, re1, im1;
134     const uint16_t *revtab = s->fft.revtab;
135     const FFTSample *tcos = s->tcos;
136     const FFTSample *tsin = s->tsin;
137     FFTComplex *x = (FFTComplex *)tmp;
138
139     n = 1 << s->nbits;
140     n2 = n >> 1;
141     n4 = n >> 2;
142     n8 = n >> 3;
143     n3 = 3 * n4;
144
145     /* pre rotation */
146     for(i=0;i<n8;i++) {
147         re = -input[2*i+3*n4] - input[n3-1-2*i];
148         im = -input[n4+2*i] + input[n4-1-2*i];
149         j = revtab[i];
150         CMUL(x[j].re, x[j].im, re, im, -tcos[i], tsin[i]);
151
152         re = input[2*i] - input[n2-1-2*i];
153         im = -(input[n2+2*i] + input[n-1-2*i]);
154         j = revtab[n8 + i];
155         CMUL(x[j].re, x[j].im, re, im, -tcos[n8 + i], tsin[n8 + i]);
156     }
157
158     ff_fft_calc(&s->fft, x);
159
160     /* post rotation */
161     for(i=0;i<n4;i++) {
162         re = x[i].re;
163         im = x[i].im;
164         CMUL(re1, im1, re, im, -tsin[i], -tcos[i]);
165         out[2*i] = im1;
166         out[n2-1-2*i] = re1;
167     }
168 }
169
170 void ff_mdct_end(MDCTContext *s)
171 {
172     av_freep(&s->tcos);
173     av_freep(&s->tsin);
174     ff_fft_end(&s->fft);
175 }