Mercurial > libavcodec.hg
comparison mdct.c @ 781:6f5e87957bcb libavcodec
new generic FFT/MDCT code for audio codecs
author | bellard |
---|---|
date | Mon, 28 Oct 2002 00:34:08 +0000 |
parents | |
children | 670009af4fc2 |
comparison
equal
deleted
inserted
replaced
780:a48bb8bc63dd | 781:6f5e87957bcb |
---|---|
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 } |