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