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