Coverage Report

Created: 2026-07-24 06:21

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/proc/self/cwd/libfaad/ps_dec.c
Line
Count
Source
1
/*
2
** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
3
** Copyright (C) 2003-2005 M. Bakker, Nero AG, http://www.nero.com
4
**
5
** This program is free software; you can redistribute it and/or modify
6
** it under the terms of the GNU General Public License as published by
7
** the Free Software Foundation; either version 2 of the License, or
8
** (at your option) any later version.
9
**
10
** This program is distributed in the hope that it will be useful,
11
** but WITHOUT ANY WARRANTY; without even the implied warranty of
12
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13
** GNU General Public License for more details.
14
**
15
** You should have received a copy of the GNU General Public License
16
** along with this program; if not, write to the Free Software
17
** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
18
**
19
** Any non-GPL usage of this software or parts of this software is strictly
20
** forbidden.
21
**
22
** The "appropriate copyright message" mentioned in section 2c of the GPLv2
23
** must read: "Code from FAAD2 is copyright (c) Nero AG, www.nero.com"
24
**
25
** Commercial non-GPL licensing of this software is possible.
26
** For more info contact Nero AG through Mpeg4AAClicense@nero.com.
27
**
28
** $Id: ps_dec.c,v 1.16 2009/01/26 22:32:31 menno Exp $
29
**/
30
31
#include "common.h"
32
33
#ifdef PS_DEC
34
35
#include <stdlib.h>
36
#include <stdio.h>
37
#include "ps_dec.h"
38
#include "ps_tables.h"
39
40
/* constants */
41
13.8M
#define NEGATE_IPD_MASK            (0x1000)
42
104k
#define DECAY_SLOPE                FRAC_CONST(0.05)
43
#define COEF_SQRT2                 COEF_CONST(1.4142135623731)
44
45
/* tables */
46
/* filters are mirrored in coef 6, second half left out */
47
static const real_t p8_13_20[7] =
48
{
49
    FRAC_CONST(0.00746082949812),
50
    FRAC_CONST(0.02270420949825),
51
    FRAC_CONST(0.04546865930473),
52
    FRAC_CONST(0.07266113929591),
53
    FRAC_CONST(0.09885108575264),
54
    FRAC_CONST(0.11793710567217),
55
    FRAC_CONST(0.125)
56
};
57
58
static const real_t p2_13_20[7] =
59
{
60
    FRAC_CONST(0.0),
61
    FRAC_CONST(0.01899487526049),
62
    FRAC_CONST(0.0),
63
    FRAC_CONST(-0.07293139167538),
64
    FRAC_CONST(0.0),
65
    FRAC_CONST(0.30596630545168),
66
    FRAC_CONST(0.5)
67
};
68
69
static const real_t p12_13_34[7] =
70
{
71
    FRAC_CONST(0.04081179924692),
72
    FRAC_CONST(0.03812810994926),
73
    FRAC_CONST(0.05144908135699),
74
    FRAC_CONST(0.06399831151592),
75
    FRAC_CONST(0.07428313801106),
76
    FRAC_CONST(0.08100347892914),
77
    FRAC_CONST(0.08333333333333)
78
};
79
80
static const real_t p8_13_34[7] =
81
{
82
    FRAC_CONST(0.01565675600122),
83
    FRAC_CONST(0.03752716391991),
84
    FRAC_CONST(0.05417891378782),
85
    FRAC_CONST(0.08417044116767),
86
    FRAC_CONST(0.10307344158036),
87
    FRAC_CONST(0.12222452249753),
88
    FRAC_CONST(0.125)
89
};
90
91
static const real_t p4_13_34[7] =
92
{
93
    FRAC_CONST(-0.05908211155639),
94
    FRAC_CONST(-0.04871498374946),
95
    FRAC_CONST(0.0),
96
    FRAC_CONST(0.07778723915851),
97
    FRAC_CONST(0.16486303567403),
98
    FRAC_CONST(0.23279856662996),
99
    FRAC_CONST(0.25)
100
};
101
102
#ifdef PARAM_32KHZ
103
static const uint8_t delay_length_d[2][NO_ALLPASS_LINKS] = {
104
    { 1, 2, 3 } /* d_24kHz */,
105
    { 3, 4, 5 } /* d_48kHz */
106
};
107
#else
108
static const uint8_t delay_length_d[NO_ALLPASS_LINKS] = {
109
    3, 4, 5 /* d_48kHz */
110
};
111
#endif
112
static const real_t filter_a[NO_ALLPASS_LINKS] = { /* a(m) = exp(-d_48kHz(m)/7) */
113
    FRAC_CONST(0.65143905753106),
114
    FRAC_CONST(0.56471812200776),
115
    FRAC_CONST(0.48954165955695)
116
};
117
118
static const uint8_t group_border20[10+12 + 1] =
119
{
120
    6, 7, 0, 1, 2, 3, /* 6 subqmf subbands */
121
    9, 8,             /* 2 subqmf subbands */
122
    10, 11,           /* 2 subqmf subbands */
123
    3, 4, 5, 6, 7, 8, 9, 11, 14, 18, 23, 35, 64
124
};
125
126
static const uint8_t group_border34[32+18 + 1] =
127
{
128
     0,  1,  2,  3,  4,  5,  6,  7,  8,  9,  10, 11, /* 12 subqmf subbands */
129
     12, 13, 14, 15, 16, 17, 18, 19,                 /*  8 subqmf subbands */
130
     20, 21, 22, 23,                                 /*  4 subqmf subbands */
131
     24, 25, 26, 27,                                 /*  4 subqmf subbands */
132
     28, 29, 30, 31,                                 /*  4 subqmf subbands */
133
     32-27, 33-27, 34-27, 35-27, 36-27, 37-27, 38-27, 40-27, 42-27, 44-27, 46-27, 48-27, 51-27, 54-27, 57-27, 60-27, 64-27, 68-27, 91-27
134
};
135
136
static const uint16_t map_group2bk20[10+12] =
137
{
138
    (NEGATE_IPD_MASK | 1), (NEGATE_IPD_MASK | 0),
139
    0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
140
};
141
142
static const uint16_t map_group2bk34[32+18] =
143
{
144
    0,  1,  2,  3,  4,  5,  6,  6,  7, (NEGATE_IPD_MASK | 2), (NEGATE_IPD_MASK | 1), (NEGATE_IPD_MASK | 0),
145
    10, 10, 4,  5,  6,  7,  8,  9,
146
    10, 11, 12, 9,
147
    14, 11, 12, 13,
148
    14, 15, 16, 13,
149
    16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33
150
};
151
152
/* type definitions */
153
typedef struct
154
{
155
    uint8_t frame_len;
156
    uint8_t resolution20[3];
157
    uint8_t resolution34[5];
158
159
    qmf_t *work;
160
    qmf_t **buffer;
161
    qmf_t **temp;
162
} hyb_info;
163
164
/* static function declarations */
165
static void ps_data_decode(ps_info *ps);
166
static hyb_info *hybrid_init(uint8_t numTimeSlotsRate);
167
static void channel_filter2(hyb_info *hyb, uint8_t frame_len, const real_t *filter,
168
                            qmf_t *buffer, qmf_t **X_hybrid);
169
static void INLINE DCT3_4_unscaled(real_t *y, real_t *x);
170
static void channel_filter8(hyb_info *hyb, uint8_t frame_len, const real_t *filter,
171
                            qmf_t *buffer, qmf_t **X_hybrid);
172
static void hybrid_analysis(hyb_info *hyb, qmf_t X[32][64], qmf_t X_hybrid[32][32],
173
                            uint8_t use34, uint8_t numTimeSlotsRate);
174
static void hybrid_synthesis(hyb_info *hyb, qmf_t X[32][64], qmf_t X_hybrid[32][32],
175
                             uint8_t use34, uint8_t numTimeSlotsRate);
176
static int8_t delta_clip(int8_t i, int8_t min, int8_t max);
177
static void delta_decode(uint8_t enable, int8_t *index, int8_t *index_prev,
178
                         uint8_t dt_flag, uint8_t nr_par, uint8_t stride,
179
                         int8_t min_index, int8_t max_index);
180
static void delta_modulo_decode(uint8_t enable, int8_t *index, int8_t *index_prev,
181
                                uint8_t dt_flag, uint8_t nr_par, uint8_t stride,
182
                                int8_t and_modulo);
183
static void map20indexto34(int8_t *index, uint8_t bins);
184
#ifdef PS_LOW_POWER
185
static void map34indexto20(int8_t *index, uint8_t bins);
186
#endif
187
static void ps_data_decode(ps_info *ps);
188
static void ps_decorrelate(ps_info *ps, qmf_t X_left[38][64], qmf_t X_right[38][64],
189
                           qmf_t X_hybrid_left[32][32], qmf_t X_hybrid_right[32][32]);
190
static void ps_mix_phase(ps_info *ps, qmf_t X_left[38][64], qmf_t X_right[38][64],
191
                         qmf_t X_hybrid_left[32][32], qmf_t X_hybrid_right[32][32]);
192
193
/*  */
194
195
196
static hyb_info *hybrid_init(uint8_t numTimeSlotsRate)
197
8.80k
{
198
8.80k
    uint8_t i;
199
200
8.80k
    hyb_info *hyb = (hyb_info*)faad_malloc(sizeof(hyb_info));
201
202
8.80k
    hyb->resolution34[0] = 12;
203
8.80k
    hyb->resolution34[1] = 8;
204
8.80k
    hyb->resolution34[2] = 4;
205
8.80k
    hyb->resolution34[3] = 4;
206
8.80k
    hyb->resolution34[4] = 4;
207
208
8.80k
    hyb->resolution20[0] = 8;
209
8.80k
    hyb->resolution20[1] = 2;
210
8.80k
    hyb->resolution20[2] = 2;
211
212
8.80k
    hyb->frame_len = numTimeSlotsRate;
213
214
8.80k
    hyb->work = (qmf_t*)faad_malloc((hyb->frame_len+12) * sizeof(qmf_t));
215
8.80k
    memset(hyb->work, 0, (hyb->frame_len+12) * sizeof(qmf_t));
216
217
8.80k
    hyb->buffer = (qmf_t**)faad_malloc(5 * sizeof(qmf_t*));
218
52.8k
    for (i = 0; i < 5; i++)
219
44.0k
    {
220
44.0k
        hyb->buffer[i] = (qmf_t*)faad_malloc(hyb->frame_len * sizeof(qmf_t));
221
44.0k
        memset(hyb->buffer[i], 0, hyb->frame_len * sizeof(qmf_t));
222
44.0k
    }
223
224
8.80k
    hyb->temp = (qmf_t**)faad_malloc(hyb->frame_len * sizeof(qmf_t*));
225
287k
    for (i = 0; i < hyb->frame_len; i++)
226
278k
    {
227
278k
        hyb->temp[i] = (qmf_t*)faad_malloc(12 /*max*/ * sizeof(qmf_t));
228
278k
    }
229
230
8.80k
    return hyb;
231
8.80k
}
232
233
static void hybrid_free(hyb_info *hyb)
234
8.80k
{
235
8.80k
    uint8_t i;
236
237
8.80k
  if (!hyb) return;
238
239
8.80k
    if (hyb->work)
240
8.80k
        faad_free(hyb->work);
241
242
52.8k
    for (i = 0; i < 5; i++)
243
44.0k
    {
244
44.0k
        if (hyb->buffer[i])
245
44.0k
            faad_free(hyb->buffer[i]);
246
44.0k
    }
247
8.80k
    if (hyb->buffer)
248
8.80k
        faad_free(hyb->buffer);
249
250
287k
    for (i = 0; i < hyb->frame_len; i++)
251
278k
    {
252
278k
        if (hyb->temp[i])
253
278k
            faad_free(hyb->temp[i]);
254
278k
    }
255
8.80k
    if (hyb->temp)
256
8.80k
        faad_free(hyb->temp);
257
258
8.80k
    faad_free(hyb);
259
8.80k
}
260
261
/* real filter, size 2 */
262
static void channel_filter2(hyb_info *hyb, uint8_t frame_len, const real_t *filter,
263
                            qmf_t *buffer, qmf_t **X_hybrid)
264
6.95k
{
265
6.95k
    uint8_t i;
266
6.95k
    (void)hyb;  /* TODO: remove parameter? */
267
268
224k
    for (i = 0; i < frame_len; i++)
269
217k
    {
270
217k
        real_t r0 = MUL_F(filter[0],(QMF_RE(buffer[0+i]) + QMF_RE(buffer[12+i])));
271
217k
        real_t r1 = MUL_F(filter[1],(QMF_RE(buffer[1+i]) + QMF_RE(buffer[11+i])));
272
217k
        real_t r2 = MUL_F(filter[2],(QMF_RE(buffer[2+i]) + QMF_RE(buffer[10+i])));
273
217k
        real_t r3 = MUL_F(filter[3],(QMF_RE(buffer[3+i]) + QMF_RE(buffer[9+i])));
274
217k
        real_t r4 = MUL_F(filter[4],(QMF_RE(buffer[4+i]) + QMF_RE(buffer[8+i])));
275
217k
        real_t r5 = MUL_F(filter[5],(QMF_RE(buffer[5+i]) + QMF_RE(buffer[7+i])));
276
217k
        real_t r6 = MUL_F(filter[6],QMF_RE(buffer[6+i]));
277
217k
        real_t i0 = MUL_F(filter[0],(QMF_IM(buffer[0+i]) + QMF_IM(buffer[12+i])));
278
217k
        real_t i1 = MUL_F(filter[1],(QMF_IM(buffer[1+i]) + QMF_IM(buffer[11+i])));
279
217k
        real_t i2 = MUL_F(filter[2],(QMF_IM(buffer[2+i]) + QMF_IM(buffer[10+i])));
280
217k
        real_t i3 = MUL_F(filter[3],(QMF_IM(buffer[3+i]) + QMF_IM(buffer[9+i])));
281
217k
        real_t i4 = MUL_F(filter[4],(QMF_IM(buffer[4+i]) + QMF_IM(buffer[8+i])));
282
217k
        real_t i5 = MUL_F(filter[5],(QMF_IM(buffer[5+i]) + QMF_IM(buffer[7+i])));
283
217k
        real_t i6 = MUL_F(filter[6],QMF_IM(buffer[6+i]));
284
285
        /* q = 0 */
286
217k
        QMF_RE(X_hybrid[i][0]) = r0 + r1 + r2 + r3 + r4 + r5 + r6;
287
217k
        QMF_IM(X_hybrid[i][0]) = i0 + i1 + i2 + i3 + i4 + i5 + i6;
288
289
        /* q = 1 */
290
217k
        QMF_RE(X_hybrid[i][1]) = r0 - r1 + r2 - r3 + r4 - r5 + r6;
291
217k
        QMF_IM(X_hybrid[i][1]) = i0 - i1 + i2 - i3 + i4 - i5 + i6;
292
217k
    }
293
6.95k
}
294
295
/* complex filter, size 4 */
296
static void channel_filter4(hyb_info *hyb, uint8_t frame_len, const real_t *filter,
297
                            qmf_t *buffer, qmf_t **X_hybrid)
