ffmpeg / libavcodec / ac3enc_fixed.c @ 7e0a284b
History  View  Annotate  Download (11.5 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(DSPContext *dsp, 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(AC3EncodeContext *s, int16_t *src, int len) 
274 
{ 
275 
int v = s>ac3dsp.ac3_max_msb_abs_int16(src, len);

276 
return av_log2(v);

277 
} 
278  
279  
280 
/**

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

282 
* @param tab input array

283 
* @param n number of values in the array

284 
* @param lshift left shift amount

285 
*/

286 
static void lshift_tab(int16_t *tab, int n, unsigned int lshift) 
287 
{ 
288 
int i;

289  
290 
if (lshift > 0) { 
291 
for (i = 0; i < n; i++) 
292 
tab[i] <<= lshift; 
293 
} 
294 
} 
295  
296  
297 
/**

298 
* Shift each value in an array by a specified amount.

299 
* @param src input array

300 
* @param n number of values in the array

301 
* @param shift shift amount (negative=right, positive=left)

302 
*/

303 
static void shift_int32(int32_t *src, int n, int shift) 
304 
{ 
305 
int i;

306  
307 
if (shift > 0) { 
308 
for (i = 0; i < n; i++) 
309 
src[i] <<= shift; 
310 
} else if (shift < 0) { 
311 
shift = shift; 
312 
for (i = 0; i < n; i++) 
313 
src[i] >>= shift; 
314 
} 
315 
} 
316  
317  
318 
/**

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

320 
* This assumes signed 16bit input samples.

321 
*

322 
* @return coefficient shift

323 
*/

324 
static int normalize_samples(AC3EncodeContext *s) 
325 
{ 
326 
int v = 14  log2_tab(s, s>windowed_samples, AC3_WINDOW_SIZE); 
327 
lshift_tab(s>windowed_samples, AC3_WINDOW_SIZE, v); 
328 
return 9  v; 
329 
} 
330  
331  
332 
/**

333 
* Scale MDCT coefficients to 24bit fixedpoint.

334 
*/

335 
static void scale_coefficients(AC3EncodeContext *s) 
336 
{ 
337 
int blk, ch;

338  
339 
for (blk = 0; blk < AC3_MAX_BLOCKS; blk++) { 
340 
AC3Block *block = &s>blocks[blk]; 
341 
for (ch = 0; ch < s>channels; ch++) { 
342 
shift_int32(block>mdct_coef[ch], AC3_MAX_COEFS, 
343 
block>coeff_shift[ch]); 
344 
} 
345 
} 
346 
} 
347  
348  
349 
#ifdef TEST

350 
/*************************************************************************/

351 
/* TEST */

352  
353 
#include "libavutil/lfg.h" 
354  
355 
#define MDCT_NBITS 9 
356 
#define MDCT_SAMPLES (1 << MDCT_NBITS) 
357 
#define FN (MDCT_SAMPLES/4) 
358  
359  
360 
static void fft_test(AC3MDCTContext *mdct, AVLFG *lfg) 
361 
{ 
362 
IComplex in[FN], in1[FN]; 
363 
int k, n, i;

364 
float sum_re, sum_im, a;

365  
366 
for (i = 0; i < FN; i++) { 
367 
in[i].re = av_lfg_get(lfg) % 65535  32767; 
368 
in[i].im = av_lfg_get(lfg) % 65535  32767; 
369 
in1[i] = in[i]; 
370 
} 
371 
fft(mdct, in, 7);

372  
373 
/* do it by hand */

374 
for (k = 0; k < FN; k++) { 
375 
sum_re = 0;

376 
sum_im = 0;

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

379 
sum_re += in1[n].re * cos(a)  in1[n].im * sin(a); 
380 
sum_im += in1[n].re * sin(a) + in1[n].im * cos(a); 
381 
} 
382 
av_log(NULL, AV_LOG_DEBUG, "%3d: %6d,%6d %6.0f,%6.0f\n", 
383 
k, in[k].re, in[k].im, sum_re / FN, sum_im / FN); 
384 
} 
385 
} 
386  
387  
388 
static void mdct_test(AC3MDCTContext *mdct, AVLFG *lfg) 
389 
{ 
390 
int16_t input[MDCT_SAMPLES]; 
391 
int32_t output[AC3_MAX_COEFS]; 
392 
float input1[MDCT_SAMPLES];

393 
float output1[AC3_MAX_COEFS];

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

395 
int i, k, n;

396  
397 
for (i = 0; i < MDCT_SAMPLES; i++) { 
398 
input[i] = (av_lfg_get(lfg) % 65535  32767) * 9 / 10; 
399 
input1[i] = input[i]; 
400 
} 
401  
402 
mdct512(mdct, output, input); 
403  
404 
/* do it by hand */

405 
for (k = 0; k < AC3_MAX_COEFS; k++) { 
406 
s = 0;

407 
for (n = 0; n < MDCT_SAMPLES; n++) { 
408 
a = (2*M_PI*(2*n+1+MDCT_SAMPLES/2)*(2*k+1) / (4 * MDCT_SAMPLES)); 
409 
s += input1[n] * cos(a); 
410 
} 
411 
output1[k] = 2 * s / MDCT_SAMPLES;

412 
} 
413  
414 
err = 0;

415 
emax = 0;

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

420 
emax = e; 
421 
err += e * e; 
422 
} 
423 
av_log(NULL, AV_LOG_DEBUG, "err2=%f emax=%f\n", err / AC3_MAX_COEFS, emax); 
424 
} 
425  
426  
427 
int main(void) 
428 
{ 
429 
AVLFG lfg; 
430 
AC3MDCTContext mdct; 
431  
432 
mdct.avctx = NULL;

433 
av_log_set_level(AV_LOG_DEBUG); 
434 
mdct_init(&mdct, 9);

435  
436 
fft_test(&mdct, &lfg); 
437 
mdct_test(&mdct, &lfg); 
438  
439 
return 0; 
440 
} 
441 
#endif /* TEST */ 
442  
443  
444 
AVCodec ff_ac3_fixed_encoder = { 
445 
"ac3_fixed",

446 
AVMEDIA_TYPE_AUDIO, 
447 
CODEC_ID_AC3, 
448 
sizeof(AC3EncodeContext),

449 
ac3_encode_init, 
450 
ac3_encode_frame, 
451 
ac3_encode_close, 
452 
NULL,

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

455 
.channel_layouts = ac3_channel_layouts, 
456 
}; 