Line data Source code
1 : /******************************************************************************
2 : * ETSI TS 103 634 V1.5.1 *
3 : * Low Complexity Communication Codec Plus (LC3plus) *
4 : * *
5 : * Copyright licence is solely granted through ETSI Intellectual Property *
6 : * Rights Policy, 3rd April 2019. No patent licence is granted by implication, *
7 : * estoppel or otherwise. *
8 : ******************************************************************************/
9 :
10 : #include "defines.h"
11 : #include "functions.h"
12 : #include "math.h" /*dbg*/
13 :
14 : /*---------------------------------------------------------------------*
15 : * Local constants
16 : *---------------------------------------------------------------------*/
17 :
18 : #define DELTA_CORR 5 /* Tuning parameter - defining range for phase correction around peak */
19 : #define DELTA_CORR_F0_INT 2 /* Tuning parameter - defining range for phase correction around peak */
20 :
21 : #define MAX_INCREASE_GRPPOW_FX 0 /* max. amplification in case of transients (in dB scale) */
22 :
23 : #define ONE_SIDED_SINE_WIDTH (4) /* expected pure sine main lobe maximum width (4+1+4) bins *62.5 hz/bin => approx 560 Hz total width */
24 : #define SIDE_LIM 12539859L /* 10^ (4.5/20.0) = 2^(a); --> x= 0.747433821 -> Lx_Q24 = round((1L<<(24))*0.747433821)) = 12539859 */
25 : #define LFHF_LIM 16719812L /* 10^ (6.0/20.0) = 2^(b); --> x= 0.996578428 -> Lx_Q24 = round((1L<<(24))*0.996578428)) = 16719812 */
26 :
27 : #define CMPLMNT_PLOC_SENS_FX 2294 /* (1.0 - p_locator_sens) in Q15 */
28 : #define FEC_HQ_ECU_ROOT2 23170 /*(0x5a83) */ /* sqrt(2) in Q14 */
29 : #define FEC_TWOTHIRDS_Q15 21845 /* round(2^15*2/3) */
30 :
31 : static void get_sin_cosQ10opt(Word16 phase, /* Q10 0..1023 i.e. 1024=2*pi */
32 : Word16 *ptrSin, /* Q15 */
33 : Word16 *ptrCos); /* Q15 */
34 :
35 : static Word16 sqrt2ndOrder(const Word16);
36 :
37 : void my_wtda_fx(const Word16 *new_audio, /* i : input audio to be windowed Q0 20 ms , OPT can be output as well */
38 : const Word16 *const win2ms_init, /* i: 2 ms initial part of pre_tda window */
39 : const Word16 *const win16ms_center, /* i: 16 ms combined part of pre_tda IWHR+MDCT-ana */
40 : Word32 * L_wtda_audio, /* o : tda audio Q16 20 ms */
41 : const Word16 L, Word8 *scratchBuffer);
42 :
43 : static void windowing_L(const Word16 *, Word32 *, const Word16 *, const Word16, const Word16);
44 : static void windowing_ola(const Word16 *, Word16 *, const Word16 *, const Word16);
45 : static void ola_add(const Word16 *, const Word16 *, Word16 *, const Word16);
46 : static void intlvW32_2_flippedW16(Word32 *L_x, const Word16 numPairs, const Word16 L_prot, Word16 *x);
47 : static void flippedW16_2_intlvW32(Word16 *x, const Word16 numPairs, const Word16 Lprot, Word32 *L_x);
48 : static Word16 imax_fx(const Word16 *, const Word16);
49 :
50 :
51 : Word16 rand_phase_fx(const Word16 seed, Word16 *sin_F, Word16 *cos_F);
52 :
53 : static Word16 imax2_jacobsen_mag_fx(const Word16 *y_re, const Word16 *y_im, const Word16 special);
54 : static void fft_spec2_sqrt_approx_fx(const Word16 x[], Word16 xMagSqrt[], const Word16 N);
55 : static Word16 sqrtMagnApprox_fx(const Word16 re, const Word16 im);
56 :
57 : static Word16 plc_phEcu_nonpure_tone_ana_fx(const Word16* plocs, const Word16 n_plocs, const Word16* X,
58 : const Word32 *L_Xavg, /* i : Frequency group amp averages for tonal tilt analysis pref. Max upshifted */
59 : const Word16 Lprot, const Word16 fs_idx);
60 :
61 0 : static void rotate_W16_fx(Word16 re_in, Word16 im_in, Word16 cosFactor, Word16 sinFactor, Word16 *re_out_ptr,
62 : Word16 *im_out_ptr)
63 : {
64 : #ifdef WMOPS
65 : push_wmops("PhECU::rotate_W16_fx");
66 : #endif
67 0 : *re_out_ptr = msu_r(L_mult(re_in, cosFactor), im_in, sinFactor); /* 2 ops no move when inlined */
68 0 : *im_out_ptr = mac_r(L_mult(re_in, sinFactor), im_in, cosFactor); /* 2 ops no move when inlined */
69 : #ifdef WMOPS
70 : pop_wmops();
71 : #endif
72 0 : return;
73 : }
74 :
75 0 : static void valley_magnitude_adj_fx(Word16 *re_ptr, Word16 *im_ptr, Word16 uniFactor, Word16 cosFactor)
76 : {
77 : Word16 scale_fx;
78 : #ifdef WMOPS
79 : push_wmops("PhECU::valley_magnitude_adj_fx");
80 : #endif
81 :
82 : /* y = 0.5*((2*rand(1,10000) + 1*cos(2*pi*x)) - 1 */ /* y will be in -1 to 1 range */
83 : /* y = 1*((1*rand(1,10000) + 0.5*cos(2*pi*x)) - 1 */ /* y will be in -1 to 1 range */
84 :
85 0 : scale_fx /*Q15*/ = mac_r(L_mult(uniFactor, 16384), cosFactor, 16384);
86 : /* make gain distribution more like N(0,1) than uniform */
87 :
88 0 : scale_fx /*Q14*/ = round_fx(L_mac((Word32)(16384L << 16), scale_fx, 4096));
89 : /* create a random gain scaling value with mean 1.0 and max 1.25 and min 0.75 */
90 0 : ASSERT(scale_fx <= (16384 + 8192));
91 0 : ASSERT(scale_fx >= (-16384 - 8192));
92 0 : *re_ptr = mult_r(scale_fx, shl_sat(*re_ptr, 1)); /* no moves , should be inlined */
93 0 : *im_ptr = mult_r(scale_fx, shl_sat(*im_ptr, 1)); /* no moves , should be inlined */
94 :
95 : #ifdef WMOPS
96 : pop_wmops();
97 : #endif
98 0 : return;
99 : }
100 :
101 : /*------------------------------------------------------------------*
102 : * rand_phase()
103 : *
104 : * randomized phase in form of sin and cos components
105 : *------------------------------------------------------------------*/
106 0 : Word16 rand_phase_fx(const Word16 seed, Word16 *sin_F, Word16 *cos_F)
107 : {
108 : /* 4x8+8 lookup scheme requiring ~40 Words of ROM freqRes 90/8 = 11.25 degrees */
109 :
110 : /* x=(0:(5*8-1))*(2*pi)/32; y=sin(x);y_int=max(-32768,min(32767,round(y*32768))), y_int/32768 */
111 :
112 0 : const Word16 *sincos_lowres_tab_cosQ15_fx = sincos_lowres_tab_sinQ15_fx + 8;
113 : /* position at 90 degrees , ptr init */
114 : Word16 seed2;
115 : Word16 seed2_shift;
116 :
117 : #ifdef DYNMEM_COUNT
118 : Dyn_Mem_In("rand_phase_fx", sizeof(struct {
119 : const Word16 *sincos_lowres_tab_cosQ15_fx; /* position at 90 degrees */
120 : Word16 seed2; /* 16 bit signed */
121 : Word16 seed2_shift;
122 : }));
123 : #endif
124 :
125 : #ifdef WMOPS
126 : push_wmops("PhECU::rand_phase_fx");
127 : #endif
128 :
129 0 : seed2 = extract_l(L_mac0(13849, seed, 31821));
130 0 : seed2_shift = lshr(seed2, 11);
131 : /* logical shift to get uniform random 5 msb bits; 0-31 , 0 degrees to 31*360/32= 348.75 */
132 0 : *sin_F = sincos_lowres_tab_sinQ15_fx[seed2_shift]; move16(); /* these moves can often be avoided by returning seed2shift and inlining */
133 0 : *cos_F = sincos_lowres_tab_cosQ15_fx[seed2_shift]; move16(); /* these moves can often be avoided by inlining */
134 : /* total WC 5 ops */
135 : #ifdef DYNMEM_COUNT
136 : Dyn_Mem_Out();
137 : #endif
138 : #ifdef WMOPS
139 : pop_wmops();
140 : #endif
141 0 : return seed2;
142 : }
143 :
144 : /*-----------------------------------------------------------------------------
145 : * trans_burst_ana_fx()
146 : *
147 : * Transient analysis
148 : *----------------------------------------------------------------------------*/
149 0 : void trans_burst_ana_fx(
150 : const Word16 *xfp, /* i : Input signal (, only used if time_offset==0) now in up_scaled *Q_spec */
151 : Word16 * mag_chg, /* o : Magnitude modification vector Q15 */
152 : Word16 * ph_dith, /* o : Phase dither, 2*PI is not included (Q15, i.e., between 0.0 and 1.0) */
153 : Word16 * mag_chg_1st, /* i/o: per band magnitude modifier for transients Q15 */
154 : const Word16 output_frame, /* i : Frame length */
155 : const Word16 time_offs, /* i : Time offset (integral multiple of output_frame) */
156 : const Word16 est_stab_content, /* i : 0.0=dynamic ... 1.0=Stable (==st->env_stab ) */
157 : Word16 * alpha, /* o : Magnitude modification factors for fade to average */
158 : Word16 * beta, /* : Magnitude modification factors for fade to average */
159 : Word16 * beta_mute, /* i/o : Factor for long-term mute */
160 : Word16 * Xavg, /* o : Frequency group average gain to fade to in same Q as X_sav */
161 : Word16 Q_spec, Word32 L_oold_xfp_w_E_fx, Word16 oold_xfp_w_E_exp_fx, Word16 oold_Ltot_exp_fx,
162 : Word16 *oold_grp_shape_fx,
163 :
164 : Word32 L_old_xfp_w_E_fx, Word16 old_xfp_w_E_exp_fx, Word16 old_Ltot_exp_fx, Word16 *old_grp_shape_fx,
165 : Word16 fadeout,
166 : Word32 * L_Xavg, /* full scale band amplitudes */
167 : Word8 *scratchBuffer /* Size = 4*4 * MAX_LTRANA + (2*4 + 1*2) * MAX_LGW + 8 */
168 : )
169 : {
170 : Word16 att_val, attDegreeFrames;
171 : Word32 * L_pGrPowLeft, *L_pGrPowRight;
172 : Word32 * L_gr_pow_left, *L_gr_pow_right;
173 : Word16 Lgw, i, k, burst_len;
174 : Word16 man, expo;
175 0 : Word16 att_always = 0; /* fixed attenuation per frequency group if set to 1 */
176 : Word16 oneOverFrame, roundEstMusContent, tmp16;
177 : Word16 burst_att_thresh;
178 : Word16 att_per_frame;
179 : Word16 * tr_dec;
180 : Word32 L_acc;
181 : Word16 fs_scale;
182 : Word16 scale_sh;
183 :
184 : Word32 L_oold_tmp, L_old_tmp;
185 : Word16 oold_exp_fx, old_exp_fx;
186 : Word16 margin_oold, margin_old;
187 : Word16 fs_idx;
188 : Word16 exp_diff;
189 : Word16 Xavg_exp_fx, Xavg_mod_exp_fx;
190 : Word16 tr_rise[MAX_LGW];
191 : Word16 tr_decay[MAX_LGW];
192 : Word16 man_in, expo_in, tmp;
193 : Word32 L_tmp, L_tmp2;
194 : Word16 thresh_tr_rise_lin_Q15;
195 : Word16 thresh_tr_decay_lin_Q15;
196 : Word16 beta_mute_thr;
197 : Word16 fade_ms_ind;
198 :
199 : #ifdef DYNMEM_COUNT
200 : Dyn_Mem_In("trans_burst_ana_fx", sizeof(struct {
201 :
202 : Word16 att_val, attDegreeFrames;
203 : Word32 * pGrPowLeft_L, *pGrPowRight_L;
204 : Word32 * L_gr_pow_left, *L_gr_pow_right;
205 : Word16 Lprot;
206 : Word16 Lgw, i, k, burst_len;
207 : Word16 man, expo;
208 : Word16 att_always; /* fixed attenuation per frequency group if set to 1 */
209 : Word16 oneOverFrame, roundEstMusContent, tmp16;
210 :
211 : Word16 burst_att_thresh;
212 : Word16 att_per_frame;
213 :
214 : Word16 * tr_dec;
215 : UWord16 lsb;
216 : Word32 L_acc;
217 : Word16 fs_scale;
218 : Word16 scale_sh;
219 :
220 : Word32 L_oold_tmp;
221 : Word32 L_old_tmp;
222 : Word16 fs_idx;
223 : Word16 shift32;
224 : Word16 margin_old;
225 : Word16 margin_oold;
226 :
227 : Word16 Xavg_exp_fx, Xavg_mod_exp_fx;
228 : Word16 tr_rise[MAX_LGW];
229 : Word16 tr_decay[MAX_LGW];
230 :
231 : Word16 man_in, expo_in, tmp;
232 : Word32 L_tmp, L_tmp2;
233 : Word16 thresh_tr_rise_lin_Q15;
234 : Word16 thresh_tr_decay_lin_Q15;
235 :
236 : Word16 beta_mute_thr;
237 : Word16 fade_ms_ind;
238 : }));
239 : #endif
240 :
241 : UNUSED(xfp);
242 : UNUSED(oold_xfp_w_E_exp_fx);
243 : UNUSED(old_xfp_w_E_exp_fx);
244 :
245 : if (time_offs == 0)
246 : {
247 : #ifdef WMOPS
248 : push_wmops("PhECU::trans_burst_ana_fx(1st)");
249 : #endif
250 : }
251 : else
252 : {
253 : #ifdef WMOPS
254 : push_wmops("PhECU::trans_burst_ana_fx(N)");
255 : #endif
256 : }
257 :
258 0 : fs_idx = mult(output_frame, (Word16)(32768.0 / 99.0)); /* truncation needed , i.e no rounding can be applied here */
259 0 : ASSERT(fs_idx == (output_frame / 100));
260 :
261 0 : L_gr_pow_left = (Word32 *)scratchAlign(scratchBuffer, 0); /* Size = 4 * MAX_LGW */ /* Size = 4 * MAX_LGW */
262 :
263 0 : L_gr_pow_right = (Word32 *)scratchAlign(L_gr_pow_left, sizeof(*L_gr_pow_left) * MAX_LGW); /* Size = 4 * MAX_LGW */
264 :
265 0 : tr_dec = (Word16 *)scratchAlign(L_gr_pow_right, sizeof(*L_gr_pow_right) * MAX_LGW); /* Size = 2bytes * MAX_LGW */
266 :
267 0 : oneOverFrame = oneOverFrameQ15Tab[fs_idx];
268 0 : Lgw = s_min(add(fs_idx, LGW8K), LGW48K); /* 4,5,6,7, (7/8) */
269 :
270 0 : burst_len = add(mult_r(time_offs, oneOverFrame), 1);
271 :
272 : UNUSED(est_stab_content);
273 : UNUSED(roundEstMusContent);
274 :
275 0 : move16();
276 0 : fade_ms_ind = (PLC2_FADEOUT_IN_MS - PLC2_FADEOUT_IN_MS_MIN) / PLC2_FADEOUT_RES; /* a shorter fading entry in fade_scheme_tab_fx */
277 0 : test();
278 0 : if (fadeout != 0)
279 : {
280 0 : fade_ms_ind = (PLC2_FADEOUT_LONG_IN_MS - PLC2_FADEOUT_IN_MS_MIN) / PLC2_FADEOUT_RES; move16();
281 : /* a long fading table entry in fade_scheme_tab */
282 : }
283 :
284 0 : move16(); move16(); move16();
285 0 : att_per_frame = fade_scheme_tab_fx[fade_ms_ind][0];
286 0 : burst_att_thresh = fade_scheme_tab_fx[fade_ms_ind][1]; /* number of 1.0 frames before muting phase */
287 : /* band gain muting can take place earlier due to a band transient */
288 0 : beta_mute_thr = fade_scheme_tab_fx[fade_ms_ind][2]; /* faster muting of added noise starts when slow main signal fadeout is over */
289 :
290 : #ifdef PLC_FADEOUT_IN_MS
291 0 : ASSERT(att_per_frame >= 1 && att_per_frame <= 12); /* table based lookup restriction */
292 : #else
293 : ASSERT(att_per_frame == 1 || att_per_frame == 2); /* table based lookup restriction */
294 : #endif
295 :
296 0 : *ph_dith = 0; /* peak scrambling, not in use */
297 :
298 0 : attDegreeFrames = 0; move16();
299 0 : IF(sub(burst_len, burst_att_thresh) > 0)
300 : {
301 0 : att_always = 1; move16();
302 : /* increase degree of attenuation */
303 :
304 : /* N.B. To facilitate the subsequent 10^(-att_degree/20) implementation
305 : * so as to use direct table-lookup,
306 : * the first (burstLen - burst_att_thresh) are NOT multiplied by "att_per_frame". */
307 0 : attDegreeFrames = sub(burst_len, burst_att_thresh); /* multiplied by 1.0 , */
308 : /* Furthermore, in order to minimize the size of the lookup-table required to
309 : * implement 10^(-att_degree/10), hard limit attDegreeFrames to (30% of 100)=30.
310 : * If attDegreeFrames is greater than 30, it means there are more than 30 successive
311 : * bad frames. */
312 0 : if (sub(attDegreeFrames, OFF_FRAMES_LIMIT) > 0)
313 : {/* Hard limit the no. of frames, for table lookup */
314 0 : attDegreeFrames = OFF_FRAMES_LIMIT; move16();
315 : }
316 : }
317 :
318 0 : plc_phEcu_initWord16(alpha, 32767, MAX_LGW);
319 0 : basop_memset(beta, 0, (MAX_LGW) * sizeof(Word16));
320 0 : IF(sub(burst_len, 1) <= 0)
321 : {
322 0 : *beta_mute = BETA_MUTE_FAC_INI; move16();
323 0 : *beta_mute = shr_pos(*beta_mute , 1); /* perceptual decrease */
324 : }
325 :
326 0 : IF(sub(burst_len, 1) <= 0)
327 : {
328 0 : L_pGrPowLeft = &L_gr_pow_left[0]; /* ptr init*/
329 0 : L_pGrPowRight = &L_gr_pow_right[0]; /* ptr init*/
330 :
331 0 : fs_scale = xfp_wE_MDCT2FFTQ11[fs_idx]; move16();
332 0 : scale_sh = 4; /* 15-11 */ move16();
333 : /* L_*old_xfp_w_E_fx, always upscaled to max from the calculating function */
334 :
335 :
336 0 : L_oold_tmp = Mpy_32_16_lc3plus(L_oold_xfp_w_E_fx, fs_scale);
337 0 : L_old_tmp = Mpy_32_16_lc3plus(L_old_xfp_w_E_fx, fs_scale);
338 :
339 0 : oold_exp_fx = add(oold_Ltot_exp_fx, scale_sh);
340 0 : old_exp_fx = add(old_Ltot_exp_fx, scale_sh);
341 :
342 : /*re-normalize L_mantissas and adjust exps */
343 0 : margin_oold = norm_l(L_oold_tmp);
344 0 : L_oold_tmp = L_shl_pos(L_oold_tmp, margin_oold);
345 0 : oold_exp_fx = sub(oold_exp_fx, margin_oold);
346 :
347 0 : margin_old = norm_l(L_old_tmp);
348 0 : L_old_tmp = L_shl_pos(L_old_tmp, margin_old);
349 0 : old_exp_fx = sub(old_exp_fx, margin_old);
350 :
351 : /* now time to analyze how the actual L_tot exponent scaling should be done */
352 : /* bring up the lowest exp to the same exp as the higher exp, and scale down the corresponding mantissa */
353 0 : exp_diff = sub(old_exp_fx, oold_exp_fx); /* energy increase from oold to old in log2 shifts */
354 :
355 : /* Overflow2 fix */
356 0 : exp_diff = s_max(-31, exp_diff);
357 0 : exp_diff = s_min(31, exp_diff);
358 0 : if (exp_diff > 0)
359 : { /* oold_exp < old_exp */
360 : /* old_exp is limiting, shift down oold mantissa */
361 0 : L_oold_tmp = L_shr_pos(L_oold_tmp, exp_diff);
362 : }
363 0 : if (exp_diff < 0)
364 : { /* oold_exp > old_exp */
365 : /* oold_exp is limiting, shift down old mantissa */
366 0 : L_old_tmp = L_shr_pos(L_old_tmp, negate(exp_diff));
367 : }
368 0 : oold_exp_fx = s_max(oold_exp_fx, old_exp_fx);
369 0 : old_exp_fx = oold_exp_fx; move16();
370 :
371 : /* safety set lowest energy to 2 , as one bit is shifted away in avg calculation */
372 0 : L_oold_tmp = L_max(L_oold_tmp, 2L);
373 0 : L_old_tmp = L_max(L_old_tmp, 2L);
374 :
375 0 : FOR(k = 0; k < Lgw; k++) /* NB Lgw may be shorter than all defined bands , e.g at at 48k */
376 : {
377 0 : L_gr_pow_left[k] = Mpy_32_16_lc3plus(L_oold_tmp, oold_grp_shape_fx[k]); move32();
378 0 : L_gr_pow_right[k] = Mpy_32_16_lc3plus(L_old_tmp, old_grp_shape_fx[k]); move32();
379 :
380 : /*Xavg[k] = sqrt(0.5f*(gr_pow_left[k]+gr_pow_right[k])/(float)(gw[k+1]-gw[k]));*/
381 0 : Xavg_exp_fx = sub(old_exp_fx, 1); /* virtual pre divide X_avg by 2 too keep precision in summation */
382 0 : L_acc = L_add(L_shr_pos(L_gr_pow_left[k], 1), L_shr_pos(L_gr_pow_right[k], 1));
383 0 : L_acc = L_shr_pos(L_acc, gw_len_inv_shift_fx[k]); /* divide by (bandwidth/2) in bins */
384 :
385 : { /* new Xavg_fx calculation */
386 0 : L_acc = L_max(L_acc, 1L);
387 0 : tmp = norm_l(L_acc);
388 0 : Xavg_exp_fx = sub(Xavg_exp_fx, tmp);
389 0 : L_acc = L_shl_pos(L_acc, tmp); /* now between 0.5 an 1.0*/
390 :
391 0 : expo_in = add(Xavg_exp_fx, 0);
392 0 : man_in = round_fx_sat(L_acc);
393 :
394 : /* now allow both positive and negative expos into sqrt */
395 0 : man = sqrt2ndOrder(man_in);
396 0 : if (s_and(expo_in, 1) != 0)
397 : {
398 0 : man = mult_r(man, FEC_HQ_ECU_ROOT2); /* odd exp operation */ /* 1/sqrt(2) */
399 : }
400 0 : expo = shr_r(expo_in, 1); /* apply even part of exp , square root operation. shr_r needed for positive side exps */
401 :
402 :
403 :
404 0 : L_acc = L_deposit_h(man);
405 0 : Xavg_exp_fx = add(expo, 0);
406 : /*Note: sqrt approximaton may overshoot e.g- sqrt(1.0) may become 1.0001 i.e. saturation is needed when eventually applying expo */
407 :
408 :
409 : /* Xavg (unscaled flt in L_acc*2^(exp-31)) needs to be saved in the same scale + Q as the stored 16bit
410 : * Xsav_fx, for use in subst_spec() */
411 : /* move Xavg fft scale to fx domain fx-fft scale*/
412 0 : L_acc = Mpy_32_16_lc3plus(L_acc, PhEcu_Xsav_Flt2FxScaleQ15[fs_idx]); /* fs fixed fractional change */
413 0 : Xavg_mod_exp_fx = sub(Xavg_exp_fx, PhEcu_Xsav_Flt2FxDnShift[fs_idx]); /* fs fixed exp change*/
414 0 : Xavg_mod_exp_fx = add(Xavg_mod_exp_fx, Q_spec); /* signal adaptive exp change*/
415 :
416 : /* :: move to Q_spec domain of Xsav , Q fixed in first BFI frames */
417 :
418 : /* extract Q0 value shift so that the mantissa is in the high part with man*2,^(0-15) */
419 0 : exp_diff = sub(15, Xavg_mod_exp_fx);
420 :
421 :
422 0 : exp_diff = s_min(31, exp_diff); /* limit to meaningfull DSP shifts as described by up to 6 bits */
423 0 : exp_diff = s_max(-32, exp_diff);
424 0 : if (exp_diff > 0)
425 : {
426 0 : L_acc = L_shr_pos(L_acc, exp_diff); /* may underflow */
427 : }
428 :
429 0 : if (exp_diff < 0)
430 : {
431 0 : L_acc = L_shr_sat(L_acc, exp_diff);
432 : }
433 0 : L_Xavg[k] = L_acc; /* export full 31 bit scale band amplitude */
434 0 : Xavg[k] = round_fx_sat(L_acc); /* extract high part */
435 :
436 : } /*end of new Xavg_fx calculation */
437 : /* internal transition detection */
438 :
439 : { /* pure percentage based transient detection */
440 0 : thresh_tr_rise_lin_Q15 = PhEcu_frac_thr_rise_lin_Q15[k];
441 0 : thresh_tr_decay_lin_Q15 = PhEcu_frac_thr_decay_lin_Q15[k];
442 :
443 : /* analyse rise */
444 : /* one of L_left or L_right should be pre-upshifted to a near max mantissa, (in one band ) */
445 0 : L_tmp2 = L_deposit_h(0);
446 0 : L_tmp = Mpy_32_16_lc3plus(*L_pGrPowRight, thresh_tr_rise_lin_Q15);
447 0 : if (L_sub(*L_pGrPowLeft, L_tmp) <= 0)
448 : {
449 0 : L_tmp2 = L_deposit_l(1);
450 : }
451 :
452 0 : if (*L_pGrPowLeft == 0) /* denominator zero special cases */
453 : {
454 : /* rise: Right/Left ; " * / 0 " --> tr_rise=1 ; "0/0" --> tr_rise = 0 */
455 0 : L_tmp2 = L_min(*L_pGrPowRight, 1L);
456 : }
457 0 : tr_rise[k] = extract_l(L_tmp2); move16();
458 :
459 : /* analyse decay */
460 0 : L_tmp2 = L_deposit_h(0);
461 0 : L_tmp = Mpy_32_16_lc3plus(*L_pGrPowLeft, thresh_tr_decay_lin_Q15);
462 0 : if (L_sub(L_tmp, *L_pGrPowRight) >= 0)
463 : {
464 0 : L_tmp2 = L_deposit_l(1);
465 : }
466 0 : if (*L_pGrPowRight == 0) /* right side no energy , special cases */
467 : {
468 : /* decay: Right/Left ; " 0 / * " --> tr_decay=0 ; "0/0" --> tr_decay = 0 */
469 0 : L_tmp2 = L_deposit_h(0);
470 : }
471 0 : tr_decay[k] = extract_l(L_tmp2); move16();
472 :
473 0 : tr_dec[k] = s_max(tr_rise[k], tr_decay[k]); move16();
474 :
475 : } /* percentage tr_dec */
476 : /* magnitude modification, calculated for decay only */
477 0 : IF(add(tr_dec[k], att_always) != 0)
478 : {
479 :
480 : #if MAX_INCREASE_GRPPOW_FX != 0
481 : #error trans_burst_ana_fx-- The following implementation is incorrect
482 : #endif
483 0 : att_val = 32767; move16();
484 0 : IF(L_sub(*L_pGrPowRight, 0) > 0)
485 : {
486 0 : IF(L_sub(*L_pGrPowRight, *L_pGrPowLeft) < 0) /* decay , i.e., (gr_pow_right/gr_pow_left) < 1.0 */
487 : {
488 : /* Compute sqrt(grp_pow_chg), where grp_pow_chg = gr_pow_right/gr_pow_left. */
489 0 : tmp16 = plc_phEcu_ratio_fx(*L_pGrPowRight, *L_pGrPowLeft, &expo); /* output tmp16 in Q14 */
490 :
491 0 : expo = sub(expo, (15 - 14)); /* Now, tmp16 is considered in Q15 */
492 0 : i = norm_s(tmp16);
493 0 : man = shl_pos(tmp16, i); /* Mandatory normalization before sqrtNthOrder(). */
494 0 : expo = add(expo, i);
495 0 : man = sqrt2ndOrder(man);
496 0 : if (s_and(expo, 1) != 0) /* Check even or odd. */
497 : {
498 0 : man = mult_r(man, FEC_HQ_ECU_ROOT2);
499 : }
500 0 : expo = shr_pos(expo, 1); /* Divided by 2-- square root operation. */
501 0 : att_val = shr(man, expo); /* Denormalize the mantissa back to Q15. */
502 : }
503 : /* ELSE
504 : {
505 : do nothing because (gr_pow_right/gr_pow_left) >= 1.0
506 : }
507 : */
508 : }
509 :
510 0 : mag_chg_1st[k] = att_val; move16();
511 0 : mag_chg[k] = att_val; move16();
512 : }
513 : ELSE
514 : {
515 0 : mag_chg_1st[k] = 32767; move16();
516 0 : mag_chg[k] = 32767; move16(); /* Set to ]1.0 in Q15 */
517 : }
518 :
519 0 : L_pGrPowLeft++;
520 0 : L_pGrPowRight++;
521 : } /* FOR band k */
522 : }
523 : ELSE /* sub(burst_len, 1) <= 0) */
524 : {
525 : /* BURST path */
526 :
527 : /* Since attDegreeFrames is discrete (integer) and hard limited to OFF_FRAMES_LIMIT,
528 : * it is easier to implement 10^(-att_degree/20.0) by a simply direct
529 : * table-lookup. Also, att_per_frame is discrete as well and can be
530 : * either ATT_PER_FRAME-1 or ATT_PER_FRAME and nothing else. This
531 : * means only 2 tables of size=(OFF_FRAMES_LIMIT+1) each are required.
532 : * To take square root into account, it is divided by 20 instead of 10. */
533 :
534 0 : if (sub(burst_len, beta_mute_thr) > 0) /* beta_mute_thr coincides/close to stronger 6dB muting phase */
535 : {
536 0 : *beta_mute = shr_pos_pos(*beta_mute, 1);
537 : }
538 :
539 0 : FOR(k = 0; k < Lgw; k++) /* Lgw may be shorter than all bands at 48k */
540 : {
541 : /* global burst attenuation */
542 : #if PLC2_FADEOUT_IN_MS != 0
543 : /* att_per_frame idx = "1:12") */
544 0 : att_val = POW_ATT_TABLES[att_per_frame][s_min(OFF_FRAMES_LIMIT, attDegreeFrames)]; move16();
545 : #else
546 : /* att_per_frame idx = "1:2") */
547 : att_val = POW_ATT_TABLES[att_per_frame][s_min(OFF_FRAMES_LIMIT, attDegreeFrames)]; move16();
548 : /* 10^(-attDegreeFrames*(att_per_frame = "1 or 2")/20) */
549 : #endif
550 :
551 0 : mag_chg[k] = mult_r(mag_chg_1st[k], att_val); /* Q15 */
552 :
553 0 : alpha[k] = mag_chg[k]; move16();
554 0 : ASSERT(beta[k] == 0); /* initialization required */
555 0 : IF(sub(alpha[k], 32766) < 0)
556 : {
557 : /* beta_pre[k] = sqrt(1.0f - SQR(alpha[k])); */
558 : /* beta[k] = beta_pre[k]* *beta_mute;*/
559 : /* (1.0-alpha.^2), in exp 1 due to L_mult0; */
560 :
561 0 : L_acc = L_sub((INT32_MAX >> 1) + 1, L_mult0(alpha[k], alpha[k]));
562 : {
563 :
564 : /* use lower complex(WMOPS/ROM) 2nd-order sqrt approximation */
565 0 : Word32 L_man, L_acc2 = L_acc;
566 : Word16 tmp, expo_in, expo2, man_in, man;
567 : /* updated code using the 2nd order approximation routine */
568 : /* form is flt=(L_acc2*2.^(-31 + 1) */
569 :
570 0 : tmp = norm_l(L_acc2); /* tmp is always 1 or higher due to Lmac0 downshift */
571 0 : man_in = round_fx_sat(L_shl_pos(L_acc2, tmp));
572 0 : expo_in = sub(1, tmp); /* 1 due to original 1 bit margin gain in L_mult0 */
573 :
574 : /* both positive and negative expos into sqrt */
575 0 : man = sqrt2ndOrder(man_in);
576 0 : if (s_and(expo_in, 1) != 0)
577 : {
578 0 : man = mult_r(man, FEC_HQ_ECU_ROOT2); /* odd exp operation */ /* sqrt(2)/2 */
579 : }
580 :
581 0 : expo2 = shr_r(expo_in, 1); /* apply square root operation. shr_r needed for pos and neg exps */
582 0 : ASSERT(expo2 <= 1);
583 :
584 0 : L_man = L_deposit_h(man);
585 0 : L_man = L_shl_sat(L_man, expo2); /* move to a zero exp , _sat needed for 1.0 input (due to approximation overshoot) */
586 :
587 0 : man = round_fx_sat(L_man); /* better perf with round here */
588 :
589 :
590 0 : beta[k] = mult_r(*beta_mute, man); move16();
591 : }
592 :
593 : /* bw Lowpass shape additive component */
594 : /* tab[LGW48K + 1] = { 1.0, ....1.0, 0.5,0.5, ... 0.1, 0.1 } */
595 :
596 0 : IF(sub(k, LGW32K - 1) >= 0)
597 : {
598 0 : beta[k] = mult_r(beta[k], 3277); /* 0.1 in Q15 */
599 : }
600 0 : ELSE IF(sub(k, LGW16K - 1) >= 0)
601 : {
602 0 : beta[k] = mult_r(beta[k], 16384); /* 0.5 in Q15 */
603 : }
604 :
605 : /*
606 : % limit Xavg noise contribution further in case of offset/tr_decay
607 : % attenuation was already active
608 :
609 : if (burst_len <= burst_att_thresh) && ( stPhECU_mag_chg_1st(k) < (32766/32768) )
610 : XavgFadeinFactor = (burst_len-1)/burst_att_thresh;
611 : XavgFadeinFactor = min(1.0, XavgFadeinFactor);
612 : beta(k) = beta(k)*XavgFadeinFactor;
613 : % limit initial Xavg noise contribution until We have reached regular burst attenuation
614 : end
615 : end
616 : */
617 0 : IF( sub(mag_chg_1st[k], 32767) <0 )
618 : { /* offset muting was started before burst muting phase */
619 : /* Xavg noise gradually increased during a short period */
620 0 : Word16 XavgFadeinFactor = 32767;
621 0 : Word16 ratio2_3_4_5tab[4][5 - 1] = {
622 : {(Word16)(.5 * 32768.0), (Word16)(1.0 * 32767.0), (Word16)(1.0 * 32767.0), (Word16)(1.0 * 32767.0)}, /* 1/2*/
623 : {(Word16)(.333 * 32768.0), (Word16)(.666 * 32768.0), (Word16)(1.0 * 32767.0), (Word16)(1.0 * 32767.0)}, /* 1/3*/
624 : {(Word16)(.25 * 32768.0), (Word16)(.5 * 32768.0), (Word16)(.75 * 32768.0), (Word16)(1.0 * 32767.0)}, /* 1/4 */
625 : {(Word16)(.2 * 32768.0), (Word16)(.4 * 32768.0), (Word16)(.6 * 32768.0), (Word16)(.8 * 32768.0)} /* 1/5 */
626 :
627 : };
628 0 : ASSERT(burst_att_thresh >= 1 && burst_att_thresh <= 5);
629 0 : ASSERT(burst_len >= 2);
630 0 : if (sub(burst_len,burst_att_thresh) <= 0)
631 : {
632 0 : ASSERT(burst_len - 2 < (5-1));
633 0 : ASSERT(burst_att_thresh-1-1 < (4));
634 0 : XavgFadeinFactor = ratio2_3_4_5tab[burst_att_thresh-1-1][burst_len - 2]; /* second bfi frame burst_len= is 2 */
635 : }
636 0 : beta[k] = mult_r(beta[k], XavgFadeinFactor); /* n Q15 */
637 : }
638 : } /* IF (sub(alpha[k], 32766) < 0) */
639 : } /* FOR k*/
640 :
641 : } /* BURST */
642 :
643 :
644 :
645 0 : IF(sub(output_frame, L_FRAME48K) == 0)
646 : { /* for 48kHz set/handle scalings of last group/band the same way as previous lower freq band(s) */
647 :
648 0 : FOR(k = Lgw; k < MAX_LGW; k++)
649 : {
650 0 : tr_dec[k] = tr_dec[k - 1]; move16(); /* only available in first bfi frame */
651 0 : Xavg[k] = Xavg[k - 1]; move16();
652 0 : mag_chg_1st[k] = mag_chg_1st[k - 1]; move16();
653 0 : mag_chg[k] = mag_chg[k - 1]; move16();
654 0 : alpha[k] = alpha[k - 1]; move16();
655 0 : beta[k] = beta[k - 1]; move16();
656 : }
657 : }
658 :
659 : #ifdef DYNMEM_COUNT
660 : Dyn_Mem_Out();
661 : #endif
662 : #ifdef WMOPS
663 : pop_wmops();
664 : #endif
665 0 : }
666 :
667 : /*-----------------------------------------------------------------------------
668 : * imax_fx()
669 : *
670 : * Get interpolated maximum position
671 : *-----------------------------------------------------------------------------*/
672 0 : static Word16 imax_fx( /* o: The location, relative to the middle of the 3 given data point, of the maximum. (Q15) */
673 : const Word16 *y, /* i: The 3 given data points. */
674 : const Word16 special /* i: -1 = left edge special case, 0 = normal, +1 = right edge special case */
675 : )
676 : {
677 : Word16 posi;
678 : Word16 man, expo, edge;
679 : const Word16 *pY;
680 : Word32 L_y1, L_y2, L_y3, L_numer, L_denom, L_sign, L_acc, L_y3_y1;
681 :
682 : #ifdef DYNMEM_COUNT
683 : Dyn_Mem_In("imax_fx", sizeof(struct {
684 : Word16 posi;
685 : Word16 man, expo, edge;
686 : const Word16 *pY;
687 : Word32 L_y1, L_y2, L_y3, L_numer, L_denom, L_sign, L_acc, L_y3_y1;
688 : }));
689 : #endif
690 :
691 : #ifdef WMOPS
692 : push_wmops("PhECU::imax_fx");
693 : #endif
694 :
695 : /* Seek the extremum of the parabola P(x) defined by 3 consecutive points
696 : so that P([-1 0 1]) = [y1 y2 y3] */
697 0 : pY = y;
698 0 : L_y1 = L_deposit_l(*pY++), L_y2 = L_deposit_l(*pY++), L_y3 = L_deposit_l(*pY);
699 :
700 : /* The extremum value:
701 : * y2i = -0.125f * SQR(y3_y1) / (y1+y3-2*y2)+y2
702 : * is not computed. Alternatively, the derivative of the parabola evaluated at y=0,
703 : * dP/dy|y=0, is used to determine whether the extremum is maximum or not.
704 : */
705 :
706 : /* Compute the extremum location: posi = (y3 - y1)/(4*y2 - 2*y1 - 2*y3). */
707 0 : L_y3_y1 = L_sub(L_y3, L_y1);
708 0 : L_acc = L_shl_pos(L_y2, 1); /* N.B. y2 is multiplied by 2 not 4. */
709 0 : L_acc = L_sub(L_acc, L_y1); /* N.B. Y1 is not multiplied by 2. */
710 0 : L_denom = L_sub(L_acc, L_y3); /* N.B. Y3 is not multiplied by 2. */
711 0 : L_sign = L_xor(L_y3_y1, L_denom); /* Preserve the sign since div_s() only takes positive arguments. */
712 0 : L_numer = L_abs(L_y3_y1);
713 0 : L_denom = L_abs(L_denom);
714 :
715 0 : test();
716 0 : IF(L_numer == 0 || L_denom == 0)
717 : {
718 0 : posi = 0; move16(); /* flat top , exit with center freq. */
719 : }
720 : ELSE
721 : {
722 :
723 0 : IF(L_sub(L_denom, L_shr_pos_pos(L_numer, 1)) > 0)
724 : {
725 : /* Although the output of ratio() is in Q14, adding the missing factor of 2 (See above)
726 : * in the denominator, the output is now considered to be in Q15. */
727 0 : man = plc_phEcu_ratio_fx(L_numer, L_denom, &expo); /* The mantissa is considered in Q15 */
728 :
729 0 : posi = shr_sat(man, expo); /* in Q15 (Due to saturation, it is automatically bound inside [-1.0,1.0].) */
730 : }
731 : ELSE
732 : {
733 0 : posi = 0x7fff; move16();
734 : }
735 :
736 0 : if (L_sign < 0) /* Restore the sign. */
737 : {
738 0 : posi = negate(posi);
739 : }
740 :
741 : /* For both edges (left and right), the extremum found above may be minimum.
742 : * It needs to reject the minimum. */
743 0 : IF(special != 0) /* Either edge special case. */
744 : {
745 0 : edge = 0x7fff; /* 1 in Q15 for the right edge special case */ move16();
746 0 : if (special < 0)
747 : {
748 0 : edge = 0; /* Left edge special case */ move16();
749 : }
750 :
751 : /* The derivative (slope) of the interpolating parabola = 2*A*y + B,
752 : * where A = (y3 + y1)/2 - y2
753 : * and B = (y3 - y1)/2.
754 : * Therefore, the slope at y=0 is simply B. Use this slope to determine
755 : * if the parabola is concave upward or downward.
756 : */
757 0 : IF(posi > 0) /* The extremum is in between the middle and the right given data points. */
758 : {
759 0 : posi = sub(0x7fff, posi); /* maximum case */
760 0 : if (L_sub(L_y3, L_y1) <= 0) /* Check the slope at y=0, i.e., at the middle given data point. */
761 : {
762 0 : posi = edge; /* minimum case */ move16();
763 : }
764 : }
765 : ELSE /* The extremum is in between the left and the middle given data points. */
766 : {
767 0 : posi = add(0x7fff, posi); /* maximum case */
768 0 : if (L_sub(L_y3, L_y1) >= 0)
769 : {
770 0 : posi = edge; /* minimum case */ move16();
771 : }
772 : }
773 : }
774 : }
775 : #ifdef DYNMEM_COUNT
776 : Dyn_Mem_Out();
777 : #endif
778 : #ifdef WMOPS
779 : pop_wmops();
780 : #endif
781 0 : return posi;
782 : /* Q15. The position either left or right relative to the index of the middle of the 3 given data points. */
783 : }
784 :
785 : /*-----------------------------------------------------------------------------
786 : * spec_ana_fx()
787 : *
788 : * Spectral analysis
789 : *-----------------------------------------------------------------------------*/
790 : /* OPT add the FB transient input flags , and skip peakfinder if fullband transient is set */
791 0 : void spec_ana_fx(Word16 * xfp, /* i/o : Input 16ms pre-upscaled time signal, output xfp utility buffer */
792 : Word16 * plocs, /* o : The indicies of the identified peaks Q0 */
793 : Word32 * L_plocsi, /* o : Interpolated positions of the identified peaks Q16 */
794 : Word16 * num_plocs, /* o : Number of identified peaks Q0 */
795 : Word16 * X_sav, /* o : Stored fft spectrum */
796 : const Word16 output_frame, /* i : Frame length Q0 */
797 : const Word16 bwidth_fx, /* i : Encoded Fs index Q0 */
798 : const Word16 *sp_ana_win, /* i : spectral analysis window Q15 */
799 : const Word16 f0hzLtpBinQ7, /* i : LTP bin frequency in normalized Hz Q7 */
800 : const Word16 norm_corrQ15_fx, /* i : correlation for lag at f0hzLtpBinQ7 */
801 : Word16 maxLprot, Word16 maxPlocs,
802 : Word8 *scratchBuffer /* Size = 4 * (MAX_LPROT + MAX_LPROT_RED + 1) + 2 * MAX_PLOCS */
803 : )
804 : {
805 : Counter n, k;
806 0 : Word16 nJacob, Lprot, hamm_len2 = 0, Lprot2, Lprot2p1;
807 : Word32 *L_xfp;
808 :
809 : Word16 *pXfp;
810 : Word16 *pPlocs;
811 : Word16 Xmax, Xmin, sens;
812 : Word16 rectLength, fraction;
813 : Word32 *pPlocsi_L;
814 : Word32 L_acc;
815 : Word16 peak_range_1;
816 : Word16 stop_band_start;
817 : Word16 stop_band_length;
818 : Word16 fft_scale;
819 : Word8 * buffer_fft;
820 : Word16 currPlocs, endPlocs;
821 : Word16 P_in_plocs;
822 : Word16 n_real_interp_tail;
823 :
824 : #ifdef WMOPS
825 : push_wmops("PhECU::spec_ana_fx(1st)");
826 : #endif
827 : #ifdef DYNMEM_COUNT
828 : Dyn_Mem_In("spec_ana_fx", sizeof(struct {
829 : Counter n, k;
830 : Word16 nJacob ,Lprot, hamm_len2, Lprot2, Lprot2p1;
831 : Word32 *L_xfp;
832 : Word32 *pXfp_L;
833 : Word16 *y_re_ptr, *y_im_ptr; /* otrs to Xsav as xfp was overwritten */
834 : Word16 *pXfp, *pXfp1, *pPlocs;
835 : Word16 Xmax, Xmin, sel;
836 : Word16 rectLength, fraction, special;
837 : Word32 *pPlocsi_L;
838 : Word32 L_acc;
839 : Word16 peak_range_1;
840 : Word16 stop_band_start;
841 : Word16 stop_band_length;
842 : Word16 fft_scale;
843 : Word8 * buffer_fft;
844 : Word16 fft_scale_by4;
845 : Word16 currPlocs, endPlocs;
846 : Word16 P_in_plocs;
847 : Word16 n_real_interp_tail;
848 : }));
849 : #endif
850 : /* Initialize for 48k to avoid warnings
851 : Lprot - length of saved prototype samples
852 : hamm_len2 - half Hamming window length
853 : pFftTbl - Table for real input FFT
854 : LprotLog2Minus1 - FFT stages for complex input FFT
855 : */
856 :
857 :
858 0 : L_xfp = (Word32 *)scratchAlign(scratchBuffer, 0); /* Size = 4 * MAX_LPROT bytes */
859 0 : buffer_fft = scratchAlign(L_xfp, sizeof(*L_xfp) * maxLprot); /* Size = 4 * (MAX_LPROT_RED + 1) + 2 * MAX_PLOCS */
860 :
861 0 : ASSERT(bwidth_fx >= 0 && bwidth_fx <= 4); /* avoid bwidth_fx variable warning */
862 :
863 0 : Lprot = LprotSzPtr[bwidth_fx];
864 0 : move16();
865 0 : hamm_len2 = DEPR_i_mult(3, mult(output_frame, (Word16)(3277) /* divBy10 + floor */)); /* 3 ms */
866 0 : fft_scale = PhEcuFftScale[bwidth_fx];
867 0 : move16(); /* 32,16,8 kHz all have fft scale 0, 24 has 8, 48 has 4 */
868 0 : Lprot2 = shr_pos(Lprot, 1);
869 0 : Lprot2p1 = add(Lprot2, 1); /* Magnitude lengths */
870 0 : rectLength = sub(Lprot, shl_pos(hamm_len2, 1));
871 : /* The length of the rectangular portion of the Hamming-Rectangular window. */
872 : {
873 : #ifdef WMOPS
874 : push_wmops("PhECU::WhrAnaWin+fft");
875 : #endif
876 :
877 : /* Apply hamming-rect window */
878 0 : windowing_L(xfp, L_xfp, sp_ana_win, rectLength, hamm_len2);
879 0 : BASOP_rfftN(L_xfp, Lprot, &fft_scale, buffer_fft);
880 : }
881 : #ifdef WMOPS
882 : pop_wmops(); /* anawin+fft */
883 : #endif
884 :
885 : /* Convert 32 Bit intlv FFT into phecu 16 bit flipped fft format */
886 : /* can not yet be an inplace operation */
887 :
888 0 : intlvW32_2_flippedW16(L_xfp, sub(Lprot2, 1), Lprot, xfp);
889 :
890 : /* Apply zeroing of non-coded FFT spectrum above 20 kHz */
891 0 : IF(sub(output_frame, ((L_FRAME48K) * 40) / 48) >= 0) /* only relevant for 48kHz in LC3plus */
892 : {
893 0 : stop_band_start = ((LPROT48K / 2) * 40) / 48; /* initial start position in real part , 320 */
894 0 : stop_band_length = ((LPROT48K * 8) / 48) - 1; /* real tail and into Im parts , 128-1 */
895 0 : stop_band_start = add(stop_band_start, 1); /* exclude DC ... */
896 :
897 0 : basop_memset(xfp + stop_band_start, 0, (stop_band_length) * sizeof(Word16));
898 : }
899 :
900 0 : peak_range_1 = s_min(Lprot2p1, MAX_LPROT_RED / 2 + 1); /* limit preliminary only active for 48k to save WMOPS */
901 :
902 0 : basop_memmove(X_sav, xfp, (Lprot) * sizeof(Word16));
903 :
904 : /* Magnitude representation */
905 0 : fft_spec2_sqrt_approx_fx(xfp, xfp, Lprot);
906 : /* inplace, i.e. [ Dc, real part of xfp ,Fs/2 ] will be replaced by magnitude(scaled by .5) */
907 :
908 : /* Find global maximum and minimum. */
909 0 : plc_phEcu_maxval_fx(xfp, peak_range_1, &Xmax);
910 0 : plc_phEcu_minval_fx(xfp, peak_range_1, &Xmin);
911 0 : sens = mult_r(sub(Xmax, Xmin), CMPLMNT_PLOC_SENS_FX);
912 :
913 :
914 0 : plc_phEcu_peak_locator_fx(xfp, peak_range_1, plocs, num_plocs, sens, Xmax, Xmin, MAX_LPROT_RED, buffer_fft);
915 :
916 :
917 : #ifdef WMOPS
918 : push_wmops("PhECU::Peaks_refine");
919 : #endif
920 :
921 : /* Refine peaks */
922 0 : pPlocsi_L = L_plocsi;
923 0 : pPlocs = plocs;
924 : /* n = sub(*num_plocs, 1); */ /* -1 so as to exclude the very last peak. */
925 0 : n = *num_plocs; /* number of peaks to process */
926 : /* Special case-- The very 1st peak if it is at 0 index position (DC) */
927 : /* With DELTA_CORR_F0_INT == 2 one needs to handle both *pPlocs==0 and *pPlocs==1 */
928 0 : logic16();
929 0 : IF((n > 0) && (sub(*pPlocs, 0) == 0)) /* Very 1st peak position possible to have a peak at 0/DC index position. */
930 : {
931 0 : fraction = imax_fx(xfp, -1); /* -1 signifies special left edge case. */
932 0 : L_acc = L_deposit_h(*pPlocs++); /* N.B., (*pPlocs) must be zero here. */
933 0 : *pPlocsi_L++ = L_mac(L_acc, fraction, 1); move32(); /* in Q16 */
934 0 : n = sub(n, 1); /* This special case is taken care of -- one less peak to go */
935 : }
936 0 : logic16();
937 0 : IF((n > 0) && (sub(*pPlocs, 1) == 0)) /* Also 2nd peak position uses DC which makes jacobsen unsuitable. */
938 : {
939 0 : fraction = imax_fx(xfp, 0); /* for parabolic this is not a special case. */
940 0 : L_acc = L_deposit_h(*pPlocs++); /* N.B., (*pPlocs) must be 1 here. */
941 0 : *pPlocsi_L++ = L_mac(L_acc, fraction, 1); move32(); /* in Q16 */
942 0 : n = sub(n, 1); /* This special case is taken care of -- one less peak to go */
943 : }
944 :
945 : /* All remaining peaks except the very last two possible integer positions */
946 0 : currPlocs = *pPlocs++; move16();
947 0 : endPlocs = sub(Lprot2p1, DELTA_CORR_F0_INT); /* last *pPlocs position for Jacobsen */
948 :
949 : /* precompute number of turns based on endpoint integer location and make into a proper FOR loop */
950 0 : IF(n > 0)
951 : {
952 0 : nJacob = n; move16();
953 :
954 : /* catch all three xxx01 , xxx10 and xxx11 ,
955 : and not only xxx01, xxx10 */
956 :
957 0 : n_real_interp_tail = 0; move16();
958 0 : if (sub(endPlocs, plocs[sub(*num_plocs, 1)]) <= 0)
959 : {
960 0 : n_real_interp_tail = add(n_real_interp_tail, 1);
961 : }
962 :
963 0 : logic16();
964 0 : if (sub(n,1) > 0 && sub(endPlocs, plocs[sub(*num_plocs, 2)]) <= 0)
965 : {
966 0 : n_real_interp_tail = add(n_real_interp_tail, 1);
967 : }
968 :
969 0 : nJacob = sub(nJacob, n_real_interp_tail);
970 :
971 0 : FOR (k = 0; k < nJacob; k++)
972 : {
973 0 : fraction = imax2_jacobsen_mag_fx(&(X_sav[currPlocs - 1]), &(X_sav[Lprot - 1 - currPlocs]), 0); /* in Q15 */ /* not endpoint */
974 0 : move16(); move16(); /* account for inloop indirect ptrs into Xsav */
975 :
976 0 : L_acc = L_deposit_h(currPlocs);
977 0 : *pPlocsi_L++ = L_mac(L_acc, fraction, 1); move32(); /* in Q16. Append the fractional part to the integral part. */
978 0 : currPlocs = *pPlocs++; move16();
979 : }
980 0 : n = sub(n, nJacob);
981 : }
982 :
983 : /* At this point there should at most two plocs left to process */
984 : /* the position before fs/2 and fs/2 both use the same magnitude points */
985 :
986 0 : IF(n > 0)
987 : {
988 : /* [ . . . . . . . ] Lprot/2+1 positions */
989 : /* | | | */
990 : /* 0 (Lprot/2-2) (Lprot/2) */
991 :
992 0 : pXfp = xfp + sub(Lprot2, 2);
993 0 : IF(sub(currPlocs, sub(Lprot2p1, DELTA_CORR_F0_INT)) == 0)
994 : /* Also 2nd last peak position uses fs/2 which makes jacobsen less suitable. */
995 : {
996 0 : fraction = imax_fx(pXfp, 0); /* for parabolic this is not a special case. */
997 :
998 0 : L_acc = L_deposit_h(currPlocs); /* N.B., (*pPlocs) must be 1 here. */
999 0 : *pPlocsi_L++ = L_mac(L_acc, fraction, 1); move32(); /* in Q16 */
1000 0 : currPlocs = *pPlocs++; move16();
1001 0 : n = sub(n, 1); /* This special case is taken care of -- one less peak to go */
1002 :
1003 0 : if (n > 0)
1004 : {
1005 : /* allow for an additional consecutive final plocs after (Fs/2-1) at Fs/2 */
1006 0 : currPlocs = *pPlocs++; move16();
1007 : }
1008 : }
1009 :
1010 : /* Here the only remaining point would be a fs/2 plocs */
1011 : /* pXfp = xfp + sub(Lprot2,1); already set just a reminder where it whould point */
1012 0 : IF(n > 0) /* fs/2 which makes special case . */
1013 : {
1014 0 : fraction = imax_fx(pXfp, 1); /* for parabolic this is a special case. */
1015 :
1016 0 : L_acc = L_deposit_h(currPlocs); /* N.B., (*pPlocs) must be 1 here. */
1017 0 : *pPlocsi_L++ = L_mac(L_acc, fraction, 1); move32(); /* in Q16 */
1018 0 : currPlocs = *pPlocs++; move16();
1019 0 : n = sub(n, 1); /* This special case is taken care of -- one less peak to go */
1020 : }
1021 : }
1022 :
1023 : /* here n should be 0 if all peaks have been processed */
1024 0 : ASSERT(n == 0);
1025 :
1026 : /* Check number of plocs within an assumed pitch range */
1027 0 : P_in_plocs = 0; move16();
1028 0 : FOR(n = 0; n < *num_plocs; n++)
1029 : {
1030 : /* count number of peaks in locations 1,2,3,4,5,6 , ~= 60 Hz ... 380 Hz */
1031 0 : fraction = s_min(1, plocs[n]); /* 0 stays zero , otherwise 1 */
1032 0 : if (sub(plocs[n], 7) < 0)
1033 : {
1034 0 : P_in_plocs = add(P_in_plocs, fraction);
1035 : }
1036 : }
1037 :
1038 : #ifdef WMOPS
1039 : pop_wmops(); /* peaks refine */
1040 : #endif
1041 :
1042 0 : logic16();
1043 0 : IF(f0hzLtpBinQ7 > 0 && P_in_plocs > 0)
1044 : {
1045 : Word16 n_plocs_in, n_plocs_out;
1046 :
1047 0 : n_plocs_in = *num_plocs; move16();
1048 :
1049 : /* NB LF peak analysis may add adjacent peaks in { plocs, L_plocsi}, (output from peakfinder did not have
1050 : * adjacent peaks ) */
1051 0 : plc_phEcu_LF_peak_analysis_fx(plocs /* i/o */, num_plocs /* i/o */, L_plocsi /* i/o */, xfp, f0hzLtpBinQ7,
1052 : norm_corrQ15_fx, 2, maxPlocs, buffer_fft);
1053 0 : n_plocs_out = *num_plocs; move16();
1054 :
1055 0 : IF(sub(n_plocs_in, n_plocs_out) == 0)
1056 : {
1057 : /* adjust first peak coinciding with LTPF0 measures if it indicates so */
1058 0 : plc_phEcu_F0_refine_first_fx(plocs /* i/o */, *num_plocs, L_plocsi /* i/o */, f0hzLtpBinQ7, norm_corrQ15_fx,
1059 : 3);
1060 : }
1061 : }
1062 : /* moved inside spec_ana_fx , to include validated pitch peak P_in_plocs*/
1063 : {
1064 0 : Word16 peak_limits_fx[5] = { 14 /*NB*/, 14 /*WB*/, 14 /*sWB*/, 14 /*SWB*/,
1065 : 14 /* FB */ }; /* to be trained for each BW */
1066 :
1067 0 : test();
1068 0 : test();
1069 0 : IF((sub(norm_corrQ15_fx, 0) > 0) && /* == 0 indicates a negative correlation, which could be likely stable */
1070 : (sub((Word16)(0.5 * 32768), norm_corrQ15_fx) > 0) && (sub(*num_plocs, peak_limits_fx[bwidth_fx]) > 0))
1071 : {
1072 0 : if (P_in_plocs > 0)
1073 : {
1074 0 : *num_plocs = 0; move16(); /*activate noise only path only if normcorr vas valid energywise */
1075 : }
1076 : }
1077 : }
1078 :
1079 : #ifdef DYNMEM_COUNT
1080 : Dyn_Mem_Out();
1081 : #endif
1082 : #ifdef WMOPS
1083 : pop_wmops();
1084 : #endif
1085 0 : }
1086 :
1087 : /*-------------------------------------------------------------------*
1088 : * subst_spec_fx()
1089 : *
1090 : * Substitution spectrum calculation
1091 : *-------------------------------------------------------------------*/
1092 :
1093 0 : void subst_spec_fx(
1094 : const Word16 *plocs, /* i : The indices of the identified peaks Q0 */
1095 : const Word32 *L_plocsi, /* i : Interpolated positions of the identified peaks Q16 */
1096 : Word16 * num_plocs, /* i/o : Number of identified peaks Q0 */
1097 : const Word16 time_offs, /* i : Time offset Q0 */
1098 : Word16 * X, /* i/o : FFT spectrum */
1099 : const Word16 *mag_chg, /* i : Magnitude modification Q15 */
1100 : const Word16 ph_dith, /* i : Phase dither, 2*PI is not included. (Q15, i.e., between 0.0 and 1.0) */
1101 : const Word16 *is_trans, /* i : (Transient) noise generation flags (either 0 or not 1 ) */
1102 : const Word16 output_frame, /* i : Frame length Q0 */
1103 : Word16 * seed, /* i/o : Random seed */
1104 : const Word16 *alpha, /* i : Magnitude modification factors for fade to average Q15 */
1105 : const Word16 *beta, /* i : Magnitude modification factors for fade to average Q15 */
1106 : const Word16 *Xavg, /* i : Frequency group averages to fade to Q0 */
1107 : const Word16 t_adv /* i : time adjustement excluding time_offs Q0 */
1108 : ,const Word16 fadeout, /* need for DC muting */
1109 : Word16 * nonpure_tone_flag_ptr, /* i/o : non-pure single tone indicator state */
1110 : const Word32 * L_Xavg /* i : Frequency group amp averages for tonal tilt analysis Max upshifted */
1111 : )
1112 : {
1113 : Word16 Xph_short;
1114 : Word32 L_corr_phase[MAX_PLOCS], L_Xph;
1115 : Word32 *pCorrPhase_L;
1116 : Word16 cos_F, sin_F, tmp;
1117 0 : Word16 peak_sin_F = 0, peak_cos_F = 0;
1118 : Word16 sin_F_fade2avg, cos_F_fade2avg;
1119 : Word16 fs_idx;
1120 : Word16 Lprot, m, i, e, im_ind, delta_corr_up, delta_corr_dn, delta_tmp;
1121 : UWord16 lsb;
1122 : Word16 j, re, im, *pReX, *pImX, lastPeak, lprotBy2Minus1, segmentLen;
1123 : Word16 pkLocation_1, pkLocation, pkLocation1;
1124 : const Word16 *pPlocs;
1125 : const Word32 *pPlocsi_L;
1126 : Word32 L_acc;
1127 : Word16 Lprot_inv;
1128 : Word16 k;
1129 : Word16 tmp2;
1130 : Word16 alpha_local;
1131 : Word32 tmp_L;
1132 : Word16 mag_chg_local;
1133 : const Word16 *gwlpr_fxPlus1;
1134 : Word16 one_peak_flag_mask;
1135 : Word16 noise_mag_scale_neg;
1136 : Word16 up_shift_adj;
1137 :
1138 : #ifdef DYNMEM_COUNT
1139 : Dyn_Mem_In("subst_spec_fx", sizeof(struct {
1140 : Word16 Xph_short;
1141 : Word32 L_corr_phase[MAX_PLOCS], L_Xph;
1142 : Word32 *pCorrPhase_L;
1143 : Word16 cos_F, sin_F, tmp;
1144 : Word16 peak_sin_F, peak_cos_F;
1145 : Word16 sin_F_fade2avg, cos_F_fade2avg;
1146 : Word16 fs_idx;
1147 : Word16 Lprot, m, i, e, im_ind, delta_corr_up, delta_corr_dn, delta_tmp;
1148 : UWord16 lsb;
1149 : Word16 j, re, im, *pReX, *pImX, lastPeak, lprotBy2Minus1, segmentLen;
1150 : Word16 pkLocation_1, pkLocation, pkLocation1;
1151 : const Word16 *pPlocs;
1152 : const Word32 *pPlocsi_L;
1153 : Word32 L_acc;
1154 : Word16 Lprot_inv;
1155 : Word16 k;
1156 : Word16 tmp2;
1157 : Word16 alpha_local;
1158 : Word16 expo;
1159 : Word32 tmp_L;
1160 : Word16 mag_chg_local;
1161 : const Word16 *gwlpr_fxPlus1;
1162 : Word16 one_peak_flag_mask;
1163 : Word16 noise_mag_scale_neg;
1164 : Word16 up_shift_adj;
1165 : }));
1166 : #endif
1167 :
1168 : if (time_offs == 0)
1169 : {
1170 : #ifdef WMOPS
1171 : push_wmops("PhECU::subst_spec_fx(1st)");
1172 : #endif
1173 : }
1174 : else
1175 : {
1176 : #ifdef WMOPS
1177 : push_wmops("PhECU::subst_spec_fx(N)");
1178 : #endif
1179 : }
1180 :
1181 :
1182 0 : gwlpr_fxPlus1 = &(gwlpr_fx[1]); /* ptr init */
1183 0 : fs_idx = mult(output_frame, (Word16)(32768.0 / 99.0)); /* truncation needed , i.e no rounding can be applied here */
1184 0 : ASSERT(fs_idx == (output_frame / 100));
1185 0 : Lprot = LprotSzPtr[fs_idx];move16();
1186 0 : Lprot_inv = InvLprot_Q22[fs_idx]; move16();
1187 :
1188 0 : tmp2 = add(mult_r(time_offs, oneOverFrameQ15Tab[fs_idx]), 1);/* save a local burst_len for securing DC and fs/2 muting */
1189 :
1190 : /* Correction/evolution phase of the identified peaks */
1191 0 : IF(s_or(is_trans[0], is_trans[1]) != 0)
1192 : {
1193 0 : *num_plocs = 0; move16();
1194 : }
1195 : ELSE
1196 : {
1197 :
1198 0 : tmp = t_adv;
1199 :
1200 0 : tmp = add_sat(tmp, time_offs);
1201 : /* NB tmp can be stored in Word16 Q0 as max used value is 684+(OFF FRAMELIMIT==60)*480 = 29484 */
1202 0 : tmp_L = L_mult0(tmp, Lprot_inv);
1203 :
1204 :
1205 0 : tmp = norm_l(tmp_L);
1206 0 : up_shift_adj = s_max(0, sub(4, tmp)); /* 48kHz : PLC frame 1..49 -> tmp_L<<4 , frame 50..60 -> tmpL<<3 */
1207 0 : tmp_L = L_shl_sat(tmp_L, sub(4, up_shift_adj));
1208 0 : tmp = round_fx_sat( tmp_L);
1209 :
1210 :
1211 0 : pPlocsi_L = L_plocsi;
1212 0 : pCorrPhase_L = L_corr_phase;
1213 0 : FOR(m = 0; m < *num_plocs; m++)
1214 : {
1215 :
1216 : /* tmp has variable resolution 10 or 9 bits */
1217 0 : ASSERT( up_shift_adj >= 0);
1218 0 : Mpy_32_16_ss(L_shl_pos(*pPlocsi_L++, up_shift_adj), tmp, &L_acc, &lsb);
1219 0 : L_acc = L_add(L_shl_pos(L_acc, 5), lshr(lsb, 11));
1220 : /* 5 lsb's actually unused though, as 6 bits are shifted out */
1221 0 : *pCorrPhase_L++ = L_acc; move32();/* in Q16. 2*PI is not included. */
1222 : }
1223 : }
1224 :
1225 0 : one_peak_flag_mask = (Word16)0xFFFF; move16(); /* all ones mask -> keep */
1226 0 : test(); logic16();
1227 0 : IF((*num_plocs > 0) && sub(*num_plocs, 3) < 0)
1228 : {
1229 0 : one_peak_flag_mask = 0x0000; move16(); /* all zeroes mask -> zeroes in valleys, single clean tone assumed */
1230 :
1231 : /* revert initial pure tone decision in some cases */
1232 0 : logic16(); logic16();
1233 0 : IF((sub(*nonpure_tone_flag_ptr, 0) < 0) &&
1234 : ((sub(fs_idx, 2) == 0)/* SemiSWB 24 kHz */ || (sub(fs_idx, 4) >= 0)) /* FB 48 kHz */
1235 : )
1236 : {
1237 : /* in the first lost frame analyze spectra and spectral bands to possibly reverse an initial pure sine assumption */
1238 0 : *nonpure_tone_flag_ptr = plc_phEcu_nonpure_tone_ana_fx(plocs, *num_plocs, X, L_Xavg, Lprot, fs_idx);
1239 :
1240 : #ifdef LOCAL_PLC2_TON_ANA_DEACTIVATE
1241 : *nonpure_tone_flag_ptr = 0; /* dbg of inactive tone analysis */
1242 : #endif
1243 : }
1244 :
1245 :
1246 0 : if (sub(*nonpure_tone_flag_ptr, 0) > 0)
1247 : {
1248 : /* actually revert single pure tone detection */ /* 0-> mute all surrounding valley bins in evolution , 0xff -> generate noise in all valleys */
1249 0 : one_peak_flag_mask = (Word16)0xFFFF; move16(); /* all ones mask -> keep */
1250 : }
1251 : }
1252 :
1253 0 : noise_mag_scale_neg = 0; move16(); /* no change of valley noise magnitude */
1254 0 : logic16();
1255 0 : if ((*num_plocs == 0) || (time_offs != 0))
1256 : { /* only adj_scale noise amplitude when we have no WC path , or in all noise frame */
1257 0 : noise_mag_scale_neg = -32768; move16(); /* all ones --> scale_noise */
1258 : }
1259 :
1260 0 : IF(*num_plocs == 0)
1261 : {
1262 0 : X[0] = 0; move16(); /* reset DC if there are no peaks */
1263 0 : X[shr_pos(Lprot, 1)] = 0; move16(); /* also reset fs/2 if there are no peaks */
1264 : }
1265 :
1266 : /* the binary selection of fadeout scheme */
1267 0 : tmp = (PLC2_FADEOUT_IN_MS - PLC2_FADEOUT_IN_MS_MIN) / PLC2_FADEOUT_RES; move16();
1268 0 : if (fadeout != 0)
1269 : {
1270 0 : tmp = (PLC2_FADEOUT_LONG_IN_MS - PLC2_FADEOUT_IN_MS_MIN) / PLC2_FADEOUT_RES; move16();
1271 : }
1272 :
1273 0 : IF(sub(tmp2, add(fade_scheme_tab_fx[tmp][1],1) ) > 0)
1274 : {
1275 : /* also start DC scaling attenuation */
1276 0 : X[0] = mult(alpha[0], X[0]); move16();
1277 : /* start fs/by2 attenuation */
1278 0 : X[shr_pos(Lprot, 1)] = mult(alpha[s_min(add(fs_idx, LGW8K), LGW48K)], X[shr_pos(Lprot, 1)]); move16();
1279 : }
1280 :
1281 :
1282 0 : lprotBy2Minus1 = sub(shr_pos(Lprot, 1), 1);
1283 :
1284 : /* for 48k the last (63+1+63) 128/2-1 values above 20 kHz are always zeroes */
1285 : { /* after a last, peak no need to evolve above 20kHz , those coeffs have been/ will be zeroed already */
1286 : /* DC, Xused,ReZeroed, fs/2, ImZeroed, ImUsed */
1287 : /*N: 1, 320 , 63 , 1 , 63 , 320 */
1288 : /*Cind: 0, 1-320, 321-383 , 384 , 385-447, 448 -767 */
1289 0 : lprotBy2Minus1 = s_min( lprotBy2Minus1, 320);
1290 : /*only process up to bin X[320] , 63 coeffs X[321] fwd should already be zeroed */
1291 : }
1292 :
1293 0 : i = 1; move16(); /* index in the X[DC, ReX part of X_Sav ,i.e *pReX == X[i] */
1294 0 : k = 0; move16();
1295 :
1296 0 : pReX = X + i;
1297 0 : im_ind = (Lprot - 1); /* ptr init */
1298 0 : pImX = X + im_ind;
1299 0 : pPlocs = plocs;
1300 0 : pCorrPhase_L = L_corr_phase;
1301 :
1302 :
1303 0 : pkLocation_1 = -4; move16(); /* dummy value to avoid Out of Array Read */
1304 0 : pkLocation = -3; move16(); /* dummy value to avoid Out of Array Read for the *num_plocs==0 case */
1305 0 : pkLocation1 = -2; move16(); /* dummy value to avoid Out of Array Read */
1306 :
1307 0 : IF (*num_plocs != 0)
1308 : {
1309 :
1310 0 : pkLocation = *pPlocs; move16(); /* N.B. No post-increment */
1311 0 : lastPeak = sub(*num_plocs, 1);
1312 0 : if (lastPeak >= 0)
1313 : {
1314 0 : pkLocation1 = *pPlocs++; move16(); /* get a next peak */
1315 : }
1316 : }
1317 :
1318 0 : FOR(m = 0; m < *num_plocs; m++)
1319 : {
1320 0 : pkLocation_1 = pkLocation; /* plocs[m - 1] */ move16();
1321 0 : pkLocation = pkLocation1; /* plocs[m] */move16();
1322 :
1323 : /*location dependent update of plocs[m + 1] */
1324 0 : if ( sub(lastPeak, m) > 0)
1325 : { /* only read additional peak when teher is a valid position to read */
1326 0 : pkLocation1 = *pPlocs++; /* plocs[m + 1] */ move16();
1327 : }
1328 :
1329 0 : delta_tmp = shr_pos(sub(sub(pkLocation, pkLocation_1), 1), 1);
1330 0 : if (m == 0)
1331 : {
1332 0 : delta_tmp = DELTA_CORR; move16(); /* first peak special case */
1333 : }
1334 0 : delta_corr_dn = s_min(delta_tmp, DELTA_CORR);
1335 :
1336 0 : delta_tmp = shr_pos(sub(sub(pkLocation1, pkLocation), 1), 1);
1337 0 : if (sub(m, lastPeak) >= 0)
1338 : {
1339 0 : delta_tmp = DELTA_CORR; move16();
1340 : }
1341 0 : delta_corr_up = s_min(delta_tmp, DELTA_CORR); /* last peak special case */
1342 :
1343 :
1344 : /* Input Xph */
1345 0 : segmentLen = sub(sub(pkLocation, delta_corr_dn), i); /* may be negative */
1346 :
1347 0 : ASSERT(pReX == &(X[i]));
1348 0 : ASSERT(*pReX == X[i]); /*before first, nth valley*/
1349 :
1350 0 : FOR(j = 0; j < segmentLen; j++) /* valley section , may be skipped for segmentlen < 0 */
1351 : {
1352 :
1353 0 : *seed = rand_phase_fx(*seed, &sin_F, &cos_F);
1354 :
1355 : /* phase scrambling */
1356 0 : rotate_W16_fx(*pReX, *pImX, cos_F, sin_F, &tmp, &im); /* tmp=real part, should be inlined */
1357 : UNUSED(re);
1358 :
1359 : /* i.e. add a bit of magnitude scrambling around 1.0 in longer bursts */
1360 0 : *seed = rand_phase_fx(*seed, &sin_F_fade2avg, &cos_F_fade2avg);
1361 0 : IF(noise_mag_scale_neg != 0)
1362 : {
1363 0 : valley_magnitude_adj_fx(&tmp, &im, *seed, cos_F);
1364 : /* use two random variables */ /*reuse cosF from regular X scrambling, + reuse *seed from fade2avg
1365 : */
1366 : }
1367 :
1368 0 : IF( beta[k] != 0 )
1369 : {
1370 0 : alpha_local = alpha[k]; /* no move as alpha_local is only needed for debug */
1371 : /* the fade2avg mixing branch, fade2avg and transient downscaling and evolution rotation */
1372 0 : tmp2 = mult_r(beta[k], Xavg[k]);
1373 :
1374 : {
1375 0 : tmp2 = s_and(tmp2, one_peak_flag_mask);
1376 0 : tmp = s_and(tmp, one_peak_flag_mask);
1377 0 : im = s_and(im, one_peak_flag_mask);
1378 : }
1379 :
1380 0 : *pReX++ = mac_r(L_mult(alpha_local, tmp), tmp2, cos_F_fade2avg); move16();
1381 0 : *pImX-- = mac_r(L_mult(alpha_local, im), tmp2, sin_F_fade2avg); move16();
1382 : }
1383 : ELSE
1384 : { /* the no fade2avg mixing branch, only transient downscaling and evolution rotation */
1385 : {
1386 0 : tmp = s_and(tmp, one_peak_flag_mask);
1387 0 : im = s_and(im, one_peak_flag_mask);
1388 : }
1389 :
1390 :
1391 0 : *pReX++ = mult_r(mag_chg[k], tmp); move16();
1392 0 : *pImX-- = mult_r(mag_chg[k], im); move16();
1393 : }
1394 0 : i = add(i, 1);
1395 0 : ASSERT(pReX == &(X[i])); /* nth Valley */
1396 0 : ASSERT(*pReX == X[i]); /* nth Valley */
1397 0 : if (sub(i, gwlpr_fxPlus1[k]) >= 0)
1398 : {
1399 0 : k = add(k, 1);
1400 : }
1401 : } /* segment_len 1st++ etc valley excluding last */
1402 :
1403 : /* peak area section */
1404 :
1405 0 : e = s_min(lprotBy2Minus1, add(pkLocation, delta_corr_up));
1406 0 : segmentLen = sub(e, sub(i, 1));
1407 :
1408 :
1409 0 : L_Xph = *pCorrPhase_L; move32();
1410 : /* rounding here, add "0.5" before extracting the 10 bits for table lookup */
1411 0 : Xph_short = s_and(extract_l(L_shr_pos_pos(L_add(L_Xph, (1L << (16 - 10 - 1))), 16 - 10)), 0x3ff);
1412 : /* 10 bits precision after radix point, for a virtual 0-1023 sin/cos table lookup */
1413 :
1414 :
1415 0 : get_sin_cosQ10opt(Xph_short, &peak_sin_F, &peak_cos_F);
1416 :
1417 :
1418 0 : ASSERT(pReX == &(X[i])); /*before peak*/
1419 0 : FOR(j = 0; j < segmentLen; j++)
1420 : {
1421 0 : mag_chg_local = mag_chg[k]; /* mag_chg_local actually only need for debugging , no move16()*/
1422 :
1423 : UNUSED(ph_dith);
1424 0 : *seed = extract_l(L_mac0(13849, *seed, 31821));
1425 0 : rotate_W16_fx(*pReX, *pImX, peak_cos_F, peak_sin_F, &tmp, &im); /* should be inlined */
1426 : UNUSED(re);
1427 :
1428 0 : *seed = rand_phase_fx(*seed, &sin_F, &cos_F);
1429 :
1430 0 : IF( beta[k] != 0 )
1431 : { /* fade2avg path alpha*X_sav + beta *Xavg */
1432 0 : alpha_local = mag_chg_local; /* no move alpha_local only needed for dbg */
1433 0 : tmp2 = mult_r(beta[k], Xavg[k]);
1434 :
1435 :
1436 0 : *pReX++ = mac_r(L_mult(alpha_local, tmp), tmp2, cos_F); move16();
1437 0 : *pImX-- = mac_r(L_mult(alpha_local, im), tmp2, sin_F); move16();
1438 : }
1439 : ELSE
1440 : {
1441 0 : *pReX++ = mult_r(mag_chg_local, tmp); move16();
1442 0 : *pImX-- = mult_r(mag_chg_local, im); move16();
1443 : }
1444 :
1445 0 : i = add(i, 1);
1446 0 : ASSERT(pReX == &(X[i]));
1447 0 : ASSERT(*pReX == X[i]);
1448 0 : if (sub(i, gwlpr_fxPlus1[k]) > 0)
1449 : {
1450 0 : k = add(k, 1);
1451 : }
1452 : } /* segment_length Peak*/
1453 0 : pCorrPhase_L++;
1454 : }
1455 :
1456 0 : segmentLen = sub(lprotBy2Minus1, sub(i, 1)); /* tail/valley noise-bins */
1457 :
1458 : /* for 48k the last 63+1+63 = values above 20 kHz */
1459 : { /* after last, peak no need to scramble above 20kHz , those coeffs have been/ will be zeroed already */
1460 : /* DC, Xused,ReZeroed, fs/2, ImZeroed, ImUsed */
1461 : /*N: 1, 320 , 63 , 1 , 63 , 320 */
1462 : /*ind: 0, 1-320, 321-383 ,384 , 385-447, 448 -767 */
1463 : /* segmentLen = sub(320, i-1); */ /*only process up to bin X[320] , 63 coeffs X[321] fwd should already be
1464 : zeroed */
1465 : /* ASSERT(i-1+ segmentLen == lprotBy2Minus1 ); */
1466 : }
1467 :
1468 0 : ASSERT(*pReX == X[i]); /* before a final valley*/
1469 0 : FOR(j = 0; j < segmentLen; j++)
1470 : {
1471 :
1472 :
1473 0 : *seed = rand_phase_fx(*seed, &sin_F, &cos_F);
1474 0 : rotate_W16_fx(*pReX, *pImX, cos_F, sin_F, &tmp, &im); /* should be inlined , tmp=real part */
1475 :
1476 0 : *seed = rand_phase_fx(*seed, &sin_F_fade2avg, &cos_F_fade2avg);
1477 : /* i.e. add a bit of magnitude scrambling around 1.0 in longer bursts */
1478 0 : IF(noise_mag_scale_neg != 0)
1479 : {
1480 0 : valley_magnitude_adj_fx(&tmp, &im, *seed, cos_F);
1481 : /*reuse cosF from regular X scrambling, + reuse *seed from fade2avg */
1482 : }
1483 :
1484 0 : tmp = s_and(tmp, one_peak_flag_mask);
1485 0 : im = s_and(im, one_peak_flag_mask);
1486 0 : IF(beta[k] != 0)
1487 : { /* fade2avg path never in first BFI frame */
1488 :
1489 0 : alpha_local = alpha[k]; /* no move as alpha_local is only needed for debug */
1490 0 : tmp2 = mult_r(beta[k], Xavg[k]);
1491 : {
1492 0 : tmp2 = s_and(tmp2, one_peak_flag_mask);
1493 : }
1494 :
1495 :
1496 0 : *pReX++ = mac_r(L_mult(alpha_local, tmp), tmp2, cos_F_fade2avg); move16();
1497 0 : *pImX-- = mac_r(L_mult(alpha_local, im), tmp2, sin_F_fade2avg); move16();
1498 : }
1499 : ELSE
1500 : {
1501 :
1502 0 : *pReX++ = mult_r(mag_chg[k], tmp); move16();
1503 0 : *pImX-- = mult_r(mag_chg[k], im); move16();
1504 : }
1505 0 : i = add(i, 1);
1506 0 : ASSERT(*pReX == X[i]);
1507 0 : if (sub(i, gwlpr_fxPlus1[k]) > 0)
1508 : {
1509 0 : k = add(k, 1);
1510 : }
1511 : } /* segment_len last valley */
1512 :
1513 :
1514 :
1515 :
1516 : #ifdef DYNMEM_COUNT
1517 : Dyn_Mem_Out();
1518 : #endif
1519 : #ifdef WMOPS
1520 : pop_wmops();
1521 : #endif
1522 0 : }
1523 :
1524 0 : void my_wtda_fx(const Word16 *new_audio, /* i : input audio to be windowed Q0 20 ms , OPT can be output as well */
1525 : const Word16 *const win2ms_init, /* i: 2 ms initial part of pre_tda window */
1526 : const Word16 *const win16ms_center, /* i: 16 ms combined part of pre_tda IWHR+MDCT-ana */
1527 : Word32 *L_wtda_audio, /* o : tda audio Q16 20 ms */
1528 : const Word16 L, Word8 *scratchBuffer) /* Size = 8 * MAX_L_FRAME */
1529 : {
1530 : Word16 i, L2; /*,L4;*/
1531 : const Word16 *pX;
1532 : const Word16 *pW;
1533 : Word32 *p1_L, *p2_L, *p3_L, *p4_L;
1534 : Word32 *L_w_audio; /*OPt 1 use input buffer , OPT2, may be shortened from 20 ms to 16 ms due to
1535 : zeroed parts*/
1536 : Word32 *pY_L, *pa1_L, *pa2_L;
1537 :
1538 : #ifdef DYNMEM_COUNT
1539 : Dyn_Mem_In("my_wtda_fx", sizeof(struct {
1540 : Word16 i, L2; /*,L4;*/
1541 : const Word16 *pX;
1542 : const Word16 *pW;
1543 : Word32 *p1_L, *p2_L, *p3_L, *p4_L;
1544 : Word32 *L_w_audio; /*OPT 1 use input buffer, OPT2, may be shortened from 20 ms to 16.25+(1.75) ms
1545 : due to zeroed parts*/
1546 : Word32 *pY_L, *pa1_L, *pa2_L;
1547 : }));
1548 : #endif
1549 :
1550 : #ifdef WMOPS
1551 : push_wmops("PhECU::my_wtda_fx");
1552 : #endif
1553 :
1554 0 : L_w_audio = (Word32 *)scratchAlign(scratchBuffer, 0); /* Size = 4 * 2 * MAX_L_FRAME */
1555 :
1556 : /* |111111|222222|333333|444444 */
1557 : /* |p1 p2| p3||p4 */
1558 :
1559 : /* Apply analysis window */
1560 :
1561 0 : pX = new_audio;
1562 0 : pY_L = L_w_audio;
1563 :
1564 0 : pW = win2ms_init;
1565 0 : FOR(i = 0; i < ((2 * L / 10)); i++) /* Loop over 2ms window MDCT-ana length */
1566 : {
1567 0 : *pY_L++ = L_mult(*pX++, *pW++); move32();
1568 : }
1569 :
1570 0 : pW = win16ms_center;
1571 : /* 20 ms - 2ms - 3.75 ms = 14.25 ms */
1572 0 : FOR(i = 0; i < (2 * L - (2 * L / 10) - ((3 * 2 * L) / 16));
1573 : i++) /* Loop over remaining 14.25ms, currently out of a 16 ms stored length */
1574 : {
1575 0 : *pY_L++ = L_mult(*pX++, *pW++); move32();
1576 : }
1577 :
1578 :
1579 : /* tda */
1580 0 : L2 = shr_pos_pos(L, 1); /* length of tda blocks */
1581 :
1582 0 : p1_L = L_w_audio; /* block 1 fwd */
1583 0 : p2_L = L_w_audio + L - 1; /* block 2 rev */
1584 0 : p3_L = L_w_audio + L + L2 - 1; /* block 3 rev */
1585 0 : p4_L = L_w_audio + L + L2; /* block 4 fwd */
1586 :
1587 0 : pa1_L = L_wtda_audio; /* first part output */
1588 0 : pa2_L = L_wtda_audio + L2; /* second part output */
1589 :
1590 0 : FOR(i = 0; i < (L / 8); i++)
1591 : {
1592 : /* first 1.25ms part of tda signal -p3_rev -p4_fwd */
1593 0 : *pa1_L++ = L_negate(L_add_sat(*p3_L--, *p4_L++)); move32();
1594 : }
1595 0 : FOR(; i < L2; i++)
1596 : {
1597 : /* first front part 3.75 ms p4_l is zeroes -p3_rev + (-p4_fwd==0) */
1598 0 : *pa1_L++ = L_negate(*p3_L--); move32();
1599 : }
1600 :
1601 0 : FOR(i = 0; i < L2; i++)
1602 : {
1603 : /* second part of tda signal p1 -p2_rev */
1604 0 : *pa2_L++ = L_sub_sat(*p1_L++, *p2_L--); move32();
1605 : }
1606 :
1607 : #ifdef DYNMEM_COUNT
1608 : Dyn_Mem_Out();
1609 : #endif
1610 : #ifdef WMOPS
1611 : pop_wmops();
1612 : #endif
1613 0 : return;
1614 : }
1615 :
1616 : /*--------------------------------------------------------------------------
1617 : * rec_wtda()
1618 : *
1619 : * Windowing and TDA of reconstructed frame
1620 : *--------------------------------------------------------------------------*/
1621 :
1622 :
1623 : static void
1624 0 : rec_wtda_fx(
1625 : Word16 *X, /* i : iFFT(X-evolved) TD-signal o: dbg 16 ms */
1626 : Word32 *L_ecu_rec, /* o : Reconstructed frame in tda domain Qx? */
1627 : const Word16 output_frame, /* i : Frame length */
1628 : const Word16 Lprot, /* i : Prototype frame length */
1629 : const Word16 *const win2ms_init, /* i: 2 ms initial part of pre_tda window */
1630 : const Word16 *const win16ms_center, /* i: 16 ms combined part of pre_tda IWHR+MDCT-ana */
1631 : const Word16 maxLen,
1632 : const Word16 *prevsynth,
1633 : const Word16 Q_psMinus1, /*i: Q prev_synth minus 1 , (-1 to match Q of Xsav in first bfi frame ) */
1634 : Word8 *scratchBuffer) /* Size = 12 * MAX_L_FRAME */
1635 : {
1636 : Word16 l, Lprot2;
1637 : Word16 *rec_buf;
1638 : Word16 *xsubst_; /*,*out_ptr;*/
1639 :
1640 : Word16 ola_old[COPY_LEN_48K + OLA_LEN_48K]; /* 3.75 ms */
1641 : Word16 work_len;
1642 : Word16 copy_len;
1643 : Word16 ola_len;
1644 :
1645 :
1646 : Word8 *buffer_wtda;
1647 :
1648 : #ifdef DYNMEM_COUNT
1649 : Dyn_Mem_In("rec_wtda_fx", sizeof(struct {
1650 : Word16 l, Lprot2;
1651 : Word16 *rec_buf;
1652 : Word16 *xsubst_; /*,*out_ptr;*/
1653 : Word8 *buffer_wtda;
1654 : }));
1655 : #endif
1656 :
1657 : #ifdef WMOPS
1658 : push_wmops("PhECU::rec_wtda_fx");
1659 : #endif
1660 :
1661 0 : rec_buf = scratchAlign(scratchBuffer, 0); /* Size = 2 * 2 * MAX_L_FRAME */
1662 0 : buffer_wtda = (Word8 *)scratchAlign(rec_buf, sizeof(*rec_buf) * (2 * maxLen)); /* Size = 4 * 2 * MAX_L_FRAME */
1663 :
1664 0 : xsubst_ = rec_buf;
1665 0 : Lprot2 = shr_pos_pos(Lprot, 1);
1666 :
1667 : /* extract reconstructed frame with ld window into rec_buf */
1668 0 : l = sub(output_frame, Lprot2);
1669 0 : basop_memmove(xsubst_ + l, X, (Lprot) * sizeof(Word16)); /* 16 ms IFFT raw output */
1670 0 : copy_len = COPY_LEN[FRAME2FS_IDX(output_frame)];
1671 0 : ola_len = OLA_LEN[FRAME2FS_IDX(output_frame)];
1672 0 : work_len = add(copy_len, ola_len);
1673 : /* Copy and scale copy part 2ms from prevsynth */
1674 0 : basop_memmove(rec_buf, &prevsynth[Lprot - work_len], (copy_len) * sizeof(Word16)); /* 2ms out of 3.75 ms copied */
1675 0 : Scale_sig_sat(rec_buf, copy_len, sub(-3, Q_psMinus1)); /* inplace scaling by 2^(-Q_ps-4) */
1676 :
1677 :
1678 : /* Copy, window and scale 1.75 ms ola part from prevsynth , into a temporary buffer ola_old */
1679 :
1680 :
1681 0 : windowing_ola(&prevsynth[Lprot - ola_len], ola_old, w_old[FRAME2FS_IDX(output_frame)], ola_len); /* 1.75 ms */
1682 :
1683 0 : Scale_sig_sat(ola_old, ola_len, sub(-3, Q_psMinus1)); /* inplace scaling by 2^(-Qps-4) */
1684 :
1685 : /* Window 1.75 ms length inplace recreated X=IFFT(prototype) signal copied into rec_buf signal */
1686 0 : windowing_ola(&rec_buf[copy_len], &rec_buf[copy_len], w_new[FRAME2FS_IDX(output_frame)], ola_len);
1687 : /* mix add the two windowed components */
1688 0 : ola_add(&rec_buf[copy_len], ola_old, &rec_buf[copy_len], ola_len); /* OPT: sat check */
1689 :
1690 0 : my_wtda_fx(rec_buf,
1691 : win2ms_init, /* i: 2 ms initial part of pre_tda window */
1692 : win16ms_center, /* i: 16 ms combined part of pre_tda IWHR+MDCT-ana */
1693 : L_ecu_rec, output_frame, buffer_wtda); /* */
1694 :
1695 :
1696 : #ifdef DYNMEM_COUNT
1697 : Dyn_Mem_Out();
1698 : #endif
1699 : #ifdef WMOPS
1700 : pop_wmops();
1701 : #endif
1702 :
1703 0 : return;
1704 : }
1705 :
1706 : /*--------------------------------------------------------------------------
1707 : * rec_frame_fx()
1708 : *
1709 : * Frame reconstruction
1710 : *--------------------------------------------------------------------------*/
1711 0 : void rec_frame_fx(Word16 * x, /* i : FFT spectrum 16 ms, o: ifft() TD debug signal 16ms */
1712 : Word32 * L_ecu_rec, /* o : Reconstructed frame in tda domain 10ms buffer */
1713 : const Word16 output_frame, /* i : Frame length */
1714 : const Word16 Q, /* i : Xsav Q */
1715 : const Word16 *const win2ms_init, /* i: 2 ms initial part of pre_tda window */
1716 : const Word16 *const win16ms_center, /* i: 16 ms combined part of pre_tda IWHR+MDCT-ana */
1717 : Word16 maxLprot,
1718 : const Word16 *prevsynth,
1719 : const Word16 Q_prevsynthMinus1, /* i : prevsynthQ-1 or xfp Q */
1720 : Word8 *scratchBuffer /* Size = 4 * MAX_LPROT + 12 * MAX_L_FRAME */
1721 : )
1722 : {
1723 : Counter i;
1724 : Word16 Lprot;
1725 : Word32 *L_x;
1726 : Word32 *pX_L;
1727 : Word16 *pX;
1728 : Word16 fft_scale;
1729 : Word8 * buffer_fft;
1730 :
1731 : #ifdef DYNMEM_COUNT
1732 : Dyn_Mem_In("rec_frame_fx", sizeof(struct {
1733 : Counter i;
1734 : Word16 Lprot;
1735 : Word32 *L_x;
1736 : Word32 *pX_L;
1737 : Word16 *pX;
1738 : Word16 fft_scale;
1739 : Word8 * buffer_fft;
1740 : }));
1741 : #endif
1742 :
1743 : #ifdef WMOPS
1744 : push_wmops("PhECU::rec_frame_fx");
1745 : #endif
1746 :
1747 0 : L_x = (Word32 *)scratchAlign(scratchBuffer, 0); /* Size = 4 * MAX_LPROT */
1748 0 : buffer_fft = (Word8 *)scratchAlign(L_x, sizeof(*L_x) * maxLprot); /* Size = 4* (2+1) * MAX_L_FRAME */
1749 :
1750 : /* Initialize to FB constants */
1751 0 : Lprot = mult(output_frame, (Word16)(32768.0 / 99.0)); /* truncation needed , i.e no rounding can be applied here */
1752 0 : ASSERT(Lprot == (output_frame / 100));
1753 0 : Lprot = LprotSzPtr[Lprot]; move16();
1754 :
1755 :
1756 : /* Convert stored 16 bit into 32bit for fft */
1757 0 : flippedW16_2_intlvW32(x, sub(shr_pos_pos(Lprot, 1), 1), Lprot, L_x);
1758 : /*scratch x now free for re-use */
1759 :
1760 0 : fft_scale = -1; move16();
1761 :
1762 : #ifdef WMOPS
1763 : push_wmops("PhECU::IFFT_fx");
1764 : #endif
1765 0 : BASOP_irfftN(L_x, Lprot, &fft_scale, buffer_fft);
1766 : #ifdef WMOPS
1767 : pop_wmops();
1768 : #endif
1769 :
1770 0 : pX_L = &L_x[0];
1771 0 : pX = &x[0]; /* scratch x reused */
1772 :
1773 : {
1774 0 : FOR(i = 0; i < Lprot; i++)
1775 : {
1776 0 : *pX++ = extract_h(L_shl_sat(*pX_L++, fft_scale)); move16();
1777 : }
1778 : }
1779 : /* Scratch L_x may be released */
1780 :
1781 :
1782 : /* Saturation possible when rescaling to Q0 if there are random bit errors in the bit stream
1783 : * One may use one guard bit to better handle this - though it should not be needed for normal
1784 : * operation.
1785 : */
1786 0 : Scale_sig_sat(x, Lprot, negate(Q)); /* scale by 2^(-Q) */
1787 :
1788 : /*scratch x is Lprot word16 */
1789 : /* scatch need for win-TDA XXX */
1790 : /* scrcatch need for TDA XXX */
1791 :
1792 0 : rec_wtda_fx(x, L_ecu_rec, output_frame, Lprot,
1793 : win2ms_init, /* i: 2 ms initial part of pre_tda window */
1794 : win16ms_center, /* i: 16 ms combined part of pre_tda IWHR+MDCT-ana */
1795 0 : s_max(output_frame, 160),
1796 : prevsynth,
1797 : Q_prevsynthMinus1, /* NB: prevsynth Q may change from prev_bfi=0 to prev_bfi==1, due to phase
1798 : of signal or PLC2-muting */
1799 : buffer_fft);
1800 :
1801 :
1802 : #ifdef DYNMEM_COUNT
1803 : Dyn_Mem_Out();
1804 : #endif
1805 : #ifdef WMOPS
1806 : pop_wmops();
1807 : #endif
1808 0 : return;
1809 : }
1810 :
1811 : /*--------------------------------------------------------------------------
1812 : * hq_phase_ecu_fx()
1813 : *
1814 : * Main routine for HQ phase ECU
1815 : *--------------------------------------------------------------------------*/
1816 0 : void hq_phase_ecu_fx(
1817 : const Word16 *prevsynth, /* i : buffer of previously synthesized signal currently 16ms */
1818 : Word32 * L_ecu_rec, /* o : reconstructed frame in tda domain , also tmp w32_fft buffer */
1819 : Word16 * time_offs, /* i/o: Sample offset for consecutive frame losses*/
1820 : Word16 * X_sav, /* i/o: Stored spectrum of prototype frame */
1821 : Word16 * Q_spec, /* o: Q value of stored spectrum */
1822 : Word16 * num_p, /* i/o: Number of identified peaks */
1823 : Word16 * plocs, /* i/o: Peak locations */
1824 : Word32 * L_plocsi, /* i/o: Interpolated peak locations Q16 */
1825 : const Word16 env_stab, /* i : Envelope stability parameter */
1826 : const Word16 f0hzLtpBinQ7, /* i: LTP bin frequency in normalized Hz Q7 */
1827 : const Word16 norm_corrQ15_fx, /*i : correlation for lag at f0hzLtpBinQ7 */
1828 : const Word16 prev_bfi, /* i : indicating burst frame error */
1829 : Word16 old_is_transient[2], /* i/o : flags indicating noise generation frames */
1830 : Word16 * mag_chg_1st, /* i/o: per band magnitude modifier for transients*/
1831 : Word16 * mag_chg_gr, /* o: per band magnitude modifier, incl burst attenuation */
1832 : Word16 * Xavg, /* i/o: Frequency group average gain to fade to */
1833 : Word16 * beta_mute, /* o : Factor for long-term mute */
1834 : const Word16 bwidth_fx, /* i : Encoded bandwidth */
1835 : const Word16 output_frame, /* i : frame length */
1836 : Word16 * seed_out_fxPtr, /* o: utput dbg NULL may be used*/
1837 : Word16 * X_out, /* o: utput dbg NUll may be used */
1838 : const Word16 t_adv, /* i : time adjustment including time_offs */
1839 : const Word16 *const win2ms_init, /* i: 2 ms initial part of pre_tda window */
1840 : const Word16 *const win16ms_center, /* i: 16 ms combined part of pre_tda IWHR+MDCT-ana */
1841 : const Word16 *sp_ana_win, /* i : whr 3+10+3 window */
1842 : Word16 q_fx_old_exp, /* i : exp of prev_synth */
1843 :
1844 : Word16 maxLprot, /* i : maz spectrum buffer size */
1845 : Word16 maxPlocs, /* i : max nb of peaks */
1846 : Word32 L_oold_xfp_w_E_fx, Word16 oold_xfp_w_E_exp_fx, /* exp of time signale */
1847 : Word16 oold_Ltot_exp_fx, /*true exp of energy */
1848 : Word16 *oold_grp_shape_fx, Word32 L_old_xfp_w_E_fx, Word16 old_xfp_w_E_exp_fx, /* exp of time signale */
1849 : Word16 old_Ltot_exp_fx, /*true exp of energy */
1850 : Word16 *old_grp_shape_fx,
1851 : Word16 margin_prevsynth,
1852 : const Word16 fadeout,
1853 : Word16 *nonpure_tone_flag_ptr, /* i/o : non-pure single tone indicator state */
1854 : Word8 *scratchBuffer /* Size = 2 * MAX_LGW + 8 * MAX_LPROT + 12 * MAX_L_FRAME */
1855 : )
1856 : {
1857 : Word16 lprot;
1858 : Word16 *mag_chg, ph_dith, *X;
1859 : Word16 *xfp;
1860 : Word16 seed;
1861 : Word16 alpha[MAX_LGW], beta[MAX_LGW] ;
1862 : Word16 prevsynth_man_upshift;
1863 : Word16 Q_prevsynthMinus1;
1864 : Word8 *buffer;
1865 : Word32 L_Xavg[MAX_LGW]; /* i/o : Frequency group amp averages for tonal tilt analysis Max upshifted */
1866 : #ifdef DYNMEM_COUNT
1867 : Dyn_Mem_In("hq_phase_ecu_fx", sizeof(struct {
1868 : Counter i;
1869 : Word16 lprot;
1870 : Word16 *mag_chg, ph_dith, *X;
1871 : Word16 seed;
1872 : Word16 alpha[MAX_LGW], beta[MAX_LGW];
1873 : Word16 prevsynth_man_upshift;
1874 : Word16 Q_prevsynthMinus1;
1875 : Word8 *buffer;
1876 : // ToDO Word32 L_Xavg[MAX_LGW];
1877 : }));
1878 : #endif
1879 :
1880 : if (!prev_bfi)
1881 : {
1882 : #ifdef WMOPS
1883 : push_wmops("PhECU::hq_phase_ecu_fx(1st)");
1884 : #endif
1885 : }
1886 : else
1887 : {
1888 : #ifdef WMOPS
1889 : push_wmops("PhECU::hq_phase_ecu_fx(N)");
1890 : #endif
1891 : }
1892 :
1893 0 : mag_chg = (Word16 *)scratchAlign(scratchBuffer, 0); /* Size = 2 * MAX_LGW */
1894 0 : X = (Word16 *)scratchAlign(mag_chg, sizeof(*mag_chg) * MAX_LGW); /* Size = 2 * MAX_LPROT == 1 Word16*MAX_LPROT */
1895 0 : xfp = (Word16 *)scratchAlign(X, sizeof(*X) * maxLprot); /* Size = 2 * MAX_LPROT == 1 Word16*MAX_LPROT */
1896 0 : buffer = (Word8 *)scratchAlign(xfp, sizeof(*xfp) * maxLprot); /* Size = 4 * MAX_LPROT + 12 * MAX_L_FRAME */
1897 :
1898 : /* buffer size = Word32 * MAX_LPROT (FFT, IFFT) DRAM
1899 : + 3*Word32 * MAX_L_FRAME */
1900 :
1901 0 : basop_memset(alpha, 0, MAX_LGW*sizeof(Word16));
1902 0 : basop_memset(beta, 0, MAX_LGW*sizeof(Word16));
1903 :
1904 0 : lprot = LprotSzPtr[bwidth_fx]; move16();
1905 :
1906 0 : test();
1907 0 : ASSERT(prev_bfi >= 0 && prev_bfi <= 1);
1908 0 : IF( prev_bfi == 0 ) /* inside PhECU we can check vs 0 */
1909 : {
1910 0 : *time_offs = 0; move16();
1911 0 : *nonpure_tone_flag_ptr = -1; move16(); /* flag nonpure tone flag for new analysis */
1912 : /* analysis made outside, up/down scaling here from static RAM to dynamic RAM */
1913 : /* prevsynth_in_flt = prev_synth_man*2.^(-15 + exp_old) */
1914 : /* X_sav_flt = X_man/2.^(Q_spec) */
1915 :
1916 : /* 1 bit headroom needed in xfp_tmp_buf for FFT processing into X_sav */
1917 : /* margin_prevsynth = actual margin in incoming 16 ms xpf segment */
1918 : /* q_fx_old_exp was computed on full 21+ ms updatepcm buffer) */
1919 :
1920 0 : ASSERT(margin_prevsynth >= 0 && (margin_prevsynth <= 16));
1921 :
1922 : prevsynth_man_upshift =
1923 0 : sub(margin_prevsynth, 1); /* 0 --> -1(==down), 1--> 0 , 2 --> 1(==up), ... */
1924 0 : ASSERT(prevsynth_man_upshift >= -16 && prevsynth_man_upshift <= 15); /* avoid Overflow in shr, shl */
1925 0 : *Q_spec = sub(15, sub(q_fx_old_exp, prevsynth_man_upshift));
1926 : /* Q_spec target is to create 1 additional bit of margin in the xfp TD buffer , before converting to
1927 : * Xsav
1928 : */
1929 :
1930 0 : if (margin_prevsynth == 0)
1931 : {
1932 0 : Word16 Qold = 15 - (q_fx_old_exp + 1);
1933 0 : Word16 tmp_man_upshift = margin_prevsynth - 1; /* 1 .. -15 */
1934 0 : Word16 Qnew = 15 - (q_fx_old_exp - tmp_man_upshift);
1935 0 : assert(Qold == Qnew);
1936 : UNUSED(Qold);
1937 : UNUSED(Qnew);
1938 : }
1939 0 : if (margin_prevsynth == 1)
1940 : {
1941 0 : Word16 Qold = 15 - (q_fx_old_exp + 0);
1942 0 : assert(*Q_spec == Qold);
1943 : UNUSED(Qold);
1944 : }
1945 0 : if (margin_prevsynth == 2)
1946 : {
1947 0 : Word16 Qold = 15 - (q_fx_old_exp - 1);
1948 0 : assert(*Q_spec == Qold);
1949 : UNUSED(Qold);
1950 : }
1951 :
1952 0 : Q_prevsynthMinus1 = sub(15, add(q_fx_old_exp, +1)); /* dbg to use non-scaled prevsynth for OLA */
1953 :
1954 : /* Q of prev_synth now separated from prev_bfi=0 and prev_bfi==1 */
1955 :
1956 : #ifdef USE_TMPXFP_IN_OLA
1957 : Q_prevsynthMinus1 = *Q_spec; /* a possibly upscaled Q */
1958 : #endif
1959 :
1960 :
1961 0 : Copy_Scale_sig(
1962 : prevsynth, xfp, lprot,
1963 : prevsynth_man_upshift); /* unscaled prevsynth is still used by rec_frame, copy to a temporary
1964 : xfp analysis buffer, and scale down to a margin of 1 bit */
1965 :
1966 :
1967 0 : trans_burst_ana_fx(
1968 : ((void *)NULL),
1969 0 : mag_chg, &ph_dith, mag_chg_1st, output_frame, *time_offs, env_stab, alpha, beta, beta_mute, Xavg,
1970 0 : (*Q_spec), L_oold_xfp_w_E_fx, oold_xfp_w_E_exp_fx, oold_Ltot_exp_fx, oold_grp_shape_fx,
1971 : L_old_xfp_w_E_fx, old_xfp_w_E_exp_fx, old_Ltot_exp_fx, old_grp_shape_fx,
1972 : fadeout,
1973 : L_Xavg, /* full scale band amplitudes in first bfi frame */
1974 : buffer);
1975 :
1976 0 : spec_ana_fx(&(xfp[0]), plocs, L_plocsi, num_p, X_sav, output_frame, bwidth_fx,
1977 : sp_ana_win, f0hzLtpBinQ7, norm_corrQ15_fx, maxLprot, maxPlocs, buffer);
1978 : }
1979 : ELSE
1980 : {
1981 : /* analysis made outside, possible up/down scaling here from static RAM to dynamic RAM */
1982 0 : ASSERT(margin_prevsynth >= 0 && (margin_prevsynth <= 15));
1983 :
1984 : /* prev_synth may have been scaled outside by update_pcm() function */
1985 : /* prevsynth_in_flt = prev_synth_man*2.^(-15 + exp_old) */
1986 : /* first bfi frame: X_sav_flt = X_man/2.^(*Q_spec) */
1987 : /* this bfi frame prev_synth(TD) is mixed with ifft signal(TD,from Q_spec) in rec_frame_fx
1988 : variable Q_prevsynthMinus1 provided to rec_frame for this TD+TD mixing
1989 : */
1990 :
1991 0 : Q_prevsynthMinus1 = sub(15, add(q_fx_old_exp, +1)); /* for now, "old" use of unscaled prevsynth for added OLA */
1992 :
1993 : #ifdef USE_TMPXFP_IN_OLA
1994 : { /* DBG target here is to instead upshift to a zero margin */
1995 : /* but Q value is maintained to match the save Xsav *Q_spec which is used with an offset of " -1" */
1996 : /* prevsynth_man_upshift = -prevsynth_margin ; */
1997 : Word16 tmp, worklen = add(hamm_len2Tab[bwidth_fx], shr(hamm_len2Tab[bwidth_fx], 2));
1998 : tmp = mult(7680, lprot); /* needs to truncate to integer */
1999 : ASSERT(tmp == worklen);
2000 : Copy_Scale_sig(&prevsynth[lprot - worklen], &xfp[lprot - worklen], worklen,
2001 : margin_prevsynth); /* only scale up last 3.75 ms */
2002 : Q_prevsynthMinus1 = sub(15, add(q_fx_old_exp, margin_prevsynth)); /*-1 does not make sense */
2003 : }
2004 : #endif
2005 :
2006 0 : *time_offs = add_sat(*time_offs, output_frame);
2007 :
2008 0 : trans_burst_ana_fx(
2009 : ((void *)NULL),
2010 0 : mag_chg, &ph_dith, mag_chg_1st, output_frame, *time_offs, env_stab, alpha, beta, beta_mute, Xavg,
2011 : (0), /* *Q_spec input only used in first bfi frames for burst analysis */
2012 : L_oold_xfp_w_E_fx, oold_xfp_w_E_exp_fx, oold_Ltot_exp_fx, oold_grp_shape_fx, L_old_xfp_w_E_fx,
2013 : old_xfp_w_E_exp_fx, old_Ltot_exp_fx, old_grp_shape_fx,
2014 : fadeout,
2015 : NULL, /* full scale band amplitudes , only used in first bfi frame */
2016 : buffer);
2017 : }
2018 : /* cpy LPROT Word16 from Static RAM Xsav to working DRAM/scratch buffer X ;*/
2019 0 : basop_memmove(X, X_sav, (lprot) * sizeof(Word16));
2020 : /* seed for own_rand2 */
2021 0 : seed = *time_offs; move16();
2022 :
2023 0 : if (mag_chg_gr != NULL) /* o: dbg per band magnitude modifier, incl. burst attenuation */
2024 : {
2025 : Counter k;
2026 0 : Word16 lgw_local = s_min(add(bwidth_fx, LGW8K), LGW48K); /* 4,5,6,7, (7/8) */
2027 0 : for (k = 0; k < lgw_local; k++)
2028 : {
2029 0 : mag_chg_gr[k] = mag_chg[k]; /* dbg SNR output */
2030 : }
2031 : }
2032 :
2033 0 : subst_spec_fx(plocs, L_plocsi, num_p, *time_offs, X, mag_chg, ph_dith, old_is_transient, output_frame,
2034 : &seed, alpha, beta,
2035 : Xavg, t_adv
2036 : ,fadeout,
2037 : nonpure_tone_flag_ptr, /* i/o : non-pure single tone indicator state */
2038 : L_Xavg /*i : only used in first bfi frame */
2039 : );
2040 :
2041 0 : if (seed_out_fxPtr != NULL)
2042 : {
2043 0 : *seed_out_fxPtr = seed; /* verify seed synch after subst_spec */
2044 : }
2045 :
2046 :
2047 0 : if (X_out != NULL)
2048 : {
2049 : Word16 ii;
2050 0 : for (ii = 0; ii < lprot; ii++)
2051 : {
2052 0 : X_out[ii] = X[ii]; /* evolve spectrum, no moves counted as this is a dbg-vector info cpy */
2053 : }
2054 : }
2055 :
2056 : /* reconstructed frame in tda domain */
2057 : /* NB *Q_spec only updated in first bfi frame */
2058 :
2059 : /*Scratch Analysis at this point */
2060 : /* X fft input needed as evolved spectrum Lprot* Word16 */
2061 : /* xfp not needed any longer Word16* Lprot */
2062 : /* IFFT operation needs Word32*Lprot in and Word32*Lprot out */
2063 : /* MDCT win operation may be inplace 2x MAX_L_FRAME */
2064 : /* TDA output is 1x MAX_L_FRAME */
2065 :
2066 0 : rec_frame_fx(X, L_ecu_rec, output_frame, *Q_spec,
2067 : win2ms_init, win16ms_center,
2068 : maxLprot,
2069 : #ifdef USE_TMPXFP_IN_OLA
2070 : xfp, /* only last 3.75 ms used, in both prevbfi=0 and prevBfi=1 i.e frames, xfp is an
2071 : upscaled prev_synth */
2072 : #else
2073 : prevsynth, /*only last 3.75 ms used in both prevbfi=0 and prevBfi=1 i.e frames */
2074 : #endif
2075 : Q_prevsynthMinus1,
2076 : buffer);
2077 :
2078 :
2079 : #ifdef DYNMEM_COUNT
2080 : Dyn_Mem_Out();
2081 : #endif
2082 : #ifdef WMOPS
2083 : pop_wmops();
2084 : #endif
2085 0 : }
2086 :
2087 :
2088 0 : static Word16 sqrt2ndOrder( /* o: in Q15 (2nd order least square approx.) */
2089 : const Word16 x /* i: x must be in between 0.5 and 1.0 (Q15). */
2090 : )
2091 : {
2092 : Word32 acc;
2093 : Word16 z;
2094 : #ifdef WMOPS
2095 : push_wmops("PhECU::sqrt2ndOrder");
2096 : #endif
2097 :
2098 0 : ASSERT(x >= 16384);
2099 :
2100 0 : acc = 1890205600L; /* 0.880195572812922 in Q31 */ move32();
2101 0 : z = mac_r(acc, x, -6506); /* -0.198537395405340 in Q15 */
2102 0 : acc = 682030261L; /* 0.317595089462249 in Q31 */ move32();
2103 0 : z = mac_r(acc, z, x); /* in Q15 */
2104 :
2105 :
2106 : #ifdef WMOPS
2107 : pop_wmops();
2108 : #endif
2109 0 : return z;
2110 : }
2111 :
2112 :
2113 : /* Modified version to produce Word32 input to FFT */
2114 0 : static void windowing_L(
2115 : const Word16 *x, /* i: Input signal */
2116 : Word32 * L_y, /* o: Windowed output */
2117 : const Word16 *win, /* i: Window coefficients */
2118 : const Word16 rectLength, /* i: Offset in between the 1st and 2nd symmetric halves of the Hamming window */
2119 : const Word16 halfLength /* i: Half of the total length of a complete Hamming window. */
2120 : )
2121 : {
2122 : Counter i;
2123 : Word32 * pY_L;
2124 : const Word16 *pX, *pW;
2125 : Word16 tmp_RL;
2126 :
2127 : #ifdef DYNMEM_COUNT
2128 : Dyn_Mem_In("windowing_L", sizeof(struct {
2129 : Counter i;
2130 : Word32 * pY_L;
2131 : const Word16 *pX, *pW;
2132 : }));
2133 : #endif
2134 :
2135 : #ifdef WMOPS
2136 : push_wmops("PhECU::windowing_L");
2137 : #endif
2138 :
2139 0 : pX = x;
2140 0 : pW = win;
2141 0 : pY_L = L_y;
2142 :
2143 0 : FOR(i = 0; i < halfLength; i++) /* 1st symmetric half of the Hamming window */
2144 : {
2145 0 : *pY_L++ = L_mult(*pX++, *pW++); move32();
2146 : }
2147 : /* Periodic filter - one more rect sample before end tapering */
2148 :
2149 0 : tmp_RL = add(rectLength, 1);
2150 :
2151 0 : if (rectLength == 0)
2152 : {
2153 0 : tmp_RL = 0; move16();
2154 : }
2155 : /* If rectLength is zero, it's a pure Hamming window; otherwise Hamming-Rectangular. */
2156 0 : FOR(i = 0; i < tmp_RL;i++)
2157 : {
2158 0 : *pY_L++ = L_deposit_h(*pX++); move32();
2159 : }
2160 0 : tmp_RL = sub(halfLength, 1);
2161 0 : if (rectLength == 0)
2162 : {
2163 0 : tmp_RL = halfLength; move16();
2164 : }
2165 0 : FOR(i = 0; i < tmp_RL; i++) /* 2nd symmetric half of the Hamming window. */
2166 : {
2167 0 : *pY_L++ = L_mult(*pX++, *(--pW)); move32();
2168 : }
2169 :
2170 : #ifdef DYNMEM_COUNT
2171 : Dyn_Mem_Out();
2172 : #endif
2173 : #ifdef WMOPS
2174 : pop_wmops();
2175 : #endif
2176 0 : }
2177 :
2178 : /* Modified version to produce Word16 input to FFT */
2179 0 : static void windowing_ola(const Word16 *x, /* i: Input signal */
2180 : Word16 * y, /* o: Windowed output */
2181 : const Word16 *win, /* i: Window coefficients */
2182 : const Word16 Length /* i: Half of the total length of a complete Hamming window. */
2183 : )
2184 : {
2185 : Counter i;
2186 : Word16 * pY;
2187 : const Word16 *pX, *pW;
2188 :
2189 : #ifdef DYNMEM_COUNT
2190 : Dyn_Mem_In("windowing_ola", sizeof(struct {
2191 : Counter i;
2192 : Word16 * pY;
2193 : const Word16 *pX, *pW;
2194 : }));
2195 : #endif
2196 :
2197 : #ifdef WMOPS
2198 : push_wmops("PhECU::windowing_ola");
2199 : #endif
2200 :
2201 0 : pX = x;
2202 0 : pW = win;
2203 0 : pY = y;
2204 :
2205 0 : FOR(i = 0; i < Length; i++) /* 1st symmetric half of the Hamming window */
2206 : {
2207 0 : *pY++ = mult_r(*pX++, *pW++); move16();
2208 : }
2209 :
2210 : #ifdef DYNMEM_COUNT
2211 : Dyn_Mem_Out();
2212 : #endif
2213 : #ifdef WMOPS
2214 : pop_wmops();
2215 : #endif
2216 0 : }
2217 :
2218 : /* Modified version to produce Word16 input to FFT */
2219 0 : static void ola_add(const Word16 *x, /* i: Input signal 1 */
2220 : const Word16 *y, /* i: Input signal 2 */
2221 : Word16 * z, /* o: Output signal */
2222 : const Word16 Length /* i: Half of the total length of a complete window. */
2223 : )
2224 : {
2225 : Counter i;
2226 : Word16 * pZ;
2227 : const Word16 *pX, *pY;
2228 :
2229 : #ifdef DYNMEM_COUNT
2230 : Dyn_Mem_In("windowing_ola", sizeof(struct {
2231 : Counter i;
2232 : Word16 * pY;
2233 : const Word16 *pX, *pW;
2234 : }));
2235 : #endif
2236 :
2237 :
2238 :
2239 0 : pX = x;
2240 0 : pY = y;
2241 0 : pZ = z;
2242 :
2243 0 : FOR(i = 0; i < Length; i++) /* 1st symmetric half of the Hamming window */
2244 : {
2245 0 : *pZ++ = add_sat(*pX++, *pY++); move16();
2246 : }
2247 :
2248 : #ifdef DYNMEM_COUNT
2249 : Dyn_Mem_Out();
2250 : #endif
2251 :
2252 0 : }
2253 :
2254 : /*-----------------------------------------------------------------------------
2255 : * magnSqrtApprox_fx()
2256 : *
2257 : * Approximation of sqrt(Square magnitude) of fft spectrum
2258 : * if min_abs <= 0.4142135*max_abs
2259 : * abs = 0.99*max_abs + 0.197*min_abs
2260 : * else
2261 : * abs = 0.84*max_abs + 0.561*min_abs
2262 : * end
2263 : *
2264 : * Note: even to handle the dynamics of sqrt(re^2+im^2) located on
2265 : * a scaled unit circle. One need to scale down the results
2266 : * with a factor 2, that is Q_out = Q_in - 1
2267 : * sqrt(32768.^2+32768.^2) results in = 23170 Q0-1,
2268 : * which corresponds to 46341 in the Q0 domain
2269 : *----------------------------------------------------------------------------*/
2270 0 : Word16 sqrtMagnApprox_fx( /* o : sqrt of magnitude square spectrum Q_in-1*/
2271 : const Word16 re, /* i : Real part Q_in */
2272 : const Word16 im /* i : Imag part Q_in */
2273 : )
2274 : {
2275 : /* Constants for Approximation of sqrt(Square magnitude) of fft spectrum
2276 : * >> num2str(round(0.4142135]*2.^15))
2277 : * ans = 13573
2278 : * >> num2str(round([0.99 0.197 0.84 0.561]*2.^14))
2279 : * ans = 16220 3228 13763 9191
2280 : */
2281 :
2282 : #define C_0p4142135_Q15 13573
2283 :
2284 : #define C_0p99_Q14 16220
2285 : #define C_0p197_Q14 3228
2286 : #define C_0p84_Q14 13763
2287 : #define C_0p561_Q14 9191
2288 : Word16 sgn_bit, re_abs, im_abs, max_abs, min_abs, sum;
2289 0 : Word16 jcoeffs[2][2] = { {C_0p99_Q14, C_0p197_Q14}, {C_0p84_Q14, C_0p561_Q14} };
2290 :
2291 : #ifdef DYNMEM_COUNT
2292 : Dyn_Mem_In("sqrtMagnApprox_fx", sizeof(struct {
2293 : Word16 sgn_bit, re_abs, im_abs, max_abs, min_abs, sum;
2294 : Word16 jcoeffs[2][2];
2295 : }));
2296 : #endif
2297 :
2298 : #ifdef WMOPS
2299 : push_wmops("PhECU::sqrtMagnApprox_fx");
2300 : #endif
2301 :
2302 : /* Get values and move pointers */
2303 0 : re_abs = abs_s(re); /* 1 cycle */
2304 0 : im_abs = abs_s(im); /* 1 cycle */
2305 :
2306 : /* Find max and min value */
2307 0 : min_abs = s_min(re_abs, im_abs); /* 1 cycle */
2308 0 : max_abs = s_max(re_abs, im_abs); /* 1 cycle */
2309 :
2310 : /* Calc approximation depending on relation */
2311 0 : sgn_bit = lshr(sub(mult(max_abs, C_0p4142135_Q15), min_abs), 15); /* 3 cycles */
2312 0 : sum = mac_r(L_mult(max_abs, jcoeffs[sgn_bit][0]), min_abs, jcoeffs[sgn_bit][1]); /* 2 cycles */
2313 : #ifdef DYNMEM_COUNT
2314 : Dyn_Mem_Out();
2315 : #endif
2316 : #ifdef WMOPS
2317 : pop_wmops();
2318 : #endif
2319 0 : return sum;
2320 : }
2321 :
2322 : /*-----------------------------------------------------------------------------
2323 : * fft_spec2_sqrt_approx_fx()
2324 : *
2325 : * Approximation of sqrt(Square magnitude) of fft spectrum
2326 : * if min_abs <= 0.4142135*max_abs
2327 : * abs = 0.99 max_abs + 0.197*min_abs
2328 : * else
2329 : * abs = 0.84 max_abs + 0.561*min_abs
2330 : * end
2331 : *
2332 : * Note: even to handle the dynamics of sqrt(re^2+im^2) located on
2333 : * a scaled unit circle. One need to scale down the results
2334 : * with a factor 2, that is Q_out = Q_in - 1
2335 : * sqrt( -32768.^2 + -32768.^2 ) results in = 23170 Qin-1,
2336 : * which corresponds to 46341 in the Qin domain
2337 : *----------------------------------------------------------------------------*/
2338 :
2339 0 : void fft_spec2_sqrt_approx_fx(const Word16 x[], /* i : Input vector: complex spectrum , Qin */
2340 : Word16 xMagSqrt[], /* o : sqrt of magnitude square spectrum Qout=Qin-1*/
2341 : const Word16 N /* i : Input vector x length */
2342 : )
2343 : {
2344 : Counter i;
2345 : Word16 l;
2346 : const Word16 *pRe, *pIm;
2347 : Word16 * pMagSqrt;
2348 :
2349 : #ifdef DYNMEM_COUNT
2350 : Dyn_Mem_In("fft_spec2_sqrt_approx_fx", sizeof(struct {
2351 : Counter i;
2352 : Word16 l;
2353 : const Word16 *pRe, *pIm;
2354 : Word16 * pMagSqrt;
2355 : }));
2356 : #endif
2357 :
2358 : #ifdef WMOPS
2359 : push_wmops("PhECU::fft_spec2_sqrt_approx_fx");
2360 : #endif
2361 :
2362 : /* Magnitude at 0. only real component */
2363 0 : pMagSqrt = &xMagSqrt[0];
2364 0 : pRe = &x[0];
2365 0 : *pMagSqrt++ = mult(abs_s(*pRe++), C_0p99_Q14); move16();
2366 :
2367 : /* From 1 to (N/2 - 1). */
2368 0 : pIm = &x[N - 1];
2369 :
2370 0 : l = sub(shr_pos_pos(N, 1), 1); /* N/2 - 1. */
2371 0 : l = s_min(l, (LPROT48K_RED / 2 - 1) + DELTA_CORR_F0_INT);
2372 : /* at 48 k the top 8 khz are always zero, and further peaks are not
2373 : located above LPROT48K_RED 32 kHz */
2374 0 : FOR(i = 0; i < l; i++)
2375 : {
2376 0 : *pMagSqrt++ = sqrtMagnApprox_fx(*pRe++, *pIm--); move16();
2377 : }
2378 :
2379 : /* The sqrt magnitude square at N/2 - only real component */
2380 0 : *pMagSqrt = mult(abs_s(*pRe), C_0p99_Q14); move16();
2381 :
2382 : #ifdef DYNMEM_COUNT
2383 : Dyn_Mem_Out();
2384 : #endif
2385 : #ifdef WMOPS
2386 : pop_wmops();
2387 : #endif
2388 0 : return;
2389 : }
2390 :
2391 : Word16
2392 0 : imax2_jacobsen_mag_fx(/* o: The location, relative to the middle of the 3 given data point, of the maximum.
2393 : (Q15) */
2394 : const Word16 *y_re, /* i: The 3 given data points. real part order -1 0 1 */
2395 : const Word16 *y_im, /* i: The 3 given data points. imag part order 1 0 -1 (from FFT)*/
2396 : const Word16
2397 : special /* i: -1 = left edge special case, 0 = normal, +1 = right edge special case */
2398 : )
2399 : {
2400 : Word16 posi;
2401 : Word16 man, expo;
2402 : const Word16 *pY;
2403 : Word16 y_m1_re, y_0_re, y_p1_re;
2404 : Word16 y_m1_im, y_0_im, y_p1_im;
2405 :
2406 : Word16 D_re, D_im, N_re, N_im;
2407 :
2408 : Word32 L_sign, L_denom, L_numer;
2409 :
2410 : #ifdef DYNMEM_COUNT
2411 : Dyn_Mem_In("imax2_jacobsen_mag_fx", sizeof(struct {
2412 : Word16 posi;
2413 : Word16 man, expo;
2414 : const Word16 *pY;
2415 : Word16 y_m1_re, y_0_re, y_p1_re;
2416 : Word16 y_m1_im, y_0_im, y_p1_im;
2417 : Word16 D_re, D_im, N_re, N_im;
2418 : Word32 L_sign, L_denom, L_numer;
2419 : }));
2420 : #endif
2421 :
2422 : #ifdef WMOPS
2423 : push_wmops("PhECU::imax2_jacobsen_mag_fx");
2424 : #endif
2425 :
2426 : /* Jacobsen estimates peak offset relative y_0 using
2427 : * X_m1 - X_p1
2428 : * d = REAL ( ------------------- ) * c_jacob
2429 : * 2*X_0 - X_m1 -Xp1
2430 : *
2431 : * Where c_jacob is a window dependent constant
2432 : */
2433 : #define C_JACOB_Q14 18725 /* c_jacob = 1.1429; % assume 0.1875 hammrect window 'periodic' */
2434 :
2435 0 : ASSERT(special == 0); /* always use other imax for edges cases */
2436 : UNUSED(special);
2437 :
2438 : /* Get the bin parameters into variables */
2439 0 : pY = y_re;
2440 0 : y_m1_re = *pY++;
2441 0 : y_0_re = *pY++;
2442 0 : y_p1_re = *pY++;
2443 :
2444 : /* Same for imaginary parts - note reverse order from FFT */
2445 0 : pY = y_im;
2446 0 : y_p1_im = *pY++;
2447 0 : y_0_im = *pY++;
2448 0 : y_m1_im = *pY++;
2449 :
2450 0 : test();
2451 0 : IF( norm_s(y_0_re) == 0 || norm_s(y_0_im) == 0 )
2452 : {
2453 : #define JACOB_MARGIN 2
2454 : /* for very high peaks the Complex denominator values may need to be downshifted 2 steps */
2455 0 : y_0_re = shr_pos(y_0_re, JACOB_MARGIN);
2456 0 : y_0_im = shr_pos(y_0_im, JACOB_MARGIN);
2457 :
2458 0 : y_m1_re = shr_pos(y_m1_re, JACOB_MARGIN);
2459 0 : y_m1_im = shr_pos(y_m1_im, JACOB_MARGIN);
2460 :
2461 0 : y_p1_re = shr_pos(y_p1_re, JACOB_MARGIN);
2462 0 : y_p1_im = shr_pos(y_p1_im, JACOB_MARGIN);
2463 : }
2464 :
2465 : /* prepare numerator real and imaginary parts*/
2466 0 : N_re = sub(y_m1_re, y_p1_re);
2467 0 : N_im = sub(y_m1_im, y_p1_im);
2468 :
2469 : /* prepare denominator real and imaginary parts */
2470 :
2471 0 : D_re = sub(sub(shl_pos(y_0_re, 1), y_m1_re), y_p1_re);
2472 0 : D_im = sub(sub(shl_pos(y_0_im, 1), y_m1_im), y_p1_im);
2473 :
2474 : /* REAL part of complex division */
2475 0 : L_numer = L_mac0(L_mult0(N_re, D_re), N_im, D_im);
2476 0 : L_denom = L_mac0(L_mult0(D_re, D_re), D_im, D_im);
2477 0 : L_sign = L_xor(L_numer, L_denom); /* Preserve the sign since div_s() only takes positive arguments. */
2478 :
2479 0 : L_numer = L_abs(L_numer);
2480 0 : L_denom = L_abs(L_denom);
2481 :
2482 0 : test();
2483 0 : IF(L_numer != 0 && L_denom != 0)
2484 : {
2485 :
2486 0 : man = plc_phEcu_ratio_fx(L_numer, L_denom, &expo); /* The mantissa is considered in Q15 */
2487 :
2488 0 : man = mult_r(man, C_JACOB_Q14);
2489 0 : posi = shr_sat( man, sub(expo, 2));
2490 : /* to Q15 (and due to saturation, it is automatically bound inside [-1.0,1.0].) */
2491 :
2492 0 : if (L_sign < 0) /* Restore the sign. */
2493 : {
2494 0 : posi = negate(posi);
2495 : }
2496 : }
2497 : ELSE
2498 : {
2499 0 : posi = 0; move16(); /* flat top, division is not possible choose center freq */
2500 : }
2501 :
2502 : #ifdef DYNMEM_COUNT
2503 : Dyn_Mem_Out();
2504 : #endif
2505 : #ifdef WMOPS
2506 : pop_wmops();
2507 : #endif
2508 0 : return posi; /* Q15. The position either left or right relative to the index of the middle of the 3 given
2509 : data points. */
2510 : }
2511 :
2512 : /* Convert 32 Bit FFT into 16 bit fft domain */
2513 0 : static void intlvW32_2_flippedW16(
2514 : Word32 * L_x, /* i : interleaved coeffs DC, Fs/2, Re(1),Im(1), Re(2),Im(2), ... ] */
2515 : const Word16 numPairs, /* i: typically (fft-size/2 -1), re/im coeffs to copy */
2516 : const Word16 Lprot, /* i: fft size , including DC+fs/2 */
2517 : Word16 * x /* o : flipped coeffs , [DC, Re(1),.. Re(N-1/2) , Fs/2, Im(N-1/2) ... Im(1) ] */
2518 : )
2519 : {
2520 :
2521 : /* reorder real and imag components, and apply fractional scale for 24/48Khz */
2522 : Counter m;
2523 : Counter numPairsLocal;
2524 0 : Word32 *pX_L = &L_x[2]; /*ptr init*/
2525 0 : Word16 *pX_Re = &x[1]; /*ptr init*/
2526 0 : Word16 *pX_Im = &x[Lprot - 1]; /*ptr init*/
2527 :
2528 : #define FHG_FFT_UPSHIFT 2
2529 :
2530 : #ifdef WMOPS
2531 : push_wmops("PhECU::intlvW32_2_flippedW16");
2532 : #endif
2533 : #ifdef DYNMEM_COUNT
2534 : Dyn_Mem_In("intlvW32_2_flippedW16", sizeof(struct {
2535 : Counter m;
2536 : Counter numPairsLocal;
2537 : Word32 *pX_L;
2538 : Word16 *pX_Re;
2539 : Word16 *pX_Im;
2540 : }));
2541 : #endif
2542 :
2543 : /* make the scaling of 8/3= 4*0.666 here for 24 kHz and 48 kHz using 16x16 instead or 32x16 ops
2544 : a limited loss of SNR */
2545 :
2546 0 : test();
2547 0 : IF(sub(numPairs, 383) == 0 || sub(numPairs, 191) == 0)
2548 : { /* 24,48 kHz , 16 ms , scale by 8/3 = .666*4 */
2549 0 : numPairsLocal = s_min(numPairs, 383 - 63); /* do not copy bins above 20 kHz */
2550 : /* for 48 kHz is to only go up to 40 kHz pairs , */
2551 0 : FOR(m = 0; m < numPairsLocal; m++)
2552 : {
2553 : /* multiply by (8/3)*(2.^FHG_FFT_UPSHIFT) */
2554 : /* note: multiplication by 2/3 need to preceed upshift , due to FFT scaling being very close to the 48/24 3
2555 : * split kHz bit margin */
2556 0 : *pX_Re++ = extract_h(L_shl_pos(Mpy_32_16_lc3plus(*pX_L++, FEC_TWOTHIRDS_Q15), FHG_FFT_UPSHIFT + 2));
2557 0 : move16();
2558 0 : *pX_Im-- = extract_h(L_shl_pos(Mpy_32_16_lc3plus(*pX_L++, FEC_TWOTHIRDS_Q15), FHG_FFT_UPSHIFT + 2));
2559 0 : move16();
2560 : }
2561 : /* Place the two real only components */
2562 0 : x[0] = extract_h(L_shl_pos(Mpy_32_16_lc3plus(L_x[0], FEC_TWOTHIRDS_Q15), FHG_FFT_UPSHIFT + 2)); /* DC */
2563 0 : move16();
2564 0 : m = shr_pos_pos(Lprot, 1);
2565 0 : x[m] = extract_h(L_shl_pos(Mpy_32_16_lc3plus(L_x[1], FEC_TWOTHIRDS_Q15), FHG_FFT_UPSHIFT + 2)); /* fs/2 */
2566 0 : move16();
2567 : }
2568 : ELSE
2569 : { /* 8,16,32 kHz 16 ms no additional scaling by 8/3 */
2570 :
2571 0 : FOR(m = 0; m < numPairs; m++)
2572 : {
2573 0 : *pX_Re++ = extract_h(L_shl_pos(*pX_L++, FHG_FFT_UPSHIFT));
2574 0 : move16();
2575 0 : *pX_Im-- = extract_h(L_shl_pos(*pX_L++, FHG_FFT_UPSHIFT));
2576 0 : move16();
2577 : }
2578 :
2579 : /* Place the two real only components */
2580 0 : x[0] = extract_h(L_shl_pos(L_x[0], FHG_FFT_UPSHIFT)); /* DC */ move16();
2581 0 : m = shr_pos_pos(Lprot, 1);
2582 0 : x[m] = extract_h(L_shl_pos(L_x[1], FHG_FFT_UPSHIFT)); /* fs/2 */ move16();
2583 : }
2584 :
2585 : #ifdef DYNMEM_COUNT
2586 : Dyn_Mem_Out();
2587 : #endif
2588 : #ifdef WMOPS
2589 : pop_wmops();
2590 : #endif
2591 0 : }
2592 :
2593 0 : static void flippedW16_2_intlvW32(
2594 : Word16 * x, /* i : flipped coeffs , [DC, Re(1),.. Re(N-1/2) , Fs/2, Im(N-1/2) ... Im(1) ] */
2595 : const Word16 numPairs, /* i: typically (fft-size/2 -1), re/im coeffs to copy */
2596 : const Word16 Lprot, /* i: fft size , including DC+fs/2 */
2597 : Word32 * L_x /* o : interleaved coeffs DC, Fs/2, Re(1),Im(1), Re(2),Im(2), ... ] */
2598 : )
2599 : {
2600 : Counter i;
2601 : Counter numPairsLocal;
2602 : Word32 *pX_L;
2603 : Word16 *pX_re, *pX_im;
2604 : #ifdef WMOPS
2605 : push_wmops("PhECU::flippedW16_2_intlvW32");
2606 : #endif
2607 : #ifdef DYNMEM_COUNT
2608 : Dyn_Mem_In("flippedW16_2_intlvW32", sizeof(struct {
2609 : Counter i;
2610 : Counter numPairsLocal;
2611 : Word32 *pX_L;
2612 : Word16 *pX_re, *pX_im;
2613 : }));
2614 : #endif
2615 :
2616 : /* Convert stored 16 bit into 32bit for fft */
2617 : /* Note during save FFT output was left shifted FHG_FFT_UPSHIFT */
2618 : /* this needs to be restored before one calls ifft to avoid overflow */
2619 :
2620 0 : pX_L = &L_x[2]; /*ptr init*/
2621 0 : pX_re = &x[1]; /*ptr init*/
2622 0 : pX_im = &x[Lprot - 1]; /*ptr init*/
2623 :
2624 0 : numPairsLocal = s_min(320, numPairs); /* 48kHz optimization */
2625 0 : FOR(i = 0; i < numPairsLocal; i++)
2626 : {
2627 0 : *pX_L++ = L_shr_pos(L_deposit_h(*pX_re++), FHG_FFT_UPSHIFT); move32();
2628 0 : *pX_L++ = L_shr_pos(L_deposit_h(*pX_im--), FHG_FFT_UPSHIFT); move32();
2629 : }
2630 : /* at 48KHz zero tail 2x63= 126 bins for Word32 IFFT input */
2631 0 : basop_memset(pX_L, 0, sizeof(Word32) * shl_pos(sub(numPairs, numPairsLocal), 1));
2632 :
2633 : /* Place the two real only components */
2634 0 : L_x[0] = L_shr_pos(L_deposit_h(x[0]), FHG_FFT_UPSHIFT); move32();
2635 0 : L_x[1] = L_shr_pos(L_deposit_h(x[Lprot / 2]), FHG_FFT_UPSHIFT); move32();
2636 :
2637 : #ifdef DYNMEM_COUNT
2638 : Dyn_Mem_Out();
2639 : #endif
2640 : #ifdef WMOPS
2641 : pop_wmops();
2642 : #endif
2643 0 : }
2644 :
2645 0 : void get_sin_cosQ10opt(Word16 phase, /* Q10 0..1023 i.e. 1024=2*pi */
2646 : Word16 * ptrSin, /* Q15 */
2647 : Word16 * ptrCos /* Q15 */
2648 : )
2649 : {
2650 : Word16 sign_val, idx, idx2, idx3;
2651 :
2652 : #ifdef DYNMEM_COUNT
2653 : Dyn_Mem_In("get_sin_cosQ10", sizeof(struct { Word16 sign_val, idx, idx2, idx3; }));
2654 : #endif
2655 :
2656 : #ifdef WMOPS
2657 : push_wmops("PhECU::get_sin_cosQ10opt");
2658 : #endif
2659 :
2660 : /* sin table has a range up to pi/2 (256+1)=257 coeffs*/
2661 :
2662 0 : sign_val = shr_pos_pos(phase, 9); /* highest bit is the sinus sign */
2663 0 : idx = s_and(phase, 0x1ff); /* mask away sign */
2664 :
2665 0 : idx2 = sub(idx, 256);
2666 0 : if (idx2 < 0)
2667 : { /*rising sine part */
2668 0 : *ptrSin = sin_quarterQ15_fx[idx]; move16();
2669 : }
2670 :
2671 0 : idx3 = sub(512, idx);
2672 0 : if (idx2 >= 0)
2673 : { /* decaying part, reverse idx */
2674 0 : *ptrSin = sin_quarterQ15_fx[idx3]; move16();
2675 : }
2676 :
2677 0 : if (sign_val != 0)
2678 : {
2679 0 : *ptrSin = negate(*ptrSin); /*no move when inlined , no sat as max in table is 32767 */
2680 : }
2681 :
2682 : /*cos*/
2683 0 : idx = add(phase, 256); /* +pi/2, i.e. move to cos phase */
2684 0 : idx = s_and(idx, 0x3ff); /* wrap on 10 bits limit */
2685 :
2686 0 : sign_val = shr_pos_pos(idx, 9); /* highest bit is the sign */
2687 0 : idx = s_and(idx, 0x1ff); /* mask away sign */
2688 :
2689 0 : idx2 = sub(idx, 256);
2690 0 : if (idx2 < 0)
2691 : { /*rising sine part */
2692 0 : *ptrCos = sin_quarterQ15_fx[idx]; move16();
2693 : }
2694 :
2695 0 : idx3 = sub(512, idx);
2696 0 : if (idx2 >= 0)
2697 : { /* decaying part, reverse idx*/
2698 0 : *ptrCos = sin_quarterQ15_fx[idx3]; move16();
2699 : }
2700 :
2701 0 : if (sign_val != 0)
2702 : {
2703 0 : *ptrCos = negate(*ptrCos); /* no move when inlined , no sat as max in table is 32767 */
2704 : }
2705 :
2706 : #ifdef DYNMEM_COUNT
2707 : Dyn_Mem_Out();
2708 : #endif
2709 : #ifdef WMOPS
2710 : pop_wmops();
2711 : #endif
2712 0 : }
2713 :
2714 0 : static Word16 plc_phEcu_nonpure_tone_ana_fx(const Word16* plocs,
2715 : const Word16 n_plocs,
2716 : const Word16* X, /*i: { DC, Re1,Re2 ,.....ReN , Fs/2, ImN... , Im2, Im1} */
2717 : const Word32 *L_Xavg, /* i : Frequency group amp averages for tonal tilt analysis Max upshifted */
2718 : const Word16 Lprot,
2719 : const Word16 fs_idx)
2720 : {
2721 : #ifdef DYNMEM_COUNT
2722 : Dyn_Mem_In("PhECU::nonpure_tone_ana", sizeof(struct {
2723 : Word16 nonpure_tone_detect;
2724 : Word16 n_ind, tone_ind, low_ind, high_ind;
2725 : Word16 peak_amp, peak_amp2, valley_amp, x_abs[(1 + 2 * ONE_SIDED_SINE_WIDTH + 2 * 1)];
2726 : Word16 sineband_ind_low, sineband_ind_high;
2727 : Word16 i, N_grp;
2728 : Word16 tmp, tmp_dB;
2729 : Word32 L_tot_inc_HF, L_tot_inc_LF;
2730 : Word16* pRe;
2731 : Word16* pIm;
2732 : Word16* gwlpr_fxPlus1;
2733 : Word32 L_nf;
2734 : Word32 L_XavgL2_fx[MAX_LGW + 1];
2735 : Word32 L_XavgUp_fx[MAX_LGW + 1];
2736 : }));
2737 : #endif
2738 :
2739 : #ifdef WMOPS
2740 : push_wmops("PhEcu::nonpure_tone_ana_fx");
2741 : #endif
2742 :
2743 : Word16 nonpure_tone_detect; /* output variable */
2744 : Word16 n_ind, tone_ind, low_ind, high_ind;
2745 : Word16 peak_amp, peak_amp2, valley_amp, x_abs[(1 + 2 * ONE_SIDED_SINE_WIDTH + 2 * 1)];
2746 :
2747 : Word16 sineband_ind_low, sineband_ind_high;
2748 : Word16 i, N_grp;
2749 : Word16 tmp;
2750 :
2751 : Word32 L_Ltot_inc_HF, L_Ltot_inc_LF;
2752 : Word32 L_tmp_dB;
2753 :
2754 : const Word16* pRe;
2755 : const Word16* pIm;
2756 : const Word16* gwlpr_fxPlus1;
2757 : Word32 L_nf;
2758 : Word32 L_XavgL2_fx[MAX_LGW + 1], L_tmp ;
2759 : Word32 L_XavgUp_fx[MAX_LGW + 1]; /* upscaled values , excluding peak band(s) */
2760 :
2761 : /* use of a compressed hearing sensitivity curve allowing more energy deviation in highest and lowest bands */
2762 : /* ROM table Word16 scATHFx[MAX_LGW - 1] */
2763 :
2764 : /* init */
2765 0 : nonpure_tone_detect = 0; move16(); /* Word16 register with decisions */
2766 :
2767 0 : peak_amp = 0; move16();
2768 0 : peak_amp2 = 0; move16();
2769 :
2770 : /* limit single sine optimization to when 2 peaks are close enough to represent a single sinusoid */
2771 0 : test();
2772 0 : if ((sub(n_plocs, 2) == 0) &&
2773 0 : (sub(sub(plocs[1], plocs[0]), ONE_SIDED_SINE_WIDTH) >= 0)
2774 : ) /* NB, plocs is an ordered vector */
2775 : {
2776 0 : nonpure_tone_detect = s_or(nonpure_tone_detect, 0x01);
2777 : }
2778 :
2779 : /* local bin wise dynamics analysis, for 1 or 2 initial local maxima peaks ,
2780 : if 2 peaks , we do the analysis based on the location of the largest abs peak */
2781 : {
2782 0 : tone_ind = 0; move16(); /* one peak only , use position plocs[tone_ind], tone_ind==0 */
2783 :
2784 0 : peak_amp = shr(abs_s(X[0]), 1); /* plocs[0]=0, simply multiply DC*0.5 to match scaling in function sqrtMagnApprox_fx */
2785 0 : test();
2786 0 : IF(plocs[0] != 0)
2787 : {
2788 0 : peak_amp = sqrtMagnApprox_fx(X[plocs[0]], X[sub(Lprot, plocs[0])]);
2789 : }
2790 :
2791 0 : IF(sub(n_plocs, 2) == 0) /* locate abs max peak */
2792 : {
2793 0 : peak_amp2 = sqrtMagnApprox_fx(X[plocs[1]], X[sub(Lprot, plocs[1])]); /* Re-part and Im-part in different ends of array X */
2794 :
2795 0 : tmp = sub(peak_amp, peak_amp2);
2796 0 : tone_ind = lshr(tmp, 15); /* mask out sign bit 0(for peak_amp0), 1( for peak_amp2), lshr --> no sign extension in shift */
2797 0 : peak_amp = s_max(peak_amp, peak_amp2);
2798 : }
2799 :
2800 0 : low_ind = s_max(1, sub(plocs[tone_ind], (ONE_SIDED_SINE_WIDTH + 1))); /* DC is not allowed in the range */
2801 0 : high_ind = s_min(sub(shr(Lprot, 1), 2), add(plocs[tone_ind], (ONE_SIDED_SINE_WIDTH + 1))); /* Fs/2 is not allowed in the range */
2802 :
2803 0 : n_ind = add(sub(high_ind, low_ind), 1);
2804 : /* find lowest amplitude around the assumed main lobe center location */
2805 :
2806 0 : pRe = &(X[low_ind]); /* ptr init */
2807 0 : pIm = &(X[Lprot - low_ind]); /* ptr init*/
2808 0 : FOR(i = 0; i < n_ind; i++)
2809 : {
2810 0 : x_abs[i] = sqrtMagnApprox_fx(*pRe++, *pIm--); move16(); /* x_abs is downscaled by 0.5 in abs(complex) approximation */
2811 : }
2812 :
2813 0 : valley_amp = peak_amp; move16();
2814 0 : FOR(i = 0; i < n_ind; i++)
2815 : {
2816 0 : valley_amp = s_min(x_abs[i], valley_amp);
2817 : }
2818 :
2819 : /* at least a localized amplitude ratio of 16 (24 dB) required to declare a pure noise-free sinusoid */
2820 0 : if (sub(shr(peak_amp, 4), valley_amp) < 0) /* 1/16 */
2821 : {
2822 0 : nonpure_tone_detect = s_or(nonpure_tone_detect, 0x02); /* not a pure tone due to too low local SNR P2V */
2823 : }
2824 : }
2825 :
2826 :
2827 : /* analyze LF/ HF bands energy dynamics vs the assumed single tone band ( for one or two peaks found cases ) */
2828 : {
2829 : /* Xavg , is a vector of rather rough MDCT based band amplitude estimates in perceptually motivated bands. from approx the last 26 ms of synthesis */
2830 :
2831 : /* eval amplitude relations for assumed tonal band vs lower and higher bands */
2832 : /*N_grp = xavg_N_grp[fs_idx];*/ /* { 4 NB , 5 WB , 6 SSWB , 7 SWB, 8 FB }; */
2833 0 : N_grp = add(fs_idx, 4);
2834 0 : assert(fs_idx < 5);
2835 :
2836 :
2837 : /* establish band(s) with assumed sinusoid tone */
2838 : /* if tone freq location is below first MDCT-band definition, use first band as band location anyway */
2839 :
2840 : /* band 0 , 1 , 2 , 3 , ...*/
2841 : /* dct-inds "c" "0"...11, 12...19, 20...35, 36 ... */
2842 : /* gwplr_fx= [ 1 , 12(750Hz), 20(1250Hz) , 36 , ... */ /* lowest lim+1 in gwplr */
2843 :
2844 : /* for-loop BASOP version */
2845 0 : tmp = 0; move16();
2846 0 : gwlpr_fxPlus1 = &(gwlpr_fx[1]); /* ptr init */
2847 0 : FOR(i = 0; i < N_grp; i++)
2848 : {
2849 0 : if (sub(plocs[tone_ind], gwlpr_fxPlus1[i]) >= 0)
2850 : {
2851 0 : tmp = add(i, 1);
2852 : }
2853 : }
2854 :
2855 0 : sineband_ind_low = tmp; move16();
2856 0 : sineband_ind_high = tmp; move16(); /* typically in the same band as low */
2857 :
2858 : /* a single tone may end up on a band border
2859 : , handle case when assumed tone is more or less right in-between two perceptual bands +/- 4*62.5 Hz */
2860 :
2861 0 : test(); logic16();
2862 0 : if ((sineband_ind_high > 0) &&
2863 0 : (sub(sub(plocs[tone_ind], ONE_SIDED_SINE_WIDTH), gwlpr_fx[add(sineband_ind_high, 1)]) >= 0)
2864 : )
2865 : {
2866 0 : sineband_ind_low = sub(sineband_ind_high, 1);
2867 : }
2868 :
2869 0 : logic16();
2870 0 : if ((sub(sineband_ind_low, sub(N_grp, 1)) < 0) &&
2871 0 : (sub(add(plocs[tone_ind], ONE_SIDED_SINE_WIDTH), gwlpr_fx[add(sineband_ind_low, 1)]) >= 0)
2872 : )
2873 : {
2874 0 : sineband_ind_high = add(sineband_ind_low, 1);
2875 : }
2876 : }
2877 :
2878 : /* intraframe(26 ms), weighted LB and HB envelope dynamics/variation analysis */
2879 : /* envelope analysis ,
2880 : require at least two HF or two LF bands in the envelope taper/roll-off analysis , otherwise skip this condition */
2881 :
2882 0 : logic16();
2883 0 : test(); logic16();
2884 0 : IF( (nonpure_tone_detect == 0) &&
2885 : ( (sub( add(sineband_ind_high, 2), N_grp) < 0 ) ||
2886 : (sub(sineband_ind_low, 2+1 ) >= 0)
2887 : )
2888 : )
2889 : {
2890 :
2891 : /* delta taper-off analysis solution, less sensitive to input bandwidth limitation and levels */
2892 :
2893 : /* test Xavg Word16 result above vs a high resolution Word32 L_Xavg */
2894 : /* strong spectral tilt causes HF to be truncated in Word16 */
2895 :
2896 0 : basop_memcpy(L_XavgUp_fx, L_Xavg, N_grp * sizeof(Word32) ) ;
2897 :
2898 : /* first remove all energy from the assumed tonal peak band(s) */
2899 0 : L_XavgUp_fx[sineband_ind_low] = L_min(L_XavgUp_fx[sineband_ind_low], 1); move32();
2900 0 : L_XavgUp_fx[sineband_ind_high] = L_min(L_XavgUp_fx[sineband_ind_high], 1); move32();
2901 :
2902 0 : tmp = getScaleFactor32_0(L_XavgUp_fx, N_grp); /* o: measured headroom in range [0..32], 32 if all L_x[i] == 0 */
2903 0 : if (sub(tmp, 32) == 0)
2904 : {
2905 0 : nonpure_tone_detect = s_or(nonpure_tone_detect, 0x100); /* also set a flag for an all zero L_Xavg coarse spectrum estimate signal */
2906 : }
2907 :
2908 :
2909 : /* add noise floor to L_Xavg before log2 function call */
2910 0 : L_nf = L_shl(1L, sub(tmp, 1));
2911 0 : L_nf = L_max(L_nf, 1L);
2912 0 : for (i = 0; i < N_grp; i++)
2913 : {
2914 0 : L_XavgUp_fx[i] = L_shl_sat(L_XavgUp_fx[i], tmp); move32(); /* maximize precision before actual log2_fx(Word32) calc call */
2915 :
2916 0 : L_tmp = L_XavgUp_fx[i]; move32();
2917 0 : test();
2918 0 : if (L_tmp <= 0)
2919 : {
2920 0 : L_tmp = L_nf; move32();
2921 : }
2922 :
2923 : /* log2(Upshifted Word32) */
2924 : /* maximize precision in BASOP Log2 function */
2925 0 : L_XavgL2_fx[i] = BASOP_Util_Log2_lc3plus(L_tmp); /* L_input 1.0 or lower --> output always negative */
2926 : /* add 31.0 to store as fractional bits of an upscaled positive Word32 integer input ) */
2927 0 : L_XavgL2_fx[i] = L_add(31L << 25, L_XavgL2_fx[i]); /* only diffs added so 31.0 is cancelled out later , in total only values between +/- 2^6 = [-64 ... -64[ are possible */
2928 : } /* band i for L_Xavg calc*/
2929 : /* verify that an assumed clean sine does not have any odd HF content indications by thresholding the accumulated delta rise in LF/HF side lobes */
2930 0 : L_Ltot_inc_HF = 0; move32();
2931 0 : for (i = (sineband_ind_high + 1); i < (N_grp - 1); i++)
2932 : {
2933 0 : L_tmp_dB = 0; move32();
2934 0 : if (L_sub(L_Xavg[i], L_Xavg[i + 1]) < 0) /* only increases are accumulated */
2935 : {
2936 0 : L_tmp_dB = L_sub(L_XavgL2_fx[i + 1], L_XavgL2_fx[i]); /* obtain ratio in log2 domain */
2937 : }
2938 0 : L_tmp = Mpy_32_16_lc3plus(L_tmp_dB, scATHFx[i]); /* scATHFx[i] is the ATH weight for band i and band i+1 */
2939 0 : L_tmp = L_shr_pos(L_tmp, 1); /* Q25 -> Q24 */
2940 :
2941 0 : L_Ltot_inc_HF = L_add(L_Ltot_inc_HF, L_tmp); /* Q24 */
2942 : }
2943 :
2944 : /* verify that an assumed clean sine does not have any odd LF content by thresholding the accumulated LF reverse up tilt */
2945 :
2946 0 : L_Ltot_inc_LF = 0; move32();
2947 0 : tmp = s_max(0, sub(sineband_ind_low, 1)); /* could also be pointer init */
2948 0 : for (i = tmp; i > 0; i--)
2949 : {
2950 0 : L_tmp_dB = 0; move32();
2951 0 : if (L_sub(L_Xavg[i - 1], L_Xavg[i]) > 0) /* only increases are accumulated */
2952 : {
2953 0 : L_tmp_dB = L_sub(L_XavgL2_fx[i - 1], L_XavgL2_fx[i]); /* obtain ratio in log2 domain */
2954 : }
2955 0 : L_tmp = Mpy_32_16_lc3plus(L_tmp_dB, scATHFx[i-1]); /* scATHFx[i-1] is the ATH weight in between band i-1 and band i */
2956 0 : L_tmp = L_shr_pos(L_tmp, 1); /* Q25 -> Q24 */
2957 0 : L_Ltot_inc_LF = L_add(L_Ltot_inc_LF, L_tmp); /* Q24 */
2958 : }
2959 :
2960 0 : if (L_sub(L_Ltot_inc_HF, SIDE_LIM) > 0) /* side_lim=round(Q24x*0.7474) */
2961 : {
2962 : /* 4.5 dB in log2 is 0.7474 */
2963 0 : nonpure_tone_detect = s_or(nonpure_tone_detect, 0x10); /* still not a pure tone, too great HF side increase */
2964 : }
2965 :
2966 0 : if (L_sub(L_Ltot_inc_LF, SIDE_LIM) > 0) /* side_lim=round(Q24x*0.7474) */
2967 : {
2968 : /* 4.5 dB in log2 is 0.7474 */
2969 0 : nonpure_tone_detect = s_or(nonpure_tone_detect, 0x20); /* still not a pure tone, too great HF side increase */
2970 : }
2971 :
2972 : /* verify that an assumed clean sine does not have any odd LF+HF content by thresholding the accumulated LF+HF unexpected tilt */
2973 0 : if (L_sub(L_add_sat(L_Ltot_inc_LF, L_Ltot_inc_HF), LFHF_LIM) > 0) /* side_lim=round(Q24x*0.7474) */
2974 : {
2975 : /* 6.0 dB in log2 is 0.996578428 */
2976 0 : nonpure_tone_detect = s_or(nonpure_tone_detect, 0x40); /* still not a pure tone, too great LF+HF side increase */
2977 : }
2978 :
2979 : } /* bands available*/
2980 :
2981 : #ifdef WMOPS
2982 : pop_wmops();
2983 : #endif
2984 :
2985 0 : return nonpure_tone_detect;
2986 : }
|