298
6.11k
{
299
6.11k
    uint8_t i;
300
6.11k
    real_t input_re1[2], input_re2[2], input_im1[2], input_im2[2];
301
6.11k
    (void)hyb;  /* TODO: remove parameter? */
302
303
193k
    for (i = 0; i < frame_len; i++)
304
187k
    {
305
187k
        input_re1[0] = -MUL_F(filter[2], (QMF_RE(buffer[i+2]) + QMF_RE(buffer[i+10]))) +
306
187k
            MUL_F(filter[6], QMF_RE(buffer[i+6]));
307
187k
        input_re1[1] = MUL_F(FRAC_CONST(-0.70710678118655),
308
187k
            (MUL_F(filter[1], (QMF_RE(buffer[i+1]) + QMF_RE(buffer[i+11]))) +
309
187k
            MUL_F(filter[3], (QMF_RE(buffer[i+3]) + QMF_RE(buffer[i+9]))) -
310
187k
            MUL_F(filter[5], (QMF_RE(buffer[i+5]) + QMF_RE(buffer[i+7])))));
311
312
187k
        input_im1[0] = MUL_F(filter[0], (QMF_IM(buffer[i+0]) - QMF_IM(buffer[i+12]))) -
313
187k
            MUL_F(filter[4], (QMF_IM(buffer[i+4]) - QMF_IM(buffer[i+8])));
314
187k
        input_im1[1] = MUL_F(FRAC_CONST(0.70710678118655),
315
187k
            (MUL_F(filter[1], (QMF_IM(buffer[i+1]) - QMF_IM(buffer[i+11]))) -
316
187k
            MUL_F(filter[3], (QMF_IM(buffer[i+3]) - QMF_IM(buffer[i+9]))) -
317
187k
            MUL_F(filter[5], (QMF_IM(buffer[i+5]) - QMF_IM(buffer[i+7])))));
318
319
187k
        input_re2[0] = MUL_F(filter[0], (QMF_RE(buffer[i+0]) - QMF_RE(buffer[i+12]))) -
320
187k
            MUL_F(filter[4], (QMF_RE(buffer[i+4]) - QMF_RE(buffer[i+8])));
321
187k
        input_re2[1] = MUL_F(FRAC_CONST(0.70710678118655),
322
187k
            (MUL_F(filter[1], (QMF_RE(buffer[i+1]) - QMF_RE(buffer[i+11]))) -
323
187k
            MUL_F(filter[3], (QMF_RE(buffer[i+3]) - QMF_RE(buffer[i+9]))) -
324
187k
            MUL_F(filter[5], (QMF_RE(buffer[i+5]) - QMF_RE(buffer[i+7])))));
325
326
187k
        input_im2[0] = -MUL_F(filter[2], (QMF_IM(buffer[i+2]) + QMF_IM(buffer[i+10]))) +
327
187k
            MUL_F(filter[6], QMF_IM(buffer[i+6]));
328
187k
        input_im2[1] = MUL_F(FRAC_CONST(-0.70710678118655),
329
187k
            (MUL_F(filter[1], (QMF_IM(buffer[i+1]) + QMF_IM(buffer[i+11]))) +
330
187k
            MUL_F(filter[3], (QMF_IM(buffer[i+3]) + QMF_IM(buffer[i+9]))) -
331
187k
            MUL_F(filter[5], (QMF_IM(buffer[i+5]) + QMF_IM(buffer[i+7])))));
332
333
        /* q == 0 */
334
187k
        QMF_RE(X_hybrid[i][0]) =  input_re1[0] + input_re1[1] + input_im1[0] + input_im1[1];
335
187k
        QMF_IM(X_hybrid[i][0]) = -input_re2[0] - input_re2[1] + input_im2[0] + input_im2[1];
336
337
        /* q == 1 */
338
187k
        QMF_RE(X_hybrid[i][1]) =  input_re1[0] - input_re1[1] - input_im1[0] + input_im1[1];
339
187k
        QMF_IM(X_hybrid[i][1]) =  input_re2[0] - input_re2[1] + input_im2[0] - input_im2[1];
340
341
        /* q == 2 */
342
187k
        QMF_RE(X_hybrid[i][2]) =  input_re1[0] - input_re1[1] + input_im1[0] - input_im1[1];
343
187k
        QMF_IM(X_hybrid[i][2]) = -input_re2[0] + input_re2[1] + input_im2[0] - input_im2[1];
344
345
        /* q == 3 */
346
187k
        QMF_RE(X_hybrid[i][3]) =  input_re1[0] + input_re1[1] - input_im1[0] - input_im1[1];
347
187k
        QMF_IM(X_hybrid[i][3]) =  input_re2[0] + input_re2[1] + input_im2[0] + input_im2[1];
348
187k
    }
349
6.11k
}
350
351
static void INLINE DCT3_4_unscaled(real_t *y, real_t *x)
352
684k
{
353
684k
    real_t f0, f1, f2, f3, f4, f5, f6, f7, f8;
354
355
684k
    f0 = MUL_F(x[2], FRAC_CONST(0.7071067811865476));
356
684k
    f1 = x[0] - f0;
357
684k
    f2 = x[0] + f0;
358
684k
    f3 = x[1] + x[3];
359
684k
    f4 = MUL_C(x[1], COEF_CONST(1.3065629648763766));
360
684k
    f5 = MUL_F(f3, FRAC_CONST(-0.9238795325112866));
361
684k
    f6 = MUL_F(x[3], FRAC_CONST(-0.5411961001461967));
362
684k
    f7 = f4 + f5;
363
684k
    f8 = f6 - f5;
364
684k
    y[3] = f2 - f8;
365
684k
    y[0] = f2 + f8;
366
684k
    y[2] = f1 - f7;
367
684k
    y[1] = f1 + f7;
368
684k
}
369
370
/* complex filter, size 8 */
371
static void channel_filter8(hyb_info *hyb, uint8_t frame_len, const real_t *filter,
372
                            qmf_t *buffer, qmf_t **X_hybrid)
373
5.51k
{
374
5.51k
    uint8_t i, n;
375
5.51k
    real_t input_re1[4], input_re2[4], input_im1[4], input_im2[4];
376
5.51k
    real_t x[4];
377
5.51k
    (void)hyb;  /* TODO: remove parameter? */
378
379
176k
    for (i = 0; i < frame_len; i++)
380
171k
    {
381
171k
        input_re1[0] =  MUL_F(filter[6],QMF_RE(buffer[6+i]));
382
171k
        input_re1[1] =  MUL_F(filter[5],(QMF_RE(buffer[5+i]) + QMF_RE(buffer[7+i])));
383
171k
        input_re1[2] = -MUL_F(filter[0],(QMF_RE(buffer[0+i]) + QMF_RE(buffer[12+i]))) + MUL_F(filter[4],(QMF_RE(buffer[4+i]) + QMF_RE(buffer[8+i])));
384
171k
        input_re1[3] = -MUL_F(filter[1],(QMF_RE(buffer[1+i]) + QMF_RE(buffer[11+i]))) + MUL_F(filter[3],(QMF_RE(buffer[3+i]) + QMF_RE(buffer[9+i])));
385
386
171k
        input_im1[0] = MUL_F(filter[5],(QMF_IM(buffer[7+i]) - QMF_IM(buffer[5+i])));
387
171k
        input_im1[1] = MUL_F(filter[0],(QMF_IM(buffer[12+i]) - QMF_IM(buffer[0+i]))) + MUL_F(filter[4],(QMF_IM(buffer[8+i]) - QMF_IM(buffer[4+i])));
388
171k
        input_im1[2] = MUL_F(filter[1],(QMF_IM(buffer[11+i]) - QMF_IM(buffer[1+i]))) + MUL_F(filter[3],(QMF_IM(buffer[9+i]) - QMF_IM(buffer[3+i])));
389
171k
        input_im1[3] = MUL_F(filter[2],(QMF_IM(buffer[10+i]) - QMF_IM(buffer[2+i])));
390
391
855k
        for (n = 0; n < 4; n++)
392
684k
        {
393
684k
            x[n] = input_re1[n] - input_im1[3-n];
394
684k
        }
395
171k
        DCT3_4_unscaled(x, x);
396
171k
        QMF_RE(X_hybrid[i][7]) = x[0];
397
171k
        QMF_RE(X_hybrid[i][5]) = x[2];
398
171k
        QMF_RE(X_hybrid[i][3]) = x[3];
399
171k
        QMF_RE(X_hybrid[i][1]) = x[1];
400
401
855k
        for (n = 0; n < 4; n++)
402
684k
        {
403
684k
            x[n] = input_re1[n] + input_im1[3-n];
404
684k
        }
405
171k
        DCT3_4_unscaled(x, x);
406
171k
        QMF_RE(X_hybrid[i][6]) = x[1];
407
171k
        QMF_RE(X_hybrid[i][4]) = x[3];
408
171k
        QMF_RE(X_hybrid[i][2]) = x[2];
409
171k
        QMF_RE(X_hybrid[i][0]) = x[0];
410
411
171k
        input_im2[0] =  MUL_F(filter[6],QMF_IM(buffer[6+i]));
412
171k
        input_im2[1] =  MUL_F(filter[5],(QMF_IM(buffer[5+i]) + QMF_IM(buffer[7+i])));
413
171k
        input_im2[2] = -MUL_F(filter[0],(QMF_IM(buffer[0+i]) + QMF_IM(buffer[12+i]))) + MUL_F(filter[4],(QMF_IM(buffer[4+i]) + QMF_IM(buffer[8+i])));
414
171k
        input_im2[3] = -MUL_F(filter[1],(QMF_IM(buffer[1+i]) + QMF_IM(buffer[11+i]))) + MUL_F(filter[3],(QMF_IM(buffer[3+i]) + QMF_IM(buffer[9+i])));
415
416
171k
        input_re2[0] = MUL_F(filter[5],(QMF_RE(buffer[7+i]) - QMF_RE(buffer[5+i])));
417
171k
        input_re2[1] = MUL_F(filter[0],(QMF_RE(buffer[12+i]) - QMF_RE(buffer[0+i]))) + MUL_F(filter[4],(QMF_RE(buffer[8+i]) - QMF_RE(buffer[4+i])));
418
171k
        input_re2[2] = MUL_F(filter[1],(QMF_RE(buffer[11+i]) - QMF_RE(buffer[1+i]))) + MUL_F(filter[3],(QMF_RE(buffer[9+i]) - QMF_RE(buffer[3+i])));
419
171k
        input_re2[3] = MUL_F(filter[2],(QMF_RE(buffer[10+i]) - QMF_RE(buffer[2+i])));
420
421
855k
        for (n = 0; n < 4; n++)
422
684k
        {
423
684k
            x[n] = input_im2[n] + input_re2[3-n];
424
684k
        }
425
171k
        DCT3_4_unscaled(x, x);
426
171k
        QMF_IM(X_hybrid[i][7]) = x[0];
427
171k
        QMF_IM(X_hybrid[i][5]) = x[2];
428
171k
        QMF_IM(X_hybrid[i][3]) = x[3];
429
171k
        QMF_IM(X_hybrid[i][1]) = x[1];
430
431
855k
        for (n = 0; n < 4; n++)
432
684k
        {
433
684k
            x[n] = input_im2[n] - input_re2[3-n];
434
684k
        }
435
171k
        DCT3_4_unscaled(x, x);
436
171k
        QMF_IM(X_hybrid[i][6]) = x[1];
437
171k
        QMF_IM(X_hybrid[i][4]) = x[3];
438
171k
        QMF_IM(X_hybrid[i][2]) = x[2];
439
171k
        QMF_IM(X_hybrid[i][0]) = x[0];
440
171k
    }
441
5.51k
}
442
443
static void INLINE DCT3_6_unscaled(real_t *y, real_t *x)
444
250k
{
445
250k
    real_t f0, f1, f2, f3, f4, f5, f6, f7;
446
447
250k
    f0 = MUL_F(x[3], FRAC_CONST(0.70710678118655));
448
250k
    f1 = x[0] + f0;
449
250k
    f2 = x[0] - f0;
450
250k
    f3 = MUL_F((x[1] - x[5]), FRAC_CONST(0.70710678118655));
451
250k
    f4 = MUL_F(x[2], FRAC_CONST(0.86602540378444)) + MUL_F(x[4], FRAC_CONST(0.5));
452
250k
    f5 = f4 - x[4];
453
250k
    f6 = MUL_F(x[1], FRAC_CONST(0.96592582628907)) + MUL_F(x[5], FRAC_CONST(0.25881904510252));
454
250k
    f7 = f6 - f3;
455
250k
    y[0] = f1 + f6 + f4;
456
250k
    y[1] = f2 + f3 - x[4];
457
250k
    y[2] = f7 + f2 - f5;
458
250k
    y[3] = f1 - f7 - f5;
459
250k
    y[4] = f1 - f3 - x[4];
460
250k
    y[5] = f2 - f6 + f4;
461
250k
}
462
463
/* complex filter, size 12 */
464
static void channel_filter12(hyb_info *hyb, uint8_t frame_len, const real_t *filter,
465
                             qmf_t *buffer, qmf_t **X_hybrid)
