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