dcaenc: update
[ffmpeg.git] / libavcodec / dcaenc.c
1 /*
2  * DCA encoder
3  * Copyright (C) 2008-2012 Alexander E. Patrakov
4  *               2010 Benjamin Larsson
5  *               2011 Xiang Wang
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 02110-1301 USA
22  */
23
24 #include "libavutil/avassert.h"
25 #include "libavutil/channel_layout.h"
26 #include "libavutil/common.h"
27 #include "avcodec.h"
28 #include "dca.h"
29 #include "dcadata.h"
30 #include "dcaenc.h"
31 #include "internal.h"
32 #include "put_bits.h"
33
34 #define MAX_CHANNELS 6
35 #define DCA_MAX_FRAME_SIZE 16384
36 #define DCA_HEADER_SIZE 13
37 #define DCA_LFE_SAMPLES 8
38
39 #define DCA_SUBBANDS 32
40 #define SUBFRAMES 1
41 #define SUBSUBFRAMES 2
42 #define SUBBAND_SAMPLES (SUBFRAMES * SUBSUBFRAMES * 8)
43 #define AUBANDS 25
44
45 typedef struct DCAContext {
46     PutBitContext pb;
47     int frame_size;
48     int frame_bits;
49     int fullband_channels;
50     int channels;
51     int lfe_channel;
52     int samplerate_index;
53     int bitrate_index;
54     int channel_config;
55     const int32_t *band_interpolation;
56     const int32_t *band_spectrum;
57     int lfe_scale_factor;
58     softfloat lfe_quant;
59     int32_t lfe_peak_cb;
60
61     int32_t history[512][MAX_CHANNELS]; /* This is a circular buffer */
62     int32_t subband[SUBBAND_SAMPLES][DCA_SUBBANDS][MAX_CHANNELS];
63     int32_t quantized[SUBBAND_SAMPLES][DCA_SUBBANDS][MAX_CHANNELS];
64     int32_t peak_cb[DCA_SUBBANDS][MAX_CHANNELS];
65     int32_t downsampled_lfe[DCA_LFE_SAMPLES];
66     int32_t masking_curve_cb[SUBSUBFRAMES][256];
67     int abits[DCA_SUBBANDS][MAX_CHANNELS];
68     int scale_factor[DCA_SUBBANDS][MAX_CHANNELS];
69     softfloat quant[DCA_SUBBANDS][MAX_CHANNELS];
70     int32_t eff_masking_curve_cb[256];
71     int32_t band_masking_cb[32];
72     int32_t worst_quantization_noise;
73     int32_t worst_noise_ever;
74     int consumed_bits;
75 } DCAContext;
76
77 static int32_t cos_table[2048];
78 static int32_t band_interpolation[2][512];
79 static int32_t band_spectrum[2][8];
80 static int32_t auf[9][AUBANDS][256];
81 static int32_t cb_to_add[256];
82 static int32_t cb_to_level[2048];
83 static int32_t lfe_fir_64i[512];
84
85 /* Transfer function of outer and middle ear, Hz -> dB */
86 static double hom(double f)
87 {
88     double f1 = f / 1000;
89
90     return -3.64 * pow(f1, -0.8)
91            + 6.8 * exp(-0.6 * (f1 - 3.4) * (f1 - 3.4))
92            - 6.0 * exp(-0.15 * (f1 - 8.7) * (f1 - 8.7))
93            - 0.0006 * (f1 * f1) * (f1 * f1);
94 }
95
96 static double gammafilter(int i, double f)
97 {
98     double h = (f - fc[i]) / erb[i];
99
100     h = 1 + h * h;
101     h = 1 / (h * h);
102     return 20 * log10(h);
103 }
104
105 static int encode_init(AVCodecContext *avctx)
106 {
107     DCAContext *c = avctx->priv_data;
108     uint64_t layout = avctx->channel_layout;
109     int i, min_frame_bits;
110
111     c->fullband_channels = c->channels = avctx->channels;
112     c->lfe_channel = (avctx->channels == 3 || avctx->channels == 6);
113     c->band_interpolation = band_interpolation[1];
114     c->band_spectrum = band_spectrum[1];
115     c->worst_quantization_noise = -2047;
116     c->worst_noise_ever = -2047;
117
118     if (!layout) {
119         av_log(avctx, AV_LOG_WARNING, "No channel layout specified. The "
120                                       "encoder will guess the layout, but it "
121                                       "might be incorrect.\n");
122         layout = av_get_default_channel_layout(avctx->channels);
123     }
124     switch (layout) {
125     case AV_CH_LAYOUT_MONO:         c->channel_config = 0; break;
126     case AV_CH_LAYOUT_STEREO:       c->channel_config = 2; break;
127     case AV_CH_LAYOUT_2_2:          c->channel_config = 8; break;
128     case AV_CH_LAYOUT_5POINT0:      c->channel_config = 9; break;
129     case AV_CH_LAYOUT_5POINT1:      c->channel_config = 9; break;
130     default:
131         av_log(avctx, AV_LOG_ERROR, "Unsupported channel layout!\n");
132         return AVERROR_PATCHWELCOME;
133     }
134
135     if (c->lfe_channel)
136         c->fullband_channels--;
137
138     for (i = 0; i < 9; i++) {
139         if (sample_rates[i] == avctx->sample_rate)
140             break;
141     }
142     if (i == 9)
143         return AVERROR(EINVAL);
144     c->samplerate_index = i;
145
146     if (avctx->bit_rate < 32000 || avctx->bit_rate > 3840000) {
147         av_log(avctx, AV_LOG_ERROR, "Bit rate %i not supported.", avctx->bit_rate);
148         return AVERROR(EINVAL);
149     }
150     for (i = 0; dca_bit_rates[i] < avctx->bit_rate; i++)
151         ;
152     c->bitrate_index = i;
153     avctx->bit_rate = dca_bit_rates[i];
154     c->frame_bits = FFALIGN((avctx->bit_rate * 512 + avctx->sample_rate - 1) / avctx->sample_rate, 32);
155     min_frame_bits = 132 + (493 + 28 * 32) * c->fullband_channels + c->lfe_channel * 72;
156     if (c->frame_bits < min_frame_bits || c->frame_bits > (DCA_MAX_FRAME_SIZE << 3))
157         return AVERROR(EINVAL);
158
159     c->frame_size = (c->frame_bits + 7) / 8;
160
161     avctx->frame_size = 32 * SUBBAND_SAMPLES;
162
163     if (!cos_table[0]) {
164         int j, k;
165
166         for (i = 0; i < 2048; i++) {
167             cos_table[i]   = (int32_t)(0x7fffffff * cos(M_PI * i / 1024));
168             cb_to_level[i] = (int32_t)(0x7fffffff * pow(10, -0.005 * i));
169         }
170
171         for (i = 0; i < 512; i++) {
172             lfe_fir_64i[i]           = (int32_t)(0x01ffffff * lfe_fir_64[i]);
173             band_interpolation[0][i] = (int32_t)(0x1000000000ULL * fir_32bands_perfect[i]);
174             band_interpolation[1][i] = (int32_t)(0x1000000000ULL * fir_32bands_nonperfect[i]);
175         }
176
177         for (i = 0; i < 9; i++) {
178             for (j = 0; j < AUBANDS; j++) {
179                 for (k = 0; k < 256; k++) {
180                     double freq = sample_rates[i] * (k + 0.5) / 512;
181
182                     auf[i][j][k] = (int32_t)(10 * (hom(freq) + gammafilter(j, freq)));
183                 }
184             }
185         }
186
187         for (i = 0; i < 256; i++) {
188             double add = 1 + pow(10, -0.01 * i);
189             cb_to_add[i] = (int32_t)(100 * log10(add));
190         }
191         for (j = 0; j < 8; j++) {
192             double accum = 0;
193             for (i = 0; i < 512; i++) {
194                 double reconst = fir_32bands_perfect[i] * ((i & 64) ? (-1) : 1);
195                 accum += reconst * cos(2 * M_PI * (i + 0.5 - 256) * (j + 0.5) / 512);
196             }
197             band_spectrum[0][j] = (int32_t)(200 * log10(accum));
198         }
199         for (j = 0; j < 8; j++) {
200             double accum = 0;
201             for (i = 0; i < 512; i++) {
202                 double reconst = fir_32bands_nonperfect[i] * ((i & 64) ? (-1) : 1);
203                 accum += reconst * cos(2 * M_PI * (i + 0.5 - 256) * (j + 0.5) / 512);
204             }
205             band_spectrum[1][j] = (int32_t)(200 * log10(accum));
206         }
207     }
208     return 0;
209 }
210
211 static inline int32_t cos_t(int x)
212 {
213     return cos_table[x & 2047];
214 }
215
216 static inline int32_t sin_t(int x)
217 {
218     return cos_t(x - 512);
219 }
220
221 static inline int32_t half32(int32_t a)
222 {
223     return (a + 1) >> 1;
224 }
225
226 static inline int32_t mul32(int32_t a, int32_t b)
227 {
228     int64_t r = (int64_t)a * b + 0x80000000ULL;
229     return r >> 32;
230 }
231
232 static void subband_transform(DCAContext *c, const int32_t *input)
233 {
234     int ch, subs, i, k, j;
235
236     for (ch = 0; ch < c->fullband_channels; ch++) {
237         /* History is copied because it is also needed for PSY */
238         int32_t hist[512];
239         int hist_start = 0;
240
241         for (i = 0; i < 512; i++)
242             hist[i] = c->history[i][ch];
243
244         for (subs = 0; subs < SUBBAND_SAMPLES; subs++) {
245             int32_t accum[64];
246             int32_t resp;
247             int band;
248
249             /* Calculate the convolutions at once */
250             for (i = 0; i < 64; i++)
251                 accum[i] = 0;
252
253             for (k = 0, i = hist_start, j = 0;
254                     i < 512; k = (k + 1) & 63, i++, j++)
255                 accum[k] += mul32(hist[i], c->band_interpolation[j]);
256             for (i = 0; i < hist_start; k = (k + 1) & 63, i++, j++)
257                 accum[k] += mul32(hist[i], c->band_interpolation[j]);
258
259             for (k = 16; k < 32; k++)
260                 accum[k] = accum[k] - accum[31 - k];
261             for (k = 32; k < 48; k++)
262                 accum[k] = accum[k] + accum[95 - k];
263
264             for (band = 0; band < 32; band++) {
265                 resp = 0;
266                 for (i = 16; i < 48; i++) {
267                     int s = (2 * band + 1) * (2 * (i + 16) + 1);
268                     resp += mul32(accum[i], cos_t(s << 3)) >> 3;
269                 }
270
271                 c->subband[subs][band][ch] = ((band + 1) & 2) ? -resp : resp;
272             }
273
274             /* Copy in 32 new samples from input */
275             for (i = 0; i < 32; i++)
276                 hist[i + hist_start] = input[(subs * 32 + i) * c->channels + ch];
277             hist_start = (hist_start + 32) & 511;
278         }
279     }
280 }
281
282 static void lfe_downsample(DCAContext *c, const int32_t *input)
283 {
284     /* FIXME: make 128x LFE downsampling possible */
285     int i, j, lfes;
286     int32_t hist[512];
287     int32_t accum;
288     int hist_start = 0;
289
290     for (i = 0; i < 512; i++)
291         hist[i] = c->history[i][c->channels - 1];
292
293     for (lfes = 0; lfes < DCA_LFE_SAMPLES; lfes++) {
294         /* Calculate the convolution */
295         accum = 0;
296
297         for (i = hist_start, j = 0; i < 512; i++, j++)
298             accum += mul32(hist[i], lfe_fir_64i[j]);
299         for (i = 0; i < hist_start; i++, j++)
300             accum += mul32(hist[i], lfe_fir_64i[j]);
301
302         c->downsampled_lfe[lfes] = accum;
303
304         /* Copy in 64 new samples from input */
305         for (i = 0; i < 64; i++)
306             hist[i + hist_start] = input[(lfes * 64 + i) * c->channels + c->channels - 1];
307
308         hist_start = (hist_start + 64) & 511;
309     }
310 }
311
312 typedef struct {
313     int32_t re;
314     int32_t im;
315 } cplx32;
316
317 static void fft(const int32_t in[2 * 256], cplx32 out[256])
318 {
319     cplx32 buf[256], rin[256], rout[256];
320     int i, j, k, l;
321
322     /* do two transforms in parallel */
323     for (i = 0; i < 256; i++) {
324         /* Apply the Hann window */
325         rin[i].re = mul32(in[2 * i], 0x3fffffff - (cos_t(8 * i + 2) >> 1));
326         rin[i].im = mul32(in[2 * i + 1], 0x3fffffff - (cos_t(8 * i + 6) >> 1));
327     }
328     /* pre-rotation */
329     for (i = 0; i < 256; i++) {
330         buf[i].re = mul32(cos_t(4 * i + 2), rin[i].re)
331                   - mul32(sin_t(4 * i + 2), rin[i].im);
332         buf[i].im = mul32(cos_t(4 * i + 2), rin[i].im)
333                   + mul32(sin_t(4 * i + 2), rin[i].re);
334     }
335
336     for (j = 256, l = 1; j != 1; j >>= 1, l <<= 1) {
337         for (k = 0; k < 256; k += j) {
338             for (i = k; i < k + j / 2; i++) {
339                 cplx32 sum, diff;
340                 int t = 8 * l * i;
341
342                 sum.re = buf[i].re + buf[i + j / 2].re;
343                 sum.im = buf[i].im + buf[i + j / 2].im;
344
345                 diff.re = buf[i].re - buf[i + j / 2].re;
346                 diff.im = buf[i].im - buf[i + j / 2].im;
347
348                 buf[i].re = half32(sum.re);
349                 buf[i].im = half32(sum.im);
350
351                 buf[i + j / 2].re = mul32(diff.re, cos_t(t))
352                                   - mul32(diff.im, sin_t(t));
353                 buf[i + j / 2].im = mul32(diff.im, cos_t(t))
354                                   + mul32(diff.re, sin_t(t));
355             }
356         }
357     }
358     /* post-rotation */
359     for (i = 0; i < 256; i++) {
360         int b = ff_reverse[i];
361         rout[i].re = mul32(buf[b].re, cos_t(4 * i))
362                    - mul32(buf[b].im, sin_t(4 * i));
363         rout[i].im = mul32(buf[b].im, cos_t(4 * i))
364                    + mul32(buf[b].re, sin_t(4 * i));
365     }
366     for (i = 0; i < 256; i++) {
367         /* separate the results of the two transforms */
368         cplx32 o1, o2;
369
370         o1.re =  rout[i].re - rout[255 - i].re;
371         o1.im =  rout[i].im + rout[255 - i].im;
372
373         o2.re =  rout[i].im - rout[255 - i].im;
374         o2.im = -rout[i].re - rout[255 - i].re;
375
376         /* combine them into one long transform */
377         out[i].re = mul32( o1.re + o2.re, cos_t(2 * i + 1))
378                   + mul32( o1.im - o2.im, sin_t(2 * i + 1));
379         out[i].im = mul32( o1.im + o2.im, cos_t(2 * i + 1))
380                   + mul32(-o1.re + o2.re, sin_t(2 * i + 1));
381     }
382 }
383
384 static int32_t get_cb(int32_t in)
385 {
386     int i, res;
387
388     res = 0;
389     if (in < 0)
390         in = -in;
391     for (i = 1024; i > 0; i >>= 1) {
392         if (cb_to_level[i + res] >= in)
393             res += i;
394     }
395     return -res;
396 }
397
398 static int32_t add_cb(int32_t a, int32_t b)
399 {
400     if (a < b)
401         FFSWAP(int32_t, a, b);
402
403     if (a - b >= 256)
404         return a;
405     return a + cb_to_add[a - b];
406 }
407
408 static void adjust_jnd(int samplerate_index,
409                        const int32_t in[512], int32_t out_cb[256])
410 {
411     int32_t power[256];
412     cplx32 out[256];
413     int32_t out_cb_unnorm[256];
414     int32_t denom;
415     const int32_t ca_cb = -1114;
416     const int32_t cs_cb = 928;
417     int i, j;
418
419     fft(in, out);
420
421     for (j = 0; j < 256; j++) {
422         power[j] = add_cb(get_cb(out[j].re), get_cb(out[j].im));
423         out_cb_unnorm[j] = -2047; /* and can only grow */
424     }
425
426     for (i = 0; i < AUBANDS; i++) {
427         denom = ca_cb; /* and can only grow */
428         for (j = 0; j < 256; j++)
429             denom = add_cb(denom, power[j] + auf[samplerate_index][i][j]);
430         for (j = 0; j < 256; j++)
431             out_cb_unnorm[j] = add_cb(out_cb_unnorm[j],
432                     -denom + auf[samplerate_index][i][j]);
433     }
434
435     for (j = 0; j < 256; j++)
436         out_cb[j] = add_cb(out_cb[j], -out_cb_unnorm[j] - ca_cb - cs_cb);
437 }
438
439 typedef void (*walk_band_t)(DCAContext *c, int band1, int band2, int f,
440                             int32_t spectrum1, int32_t spectrum2, int channel,
441                             int32_t * arg);
442
443 static void walk_band_low(DCAContext *c, int band, int channel,
444                           walk_band_t walk, int32_t *arg)
445 {
446     int f;
447
448     if (band == 0) {
449         for (f = 0; f < 4; f++)
450             walk(c, 0, 0, f, 0, -2047, channel, arg);
451     } else {
452         for (f = 0; f < 8; f++)
453             walk(c, band, band - 1, 8 * band - 4 + f,
454                     c->band_spectrum[7 - f], c->band_spectrum[f], channel, arg);
455     }
456 }
457
458 static void walk_band_high(DCAContext *c, int band, int channel,
459                            walk_band_t walk, int32_t *arg)
460 {
461     int f;
462
463     if (band == 31) {
464         for (f = 0; f < 4; f++)
465             walk(c, 31, 31, 256 - 4 + f, 0, -2047, channel, arg);
466     } else {
467         for (f = 0; f < 8; f++)
468             walk(c, band, band + 1, 8 * band + 4 + f,
469                     c->band_spectrum[f], c->band_spectrum[7 - f], channel, arg);
470     }
471 }
472
473 static void update_band_masking(DCAContext *c, int band1, int band2,
474                                 int f, int32_t spectrum1, int32_t spectrum2,
475                                 int channel, int32_t * arg)
476 {
477     int32_t value = c->eff_masking_curve_cb[f] - spectrum1;
478
479     if (value < c->band_masking_cb[band1])
480         c->band_masking_cb[band1] = value;
481 }
482
483 static void calc_masking(DCAContext *c, const int32_t *input)
484 {
485     int i, k, band, ch, ssf;
486     int32_t data[512];
487
488     for (i = 0; i < 256; i++)
489         for (ssf = 0; ssf < SUBSUBFRAMES; ssf++)
490             c->masking_curve_cb[ssf][i] = -2047;
491
492     for (ssf = 0; ssf < SUBSUBFRAMES; ssf++)
493         for (ch = 0; ch < c->fullband_channels; ch++) {
494             for (i = 0, k = 128 + 256 * ssf; k < 512; i++, k++)
495                 data[i] = c->history[k][ch];
496             for (k -= 512; i < 512; i++, k++)
497                 data[i] = input[k * c->channels + ch];
498             adjust_jnd(c->samplerate_index, data, c->masking_curve_cb[ssf]);
499         }
500     for (i = 0; i < 256; i++) {
501         int32_t m = 2048;
502
503         for (ssf = 0; ssf < SUBSUBFRAMES; ssf++)
504             if (c->masking_curve_cb[ssf][i] < m)
505                 m = c->masking_curve_cb[ssf][i];
506         c->eff_masking_curve_cb[i] = m;
507     }
508
509     for (band = 0; band < 32; band++) {
510         c->band_masking_cb[band] = 2048;
511         walk_band_low(c, band, 0, update_band_masking, NULL);
512         walk_band_high(c, band, 0, update_band_masking, NULL);
513     }
514 }
515
516 static void find_peaks(DCAContext *c)
517 {
518     int band, ch;
519
520     for (band = 0; band < 32; band++)
521         for (ch = 0; ch < c->fullband_channels; ch++) {
522             int sample;
523             int32_t m = 0;
524
525             for (sample = 0; sample < SUBBAND_SAMPLES; sample++) {
526                 int32_t s = abs(c->subband[sample][band][ch]);
527                 if (m < s)
528                     m = s;
529             }
530             c->peak_cb[band][ch] = get_cb(m);
531         }
532
533     if (c->lfe_channel) {
534         int sample;
535         int32_t m = 0;
536
537         for (sample = 0; sample < DCA_LFE_SAMPLES; sample++)
538             if (m < abs(c->downsampled_lfe[sample]))
539                 m = abs(c->downsampled_lfe[sample]);
540         c->lfe_peak_cb = get_cb(m);
541     }
542 }
543
544 static const int snr_fudge = 128;
545 #define USED_1ABITS 1
546 #define USED_NABITS 2
547 #define USED_26ABITS 4
548
549 static int init_quantization_noise(DCAContext *c, int noise)
550 {
551     int ch, band, ret = 0;
552
553     c->consumed_bits = 132 + 493 * c->fullband_channels;
554     if (c->lfe_channel)
555         c->consumed_bits += 72;
556
557     /* attempt to guess the bit distribution based on the prevoius frame */
558     for (ch = 0; ch < c->fullband_channels; ch++) {
559         for (band = 0; band < 32; band++) {
560             int snr_cb = c->peak_cb[band][ch] - c->band_masking_cb[band] - noise;
561
562             if (snr_cb >= 1312) {
563                 c->abits[band][ch] = 26;
564                 ret |= USED_26ABITS;
565             } else if (snr_cb >= 222) {
566                 c->abits[band][ch] = 8 + mul32(snr_cb - 222, 69000000);
567                 ret |= USED_NABITS;
568             } else if (snr_cb >= 0) {
569                 c->abits[band][ch] = 2 + mul32(snr_cb, 106000000);
570                 ret |= USED_NABITS;
571             } else {
572                 c->abits[band][ch] = 1;
573                 ret |= USED_1ABITS;
574             }
575         }
576     }
577
578     for (band = 0; band < 32; band++)
579         for (ch = 0; ch < c->fullband_channels; ch++) {
580             c->consumed_bits += bit_consumption[c->abits[band][ch]];
581         }
582
583     return ret;
584 }
585
586 static void assign_bits(DCAContext *c)
587 {
588     /* Find the bounds where the binary search should work */
589     int low, high, down;
590     int used_abits = 0;
591
592     init_quantization_noise(c, c->worst_quantization_noise);
593     low = high = c->worst_quantization_noise;
594     if (c->consumed_bits > c->frame_bits) {
595         while (c->consumed_bits > c->frame_bits) {
596             av_assert0(used_abits != USED_1ABITS);
597             low = high;
598             high += snr_fudge;
599             used_abits = init_quantization_noise(c, high);
600         }
601     } else {
602         while (c->consumed_bits <= c->frame_bits) {
603             high = low;
604             if (used_abits == USED_26ABITS)
605                 goto out; /* The requested bitrate is too high, pad with zeros */
606             low -= snr_fudge;
607             used_abits = init_quantization_noise(c, low);
608         }
609     }
610
611     /* Now do a binary search between low and high to see what fits */
612     for (down = snr_fudge >> 1; down; down >>= 1) {
613         init_quantization_noise(c, high - down);
614         if (c->consumed_bits <= c->frame_bits)
615             high -= down;
616     }
617     init_quantization_noise(c, high);
618 out:
619     c->worst_quantization_noise = high;
620     if (high > c->worst_noise_ever)
621         c->worst_noise_ever = high;
622 }
623
624 static void shift_history(DCAContext *c, const int32_t *input)
625 {
626     int k, ch;
627
628     for (k = 0; k < 512; k++)
629         for (ch = 0; ch < c->channels; ch++)
630             c->history[k][ch] = input[k * c->channels + ch];
631 }
632
633 static int32_t quantize_value(int32_t value, softfloat quant)
634 {
635     int32_t offset = 1 << (quant.e - 1);
636
637     value = mul32(value, quant.m) + offset;
638     value = value >> quant.e;
639     return value;
640 }
641
642 static int calc_one_scale(int32_t peak_cb, int abits, softfloat *quant)
643 {
644     int32_t peak;
645     int our_nscale, try_remove;
646     softfloat our_quant;
647
648     av_assert0(peak_cb <= 0);
649     av_assert0(peak_cb >= -2047);
650
651     our_nscale = 127;
652     peak = cb_to_level[-peak_cb];
653
654     for (try_remove = 64; try_remove > 0; try_remove >>= 1) {
655         if (scalefactor_inv[our_nscale - try_remove].e + stepsize_inv[abits].e <= 17)
656             continue;
657         our_quant.m = mul32(scalefactor_inv[our_nscale - try_remove].m, stepsize_inv[abits].m);
658         our_quant.e = scalefactor_inv[our_nscale - try_remove].e + stepsize_inv[abits].e - 17;
659         if ((quant_levels[abits] - 1) / 2 < quantize_value(peak, our_quant))
660             continue;
661         our_nscale -= try_remove;
662     }
663
664     if (our_nscale >= 125)
665         our_nscale = 124;
666
667     quant->m = mul32(scalefactor_inv[our_nscale].m, stepsize_inv[abits].m);
668     quant->e = scalefactor_inv[our_nscale].e + stepsize_inv[abits].e - 17;
669     av_assert0((quant_levels[abits] - 1) / 2 >= quantize_value(peak, *quant));
670
671     return our_nscale;
672 }
673
674 static void calc_scales(DCAContext *c)
675 {
676     int band, ch;
677
678     for (band = 0; band < 32; band++)
679         for (ch = 0; ch < c->fullband_channels; ch++)
680             c->scale_factor[band][ch] = calc_one_scale(c->peak_cb[band][ch],
681                                                        c->abits[band][ch],
682                                                        &c->quant[band][ch]);
683
684     if (c->lfe_channel)
685         c->lfe_scale_factor = calc_one_scale(c->lfe_peak_cb, 11, &c->lfe_quant);
686 }
687
688 static void quantize_all(DCAContext *c)
689 {
690     int sample, band, ch;
691
692     for (sample = 0; sample < SUBBAND_SAMPLES; sample++)
693         for (band = 0; band < 32; band++)
694             for (ch = 0; ch < c->fullband_channels; ch++)
695                 c->quantized[sample][band][ch] = quantize_value(c->subband[sample][band][ch], c->quant[band][ch]);
696 }
697
698 static void put_frame_header(DCAContext *c)
699 {
700     /* SYNC */
701     put_bits(&c->pb, 16, 0x7ffe);
702     put_bits(&c->pb, 16, 0x8001);
703
704     /* Frame type: normal */
705     put_bits(&c->pb, 1, 1);
706
707     /* Deficit sample count: none */
708     put_bits(&c->pb, 5, 31);
709
710     /* CRC is not present */
711     put_bits(&c->pb, 1, 0);
712
713     /* Number of PCM sample blocks */
714     put_bits(&c->pb, 7, SUBBAND_SAMPLES - 1);
715
716     /* Primary frame byte size */
717     put_bits(&c->pb, 14, c->frame_size - 1);
718
719     /* Audio channel arrangement */
720     put_bits(&c->pb, 6, c->channel_config);
721
722     /* Core audio sampling frequency */
723     put_bits(&c->pb, 4, bitstream_sfreq[c->samplerate_index]);
724
725     /* Transmission bit rate */
726     put_bits(&c->pb, 5, c->bitrate_index);
727
728     /* Embedded down mix: disabled */
729     put_bits(&c->pb, 1, 0);
730
731     /* Embedded dynamic range flag: not present */
732     put_bits(&c->pb, 1, 0);
733
734     /* Embedded time stamp flag: not present */
735     put_bits(&c->pb, 1, 0);
736
737     /* Auxiliary data flag: not present */
738     put_bits(&c->pb, 1, 0);
739
740     /* HDCD source: no */
741     put_bits(&c->pb, 1, 0);
742
743     /* Extension audio ID: N/A */
744     put_bits(&c->pb, 3, 0);
745
746     /* Extended audio data: not present */
747     put_bits(&c->pb, 1, 0);
748
749     /* Audio sync word insertion flag: after each sub-frame */
750     put_bits(&c->pb, 1, 0);
751
752     /* Low frequency effects flag: not present or 64x subsampling */
753     put_bits(&c->pb, 2, c->lfe_channel ? 2 : 0);
754
755     /* Predictor history switch flag: on */
756     put_bits(&c->pb, 1, 1);
757
758     /* No CRC */
759     /* Multirate interpolator switch: non-perfect reconstruction */
760     put_bits(&c->pb, 1, 0);
761
762     /* Encoder software revision: 7 */
763     put_bits(&c->pb, 4, 7);
764
765     /* Copy history: 0 */
766     put_bits(&c->pb, 2, 0);
767
768     /* Source PCM resolution: 16 bits, not DTS ES */
769     put_bits(&c->pb, 3, 0);
770
771     /* Front sum/difference coding: no */
772     put_bits(&c->pb, 1, 0);
773
774     /* Surrounds sum/difference coding: no */
775     put_bits(&c->pb, 1, 0);
776
777     /* Dialog normalization: 0 dB */
778     put_bits(&c->pb, 4, 0);
779 }
780
781 static void put_primary_audio_header(DCAContext *c)
782 {
783     static const int bitlen[11] = { 0, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 };
784     static const int thr[11]    = { 0, 1, 3, 3, 3, 3, 7, 7, 7, 7, 7 };
785
786     int ch, i;
787     /* Number of subframes */
788     put_bits(&c->pb, 4, SUBFRAMES - 1);
789
790     /* Number of primary audio channels */
791     put_bits(&c->pb, 3, c->fullband_channels - 1);
792
793     /* Subband activity count */
794     for (ch = 0; ch < c->fullband_channels; ch++)
795         put_bits(&c->pb, 5, DCA_SUBBANDS - 2);
796
797     /* High frequency VQ start subband */
798     for (ch = 0; ch < c->fullband_channels; ch++)
799         put_bits(&c->pb, 5, DCA_SUBBANDS - 1);
800
801     /* Joint intensity coding index: 0, 0 */
802     for (ch = 0; ch < c->fullband_channels; ch++)
803         put_bits(&c->pb, 3, 0);
804
805     /* Transient mode codebook: A4, A4 (arbitrary) */
806     for (ch = 0; ch < c->fullband_channels; ch++)
807         put_bits(&c->pb, 2, 0);
808
809     /* Scale factor code book: 7 bit linear, 7-bit sqrt table (for each channel) */
810     for (ch = 0; ch < c->fullband_channels; ch++)
811         put_bits(&c->pb, 3, 6);
812
813     /* Bit allocation quantizer select: linear 5-bit */
814     for (ch = 0; ch < c->fullband_channels; ch++)
815         put_bits(&c->pb, 3, 6);
816
817     /* Quantization index codebook select: dummy data
818        to avoid transmission of scale factor adjustment */
819     for (i = 1; i < 11; i++)
820         for (ch = 0; ch < c->fullband_channels; ch++)
821             put_bits(&c->pb, bitlen[i], thr[i]);
822
823     /* Scale factor adjustment index: not transmitted */
824     /* Audio header CRC check word: not transmitted */
825 }
826
827 static void put_subframe_samples(DCAContext *c, int ss, int band, int ch)
828 {
829     if (c->abits[band][ch] <= 7) {
830         int sum, i, j;
831         for (i = 0; i < 8; i += 4) {
832             sum = 0;
833             for (j = 3; j >= 0; j--) {
834                 sum *= quant_levels[c->abits[band][ch]];
835                 sum += c->quantized[ss * 8 + i + j][band][ch];
836                 sum += (quant_levels[c->abits[band][ch]] - 1) / 2;
837             }
838             put_bits(&c->pb, bit_consumption[c->abits[band][ch]] / 4, sum);
839         }
840     } else {
841         int i;
842         for (i = 0; i < 8; i++) {
843             int bits = bit_consumption[c->abits[band][ch]] / 16;
844             int32_t mask = (1 << bits) - 1;
845             put_bits(&c->pb, bits, c->quantized[ss * 8 + i][band][ch] & mask);
846         }
847     }
848 }
849
850 static void put_subframe(DCAContext *c, int subframe)
851 {
852     int i, band, ss, ch;
853
854     /* Subsubframes count */
855     put_bits(&c->pb, 2, SUBSUBFRAMES -1);
856
857     /* Partial subsubframe sample count: dummy */
858     put_bits(&c->pb, 3, 0);
859
860     /* Prediction mode: no ADPCM, in each channel and subband */
861     for (ch = 0; ch < c->fullband_channels; ch++)
862         for (band = 0; band < DCA_SUBBANDS; band++)
863             put_bits(&c->pb, 1, 0);
864
865     /* Prediction VQ addres: not transmitted */
866     /* Bit allocation index */
867     for (ch = 0; ch < c->fullband_channels; ch++)
868         for (band = 0; band < DCA_SUBBANDS; band++)
869             put_bits(&c->pb, 5, c->abits[band][ch]);
870
871     if (SUBSUBFRAMES > 1) {
872         /* Transition mode: none for each channel and subband */
873         for (ch = 0; ch < c->fullband_channels; ch++)
874             for (band = 0; band < DCA_SUBBANDS; band++)
875                 put_bits(&c->pb, 1, 0); /* codebook A4 */
876     }
877
878     /* Scale factors */
879     for (ch = 0; ch < c->fullband_channels; ch++)
880         for (band = 0; band < DCA_SUBBANDS; band++)
881             put_bits(&c->pb, 7, c->scale_factor[band][ch]);
882
883     /* Joint subband scale factor codebook select: not transmitted */
884     /* Scale factors for joint subband coding: not transmitted */
885     /* Stereo down-mix coefficients: not transmitted */
886     /* Dynamic range coefficient: not transmitted */
887     /* Stde information CRC check word: not transmitted */
888     /* VQ encoded high frequency subbands: not transmitted */
889
890     /* LFE data: 8 samples and scalefactor */
891     if (c->lfe_channel) {
892         for (i = 0; i < DCA_LFE_SAMPLES; i++)
893             put_bits(&c->pb, 8, quantize_value(c->downsampled_lfe[i], c->lfe_quant) & 0xff);
894         put_bits(&c->pb, 8, c->lfe_scale_factor);
895     }
896
897     /* Audio data (subsubframes) */
898     for (ss = 0; ss < SUBSUBFRAMES ; ss++)
899         for (ch = 0; ch < c->fullband_channels; ch++)
900             for (band = 0; band < DCA_SUBBANDS; band++)
901                     put_subframe_samples(c, ss, band, ch);
902
903     /* DSYNC */
904     put_bits(&c->pb, 16, 0xffff);
905 }
906
907 static int encode_frame(AVCodecContext *avctx, AVPacket *avpkt,
908                         const AVFrame *frame, int *got_packet_ptr)
909 {
910     DCAContext *c = avctx->priv_data;
911     const int32_t *samples;
912     int ret, i;
913
914     if ((ret = ff_alloc_packet2(avctx, avpkt, c->frame_size )) < 0)
915         return ret;
916
917     samples = (const int32_t *)frame->data[0];
918
919     subband_transform(c, samples);
920     if (c->lfe_channel)
921         lfe_downsample(c, samples);
922
923     calc_masking(c, samples);
924     find_peaks(c);
925     assign_bits(c);
926     calc_scales(c);
927     quantize_all(c);
928     shift_history(c, samples);
929
930     init_put_bits(&c->pb, avpkt->data, avpkt->size);
931     put_frame_header(c);
932     put_primary_audio_header(c);
933     for (i = 0; i < SUBFRAMES; i++)
934         put_subframe(c, i);
935
936     flush_put_bits(&c->pb);
937
938     avpkt->pts      = frame->pts;
939     avpkt->duration = ff_samples_to_time_base(avctx, frame->nb_samples);
940     avpkt->size     = c->frame_size + 1;
941     *got_packet_ptr = 1;
942     return 0;
943 }
944
945 static const AVCodecDefault defaults[] = {
946     { "b",          "1411200" },
947     { NULL },
948 };
949
950 AVCodec ff_dca_encoder = {
951     .name                  = "dca",
952     .type                  = AVMEDIA_TYPE_AUDIO,
953     .id                    = AV_CODEC_ID_DTS,
954     .priv_data_size        = sizeof(DCAContext),
955     .init                  = encode_init,
956     .encode2               = encode_frame,
957     .capabilities          = CODEC_CAP_EXPERIMENTAL,
958     .sample_fmts           = (const enum AVSampleFormat[]){ AV_SAMPLE_FMT_S32,
959                                                             AV_SAMPLE_FMT_NONE },
960     .long_name             = NULL_IF_CONFIG_SMALL("DCA (DTS Coherent Acoustics)"),
961     .supported_samplerates = sample_rates,
962     .channel_layouts       = (const uint64_t[]) { AV_CH_LAYOUT_MONO,
963                                                   AV_CH_LAYOUT_STEREO,
964                                                   AV_CH_LAYOUT_2_2,
965                                                   AV_CH_LAYOUT_5POINT0,
966                                                   AV_CH_LAYOUT_5POINT1,
967                                                   0 },
968     .defaults              = defaults,
969 };