ffmpeg / libavcodec / ac3enc_fixed.c @ 6fd96d1a
History  View  Annotate  Download (10.9 KB)
1 
/*


2 
* The simplest AC3 encoder

3 
* Copyright (c) 2000 Fabrice Bellard

4 
* Copyright (c) 20062010 Justin Ruggles <justin.ruggles@gmail.com>

5 
* Copyright (c) 20062010 Prakash Punnoor <prakash@punnoor.de>

6 
*

7 
* This file is part of FFmpeg.

8 
*

9 
* FFmpeg is free software; you can redistribute it and/or

10 
* modify it under the terms of the GNU Lesser General Public

11 
* License as published by the Free Software Foundation; either

12 
* version 2.1 of the License, or (at your option) any later version.

13 
*

14 
* FFmpeg is distributed in the hope that it will be useful,

15 
* but WITHOUT ANY WARRANTY; without even the implied warranty of

16 
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU

17 
* Lesser General Public License for more details.

18 
*

19 
* You should have received a copy of the GNU Lesser General Public

20 
* License along with FFmpeg; if not, write to the Free Software

21 
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA

22 
*/

23  
24 
/**

25 
* @file

26 
* fixedpoint AC3 encoder.

27 
*/

28  
29 
#undef CONFIG_AC3ENC_FLOAT

30 
#include "ac3enc.c" 
31  
32  
33 
/** Scale a float value by 2^15, convert to an integer, and clip to range 32767..32767. */

34 
#define FIX15(a) av_clip(SCALE_FLOAT(a, 15), 32767, 32767) 
35  
36  
37 
/**

38 
* Finalize MDCT and free allocated memory.

39 
*/

40 
static av_cold void mdct_end(AC3MDCTContext *mdct) 
41 
{ 
42 
mdct>nbits = 0;

43 
av_freep(&mdct>costab); 
44 
av_freep(&mdct>sintab); 
45 
av_freep(&mdct>xcos1); 
46 
av_freep(&mdct>xsin1); 
47 
av_freep(&mdct>rot_tmp); 
48 
av_freep(&mdct>cplx_tmp); 
49 
} 
50  
51  
52 
/**

53 
* Initialize FFT tables.

54 
* @param ln log2(FFT size)

55 
*/

56 
static av_cold int fft_init(AVCodecContext *avctx, AC3MDCTContext *mdct, int ln) 
57 
{ 
58 
int i, n, n2;

59 
float alpha;

60  
61 
n = 1 << ln;

62 
n2 = n >> 1;

63  
64 
FF_ALLOC_OR_GOTO(avctx, mdct>costab, n2 * sizeof(*mdct>costab), fft_alloc_fail);

65 
FF_ALLOC_OR_GOTO(avctx, mdct>sintab, n2 * sizeof(*mdct>sintab), fft_alloc_fail);

66  
67 
for (i = 0; i < n2; i++) { 
68 
alpha = 2.0 * M_PI * i / n; 
69 
mdct>costab[i] = FIX15(cos(alpha)); 
70 
mdct>sintab[i] = FIX15(sin(alpha)); 
71 
} 
72  
73 
return 0; 
74 
fft_alloc_fail:

75 
mdct_end(mdct); 
76 
return AVERROR(ENOMEM);

77 
} 
78  
79  
80 
/**

81 
* Initialize MDCT tables.

82 
* @param nbits log2(MDCT size)

83 
*/

84 
static av_cold int mdct_init(AVCodecContext *avctx, AC3MDCTContext *mdct, 
85 
int nbits)

