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 3667821 : 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 3667821 : 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 3667821 : 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 48432750 : 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 3667821 : 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 3667821 : fac_num = cfg_ptr->fac_num;
119 3667821 : fac_den = (int16_t) ( ( cfg_ptr->fin * fac_num ) / cfg_ptr->fout );
120 3667821 : lg_out = ( lg * fac_num ) / fac_den;
121 3667821 : filt_len = cfg_ptr->filt_len;
122 :
123 3667821 : mem_len = 2 * filt_len;
124 3667821 : plus_sample_in = 0; /* default, regular delay */
125 3667821 : frac = 0;
126 :
127 3667821 : if ( fin == 8000 && fout == 12800 )
128 : {
129 0 : plus_sample_in = 7;
130 0 : frac = 4;
131 : }
132 :
133 3667821 : signal = signal_tab + 2 * L_FILT_MAX + L_FRAME48k - mem_len - lg;
134 3667821 : signal_ana = signal;
135 3667821 : 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 302126709 : for ( i = 0; i < 2 * filt_len; i++ )
145 : {
146 298458888 : signal[i] = mem[i];
147 : }
148 :
149 1608879073 : for ( i = 0; i < lg; i++ )
150 : {
151 1605211252 : signal[i + ( 2 * filt_len )] = sigIn[i];
152 : }
153 :
154 3667821 : if ( plus_sample_in > 0 )
155 : {
156 0 : autocorr( signal_ana + mem_len_ana + lg - LEN_WIN_SSS, r, 1, LEN_WIN_SSS, wind_sss, 0, 0, 0 );
157 :
158 0 : mu_preemph = r[1] / r[0];
159 0 : mem_preemph = signal_ana[mem_len_ana + lg - LEN_WIN_SSS - 1];
160 0 : preemph( signal_ana + mem_len_ana + lg - LEN_WIN_SSS, mu_preemph, LEN_WIN_SSS, &mem_preemph );
161 :
162 : /* Autocorrelations */
163 0 : autocorr( signal_ana + mem_len_ana + lg - LEN_WIN_SSS, r, M, LEN_WIN_SSS, wind_sss, 0, 0, 0 );
164 :
165 0 : lag_wind( r, M, fin, LAGW_STRONG );
166 :
167 : /* Levinson-Durbin */
168 0 : lev_dur( A, r, M, epsP );
169 :
170 0 : for ( i = 0; i < plus_sample_in; i++ )
171 : {
172 0 : val = 0;
173 0 : for ( j = 1; j <= M; j++ )
174 : {
175 0 : val -= signal[i + lg + mem_len - j] * A[j];
176 : }
177 0 : signal[i + lg + mem_len] = val; /* AZ ringing padding */
178 : }
179 :
180 0 : mem_preemph = signal[mem_len + lg - LEN_WIN_SSS - 1];
181 0 : deemph( signal + mem_len + lg - LEN_WIN_SSS, mu_preemph, LEN_WIN_SSS + plus_sample_in, &mem_preemph );
182 : }
183 :
184 : /* interpolation */
185 3667821 : datastep = fac_den / fac_num;
186 3667821 : fracstep = fac_den - datastep * fac_num; /* equivalent to datastep = fac_den % fac_num */
187 :
188 3667821 : sigIn_ptr = signal + filt_len + plus_sample_in;
189 :
190 3667821 : filt_len_tmp = filt_len;
191 3667821 : if ( flag_low_order )
192 : {
193 0 : filt_len_tmp = ( filt_len + 1 ) >> 1;
194 : }
195 :
196 538532309 : for ( i = 0; i < lg_out; i++ )
197 : {
198 534864488 : sigOut[i] = interpolation( sigIn_ptr, cfg_ptr->filter, frac, fac_num, filt_len_tmp );
199 :
200 534864488 : frac = frac + fracstep;
201 534864488 : if ( frac >= fac_num )
202 : {
203 289408259 : frac = frac - fac_num;
204 289408259 : sigIn_ptr++;
205 : }
206 :
207 534864488 : sigIn_ptr += datastep;
208 : }
209 :
210 : /* rescaling */
211 3667821 : if ( ( fac_num > fac_den ) == ( ( cfg_ptr->flags & RS_INV_FAC ) != 0 ) )
212 : {
213 3662989 : num_den = (float) fac_num / fac_den;
214 :
215 536057237 : for ( i = 0; i < lg_out; i++ )
216 : {
217 532394248 : sigOut[i] *= num_den;
218 : }
219 : }
220 :
221 : /* update the filter memory */
222 302126709 : for ( i = 0; i < 2 * filt_len; i++ )
223 : {
224 298458888 : mem[i] = signal[i + lg];
225 : }
226 :
227 3667821 : return lg_out;
228 : }
229 :
230 : /*-------------------------------------------------------------------*
231 : * modify_Fs_intcub3m_sup()
232 : *
233 : *
234 : *-------------------------------------------------------------------*/
235 :
236 : /*! r: length of output */
237 75 : 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 75 : const double( *cu )[3] = 0;
251 :
252 : /*-------------------------------------------------------------------*
253 : * Find the resampling configuration
254 : *-------------------------------------------------------------------*/
255 :
256 : /* check if fin and fout are the same */
257 75 : 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 75 : 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 75 : cc[2][0] = sigIn[-1] / 3;
282 75 : cc[2][2] = ( sigIn[-2] + sigIn[0] ) / 2 - sigIn[-1];
283 75 : cc[2][3] = ( sigIn[-2] + sigIn[1] - sigIn[-1] - sigIn[0] - 4 * cc[2][2] ) / 6;
284 75 : cc[2][1] = sigIn[0] - sigIn[-1] - cc[2][3] - cc[2][2];
285 :
286 : /* cc[3][] : indexes -1 0 1 2 */
287 75 : cc[3][0] = sigIn[0] / 3;
288 75 : cc[3][2] = ( sigIn[-1] + sigIn[1] ) / 2 - sigIn[0];
289 75 : cc[3][3] = ( sigIn[-1] + sigIn[2] - sigIn[0] - sigIn[1] - 4 * cc[3][2] ) / 6;
290 75 : cc[3][1] = sigIn[1] - sigIn[0] - cc[3][3] - cc[3][2];
291 75 : j = 0;
292 :
293 75 : if ( fin == 12800 )
294 : {
295 15 : if ( fout == 8000 )
296 : {
297 0 : cind = 0;
298 : }
299 15 : else if ( fout == 16000 )
300 : {
301 0 : cind = 1;
302 : }
303 15 : else if ( fout == 32000 )
304 : {
305 15 : cind = 2;
306 : }
307 0 : else if ( fout == 48000 )
308 : {
309 0 : 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 60 : else if ( fin == 16000 )
318 : {
319 60 : if ( fout == 12800 )
320 : {
321 0 : cind = 4;
322 : }
323 60 : else if ( fout == 32000 )
324 : {
325 0 : cind = 5;
326 : }
327 60 : else if ( fout == 48000 )
328 : {
329 60 : 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 75 : *delayout = ct2[cind][9];
344 :
345 75 : if ( ct2[cind][12] == 15 )
346 : {
347 75 : cu = cu15;
348 : }
349 :
350 75 : if ( ct2[cind][12] == 4 )
351 : {
352 0 : cu = cu4;
353 : }
354 :
355 75 : fk1 = 2 * ct2[cind][12] - 2;
356 75 : k2d = fk1 / 2; /* shift of index in cu with respect to the next sample (ex 1.25 -> 0.25 ) */
357 75 : k3d = fk1 - 1; /* to compurte index in cu with respect to the last sample with - sign (ex 1.25 -> -0.75 ) */
358 :
359 75 : kk = 0;
360 1125 : for ( i = 0; i < lg - ct2[cind][11]; )
361 : {
362 1050 : sigOut[j++] = sigIn[i];
363 2190 : for ( k = 0; k < ct2[cind][10]; k++ )
364 : {
365 1140 : cc[kk][0] = sigIn[i + 1] / 3;
366 1140 : cc[kk][2] = ( sigIn[i] + sigIn[i + 2] ) / 2 - sigIn[i + 1];
367 1140 : cc[kk][3] = ( sigIn[i] + sigIn[i + 3] - sigIn[i + 1] - sigIn[i + 2] - 4 * cc[kk][2] ) / 6;
368 1140 : cc[kk][1] = sigIn[i + 2] - sigIn[i + 1] - cc[kk][3] - cc[kk][2];
369 1140 : i++;
370 :
371 1140 : i2 = kk - 2;
372 1140 : i1 = kk - 1;
373 1140 : if ( i1 < 0 )
374 : {
375 285 : i1 += 4;
376 : }
377 :
378 1140 : if ( i2 < 0 )
379 : {
380 570 : i2 += 4;
381 : }
382 :
383 3420 : for ( k1 = ct2[cind][k]; k1 < fk1; k1 += ct2[cind][8] )
384 : {
385 2280 : k2 = k1 - k2d;
386 2280 : k3 = k3d - k1;
387 2280 : vv = (float) ( cu[k1][2] * cc[i2][3] + cu[k1][1] * cc[i2][2] + cu[k1][0] * cc[i2][1] + cc[i2][0] );
388 2280 : vv += (float) ( cu[k2][2] * cc[i1][3] + cu[k2][1] * cc[i1][2] + cu[k2][0] * cc[i1][1] + cc[i1][0] );
389 2280 : vv += (float) ( -cu[k3][2] * cc[kk][3] + cu[k3][1] * cc[kk][2] - cu[k3][0] * cc[kk][1] + cc[kk][0] );
390 2280 : sigOut[j++] = vv;
391 : }
392 :
393 1140 : kk++;
394 1140 : if ( kk == 4 )
395 : {
396 285 : kk = 0;
397 : }
398 : }
399 : }
400 :
401 75 : sigOut[j++] = sigIn[i];
402 :
403 150 : for ( k = 0; k < ct2[cind][11] - 3; k++ )
404 : {
405 75 : cc[kk][0] = sigIn[i + 1] / 3;
406 75 : cc[kk][2] = ( sigIn[i] + sigIn[i + 2] ) / 2 - sigIn[i + 1];
407 75 : cc[kk][3] = ( sigIn[i] + sigIn[i + 3] - sigIn[i + 1] - sigIn[i + 2] - 4 * cc[kk][2] ) / 6;
408 75 : cc[kk][1] = sigIn[i + 2] - sigIn[i + 1] - cc[kk][3] - cc[kk][2];
409 75 : i++;
410 :
411 75 : i2 = kk - 2;
412 75 : i1 = kk - 1;
413 75 : if ( i1 < 0 )
414 : {
415 75 : i1 += 4;
416 : }
417 :
418 75 : if ( i2 < 0 )
419 : {
420 75 : i2 += 4;
421 : }
422 :
423 225 : for ( k1 = ct2[cind][k]; k1 < fk1; k1 += ct2[cind][8] )
424 : {
425 150 : k2 = k1 - k2d;
426 150 : k3 = k3d - k1;
427 150 : vv = (float) ( cu[k1][2] * cc[i2][3] + cu[k1][1] * cc[i2][2] + cu[k1][0] * cc[i2][1] + cc[i2][0] );
428 150 : vv += (float) ( cu[k2][2] * cc[i1][3] + cu[k2][1] * cc[i1][2] + cu[k2][0] * cc[i1][1] + cc[i1][0] );
429 150 : vv += (float) ( -cu[k3][2] * cc[kk][3] + cu[k3][1] * cc[kk][2] - cu[k3][0] * cc[kk][1] + cc[kk][0] );
430 150 : sigOut[j++] = vv;
431 : }
432 :
433 75 : kk++;
434 :
435 75 : if ( kk == 4 )
436 : {
437 0 : kk = 0;
438 : }
439 : }
440 :
441 75 : kk--;
442 75 : if ( kk == -1 )
443 : {
444 0 : kk = 3;
445 : }
446 :
447 75 : if ( ct2[cind][10] == 1 )
448 : {
449 60 : sigOut[j++] = sigIn[i];
450 : }
451 :
452 225 : for ( k1 = ct2[cind][k]; k1 < fk1; k1 += ct2[cind][8] )
453 : {
454 150 : k2 = k1 - k2d;
455 150 : vv = (float) ( cu[k2][2] * cc[kk][3] + cu[k2][1] * cc[kk][2] + cu[k2][0] * cc[kk][1] + cc[kk][0] );
456 150 : sigOut[j++] = vv * 3;
457 : }
458 :
459 75 : if ( ct2[cind][10] < 3 )
460 : {
461 75 : sigOut[j++] = sigIn[i + 1];
462 : }
463 :
464 225 : for ( k1 = ct2[cind][k + 1]; k1 < fk1; k1 += ct2[cind][8] )
465 : {
466 150 : vv = (float) ( cu[k1][2] * cc[kk][3] + cu[k1][1] * cc[kk][2] + cu[k1][0] * cc[kk][1] + cc[kk][0] );
467 150 : sigOut[j++] = vv * 3;
468 : }
469 :
470 75 : if ( ct2[cind][10] == 1 )
471 : {
472 60 : sigOut[j++] = sigIn[i + 2];
473 : }
474 : }
475 :
476 75 : return lg_out;
477 : }
478 :
479 : /*-------------------------------------------------------------------*
480 : * Interpolate_allpass_steep()
481 : *
482 : * Interpolation by a factor 2
483 : *-------------------------------------------------------------------*/
484 :
485 458871 : 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 128449431 : for ( k = 0; k < N; k++ )
497 : {
498 127990560 : temp[0] = mem[0] + AP2_STEEP[0] * in[k];
499 127990560 : mem[0] = in[k] - AP2_STEEP[0] * temp[0];
500 :
501 : /* for better performance, unroll this loop */
502 255981120 : for ( n = 1; n < ALLPASSSECTIONS_STEEP - 1; n++ )
503 : {
504 127990560 : temp[n] = mem[n] + AP2_STEEP[n] * temp[n - 1];
505 127990560 : mem[n] = temp[n - 1] - AP2_STEEP[n] * temp[n];
506 : }
507 :
508 127990560 : out[2 * k + 1] = mem[ALLPASSSECTIONS_STEEP - 1] + AP2_STEEP[ALLPASSSECTIONS_STEEP - 1] * temp[ALLPASSSECTIONS_STEEP - 2];
509 127990560 : 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 128449431 : for ( k = 0; k < N; k++ )
514 : {
515 127990560 : temp[0] = mem[ALLPASSSECTIONS_STEEP] + AP1_STEEP[0] * in[k];
516 127990560 : mem[ALLPASSSECTIONS_STEEP] = in[k] - AP1_STEEP[0] * temp[0];
517 :
518 : /* for better performance, unroll this loop */
519 255981120 : for ( n = 1; n < ALLPASSSECTIONS_STEEP - 1; n++ )
520 : {
521 127990560 : temp[n] = mem[ALLPASSSECTIONS_STEEP + n] + AP1_STEEP[n] * temp[n - 1];
522 127990560 : mem[ALLPASSSECTIONS_STEEP + n] = temp[n - 1] - AP1_STEEP[n] * temp[n];
523 : }
524 :
525 127990560 : out[2 * k] = mem[2 * ALLPASSSECTIONS_STEEP - 1] + AP1_STEEP[ALLPASSSECTIONS_STEEP - 1] * temp[ALLPASSSECTIONS_STEEP - 2];
526 127990560 : mem[2 * ALLPASSSECTIONS_STEEP - 1] = temp[ALLPASSSECTIONS_STEEP - 2] - AP1_STEEP[ALLPASSSECTIONS_STEEP - 1] * out[2 * k];
527 : }
528 :
529 458871 : return;
530 : }
531 :
532 : /*-------------------------------------------------------------------*
533 : * Decimate_allpass_steep()
534 : *
535 : * Decimation by a factor 2
536 : *-------------------------------------------------------------------*/
537 :
538 881902 : 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 250911881 : for ( k = 0; k < N / 2; k++ )
550 : {
551 250029979 : temp[0] = mem[0] + AP1_STEEP[0] * in[2 * k];
552 250029979 : mem[0] = in[2 * k] - AP1_STEEP[0] * temp[0];
553 :
554 : /* for better performance, unroll this loop */
555 500059958 : for ( n = 1; n < ALLPASSSECTIONS_STEEP - 1; n++ )
556 : {
557 250029979 : temp[n] = mem[n] + AP1_STEEP[n] * temp[n - 1];
558 250029979 : if ( fabs( temp[n] ) < 1e-12 )
559 : {
560 978835 : temp[n] = sign( temp[n] ) * 1e-12f;
561 : }
562 250029979 : mem[n] = temp[n - 1] - AP1_STEEP[n] * temp[n];
563 : }
564 :
565 250029979 : out[k] = mem[ALLPASSSECTIONS_STEEP - 1] + AP1_STEEP[ALLPASSSECTIONS_STEEP - 1] * temp[ALLPASSSECTIONS_STEEP - 2];
566 250029979 : mem[ALLPASSSECTIONS_STEEP - 1] = temp[ALLPASSSECTIONS_STEEP - 2] - AP1_STEEP[ALLPASSSECTIONS_STEEP - 1] * out[k];
567 : }
568 :
569 : /* lower allpass filter chain */
570 881902 : temp[0] = mem[ALLPASSSECTIONS_STEEP] + AP2_STEEP[0] * mem[2 * ALLPASSSECTIONS_STEEP];
571 881902 : mem[ALLPASSSECTIONS_STEEP] = mem[2 * ALLPASSSECTIONS_STEEP] - AP2_STEEP[0] * temp[0];
572 :
573 : /* for better performance, unroll this loop */
574 1763804 : for ( n = 1; n < ALLPASSSECTIONS_STEEP - 1; n++ )
575 : {
576 881902 : temp[n] = mem[ALLPASSSECTIONS_STEEP + n] + AP2_STEEP[n] * temp[n - 1];
577 881902 : if ( fabs( temp[n] ) < 1e-12 )
578 : {
579 273723 : temp[n] = sign( temp[n] ) * 1e-12f;
580 : }
581 881902 : mem[ALLPASSSECTIONS_STEEP + n] = temp[n - 1] - AP2_STEEP[n] * temp[n];
582 : }
583 :
584 881902 : temp[ALLPASSSECTIONS_STEEP - 1] = mem[2 * ALLPASSSECTIONS_STEEP - 1] + AP2_STEEP[ALLPASSSECTIONS_STEEP - 1] *
585 881902 : temp[ALLPASSSECTIONS_STEEP - 2];
586 :
587 881902 : mem[2 * ALLPASSSECTIONS_STEEP - 1] = temp[ALLPASSSECTIONS_STEEP - 2] - AP2_STEEP[ALLPASSSECTIONS_STEEP - 1] *
588 881902 : temp[ALLPASSSECTIONS_STEEP - 1];
589 881902 : out[0] = (float) ( ( out[0] + temp[ALLPASSSECTIONS_STEEP - 1] ) * 0.5 );
590 :
591 250029979 : for ( k = 1; k < N / 2; k++ )
592 : {
593 249148077 : temp[0] = mem[ALLPASSSECTIONS_STEEP] + AP2_STEEP[0] * in[2 * k - 1];
594 249148077 : mem[ALLPASSSECTIONS_STEEP] = in[2 * k - 1] - AP2_STEEP[0] * temp[0];
595 :
596 : /* for better performance, unroll this loop */
597 498296154 : for ( n = 1; n < ALLPASSSECTIONS_STEEP - 1; n++ )
598 : {
599 249148077 : temp[n] = mem[ALLPASSSECTIONS_STEEP + n] + AP2_STEEP[n] * temp[n - 1];
600 249148077 : if ( fabs( temp[n] ) < 1e-12 )
601 : {
602 873627 : temp[n] = sign( temp[n] ) * 1e-12f;
603 : }
604 249148077 : mem[ALLPASSSECTIONS_STEEP + n] = temp[n - 1] - AP2_STEEP[n] * temp[n];
605 : }
606 :
607 249148077 : temp[ALLPASSSECTIONS_STEEP - 1] = mem[2 * ALLPASSSECTIONS_STEEP - 1] + AP2_STEEP[ALLPASSSECTIONS_STEEP - 1] *
608 249148077 : temp[ALLPASSSECTIONS_STEEP - 2];
609 249148077 : mem[2 * ALLPASSSECTIONS_STEEP - 1] = temp[ALLPASSSECTIONS_STEEP - 2] - AP2_STEEP[ALLPASSSECTIONS_STEEP - 1] *
610 249148077 : temp[ALLPASSSECTIONS_STEEP - 1];
611 249148077 : out[k] = (float) ( ( out[k] + temp[ALLPASSSECTIONS_STEEP - 1] ) * 0.5 );
612 : }
613 :
614 : /* z^(-1) */
615 881902 : mem[2 * ALLPASSSECTIONS_STEEP] = in[N - 1];
616 :
617 881902 : 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 280968 : 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 280968 : const float *filt_coeff = allpass_poles_3_ov_2;
639 :
640 280968 : out1 = out1_buff;
641 :
642 173979714 : for ( i = 0; i < len; i++ )
643 : {
644 : /* Upper branch */
645 173698746 : Vu[0] = mem[0] + filt_coeff[0] * ( input[i] - mem[1] );
646 173698746 : Vu[1] = mem[1] + filt_coeff[1] * ( Vu[0] - mem[2] );
647 173698746 : mem[3] = mem[2] + filt_coeff[2] * ( Vu[1] - mem[3] );
648 :
649 173698746 : mem[1] = Vu[0];
650 173698746 : mem[2] = Vu[1];
651 173698746 : *out1++ = mem[3];
652 :
653 : /* Middle branch */
654 173698746 : Vm[0] = mem[0] + filt_coeff[3] * ( input[i] - mem[4] );
655 173698746 : Vm[1] = mem[4] + filt_coeff[4] * ( Vm[0] - mem[5] );
656 173698746 : mem[6] = mem[5] + filt_coeff[5] * ( Vm[1] - mem[6] );
657 :
658 173698746 : mem[4] = Vm[0];
659 173698746 : mem[5] = Vm[1];
660 173698746 : *out1++ = mem[6];
661 :
662 : /* Lower branch */
663 173698746 : Vl[0] = mem[0] + filt_coeff[6] * ( input[i] - mem[7] );
664 173698746 : Vl[1] = mem[7] + filt_coeff[7] * ( Vl[0] - mem[8] );
665 173698746 : mem[9] = mem[8] + filt_coeff[8] * ( Vl[1] - mem[9] );
666 :
667 173698746 : mem[0] = input[i];
668 173698746 : mem[7] = Vl[0];
669 173698746 : mem[8] = Vl[1];
670 173698746 : *out1++ = mem[9];
671 : }
672 :
673 280968 : loop_len = len * 3 / 2;
674 :
675 : /*decimate by 2 and LPF*/
676 260829087 : for ( i = 0; i < loop_len; i++ )
677 : {
678 260548119 : mem_temp = out1_buff[2 * i];
679 260548119 : out[i] = ( ( ( 0.0473147f ) * ( mem_temp + mem[10] ) ) + ( ( -0.151521f ) * ( mem[11] + mem[14] ) ) );
680 260548119 : out[i] = ( out[i] + ( ( 0.614152f ) * ( mem[12] + mem[13] ) ) );
681 260548119 : mem[10] = mem[11];
682 260548119 : mem[11] = mem[12];
683 260548119 : mem[12] = mem[13];
684 260548119 : mem[13] = mem[14];
685 260548119 : mem[14] = mem_temp;
686 : }
687 :
688 280968 : return;
689 : }
690 :
691 : /*-------------------------------------------------------------------*
692 : * decimate_2_over_3_allpass()
693 : *
694 : * Decimate 2/3 using allpass iir polyphase filter.
695 : *-------------------------------------------------------------------*/
696 :
697 128006 : 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 128006 : const float *filt_coeff = allpass_poles_3_ov_2;
712 128006 : const float *lp_num = decimate_3_ov_2_lowpass_num;
713 128006 : 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 128006 : out1 = out1_buff;
717 :
718 128006 : *out1++ = lp_num[0] * ( input[0] + lp_mem[0] ) - lp_den[2] * lp_mem[2];
719 128006 : *out1++ = lp_num[1] * input[0] - lp_den[2] * lp_mem[1];
720 :
721 112423344 : for ( i = 1; i < len; i++ )
722 : {
723 112295338 : tmp = lp_num[0] * ( input[i] + input[i - 1] ) - lp_den[2] * out1[-2];
724 112295338 : *out1++ = tmp;
725 112295338 : tmp = lp_num[1] * input[i] - lp_den[2] * out1[-2];
726 112295338 : *out1++ = tmp;
727 : }
728 128006 : lp_mem[0] = input[len - 1];
729 128006 : lp_mem[1] = out1[-1];
730 128006 : lp_mem[2] = out1[-2];
731 :
732 : /* do the all pass polyphase filter with pi/3 cutoff */
733 128006 : out1 = out;
734 128006 : in = out1_buff;
735 128006 : loop_len = (int16_t) len * 2 / 3;
736 :
737 75076902 : for ( i = 0; i < loop_len; i++ )
738 : {
739 : /* Lower branch */
740 74948896 : Vl[0] = mem[8] + filt_coeff[6] * ( *in - mem[9] );
741 74948896 : Vl[1] = mem[9] + filt_coeff[7] * ( Vl[0] - mem[10] );
742 74948896 : mem[11] = mem[10] + filt_coeff[8] * ( Vl[1] - mem[11] );
743 :
744 74948896 : mem[8] = *in++;
745 74948896 : mem[9] = Vl[0];
746 74948896 : mem[10] = Vl[1];
747 74948896 : *out1 = mem[11];
748 :
749 : /* Middle branch */
750 74948896 : Vm[0] = mem[4] + filt_coeff[3] * ( *in - mem[5] );
751 74948896 : Vm[1] = mem[5] + filt_coeff[4] * ( Vm[0] - mem[6] );
752 74948896 : mem[7] = mem[6] + filt_coeff[5] * ( Vm[1] - mem[7] );
753 :
754 74948896 : mem[4] = *in++;
755 74948896 : mem[5] = Vm[0];
756 74948896 : mem[6] = Vm[1];
757 74948896 : *out1 += mem[7];
758 :
759 : /* Upper branch */
760 74948896 : Vu[0] = mem[0] + filt_coeff[0] * ( *in - mem[1] );
761 74948896 : Vu[1] = mem[1] + filt_coeff[1] * ( Vu[0] - mem[2] );
762 74948896 : mem[3] = mem[2] + filt_coeff[2] * ( Vu[1] - mem[3] );
763 :
764 74948896 : mem[0] = *in++;
765 74948896 : mem[1] = Vu[0];
766 74948896 : mem[2] = Vu[1];
767 74948896 : *out1++ += mem[3];
768 : }
769 :
770 128006 : 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 168624 : 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 168624 : const float *filt_coeff = allpass_poles_3_ov_2;
791 :
792 168624 : out1 = &out[0];
793 :
794 16060464 : for ( i = 0; i < len; i++ )
795 : {
796 : /* Upper branch */
797 15891840 : Vu[0] = mem[0] + filt_coeff[0] * ( input[i] - mem[1] );
798 15891840 : Vu[1] = mem[1] + filt_coeff[1] * ( Vu[0] - mem[2] );
799 15891840 : mem[3] = mem[2] + filt_coeff[2] * ( Vu[1] - mem[3] );
800 :
801 15891840 : mem[1] = Vu[0];
802 15891840 : mem[2] = Vu[1];
803 15891840 : *out1++ = mem[3];
804 :
805 : /* Middle branch */
806 15891840 : Vm[0] = mem[0] + filt_coeff[3] * ( input[i] - mem[4] );
807 15891840 : Vm[1] = mem[4] + filt_coeff[4] * ( Vm[0] - mem[5] );
808 15891840 : mem[6] = mem[5] + filt_coeff[5] * ( Vm[1] - mem[6] );
809 :
810 15891840 : mem[4] = Vm[0];
811 15891840 : mem[5] = Vm[1];
812 15891840 : *out1++ = mem[6];
813 :
814 : /* Lower branch */
815 15891840 : Vl[0] = mem[0] + filt_coeff[6] * ( input[i] - mem[7] );
816 15891840 : Vl[1] = mem[7] + filt_coeff[7] * ( Vl[0] - mem[8] );
817 15891840 : mem[9] = mem[8] + filt_coeff[8] * ( Vl[1] - mem[9] );
818 :
819 15891840 : mem[0] = input[i];
820 15891840 : mem[7] = Vl[0];
821 15891840 : mem[8] = Vl[1];
822 15891840 : *out1++ = mem[9];
823 : }
824 :
825 : /*LPF*/
826 47844144 : for ( i = 0; i < len * 3; i++ )
827 : {
828 47675520 : mem_temp = out[i];
829 47675520 : out[i] = ( ( ( 0.572769f ) * ( mem[12] + mem[11] ) ) - ( ( 0.074005f ) * ( mem_temp + mem[10] ) ) );
830 47675520 : mem[10] = mem[11];
831 47675520 : mem[11] = mem[12];
832 47675520 : mem[12] = mem_temp;
833 : }
834 :
835 168624 : return;
836 : }
837 :
838 :
839 : /*-------------------------------------------------------------------*
840 : * retro_interp4_5()
841 : *
842 : *
843 : *-------------------------------------------------------------------*/
844 :
845 1674 : 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 1674 : pf4 = (float *) &pst_old_syn[58];
854 1674 : pf5 = (float *) pst_old_syn;
855 97092 : for ( c = 0; c < 57; c++ )
856 : {
857 95418 : *pf5++ = pf4[0];
858 95418 : *pf5++ = 0.2f * pf4[0] + 0.8f * pf4[1];
859 95418 : *pf5++ = 0.4f * pf4[1] + 0.6f * pf4[2];
860 95418 : *pf5++ = 0.6f * pf4[2] + 0.4f * pf4[3];
861 95418 : *pf5++ = 0.8f * pf4[3] + 0.2f * pf4[4];
862 95418 : pf4 += 4;
863 : }
864 1674 : *pf5++ = pf4[0];
865 1674 : *pf5++ = 0.2f * pf4[0] + 0.8f * pf4[1];
866 1674 : *pf5++ = 0.4f * pf4[1] + 0.6f * pf4[2];
867 1674 : *pf5++ = 0.6f * pf4[2] + 0.4f * pf4[3];
868 1674 : *pf5++ = 0.8f * pf4[3] + 0.2f * syn[0];
869 : /* all samples processed: NBPSF_PIT_MAX = 290 = (58*5) */
870 :
871 1674 : return;
872 : }
873 :
874 :
875 : /*-------------------------------------------------------------------*
876 : * retro_interp5_4()
877 : *
878 : *
879 : *-------------------------------------------------------------------*/
880 :
881 3714 : 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 3714 : pf4 = (float *) &pst_old_syn[NBPSF_PIT_MAX - 1];
889 3714 : pf5 = pf4;
890 219126 : for ( c = 0; c < 58; c++ )
891 : {
892 215412 : *pf4-- = 0.75f * pf5[0] + 0.25f * pf5[-1];
893 215412 : *pf4-- = 0.50f * pf5[-1] + 0.50f * pf5[-2];
894 215412 : *pf4-- = 0.25f * pf5[-2] + 0.75f * pf5[-3];
895 215412 : *pf4-- = pf5[-4];
896 215412 : pf5 -= 5;
897 : }
898 : /* all samples processed: NBPSF_PIT_MAX = 290 = (58*5) */
899 :
900 3714 : return;
901 : }
|
