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