466
2.03k
{
467
2.03k
    uint8_t i, n;
468
2.03k
    real_t input_re1[6], input_re2[6], input_im1[6], input_im2[6];
469
2.03k
    real_t out_re1[6], out_re2[6], out_im1[6], out_im2[6];
470
2.03k
    (void)hyb;  /* TODO: remove parameter? */
471
472
64.5k
    for (i = 0; i < frame_len; i++)
473
62.5k
    {
474
437k
        for (n = 0; n < 6; n++)
475
375k
        {
476
375k
            if (n == 0)
477
62.5k
            {
478
62.5k
                input_re1[0] = MUL_F(QMF_RE(buffer[6+i]), filter[6]);
479
62.5k
                input_re2[0] = MUL_F(QMF_IM(buffer[6+i]), filter[6]);
480
312k
            } else {
481
312k
                input_re1[6-n] = MUL_F((QMF_RE(buffer[n+i]) + QMF_RE(buffer[12-n+i])), filter[n]);
482
312k
                input_re2[6-n] = MUL_F((QMF_IM(buffer[n+i]) + QMF_IM(buffer[12-n+i])), filter[n]);
483
312k
            }
484
375k
            input_im2[n] = MUL_F((QMF_RE(buffer[n+i]) - QMF_RE(buffer[12-n+i])), filter[n]);
485
375k
            input_im1[n] = MUL_F((QMF_IM(buffer[n+i]) - QMF_IM(buffer[12-n+i])), filter[n]);
486
375k
        }
487
488
62.5k
        DCT3_6_unscaled(out_re1, input_re1);
489
62.5k
        DCT3_6_unscaled(out_re2, input_re2);
490
491
62.5k
        DCT3_6_unscaled(out_im1, input_im1);
492
62.5k
        DCT3_6_unscaled(out_im2, input_im2);
493
494
250k
        for (n = 0; n < 6; n += 2)
495
187k
        {
496
187k
            QMF_RE(X_hybrid[i][n]) = out_re1[n] - out_im1[n];
497
187k
            QMF_IM(X_hybrid[i][n]) = out_re2[n] + out_im2[n];
498
187k
            QMF_RE(X_hybrid[i][n+1]) = out_re1[n+1] + out_im1[n+1];
499
187k
            QMF_IM(X_hybrid[i][n+1]) = out_re2[n+1] - out_im2[n+1];
500
501
187k
            QMF_RE(X_hybrid[i][10-n]) = out_re1[n+1] - out_im1[n+1];
502
187k
            QMF_IM(X_hybrid[i][10-n]) = out_re2[n+1] + out_im2[n+1];
503
187k
            QMF_RE(X_hybrid[i][11-n]) = out_re1[n] + out_im1[n];
504
187k
            QMF_IM(X_hybrid[i][11-n]) = out_re2[n] - out_im2[n];
505
187k
        }
506
62.5k
    }
507
2.03k
}
508
509
/* Hybrid analysis: further split up QMF subbands
510
 * to improve frequency resolution
511
 */
512
static void hybrid_analysis(hyb_info *hyb, qmf_t X[32][64], qmf_t X_hybrid[32][32],
513
                            uint8_t use34, uint8_t numTimeSlotsRate)
514
5.51k
{
515
5.51k
    uint8_t k, n, band;
516
5.51k
    uint8_t offset = 0;
517
5.51k
    uint8_t qmf_bands = (use34) ? 5 : 3;
518
5.51k
    uint8_t *resolution = (use34) ? hyb->resolution34 : hyb->resolution20;
519
520
26.1k
    for (band = 0; band < qmf_bands; band++)
521
20.6k
    {
522
        /* build working buffer */
523
20.6k
        memcpy(hyb->work, hyb->buffer[band], 12 * sizeof(qmf_t));
524
525
        /* add new samples */
526
659k
        for (n = 0; n < hyb->frame_len; n++)
527
638k
        {
528
638k
            QMF_RE(hyb->work[12 + n]) = QMF_RE(X[n + 6 /*delay*/][band]);
529
638k
            QMF_IM(hyb->work[12 + n]) = QMF_IM(X[n + 6 /*delay*/][band]);
530
638k
        }
531
532
        /* store samples */
533
20.6k
        memcpy(hyb->buffer[band], hyb->work + hyb->frame_len, 12 * sizeof(qmf_t));
534
535
536
20.6k
        switch(resolution[band])
537
20.6k
        {
538
6.95k
        case 2:
539
            /* Type B real filter, Q[p] = 2 */
540
6.95k
            channel_filter2(hyb, hyb->frame_len, p2_13_20, hyb->work, hyb->temp);
541
6.95k
            break;
542
6.11k
        case 4:
543
            /* Type A complex filter, Q[p] = 4 */
544
6.11k
            channel_filter4(hyb, hyb->frame_len, p4_13_34, hyb->work, hyb->temp);
545
6.11k
            break;
546
5.51k
        case 8:
547
            /* Type A complex filter, Q[p] = 8 */
548
5.51k
            channel_filter8(hyb, hyb->frame_len, (use34) ? p8_13_34 : p8_13_20,
549
5.51k
                hyb->work, hyb->temp);
550
5.51k
            break;
551
2.03k
        case 12:
552
            /* Type A complex filter, Q[p] = 12 */
553
2.03k
            channel_filter12(hyb, hyb->frame_len, p12_13_34, hyb->work, hyb->temp);
554
2.03k
            break;
555
20.6k
        }
556
557
659k
        for (n = 0; n < hyb->frame_len; n++)
558
638k
        {
559
3.94M
            for (k = 0; k < resolution[band]; k++)
560
3.30M
            {
561
3.30M
                QMF_RE(X_hybrid[n][offset + k]) = QMF_RE(hyb->temp[n][k]);
562
3.30M
                QMF_IM(X_hybrid[n][offset + k]) = QMF_IM(hyb->temp[n][k]);
563
3.30M
            }
564
638k
        }
565
20.6k
        offset += resolution[band];
566
20.6k
    }
567
568
    /* group hybrid channels */
569
5.51k
    if (!use34)
570
3.47k
    {
571
112k
        for (n = 0; n < numTimeSlotsRate; n++)
572
108k
        {
573
108k
            QMF_RE(X_hybrid[n][3]) += QMF_RE(X_hybrid[n][4]);
574
108k
            QMF_IM(X_hybrid[n][3]) += QMF_IM(X_hybrid[n][4]);
575
108k
            QMF_RE(X_hybrid[n][4]) = 0;
576
108k
            QMF_IM(X_hybrid[n][4]) = 0;
577
578
108k
            QMF_RE(X_hybrid[n][2]) += QMF_RE(X_hybrid[n][5]);
579
108k
            QMF_IM(X_hybrid[n][2]) += QMF_IM(X_hybrid[n][5]);
580
108k
            QMF_RE(X_hybrid[n][5]) = 0;
581
108k
            QMF_IM(X_hybrid[n][5]) = 0;
582
108k
        }
583
3.47k
    }
584
5.51k
}
585
586
static void hybrid_synthesis(hyb_info *hyb, qmf_t X[32][64], qmf_t X_hybrid[32][32],
587
                             uint8_t use34, uint8_t numTimeSlotsRate)
588
11.0k
{
589
11.0k
    uint8_t k, n, band;
590
11.0k
    uint8_t offset = 0;
591
11.0k
    uint8_t qmf_bands = (use34) ? 5 : 3;
592
11.0k
    uint8_t *resolution = (use34) ? hyb->resolution34 : hyb->resolution20;
593
11.0k
    (void)numTimeSlotsRate;  /* TODO: remove parameter? */
594
595
52.2k
    for(band = 0; band < qmf_bands; band++)
596
41.2k
    {
597
1.31M
        for (n = 0; n < hyb->frame_len; n++)
598
1.27M
        {
599
1.27M
            QMF_RE(X[n][band]) = 0;
600
1.27M
            QMF_IM(X[n][band]) = 0;
601
602
7.88M
            for (k = 0; k < resolution[band]; k++)
603
6.60M
            {
604
6.60M
                QMF_RE(X[n][band]) += QMF_RE(X_hybrid[n][offset + k]);
605
6.60M
                QMF_IM(X[n][band]) += QMF_IM(X_hybrid[n][offset + k]);
606
6.60M
            }
607
1.27M
        }
608
41.2k
        offset += resolution[band];
609
41.2k
    }
610
11.0k
}
611
612
/* limits the value i to the range [min,max] */
613
static int8_t delta_clip(int8_t i, int8_t min, int8_t max)
614
69.0k
{
615
69.0k
    if (i < min)
616
8.79k
        return min;
617
60.2k
    else if (i > max)
618
1.23k
        return max;
619
59.0k
    else
620
59.0k
        return i;
621
69.0k
}
622
623
//int iid = 0;
624
625
/* delta decode array */
626
static void delta_decode(uint8_t enable, int8_t *index, int8_t *index_prev,
627
                         uint8_t dt_flag, uint8_t nr_par, uint8_t stride,
628
                         int8_t min_index, int8_t max_index)
629
15.4k
{
630
15.4k
    int8_t i;
631
632
15.4k
    if (enable == 1)
633
7.30k
    {
634
7.30k
        if (dt_flag == 0)
635
4.79k
        {
636
            /* delta coded in frequency direction */
637
4.79k
            index[0] = 0 + index[0];
638
4.79k
            index[0] = delta_clip(index[0], min_index, max_index);
639
640
50.7k
            for (i = 1; i < nr_par; i++)
641
45.9k
            {
642
45.9k
                index[i] = index[i-1] + index[i];
643
45.9k
                index[i] = delta_clip(index[i], min_index, max_index);
644
45.9k
            }
645
4.79k
        } else {
646
            /* delta coded in time direction */
647
20.8k
            for (i = 0; i < nr_par; i++)
648
18.3k
            {
649
                //int8_t tmp2;
650
                //int8_t tmp = index[i];
651
652
                //printf("%d %d\n", index_prev[i*stride], index[i]);
653
                //printf("%d\n", index[i]);
654
655
18.3k
                index[i] = index_prev[i*stride] + index[i];
656
                //tmp2 = index[i];
657
18.3k
                index[i] = delta_clip(index[i], min_index, max_index);
658
659
                //if (iid)
660
                //{
661
                //    if (index[i] == 7)
662
                //    {
663
                //        printf("%d %d %d\n", index_prev[i*stride], tmp, tmp2);
664
                //    }
665
                //}
666
18.3k
            }
667
2.51k
        }
668
8.15k
    } else {
669
        /* set indices to zero */
670
16.2k
        for (i = 0; i < nr_par; i++)
671
8.10k
        {
672
8.10k
            index[i] = 0;
673
8.10k
        }
674
8.15k
    }
675
676
    /* coarse */
677
15.4k
    if (stride == 2)
678
8.86k
    {
679
29.7k
        for (i = (nr_par<<1)-1; i > 0; i--)
680
20.8k
        {
681
20.8k
            index[i] = index[i>>1];
682
20.8k
        }
683
8.86k
    }
684
15.4k
}
685
686
/* delta modulo decode array */
687
/* in: log2 value of the modulo value to allow using AND instead of MOD */
688
static void delta_modulo_decode(uint8_t enable, int8_t *index, int8_t *index_prev,
689
                                uint8_t dt_flag, uint8_t nr_par, uint8_t stride,
690
                                int8_t and_modulo)
