ffmpeg / libavcodec / mdct.c @ 979395bb
History  View  Annotate  Download (5.94 KB)
1 
/*


2 
* MDCT/IMDCT transforms

3 
* Copyright (c) 2002 Fabrice Bellard

4 
*

5 
* This file is part of FFmpeg.

6 
*

7 
* FFmpeg 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 
* FFmpeg 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 FFmpeg; if not, write to the Free Software

19 
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 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 
// Generate a KaiserBessel Derived Window.

34 
#define BESSEL_I0_ITER 50 // default: 50 iterations of Bessel I0 approximation 
35 
av_cold void ff_kbd_window_init(float *window, float alpha, int n) 
36 
{ 
37 
int i, j;

38 
double sum = 0.0, bessel, tmp; 
39 
double local_window[FF_KBD_WINDOW_MAX];

40 
double alpha2 = (alpha * M_PI / n) * (alpha * M_PI / n);

41  
42 
assert(n <= FF_KBD_WINDOW_MAX); 
43  
44 
for (i = 0; i < n; i++) { 
45 
tmp = i * (n  i) * alpha2; 
46 
bessel = 1.0; 
47 
for (j = BESSEL_I0_ITER; j > 0; j) 
48 
bessel = bessel * tmp / (j * j) + 1;

49 
sum += bessel; 
50 
local_window[i] = sum; 
51 
} 
52  
53 
sum++; 
54 
for (i = 0; i < n; i++) 
55 
window[i] = sqrt(local_window[i] / sum); 
56 
} 
57  
58 
#include "mdct_tablegen.h" 
59  
60 
/**

61 
* init MDCT or IMDCT computation.

62 
*/

63 
av_cold int ff_mdct_init(FFTContext *s, int nbits, int inverse, double scale) 
64 
{ 
65 
int n, n4, i;

66 
double alpha, theta;

67 
int tstep;

68  
69 
memset(s, 0, sizeof(*s)); 
70 
n = 1 << nbits;

71 
s>mdct_bits = nbits; 
72 
s>mdct_size = n; 
73 
n4 = n >> 2;

74 
s>mdct_permutation = FF_MDCT_PERM_NONE; 
75  
76 
if (ff_fft_init(s, s>mdct_bits  2, inverse) < 0) 
77 
goto fail;

78  
79 
s>tcos = av_malloc(n/2 * sizeof(FFTSample)); 
80 
if (!s>tcos)

81 
goto fail;

82  
83 
switch (s>mdct_permutation) {

84 
case FF_MDCT_PERM_NONE:

85 
s>tsin = s>tcos + n4; 
86 
tstep = 1;

87 
break;

88 
case FF_MDCT_PERM_INTERLEAVE:

89 
s>tsin = s>tcos + 1;

90 
tstep = 2;

91 
break;

92 
default:

93 
goto fail;

94 
} 
95  
96 
theta = 1.0 / 8.0 + (scale < 0 ? n4 : 0); 
97 
scale = sqrt(fabs(scale)); 
98 
for(i=0;i<n4;i++) { 
99 
alpha = 2 * M_PI * (i + theta) / n;

100 
s>tcos[i*tstep] = cos(alpha) * scale; 
101 
s>tsin[i*tstep] = sin(alpha) * scale; 
102 
} 
103 
return 0; 
104 
fail:

105 
ff_mdct_end(s); 
106 
return 1; 
107 
} 
108  
109 
/* complex multiplication: p = a * b */

110 
#define CMUL(pre, pim, are, aim, bre, bim) \

111 
{\ 
112 
FFTSample _are = (are);\ 
113 
FFTSample _aim = (aim);\ 
114 
FFTSample _bre = (bre);\ 
115 
FFTSample _bim = (bim);\ 
116 
(pre) = _are * _bre  _aim * _bim;\ 
117 
(pim) = _are * _bim + _aim * _bre;\ 
118 
} 
119  
120 
/**

121 
* Compute the middle half of the inverse MDCT of size N = 2^nbits,

122 
* thus excluding the parts that can be derived by symmetry

123 
* @param output N/2 samples

124 
* @param input N/2 samples

125 
*/