86 
{ 
87 
int i, n, n4, ret;

88  
89 
n = 1 << nbits;

90 
n4 = n >> 2;

91  
92 
mdct>nbits = nbits; 
93  
94 
ret = fft_init(avctx, mdct, nbits  2);

95 
if (ret)

96 
return ret;

97  
98 
mdct>window = ff_ac3_window; 
99  
100 
FF_ALLOC_OR_GOTO(avctx, mdct>xcos1, n4 * sizeof(*mdct>xcos1), mdct_alloc_fail);

101 
FF_ALLOC_OR_GOTO(avctx, mdct>xsin1, n4 * sizeof(*mdct>xsin1), mdct_alloc_fail);

102 
FF_ALLOC_OR_GOTO(avctx, mdct>rot_tmp, n * sizeof(*mdct>rot_tmp), mdct_alloc_fail);

103 
FF_ALLOC_OR_GOTO(avctx, mdct>cplx_tmp, n4 * sizeof(*mdct>cplx_tmp), mdct_alloc_fail);

104  
105 
for (i = 0; i < n4; i++) { 
106 
float alpha = 2.0 * M_PI * (i + 1.0 / 8.0) / n; 
107 
mdct>xcos1[i] = FIX15(cos(alpha)); 
108 
mdct>xsin1[i] = FIX15(sin(alpha)); 
109 
} 
110  
111 
return 0; 
112 
mdct_alloc_fail:

113 
mdct_end(mdct); 
114 
return AVERROR(ENOMEM);

115 
} 
116  
117  
118 
/** Butterfly op */

119 
#define BF(pre, pim, qre, qim, pre1, pim1, qre1, qim1) \

120 
{ \ 
121 
int ax, ay, bx, by; \

122 
bx = pre1; \ 
123 
by = pim1; \ 
124 
ax = qre1; \ 
125 
ay = qim1; \ 
126 
pre = (bx + ax) >> 1; \

127 
pim = (by + ay) >> 1; \

128 
qre = (bx  ax) >> 1; \

129 
qim = (by  ay) >> 1; \

130 
} 
131  
132  
133 
/** Complex multiply */

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

135 
{ \ 
136 
pre = (MUL16(are, bre)  MUL16(aim, bim)) >> 15; \

137 
pim = (MUL16(are, bim) + MUL16(bre, aim)) >> 15; \

138 
} 
139  
140  
141 
/**

142 
* Calculate a 2^n point complex FFT on 2^ln points.

143 
* @param z complex input/output samples

144 
* @param ln log2(FFT size)

145 
*/

146 
static void fft(AC3MDCTContext *mdct, IComplex *z, int ln) 
147 
{ 
148 
int j, l, np, np2;

149 
int nblocks, nloops;

150 
register IComplex *p,*q;

151 
int tmp_re, tmp_im;

152  
153 
np = 1 << ln;

154  
155 
/* reverse */

156 
for (j = 0; j < np; j++) { 
157 
int k = av_reverse[j] >> (8  ln); 
158 
if (k < j)

159 
FFSWAP(IComplex, z[k], z[j]); 
160 
} 
161  
162 
/* pass 0 */

163  
164 
p = &z[0];

165 
j = np >> 1;

166 
do {

167 
BF(p[0].re, p[0].im, p[1].re, p[1].im, 
168 
p[0].re, p[0].im, p[1].re, p[1].im); 
169 
p += 2;

170 
} while (j);

171  
172 
/* pass 1 */

173  
174 
p = &z[0];

175 
j = np >> 2;

176 
do {

177 
BF(p[0].re, p[0].im, p[2].re, p[2].im, 
178 
p[0].re, p[0].im, p[2].re, p[2].im); 
179 
BF(p[1].re, p[1].im, p[3].re, p[3].im, 
180 
p[1].re, p[1].im, p[3].im, p[3].re); 
181 
p+=4;

182 
} while (j);

183  
184 
/* pass 2 .. ln1 */

185  
186 
nblocks = np >> 3;

187 
nloops = 1 << 2; 
188 
np2 = np >> 1;

189 
do {

190 
p = z; 
191 
q = z + nloops; 
192 
for (j = 0; j < nblocks; j++) { 
193 
BF(p>re, p>im, q>re, q>im, 
194 
p>re, p>im, q>re, q>im); 
195 
p++; 
196 
q++; 
197 
for(l = nblocks; l < np2; l += nblocks) {

198 
CMUL(tmp_re, tmp_im, mdct>costab[l], mdct>sintab[l], q>re, q>im); 
199 
BF(p>re, p>im, q>re, q>im, 
200 
p>re, p>im, tmp_re, tmp_im); 
201 
p++; 
202 
q++; 
203 
} 
204 
p += nloops; 
205 
q += nloops; 
206 
} 
207 
nblocks = nblocks >> 1;

208 
nloops = nloops << 1;

209 
} while (nblocks);

210 
} 
211  
212  
213 
/**

214 
* Calculate a 512point MDCT

215 
* @param out 256 output frequency coefficients

216 
* @param in 512 windowed input audio samples

217 
*/