691
15.4k
{
692
15.4k
    int8_t i;
693
694
15.4k
    if (enable == 1)
695
3.85k
    {
696
3.85k
        if (dt_flag == 0)
697
1.88k
        {
698
            /* delta coded in frequency direction */
699
1.88k
            index[0] = 0 + index[0];
700
1.88k
            index[0] &= and_modulo;
701
702
12.0k
            for (i = 1; i < nr_par; i++)
703
10.1k
            {
704
10.1k
                index[i] = index[i-1] + index[i];
705
10.1k
                index[i] &= and_modulo;
706
10.1k
            }
707
1.97k
        } else {
708
            /* delta coded in time direction */
709
5.38k
            for (i = 0; i < nr_par; i++)
710
3.41k
            {
711
3.41k
                index[i] = index_prev[i*stride] + index[i];
712
3.41k
                index[i] &= and_modulo;
713
3.41k
            }
714
1.97k
        }
715
11.6k
    } else {
716
        /* set indices to zero */
717
46.5k
        for (i = 0; i < nr_par; i++)
718
34.9k
        {
719
34.9k
            index[i] = 0;
720
34.9k
        }
721
11.6k
    }
722
723
    /* coarse */
724
15.4k
    if (stride == 2)
725
0
    {
726
0
        index[0] = 0;
727
0
        for (i = (nr_par<<1)-1; i > 0; i--)
728
0
        {
729
0
            index[i] = index[i>>1];
730
0
        }
731
0
    }
732
15.4k
}
733
734
#ifdef PS_LOW_POWER
735
static void map34indexto20(int8_t *index, uint8_t bins)
736
{
737
    index[0] = (2*index[0]+index[1])/3;
738
    index[1] = (index[1]+2*index[2])/3;
739
    index[2] = (2*index[3]+index[4])/3;
740
    index[3] = (index[4]+2*index[5])/3;
741
    index[4] = (index[6]+index[7])/2;
742
    index[5] = (index[8]+index[9])/2;
743
    index[6] = index[10];
744
    index[7] = index[11];
745
    index[8] = (index[12]+index[13])/2;
746
    index[9] = (index[14]+index[15])/2;
747
    index[10] = index[16];
748
749
    if (bins == 34)
750
    {
751
        index[11] = index[17];
752
        index[12] = index[18];
753
        index[13] = index[19];
754
        index[14] = (index[20]+index[21])/2;
755
        index[15] = (index[22]+index[23])/2;
756
        index[16] = (index[24]+index[25])/2;
757
        index[17] = (index[26]+index[27])/2;
758
        index[18] = (index[28]+index[29]+index[30]+index[31])/4;
759
        index[19] = (index[32]+index[33])/2;
760
    }
761
}
762
#endif
763
764
static void map20indexto34(int8_t *index, uint8_t bins)
765
5.86k
{
766
5.86k
    index[0] = index[0];
767
5.86k
    index[1] = (index[0] + index[1])/2;
768
5.86k
    index[2] = index[1];
769
5.86k
    index[3] = index[2];
770
5.86k
    index[4] = (index[2] + index[3])/2;
771
5.86k
    index[5] = index[3];
772
5.86k
    index[6] = index[4];
773
5.86k
    index[7] = index[4];
774
5.86k
    index[8] = index[5];
775
5.86k
    index[9] = index[5];
776
5.86k
    index[10] = index[6];
777
5.86k
    index[11] = index[7];
778
5.86k
    index[12] = index[8];
779
5.86k
    index[13] = index[8];
780
5.86k
    index[14] = index[9];
781
5.86k
    index[15] = index[9];
782
5.86k
    index[16] = index[10];
783
784
5.86k
    if (bins == 34)
785
3.15k
    {
786
3.15k
        index[17] = index[11];
787
3.15k
        index[18] = index[12];
788
3.15k
        index[19] = index[13];
789
3.15k
        index[20] = index[14];
790
3.15k
        index[21] = index[14];
791
3.15k
        index[22] = index[15];
792
3.15k
        index[23] = index[15];
793
3.15k
        index[24] = index[16];
794
3.15k
        index[25] = index[16];
795
3.15k
        index[26] = index[17];
796
3.15k
        index[27] = index[17];
797
3.15k
        index[28] = index[18];
798
3.15k
        index[29] = index[18];
799
3.15k
        index[30] = index[18];
800
3.15k
        index[31] = index[18];
801
3.15k
        index[32] = index[19];
802
3.15k
        index[33] = index[19];
803
3.15k
    }
804
5.86k
}
805
806
/* parse the bitstream data decoded in ps_data() */
807
static void ps_data_decode(ps_info *ps)
808
5.51k
{
809
5.51k
    uint8_t env, bin;
810
811
    /* ps data not available, use data from previous frame */
812
5.51k
    if (ps->ps_data_available == 0)
813
1.57k
    {
814
1.57k
        ps->num_env = 0;
815
1.57k
    }
816
817
13.2k
    for (env = 0; env < ps->num_env; env++)
818
7.73k
    {
819
7.73k
        int8_t *iid_index_prev;
820
7.73k
        int8_t *icc_index_prev;
821
7.73k
        int8_t *ipd_index_prev;
822
7.73k
        int8_t *opd_index_prev;
823
824
7.73k
        int8_t num_iid_steps = (ps->iid_mode < 3) ? 7 : 15 /*fine quant*/;
825
826
7.73k
        if (env == 0)
827
2.47k
        {
828
            /* take last envelope from previous frame */
829
2.47k
            iid_index_prev = ps->iid_index_prev;
830
2.47k
            icc_index_prev = ps->icc_index_prev;
831
2.47k
            ipd_index_prev = ps->ipd_index_prev;
832
2.47k
            opd_index_prev = ps->opd_index_prev;
833
5.26k
        } else {
834
            /* take index values from previous envelope */
835
5.26k
            iid_index_prev = ps->iid_index[env - 1];
836
5.26k
            icc_index_prev = ps->icc_index[env - 1];
837
5.26k
            ipd_index_prev = ps->ipd_index[env - 1];
838
5.26k
            opd_index_prev = ps->opd_index[env - 1];
839
5.26k
        }
840
841
//        iid = 1;
842
        /* delta decode iid parameters */
843
7.73k
        delta_decode(ps->enable_iid, ps->iid_index[env], iid_index_prev,
844
7.73k
            ps->iid_dt[env], ps->nr_iid_par,
845
7.73k
            (ps->iid_mode == 0 || ps->iid_mode == 3) ? 2 : 1,
846
7.73k
            -num_iid_steps, num_iid_steps);
847
//        iid = 0;
848
849
        /* delta decode icc parameters */
850
7.73k
        delta_decode(ps->enable_icc, ps->icc_index[env], icc_index_prev,
851
7.73k
            ps->icc_dt[env], ps->nr_icc_par,
852
7.73k
            (ps->icc_mode == 0 || ps->icc_mode == 3) ? 2 : 1,
853
7.73k
            0, 7);
854
855
        /* delta modulo decode ipd parameters */
856
7.73k
        delta_modulo_decode(ps->enable_ipdopd, ps->ipd_index[env], ipd_index_prev,
857
7.73k
            ps->ipd_dt[env], ps->nr_ipdopd_par, 1, 7);
858
859
        /* delta modulo decode opd parameters */
860
7.73k
        delta_modulo_decode(ps->enable_ipdopd, ps->opd_index[env], opd_index_prev,
861
7.73k
            ps->opd_dt[env], ps->nr_ipdopd_par, 1, 7);
862
7.73k
    }
863
864
    /* handle error case */
865
5.51k
    if (ps->num_env == 0)
866
3.04k
    {
867
        /* force to 1 */
868
3.04k
        ps->num_env = 1;
869
870
3.04k
        if (ps->enable_iid)
871
2.12k
        {
872
74.3k
            for (bin = 0; bin < 34; bin++)
873
72.1k
                ps->iid_index[0][bin] = ps->iid_index_prev[bin];
874
2.12k
        } else {
875
32.1k
            for (bin = 0; bin < 34; bin++)
876
31.2k
                ps->iid_index[0][bin] = 0;
877
919
        }
878
879
3.04k
        if (ps->enable_icc)
880
1.53k
        {
881
53.6k
            for (bin = 0; bin < 34; bin++)
882
52.1k
                ps->icc_index[0][bin] = ps->icc_index_prev[bin];
883
1.53k
        } else {
884
52.7k
            for (bin = 0; bin < 34; bin++)
885
51.2k
                ps->icc_index[0][bin] = 0;
886
1.50k
        }
887
888
3.04k
        if (ps->enable_ipdopd)
889
537
        {
890
9.66k
            for (bin = 0; bin < 17; bin++)
891
9.12k
            {
892
9.12k
                ps->ipd_index[0][bin] = ps->ipd_index_prev[bin];
893
9.12k
                ps->opd_index[0][bin] = ps->opd_index_prev[bin];
894
9.12k
            }
895
2.50k
        } else {
896
45.0k
            for (bin = 0; bin < 17; bin++)
897
42.5k
            {
898
42.5k
                ps->ipd_index[0][bin] = 0;
899
42.5k
                ps->opd_index[0][bin] = 0;
900
42.5k
            }
901
2.50k
        }
902
3.04k
    }
903
904
    /* update previous indices */
905
193k
    for (bin = 0; bin < 34; bin++)
906
187k
        ps->iid_index_prev[bin] = ps->iid_index[ps->num_env-1][bin];
907
193k
    for (bin = 0; bin < 34; bin++)
908
187k
        ps->icc_index_prev[bin] = ps->icc_index[ps->num_env-1][bin];
909
99.2k
    for (bin = 0; bin < 17; bin++)
910
93.7k
    {
911
93.7k
        ps->ipd_index_prev[bin] = ps->ipd_index[ps->num_env-1][bin];
912
93.7k
        ps->opd_index_prev[bin] = ps->opd_index[ps->num_env-1][bin];
913
93.7k
    }
914
915
5.51k
    ps->ps_data_available = 0;
916
917
5.51k
    if (ps->frame_class == 0)
918
3.72k
    {
919
3.72k
        ps->border_position[0] = 0;
920
6.15k
        for (env = 1; env < ps->num_env; env++)
921
2.43k
        {
922
2.43k
            ps->border_position[env] = (env * ps->numTimeSlotsRate) / ps->num_env;
923
2.43k
        }
924
3.72k
        ps->border_position[ps->num_env] = ps->numTimeSlotsRate;
925
3.72k
    } else {
926
1.79k
        ps->border_position[0] = 0;
927
928
1.79k
        if (ps->border_position[ps->num_env] < ps->numTimeSlotsRate)
929
1.49k
        {
930
52.3k
            for (bin = 0; bin < 34; bin++)
931
50.8k
            {
932
50.8k
                ps->iid_index[ps->num_env][bin] = ps->iid_index[ps->num_env-1][bin];
933
50.8k
                ps->icc_index[ps->num_env][bin] = ps->icc_index[ps->num_env-1][bin];
934
50.8k
            }
935
26.9k
            for (bin = 0; bin < 17; bin++)
936
25.4k
            {
937
25.4k
                ps->ipd_index[ps->num_env][bin] = ps->ipd_index[ps->num_env-1][bin];
938
25.4k
                ps->opd_index[ps->num_env][bin] = ps->opd_index[ps->num_env-1][bin];
939
25.4k
            }
940
1.49k
            ps->num_env++;
941
1.49k
            ps->border_position[ps->num_env] = ps->numTimeSlotsRate;
942
1.49k
        }
943
944
6.11k
        for (env = 1; env < ps->num_env; env++)
945
4.31k
        {
946
4.31k
            int8_t thr = ps->numTimeSlotsRate - (ps->num_env - env);
947
948
4.31k
            if (ps->border_position[env] > thr)
949
658
            {
950
658
                ps->border_position[env] = thr;
951
3.66k
            } else {
952
3.66k
                thr = ps->border_position[env-1]+1;
953
3.66k
                if (ps->border_position[env] < thr)
954
1.57k
                {
955
1.57k
                    ps->border_position[env] = thr;
956
1.57k
                }
957
3.66k
            }
958
4.31k
        }
959
1.79k
    }
960
961
    /* make sure that the indices of all parameters can be mapped
962
     * to the same hybrid synthesis filterbank
963
     */
964
#ifdef PS_LOW_POWER
965
    for (env = 0; env < ps->num_env; env++)
966
    {
967
        if (ps->iid_mode == 2 || ps->iid_mode == 5)
968
            map34indexto20(ps->iid_index[env], 34);
969
        if (ps->icc_mode == 2 || ps->icc_mode == 5)
970
            map34indexto20(ps->icc_index[env], 34);
971
972
        /* disable ipd/opd */
973
        for (bin = 0; bin < 17; bin++)
974
        {
975
            ps->aaIpdIndex[env][bin] = 0;
976
            ps->aaOpdIndex[env][bin] = 0;
977
        }
978
    }
979
#else
980
5.51k
    if (ps->use34hybrid_bands)
981
2.03k
    {
982
5.32k
        for (env = 0; env < ps->num_env; env++)
983
3.28k
        {
984
3.28k
            if (ps->iid_mode != 2 && ps->iid_mode != 5)
985
1.35k
                map20indexto34(ps->iid_index[env], 34);
986
3.28k
            if (ps->icc_mode != 2 && ps->icc_mode != 5)
987
1.79k
                map20indexto34(ps->icc_index[env], 34);
988
3.28k
            if (ps->ipd_mode != 2 && ps->ipd_mode != 5)
989
1.35k
            {
990
1.35k
                map20indexto34(ps->ipd_index[env], 17);
991
1.35k
                map20indexto34(ps->opd_index[env], 17);
992
1.35k
            }
993
3.28k
        }
994
2.03k
    }
995
5.51k
#endif
996
997
#if 0
998
    for (env = 0; env < ps->num_env; env++)
999
    {
1000
        printf("iid[env:%d]:", env);
1001
        for (bin = 0; bin < 34; bin++)
1002
        {
1003
            printf(" %d", ps->iid_index[env][bin]);
1004
        }
1005
        printf("\n");
1006
    }
1007
    for (env = 0; env < ps->num_env; env++)
1008
    {
1009
        printf("icc[env:%d]:", env);
1010
        for (bin = 0; bin < 34; bin++)
1011
        {
1012
            printf(" %d", ps->icc_index[env][bin]);
1013
        }
1014
        printf("\n");
1015
    }
1016
    for (env = 0; env < ps->num_env; env++)
1017
    {
1018
        printf("ipd[env:%d]:", env);
1019
        for (bin = 0; bin < 17; bin++)
1020
        {
1021
            printf(" %d", ps->ipd_index[env][bin]);
1022
        }
1023
        printf("\n");
1024
    }
1025
    for (env = 0; env < ps->num_env; env++)
1026
    {
1027
        printf("opd[env:%d]:", env);
1028
        for (bin = 0; bin < 17; bin++)
1029
        {
1030
            printf(" %d", ps->opd_index[env][bin]);
1031
        }
1032
        printf("\n");
1033
    }
1034
    printf("\n");
1035
#endif
1036
5.51k
}
1037
1038
/* decorrelate the mono signal using an allpass filter */
1039
static void ps_decorrelate(ps_info *ps, qmf_t X_left[38][64], qmf_t X_right[38][64],
1040
                           qmf_t X_hybrid_left[32][32], qmf_t X_hybrid_right[32][32])
