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