218 
static void mdct512(AC3MDCTContext *mdct, int32_t *out, int16_t *in) 
219 
{ 
220 
int i, re, im, n, n2, n4;

221 
int16_t *rot = mdct>rot_tmp; 
222 
IComplex *x = mdct>cplx_tmp; 
223  
224 
n = 1 << mdct>nbits;

225 
n2 = n >> 1;

226 
n4 = n >> 2;

227  
228 
/* shift to simplify computations */

229 
for (i = 0; i <n4; i++) 
230 
rot[i] = in[i + 3*n4];

231 
memcpy(&rot[n4], &in[0], 3*n4*sizeof(*in)); 
232  
233 
/* pre rotation */

234 
for (i = 0; i < n4; i++) { 
235 
re = ((int)rot[ 2*i]  (int)rot[ n12*i]) >> 1; 
236 
im = ((int)rot[n2+2*i]  (int)rot[n212*i]) >> 1; 
237 
CMUL(x[i].re, x[i].im, re, im, mdct>xcos1[i], mdct>xsin1[i]); 
238 
} 
239  
240 
fft(mdct, x, mdct>nbits  2);

241  
242 
/* post rotation */

243 
for (i = 0; i < n4; i++) { 
244 
re = x[i].re; 
245 
im = x[i].im; 
246 
CMUL(out[n212*i], out[2*i], re, im, mdct>xsin1[i], mdct>xcos1[i]); 
247 
} 
248 
} 
249  
250  
251 
/**

252 
* Apply KBD window to input samples prior to MDCT.

253 
*/

254 
static void apply_window(int16_t *output, const int16_t *input, 
255 
const int16_t *window, int n) 
256 
{ 
257 
int i;

258 
int n2 = n >> 1; 
259  
260 
for (i = 0; i < n2; i++) { 
261 
output[i] = MUL16(input[i], window[i]) >> 15;

262 
output[ni1] = MUL16(input[ni1], window[i]) >> 15; 
263 
} 
264 
} 
265  
266  
267 
/**

268 
* Calculate the log2() of the maximum absolute value in an array.

269 
* @param tab input array

270 
* @param n number of values in the array

271 
* @return log2(max(abs(tab[])))

272 
*/

273 
static int log2_tab(int16_t *tab, int n) 
274 
{ 
275 
int i, v;

276  
277 
v = 0;

278 
for (i = 0; i < n; i++) 
279 
v = abs(tab[i]); 
280  
281 
return av_log2(v);

282 
} 
283  
284  
285 
/**

286 
* Leftshift each value in an array by a specified amount.

287 
* @param tab input array

288 
* @param n number of values in the array

289 
* @param lshift left shift amount. a negative value means right shift.

290 
*/

291 
static void lshift_tab(int16_t *tab, int n, int lshift) 
292 
{ 
293 
int i;

294  
295 
if (lshift > 0) { 
296 
for (i = 0; i < n; i++) 
297 
tab[i] <<= lshift; 
298 
} else if (lshift < 0) { 
299 
lshift = lshift; 
300 
for (i = 0; i < n; i++) 
301 
tab[i] >>= lshift; 
302 
} 
303 
} 
304  
305  
306 
/**

307 
* Normalize the input samples to use the maximum available precision.

308 
* This assumes signed 16bit input samples. Exponents are reduced by 9 to

309 
* match the 24bit internal precision for MDCT coefficients.

310 
*

311 
* @return exponent shift

312 
*/

313 
static int normalize_samples(AC3EncodeContext *s) 
314 
{ 
315 
int v = 14  log2_tab(s>windowed_samples, AC3_WINDOW_SIZE); 
316 
v = FFMAX(0, v);

317 
lshift_tab(s>windowed_samples, AC3_WINDOW_SIZE, v); 
318 
return v  9; 
319 
} 
320  
321  
322 
#ifdef TEST

323 
/*************************************************************************/

324 
/* TEST */

