aaccoder: add a new perceptual noise substitution implementation
[ffmpeg.git] / libavcodec / aaccoder.c
1 /*
2  * AAC coefficients encoder
3  * Copyright (C) 2008-2009 Konstantin Shishkov
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 02110-1301 USA
20  */
21
22 /**
23  * @file
24  * AAC coefficients encoder
25  */
26
27 /***********************************
28  *              TODOs:
29  * speedup quantizer selection
30  * add sane pulse detection
31  ***********************************/
32
33 #include "libavutil/libm.h" // brought forward to work around cygwin header breakage
34
35 #include <float.h>
36 #include "libavutil/mathematics.h"
37 #include "avcodec.h"
38 #include "put_bits.h"
39 #include "aac.h"
40 #include "aacenc.h"
41 #include "aactab.h"
42
43 /** Frequency in Hz for lower limit of noise substitution **/
44 #define NOISE_LOW_LIMIT 4500
45
46 /* Energy spread threshold value below which no PNS is used, this corresponds to
47  * typically around 17Khz, after which PNS usage decays ending at 19Khz */
48 #define NOISE_SPREAD_THRESHOLD 0.5f
49
50 /* This constant gets divided by lambda to return ~1.65 which when multiplied
51  * by the band->threshold and compared to band->energy is the boundary between
52  * excessive PNS and little PNS usage. */
53 #define NOISE_LAMBDA_NUMERATOR 252.1f
54
55 /** Total number of usable codebooks **/
56 #define CB_TOT 12
57
58 /** Total number of codebooks, including special ones **/
59 #define CB_TOT_ALL 15
60
61 /** bits needed to code codebook run value for long windows */
62 static const uint8_t run_value_bits_long[64] = {
63      5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5,
64      5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5,  5, 10,
65     10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
66     10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 15
67 };
68
69 /** bits needed to code codebook run value for short windows */
70 static const uint8_t run_value_bits_short[16] = {
71     3, 3, 3, 3, 3, 3, 3, 6, 6, 6, 6, 6, 6, 6, 6, 9
72 };
73
74 static const uint8_t * const run_value_bits[2] = {
75     run_value_bits_long, run_value_bits_short
76 };
77
78 /** Map to convert values from BandCodingPath index to a codebook index **/
79 static const uint8_t aac_cb_out_map[CB_TOT_ALL]  = {0,1,2,3,4,5,6,7,8,9,10,11,13,14,15};
80 /** Inverse map to convert from codebooks to BandCodingPath indices **/
81 static const uint8_t aac_cb_in_map[CB_TOT_ALL+1] = {0,1,2,3,4,5,6,7,8,9,10,11,0,12,13,14};
82
83 /**
84  * Quantize one coefficient.
85  * @return absolute value of the quantized coefficient
86  * @see 3GPP TS26.403 5.6.2 "Scalefactor determination"
87  */
88 static av_always_inline int quant(float coef, const float Q)
89 {
90     float a = coef * Q;
91     return sqrtf(a * sqrtf(a)) + 0.4054;
92 }
93
94 static void quantize_bands(int *out, const float *in, const float *scaled,
95                            int size, float Q34, int is_signed, int maxval)
96 {
97     int i;
98     double qc;
99     for (i = 0; i < size; i++) {
100         qc = scaled[i] * Q34;
101         out[i] = (int)FFMIN(qc + 0.4054, (double)maxval);
102         if (is_signed && in[i] < 0.0f) {
103             out[i] = -out[i];
104         }
105     }
106 }
107
108 static void abs_pow34_v(float *out, const float *in, const int size)
109 {
110 #ifndef USE_REALLY_FULL_SEARCH
111     int i;
112     for (i = 0; i < size; i++) {
113         float a = fabsf(in[i]);
114         out[i] = sqrtf(a * sqrtf(a));
115     }
116 #endif /* USE_REALLY_FULL_SEARCH */
117 }
118
119 static const uint8_t aac_cb_range [12] = {0, 3, 3, 3, 3, 9, 9, 8, 8, 13, 13, 17};
120 static const uint8_t aac_cb_maxval[12] = {0, 1, 1, 2, 2, 4, 4, 7, 7, 12, 12, 16};
121
122 /**
123  * Calculate rate distortion cost for quantizing with given codebook
124  *
125  * @return quantization distortion
126  */
127 static av_always_inline float quantize_and_encode_band_cost_template(
128                                 struct AACEncContext *s,
129                                 PutBitContext *pb, const float *in,
130                                 const float *scaled, int size, int scale_idx,
131                                 int cb, const float lambda, const float uplim,
132                                 int *bits, int BT_ZERO, int BT_UNSIGNED,
133                                 int BT_PAIR, int BT_ESC, int BT_NOISE, int BT_STEREO)
134 {
135     const int q_idx = POW_SF2_ZERO - scale_idx + SCALE_ONE_POS - SCALE_DIV_512;
136     const float Q   = ff_aac_pow2sf_tab [q_idx];
137     const float Q34 = ff_aac_pow34sf_tab[q_idx];
138     const float IQ  = ff_aac_pow2sf_tab [POW_SF2_ZERO + scale_idx - SCALE_ONE_POS + SCALE_DIV_512];
139     const float CLIPPED_ESCAPE = 165140.0f*IQ;
140     int i, j;
141     float cost = 0;
142     const int dim = BT_PAIR ? 2 : 4;
143     int resbits = 0;
144     int off;
145
146     if (BT_ZERO || BT_NOISE || BT_STEREO) {
147         for (i = 0; i < size; i++)
148             cost += in[i]*in[i];
149         if (bits)
150             *bits = 0;
151         return cost * lambda;
152     }
153     if (!scaled) {
154         abs_pow34_v(s->scoefs, in, size);
155         scaled = s->scoefs;
156     }
157     quantize_bands(s->qcoefs, in, scaled, size, Q34, !BT_UNSIGNED, aac_cb_maxval[cb]);
158     if (BT_UNSIGNED) {
159         off = 0;
160     } else {
161         off = aac_cb_maxval[cb];
162     }
163     for (i = 0; i < size; i += dim) {
164         const float *vec;
165         int *quants = s->qcoefs + i;
166         int curidx = 0;
167         int curbits;
168         float rd = 0.0f;
169         for (j = 0; j < dim; j++) {
170             curidx *= aac_cb_range[cb];
171             curidx += quants[j] + off;
172         }
173         curbits =  ff_aac_spectral_bits[cb-1][curidx];
174         vec     = &ff_aac_codebook_vectors[cb-1][curidx*dim];
175         if (BT_UNSIGNED) {
176             for (j = 0; j < dim; j++) {
177                 float t = fabsf(in[i+j]);
178                 float di;
179                 if (BT_ESC && vec[j] == 64.0f) { //FIXME: slow
180                     if (t >= CLIPPED_ESCAPE) {
181                         di = t - CLIPPED_ESCAPE;
182                         curbits += 21;
183                     } else {
184                         int c = av_clip_uintp2(quant(t, Q), 13);
185                         di = t - c*cbrtf(c)*IQ;
186                         curbits += av_log2(c)*2 - 4 + 1;
187                     }
188                 } else {
189                     di = t - vec[j]*IQ;
190                 }
191                 if (vec[j] != 0.0f)
192                     curbits++;
193                 rd += di*di;
194             }
195         } else {
196             for (j = 0; j < dim; j++) {
197                 float di = in[i+j] - vec[j]*IQ;
198                 rd += di*di;
199             }
200         }
201         cost    += rd * lambda + curbits;
202         resbits += curbits;
203         if (cost >= uplim)
204             return uplim;
205         if (pb) {
206             put_bits(pb, ff_aac_spectral_bits[cb-1][curidx], ff_aac_spectral_codes[cb-1][curidx]);
207             if (BT_UNSIGNED)
208                 for (j = 0; j < dim; j++)
209                     if (ff_aac_codebook_vectors[cb-1][curidx*dim+j] != 0.0f)
210                         put_bits(pb, 1, in[i+j] < 0.0f);
211             if (BT_ESC) {
212                 for (j = 0; j < 2; j++) {
213                     if (ff_aac_codebook_vectors[cb-1][curidx*2+j] == 64.0f) {
214                         int coef = av_clip_uintp2(quant(fabsf(in[i+j]), Q), 13);
215                         int len = av_log2(coef);
216
217                         put_bits(pb, len - 4 + 1, (1 << (len - 4 + 1)) - 2);
218                         put_sbits(pb, len, coef);
219                     }
220                 }
221             }
222         }
223     }
224
225     if (bits)
226         *bits = resbits;
227     return cost;
228 }
229
230 static float quantize_and_encode_band_cost_NONE(struct AACEncContext *s, PutBitContext *pb,
231                                                 const float *in, const float *scaled,
232                                                 int size, int scale_idx, int cb,
233                                                 const float lambda, const float uplim,
234                                                 int *bits) {
235     av_assert0(0);
236     return 0.0f;
237 }
238
239 #define QUANTIZE_AND_ENCODE_BAND_COST_FUNC(NAME, BT_ZERO, BT_UNSIGNED, BT_PAIR, BT_ESC, BT_NOISE, BT_STEREO) \
240 static float quantize_and_encode_band_cost_ ## NAME(                                         \
241                                 struct AACEncContext *s,                                     \
242                                 PutBitContext *pb, const float *in,                          \
243                                 const float *scaled, int size, int scale_idx,                \
244                                 int cb, const float lambda, const float uplim,               \
245                                 int *bits) {                                                 \
246     return quantize_and_encode_band_cost_template(                                           \
247                                 s, pb, in, scaled, size, scale_idx,                          \
248                                 BT_ESC ? ESC_BT : cb, lambda, uplim, bits,                   \
249                                 BT_ZERO, BT_UNSIGNED, BT_PAIR, BT_ESC, BT_NOISE, BT_STEREO); \
250 }
251
252 QUANTIZE_AND_ENCODE_BAND_COST_FUNC(ZERO,  1, 0, 0, 0, 0, 0)
253 QUANTIZE_AND_ENCODE_BAND_COST_FUNC(SQUAD, 0, 0, 0, 0, 0, 0)
254 QUANTIZE_AND_ENCODE_BAND_COST_FUNC(UQUAD, 0, 1, 0, 0, 0, 0)
255 QUANTIZE_AND_ENCODE_BAND_COST_FUNC(SPAIR, 0, 0, 1, 0, 0, 0)
256 QUANTIZE_AND_ENCODE_BAND_COST_FUNC(UPAIR, 0, 1, 1, 0, 0, 0)
257 QUANTIZE_AND_ENCODE_BAND_COST_FUNC(ESC,   0, 1, 1, 1, 0, 0)
258 QUANTIZE_AND_ENCODE_BAND_COST_FUNC(NOISE, 0, 0, 0, 0, 1, 0)
259 QUANTIZE_AND_ENCODE_BAND_COST_FUNC(STEREO,0, 0, 0, 0, 0, 1)
260
261 static float (*const quantize_and_encode_band_cost_arr[])(
262                                 struct AACEncContext *s,
263                                 PutBitContext *pb, const float *in,
264                                 const float *scaled, int size, int scale_idx,
265                                 int cb, const float lambda, const float uplim,
266                                 int *bits) = {
267     quantize_and_encode_band_cost_ZERO,
268     quantize_and_encode_band_cost_SQUAD,
269     quantize_and_encode_band_cost_SQUAD,
270     quantize_and_encode_band_cost_UQUAD,
271     quantize_and_encode_band_cost_UQUAD,
272     quantize_and_encode_band_cost_SPAIR,
273     quantize_and_encode_band_cost_SPAIR,
274     quantize_and_encode_band_cost_UPAIR,
275     quantize_and_encode_band_cost_UPAIR,
276     quantize_and_encode_band_cost_UPAIR,
277     quantize_and_encode_band_cost_UPAIR,
278     quantize_and_encode_band_cost_ESC,
279     quantize_and_encode_band_cost_NONE,     /* CB 12 doesn't exist */
280     quantize_and_encode_band_cost_NOISE,
281     quantize_and_encode_band_cost_STEREO,
282     quantize_and_encode_band_cost_STEREO,
283 };
284
285 #define quantize_and_encode_band_cost(                                  \
286                                 s, pb, in, scaled, size, scale_idx, cb, \
287                                 lambda, uplim, bits)                    \
288     quantize_and_encode_band_cost_arr[cb](                              \
289                                 s, pb, in, scaled, size, scale_idx, cb, \
290                                 lambda, uplim, bits)
291
292 static float quantize_band_cost(struct AACEncContext *s, const float *in,
293                                 const float *scaled, int size, int scale_idx,
294                                 int cb, const float lambda, const float uplim,
295                                 int *bits)
296 {
297     return quantize_and_encode_band_cost(s, NULL, in, scaled, size, scale_idx,
298                                          cb, lambda, uplim, bits);
299 }
300
301 static void quantize_and_encode_band(struct AACEncContext *s, PutBitContext *pb,
302                                      const float *in, int size, int scale_idx,
303                                      int cb, const float lambda)
304 {
305     quantize_and_encode_band_cost(s, pb, in, NULL, size, scale_idx, cb, lambda,
306                                   INFINITY, NULL);
307 }
308
309 static float find_max_val(int group_len, int swb_size, const float *scaled) {
310     float maxval = 0.0f;
311     int w2, i;
312     for (w2 = 0; w2 < group_len; w2++) {
313         for (i = 0; i < swb_size; i++) {
314             maxval = FFMAX(maxval, scaled[w2*128+i]);
315         }
316     }
317     return maxval;
318 }
319
320 static int find_min_book(float maxval, int sf) {
321     float Q = ff_aac_pow2sf_tab[POW_SF2_ZERO - sf + SCALE_ONE_POS - SCALE_DIV_512];
322     float Q34 = sqrtf(Q * sqrtf(Q));
323     int qmaxval, cb;
324     qmaxval = maxval * Q34 + 0.4054f;
325     if      (qmaxval ==  0) cb = 0;
326     else if (qmaxval ==  1) cb = 1;
327     else if (qmaxval ==  2) cb = 3;
328     else if (qmaxval <=  4) cb = 5;
329     else if (qmaxval <=  7) cb = 7;
330     else if (qmaxval <= 12) cb = 9;
331     else                    cb = 11;
332     return cb;
333 }
334
335 /**
336  * structure used in optimal codebook search
337  */
338 typedef struct BandCodingPath {
339     int prev_idx; ///< pointer to the previous path point
340     float cost;   ///< path cost
341     int run;
342 } BandCodingPath;
343
344 /**
345  * Encode band info for single window group bands.
346  */
347 static void encode_window_bands_info(AACEncContext *s, SingleChannelElement *sce,
348                                      int win, int group_len, const float lambda)
349 {
350     BandCodingPath path[120][CB_TOT_ALL];
351     int w, swb, cb, start, size;
352     int i, j;
353     const int max_sfb  = sce->ics.max_sfb;
354     const int run_bits = sce->ics.num_windows == 1 ? 5 : 3;
355     const int run_esc  = (1 << run_bits) - 1;
356     int idx, ppos, count;
357     int stackrun[120], stackcb[120], stack_len;
358     float next_minrd = INFINITY;
359     int next_mincb = 0;
360
361     abs_pow34_v(s->scoefs, sce->coeffs, 1024);
362     start = win*128;
363     for (cb = 0; cb < CB_TOT_ALL; cb++) {
364         path[0][cb].cost     = 0.0f;
365         path[0][cb].prev_idx = -1;
366         path[0][cb].run      = 0;
367     }
368     for (swb = 0; swb < max_sfb; swb++) {
369         size = sce->ics.swb_sizes[swb];
370         if (sce->zeroes[win*16 + swb]) {
371             for (cb = 0; cb < CB_TOT_ALL; cb++) {
372                 path[swb+1][cb].prev_idx = cb;
373                 path[swb+1][cb].cost     = path[swb][cb].cost;
374                 path[swb+1][cb].run      = path[swb][cb].run + 1;
375             }
376         } else {
377             float minrd = next_minrd;
378             int mincb = next_mincb;
379             next_minrd = INFINITY;
380             next_mincb = 0;
381             for (cb = 0; cb < CB_TOT_ALL; cb++) {
382                 float cost_stay_here, cost_get_here;
383                 float rd = 0.0f;
384                 if (cb >= 12 && sce->band_type[win*16+swb] < aac_cb_out_map[cb] ||
385                     cb  < aac_cb_in_map[sce->band_type[win*16+swb]] && sce->band_type[win*16+swb] > aac_cb_out_map[cb]) {
386                     path[swb+1][cb].prev_idx = -1;
387                     path[swb+1][cb].cost     = INFINITY;
388                     path[swb+1][cb].run      = path[swb][cb].run + 1;
389                     continue;
390                 }
391                 for (w = 0; w < group_len; w++) {
392                     FFPsyBand *band = &s->psy.ch[s->cur_channel].psy_bands[(win+w)*16+swb];
393                     rd += quantize_band_cost(s, sce->coeffs + start + w*128,
394                                              s->scoefs + start + w*128, size,
395                                              sce->sf_idx[(win+w)*16+swb], aac_cb_out_map[cb],
396                                              lambda / band->threshold, INFINITY, NULL);
397                 }
398                 cost_stay_here = path[swb][cb].cost + rd;
399                 cost_get_here  = minrd              + rd + run_bits + 4;
400                 if (   run_value_bits[sce->ics.num_windows == 8][path[swb][cb].run]
401                     != run_value_bits[sce->ics.num_windows == 8][path[swb][cb].run+1])
402                     cost_stay_here += run_bits;
403                 if (cost_get_here < cost_stay_here) {
404                     path[swb+1][cb].prev_idx = mincb;
405                     path[swb+1][cb].cost     = cost_get_here;
406                     path[swb+1][cb].run      = 1;
407                 } else {
408                     path[swb+1][cb].prev_idx = cb;
409                     path[swb+1][cb].cost     = cost_stay_here;
410                     path[swb+1][cb].run      = path[swb][cb].run + 1;
411                 }
412                 if (path[swb+1][cb].cost < next_minrd) {
413                     next_minrd = path[swb+1][cb].cost;
414                     next_mincb = cb;
415                 }
416             }
417         }
418         start += sce->ics.swb_sizes[swb];
419     }
420
421     //convert resulting path from backward-linked list
422     stack_len = 0;
423     idx       = 0;
424     for (cb = 1; cb < CB_TOT_ALL; cb++)
425         if (path[max_sfb][cb].cost < path[max_sfb][idx].cost)
426             idx = cb;
427     ppos = max_sfb;
428     while (ppos > 0) {
429         av_assert1(idx >= 0);
430         cb = idx;
431         stackrun[stack_len] = path[ppos][cb].run;
432         stackcb [stack_len] = cb;
433         idx = path[ppos-path[ppos][cb].run+1][cb].prev_idx;
434         ppos -= path[ppos][cb].run;
435         stack_len++;
436     }
437     //perform actual band info encoding
438     start = 0;
439     for (i = stack_len - 1; i >= 0; i--) {
440         cb = aac_cb_out_map[stackcb[i]];
441         put_bits(&s->pb, 4, cb);
442         count = stackrun[i];
443         memset(sce->zeroes + win*16 + start, !cb, count);
444         //XXX: memset when band_type is also uint8_t
445         for (j = 0; j < count; j++) {
446             sce->band_type[win*16 + start] = cb;
447             start++;
448         }
449         while (count >= run_esc) {
450             put_bits(&s->pb, run_bits, run_esc);
451             count -= run_esc;
452         }
453         put_bits(&s->pb, run_bits, count);
454     }
455 }
456
457 static void codebook_trellis_rate(AACEncContext *s, SingleChannelElement *sce,
458                                   int win, int group_len, const float lambda)
459 {
460     BandCodingPath path[120][CB_TOT_ALL];
461     int w, swb, cb, start, size;
462     int i, j;
463     const int max_sfb  = sce->ics.max_sfb;
464     const int run_bits = sce->ics.num_windows == 1 ? 5 : 3;
465     const int run_esc  = (1 << run_bits) - 1;
466     int idx, ppos, count;
467     int stackrun[120], stackcb[120], stack_len;
468     float next_minbits = INFINITY;
469     int next_mincb = 0;
470
471     abs_pow34_v(s->scoefs, sce->coeffs, 1024);
472     start = win*128;
473     for (cb = 0; cb < CB_TOT_ALL; cb++) {
474         path[0][cb].cost     = run_bits+4;
475         path[0][cb].prev_idx = -1;
476         path[0][cb].run      = 0;
477     }
478     for (swb = 0; swb < max_sfb; swb++) {
479         size = sce->ics.swb_sizes[swb];
480         if (sce->zeroes[win*16 + swb]) {
481             float cost_stay_here = path[swb][0].cost;
482             float cost_get_here  = next_minbits + run_bits + 4;
483             if (   run_value_bits[sce->ics.num_windows == 8][path[swb][0].run]
484                 != run_value_bits[sce->ics.num_windows == 8][path[swb][0].run+1])
485                 cost_stay_here += run_bits;
486             if (cost_get_here < cost_stay_here) {
487                 path[swb+1][0].prev_idx = next_mincb;
488                 path[swb+1][0].cost     = cost_get_here;
489                 path[swb+1][0].run      = 1;
490             } else {
491                 path[swb+1][0].prev_idx = 0;
492                 path[swb+1][0].cost     = cost_stay_here;
493                 path[swb+1][0].run      = path[swb][0].run + 1;
494             }
495             next_minbits = path[swb+1][0].cost;
496             next_mincb = 0;
497             for (cb = 1; cb < CB_TOT_ALL; cb++) {
498                 path[swb+1][cb].cost = 61450;
499                 path[swb+1][cb].prev_idx = -1;
500                 path[swb+1][cb].run = 0;
501             }
502         } else {
503             float minbits = next_minbits;
504             int mincb = next_mincb;
505             int startcb = sce->band_type[win*16+swb];
506             startcb = aac_cb_in_map[startcb];
507             next_minbits = INFINITY;
508             next_mincb = 0;
509             for (cb = 0; cb < startcb; cb++) {
510                 path[swb+1][cb].cost = 61450;
511                 path[swb+1][cb].prev_idx = -1;
512                 path[swb+1][cb].run = 0;
513             }
514             for (cb = startcb; cb < CB_TOT_ALL; cb++) {
515                 float cost_stay_here, cost_get_here;
516                 float bits = 0.0f;
517                 if (cb >= 12 && sce->band_type[win*16+swb] != aac_cb_out_map[cb]) {
518                     path[swb+1][cb].cost = 61450;
519                     path[swb+1][cb].prev_idx = -1;
520                     path[swb+1][cb].run = 0;
521                     continue;
522                 }
523                 for (w = 0; w < group_len; w++) {
524                     bits += quantize_band_cost(s, sce->coeffs + start + w*128,
525                                                s->scoefs + start + w*128, size,
526                                                sce->sf_idx[(win+w)*16+swb],
527                                                aac_cb_out_map[cb],
528                                                0, INFINITY, NULL);
529                 }
530                 cost_stay_here = path[swb][cb].cost + bits;
531                 cost_get_here  = minbits            + bits + run_bits + 4;
532                 if (   run_value_bits[sce->ics.num_windows == 8][path[swb][cb].run]
533                     != run_value_bits[sce->ics.num_windows == 8][path[swb][cb].run+1])
534                     cost_stay_here += run_bits;
535                 if (cost_get_here < cost_stay_here) {
536                     path[swb+1][cb].prev_idx = mincb;
537                     path[swb+1][cb].cost     = cost_get_here;
538                     path[swb+1][cb].run      = 1;
539                 } else {
540                     path[swb+1][cb].prev_idx = cb;
541                     path[swb+1][cb].cost     = cost_stay_here;
542                     path[swb+1][cb].run      = path[swb][cb].run + 1;
543                 }
544                 if (path[swb+1][cb].cost < next_minbits) {
545                     next_minbits = path[swb+1][cb].cost;
546                     next_mincb = cb;
547                 }
548             }
549         }
550         start += sce->ics.swb_sizes[swb];
551     }
552
553     //convert resulting path from backward-linked list
554     stack_len = 0;
555     idx       = 0;
556     for (cb = 1; cb < CB_TOT_ALL; cb++)
557         if (path[max_sfb][cb].cost < path[max_sfb][idx].cost)
558             idx = cb;
559     ppos = max_sfb;
560     while (ppos > 0) {
561         av_assert1(idx >= 0);
562         cb = idx;
563         stackrun[stack_len] = path[ppos][cb].run;
564         stackcb [stack_len] = cb;
565         idx = path[ppos-path[ppos][cb].run+1][cb].prev_idx;
566         ppos -= path[ppos][cb].run;
567         stack_len++;
568     }
569     //perform actual band info encoding
570     start = 0;
571     for (i = stack_len - 1; i >= 0; i--) {
572         cb = aac_cb_out_map[stackcb[i]];
573         put_bits(&s->pb, 4, cb);
574         count = stackrun[i];
575         memset(sce->zeroes + win*16 + start, !cb, count);
576         //XXX: memset when band_type is also uint8_t
577         for (j = 0; j < count; j++) {
578             sce->band_type[win*16 + start] = cb;
579             start++;
580         }
581         while (count >= run_esc) {
582             put_bits(&s->pb, run_bits, run_esc);
583             count -= run_esc;
584         }
585         put_bits(&s->pb, run_bits, count);
586     }
587 }
588
589 /** Return the minimum scalefactor where the quantized coef does not clip. */
590 static av_always_inline uint8_t coef2minsf(float coef) {
591     return av_clip_uint8(log2f(coef)*4 - 69 + SCALE_ONE_POS - SCALE_DIV_512);
592 }
593
594 /** Return the maximum scalefactor where the quantized coef is not zero. */
595 static av_always_inline uint8_t coef2maxsf(float coef) {
596     return av_clip_uint8(log2f(coef)*4 +  6 + SCALE_ONE_POS - SCALE_DIV_512);
597 }
598
599 typedef struct TrellisPath {
600     float cost;
601     int prev;
602 } TrellisPath;
603
604 #define TRELLIS_STAGES 121
605 #define TRELLIS_STATES (SCALE_MAX_DIFF+1)
606
607 static void set_special_band_scalefactors(AACEncContext *s, SingleChannelElement *sce)
608 {
609     int w, g, start = 0;
610     int minscaler_n = sce->sf_idx[0], minscaler_i = sce->sf_idx[0];
611     int bands = 0;
612
613     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
614         start = 0;
615         for (g = 0;  g < sce->ics.num_swb; g++) {
616             if (sce->band_type[w*16+g] == INTENSITY_BT || sce->band_type[w*16+g] == INTENSITY_BT2) {
617                 sce->sf_idx[w*16+g] = av_clip(ceilf(log2f(sce->is_ener[w*16+g])*2), -155, 100);
618                 minscaler_i = FFMIN(minscaler_i, sce->sf_idx[w*16+g]);
619                 bands++;
620             } else if (sce->band_type[w*16+g] == NOISE_BT) {
621                 sce->sf_idx[w*16+g] = av_clip(4+log2f(sce->pns_ener[w*16+g])*2, -100, 155);
622                 minscaler_n = FFMIN(minscaler_n, sce->sf_idx[w*16+g]);
623                 bands++;
624             }
625             start += sce->ics.swb_sizes[g];
626         }
627     }
628
629     if (!bands)
630         return;
631
632     /* Clip the scalefactor indices */
633     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
634         for (g = 0;  g < sce->ics.num_swb; g++) {
635             if (sce->band_type[w*16+g] == INTENSITY_BT || sce->band_type[w*16+g] == INTENSITY_BT2) {
636                 sce->sf_idx[w*16+g] = av_clip(sce->sf_idx[w*16+g], minscaler_i, minscaler_i + SCALE_MAX_DIFF);
637             } else if (sce->band_type[w*16+g] == NOISE_BT) {
638                 sce->sf_idx[w*16+g] = av_clip(sce->sf_idx[w*16+g], minscaler_n, minscaler_n + SCALE_MAX_DIFF);
639             }
640         }
641     }
642 }
643
644 static void search_for_quantizers_anmr(AVCodecContext *avctx, AACEncContext *s,
645                                        SingleChannelElement *sce,
646                                        const float lambda)
647 {
648     int q, w, w2, g, start = 0;
649     int i, j;
650     int idx;
651     TrellisPath paths[TRELLIS_STAGES][TRELLIS_STATES];
652     int bandaddr[TRELLIS_STAGES];
653     int minq;
654     float mincost;
655     float q0f = FLT_MAX, q1f = 0.0f, qnrgf = 0.0f;
656     int q0, q1, qcnt = 0;
657
658     for (i = 0; i < 1024; i++) {
659         float t = fabsf(sce->coeffs[i]);
660         if (t > 0.0f) {
661             q0f = FFMIN(q0f, t);
662             q1f = FFMAX(q1f, t);
663             qnrgf += t*t;
664             qcnt++;
665         }
666     }
667
668     if (!qcnt) {
669         memset(sce->sf_idx, 0, sizeof(sce->sf_idx));
670         memset(sce->zeroes, 1, sizeof(sce->zeroes));
671         return;
672     }
673
674     //minimum scalefactor index is when minimum nonzero coefficient after quantizing is not clipped
675     q0 = coef2minsf(q0f);
676     //maximum scalefactor index is when maximum coefficient after quantizing is still not zero
677     q1 = coef2maxsf(q1f);
678     if (q1 - q0 > 60) {
679         int q0low  = q0;
680         int q1high = q1;
681         //minimum scalefactor index is when maximum nonzero coefficient after quantizing is not clipped
682         int qnrg = av_clip_uint8(log2f(sqrtf(qnrgf/qcnt))*4 - 31 + SCALE_ONE_POS - SCALE_DIV_512);
683         q1 = qnrg + 30;
684         q0 = qnrg - 30;
685         if (q0 < q0low) {
686             q1 += q0low - q0;
687             q0  = q0low;
688         } else if (q1 > q1high) {
689             q0 -= q1 - q1high;
690             q1  = q1high;
691         }
692     }
693
694     for (i = 0; i < TRELLIS_STATES; i++) {
695         paths[0][i].cost    = 0.0f;
696         paths[0][i].prev    = -1;
697     }
698     for (j = 1; j < TRELLIS_STAGES; j++) {
699         for (i = 0; i < TRELLIS_STATES; i++) {
700             paths[j][i].cost    = INFINITY;
701             paths[j][i].prev    = -2;
702         }
703     }
704     idx = 1;
705     abs_pow34_v(s->scoefs, sce->coeffs, 1024);
706     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
707         start = w*128;
708         for (g = 0; g < sce->ics.num_swb; g++) {
709             const float *coefs = sce->coeffs + start;
710             float qmin, qmax;
711             int nz = 0;
712
713             bandaddr[idx] = w * 16 + g;
714             qmin = INT_MAX;
715             qmax = 0.0f;
716             for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
717                 FFPsyBand *band = &s->psy.ch[s->cur_channel].psy_bands[(w+w2)*16+g];
718                 if (band->energy <= band->threshold || band->threshold == 0.0f) {
719                     sce->zeroes[(w+w2)*16+g] = 1;
720                     continue;
721                 }
722                 sce->zeroes[(w+w2)*16+g] = 0;
723                 nz = 1;
724                 for (i = 0; i < sce->ics.swb_sizes[g]; i++) {
725                     float t = fabsf(coefs[w2*128+i]);
726                     if (t > 0.0f)
727                         qmin = FFMIN(qmin, t);
728                     qmax = FFMAX(qmax, t);
729                 }
730             }
731             if (nz) {
732                 int minscale, maxscale;
733                 float minrd = INFINITY;
734                 float maxval;
735                 //minimum scalefactor index is when minimum nonzero coefficient after quantizing is not clipped
736                 minscale = coef2minsf(qmin);
737                 //maximum scalefactor index is when maximum coefficient after quantizing is still not zero
738                 maxscale = coef2maxsf(qmax);
739                 minscale = av_clip(minscale - q0, 0, TRELLIS_STATES - 1);
740                 maxscale = av_clip(maxscale - q0, 0, TRELLIS_STATES);
741                 maxval = find_max_val(sce->ics.group_len[w], sce->ics.swb_sizes[g], s->scoefs+start);
742                 for (q = minscale; q < maxscale; q++) {
743                     float dist = 0;
744                     int cb = find_min_book(maxval, sce->sf_idx[w*16+g]);
745                     for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
746                         FFPsyBand *band = &s->psy.ch[s->cur_channel].psy_bands[(w+w2)*16+g];
747                         dist += quantize_band_cost(s, coefs + w2*128, s->scoefs + start + w2*128, sce->ics.swb_sizes[g],
748                                                    q + q0, cb, lambda / band->threshold, INFINITY, NULL);
749                     }
750                     minrd = FFMIN(minrd, dist);
751
752                     for (i = 0; i < q1 - q0; i++) {
753                         float cost;
754                         cost = paths[idx - 1][i].cost + dist
755                                + ff_aac_scalefactor_bits[q - i + SCALE_DIFF_ZERO];
756                         if (cost < paths[idx][q].cost) {
757                             paths[idx][q].cost    = cost;
758                             paths[idx][q].prev    = i;
759                         }
760                     }
761                 }
762             } else {
763                 for (q = 0; q < q1 - q0; q++) {
764                     paths[idx][q].cost = paths[idx - 1][q].cost + 1;
765                     paths[idx][q].prev = q;
766                 }
767             }
768             sce->zeroes[w*16+g] = !nz;
769             start += sce->ics.swb_sizes[g];
770             idx++;
771         }
772     }
773     idx--;
774     mincost = paths[idx][0].cost;
775     minq    = 0;
776     for (i = 1; i < TRELLIS_STATES; i++) {
777         if (paths[idx][i].cost < mincost) {
778             mincost = paths[idx][i].cost;
779             minq = i;
780         }
781     }
782     while (idx) {
783         sce->sf_idx[bandaddr[idx]] = minq + q0;
784         minq = paths[idx][minq].prev;
785         idx--;
786     }
787     //set the same quantizers inside window groups
788     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w])
789         for (g = 0;  g < sce->ics.num_swb; g++)
790             for (w2 = 1; w2 < sce->ics.group_len[w]; w2++)
791                 sce->sf_idx[(w+w2)*16+g] = sce->sf_idx[w*16+g];
792 }
793
794 /**
795  * two-loop quantizers search taken from ISO 13818-7 Appendix C
796  */
797 static void search_for_quantizers_twoloop(AVCodecContext *avctx,
798                                           AACEncContext *s,
799                                           SingleChannelElement *sce,
800                                           const float lambda)
801 {
802     int start = 0, i, w, w2, g;
803     int destbits = avctx->bit_rate * 1024.0 / avctx->sample_rate / avctx->channels * (lambda / 120.f);
804     float dists[128] = { 0 }, uplims[128] = { 0 };
805     float maxvals[128];
806     int fflag, minscaler;
807     int its  = 0;
808     int allz = 0;
809     float minthr = INFINITY;
810
811     // for values above this the decoder might end up in an endless loop
812     // due to always having more bits than what can be encoded.
813     destbits = FFMIN(destbits, 5800);
814     //XXX: some heuristic to determine initial quantizers will reduce search time
815     //determine zero bands and upper limits
816     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
817         for (g = 0;  g < sce->ics.num_swb; g++) {
818             int nz = 0;
819             float uplim = 0.0f, energy = 0.0f;
820             for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
821                 FFPsyBand *band = &s->psy.ch[s->cur_channel].psy_bands[(w+w2)*16+g];
822                 uplim  += band->threshold;
823                 energy += band->energy;
824                 if (band->energy <= band->threshold || band->threshold == 0.0f) {
825                     sce->zeroes[(w+w2)*16+g] = 1;
826                     continue;
827                 }
828                 nz = 1;
829             }
830             uplims[w*16+g] = uplim *512;
831             sce->zeroes[w*16+g] = !nz;
832             if (nz)
833                 minthr = FFMIN(minthr, uplim);
834             allz |= nz;
835         }
836     }
837     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
838         for (g = 0;  g < sce->ics.num_swb; g++) {
839             if (sce->zeroes[w*16+g]) {
840                 sce->sf_idx[w*16+g] = SCALE_ONE_POS;
841                 continue;
842             }
843             sce->sf_idx[w*16+g] = SCALE_ONE_POS + FFMIN(log2f(uplims[w*16+g]/minthr)*4,59);
844         }
845     }
846
847     if (!allz)
848         return;
849     abs_pow34_v(s->scoefs, sce->coeffs, 1024);
850
851     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
852         start = w*128;
853         for (g = 0;  g < sce->ics.num_swb; g++) {
854             const float *scaled = s->scoefs + start;
855             maxvals[w*16+g] = find_max_val(sce->ics.group_len[w], sce->ics.swb_sizes[g], scaled);
856             start += sce->ics.swb_sizes[g];
857         }
858     }
859
860     //perform two-loop search
861     //outer loop - improve quality
862     do {
863         int tbits, qstep;
864         minscaler = sce->sf_idx[0];
865         //inner loop - quantize spectrum to fit into given number of bits
866         qstep = its ? 1 : 32;
867         do {
868             int prev = -1;
869             tbits = 0;
870             for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
871                 start = w*128;
872                 for (g = 0;  g < sce->ics.num_swb; g++) {
873                     const float *coefs = sce->coeffs + start;
874                     const float *scaled = s->scoefs + start;
875                     int bits = 0;
876                     int cb;
877                     float dist = 0.0f;
878
879                     if (sce->zeroes[w*16+g] || sce->sf_idx[w*16+g] >= 218) {
880                         start += sce->ics.swb_sizes[g];
881                         continue;
882                     }
883                     minscaler = FFMIN(minscaler, sce->sf_idx[w*16+g]);
884                     cb = find_min_book(maxvals[w*16+g], sce->sf_idx[w*16+g]);
885                     for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
886                         int b;
887                         dist += quantize_band_cost(s, coefs + w2*128,
888                                                    scaled + w2*128,
889                                                    sce->ics.swb_sizes[g],
890                                                    sce->sf_idx[w*16+g],
891                                                    cb,
892                                                    1.0f,
893                                                    INFINITY,
894                                                    &b);
895                         bits += b;
896                     }
897                     dists[w*16+g] = dist - bits;
898                     if (prev != -1) {
899                         bits += ff_aac_scalefactor_bits[sce->sf_idx[w*16+g] - prev + SCALE_DIFF_ZERO];
900                     }
901                     tbits += bits;
902                     start += sce->ics.swb_sizes[g];
903                     prev = sce->sf_idx[w*16+g];
904                 }
905             }
906             if (tbits > destbits) {
907                 for (i = 0; i < 128; i++)
908                     if (sce->sf_idx[i] < 218 - qstep)
909                         sce->sf_idx[i] += qstep;
910             } else {
911                 for (i = 0; i < 128; i++)
912                     if (sce->sf_idx[i] > 60 - qstep)
913                         sce->sf_idx[i] -= qstep;
914             }
915             qstep >>= 1;
916             if (!qstep && tbits > destbits*1.02 && sce->sf_idx[0] < 217)
917                 qstep = 1;
918         } while (qstep);
919
920         fflag = 0;
921         minscaler = av_clip(minscaler, 60, 255 - SCALE_MAX_DIFF);
922
923         for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
924             for (g = 0; g < sce->ics.num_swb; g++) {
925                 int prevsc = sce->sf_idx[w*16+g];
926                 if (dists[w*16+g] > uplims[w*16+g] && sce->sf_idx[w*16+g] > 60) {
927                     if (find_min_book(maxvals[w*16+g], sce->sf_idx[w*16+g]-1))
928                         sce->sf_idx[w*16+g]--;
929                     else //Try to make sure there is some energy in every band
930                         sce->sf_idx[w*16+g]-=2;
931                 }
932                 sce->sf_idx[w*16+g] = av_clip(sce->sf_idx[w*16+g], minscaler, minscaler + SCALE_MAX_DIFF);
933                 sce->sf_idx[w*16+g] = FFMIN(sce->sf_idx[w*16+g], 219);
934                 if (sce->sf_idx[w*16+g] != prevsc)
935                     fflag = 1;
936                 sce->band_type[w*16+g] = find_min_book(maxvals[w*16+g], sce->sf_idx[w*16+g]);
937             }
938         }
939         its++;
940     } while (fflag && its < 10);
941 }
942
943 static void search_for_quantizers_faac(AVCodecContext *avctx, AACEncContext *s,
944                                        SingleChannelElement *sce,
945                                        const float lambda)
946 {
947     int start = 0, i, w, w2, g;
948     float uplim[128], maxq[128];
949     int minq, maxsf;
950     float distfact = ((sce->ics.num_windows > 1) ? 85.80 : 147.84) / lambda;
951     int last = 0, lastband = 0, curband = 0;
952     float avg_energy = 0.0;
953     if (sce->ics.num_windows == 1) {
954         start = 0;
955         for (i = 0; i < 1024; i++) {
956             if (i - start >= sce->ics.swb_sizes[curband]) {
957                 start += sce->ics.swb_sizes[curband];
958                 curband++;
959             }
960             if (sce->coeffs[i]) {
961                 avg_energy += sce->coeffs[i] * sce->coeffs[i];
962                 last = i;
963                 lastband = curband;
964             }
965         }
966     } else {
967         for (w = 0; w < 8; w++) {
968             const float *coeffs = sce->coeffs + w*128;
969             curband = start = 0;
970             for (i = 0; i < 128; i++) {
971                 if (i - start >= sce->ics.swb_sizes[curband]) {
972                     start += sce->ics.swb_sizes[curband];
973                     curband++;
974                 }
975                 if (coeffs[i]) {
976                     avg_energy += coeffs[i] * coeffs[i];
977                     last = FFMAX(last, i);
978                     lastband = FFMAX(lastband, curband);
979                 }
980             }
981         }
982     }
983     last++;
984     avg_energy /= last;
985     if (avg_energy == 0.0f) {
986         for (i = 0; i < FF_ARRAY_ELEMS(sce->sf_idx); i++)
987             sce->sf_idx[i] = SCALE_ONE_POS;
988         return;
989     }
990     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
991         start = w*128;
992         for (g = 0; g < sce->ics.num_swb; g++) {
993             float *coefs   = sce->coeffs + start;
994             const int size = sce->ics.swb_sizes[g];
995             int start2 = start, end2 = start + size, peakpos = start;
996             float maxval = -1, thr = 0.0f, t;
997             maxq[w*16+g] = 0.0f;
998             if (g > lastband) {
999                 maxq[w*16+g] = 0.0f;
1000                 start += size;
1001                 for (w2 = 0; w2 < sce->ics.group_len[w]; w2++)
1002                     memset(coefs + w2*128, 0, sizeof(coefs[0])*size);
1003                 continue;
1004             }
1005             for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
1006                 for (i = 0; i < size; i++) {
1007                     float t = coefs[w2*128+i]*coefs[w2*128+i];
1008                     maxq[w*16+g] = FFMAX(maxq[w*16+g], fabsf(coefs[w2*128 + i]));
1009                     thr += t;
1010                     if (sce->ics.num_windows == 1 && maxval < t) {
1011                         maxval  = t;
1012                         peakpos = start+i;
1013                     }
1014                 }
1015             }
1016             if (sce->ics.num_windows == 1) {
1017                 start2 = FFMAX(peakpos - 2, start2);
1018                 end2   = FFMIN(peakpos + 3, end2);
1019             } else {
1020                 start2 -= start;
1021                 end2   -= start;
1022             }
1023             start += size;
1024             thr = pow(thr / (avg_energy * (end2 - start2)), 0.3 + 0.1*(lastband - g) / lastband);
1025             t   = 1.0 - (1.0 * start2 / last);
1026             uplim[w*16+g] = distfact / (1.4 * thr + t*t*t + 0.075);
1027         }
1028     }
1029     memset(sce->sf_idx, 0, sizeof(sce->sf_idx));
1030     abs_pow34_v(s->scoefs, sce->coeffs, 1024);
1031     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
1032         start = w*128;
1033         for (g = 0;  g < sce->ics.num_swb; g++) {
1034             const float *coefs  = sce->coeffs + start;
1035             const float *scaled = s->scoefs   + start;
1036             const int size      = sce->ics.swb_sizes[g];
1037             int scf, prev_scf, step;
1038             int min_scf = -1, max_scf = 256;
1039             float curdiff;
1040             if (maxq[w*16+g] < 21.544) {
1041                 sce->zeroes[w*16+g] = 1;
1042                 start += size;
1043                 continue;
1044             }
1045             sce->zeroes[w*16+g] = 0;
1046             scf  = prev_scf = av_clip(SCALE_ONE_POS - SCALE_DIV_512 - log2f(1/maxq[w*16+g])*16/3, 60, 218);
1047             for (;;) {
1048                 float dist = 0.0f;
1049                 int quant_max;
1050
1051                 for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
1052                     int b;
1053                     dist += quantize_band_cost(s, coefs + w2*128,
1054                                                scaled + w2*128,
1055                                                sce->ics.swb_sizes[g],
1056                                                scf,
1057                                                ESC_BT,
1058                                                lambda,
1059                                                INFINITY,
1060                                                &b);
1061                     dist -= b;
1062                 }
1063                 dist *= 1.0f / 512.0f / lambda;
1064                 quant_max = quant(maxq[w*16+g], ff_aac_pow2sf_tab[POW_SF2_ZERO - scf + SCALE_ONE_POS - SCALE_DIV_512]);
1065                 if (quant_max >= 8191) { // too much, return to the previous quantizer
1066                     sce->sf_idx[w*16+g] = prev_scf;
1067                     break;
1068                 }
1069                 prev_scf = scf;
1070                 curdiff = fabsf(dist - uplim[w*16+g]);
1071                 if (curdiff <= 1.0f)
1072                     step = 0;
1073                 else
1074                     step = log2f(curdiff);
1075                 if (dist > uplim[w*16+g])
1076                     step = -step;
1077                 scf += step;
1078                 scf = av_clip_uint8(scf);
1079                 step = scf - prev_scf;
1080                 if (FFABS(step) <= 1 || (step > 0 && scf >= max_scf) || (step < 0 && scf <= min_scf)) {
1081                     sce->sf_idx[w*16+g] = av_clip(scf, min_scf, max_scf);
1082                     break;
1083                 }
1084                 if (step > 0)
1085                     min_scf = prev_scf;
1086                 else
1087                     max_scf = prev_scf;
1088             }
1089             start += size;
1090         }
1091     }
1092     minq = sce->sf_idx[0] ? sce->sf_idx[0] : INT_MAX;
1093     for (i = 1; i < 128; i++) {
1094         if (!sce->sf_idx[i])
1095             sce->sf_idx[i] = sce->sf_idx[i-1];
1096         else
1097             minq = FFMIN(minq, sce->sf_idx[i]);
1098     }
1099     if (minq == INT_MAX)
1100         minq = 0;
1101     minq = FFMIN(minq, SCALE_MAX_POS);
1102     maxsf = FFMIN(minq + SCALE_MAX_DIFF, SCALE_MAX_POS);
1103     for (i = 126; i >= 0; i--) {
1104         if (!sce->sf_idx[i])
1105             sce->sf_idx[i] = sce->sf_idx[i+1];
1106         sce->sf_idx[i] = av_clip(sce->sf_idx[i], minq, maxsf);
1107     }
1108 }
1109
1110 static void search_for_quantizers_fast(AVCodecContext *avctx, AACEncContext *s,
1111                                        SingleChannelElement *sce,
1112                                        const float lambda)
1113 {
1114     int i, w, w2, g;
1115     int minq = 255;
1116
1117     memset(sce->sf_idx, 0, sizeof(sce->sf_idx));
1118     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
1119         for (g = 0; g < sce->ics.num_swb; g++) {
1120             for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
1121                 FFPsyBand *band = &s->psy.ch[s->cur_channel].psy_bands[(w+w2)*16+g];
1122                 if (band->energy <= band->threshold) {
1123                     sce->sf_idx[(w+w2)*16+g] = 218;
1124                     sce->zeroes[(w+w2)*16+g] = 1;
1125                 } else {
1126                     sce->sf_idx[(w+w2)*16+g] = av_clip(SCALE_ONE_POS - SCALE_DIV_512 + log2f(band->threshold), 80, 218);
1127                     sce->zeroes[(w+w2)*16+g] = 0;
1128                 }
1129                 minq = FFMIN(minq, sce->sf_idx[(w+w2)*16+g]);
1130             }
1131         }
1132     }
1133     for (i = 0; i < 128; i++) {
1134         sce->sf_idx[i] = 140;
1135         //av_clip(sce->sf_idx[i], minq, minq + SCALE_MAX_DIFF - 1);
1136     }
1137     //set the same quantizers inside window groups
1138     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w])
1139         for (g = 0;  g < sce->ics.num_swb; g++)
1140             for (w2 = 1; w2 < sce->ics.group_len[w]; w2++)
1141                 sce->sf_idx[(w+w2)*16+g] = sce->sf_idx[w*16+g];
1142 }
1143
1144 static void search_for_pns(AACEncContext *s, AVCodecContext *avctx, SingleChannelElement *sce,
1145                            const float lambda)
1146 {
1147     int start = 0, w, w2, g;
1148     const float freq_mult = avctx->sample_rate/(1024.0f/sce->ics.num_windows)/2.0f;
1149     const float spread_threshold = NOISE_SPREAD_THRESHOLD*(lambda/120.f);
1150     const float thr_mult = NOISE_LAMBDA_NUMERATOR/lambda;
1151
1152     /* Coders !twoloop don't reset the band_types */
1153     for (w = 0; w < 128; w++)
1154         if (sce->band_type[w] == NOISE_BT)
1155             sce->band_type[w] = 0;
1156
1157     for (w = 0; w < sce->ics.num_windows; w += sce->ics.group_len[w]) {
1158         start = 0;
1159         for (g = 0;  g < sce->ics.num_swb; g++) {
1160             if (start*freq_mult > NOISE_LOW_LIMIT*(lambda/170.0f)) {
1161                 float energy = 0.0f, threshold = 0.0f, spread = 0.0f;
1162                 for (w2 = 0; w2 < sce->ics.group_len[w]; w2++) {
1163                     FFPsyBand *band = &s->psy.ch[s->cur_channel+0].psy_bands[(w+w2)*16+g];
1164                     energy += band->energy;
1165                     threshold += band->threshold;
1166                     spread += band->spread;
1167                 }
1168                 if (spread > spread_threshold*sce->ics.group_len[w] &&
1169                     ((sce->zeroes[w*16+g] && energy >= threshold) ||
1170                     energy < threshold*thr_mult*sce->ics.group_len[w])) {
1171                     sce->band_type[w*16+g] = NOISE_BT;
1172                     sce->pns_ener[w*16+g] = energy / sce->ics.group_len[w];
1173                     sce->zeroes[w*16+g] = 0;
1174                 }
1175             }
1176             start += sce->ics.swb_sizes[g];
1177         }
1178     }
1179 }
1180
1181 static void search_for_ms(AACEncContext *s, ChannelElement *cpe,
1182                           const float lambda)
1183 {
1184     int start = 0, i, w, w2, g;
1185     float M[128], S[128];
1186     float *L34 = s->scoefs, *R34 = s->scoefs + 128, *M34 = s->scoefs + 128*2, *S34 = s->scoefs + 128*3;
1187     SingleChannelElement *sce0 = &cpe->ch[0];
1188     SingleChannelElement *sce1 = &cpe->ch[1];
1189     if (!cpe->common_window)
1190         return;
1191     for (w = 0; w < sce0->ics.num_windows; w += sce0->ics.group_len[w]) {
1192         start = 0;
1193         for (g = 0;  g < sce0->ics.num_swb; g++) {
1194             if (!cpe->ch[0].zeroes[w*16+g] && !cpe->ch[1].zeroes[w*16+g]) {
1195                 float dist1 = 0.0f, dist2 = 0.0f;
1196                 for (w2 = 0; w2 < sce0->ics.group_len[w]; w2++) {
1197                     FFPsyBand *band0 = &s->psy.ch[s->cur_channel+0].psy_bands[(w+w2)*16+g];
1198                     FFPsyBand *band1 = &s->psy.ch[s->cur_channel+1].psy_bands[(w+w2)*16+g];
1199                     float minthr = FFMIN(band0->threshold, band1->threshold);
1200                     float maxthr = FFMAX(band0->threshold, band1->threshold);
1201                     for (i = 0; i < sce0->ics.swb_sizes[g]; i++) {
1202                         M[i] = (sce0->pcoeffs[start+(w+w2)*128+i]
1203                               + sce1->pcoeffs[start+(w+w2)*128+i]) * 0.5;
1204                         S[i] =  M[i]
1205                               - sce1->pcoeffs[start+(w+w2)*128+i];
1206                     }
1207                     abs_pow34_v(L34, sce0->coeffs+start+(w+w2)*128, sce0->ics.swb_sizes[g]);
1208                     abs_pow34_v(R34, sce1->coeffs+start+(w+w2)*128, sce0->ics.swb_sizes[g]);
1209                     abs_pow34_v(M34, M,                         sce0->ics.swb_sizes[g]);
1210                     abs_pow34_v(S34, S,                         sce0->ics.swb_sizes[g]);
1211                     dist1 += quantize_band_cost(s, sce0->coeffs + start + (w+w2)*128,
1212                                                 L34,
1213                                                 sce0->ics.swb_sizes[g],
1214                                                 sce0->sf_idx[(w+w2)*16+g],
1215                                                 sce0->band_type[(w+w2)*16+g],
1216                                                 lambda / band0->threshold, INFINITY, NULL);
1217                     dist1 += quantize_band_cost(s, sce1->coeffs + start + (w+w2)*128,
1218                                                 R34,
1219                                                 sce1->ics.swb_sizes[g],
1220                                                 sce1->sf_idx[(w+w2)*16+g],
1221                                                 sce1->band_type[(w+w2)*16+g],
1222                                                 lambda / band1->threshold, INFINITY, NULL);
1223                     dist2 += quantize_band_cost(s, M,
1224                                                 M34,
1225                                                 sce0->ics.swb_sizes[g],
1226                                                 sce0->sf_idx[(w+w2)*16+g],
1227                                                 sce0->band_type[(w+w2)*16+g],
1228                                                 lambda / maxthr, INFINITY, NULL);
1229                     dist2 += quantize_band_cost(s, S,
1230                                                 S34,
1231                                                 sce1->ics.swb_sizes[g],
1232                                                 sce1->sf_idx[(w+w2)*16+g],
1233                                                 sce1->band_type[(w+w2)*16+g],
1234                                                 lambda / minthr, INFINITY, NULL);
1235                 }
1236                 cpe->ms_mask[w*16+g] = dist2 < dist1;
1237             }
1238             start += sce0->ics.swb_sizes[g];
1239         }
1240     }
1241 }
1242
1243 AACCoefficientsEncoder ff_aac_coders[AAC_CODER_NB] = {
1244     [AAC_CODER_FAAC] = {
1245         search_for_quantizers_faac,
1246         encode_window_bands_info,
1247         quantize_and_encode_band,
1248         set_special_band_scalefactors,
1249         search_for_pns,
1250         search_for_ms,
1251     },
1252     [AAC_CODER_ANMR] = {
1253         search_for_quantizers_anmr,
1254         encode_window_bands_info,
1255         quantize_and_encode_band,
1256         set_special_band_scalefactors,
1257         search_for_pns,
1258         search_for_ms,
1259     },
1260     [AAC_CODER_TWOLOOP] = {
1261         search_for_quantizers_twoloop,
1262         codebook_trellis_rate,
1263         quantize_and_encode_band,
1264         set_special_band_scalefactors,
1265         search_for_pns,
1266         search_for_ms,
1267     },
1268     [AAC_CODER_FAST] = {
1269         search_for_quantizers_fast,
1270         encode_window_bands_info,
1271         quantize_and_encode_band,
1272         set_special_band_scalefactors,
1273         search_for_pns,
1274         search_for_ms,
1275     },
1276 };