1041
5.51k
{
1042
5.51k
    uint8_t gr, n, bk;
1043
5.51k
    uint8_t temp_delay = 0;
1044
5.51k
    uint8_t sb, maxsb;
1045
5.51k
    const complex_t *Phi_Fract_SubQmf;
1046
5.51k
    uint8_t temp_delay_ser[NO_ALLPASS_LINKS];
1047
5.51k
    real_t P_SmoothPeakDecayDiffNrg, nrg;
1048
5.51k
    real_t P[32][34];
1049
5.51k
    real_t G_TransientRatio[32][34] = {{0}};
1050
5.51k
    complex_t inputLeft;
1051
1052
1053
    /* chose hybrid filterbank: 20 or 34 band case */
1054
5.51k
    if (ps->use34hybrid_bands)
1055
2.03k
    {
1056
2.03k
        Phi_Fract_SubQmf = Phi_Fract_SubQmf34;
1057
3.47k
    } else{
1058
3.47k
        Phi_Fract_SubQmf = Phi_Fract_SubQmf20;
1059
3.47k
    }
1060
1061
    /* clear the energy values */
1062
181k
    for (n = 0; n < 32; n++)
1063
176k
    {
1064
6.17M
        for (bk = 0; bk < 34; bk++)
1065
6.00M
        {
1066
6.00M
            P[n][bk] = 0;
1067
6.00M
        }
1068
176k
    }
1069
1070
    /* calculate the energy in each parameter band b(k) */
1071
183k
    for (gr = 0; gr < ps->num_groups; gr++)
1072
178k
    {
1073
        /* select the parameter index b(k) to which this group belongs */
1074
178k
        bk = (~NEGATE_IPD_MASK) & ps->map_group2bk[gr];
1075
1076
        /* select the upper subband border for this group */
1077
178k
        maxsb = (gr < ps->num_hybrid_groups) ? ps->group_border[gr]+1 : ps->group_border[gr+1];
1078
1079
610k
        for (sb = ps->group_border[gr]; sb < maxsb; sb++)
1080
432k
        {
1081
13.8M
            for (n = ps->border_position[0]; n < ps->border_position[ps->num_env]; n++)
1082
13.4M
            {
1083
#ifdef FIXED_POINT
1084
                uint32_t in_re, in_im;
1085
#endif
1086
1087
                /* input from hybrid subbands or QMF subbands */
1088
13.4M
                if (gr < ps->num_hybrid_groups)
1089
3.09M
                {
1090
3.09M
                    RE(inputLeft) = QMF_RE(X_hybrid_left[n][sb]);
1091
3.09M
                    IM(inputLeft) = QMF_IM(X_hybrid_left[n][sb]);
1092
10.3M
                } else {
1093
10.3M
                    RE(inputLeft) = QMF_RE(X_left[n][sb]);
1094
10.3M
                    IM(inputLeft) = QMF_IM(X_left[n][sb]);
1095
10.3M
                }
1096
1097
                /* accumulate energy */
1098
#ifdef FIXED_POINT
1099
                /* NOTE: all input is scaled by 2^(-5) because of fixed point QMF
1100
                 * meaning that P will be scaled by 2^(-10) compared to floating point version
1101
                 */
1102
                in_re = ((abs(RE(inputLeft))+(1<<(REAL_BITS-1)))>>REAL_BITS);
1103
                in_im = ((abs(IM(inputLeft))+(1<<(REAL_BITS-1)))>>REAL_BITS);
1104
                P[n][bk] += in_re*in_re + in_im*in_im;
1105
#else
1106
13.4M
                P[n][bk] += MUL_R(RE(inputLeft),RE(inputLeft)) + MUL_R(IM(inputLeft),IM(inputLeft));
1107
13.4M
#endif
1108
13.4M
            }
1109
432k
        }
1110
178k
    }
1111
1112
#if 0
1113
    for (n = 0; n < 32; n++)
1114
    {
1115
        for (bk = 0; bk < 34; bk++)
1116
        {
1117
#ifdef FIXED_POINT
1118
            printf("%d %d: %d\n", n, bk, P[n][bk] /*/(float)REAL_PRECISION*/);
1119
#else
1120
            printf("%d %d: %f\n", n, bk, P[n][bk]/1024.0);
1121
#endif
1122
        }
1123
    }
1124
#endif
1125
1126
    /* calculate transient reduction ratio for each parameter band b(k) */
1127
144k
    for (bk = 0; bk < ps->nr_par_bands; bk++)
1128
138k
    {
1129
4.44M
        for (n = ps->border_position[0]; n < ps->border_position[ps->num_env]; n++)
1130
4.30M
        {
1131
4.30M
            const real_t gamma = COEF_CONST(1.5);
1132
1133
4.30M
            ps->P_PeakDecayNrg[bk] = MUL_F(ps->P_PeakDecayNrg[bk], ps->alpha_decay);
1134
4.30M
            if (ps->P_PeakDecayNrg[bk] < P[n][bk])
1135
46.6k
                ps->P_PeakDecayNrg[bk] = P[n][bk];
1136
1137
            /* apply smoothing filter to peak decay energy */
1138
4.30M
            P_SmoothPeakDecayDiffNrg = ps->P_SmoothPeakDecayDiffNrg_prev[bk];
1139
4.30M
            P_SmoothPeakDecayDiffNrg += MUL_F((ps->P_PeakDecayNrg[bk] - P[n][bk] - ps->P_SmoothPeakDecayDiffNrg_prev[bk]), ps->alpha_smooth);
1140
4.30M
            ps->P_SmoothPeakDecayDiffNrg_prev[bk] = P_SmoothPeakDecayDiffNrg;
1141
1142
            /* apply smoothing filter to energy */
1143
4.30M
            nrg = ps->P_prev[bk];
1144
4.30M
            nrg += MUL_F((P[n][bk] - ps->P_prev[bk]), ps->alpha_smooth);
1145
4.30M
            ps->P_prev[bk] = nrg;
1146
1147
            /* calculate transient ratio */
1148
4.30M
            if (MUL_C(P_SmoothPeakDecayDiffNrg, gamma) <= nrg)
1149
4.27M
            {
1150
4.27M
                G_TransientRatio[n][bk] = REAL_CONST(1.0);
1151
4.27M
            } else {
1152
36.2k
                G_TransientRatio[n][bk] = DIV_R(nrg, (MUL_C(P_SmoothPeakDecayDiffNrg, gamma)));
1153
36.2k
            }
1154
4.30M
        }
1155
138k
    }
1156
1157
#if 0
1158
    for (n = 0; n < 32; n++)
1159
    {
1160
        for (bk = 0; bk < 34; bk++)
1161
        {
1162
#ifdef FIXED_POINT
1163
            printf("%d %d: %f\n", n, bk, G_TransientRatio[n][bk]/(float)REAL_PRECISION);
1164
#else
1165
            printf("%d %d: %f\n", n, bk, G_TransientRatio[n][bk]);
1166
#endif
1167
        }
1168
    }
1169
#endif
1170
1171
    /* apply stereo decorrelation filter to the signal */
1172
183k
    for (gr = 0; gr < ps->num_groups; gr++)
1173
178k
    {
1174
178k
        if (gr < ps->num_hybrid_groups)
1175
99.9k
            maxsb = ps->group_border[gr] + 1;
1176
78.4k
        else
1177
78.4k
            maxsb = ps->group_border[gr + 1];
1178
1179
        /* QMF channel */
1180
610k
        for (sb = ps->group_border[gr]; sb < maxsb; sb++)
1181
432k
        {
1182
432k
            real_t g_DecaySlope;
1183
432k
            real_t g_DecaySlope_filt[NO_ALLPASS_LINKS];
1184
1185
            /* g_DecaySlope: [0..1] */
1186
432k
            if (gr < ps->num_hybrid_groups || sb <= ps->decay_cutoff)
1187
105k
            {
1188
105k
                g_DecaySlope = FRAC_CONST(1.0);
1189
326k
            } else {
1190
326k
                int8_t decay = ps->decay_cutoff - sb;
1191
326k
                if (decay <= -20 /* -1/DECAY_SLOPE */)
1192
222k
                {
1193
222k
                    g_DecaySlope = 0;
1194
222k
                } else {
1195
                    /* decay(int)*decay_slope(frac) = g_DecaySlope(frac) */
1196
104k
                    g_DecaySlope = FRAC_CONST(1.0) + DECAY_SLOPE * decay;
1197
104k
                }
1198
326k
            }
1199
1200
            /* calculate g_DecaySlope_filt for every n multiplied by filter_a[n] */
1201
1.72M
            for (n = 0; n < NO_ALLPASS_LINKS; n++)
1202
1.29M
            {
1203
1.29M
                g_DecaySlope_filt[n] = MUL_F(g_DecaySlope, filter_a[n]);
1204
1.29M
            }
1205
1206
1207
            /* set delay indices */
1208
432k
            temp_delay = ps->saved_delay;
1209
1.72M
            for (n = 0; n < NO_ALLPASS_LINKS; n++)
1210
1.29M
                temp_delay_ser[n] = ps->delay_buf_index_ser[n];
1211
1212
13.8M
            for (n = ps->border_position[0]; n < ps->border_position[ps->num_env]; n++)
1213
13.4M
            {
1214
13.4M
                complex_t tmp, tmp0, R0;
1215
13.4M
                uint8_t m;
1216
1217
13.4M
                if (gr < ps->num_hybrid_groups)
1218
3.09M
                {
1219
                    /* hybrid filterbank input */
1220
3.09M
                    RE(inputLeft) = QMF_RE(X_hybrid_left[n][sb]);
1221
3.09M
                    IM(inputLeft) = QMF_IM(X_hybrid_left[n][sb]);
1222
10.3M
                } else {
1223
                    /* QMF filterbank input */
1224
10.3M
                    RE(inputLeft) = QMF_RE(X_left[n][sb]);
1225
10.3M
                    IM(inputLeft) = QMF_IM(X_left[n][sb]);
1226
10.3M
                }
1227
1228
13.4M
                if (sb > ps->nr_allpass_bands && gr >= ps->num_hybrid_groups)
1229
7.03M
                {
1230
                    /* delay */
1231
1232
                    /* never hybrid subbands here, always QMF subbands */
1233
7.03M
                    RE(tmp) = RE(ps->delay_Qmf[ps->delay_buf_index_delay[sb]][sb]);
1234
7.03M
                    IM(tmp) = IM(ps->delay_Qmf[ps->delay_buf_index_delay[sb]][sb]);
1235
7.03M
                    RE(R0) = RE(tmp);
1236
7.03M
                    IM(R0) = IM(tmp);
1237
7.03M
                    RE(ps->delay_Qmf[ps->delay_buf_index_delay[sb]][sb]) = RE(inputLeft);
1238
7.03M
                    IM(ps->delay_Qmf[ps->delay_buf_index_delay[sb]][sb]) = IM(inputLeft);
1239
7.03M
                } else {
1240
                    /* allpass filter */
1241
6.39M
                    complex_t Phi_Fract;
1242
1243
                    /* fetch parameters */
1244
6.39M
                    if (gr < ps->num_hybrid_groups)
1245
3.09M
                    {
1246
                        /* select data from the hybrid subbands */
1247
3.09M
                        RE(tmp0) = RE(ps->delay_SubQmf[temp_delay][sb]);
1248
3.09M
                        IM(tmp0) = IM(ps->delay_SubQmf[temp_delay][sb]);
1249
1250
3.09M
                        RE(ps->delay_SubQmf[temp_delay][sb]) = RE(inputLeft);
1251
3.09M
                        IM(ps->delay_SubQmf[temp_delay][sb]) = IM(inputLeft);
1252
1253
3.09M
                        RE(Phi_Fract) = RE(Phi_Fract_SubQmf[sb]);
1254
3.09M
                        IM(Phi_Fract) = IM(Phi_Fract_SubQmf[sb]);
1255
3.30M
                    } else {
1256
                        /* select data from the QMF subbands */
1257
3.30M
                        RE(tmp0) = RE(ps->delay_Qmf[temp_delay][sb]);
1258
3.30M
                        IM(tmp0) = IM(ps->delay_Qmf[temp_delay][sb]);
1259
1260
3.30M
                        RE(ps->delay_Qmf[temp_delay][sb]) = RE(inputLeft);
1261
3.30M
                        IM(ps->delay_Qmf[temp_delay][sb]) = IM(inputLeft);
1262
1263
3.30M
                        RE(Phi_Fract) = RE(Phi_Fract_Qmf[sb]);
1264
3.30M
                        IM(Phi_Fract) = IM(Phi_Fract_Qmf[sb]);
1265
3.30M
                    }
1266
1267
                    /* z^(-2) * Phi_Fract[k] */
1268
6.39M
                    ComplexMult(&RE(tmp), &IM(tmp), RE(tmp0), IM(tmp0), RE(Phi_Fract), IM(Phi_Fract));
1269
1270
6.39M
                    RE(R0) = RE(tmp);
1271
6.39M
                    IM(R0) = IM(tmp);
1272
25.5M
                    for (m = 0; m < NO_ALLPASS_LINKS; m++)
1273
19.1M
                    {
1274
19.1M
                        complex_t Q_Fract_allpass, tmp2;
1275
1276
                        /* fetch parameters */
1277
19.1M
                        if (gr < ps->num_hybrid_groups)
1278
9.27M
                        {
1279
                            /* select data from the hybrid subbands */
1280
9.27M
                            RE(tmp0) = RE(ps->delay_SubQmf_ser[m][temp_delay_ser[m]][sb]);
1281
9.27M
                            IM(tmp0) = IM(ps->delay_SubQmf_ser[m][temp_delay_ser[m]][sb]);
1282
1283
9.27M
                            if (ps->use34hybrid_bands)
1284
6.00M
                            {
1285
6.00M
                                RE(Q_Fract_allpass) = RE(Q_Fract_allpass_SubQmf34[sb][m]);
1286
6.00M
                                IM(Q_Fract_allpass) = IM(Q_Fract_allpass_SubQmf34[sb][m]);
1287
6.00M
                            } else {
1288
3.26M
                                RE(Q_Fract_allpass) = RE(Q_Fract_allpass_SubQmf20[sb][m]);
1289
3.26M
                                IM(Q_Fract_allpass) = IM(Q_Fract_allpass_SubQmf20[sb][m]);
1290
3.26M
                            }
1291
9.91M
                        } else {
1292
                            /* select data from the QMF subbands */
1293
9.91M
                            RE(tmp0) = RE(ps->delay_Qmf_ser[m][temp_delay_ser[m]][sb]);
1294
9.91M
                            IM(tmp0) = IM(ps->delay_Qmf_ser[m][temp_delay_ser[m]][sb]);
1295
1296
9.91M
                            RE(Q_Fract_allpass) = RE(Q_Fract_allpass_Qmf[sb][m]);
1297
9.91M
                            IM(Q_Fract_allpass) = IM(Q_Fract_allpass_Qmf[sb][m]);
1298
9.91M
                        }
1299
1300
                        /* delay by a fraction */
1301
                        /* z^(-d(m)) * Q_Fract_allpass[k,m] */
1302
19.1M
                        ComplexMult(&RE(tmp), &IM(tmp), RE(tmp0), IM(tmp0), RE(Q_Fract_allpass), IM(Q_Fract_allpass));
1303
1304
                        /* -a(m) * g_DecaySlope[k] */
1305
19.1M
                        RE(tmp) += -MUL_F(g_DecaySlope_filt[m], RE(R0));
1306
19.1M
                        IM(tmp) += -MUL_F(g_DecaySlope_filt[m], IM(R0));
1307
1308
                        /* -a(m) * g_DecaySlope[k] * Q_Fract_allpass[k,m] * z^(-d(m)) */
1309
19.1M
                        RE(tmp2) = RE(R0) + MUL_F(g_DecaySlope_filt[m], RE(tmp));
1310
19.1M
                        IM(tmp2) = IM(R0) + MUL_F(g_DecaySlope_filt[m], IM(tmp));
1311
1312
                        /* store sample */
1313
19.1M
                        if (gr < ps->num_hybrid_groups)
1314
9.27M
                        {
1315
9.27M
                            RE(ps->delay_SubQmf_ser[m][temp_delay_ser[m]][sb]) = RE(tmp2);
1316
9.27M
                            IM(ps->delay_SubQmf_ser[m][temp_delay_ser[m]][sb]) = IM(tmp2);
1317
9.91M
                        } else {
1318
9.91M
                            RE(ps->delay_Qmf_ser[m][temp_delay_ser[m]][sb]) = RE(tmp2);
1319
9.91M
                            IM(ps->delay_Qmf_ser[m][temp_delay_ser[m]][sb]) = IM(tmp2);
1320
9.91M
                        }
1321
1322
                        /* store for next iteration (or as output value if last iteration) */
1323
19.1M
                        RE(R0) = RE(tmp);
1324
19.1M
                        IM(R0) = IM(tmp);
1325
19.1M
                    }
1326
6.39M
                }
1327
1328
                /* select b(k) for reading the transient ratio */
1329
13.4M
                bk = (~NEGATE_IPD_MASK) & ps->map_group2bk[gr];
1330
1331
                /* duck if a past transient is found */
1332
13.4M
                RE(R0) = MUL_R(G_TransientRatio[n][bk], RE(R0));
1333
13.4M
                IM(R0) = MUL_R(G_TransientRatio[n][bk], IM(R0));
1334
1335
13.4M
                if (gr < ps->num_hybrid_groups)
1336
3.09M
                {
1337
                    /* hybrid */
1338
3.09M
                    QMF_RE(X_hybrid_right[n][sb]) = RE(R0);
1339
3.09M
                    QMF_IM(X_hybrid_right[n][sb]) = IM(R0);
1340
10.3M
                } else {
1341
                    /* QMF */
1342
10.3M
                    QMF_RE(X_right[n][sb]) = RE(R0);
1343
10.3M
                    QMF_IM(X_right[n][sb]) = IM(R0);
1344
10.3M
                }
1345
1346
                /* Update delay buffer index */
1347
13.4M
                if (++temp_delay >= 2)
1348
6.71M
                {
1349
6.71M
                    temp_delay = 0;
1350
6.71M
                }
1351
1352
                /* update delay indices */
1353
13.4M
                if (sb > ps->nr_allpass_bands && gr >= ps->num_hybrid_groups)
1354
7.03M
                {
1355
                    /* delay_D depends on the samplerate, it can hold the values 14 and 1 */
1356
7.03M
                    if (++ps->delay_buf_index_delay[sb] >= ps->delay_D[sb])
1357
5.10M
                    {
1358
5.10M
                        ps->delay_buf_index_delay[sb] = 0;
1359
5.10M
                    }
1360
7.03M
                }
1361
1362
53.7M
                for (m = 0; m < NO_ALLPASS_LINKS; m++)
1363
40.2M
                {
1364
40.2M
                    if (++temp_delay_ser[m] >= ps->num_sample_delay_ser[m])
1365
10.3M
                    {
1366
10.3M
                        temp_delay_ser[m] = 0;
1367
10.3M
                    }
1368
40.2M
                }
1369
13.4M
            }
1370
432k
        }
1371
178k
    }
1372
1373
    /* update delay indices */
1374
5.51k
    ps->saved_delay = temp_delay;
1375
22.0k
    for (n = 0; n < NO_ALLPASS_LINKS; n++)
1376
16.5k
        ps->delay_buf_index_ser[n] = temp_delay_ser[n];
1377
5.51k
}
1378
1379
#if 0
1380
#ifdef FIXED_POINT
1381
#define step(shift) \
1382
    if ((0x40000000l >> shift) + root <= value)       \
