new generic FFT/MDCT code for audio codecs
[libav.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., 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 int mdct_init(MDCTContext *s, int nbits, int inverse)
25 {
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 void imdct_calc(MDCTContext *s, FFTSample *output,
73 const FFTSample *input, FFTSample *tmp)
74 {
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 void mdct_calc(MDCTContext *s, FFTSample *out,
125 const FFTSample *input, FFTSample *tmp)
126 {
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 void mdct_end(MDCTContext *s)
166 {
167 av_freep(&s->tcos);
168 av_freep(&s->tsin);
169 fft_end(&s->fft);
170 }