ffmpeg / libavcodec / ac3enc_fixed.c @ e6e98234
History  View  Annotate  Download (10.7 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 Libav.

8 
*

9 
* Libav 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 
* Libav 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 Libav; 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, rshift) \

135 
{ \ 
136 
pre = (MUL16(are, bre)  MUL16(aim, bim)) >> rshift; \ 
137 
pim = (MUL16(are, bim) + MUL16(bre, aim)) >> rshift; \ 
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, 15);

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], 15);

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], 0); 
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, unsigned int len) 
256 
{ 
257 
dsp>apply_window_int16(output, input, window, len); 
258 
} 
259  
260  
261 
/**

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

263 
* @param tab input array

264 
* @param n number of values in the array

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

266 
*/

267 
static int log2_tab(AC3EncodeContext *s, int16_t *src, int len) 
268 
{ 
269 
int v = s>ac3dsp.ac3_max_msb_abs_int16(src, len);

270 
return av_log2(v);

271 
} 
272  
273  
274 
/**

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

276 
* This assumes signed 16bit input samples.

277 
*

278 
* @return exponent shift

279 
*/

280 
static int normalize_samples(AC3EncodeContext *s) 
281 
{ 
282 
int v = 14  log2_tab(s, s>windowed_samples, AC3_WINDOW_SIZE); 
283 
if (v > 0) 
284 
s>ac3dsp.ac3_lshift_int16(s>windowed_samples, AC3_WINDOW_SIZE, v); 
285 
/* +6 to rightshift from 31bit to 25bit */

286 
return v + 6; 
287 
} 
288  
289  
290 
/**

291 
* Scale MDCT coefficients to 25bit signed fixedpoint.

292 
*/

293 
static void scale_coefficients(AC3EncodeContext *s) 
294 
{ 
295 
int blk, ch;

296  
297 
for (blk = 0; blk < AC3_MAX_BLOCKS; blk++) { 
298 
AC3Block *block = &s>blocks[blk]; 
299 
for (ch = 0; ch < s>channels; ch++) { 
300 
s>ac3dsp.ac3_rshift_int32(block>mdct_coef[ch], AC3_MAX_COEFS, 
301 
block>coeff_shift[ch]); 
302 
} 
303 
} 
304 
} 
305  
306  
307 
#ifdef TEST

308 
/*************************************************************************/

309 
/* TEST */

310  
311 
#include "libavutil/lfg.h" 
312  
313 
#define MDCT_NBITS 9 
314 
#define MDCT_SAMPLES (1 << MDCT_NBITS) 
315 
#define FN (MDCT_SAMPLES/4) 
316  
317  
318 
static void fft_test(AC3MDCTContext *mdct, AVLFG *lfg) 
319 
{ 
320 
IComplex in[FN], in1[FN]; 
321 
int k, n, i;

322 
float sum_re, sum_im, a;

323  
324 
for (i = 0; i < FN; i++) { 
325 
in[i].re = av_lfg_get(lfg) % 65535  32767; 
326 
in[i].im = av_lfg_get(lfg) % 65535  32767; 
327 
in1[i] = in[i]; 
328 
} 
329 
fft(mdct, in, 7);

330  
331 
/* do it by hand */

332 
for (k = 0; k < FN; k++) { 
333 
sum_re = 0;

334 
sum_im = 0;

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

337 
sum_re += in1[n].re * cos(a)  in1[n].im * sin(a); 
338 
sum_im += in1[n].re * sin(a) + in1[n].im * cos(a); 
339 
} 
340 
av_log(NULL, AV_LOG_DEBUG, "%3d: %6d,%6d %6.0f,%6.0f\n", 
341 
k, in[k].re, in[k].im, sum_re / FN, sum_im / FN); 
342 
} 
343 
} 
344  
345  
346 
static void mdct_test(AC3MDCTContext *mdct, AVLFG *lfg) 
347 
{ 
348 
int16_t input[MDCT_SAMPLES]; 
349 
int32_t output[AC3_MAX_COEFS]; 
350 
float input1[MDCT_SAMPLES];

351 
float output1[AC3_MAX_COEFS];

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

353 
int i, k, n;

354  
355 
for (i = 0; i < MDCT_SAMPLES; i++) { 
356 
input[i] = (av_lfg_get(lfg) % 65535  32767) * 9 / 10; 
357 
input1[i] = input[i]; 
358 
} 
359  
360 
mdct512(mdct, output, input); 
361  
362 
/* do it by hand */

363 
for (k = 0; k < AC3_MAX_COEFS; k++) { 
364 
s = 0;

365 
for (n = 0; n < MDCT_SAMPLES; n++) { 
366 
a = (2*M_PI*(2*n+1+MDCT_SAMPLES/2)*(2*k+1) / (4 * MDCT_SAMPLES)); 
367 
s += input1[n] * cos(a); 
368 
} 
369 
output1[k] = 2 * s / MDCT_SAMPLES;

370 
} 
371  
372 
err = 0;

373 
emax = 0;

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

378 
emax = e; 
379 
err += e * e; 
380 
} 
381 
av_log(NULL, AV_LOG_DEBUG, "err2=%f emax=%f\n", err / AC3_MAX_COEFS, emax); 
382 
} 
383  
384  
385 
int main(void) 
386 
{ 
387 
AVLFG lfg; 
388 
AC3MDCTContext mdct; 
389  
390 
mdct.avctx = NULL;

391 
av_log_set_level(AV_LOG_DEBUG); 
392 
mdct_init(&mdct, 9);

393  
394 
fft_test(&mdct, &lfg); 
395 
mdct_test(&mdct, &lfg); 
396  
397 
return 0; 
398 
} 
399 
#endif /* TEST */ 
400  
401  
402 
AVCodec ff_ac3_fixed_encoder = { 
403 
"ac3_fixed",

404 
AVMEDIA_TYPE_AUDIO, 
405 
CODEC_ID_AC3, 
406 
sizeof(AC3EncodeContext),

407 
ac3_encode_init, 
408 
ac3_encode_frame, 
409 
ac3_encode_close, 
410 
NULL,

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

413 
.channel_layouts = ac3_channel_layouts, 
414 
}; 