1383
    {                                                 \
1384
        value -= (0x40000000l >> shift) + root;       \
1385
        root = (root >> 1) | (0x40000000l >> shift);  \
1386
    } else {                                          \
1387
        root = root >> 1;                             \
1388
    }
1389
1390
/* fixed point square root approximation */
1391
static real_t ps_sqrt(real_t value)
1392
{
1393
    real_t root = 0;
1394
1395
    step( 0); step( 2); step( 4); step( 6);
1396
    step( 8); step(10); step(12); step(14);
1397
    step(16); step(18); step(20); step(22);
1398
    step(24); step(26); step(28); step(30);
1399
1400
    if (root < value)
1401
        ++root;
1402
1403
    root <<= (REAL_BITS/2);
1404
1405
    return root;
1406
}
1407
#else
1408
#define ps_sqrt(A) sqrt(A)
1409
#endif
1410
#endif
1411
1412
static const real_t ipdopd_cos_tab[] = {
1413
    FRAC_CONST(1.000000000000000),
1414
    FRAC_CONST(0.707106781186548),
1415
    FRAC_CONST(0.000000000000000),
1416
    FRAC_CONST(-0.707106781186547),
1417
    FRAC_CONST(-1.000000000000000),
1418
    FRAC_CONST(-0.707106781186548),
1419
    FRAC_CONST(-0.000000000000000),
1420
    FRAC_CONST(0.707106781186547),
1421
    FRAC_CONST(1.000000000000000)
1422
};
1423
1424
static const real_t ipdopd_sin_tab[] = {
1425
    FRAC_CONST(0.000000000000000),
1426
    FRAC_CONST(0.707106781186547),
1427
    FRAC_CONST(1.000000000000000),
1428
    FRAC_CONST(0.707106781186548),
1429
    FRAC_CONST(0.000000000000000),
1430
    FRAC_CONST(-0.707106781186547),
1431
    FRAC_CONST(-1.000000000000000),
1432
    FRAC_CONST(-0.707106781186548),
1433
    FRAC_CONST(-0.000000000000000)
1434
};
1435
1436
static real_t magnitude_c(complex_t c)
1437
92.7k
{
1438
#ifdef FIXED_POINT
1439
#define ps_abs(A) (((A) > 0) ? (A) : (-(A)))
1440
#define ALPHA FRAC_CONST(0.948059448969)
1441
#define BETA  FRAC_CONST(0.392699081699)
1442
1443
    real_t abs_inphase = ps_abs(RE(c));
1444
    real_t abs_quadrature = ps_abs(IM(c));
1445
1446
    if (abs_inphase > abs_quadrature) {
1447
        return MUL_F(abs_inphase, ALPHA) + MUL_F(abs_quadrature, BETA);
1448
    } else {
1449
        return MUL_F(abs_quadrature, ALPHA) + MUL_F(abs_inphase, BETA);
1450
    }
1451
#else
1452
92.7k
    return sqrt(RE(c)*RE(c) + IM(c)*IM(c));
1453
92.7k
#endif
1454
92.7k
}
1455
1456
static void ps_mix_phase(ps_info *ps, qmf_t X_left[38][64], qmf_t X_right[38][64],
1457
                         qmf_t X_hybrid_left[32][32], qmf_t X_hybrid_right[32][32])
