Line data Source code
1 : /******************************************************************************************************
2 :
3 : (C) 2022-2025 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB,
4 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
5 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
6 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
7 : contributors to this repository. All Rights Reserved.
8 :
9 : This software is protected by copyright law and by international treaties.
10 : The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB,
11 : Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,
12 : Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,
13 : Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other
14 : contributors to this repository retain full ownership rights in their respective contributions in
15 : the software. This notice grants no license of any kind, including but not limited to patent
16 : license, nor is any license granted by implication, estoppel or otherwise.
17 :
18 : Contributors are required to enter into the IVAS codec Public Collaboration agreement before making
19 : contributions.
20 :
21 : This software is provided "AS IS", without any express or implied warranties. The software is in the
22 : development stage. It is intended exclusively for experts who have experience with such software and
23 : solely for the purpose of inspection. All implied warranties of non-infringement, merchantability
24 : and fitness for a particular purpose are hereby disclaimed and excluded.
25 :
26 : Any dispute, controversy or claim arising under or in relation to providing this software shall be
27 : submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in
28 : accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and
29 : the United Nations Convention on Contracts on the International Sales of Goods.
30 :
31 : *******************************************************************************************************/
32 :
33 : /*====================================================================================
34 : EVS Codec 3GPP TS26.443 Nov 04, 2021. Version 12.14.0 / 13.10.0 / 14.6.0 / 15.4.0 / 16.3.0
35 : ====================================================================================*/
36 :
37 : #include <stdint.h>
38 : #include "options.h"
39 : #ifdef DEBUGGING
40 : #include "debug.h"
41 : #endif
42 : #include <math.h>
43 : #include "cnst.h"
44 : #include "prot.h"
45 : #include "rom_com.h"
46 : #include "wmc_auto.h"
47 :
48 : /*-------------------------------------------------------------------*
49 : * modify_Fs()
50 : *
51 : * Function for resampling of signals
52 : *-------------------------------------------------------------------*/
53 :
54 : /*! r: length of output */
55 49002992 : int16_t modify_Fs(
56 : const float sigIn[], /* i : signal to decimate */
57 : const int16_t lg, /* i : length of input + 0 delay signaling */
58 : const int32_t fin, /* i : frequency of input */
59 : float sigOut[], /* o : decimated signal */
60 : const int32_t fout, /* i : frequency of output */
61 : float mem[], /* i/o: filter memory */
62 : const int16_t nblp /* i : flag indicating if NB low-pass is applied */
63 : )
64 : {
65 : int16_t i;
66 : int16_t lg_out, fac_num, fac_den, filt_len, frac, mem_len;
67 : float num_den;
68 : int16_t datastep, fracstep;
69 : float *sigIn_ptr;
70 : float signal_tab[3 * L_FILT_MAX + L_FRAME48k], *signal, *signal_ana; /* 3* as 2* for memory and 1* for future prediction */
71 : float A[M + 1], r[M + 1], epsP[M + 1], val;
72 : int16_t mem_len_ana;
73 : int16_t plus_sample_in;
74 : int16_t j;
75 : float mu_preemph;
76 : float mem_preemph;
77 : const Resampling_cfg *cfg_ptr;
78 49002992 : int16_t flag_low_order = 0;
79 : int16_t filt_len_tmp;
80 :
81 : /*-------------------------------------------------------------------*
82 : * IIR filters for resampling to/from 8 kHz
83 : *-------------------------------------------------------------------*/
84 :
85 : /*-------------------------------------------------------------------*
86 : * Find the resampling configuration
87 : *-------------------------------------------------------------------*/
88 :
89 49002992 : if ( fin == fout )
90 : {
91 : /* just copy the signal and quit */
92 0 : for ( i = 0; i < lg; i++ )
93 : {
94 0 : sigOut[i] = sigIn[i];
95 : }
96 :
97 0 : return lg;
98 : }
99 : else
100 : {
101 : /* find the resampling configuration in the lookup table */
102 632588126 : for ( cfg_ptr = &resampling_cfg_tbl[0]; ( cfg_ptr->fin != 0 ) && !( cfg_ptr->fin == fin && cfg_ptr->fout == fout ); cfg_ptr++ )
103 : {
104 : }
105 :
106 :
107 : /* find config with NB 4kHz low-pass */
108 49002992 : if ( nblp && ( fin > 8000 ) && ( fout == 12800 ) )
109 : {
110 0 : flag_low_order = 1;
111 0 : for ( cfg_ptr++; ( cfg_ptr->fin != 0 ) && !( cfg_ptr->fin == fin && cfg_ptr->fout == fout ); cfg_ptr++ )
112 : {
113 : }
114 : }
115 :
116 : #define WMC_TOOL_SKIP
117 : /* Retrieve and/or calculate the resampling parameters */
118 49002992 : fac_num = cfg_ptr->fac_num;
119 49002992 : fac_den = (int16_t) ( ( cfg_ptr->fin * fac_num ) / cfg_ptr->fout );
120 49002992 : lg_out = ( lg * fac_num ) / fac_den;
121 49002992 : filt_len = cfg_ptr->filt_len;
122 :
123 49002992 : mem_len = 2 * filt_len;
124 49002992 : plus_sample_in = 0; /* default, regular delay */
125 49002992 : frac = 0;
126 :
127 49002992 : if ( fin == 8000 && fout == 12800 )
128 : {
129 7710 : plus_sample_in = 7;
130 7710 : frac = 4;
131 : }
132 :
133 49002992 : signal = signal_tab + 2 * L_FILT_MAX + L_FRAME48k - mem_len - lg;
134 49002992 : signal_ana = signal;
135 49002992 : mem_len_ana = mem_len;
136 : #undef WMC_TOOL_SKIP
137 : }
138 :
139 : /*-------------------------------------------------------------------*
140 : * FIR filters for resampling to/from 12.8, 16, 32, 48 kHz
141 : *-------------------------------------------------------------------*/
142 :
143 : /* append filter memory */
144 3951524338 : for ( i = 0; i < 2 * filt_len; i++ )
145 : {
146 3902521346 : signal[i] = mem[i];
147 : }
148 :
149 20905019738 : for ( i = 0; i < lg; i++ )
150 : {
151 20856016746 : signal[i + ( 2 * filt_len )] = sigIn[i];
152 : }
153 :
154 49002992 : if ( plus_sample_in > 0 )
155 : {
156 7710 : autocorr( signal_ana + mem_len_ana + lg - LEN_WIN_SSS, r, 1, LEN_WIN_SSS, wind_sss, 0, 0, 0 );
157 :
158 7710 : mu_preemph = r[1] / r[0];
159 7710 : mem_preemph = signal_ana[mem_len_ana + lg - LEN_WIN_SSS - 1];
160 7710 : preemph( signal_ana + mem_len_ana + lg - LEN_WIN_SSS, mu_preemph, LEN_WIN_SSS, &mem_preemph );
161 :
162 : /* Autocorrelations */
163 7710 : autocorr( signal_ana + mem_len_ana + lg - LEN_WIN_SSS, r, M, LEN_WIN_SSS, wind_sss, 0, 0, 0 );
164 :
165 7710 : lag_wind( r, M, fin, LAGW_STRONG );
166 :
167 : /* Levinson-Durbin */
168 7710 : lev_dur( A, r, M, epsP );
169 :
170 61680 : for ( i = 0; i < plus_sample_in; i++ )
171 : {
172 53970 : val = 0;
173 917490 : for ( j = 1; j <= M; j++ )
174 : {
175 863520 : val -= signal[i + lg + mem_len - j] * A[j];
176 : }
177 53970 : signal[i + lg + mem_len] = val; /* AZ ringing padding */
178 : }
179 :
180 7710 : mem_preemph = signal[mem_len + lg - LEN_WIN_SSS - 1];
181 7710 : deemph( signal + mem_len + lg - LEN_WIN_SSS, mu_preemph, LEN_WIN_SSS + plus_sample_in, &mem_preemph );
182 : }
183 :
184 : /* interpolation */
185 49002992 : datastep = fac_den / fac_num;
186 49002992 : fracstep = fac_den - datastep * fac_num; /* equivalent to datastep = fac_den % fac_num */
187 :
188 49002992 : sigIn_ptr = signal + filt_len + plus_sample_in;
189 :
190 49002992 : filt_len_tmp = filt_len;
191 49002992 : if ( flag_low_order )
192 : {
193 0 : filt_len_tmp = ( filt_len + 1 ) >> 1;
194 : }
195 :
196 6904378206 : for ( i = 0; i < lg_out; i++ )
197 : {
198 6855375214 : sigOut[i] = interpolation( sigIn_ptr, cfg_ptr->filter, frac, fac_num, filt_len_tmp );
199 :
200 6855375214 : frac = frac + fracstep;
201 6855375214 : if ( frac >= fac_num )
202 : {
203 3790070850 : frac = frac - fac_num;
204 3790070850 : sigIn_ptr++;
205 : }
206 :
207 6855375214 : sigIn_ptr += datastep;
208 : }
209 :
210 : /* rescaling */
211 49002992 : if ( ( fac_num > fac_den ) == ( ( cfg_ptr->flags & RS_INV_FAC ) != 0 ) )
212 : {
213 48954635 : num_den = (float) fac_num / fac_den;
214 :
215 6861935009 : for ( i = 0; i < lg_out; i++ )
216 : {
217 6812980374 : sigOut[i] *= num_den;
218 : }
219 : }
220 :
221 : /* update the filter memory */
222 3951524338 : for ( i = 0; i < 2 * filt_len; i++ )
223 : {
224 3902521346 : mem[i] = signal[i + lg];
225 : }
226 :
227 49002992 : return lg_out;
228 : }
229 :
230 : /*-------------------------------------------------------------------*
231 : * modify_Fs_intcub3m_sup()
232 : *
233 : *
234 : *-------------------------------------------------------------------*/
235 :
236 : /*! r: length of output */
237 322 : int16_t modify_Fs_intcub3m_sup(
238 : const float sigIn[], /* i : signal to decimate with memory of 2 samples (indexes -2 & -1) */
239 : const int16_t lg, /* i : length of input (suppose that lg is such that lg_out is integer, ex multiple of 5 in case of 16kHz to 12.8 kHz) */
240 : const int32_t fin, /* i : frequency of input */
241 : float sigOut[], /* o : decimated signal */
242 : const int32_t fout, /* i : frequency of output */
243 : int16_t *delayout /* o : delay of output */
244 : )
245 : {
246 : int16_t i, j, k, i1, i2, k1, k2, k3, kk, cind;
247 : int16_t lg_out, fk1, k2d, k3d;
248 : float cc[4][4];
249 : float vv;
250 322 : const double( *cu )[3] = 0;
251 :
252 : /*-------------------------------------------------------------------*
253 : * Find the resampling configuration
254 : *-------------------------------------------------------------------*/
255 :
256 : /* check if fin and fout are the same */
257 322 : if ( fin == fout )
258 : {
259 : /* just copy the signal and quit */
260 0 : for ( i = 0; i < lg; i++ )
261 : {
262 0 : sigOut[i] = sigIn[i];
263 : }
264 :
265 0 : *delayout = 0;
266 0 : return lg;
267 : }
268 : else
269 : {
270 : /* length of the interpolated signal */
271 322 : lg_out = (int16_t) ( lg * fout / fin );
272 :
273 : /* cc[x][3]*s*s*s + cc[x][2]*s*s + cc[x][1]*s + cc[x][0]; indexes relatives of s : -1 0 1 2 */
274 : /* d : cc[x][0] = s[0] */
275 : /* b : cc[x][2] =(s[-1]+s[1])/2-s[0] */
276 : /* a : cc[x][3] = (s[-1]+s[2]-s[0]-s[1]-4*cc[x][2]) / 6 */
277 : /* c : cc[x][1] = s[1]-s[0]-cc[x][3]-cc[x][2] */
278 :
279 : /* coef inits using memory (indexes < 0) */
280 : /* cc[2][] : indexes -2 -1 0 1 */
281 322 : cc[2][0] = sigIn[-1] / 3;
282 322 : cc[2][2] = ( sigIn[-2] + sigIn[0] ) / 2 - sigIn[-1];
283 322 : cc[2][3] = ( sigIn[-2] + sigIn[1] - sigIn[-1] - sigIn[0] - 4 * cc[2][2] ) / 6;
284 322 : cc[2][1] = sigIn[0] - sigIn[-1] - cc[2][3] - cc[2][2];
285 :
286 : /* cc[3][] : indexes -1 0 1 2 */
287 322 : cc[3][0] = sigIn[0] / 3;
288 322 : cc[3][2] = ( sigIn[-1] + sigIn[1] ) / 2 - sigIn[0];
289 322 : cc[3][3] = ( sigIn[-1] + sigIn[2] - sigIn[0] - sigIn[1] - 4 * cc[3][2] ) / 6;
290 322 : cc[3][1] = sigIn[1] - sigIn[0] - cc[3][3] - cc[3][2];
291 322 : j = 0;
292 :
293 322 : if ( fin == 12800 )
294 : {
295 67 : if ( fout == 8000 )
296 : {
297 0 : cind = 0;
298 : }
299 67 : else if ( fout == 16000 )
300 : {
301 2 : cind = 1;
302 : }
303 65 : else if ( fout == 32000 )
304 : {
305 60 : cind = 2;
306 : }
307 5 : else if ( fout == 48000 )
308 : {
309 5 : cind = 3;
310 : }
311 : else
312 : {
313 0 : printf( "warning, output sampling frequency %d not implemented for input %d", fout, fin );
314 0 : return ( -1 );
315 : }
316 : }
317 255 : else if ( fin == 16000 )
318 : {
319 255 : if ( fout == 12800 )
320 : {
321 0 : cind = 4;
322 : }
323 255 : else if ( fout == 32000 )
324 : {
325 58 : cind = 5;
326 : }
327 197 : else if ( fout == 48000 )
328 : {
329 197 : cind = 6;
330 : }
331 : else
332 : {
333 0 : printf( "warning, output sampling frequency %d not implemented for input %d", fout, fin );
334 0 : return ( -1 );
335 : }
336 : }
337 : else
338 : {
339 0 : printf( "warning, input sampling frequency %d not implemented", fin );
340 0 : return ( -1 );
341 : }
342 :
343 322 : *delayout = ct2[cind][9];
344 :
345 322 : if ( ct2[cind][12] == 15 )
346 : {
347 264 : cu = cu15;
348 : }
349 :
350 322 : if ( ct2[cind][12] == 4 )
351 : {
352 58 : cu = cu4;
353 : }
354 :
355 322 : fk1 = 2 * ct2[cind][12] - 2;
356 322 : k2d = fk1 / 2; /* shift of index in cu with respect to the next sample (ex 1.25 -> 0.25 ) */
357 322 : k3d = fk1 - 1; /* to compurte index in cu with respect to the last sample with - sign (ex 1.25 -> -0.75 ) */
358 :
359 322 : kk = 0;
360 4783 : for ( i = 0; i < lg - ct2[cind][11]; )
361 : {
362 4461 : sigOut[j++] = sigIn[i];
363 9345 : for ( k = 0; k < ct2[cind][10]; k++ )
364 : {
365 4884 : cc[kk][0] = sigIn[i + 1] / 3;
366 4884 : cc[kk][2] = ( sigIn[i] + sigIn[i + 2] ) / 2 - sigIn[i + 1];
367 4884 : cc[kk][3] = ( sigIn[i] + sigIn[i + 3] - sigIn[i + 1] - sigIn[i + 2] - 4 * cc[kk][2] ) / 6;
368 4884 : cc[kk][1] = sigIn[i + 2] - sigIn[i + 1] - cc[kk][3] - cc[kk][2];
369 4884 : i++;
370 :
371 4884 : i2 = kk - 2;
372 4884 : i1 = kk - 1;
373 4884 : if ( i1 < 0 )
374 : {
375 1221 : i1 += 4;
376 : }
377 :
378 4884 : if ( i2 < 0 )
379 : {
380 2442 : i2 += 4;
381 : }
382 :
383 13790 : for ( k1 = ct2[cind][k]; k1 < fk1; k1 += ct2[cind][8] )
384 : {
385 8906 : k2 = k1 - k2d;
386 8906 : k3 = k3d - k1;
387 8906 : vv = (float) ( cu[k1][2] * cc[i2][3] + cu[k1][1] * cc[i2][2] + cu[k1][0] * cc[i2][1] + cc[i2][0] );
388 8906 : vv += (float) ( cu[k2][2] * cc[i1][3] + cu[k2][1] * cc[i1][2] + cu[k2][0] * cc[i1][1] + cc[i1][0] );
389 8906 : vv += (float) ( -cu[k3][2] * cc[kk][3] + cu[k3][1] * cc[kk][2] - cu[k3][0] * cc[kk][1] + cc[kk][0] );
390 8906 : sigOut[j++] = vv;
391 : }
392 :
393 4884 : kk++;
394 4884 : if ( kk == 4 )
395 : {
396 1221 : kk = 0;
397 : }
398 : }
399 : }
400 :
401 322 : sigOut[j++] = sigIn[i];
402 :
403 644 : for ( k = 0; k < ct2[cind][11] - 3; k++ )
404 : {
405 322 : cc[kk][0] = sigIn[i + 1] / 3;
406 322 : cc[kk][2] = ( sigIn[i] + sigIn[i + 2] ) / 2 - sigIn[i + 1];
407 322 : cc[kk][3] = ( sigIn[i] + sigIn[i + 3] - sigIn[i + 1] - sigIn[i + 2] - 4 * cc[kk][2] ) / 6;
408 322 : cc[kk][1] = sigIn[i + 2] - sigIn[i + 1] - cc[kk][3] - cc[kk][2];
409 322 : i++;
410 :
411 322 : i2 = kk - 2;
412 322 : i1 = kk - 1;
413 322 : if ( i1 < 0 )
414 : {
415 322 : i1 += 4;
416 : }
417 :
418 322 : if ( i2 < 0 )
419 : {
420 322 : i2 += 4;
421 : }
422 :
423 911 : for ( k1 = ct2[cind][k]; k1 < fk1; k1 += ct2[cind][8] )
424 : {
425 589 : k2 = k1 - k2d;
426 589 : k3 = k3d - k1;
427 589 : vv = (float) ( cu[k1][2] * cc[i2][3] + cu[k1][1] * cc[i2][2] + cu[k1][0] * cc[i2][1] + cc[i2][0] );
428 589 : vv += (float) ( cu[k2][2] * cc[i1][3] + cu[k2][1] * cc[i1][2] + cu[k2][0] * cc[i1][1] + cc[i1][0] );
429 589 : vv += (float) ( -cu[k3][2] * cc[kk][3] + cu[k3][1] * cc[kk][2] - cu[k3][0] * cc[kk][1] + cc[kk][0] );
430 589 : sigOut[j++] = vv;
431 : }
432 :
433 322 : kk++;
434 :
435 322 : if ( kk == 4 )
436 : {
437 0 : kk = 0;
438 : }
439 : }
440 :
441 322 : kk--;
442 322 : if ( kk == -1 )
443 : {
444 0 : kk = 3;
445 : }
446 :
447 322 : if ( ct2[cind][10] == 1 )
448 : {
449 255 : sigOut[j++] = sigIn[i];
450 : }
451 :
452 916 : for ( k1 = ct2[cind][k]; k1 < fk1; k1 += ct2[cind][8] )
453 : {
454 594 : k2 = k1 - k2d;
455 594 : vv = (float) ( cu[k2][2] * cc[kk][3] + cu[k2][1] * cc[kk][2] + cu[k2][0] * cc[kk][1] + cc[kk][0] );
456 594 : sigOut[j++] = vv * 3;
457 : }
458 :
459 322 : if ( ct2[cind][10] < 3 )
460 : {
461 315 : sigOut[j++] = sigIn[i + 1];
462 : }
463 :
464 916 : for ( k1 = ct2[cind][k + 1]; k1 < fk1; k1 += ct2[cind][8] )
465 : {
466 594 : vv = (float) ( cu[k1][2] * cc[kk][3] + cu[k1][1] * cc[kk][2] + cu[k1][0] * cc[kk][1] + cc[kk][0] );
467 594 : sigOut[j++] = vv * 3;
468 : }
469 :
470 322 : if ( ct2[cind][10] == 1 )
471 : {
472 255 : sigOut[j++] = sigIn[i + 2];
473 : }
474 : }
475 :
476 322 : return lg_out;
477 : }
478 :
479 : /*-------------------------------------------------------------------*
480 : * Interpolate_allpass_steep()
481 : *
482 : * Interpolation by a factor 2
483 : *-------------------------------------------------------------------*/
484 :
485 3055387 : void Interpolate_allpass_steep(
486 : const float *in, /* i : input array of size N */
487 : float *mem, /* i/o: memory */
488 : const int16_t N, /* i : number of input samples */
489 : float *out /* o : output array of size 2*N */
490 : )
491 : {
492 : int16_t n, k;
493 : float temp[ALLPASSSECTIONS_STEEP - 1];
494 :
495 : /* upper allpass filter chain */
496 775361767 : for ( k = 0; k < N; k++ )
497 : {
498 772306380 : temp[0] = mem[0] + AP2_STEEP[0] * in[k];
499 772306380 : mem[0] = in[k] - AP2_STEEP[0] * temp[0];
500 :
501 : /* for better performance, unroll this loop */
502 1544612760 : for ( n = 1; n < ALLPASSSECTIONS_STEEP - 1; n++ )
503 : {
504 772306380 : temp[n] = mem[n] + AP2_STEEP[n] * temp[n - 1];
505 772306380 : mem[n] = temp[n - 1] - AP2_STEEP[n] * temp[n];
506 : }
507 :
508 772306380 : out[2 * k + 1] = mem[ALLPASSSECTIONS_STEEP - 1] + AP2_STEEP[ALLPASSSECTIONS_STEEP - 1] * temp[ALLPASSSECTIONS_STEEP - 2];
509 772306380 : mem[ALLPASSSECTIONS_STEEP - 1] = temp[ALLPASSSECTIONS_STEEP - 2] - AP2_STEEP[ALLPASSSECTIONS_STEEP - 1] * out[2 * k + 1];
510 : }
511 :
512 : /* lower allpass filter chain */
513 775361767 : for ( k = 0; k < N; k++ )
514 : {
515 772306380 : temp[0] = mem[ALLPASSSECTIONS_STEEP] + AP1_STEEP[0] * in[k];
516 772306380 : mem[ALLPASSSECTIONS_STEEP] = in[k] - AP1_STEEP[0] * temp[0];
517 :
518 : /* for better performance, unroll this loop */
519 1544612760 : for ( n = 1; n < ALLPASSSECTIONS_STEEP - 1; n++ )
520 : {
521 772306380 : temp[n] = mem[ALLPASSSECTIONS_STEEP + n] + AP1_STEEP[n] * temp[n - 1];
522 772306380 : mem[ALLPASSSECTIONS_STEEP + n] = temp[n - 1] - AP1_STEEP[n] * temp[n];
523 : }
524 :
525 772306380 : out[2 * k] = mem[2 * ALLPASSSECTIONS_STEEP - 1] + AP1_STEEP[ALLPASSSECTIONS_STEEP - 1] * temp[ALLPASSSECTIONS_STEEP - 2];
526 772306380 : mem[2 * ALLPASSSECTIONS_STEEP - 1] = temp[ALLPASSSECTIONS_STEEP - 2] - AP1_STEEP[ALLPASSSECTIONS_STEEP - 1] * out[2 * k];
527 : }
528 :
529 3055387 : return;
530 : }
531 :
532 : /*-------------------------------------------------------------------*
533 : * Decimate_allpass_steep()
534 : *
535 : * Decimation by a factor 2
536 : *-------------------------------------------------------------------*/
537 :
538 8534732 : void Decimate_allpass_steep(
539 : const float *in, /* i : input array of size N */
540 : float *mem, /* i/o: memory */
541 : const int16_t N, /* i : number of input samples */
542 : float *out /* o : output array of size N/2 */
543 : )
544 : {
545 : int16_t n, k;
546 : float temp[ALLPASSSECTIONS_STEEP];
547 :
548 : /* upper allpass filter chain */
549 2175211797 : for ( k = 0; k < N / 2; k++ )
550 : {
551 2166677065 : temp[0] = mem[0] + AP1_STEEP[0] * in[2 * k];
552 2166677065 : mem[0] = in[2 * k] - AP1_STEEP[0] * temp[0];
553 :
554 : /* for better performance, unroll this loop */
555 4333354130 : for ( n = 1; n < ALLPASSSECTIONS_STEEP - 1; n++ )
556 : {
557 2166677065 : temp[n] = mem[n] + AP1_STEEP[n] * temp[n - 1];
558 2166677065 : if ( fabs( temp[n] ) < 1e-12 )
559 : {
560 46430725 : temp[n] = sign( temp[n] ) * 1e-12f;
561 : }
562 2166677065 : mem[n] = temp[n - 1] - AP1_STEEP[n] * temp[n];
563 : }
564 :
565 2166677065 : out[k] = mem[ALLPASSSECTIONS_STEEP - 1] + AP1_STEEP[ALLPASSSECTIONS_STEEP - 1] * temp[ALLPASSSECTIONS_STEEP - 2];
566 2166677065 : mem[ALLPASSSECTIONS_STEEP - 1] = temp[ALLPASSSECTIONS_STEEP - 2] - AP1_STEEP[ALLPASSSECTIONS_STEEP - 1] * out[k];
567 : }
568 :
569 : /* lower allpass filter chain */
570 8534732 : temp[0] = mem[ALLPASSSECTIONS_STEEP] + AP2_STEEP[0] * mem[2 * ALLPASSSECTIONS_STEEP];
571 8534732 : mem[ALLPASSSECTIONS_STEEP] = mem[2 * ALLPASSSECTIONS_STEEP] - AP2_STEEP[0] * temp[0];
572 :
573 : /* for better performance, unroll this loop */
574 17069464 : for ( n = 1; n < ALLPASSSECTIONS_STEEP - 1; n++ )
575 : {
576 8534732 : temp[n] = mem[ALLPASSSECTIONS_STEEP + n] + AP2_STEEP[n] * temp[n - 1];
577 8534732 : if ( fabs( temp[n] ) < 1e-12 )
578 : {
579 2930470 : temp[n] = sign( temp[n] ) * 1e-12f;
580 : }
581 8534732 : mem[ALLPASSSECTIONS_STEEP + n] = temp[n - 1] - AP2_STEEP[n] * temp[n];
582 : }
583 :
584 8534732 : temp[ALLPASSSECTIONS_STEEP - 1] = mem[2 * ALLPASSSECTIONS_STEEP - 1] + AP2_STEEP[ALLPASSSECTIONS_STEEP - 1] *
585 8534732 : temp[ALLPASSSECTIONS_STEEP - 2];
586 :
587 8534732 : mem[2 * ALLPASSSECTIONS_STEEP - 1] = temp[ALLPASSSECTIONS_STEEP - 2] - AP2_STEEP[ALLPASSSECTIONS_STEEP - 1] *
588 8534732 : temp[ALLPASSSECTIONS_STEEP - 1];
589 8534732 : out[0] = (float) ( ( out[0] + temp[ALLPASSSECTIONS_STEEP - 1] ) * 0.5 );
590 :
591 2166677065 : for ( k = 1; k < N / 2; k++ )
592 : {
593 2158142333 : temp[0] = mem[ALLPASSSECTIONS_STEEP] + AP2_STEEP[0] * in[2 * k - 1];
594 2158142333 : mem[ALLPASSSECTIONS_STEEP] = in[2 * k - 1] - AP2_STEEP[0] * temp[0];
595 :
596 : /* for better performance, unroll this loop */
597 4316284666 : for ( n = 1; n < ALLPASSSECTIONS_STEEP - 1; n++ )
598 : {
599 2158142333 : temp[n] = mem[ALLPASSSECTIONS_STEEP + n] + AP2_STEEP[n] * temp[n - 1];
600 2158142333 : if ( fabs( temp[n] ) < 1e-12 )
601 : {
602 45355156 : temp[n] = sign( temp[n] ) * 1e-12f;
603 : }
604 2158142333 : mem[ALLPASSSECTIONS_STEEP + n] = temp[n - 1] - AP2_STEEP[n] * temp[n];
605 : }
606 :
607 2158142333 : temp[ALLPASSSECTIONS_STEEP - 1] = mem[2 * ALLPASSSECTIONS_STEEP - 1] + AP2_STEEP[ALLPASSSECTIONS_STEEP - 1] *
608 2158142333 : temp[ALLPASSSECTIONS_STEEP - 2];
609 2158142333 : mem[2 * ALLPASSSECTIONS_STEEP - 1] = temp[ALLPASSSECTIONS_STEEP - 2] - AP2_STEEP[ALLPASSSECTIONS_STEEP - 1] *
610 2158142333 : temp[ALLPASSSECTIONS_STEEP - 1];
611 2158142333 : out[k] = (float) ( ( out[k] + temp[ALLPASSSECTIONS_STEEP - 1] ) * 0.5 );
612 : }
613 :
614 : /* z^(-1) */
615 8534732 : mem[2 * ALLPASSSECTIONS_STEEP] = in[N - 1];
616 :
617 8534732 : return;
618 : }
619 :
620 : /*-------------------------------------------------------------------*
621 : * interpolate_3_over_2_allpass()
622 : *
623 : * Interpolate 3/2 using allpass iir polyphase filter. Delay 4 samples @48k
624 : *-------------------------------------------------------------------*/
625 :
626 1351966 : void interpolate_3_over_2_allpass(
627 : const float *input, /* i : input signal */
628 : const int16_t len, /* i : number of input samples */
629 : float *out, /* o : output signal */
630 : float *mem /* i/o: memory */
631 : )
632 : {
633 : int16_t i, loop_len;
634 : float Vu[2], Vm[2], Vl[2]; /* Outputs of three cascaded allpass stages (upper, middle, and lower) */
635 : float out1_buff[L_FRAME32k * 3];
636 : float *out1;
637 : float mem_temp;
638 1351966 : const float *filt_coeff = allpass_poles_3_ov_2;
639 :
640 1351966 : out1 = out1_buff;
641 :
642 830070318 : for ( i = 0; i < len; i++ )
643 : {
644 : /* Upper branch */
645 828718352 : Vu[0] = mem[0] + filt_coeff[0] * ( input[i] - mem[1] );
646 828718352 : Vu[1] = mem[1] + filt_coeff[1] * ( Vu[0] - mem[2] );
647 828718352 : mem[3] = mem[2] + filt_coeff[2] * ( Vu[1] - mem[3] );
648 :
649 828718352 : mem[1] = Vu[0];
650 828718352 : mem[2] = Vu[1];
651 828718352 : *out1++ = mem[3];
652 :
653 : /* Middle branch */
654 828718352 : Vm[0] = mem[0] + filt_coeff[3] * ( input[i] - mem[4] );
655 828718352 : Vm[1] = mem[4] + filt_coeff[4] * ( Vm[0] - mem[5] );
656 828718352 : mem[6] = mem[5] + filt_coeff[5] * ( Vm[1] - mem[6] );
657 :
658 828718352 : mem[4] = Vm[0];
659 828718352 : mem[5] = Vm[1];
660 828718352 : *out1++ = mem[6];
661 :
662 : /* Lower branch */
663 828718352 : Vl[0] = mem[0] + filt_coeff[6] * ( input[i] - mem[7] );
664 828718352 : Vl[1] = mem[7] + filt_coeff[7] * ( Vl[0] - mem[8] );
665 828718352 : mem[9] = mem[8] + filt_coeff[8] * ( Vl[1] - mem[9] );
666 :
667 828718352 : mem[0] = input[i];
668 828718352 : mem[7] = Vl[0];
669 828718352 : mem[8] = Vl[1];
670 828718352 : *out1++ = mem[9];
671 : }
672 :
673 1351966 : loop_len = len * 3 / 2;
674 :
675 : /*decimate by 2 and LPF*/
676 1244429494 : for ( i = 0; i < loop_len; i++ )
677 : {
678 1243077528 : mem_temp = out1_buff[2 * i];
679 1243077528 : out[i] = ( ( ( 0.0473147f ) * ( mem_temp + mem[10] ) ) + ( ( -0.151521f ) * ( mem[11] + mem[14] ) ) );
680 1243077528 : out[i] = ( out[i] + ( ( 0.614152f ) * ( mem[12] + mem[13] ) ) );
681 1243077528 : mem[10] = mem[11];
682 1243077528 : mem[11] = mem[12];
683 1243077528 : mem[12] = mem[13];
684 1243077528 : mem[13] = mem[14];
685 1243077528 : mem[14] = mem_temp;
686 : }
687 :
688 1351966 : return;
689 : }
690 :
691 : /*-------------------------------------------------------------------*
692 : * decimate_2_over_3_allpass()
693 : *
694 : * Decimate 2/3 using allpass iir polyphase filter.
695 : *-------------------------------------------------------------------*/
696 :
697 1435438 : void decimate_2_over_3_allpass(
698 : const float *input, /* i : input signal */
699 : const int16_t len, /* i : number of input samples */
700 : float *out, /* o : output signal */
701 : float *mem, /* i/o: memory */
702 : float *lp_mem /* i/o: memory */
703 : )
704 : {
705 : int16_t i, loop_len;
706 : float Vu[2], Vm[2], Vl[2]; /* Outputs of three cascaded allpass stages (upper, middle, and lower) */
707 : float *out1;
708 : float *in;
709 : float out1_buff[L_FRAME48k * 2];
710 : float tmp;
711 1435438 : const float *filt_coeff = allpass_poles_3_ov_2;
712 1435438 : const float *lp_num = decimate_3_ov_2_lowpass_num;
713 1435438 : const float *lp_den = decimate_3_ov_2_lowpass_den;
714 :
715 : /* Combine the 2nd order iir lpf with the decimation by 2 to improve the efficiency*/
716 1435438 : out1 = out1_buff;
717 :
718 1435438 : *out1++ = lp_num[0] * ( input[0] + lp_mem[0] ) - lp_den[2] * lp_mem[2];
719 1435438 : *out1++ = lp_num[1] * input[0] - lp_den[2] * lp_mem[1];
720 :
721 1069611300 : for ( i = 1; i < len; i++ )
722 : {
723 1068175862 : tmp = lp_num[0] * ( input[i] + input[i - 1] ) - lp_den[2] * out1[-2];
724 1068175862 : *out1++ = tmp;
725 1068175862 : tmp = lp_num[1] * input[i] - lp_den[2] * out1[-2];
726 1068175862 : *out1++ = tmp;
727 : }
728 1435438 : lp_mem[0] = input[len - 1];
729 1435438 : lp_mem[1] = out1[-1];
730 1435438 : lp_mem[2] = out1[-2];
731 :
732 : /* do the all pass polyphase filter with pi/3 cutoff */
733 1435438 : out1 = out;
734 1435438 : in = out1_buff;
735 1435438 : loop_len = (int16_t) len * 2 / 3;
736 :
737 714509638 : for ( i = 0; i < loop_len; i++ )
738 : {
739 : /* Lower branch */
740 713074200 : Vl[0] = mem[8] + filt_coeff[6] * ( *in - mem[9] );
741 713074200 : Vl[1] = mem[9] + filt_coeff[7] * ( Vl[0] - mem[10] );
742 713074200 : mem[11] = mem[10] + filt_coeff[8] * ( Vl[1] - mem[11] );
743 :
744 713074200 : mem[8] = *in++;
745 713074200 : mem[9] = Vl[0];
746 713074200 : mem[10] = Vl[1];
747 713074200 : *out1 = mem[11];
748 :
749 : /* Middle branch */
750 713074200 : Vm[0] = mem[4] + filt_coeff[3] * ( *in - mem[5] );
751 713074200 : Vm[1] = mem[5] + filt_coeff[4] * ( Vm[0] - mem[6] );
752 713074200 : mem[7] = mem[6] + filt_coeff[5] * ( Vm[1] - mem[7] );
753 :
754 713074200 : mem[4] = *in++;
755 713074200 : mem[5] = Vm[0];
756 713074200 : mem[6] = Vm[1];
757 713074200 : *out1 += mem[7];
758 :
759 : /* Upper branch */
760 713074200 : Vu[0] = mem[0] + filt_coeff[0] * ( *in - mem[1] );
761 713074200 : Vu[1] = mem[1] + filt_coeff[1] * ( Vu[0] - mem[2] );
762 713074200 : mem[3] = mem[2] + filt_coeff[2] * ( Vu[1] - mem[3] );
763 :
764 713074200 : mem[0] = *in++;
765 713074200 : mem[1] = Vu[0];
766 713074200 : mem[2] = Vu[1];
767 713074200 : *out1++ += mem[3];
768 : }
769 :
770 1435438 : return;
771 : }
772 :
773 : /*-------------------------------------------------------------------*
774 : * interpolate_3_over_1_allpass()
775 : *
776 : * Interpolate 3/1 using allpass iir polyphase filter. Delay 4 samples @48k
777 : *-------------------------------------------------------------------*/
778 :
779 1244016 : void interpolate_3_over_1_allpass(
780 : const float *input, /* i : input signal */
781 : const int16_t len, /* i : number of input samples */
782 : float *out, /* o : output signal */
783 : float *mem /* i/o: memory */
784 : )
785 : {
786 : int16_t i;
787 : float Vu[2], Vm[2], Vl[2]; /* Outputs of three cascaded allpass stages (upper, middle, and lower) */
788 : float *out1;
789 : float mem_temp;
790 1244016 : const float *filt_coeff = allpass_poles_3_ov_2;
791 :
792 1244016 : out1 = &out[0];
793 :
794 108801456 : for ( i = 0; i < len; i++ )
795 : {
796 : /* Upper branch */
797 107557440 : Vu[0] = mem[0] + filt_coeff[0] * ( input[i] - mem[1] );
798 107557440 : Vu[1] = mem[1] + filt_coeff[1] * ( Vu[0] - mem[2] );
799 107557440 : mem[3] = mem[2] + filt_coeff[2] * ( Vu[1] - mem[3] );
800 :
801 107557440 : mem[1] = Vu[0];
802 107557440 : mem[2] = Vu[1];
803 107557440 : *out1++ = mem[3];
804 :
805 : /* Middle branch */
806 107557440 : Vm[0] = mem[0] + filt_coeff[3] * ( input[i] - mem[4] );
807 107557440 : Vm[1] = mem[4] + filt_coeff[4] * ( Vm[0] - mem[5] );
808 107557440 : mem[6] = mem[5] + filt_coeff[5] * ( Vm[1] - mem[6] );
809 :
810 107557440 : mem[4] = Vm[0];
811 107557440 : mem[5] = Vm[1];
812 107557440 : *out1++ = mem[6];
813 :
814 : /* Lower branch */
815 107557440 : Vl[0] = mem[0] + filt_coeff[6] * ( input[i] - mem[7] );
816 107557440 : Vl[1] = mem[7] + filt_coeff[7] * ( Vl[0] - mem[8] );
817 107557440 : mem[9] = mem[8] + filt_coeff[8] * ( Vl[1] - mem[9] );
818 :
819 107557440 : mem[0] = input[i];
820 107557440 : mem[7] = Vl[0];
821 107557440 : mem[8] = Vl[1];
822 107557440 : *out1++ = mem[9];
823 : }
824 :
825 : /*LPF*/
826 323916336 : for ( i = 0; i < len * 3; i++ )
827 : {
828 322672320 : mem_temp = out[i];
829 322672320 : out[i] = ( ( ( 0.572769f ) * ( mem[12] + mem[11] ) ) - ( ( 0.074005f ) * ( mem_temp + mem[10] ) ) );
830 322672320 : mem[10] = mem[11];
831 322672320 : mem[11] = mem[12];
832 322672320 : mem[12] = mem_temp;
833 : }
834 :
835 1244016 : return;
836 : }
837 :
838 :
839 : /*-------------------------------------------------------------------*
840 : * retro_interp4_5()
841 : *
842 : *
843 : *-------------------------------------------------------------------*/
844 :
845 11388 : void retro_interp4_5(
846 : const float *syn,
847 : float *pst_old_syn )
848 : {
849 : float *pf5, *pf4;
850 : int16_t c;
851 :
852 : /* resample st->pst_old_syn in a reverse way to preserve time-alignment */
853 11388 : pf4 = (float *) &pst_old_syn[58];
854 11388 : pf5 = (float *) pst_old_syn;
855 660504 : for ( c = 0; c < 57; c++ )
856 : {
857 649116 : *pf5++ = pf4[0];
858 649116 : *pf5++ = 0.2f * pf4[0] + 0.8f * pf4[1];
859 649116 : *pf5++ = 0.4f * pf4[1] + 0.6f * pf4[2];
860 649116 : *pf5++ = 0.6f * pf4[2] + 0.4f * pf4[3];
861 649116 : *pf5++ = 0.8f * pf4[3] + 0.2f * pf4[4];
862 649116 : pf4 += 4;
863 : }
864 11388 : *pf5++ = pf4[0];
865 11388 : *pf5++ = 0.2f * pf4[0] + 0.8f * pf4[1];
866 11388 : *pf5++ = 0.4f * pf4[1] + 0.6f * pf4[2];
867 11388 : *pf5++ = 0.6f * pf4[2] + 0.4f * pf4[3];
868 11388 : *pf5++ = 0.8f * pf4[3] + 0.2f * syn[0];
869 : /* all samples processed: NBPSF_PIT_MAX = 290 = (58*5) */
870 :
871 11388 : return;
872 : }
873 :
874 :
875 : /*-------------------------------------------------------------------*
876 : * retro_interp5_4()
877 : *
878 : *
879 : *-------------------------------------------------------------------*/
880 :
881 24446 : void retro_interp5_4(
882 : float *pst_old_syn )
883 : {
884 : float *pf5, *pf4;
885 : int16_t c;
886 :
887 : /* resample st->pst_old_syn in a reverse way to preserve time-alignment */
888 24446 : pf4 = (float *) &pst_old_syn[NBPSF_PIT_MAX - 1];
889 24446 : pf5 = pf4;
890 1442314 : for ( c = 0; c < 58; c++ )
891 : {
892 1417868 : *pf4-- = 0.75f * pf5[0] + 0.25f * pf5[-1];
893 1417868 : *pf4-- = 0.50f * pf5[-1] + 0.50f * pf5[-2];
894 1417868 : *pf4-- = 0.25f * pf5[-2] + 0.75f * pf5[-3];
895 1417868 : *pf4-- = pf5[-4];
896 1417868 : pf5 -= 5;
897 : }
898 : /* all samples processed: NBPSF_PIT_MAX = 290 = (58*5) */
899 :
900 24446 : return;
901 : }
|
