Mercurial > mplayer.hg

comparison libfaad2/ps_dec.c @ 13453:6d50ef45a058

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