1458
5.51k
{
1459
5.51k
    uint8_t n;
1460
5.51k
    uint8_t gr;
1461
5.51k
    uint8_t bk = 0;
1462
5.51k
    uint8_t sb, maxsb;
1463
5.51k
    uint8_t env;
1464
5.51k
    uint8_t nr_ipdopd_par;
1465
5.51k
    complex_t h11, h12, h21, h22;  // COEF
1466
5.51k
    complex_t H11, H12, H21, H22;  // COEF
1467
5.51k
    complex_t deltaH11, deltaH12, deltaH21, deltaH22;  // COEF
1468
5.51k
    complex_t tempLeft, tempRight; // FRAC
1469
5.51k
    complex_t phaseLeft, phaseRight; // FRAC
1470
5.51k
    real_t L;
1471
5.51k
    const real_t *sf_iid;
1472
5.51k
    uint8_t no_iid_steps;
1473
1474
5.51k
    if (ps->iid_mode >= 3)
1475
2.45k
    {
1476
2.45k
        no_iid_steps = 15;
1477
2.45k
        sf_iid = sf_iid_fine;
1478
3.05k
    } else {
1479
3.05k
        no_iid_steps = 7;
1480
3.05k
        sf_iid = sf_iid_normal;
1481
3.05k
    }
1482
1483
5.51k
    if (ps->ipd_mode == 0 || ps->ipd_mode == 3)
1484
2.61k
    {
1485
2.61k
        nr_ipdopd_par = 11; /* resolution */
1486
2.89k
    } else {
1487
2.89k
        nr_ipdopd_par = ps->nr_ipdopd_par;
1488
2.89k
    }
1489
1490
183k
    for (gr = 0; gr < ps->num_groups; gr++)
1491
178k
    {
1492
178k
        bk = (~NEGATE_IPD_MASK) & ps->map_group2bk[gr];
1493
1494
        /* use one channel per group in the subqmf domain */
1495
178k
        maxsb = (gr < ps->num_hybrid_groups) ? ps->group_border[gr] + 1 : ps->group_border[gr + 1];
1496
1497
540k
        for (env = 0; env < ps->num_env; env++)
1498
361k
        {
1499
361k
            uint8_t abs_iid = (uint8_t)abs(ps->iid_index[env][bk]);
1500
            /* index range is supposed to be -7...7 or -15...15 depending on iid_mode
1501
                (Table 8.24, ISO/IEC 14496-3:2005).
1502
                if it is outside these boundaries, this is most likely an error. sanitize
1503
                it and try to process further. */
1504
361k
            if (ps->iid_index[env][bk] < -no_iid_steps) {
1505
122
                fprintf(stderr, "Warning: invalid iid_index: %d < %d\n", ps->iid_index[env][bk],
1506
122
                    -no_iid_steps);
1507
122
                ps->iid_index[env][bk] = -no_iid_steps;
1508
122
                abs_iid = no_iid_steps;
1509
361k
            } else if (ps->iid_index[env][bk] > no_iid_steps) {
1510
78
                fprintf(stderr, "Warning: invalid iid_index: %d > %d\n", ps->iid_index[env][bk],
1511
78
                    no_iid_steps);
1512
78
                ps->iid_index[env][bk] = no_iid_steps;
1513
78
                abs_iid = no_iid_steps;
1514
78
            }
1515
361k
            if (ps->icc_index[env][bk] < 0) {
1516
262
                fprintf(stderr, "Warning: invalid icc_index: %d < 0\n", ps->icc_index[env][bk]);
1517
262
                ps->icc_index[env][bk] = 0;
1518
361k
            } else if (ps->icc_index[env][bk] > 7) {
1519
0
                fprintf(stderr, "Warning: invalid icc_index: %d > 7\n", ps->icc_index[env][bk]);
1520
0
                ps->icc_index[env][bk] = 7;
1521
0
            }
1522
1523
361k
            if (ps->icc_mode < 3)
1524
219k
            {
1525
                /* type 'A' mixing as described in 8.6.4.6.2.1 */
1526
219k
                real_t c_1, c_2;  // COEF
1527
219k
                real_t cosa, sina;  // COEF
1528
219k
                real_t cosb, sinb;  // COEF
1529
219k
                real_t ab1, ab2;  // COEF
1530
219k
                real_t ab3, ab4;  // COEF
1531
1532
                /*
1533
                c_1 = sqrt(2.0 / (1.0 + pow(10.0, quant_iid[no_iid_steps + iid_index] / 10.0)));
1534
                c_2 = sqrt(2.0 / (1.0 + pow(10.0, quant_iid[no_iid_steps - iid_index] / 10.0)));
1535
                alpha = 0.5 * acos(quant_rho[icc_index]);
1536
                beta = alpha * ( c_1 - c_2 ) / sqrt(2.0);
1537
                */
1538
1539
                //printf("%d\n", ps->iid_index[env][bk]);
1540
1541
                /* calculate the scalefactors c_1 and c_2 from the intensity differences */
1542
219k
                c_1 = sf_iid[no_iid_steps + ps->iid_index[env][bk]];
1543
219k
                c_2 = sf_iid[no_iid_steps - ps->iid_index[env][bk]];
1544
1545
                /* calculate alpha and beta using the ICC parameters */
1546
219k
                cosa = cos_alphas[ps->icc_index[env][bk]];
1547
219k
                sina = sin_alphas[ps->icc_index[env][bk]];
1548
1549
219k
                if (ps->iid_mode >= 3)
1550
81.2k
                {
1551
81.2k
                    cosb = cos_betas_fine[abs_iid][ps->icc_index[env][bk]];
1552
81.2k
                    sinb = sin_betas_fine[abs_iid][ps->icc_index[env][bk]];
1553
138k
                } else {
1554
138k
                    cosb = cos_betas_normal[abs_iid][ps->icc_index[env][bk]];
1555
138k
                    sinb = sin_betas_normal[abs_iid][ps->icc_index[env][bk]];
1556
138k
                }
1557
1558
219k
                ab1 = MUL_C(cosb, cosa);
1559
219k
                ab2 = MUL_C(sinb, sina);
1560
219k
                ab3 = MUL_C(sinb, cosa);
1561
219k
                ab4 = MUL_C(cosb, sina);
1562
1563
                /* h_xy: COEF */
1564
219k
                RE(h11) = MUL_C(c_2, (ab1 - ab2));
1565
219k
                RE(h12) = MUL_C(c_1, (ab1 + ab2));
1566
219k
                RE(h21) = MUL_C(c_2, (ab3 + ab4));
1567
219k
                RE(h22) = MUL_C(c_1, (ab3 - ab4));
1568
219k
            } else {
1569
                /* type 'B' mixing as described in 8.6.4.6.2.2 */
1570
142k
                real_t sina, cosa;  // COEF
1571
142k
                real_t cosg, sing;  // COEF
1572
1573
                /*
1574
                real_t c, rho, mu, alpha, gamma;
1575
                uint8_t i;
1576
1577
                i = ps->iid_index[env][bk];
1578
                c = (real_t)pow(10.0, ((i)?(((i>0)?1:-1)*quant_iid[((i>0)?i:-i)-1]):0.)/20.0);
1579
                rho = quant_rho[ps->icc_index[env][bk]];
1580
1581
                if (rho == 0.0f && c == 1.)
1582
                {
1583
                    alpha = (real_t)M_PI/4.0f;
1584
                    rho = 0.05f;
1585
                } else {
1586
                    if (rho <= 0.05f)
1587
                    {
1588
                        rho = 0.05f;
1589
                    }
1590
                    alpha = 0.5f*(real_t)atan( (2.0f*c*rho) / (c*c-1.0f) );
1591
1592
                    if (alpha < 0.)
1593
                    {
1594
                        alpha += (real_t)M_PI/2.0f;
1595
                    }
1596
                    if (rho < 0.)
1597
                    {
1598
                        alpha += (real_t)M_PI;
1599
                    }
1600
                }
1601
                mu = c+1.0f/c;
1602
                mu = 1+(4.0f*rho*rho-4.0f)/(mu*mu);
1603
                gamma = (real_t)atan(sqrt((1.0f-sqrt(mu))/(1.0f+sqrt(mu))));
1604
                */
1605
1606
142k
                if (ps->iid_mode >= 3)
1607
70.4k
                {
1608
70.4k
                    cosa = sincos_alphas_B_fine[no_iid_steps + ps->iid_index[env][bk]][ps->icc_index[env][bk]];
1609
70.4k
                    sina = sincos_alphas_B_fine[30 - (no_iid_steps + ps->iid_index[env][bk])][ps->icc_index[env][bk]];
1610
70.4k
                    cosg = cos_gammas_fine[abs_iid][ps->icc_index[env][bk]];
1611
70.4k
                    sing = sin_gammas_fine[abs_iid][ps->icc_index[env][bk]];
1612
71.9k
                } else {
1613
71.9k
                    cosa = sincos_alphas_B_normal[no_iid_steps + ps->iid_index[env][bk]][ps->icc_index[env][bk]];
1614
71.9k
                    sina = sincos_alphas_B_normal[14 - (no_iid_steps + ps->iid_index[env][bk])][ps->icc_index[env][bk]];
1615
71.9k
                    cosg = cos_gammas_normal[abs_iid][ps->icc_index[env][bk]];
1616
71.9k
                    sing = sin_gammas_normal[abs_iid][ps->icc_index[env][bk]];
1617
71.9k
                }
1618
1619
142k
                RE(h11) = MUL_C(COEF_SQRT2, MUL_C(cosa, cosg));
1620
142k
                RE(h12) = MUL_C(COEF_SQRT2, MUL_C(sina, cosg));
1621
142k
                RE(h21) = MUL_C(COEF_SQRT2, MUL_C(-cosa, sing));
1622
142k
                RE(h22) = MUL_C(COEF_SQRT2, MUL_C(sina, sing));
1623
142k
            }
1624
361k
            IM(h11) = IM(h12) = IM(h21) = IM(h22) = 0;
1625
1626
            /* calculate phase rotation parameters H_xy */
1627
            /* note that the imaginary part of these parameters are only calculated when
1628
               IPD and OPD are enabled
1629
             */
1630
361k
            if ((ps->enable_ipdopd) && (bk < nr_ipdopd_par))
1631
46.3k
            {
1632
46.3k
                int8_t i;
1633
46.3k
                real_t xy, pq, xypq;  // FRAC
1634
1635
                /* ringbuffer index */
1636
46.3k
                i = ps->phase_hist;
1637
1638
                /* previous value */
1639
#ifdef FIXED_POINT
1640
                /* divide by 4*2, shift right 3 bits;
1641
                   extra halving to avoid overflows; it is ok, because result is normalized */
1642
                RE(tempLeft)  = RE(ps->ipd_prev[bk][i]) >> 3;
1643
                IM(tempLeft)  = IM(ps->ipd_prev[bk][i]) >> 3;
1644
                RE(tempRight) = RE(ps->opd_prev[bk][i]) >> 3;
1645
                IM(tempRight) = IM(ps->opd_prev[bk][i]) >> 3;
1646
#else
1647
46.3k
                RE(tempLeft)  = MUL_F(RE(ps->ipd_prev[bk][i]), FRAC_CONST(0.25));
1648
46.3k
                IM(tempLeft)  = MUL_F(IM(ps->ipd_prev[bk][i]), FRAC_CONST(0.25));
1649
46.3k
                RE(tempRight) = MUL_F(RE(ps->opd_prev[bk][i]), FRAC_CONST(0.25));
1650
46.3k
                IM(tempRight) = MUL_F(IM(ps->opd_prev[bk][i]), FRAC_CONST(0.25));
1651
46.3k
#endif
1652
1653
                /* save current value */
1654
46.3k
                RE(ps->ipd_prev[bk][i]) = ipdopd_cos_tab[abs(ps->ipd_index[env][bk])];
1655
46.3k
                IM(ps->ipd_prev[bk][i]) = ipdopd_sin_tab[abs(ps->ipd_index[env][bk])];
1656
46.3k
                RE(ps->opd_prev[bk][i]) = ipdopd_cos_tab[abs(ps->opd_index[env][bk])];
1657
46.3k
                IM(ps->opd_prev[bk][i]) = ipdopd_sin_tab[abs(ps->opd_index[env][bk])];
1658
1659
                /* add current value */
1660
#ifdef FIXED_POINT
1661
                /* extra halving to avoid overflows */
1662
                RE(tempLeft)  += RE(ps->ipd_prev[bk][i]) >> 1;
1663
                IM(tempLeft)  += IM(ps->ipd_prev[bk][i]) >> 1;
1664
                RE(tempRight) += RE(ps->opd_prev[bk][i]) >> 1;
1665
                IM(tempRight) += IM(ps->opd_prev[bk][i]) >> 1;
1666
#else
1667
46.3k
                RE(tempLeft)  += RE(ps->ipd_prev[bk][i]);
1668
46.3k
                IM(tempLeft)  += IM(ps->ipd_prev[bk][i]);
1669
46.3k
                RE(tempRight) += RE(ps->opd_prev[bk][i]);
1670
46.3k
                IM(tempRight) += IM(ps->opd_prev[bk][i]);
1671
46.3k
#endif
1672
1673
                /* ringbuffer index */
1674
46.3k
                if (i == 0)
1675
23.4k
                {
1676
23.4k
                    i = 2;
1677
23.4k
                }
1678
46.3k
                i--;
1679
1680
                /* get value before previous */
1681
#ifdef FIXED_POINT
1682
                /* dividing by 2*2, shift right 2 bits; extra halving to avoid overflows */
1683
                RE(tempLeft)  += (RE(ps->ipd_prev[bk][i]) >> 2);
1684
                IM(tempLeft)  += (IM(ps->ipd_prev[bk][i]) >> 2);
1685
                RE(tempRight) += (RE(ps->opd_prev[bk][i]) >> 2);
1686
                IM(tempRight) += (IM(ps->opd_prev[bk][i]) >> 2);
1687
#else
1688
46.3k
                RE(tempLeft)  += MUL_F(RE(ps->ipd_prev[bk][i]), FRAC_CONST(0.5));
1689
46.3k
                IM(tempLeft)  += MUL_F(IM(ps->ipd_prev[bk][i]), FRAC_CONST(0.5));
1690
46.3k
                RE(tempRight) += MUL_F(RE(ps->opd_prev[bk][i]), FRAC_CONST(0.5));
1691
46.3k
                IM(tempRight) += MUL_F(IM(ps->opd_prev[bk][i]), FRAC_CONST(0.5));
1692
46.3k
#endif
1693
1694
#if 0 /* original code */
1695
                ipd = (float)atan2(IM(tempLeft), RE(tempLeft));
1696
                opd = (float)atan2(IM(tempRight), RE(tempRight));
1697
1698
                /* phase rotation */
1699
                RE(phaseLeft) = (float)cos(opd);
1700
                IM(phaseLeft) = (float)sin(opd);
1701
                opd -= ipd;
1702
                RE(phaseRight) = (float)cos(opd);
1703
                IM(phaseRight) = (float)sin(opd);
1704
#else
1705
1706
                // x = IM(tempLeft)
1707
                // y = RE(tempLeft)
1708
                // p = IM(tempRight)
1709
                // q = RE(tempRight)
1710
                // cos(atan2(x,y)) = y/sqrt((x*x) + (y*y))
1711
                // sin(atan2(x,y)) = x/sqrt((x*x) + (y*y))
1712
                // cos(atan2(x,y)-atan2(p,q)) = (y*q + x*p) / ( sqrt((x*x) + (y*y)) * sqrt((p*p) + (q*q)) );
1713
                // sin(atan2(x,y)-atan2(p,q)) = (x*q - y*p) / ( sqrt((x*x) + (y*y)) * sqrt((p*p) + (q*q)) );
1714
1715
46.3k
                xy = magnitude_c(tempRight);
1716
46.3k
                pq = magnitude_c(tempLeft);
1717
1718
46.3k
                if (xy != 0)
1719
46.3k
                {
1720
46.3k
                    RE(phaseLeft) = DIV_F(RE(tempRight), xy);
1721
46.3k
                    IM(phaseLeft) = DIV_F(IM(tempRight), xy);
1722
46.3k
                } else {
1723
0
                    RE(phaseLeft) = 0;
1724
0
                    IM(phaseLeft) = 0;
1725
0
                }
1726
1727
46.3k
                xypq = MUL_F(xy, pq);
1728
1729
46.3k
                if (xypq != 0)
1730
46.3k
                {
1731
46.3k
                    real_t tmp1 = MUL_F(RE(tempRight), RE(tempLeft)) + MUL_F(IM(tempRight), IM(tempLeft));
1732
46.3k
                    real_t tmp2 = MUL_F(IM(tempRight), RE(tempLeft)) - MUL_F(RE(tempRight), IM(tempLeft));
1733
1734
46.3k
                    RE(phaseRight) = DIV_F(tmp1, xypq);
1735
46.3k
                    IM(phaseRight) = DIV_F(tmp2, xypq);
1736
46.3k
                } else {
1737
0
                    RE(phaseRight) = 0;
1738
0
                    IM(phaseRight) = 0;
1739
0
                }
1740
1741
46.3k
#endif
1742
1743
                /* MUL_F(COEF, REAL) = COEF */
1744
46.3k
                IM(h11) = MUL_F(RE(h11), IM(phaseLeft));
1745
46.3k
                IM(h12) = MUL_F(RE(h12), IM(phaseRight));
1746
46.3k
                IM(h21) = MUL_F(RE(h21), IM(phaseLeft));
1747
46.3k
                IM(h22) = MUL_F(RE(h22), IM(phaseRight));
1748
1749
46.3k
                RE(h11) = MUL_F(RE(h11), RE(phaseLeft));
1750
46.3k
                RE(h12) = MUL_F(RE(h12), RE(phaseRight));
1751
46.3k
                RE(h21) = MUL_F(RE(h21), RE(phaseLeft));
1752
46.3k
                RE(h22) = MUL_F(RE(h22), RE(phaseRight));
1753
46.3k
            }
1754
1755
            /* length of the envelope n_e+1 - n_e (in time samples) */
1756
            /* 0 < L <= 32: integer */
1757
361k
            L = (real_t)(ps->border_position[env + 1] - ps->border_position[env]);
1758
1759
            /* obtain final H_xy by means of linear interpolation */
1760
361k
            RE(deltaH11) = (RE(h11) - RE(ps->h11_prev[gr])) / L;
1761
361k
            RE(deltaH12) = (RE(h12) - RE(ps->h12_prev[gr])) / L;
1762
361k
            RE(deltaH21) = (RE(h21) - RE(ps->h21_prev[gr])) / L;
1763
361k
            RE(deltaH22) = (RE(h22) - RE(ps->h22_prev[gr])) / L;
1764
1765
361k
            RE(H11) = RE(ps->h11_prev[gr]);
1766
361k
            RE(H12) = RE(ps->h12_prev[gr]);
1767
361k
            RE(H21) = RE(ps->h21_prev[gr]);
1768
361k
            RE(H22) = RE(ps->h22_prev[gr]);
1769
361k
            IM(H11) = IM(H12) = IM(H21) = IM(H22) = 0;
1770
1771
361k
            RE(ps->h11_prev[gr]) = RE(h11);
1772
361k
            RE(ps->h12_prev[gr]) = RE(h12);
1773
361k
            RE(ps->h21_prev[gr]) = RE(h21);
1774
361k
            RE(ps->h22_prev[gr]) = RE(h22);
1775
1776
            /* only calculate imaginary part when needed */
1777
361k
            if ((ps->enable_ipdopd) && (bk < nr_ipdopd_par))
1778
46.3k
            {
1779
                /* obtain final H_xy by means of linear interpolation */
1780
46.3k
                IM(deltaH11) = (IM(h11) - IM(ps->h11_prev[gr])) / L;
1781
46.3k
                IM(deltaH12) = (IM(h12) - IM(ps->h12_prev[gr])) / L;
1782
46.3k
                IM(deltaH21) = (IM(h21) - IM(ps->h21_prev[gr])) / L;
1783
46.3k
                IM(deltaH22) = (IM(h22) - IM(ps->h22_prev[gr])) / L;
1784
1785
46.3k
                IM(H11) = IM(ps->h11_prev[gr]);
1786
46.3k
                IM(H12) = IM(ps->h12_prev[gr]);
1787
46.3k
                IM(H21) = IM(ps->h21_prev[gr]);
1788
46.3k
                IM(H22) = IM(ps->h22_prev[gr]);
1789
1790
46.3k
                if ((NEGATE_IPD_MASK & ps->map_group2bk[gr]) != 0)
1791
5.74k
                {
1792
5.74k
                    IM(deltaH11) = -IM(deltaH11);
1793
5.74k
                    IM(deltaH12) = -IM(deltaH12);
1794
5.74k
                    IM(deltaH21) = -IM(deltaH21);
1795
5.74k
                    IM(deltaH22) = -IM(deltaH22);
1796
1797
5.74k
                    IM(H11) = -IM(H11);
1798
5.74k
                    IM(H12) = -IM(H12);
1799
5.74k
                    IM(H21) = -IM(H21);
1800
5.74k
                    IM(H22) = -IM(H22);
1801
5.74k
                }
1802
1803
46.3k
                IM(ps->h11_prev[gr]) = IM(h11);
1804
46.3k
                IM(ps->h12_prev[gr]) = IM(h12);
1805
46.3k
                IM(ps->h21_prev[gr]) = IM(h21);
1806
46.3k
                IM(ps->h22_prev[gr]) = IM(h22);
1807
46.3k
            }
1808
1809
            /* apply H_xy to the current envelope band of the decorrelated subband */
1810
5.88M
            for (n = ps->border_position[env]; n < ps->border_position[env + 1]; n++)
1811
5.52M
            {
1812
                /* addition finalises the interpolation over every n */
1813
5.52M
                RE(H11) += RE(deltaH11);
1814
5.52M
                RE(H12) += RE(deltaH12);
1815
5.52M
                RE(H21) += RE(deltaH21);
1816
5.52M
                RE(H22) += RE(deltaH22);
1817
5.52M
                if ((ps->enable_ipdopd) && (bk < nr_ipdopd_par))
1818
692k
                {
1819
692k
                    IM(H11) += IM(deltaH11);
1820
692k
                    IM(H12) += IM(deltaH12);
1821
692k
                    IM(H21) += IM(deltaH21);
1822
692k
                    IM(H22) += IM(deltaH22);
1823
692k
                }
1824
1825
                /* channel is an alias to the subband */
1826
18.9M
                for (sb = ps->group_border[gr]; sb < maxsb; sb++)
1827
13.4M
                {
1828
13.4M
                    complex_t inLeft, inRight;  // precision_of in(Left|Right) == precision_of X_(left|right)
1829
1830
                    /* load decorrelated samples */
1831
13.4M
                    if (gr < ps->num_hybrid_groups)
1832
3.09M
                    {
1833
3.09M
                        RE(inLeft) =  RE(X_hybrid_left[n][sb]);
1834
3.09M
                        IM(inLeft) =  IM(X_hybrid_left[n][sb]);
1835
3.09M
                        RE(inRight) = RE(X_hybrid_right[n][sb]);
1836
3.09M
                        IM(inRight) = IM(X_hybrid_right[n][sb]);
1837
10.3M
                    } else {
1838
10.3M
                        RE(inLeft) =  RE(X_left[n][sb]);
1839
10.3M
                        IM(inLeft) =  IM(X_left[n][sb]);
1840
10.3M
                        RE(inRight) = RE(X_right[n][sb]);
1841
10.3M
                        IM(inRight) = IM(X_right[n][sb]);
1842
10.3M
                    }
1843
1844
                    /* precision_of temp(Left|Right) == precision_of X_(left|right) */
1845
1846
                    /* apply mixing */
1847
13.4M
                    RE(tempLeft) =  MUL_C(RE(H11), RE(inLeft)) + MUL_C(RE(H21), RE(inRight));
1848
13.4M
                    IM(tempLeft) =  MUL_C(RE(H11), IM(inLeft)) + MUL_C(RE(H21), IM(inRight));
1849
13.4M
                    RE(tempRight) = MUL_C(RE(H12), RE(inLeft)) + MUL_C(RE(H22), RE(inRight));
1850
13.4M
                    IM(tempRight) = MUL_C(RE(H12), IM(inLeft)) + MUL_C(RE(H22), IM(inRight));
1851
1852
                    /* only perform imaginary operations when needed */
1853
13.4M
                    if ((ps->enable_ipdopd) && (bk < nr_ipdopd_par))
1854
697k
                    {
1855
                        /* apply rotation */
1856
697k
                        RE(tempLeft)  -= MUL_C(IM(H11), IM(inLeft)) + MUL_C(IM(H21), IM(inRight));
1857
697k
                        IM(tempLeft)  += MUL_C(IM(H11), RE(inLeft)) + MUL_C(IM(H21), RE(inRight));
1858
697k
                        RE(tempRight) -= MUL_C(IM(H12), IM(inLeft)) + MUL_C(IM(H22), IM(inRight));
1859
697k
                        IM(tempRight) += MUL_C(IM(H12), RE(inLeft)) + MUL_C(IM(H22), RE(inRight));
1860
697k
                    }
1861
1862
                    /* store final samples */
1863
13.4M
                    if (gr < ps->num_hybrid_groups)
1864
3.09M
                    {
1865
3.09M
                        RE(X_hybrid_left[n][sb])  = RE(tempLeft);
1866
3.09M
                        IM(X_hybrid_left[n][sb])  = IM(tempLeft);
1867
3.09M
                        RE(X_hybrid_right[n][sb]) = RE(tempRight);
1868
3.09M
                        IM(X_hybrid_right[n][sb]) = IM(tempRight);
1869
10.3M
                    } else {
1870
10.3M
                        RE(X_left[n][sb])  = RE(tempLeft);
1871
10.3M
                        IM(X_left[n][sb])  = IM(tempLeft);
1872
10.3M
                        RE(X_right[n][sb]) = RE(tempRight);
1873
10.3M
                        IM(X_right[n][sb]) = IM(tempRight);
1874
10.3M
                    }
1875
13.4M
                }
1876
5.52M
            }
1877
1878
            /* shift phase smoother's circular buffer index */
1879
361k
            ps->phase_hist++;
1880
361k
            if (ps->phase_hist == 2)
1881
180k
            {
1882
180k
                ps->phase_hist = 0;
1883
180k
            }
1884
361k
        }
1885
178k
    }
1886
5.51k
}
1887
1888
void ps_free(ps_info *ps)
1889
8.80k
{
1890
    /* free hybrid filterbank structures */
1891
8.80k
    hybrid_free(ps->hyb);
1892
1893
8.80k
    faad_free(ps);
1894
8.80k
}
1895
1896
ps_info *ps_init(uint8_t sr_index, uint8_t numTimeSlotsRate)
1897
8.80k
{
1898
8.80k
    uint8_t i;
1899
8.80k
    uint8_t short_delay_band;
1900
1901
8.80k
    ps_info *ps = (ps_info*)faad_malloc(sizeof(ps_info));
1902
8.80k
    memset(ps, 0, sizeof(ps_info));
1903
1904
8.80k
    ps->hyb = hybrid_init(numTimeSlotsRate);
1905
8.80k
    ps->numTimeSlotsRate = numTimeSlotsRate;
1906
1907
8.80k
    ps->ps_data_available = 0;
1908
1909
    /* delay stuff*/
1910
8.80k
    ps->saved_delay = 0;
1911
1912
572k
    for (i = 0; i < 64; i++)
1913
563k
    {
1914
563k
        ps->delay_buf_index_delay[i] = 0;
1915
563k
    }
1916
1917
35.2k
    for (i = 0; i < NO_ALLPASS_LINKS; i++)
1918
26.4k
    {
1919
26.4k
        ps->delay_buf_index_ser[i] = 0;
1920
#ifdef PARAM_32KHZ
1921
        if (sr_index <= 5) /* >= 32 kHz*/
1922
        {
1923
            ps->num_sample_delay_ser[i] = delay_length_d[1][i];
1924
        } else {
1925
            ps->num_sample_delay_ser[i] = delay_length_d[0][i];
1926
        }
1927
#else
1928
26.4k
        (void)sr_index;
1929
        /* THESE ARE CONSTANTS NOW */
1930
26.4k
        ps->num_sample_delay_ser[i] = delay_length_d[i];
1931
26.4k
#endif
1932
26.4k
    }
1933
1934
#ifdef PARAM_32KHZ
1935
    if (sr_index <= 5) /* >= 32 kHz*/
1936
    {
1937
        short_delay_band = 35;
1938
        ps->nr_allpass_bands = 22;
1939
        ps->alpha_decay = FRAC_CONST(0.76592833836465);
1940
        ps->alpha_smooth = FRAC_CONST(0.25);
1941
    } else {
1942
        short_delay_band = 64;
1943
        ps->nr_allpass_bands = 45;
1944
        ps->alpha_decay = FRAC_CONST(0.58664621951003);
1945
        ps->alpha_smooth = FRAC_CONST(0.6);
1946
    }
1947
#else
1948
    /* THESE ARE CONSTANTS NOW */
1949
8.80k
    short_delay_band = 35;
1950
8.80k
    ps->nr_allpass_bands = 22;
1951
8.80k
    ps->alpha_decay = FRAC_CONST(0.76592833836465);
1952
8.80k
    ps->alpha_smooth = FRAC_CONST(0.25);
1953
8.80k
#endif
1954
1955
    /* THESE ARE CONSTANT NOW IF PS IS INDEPENDANT OF SAMPLERATE */
1956
316k
    for (i = 0; i < short_delay_band; i++)
1957
308k
    {
1958
308k
        ps->delay_D[i] = 14;
1959
308k
    }
1960
264k
    for (i = short_delay_band; i < 64; i++)
1961
255k
    {
1962
255k
        ps->delay_D[i] = 1;
1963
255k
    }
1964
1965
    /* mixing and phase */
1966
448k
    for (i = 0; i < 50; i++)
1967
440k
    {
1968
440k
        RE(ps->h11_prev[i]) = 1;
1969
440k
        IM(ps->h11_prev[i]) = 1;
1970
440k
        RE(ps->h12_prev[i]) = 1;
1971
440k
        IM(ps->h12_prev[i]) = 1;
1972
440k
    }
1973
1974
8.80k
    ps->phase_hist = 0;
1975
1976
184k
    for (i = 0; i < 20; i++)
1977
176k
    {
1978
176k
        RE(ps->ipd_prev[i][0]) = 0;
1979
176k
        IM(ps->ipd_prev[i][0]) = 0;
1980
176k
        RE(ps->ipd_prev[i][1]) = 0;
1981
176k
        IM(ps->ipd_prev[i][1]) = 0;
1982
176k
        RE(ps->opd_prev[i][0]) = 0;
1983
176k
        IM(ps->opd_prev[i][0]) = 0;
1984
176k
        RE(ps->opd_prev[i][1]) = 0;
1985
176k
        IM(ps->opd_prev[i][1]) = 0;
1986
176k
    }
1987
1988
8.80k
    return ps;
1989
8.80k
}
1990
1991
/* main Parametric Stereo decoding function */
1992
uint8_t ps_decode(ps_info *ps, qmf_t X_left[38][64], qmf_t X_right[38][64])
1993
5.51k
{
1994
5.51k
    qmf_t X_hybrid_left[32][32] = {{{0}}};
1995
5.51k
    qmf_t X_hybrid_right[32][32] = {{{0}}};
1996
1997
    /* delta decoding of the bitstream data */
1998
5.51k
    ps_data_decode(ps);
1999
2000
    /* set up some parameters depending on filterbank type */
2001
5.51k
    if (ps->use34hybrid_bands)
2002
2.03k
    {
2003
2.03k
        ps->group_border = (uint8_t*)group_border34;
2004
2.03k
        ps->map_group2bk = (uint16_t*)map_group2bk34;
2005
2.03k
        ps->num_groups = 32+18;
2006
2.03k
        ps->num_hybrid_groups = 32;
2007
2.03k
        ps->nr_par_bands = 34;
2008
2.03k
        ps->decay_cutoff = 5;
2009
3.47k
    } else {
2010
3.47k
        ps->group_border = (uint8_t*)group_border20;
2011
3.47k
        ps->map_group2bk = (uint16_t*)map_group2bk20;
2012
3.47k
        ps->num_groups = 10+12;
2013
3.47k
        ps->num_hybrid_groups = 10;
2014
3.47k
        ps->nr_par_bands = 20;
2015
3.47k
        ps->decay_cutoff = 3;
2016
3.47k
    }
2017
2018
    /* Perform further analysis on the lowest subbands to get a higher
2019
     * frequency resolution
2020
     */
2021
5.51k
    hybrid_analysis((hyb_info*)ps->hyb, X_left, X_hybrid_left,
2022
5.51k
        ps->use34hybrid_bands, ps->numTimeSlotsRate);
2023
2024
    /* decorrelate mono signal */
2025
5.51k
    ps_decorrelate(ps, X_left, X_right, X_hybrid_left, X_hybrid_right);
2026
2027
    /* apply mixing and phase parameters */
2028
5.51k
    ps_mix_phase(ps, X_left, X_right, X_hybrid_left, X_hybrid_right);
2029
2030
    /* hybrid synthesis, to rebuild the SBR QMF matrices */
2031
5.51k
    hybrid_synthesis((hyb_info*)ps->hyb, X_left, X_hybrid_left,
2032
5.51k
        ps->use34hybrid_bands, ps->numTimeSlotsRate);
2033
2034
5.51k
    hybrid_synthesis((hyb_info*)ps->hyb, X_right, X_hybrid_right,
2035
5.51k
        ps->use34hybrid_bands, ps->numTimeSlotsRate);
2036
2037
5.51k
    return 0;
2038
5.51k
}
2039
2040
#endif