325  
326 
#include "libavutil/lfg.h" 
327  
328 
#define MDCT_NBITS 9 
329 
#define MDCT_SAMPLES (1 << MDCT_NBITS) 
330 
#define FN (MDCT_SAMPLES/4) 
331  
332  
333 
static void fft_test(AC3MDCTContext *mdct, AVLFG *lfg) 
334 
{ 
335 
IComplex in[FN], in1[FN]; 
336 
int k, n, i;

337 
float sum_re, sum_im, a;

338  
339 
for (i = 0; i < FN; i++) { 
340 
in[i].re = av_lfg_get(lfg) % 65535  32767; 
341 
in[i].im = av_lfg_get(lfg) % 65535  32767; 
342 
in1[i] = in[i]; 
343 
} 
344 
fft(mdct, in, 7);

345  
346 
/* do it by hand */

347 
for (k = 0; k < FN; k++) { 
348 
sum_re = 0;

349 
sum_im = 0;

350 
for (n = 0; n < FN; n++) { 
351 
a = 2 * M_PI * (n * k) / FN;

352 
sum_re += in1[n].re * cos(a)  in1[n].im * sin(a); 
353 
sum_im += in1[n].re * sin(a) + in1[n].im * cos(a); 
354 
} 
355 
av_log(NULL, AV_LOG_DEBUG, "%3d: %6d,%6d %6.0f,%6.0f\n", 
356 
k, in[k].re, in[k].im, sum_re / FN, sum_im / FN); 
357 
} 
358 
} 
359  
360  
361 
static void mdct_test(AC3MDCTContext *mdct, AVLFG *lfg) 
362 
{ 
363 
int16_t input[MDCT_SAMPLES]; 
364 
int32_t output[AC3_MAX_COEFS]; 
365 
float input1[MDCT_SAMPLES];

366 
float output1[AC3_MAX_COEFS];

367 
float s, a, err, e, emax;

368 
int i, k, n;

369  
370 
for (i = 0; i < MDCT_SAMPLES; i++) { 
371 
input[i] = (av_lfg_get(lfg) % 65535  32767) * 9 / 10; 
372 
input1[i] = input[i]; 
373 
} 
374  
375 
mdct512(mdct, output, input); 
376  
377 
/* do it by hand */

378 
for (k = 0; k < AC3_MAX_COEFS; k++) { 
379 
s = 0;

380 
for (n = 0; n < MDCT_SAMPLES; n++) { 
381 
a = (2*M_PI*(2*n+1+MDCT_SAMPLES/2)*(2*k+1) / (4 * MDCT_SAMPLES)); 
382 
s += input1[n] * cos(a); 
383 
} 
384 
output1[k] = 2 * s / MDCT_SAMPLES;

385 
} 
386  
387 
err = 0;

388 
emax = 0;

389 
for (i = 0; i < AC3_MAX_COEFS; i++) { 
390 
av_log(NULL, AV_LOG_DEBUG, "%3d: %7d %7.0f\n", i, output[i], output1[i]); 
391 
e = output[i]  output1[i]; 
392 
if (e > emax)

393 
emax = e; 
394 
err += e * e; 
395 
} 
396 
av_log(NULL, AV_LOG_DEBUG, "err2=%f emax=%f\n", err / AC3_MAX_COEFS, emax); 
397 
} 
398  
399  
400 
int main(void) 
401 
{ 
402 
AVLFG lfg; 
403 
AC3MDCTContext mdct; 
404  
405 
mdct.avctx = NULL;

406 
av_log_set_level(AV_LOG_DEBUG); 
407 
mdct_init(&mdct, 9);

408  
409 
fft_test(&mdct, &lfg); 
410 
mdct_test(&mdct, &lfg); 
411  
412 
return 0; 
413 
} 
414 
#endif /* TEST */ 
415  
416  
417 
AVCodec ac3_fixed_encoder = { 
418 
"ac3_fixed",

419 
AVMEDIA_TYPE_AUDIO, 
420 
CODEC_ID_AC3, 
421 
sizeof(AC3EncodeContext),

422 
ac3_encode_init, 
423 
ac3_encode_frame, 
424 
ac3_encode_close, 
425 
NULL,

426 
.sample_fmts = (const enum AVSampleFormat[]){AV_SAMPLE_FMT_S16,AV_SAMPLE_FMT_NONE}, 
427 
.long_name = NULL_IF_CONFIG_SMALL("ATSC A/52A (AC3)"),

428 
.channel_layouts = ac3_channel_layouts, 
429 
}; 