126 
void ff_imdct_half_c(FFTContext *s, FFTSample *output, const FFTSample *input) 
127 
{ 
128 
int k, n8, n4, n2, n, j;

129 
const uint16_t *revtab = s>revtab;

130 
const FFTSample *tcos = s>tcos;

131 
const FFTSample *tsin = s>tsin;

132 
const FFTSample *in1, *in2;

133 
FFTComplex *z = (FFTComplex *)output; 
134  
135 
n = 1 << s>mdct_bits;

136 
n2 = n >> 1;

137 
n4 = n >> 2;

138 
n8 = n >> 3;

139  
140 
/* pre rotation */

141 
in1 = input; 
142 
in2 = input + n2  1;

143 
for(k = 0; k < n4; k++) { 
144 
j=revtab[k]; 
145 
CMUL(z[j].re, z[j].im, *in2, *in1, tcos[k], tsin[k]); 
146 
in1 += 2;

147 
in2 = 2;

148 
} 
149 
ff_fft_calc(s, z); 
150  
151 
/* post rotation + reordering */

152 
for(k = 0; k < n8; k++) { 
153 
FFTSample r0, i0, r1, i1; 
154 
CMUL(r0, i1, z[n8k1].im, z[n8k1].re, tsin[n8k1], tcos[n8k1]); 
155 
CMUL(r1, i0, z[n8+k ].im, z[n8+k ].re, tsin[n8+k ], tcos[n8+k ]); 
156 
z[n8k1].re = r0;

157 
z[n8k1].im = i0;

158 
z[n8+k ].re = r1; 
159 
z[n8+k ].im = i1; 
160 
} 
161 
} 
162  
163 
/**

164 
* Compute inverse MDCT of size N = 2^nbits

165 
* @param output N samples

166 
* @param input N/2 samples

167 
*/

168 
void ff_imdct_calc_c(FFTContext *s, FFTSample *output, const FFTSample *input) 
169 
{ 
170 
int k;

171 
int n = 1 << s>mdct_bits; 
172 
int n2 = n >> 1; 
173 
int n4 = n >> 2; 
174  
175 
ff_imdct_half_c(s, output+n4, input); 
176  
177 
for(k = 0; k < n4; k++) { 
178 
output[k] = output[n2k1];

179 
output[nk1] = output[n2+k];

180 
} 
181 
} 
182  
183 
/**

184 
* Compute MDCT of size N = 2^nbits

185 
* @param input N samples

186 
* @param out N/2 samples

187 
*/

188 
void ff_mdct_calc_c(FFTContext *s, FFTSample *out, const FFTSample *input) 
189 
{ 
190 
int i, j, n, n8, n4, n2, n3;

191 
FFTSample re, im; 
192 
const uint16_t *revtab = s>revtab;

193 
const FFTSample *tcos = s>tcos;

194 
const FFTSample *tsin = s>tsin;

195 
FFTComplex *x = (FFTComplex *)out; 
196  
197 
n = 1 << s>mdct_bits;

198 
n2 = n >> 1;

199 
n4 = n >> 2;

200 
n8 = n >> 3;

201 
n3 = 3 * n4;

202  
203 
/* pre rotation */

204 
for(i=0;i<n8;i++) { 
205 
re = input[2*i+n3]  input[n312*i]; 
206 
im = input[n4+2*i] + input[n412*i]; 
207 
j = revtab[i]; 
208 
CMUL(x[j].re, x[j].im, re, im, tcos[i], tsin[i]); 
209  
210 
re = input[2*i]  input[n212*i]; 
211 
im = (input[n2+2*i] + input[n12*i]); 
212 
j = revtab[n8 + i]; 
213 
CMUL(x[j].re, x[j].im, re, im, tcos[n8 + i], tsin[n8 + i]); 
214 
} 
215  
216 
ff_fft_calc(s, x); 
217  
218 
/* post rotation */

219 
for(i=0;i<n8;i++) { 
220 
FFTSample r0, i0, r1, i1; 
221 
CMUL(i1, r0, x[n8i1].re, x[n8i1].im, tsin[n8i1], tcos[n8i1]); 
222 
CMUL(i0, r1, x[n8+i ].re, x[n8+i ].im, tsin[n8+i ], tcos[n8+i ]); 
223 
x[n8i1].re = r0;

224 
x[n8i1].im = i0;

225 
x[n8+i ].re = r1; 
226 
x[n8+i ].im = i1; 
227 
} 
228 
} 
229  
230 
av_cold void ff_mdct_end(FFTContext *s)

231 
{ 
232 
av_freep(&s>tcos); 
233 
ff_fft_end(s); 
234 
} 