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