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 : #define _USE_MATH_DEFINES
38 :
39 : #include <assert.h>
40 : #include <stdint.h>
41 : #include "options.h"
42 : #ifdef DEBUGGING
43 : #include "debug.h"
44 : #endif
45 : #include <math.h>
46 : #include "prot.h"
47 : #include "cnst.h"
48 : #include "stat_com.h"
49 : #include "wmc_auto.h"
50 : #include "ivas_prot.h"
51 :
52 :
53 : /*---------------------------------------------------------------------*
54 : * Local function prototypes
55 : *---------------------------------------------------------------------*/
56 :
57 : static void findStrongestHarmonics( const uint16_t nSamples, const float *powerSpectrum, const float F0, const uint16_t nTotalHarmonics, uint16_t *pHarmonicIndexes, int16_t *pnHarmonics );
58 : static void CorrectF0( const uint16_t *pHarmonicIndexes, const uint16_t nHarmonics, float *pF0 );
59 : static void findCandidates( const uint16_t nSamples, const float *MDCTSpectrum, float *thresholdModificationNew, float floorPowerSpectrum );
60 : static void RefineThresholdsUsingPitch( const uint16_t nSamples, const uint16_t nSamplesCore, const float *powerSpectrum, const float lastPitchLag, const float currentPitchLag, float *pF0, float *thresholdModification );
61 : static void findTonalComponents( uint16_t *indexOfTonalPeak, uint16_t *lowerIndex, uint16_t *upperIndex, uint16_t *numIndexes, const uint16_t nSamples, const float *powerSpectrum, const float F0, const float *thresholdModification );
62 :
63 :
64 : /*-------------------------------------------------------------------*
65 : * DetectTonalComponents()
66 : *
67 : * Detect tonal components in the lastMDCTSpectrum, use
68 : * secondLastPowerSpectrum for the precise location of the peaks and
69 : * store them in indexOfTonalPeak. Updates lowerIndex, upperIndex,
70 : * pNumIndexes accordingly.
71 : *-------------------------------------------------------------------*/
72 :
73 6066 : void DetectTonalComponents(
74 : uint16_t indexOfTonalPeak[],
75 : uint16_t lowerIndex[],
76 : uint16_t upperIndex[],
77 : uint16_t *pNumIndexes,
78 : const float lastPitchLag,
79 : const float currentPitchLag,
80 : const float lastMDCTSpectrum[],
81 : const float scaleFactors[],
82 : const float secondLastPowerSpectrum[],
83 : const uint16_t nSamples,
84 : const uint16_t nSamplesCore,
85 : float floorPowerSpectrum, /* i : lower limit for power spectrum bins */
86 : const PsychoacousticParameters *psychParamsCurrent )
87 : {
88 : float F0;
89 : float thresholdModification[L_FRAME_MAX];
90 6066 : float *pScaledMdctSpectrum = thresholdModification;
91 : int16_t nBands;
92 :
93 : /* Convert from 16 bit to 32 bit */
94 6066 : mvr2r( lastMDCTSpectrum, pScaledMdctSpectrum, nSamples );
95 6066 : if ( psychParamsCurrent == NULL )
96 : {
97 3135 : nBands = FDNS_NPTS;
98 3135 : mdct_noiseShaping( pScaledMdctSpectrum, nSamplesCore, scaleFactors, nBands );
99 : }
100 : else
101 : {
102 2931 : sns_shape_spectrum( pScaledMdctSpectrum, psychParamsCurrent, scaleFactors, nSamplesCore );
103 2931 : nBands = psychParamsCurrent->nBands;
104 : }
105 :
106 6066 : v_multc( pScaledMdctSpectrum + nSamplesCore, scaleFactors[nBands - 1], pScaledMdctSpectrum + nSamplesCore, nSamples - nSamplesCore );
107 :
108 : /* Find peak candidates in the last frame. */
109 6066 : findCandidates( nSamples, pScaledMdctSpectrum, thresholdModification, floorPowerSpectrum );
110 :
111 : /* Refine peak candidates using the pitch information */
112 6066 : RefineThresholdsUsingPitch( nSamples, nSamplesCore, secondLastPowerSpectrum, lastPitchLag, currentPitchLag, &F0, thresholdModification );
113 :
114 : /* Find peaks in the second last frame */
115 6066 : findTonalComponents( indexOfTonalPeak, lowerIndex, upperIndex, pNumIndexes, nSamples, secondLastPowerSpectrum, F0, thresholdModification );
116 :
117 6066 : return;
118 : }
119 :
120 :
121 : /*-------------------------------------------------------------------*
122 : * RefineTonalComponents()
123 : *
124 : *-------------------------------------------------------------------*/
125 :
126 : /* When called, the tonal components are already stored in
127 : * indexOfTonalPeak. Detect tonal components in the lastMDCTSpectrum,
128 : * use secondLastPowerSpectrum for the precise location of the peaks and
129 : * then keep in indexOfTonalPeak only the tonal components that are
130 : * again detected Updates indexOfTonalPeak, lowerIndex, upperIndex,
131 : * phaseDiff, phases, pNumIndexes accordingly. */
132 135 : void RefineTonalComponents(
133 : uint16_t indexOfTonalPeak[],
134 : uint16_t lowerIndex[],
135 : uint16_t upperIndex[],
136 : float phaseDiff[],
137 : float phases[],
138 : uint16_t *pNumIndexes,
139 : const float lastPitchLag,
140 : const float currentPitchLag,
141 : const float lastMDCTSpectrum[],
142 : const float scaleFactors[],
143 : const float secondLastPowerSpectrum[],
144 : const uint16_t nSamples,
145 : const uint16_t nSamplesCore,
146 : float floorPowerSpectrum, /* i : lower limit for power spectrum bins */
147 : const PsychoacousticParameters *psychParamsCurrent )
148 : {
149 : uint16_t newIndexOfTonalPeak[MAX_NUMBER_OF_IDX];
150 : uint16_t newLowerIndex[MAX_NUMBER_OF_IDX];
151 : uint16_t newUpperIndex[MAX_NUMBER_OF_IDX];
152 : uint16_t newNumIndexes, nPreservedPeaks;
153 : uint16_t iNew, iOld, j;
154 : float *pOldPhase, *pNewPhase;
155 :
156 135 : DetectTonalComponents( newIndexOfTonalPeak, newLowerIndex, newUpperIndex, &newNumIndexes, lastPitchLag, currentPitchLag, lastMDCTSpectrum, scaleFactors, secondLastPowerSpectrum,
157 : nSamples, nSamplesCore, floorPowerSpectrum, psychParamsCurrent );
158 :
159 135 : nPreservedPeaks = 0;
160 135 : iNew = 0;
161 135 : pOldPhase = phases;
162 135 : pNewPhase = phases;
163 948 : for ( iOld = 0; iOld < *pNumIndexes; iOld++ )
164 : {
165 : /* We don't want that the old peak index is at the border of the new peak region, that is why >= newUpperIndex and > newLowerIndex */
166 1380 : while ( ( iNew < newNumIndexes ) && ( indexOfTonalPeak[iOld] >= newUpperIndex[iNew] ) )
167 : {
168 567 : ++iNew;
169 : }
170 :
171 813 : if ( ( iNew < newNumIndexes ) && ( indexOfTonalPeak[iOld] > newLowerIndex[iNew] ) )
172 : {
173 597 : newIndexOfTonalPeak[nPreservedPeaks] = indexOfTonalPeak[iOld];
174 597 : newLowerIndex[nPreservedPeaks] = lowerIndex[iOld];
175 597 : newUpperIndex[nPreservedPeaks] = upperIndex[iOld];
176 597 : phaseDiff[nPreservedPeaks] = phaseDiff[iOld];
177 :
178 4776 : for ( j = lowerIndex[iOld]; j <= upperIndex[iOld]; j++ )
179 : {
180 4179 : *pNewPhase++ = *pOldPhase++;
181 : }
182 :
183 597 : ++nPreservedPeaks;
184 : }
185 : else
186 : {
187 216 : pOldPhase += upperIndex[iOld] - lowerIndex[iOld] + 1;
188 : }
189 : }
190 :
191 732 : for ( iNew = 0; iNew < nPreservedPeaks; iNew++ )
192 : {
193 597 : indexOfTonalPeak[iNew] = newIndexOfTonalPeak[iNew];
194 597 : lowerIndex[iNew] = newLowerIndex[iNew];
195 597 : upperIndex[iNew] = newUpperIndex[iNew];
196 : }
197 :
198 135 : *pNumIndexes = nPreservedPeaks;
199 :
200 135 : return;
201 : }
202 :
203 :
204 : /*-------------------------------------------------------------------*
205 : * Local functions
206 : *-------------------------------------------------------------------*/
207 :
208 6066 : static void calcPseudoSpec(
209 : const float *mdctSpec,
210 : const uint16_t nSamples,
211 : float floorPowerSpectrum,
212 : float *powerSpec )
213 : {
214 : int16_t k;
215 : float x;
216 :
217 3019473 : for ( k = 1; k <= nSamples - 2; k++ )
218 : {
219 3013407 : x = ( mdctSpec[k + 1] - mdctSpec[k - 1] ); /* An MDST estimate */
220 3013407 : x = mdctSpec[k] * mdctSpec[k] + x * x;
221 3013407 : powerSpec[k] = max( floorPowerSpectrum, x );
222 : }
223 :
224 6066 : powerSpec[0] = 0.5f * powerSpec[1];
225 6066 : powerSpec[nSamples - 1] = 0.5f * powerSpec[nSamples - 2];
226 :
227 6066 : return;
228 : }
229 :
230 12132 : static void getEnvelope(
231 : const int16_t nSamples,
232 : const float *powerSpec,
233 : float F0,
234 : float *envelope,
235 : float *smoothedSpectrum )
236 : {
237 : int16_t nFilterLength, nHalfFilterLength, nSecondHalfFilterLength, n1, n2;
238 : int16_t i;
239 : double sum; /* double is required to avoid precision problems in the optimized version. */
240 :
241 12132 : if ( F0 == 0 )
242 : {
243 10698 : nFilterLength = 15;
244 : }
245 1434 : else if ( F0 <= 10.0f )
246 : {
247 540 : nFilterLength = 11;
248 : }
249 894 : else if ( F0 >= 22.0f )
250 : {
251 : /* For F0 >= 22 peak is isolated well enough with the filter length of 23.
252 : This case is however not triggered due to the limit of pit_min,
253 : but the line is left for security reasons. */
254 21 : nFilterLength = 23;
255 : }
256 : else
257 : {
258 873 : nFilterLength = 1 + 2 * (int16_t) ( F0 / 2 );
259 : }
260 :
261 12132 : nHalfFilterLength = nFilterLength / 2;
262 12132 : n1 = nHalfFilterLength + 1;
263 12132 : nSecondHalfFilterLength = nFilterLength - nHalfFilterLength;
264 12132 : n2 = nSecondHalfFilterLength - 1;
265 :
266 12132 : assert( ( nFilterLength >= 7 ) && ( nFilterLength <= 23 ) && ( nFilterLength % 2 == 1 ) );
267 :
268 12132 : sum = 0;
269 95625 : for ( i = 0; i < n2; i++ )
270 : {
271 83493 : sum += LEVEL_ABOVE_ENVELOPE * powerSpec[i];
272 : }
273 :
274 : /* No need for PTR_INIT for powerSpec[i+n2] as we continue from the previous loop */
275 107757 : for ( i = 0; i < n1; i++ )
276 : {
277 95625 : sum += LEVEL_ABOVE_ENVELOPE * powerSpec[i + n2];
278 : /* 1/(i+nSecondHalfFilterLength) needs to be stored in a table for each filter length */
279 95625 : envelope[i] = (float) sum / ( i + nSecondHalfFilterLength );
280 : }
281 12132 : sum /= nFilterLength; /* 1/nFilterLength is from a table */
282 :
283 : /* No need for PTR_INIT for powerSpec[i+n2] as we continue from the previous loop */
284 5884092 : for ( i = n1; i < nSamples - n2; i++ )
285 : {
286 5871960 : sum += LEVEL_ABOVE_ENVELOPE * ( powerSpec[i + n2] - powerSpec[i - n1] ) / nFilterLength; /* LEVEL_ABOVE_ENVELOPE/nFilterLength is from a table */
287 5871960 : envelope[i] = (float) sum;
288 : }
289 12132 : sum *= nFilterLength;
290 :
291 : /* No need for PTR_INIT for powerSpec[i-n1] as we continue from the previous loop */
292 95625 : for ( i = nSamples - n2; i < nSamples; i++ )
293 : {
294 83493 : sum -= LEVEL_ABOVE_ENVELOPE * powerSpec[i - n1];
295 : /* 1/(nSamples-(i-nHalfFilterLength)) needs to be stored in a table for each filter length */
296 83493 : envelope[i] = (float) sum / ( nSamples - ( i - nHalfFilterLength ) );
297 : }
298 :
299 6038946 : for ( i = 1; i < nSamples - 1; i++ )
300 : {
301 6026814 : smoothedSpectrum[i] = 0.75f * powerSpec[i - 1] + powerSpec[i] + 0.75f * powerSpec[i + 1];
302 : }
303 :
304 12132 : smoothedSpectrum[0] = powerSpec[0] + 0.75f * powerSpec[1];
305 12132 : smoothedSpectrum[nSamples - 1] = 0.75f * powerSpec[nSamples - 2] + powerSpec[nSamples - 1];
306 :
307 12132 : return;
308 : }
309 :
310 6066 : static void GetF0(
311 : const uint16_t nSamples, /* i*/
312 : const uint16_t nSamplesCore, /* i*/
313 : const float *powerSpectrum, /* i*/
314 : const float pitchLag, /* i*/
315 : float *pOrigF0, /* i/o*/
316 : float *pF0 /* i/o*/
317 : )
318 : {
319 : float halfPitchLag;
320 : uint16_t rgiStrongHarmonics[MAX_PEAKS_FROM_PITCH];
321 : int16_t nTotalHarmonics, nStrongHarmonics;
322 :
323 6066 : assert( LAST_HARMONIC_POS_TO_CHECK <= nSamplesCore );
324 :
325 : /* Use only F0 >= 100 Hz */
326 6066 : if ( ( pitchLag > 0 ) && ( pitchLag <= 0.5f * nSamplesCore ) )
327 : {
328 1497 : halfPitchLag = 0.5f * pitchLag;
329 1497 : *pF0 = nSamplesCore / halfPitchLag;
330 1497 : *pOrigF0 = *pF0;
331 1497 : if ( nSamples < 2 * LAST_HARMONIC_POS_TO_CHECK )
332 : {
333 108 : nTotalHarmonics = (int16_t) ( nSamples / ( *pF0 ) );
334 : }
335 : else
336 : {
337 1389 : nTotalHarmonics = (int16_t) ( 2 * LAST_HARMONIC_POS_TO_CHECK / ( *pF0 ) ); /* For correcting F0 we go 2 times the last harmonic position that will be used */
338 : }
339 :
340 : /* Get in rgiStrongHarmonics all i for which i*F0 are the strongest harmonics */
341 1497 : findStrongestHarmonics( nSamples, powerSpectrum, *pF0, nTotalHarmonics, rgiStrongHarmonics, &nStrongHarmonics );
342 1497 : CorrectF0( rgiStrongHarmonics, nStrongHarmonics, pF0 );
343 : }
344 : else
345 : {
346 4569 : *pF0 = 0;
347 4569 : *pOrigF0 = 0;
348 : }
349 :
350 6066 : return;
351 : }
352 :
353 1497 : static void findStrongestHarmonics(
354 : const uint16_t nSamples,
355 : const float *powerSpectrum,
356 : const float F0,
357 : const uint16_t nTotalHarmonics,
358 : uint16_t *pHarmonicIndexes,
359 : int16_t *pnHarmonics )
360 : {
361 : float peaks[MAX_PEAKS_FROM_PITCH], smallestPeak;
362 : uint16_t nPeaksToCheck, nPeaks, iSmallestPeak;
363 : int16_t i, l, k;
364 :
365 1497 : nPeaks = 0;
366 1497 : iSmallestPeak = 0;
367 1497 : smallestPeak = FLT_MAX;
368 1497 : nPeaksToCheck = min( nTotalHarmonics, MAX_PEAKS_FROM_PITCH + 1 );
369 :
370 16443 : for ( i = 1; i < nPeaksToCheck; i++ )
371 : {
372 : float newPeak;
373 :
374 14946 : k = (int16_t) ( i * F0 );
375 14946 : assert( k > 0 && k < 2 * LAST_HARMONIC_POS_TO_CHECK && k < nSamples );
376 :
377 14946 : newPeak = powerSpectrum[k];
378 14946 : peaks[nPeaks] = newPeak;
379 14946 : pHarmonicIndexes[nPeaks] = i;
380 14946 : if ( newPeak <= smallestPeak )
381 : {
382 7626 : iSmallestPeak = nPeaks;
383 7626 : smallestPeak = newPeak;
384 : }
385 14946 : ++nPeaks;
386 : }
387 :
388 25413 : for ( ; i < nTotalHarmonics; i++ )
389 : {
390 : float newPeak;
391 :
392 23916 : k = (int16_t) ( i * F0 );
393 23916 : assert( k > 0 && k < 2 * LAST_HARMONIC_POS_TO_CHECK && k < nSamples );
394 23916 : newPeak = powerSpectrum[k];
395 :
396 23916 : if ( newPeak > smallestPeak )
397 : {
398 3717 : peaks[iSmallestPeak] = newPeak;
399 3717 : pHarmonicIndexes[iSmallestPeak] = i;
400 3717 : smallestPeak = newPeak;
401 40887 : for ( l = 0; l < MAX_PEAKS_FROM_PITCH; l++ )
402 : {
403 37170 : if ( peaks[l] <= smallestPeak )
404 : {
405 7578 : iSmallestPeak = l;
406 7578 : smallestPeak = peaks[l];
407 : }
408 : }
409 : }
410 : }
411 1497 : sort( pHarmonicIndexes, nPeaks );
412 1497 : *pnHarmonics = nPeaks;
413 :
414 1497 : return;
415 : }
416 :
417 : /* Use new F0, for which harmonics are most common in pHarmonicIndexes */
418 1497 : static void CorrectF0(
419 : const uint16_t *pHarmonicIndexes,
420 : const uint16_t nHarmonics,
421 : float *pF0 )
422 : {
423 : int16_t i;
424 : float F0;
425 : uint16_t diff[MAX_PEAKS_FROM_PITCH - 1], sortedDiff[MAX_PEAKS_FROM_PITCH - 1];
426 : uint16_t iMostCommonDiff, nMostCommonDiff, nSameDiff, iMult;
427 :
428 14970 : for ( i = 0; i < MAX_PEAKS_FROM_PITCH - 1; i++ )
429 : {
430 13473 : diff[i] = 0;
431 13473 : sortedDiff[i] = 0;
432 : }
433 :
434 1497 : F0 = *pF0;
435 1497 : if ( F0 > 0 && nHarmonics > 0 )
436 : {
437 14946 : for ( i = 0; i < nHarmonics - 1; i++ )
438 : {
439 13449 : diff[i] = pHarmonicIndexes[i + 1] - pHarmonicIndexes[i];
440 13449 : sortedDiff[i] = diff[i];
441 : }
442 1497 : sort( sortedDiff, nHarmonics - 1 );
443 1497 : iMostCommonDiff = sortedDiff[0];
444 1497 : nSameDiff = 1;
445 1497 : i = 1;
446 1497 : if ( sortedDiff[0] * pHarmonicIndexes[0] == 1 )
447 : {
448 : /* Find how many distances between peaks have length 1 */
449 11640 : for ( ; i < nHarmonics - 1; i++ )
450 : {
451 10344 : if ( sortedDiff[i] == 1 )
452 : {
453 8127 : ++nSameDiff;
454 : }
455 : }
456 : }
457 1497 : nMostCommonDiff = nSameDiff;
458 :
459 : /* If there are at least 3 distances between peaks with length 1 and if the 1st harmonic is in pHarmonicIndexes then keep the original F0 */
460 : /* Otherwise find the most common distance between peaks */
461 1497 : if ( nSameDiff < 3 )
462 : {
463 : /* Find the most common difference */
464 2118 : for ( i = nSameDiff; i < nHarmonics - 1; i++ )
465 : {
466 1878 : if ( sortedDiff[i] == sortedDiff[i - 1] )
467 : {
468 1350 : ++nSameDiff;
469 : }
470 : else
471 : {
472 528 : if ( nSameDiff > nMostCommonDiff )
473 : {
474 240 : nMostCommonDiff = nSameDiff;
475 240 : iMostCommonDiff = sortedDiff[i - 1];
476 : }
477 : else
478 : {
479 288 : if ( ( nSameDiff == nMostCommonDiff ) && ( abs( (int16_t) iMostCommonDiff - (int16_t) pHarmonicIndexes[0] ) > abs( (int16_t) sortedDiff[i - 1] - (int16_t) pHarmonicIndexes[0] ) ) )
480 : {
481 18 : nMostCommonDiff = nSameDiff;
482 18 : iMostCommonDiff = sortedDiff[i - 1];
483 : }
484 : }
485 528 : nSameDiff = 1;
486 : }
487 : }
488 :
489 240 : if ( nSameDiff > nMostCommonDiff )
490 : {
491 69 : nMostCommonDiff = nSameDiff;
492 69 : iMostCommonDiff = sortedDiff[nHarmonics - 2];
493 : }
494 : }
495 :
496 : /* If there are enough peaks at the same distance */
497 1497 : if ( nMostCommonDiff >= MAX_PEAKS_FROM_PITCH / 2 )
498 : {
499 1368 : iMult = 1;
500 1407 : for ( i = 0; i < nHarmonics - 1; i++ )
501 : {
502 1407 : if ( diff[i] == iMostCommonDiff )
503 : {
504 1311 : iMult = pHarmonicIndexes[i];
505 1311 : break;
506 : }
507 :
508 : /* for rare cases of octave mismatch or missing harmonics */
509 96 : if ( ( i < nHarmonics - 2 ) && ( diff[i] == diff[i + 1] ) && ( diff[i] + diff[i + 1] == iMostCommonDiff ) )
510 : {
511 57 : iMult = pHarmonicIndexes[i];
512 57 : break;
513 : }
514 : }
515 :
516 : /* If the real F0 is much higher than the original F0 from the pitch */
517 1368 : if ( iMult <= 3 )
518 : {
519 : /* Use iMostCommonDiff, because the lowest pHarmonicIndexes[i] (which is equal to iMult) may not correspond to the new F0, but to it's multiple */
520 1353 : F0 = iMostCommonDiff * F0;
521 : }
522 : else
523 : {
524 15 : F0 = 0;
525 : }
526 : }
527 : /* Otherwise if there are at least 3 distances between peaks with length 1 and if the 1st harmonic is in pHarmonicIndexes then keep the original F0 */
528 : /* Otherwise don't use F0 */
529 129 : else if ( ( iMostCommonDiff > 1 ) || ( nMostCommonDiff < 3 ) )
530 : {
531 : /* Not enough peaks at the same distance => don't use the pitch. */
532 48 : F0 = 0;
533 : }
534 1497 : *pF0 = F0;
535 : }
536 :
537 1497 : return;
538 : }
539 :
540 15168 : static void modifyThreshold(
541 : const int16_t i,
542 : float F0,
543 : float threshold,
544 : float *thresholdModification )
545 : {
546 : float harmonic, fractional, twoTimesFract;
547 : int16_t k;
548 :
549 15168 : harmonic = i * F0;
550 15168 : k = (int16_t) harmonic;
551 15168 : fractional = harmonic - k; /* Fractional part of the i*F0 */
552 15168 : twoTimesFract = 2.0f * fractional; /* threshold if the centar of the peek is between k-1 and k, threshold+2 if the centar of the peek is between k and k+1 */
553 15168 : thresholdModification[k] = threshold;
554 15168 : thresholdModification[k - 1] = threshold + twoTimesFract;
555 15168 : thresholdModification[k + 1] = threshold + 2.0f - twoTimesFract;
556 :
557 15168 : return;
558 : }
559 :
560 6066 : static void modifyThresholds(
561 : float F0,
562 : float origF0,
563 : float *thresholdModification )
564 : {
565 : int16_t i, nHarmonics;
566 :
567 6066 : if ( ( F0 == 0 ) && ( origF0 > 0 ) )
568 : {
569 63 : nHarmonics = min( MAX_PEAKS_FROM_PITCH, (int16_t) ( LAST_HARMONIC_POS_TO_CHECK / origF0 ) );
570 693 : for ( i = 1; i <= nHarmonics; i++ )
571 : {
572 630 : modifyThreshold( i, origF0, 0.7f, thresholdModification );
573 : }
574 : }
575 6003 : else if ( ( F0 > 0 ) && ( origF0 > 0 ) )
576 : {
577 1434 : nHarmonics = min( MAX_PEAKS_FROM_PITCH, (int16_t) ( LAST_HARMONIC_POS_TO_CHECK / F0 ) );
578 :
579 3003 : for ( i = (int16_t) ( F0 / origF0 + 0.5f ); i > 0; i-- )
580 : {
581 1569 : modifyThreshold( i, origF0, 0.35f, thresholdModification );
582 : }
583 14403 : for ( i = 1; i <= nHarmonics; i++ )
584 : {
585 12969 : modifyThreshold( i, F0, 0.35f, thresholdModification );
586 : }
587 : }
588 :
589 6066 : return;
590 : }
591 :
592 6066 : static void findCandidates(
593 : const uint16_t nSamples,
594 : const float *MDCTSpectrum,
595 : float *thresholdModificationNew,
596 : float floorPowerSpectrum /* i : lower limit for power spectrum bins */
597 : )
598 : {
599 : float powerSpectrum[L_FRAME_MAX];
600 : float envelope[L_FRAME_MAX];
601 : float smoothedSpectrum[L_FRAME_MAX];
602 : uint16_t upperIdx, lowerIdx;
603 : int16_t k, j;
604 :
605 6066 : calcPseudoSpec( MDCTSpectrum, nSamples, floorPowerSpectrum, powerSpectrum );
606 6066 : getEnvelope( nSamples, powerSpectrum, 0, envelope, smoothedSpectrum );
607 :
608 6066 : set_f( thresholdModificationNew, UNREACHABLE_THRESHOLD, nSamples );
609 2840973 : for ( k = GROUP_LENGTH / 2; k <= nSamples - ( GROUP_LENGTH - GROUP_LENGTH / 2 ); k++ )
610 : {
611 2834907 : if ( smoothedSpectrum[k] > envelope[k] )
612 : {
613 : /* The check that bin at k is bigger than bins at k-1 and k+1 is needed to avoid deadlocks when the thresholds are low. */
614 : /* It removes some true peaks, especially if non weighted sum is used for the smoothed spectrum. */
615 : float biggerNeighbor;
616 27627 : biggerNeighbor = max( powerSpectrum[k - 1], powerSpectrum[k + 1] );
617 27627 : if ( powerSpectrum[k] >= biggerNeighbor )
618 : {
619 : /* Find the right foot */
620 152778 : for ( upperIdx = k + 1; upperIdx < nSamples - 1; upperIdx++ )
621 : {
622 152610 : if ( powerSpectrum[upperIdx] < powerSpectrum[upperIdx + 1] )
623 : {
624 : /* Side lobes may increase for certain ammount */
625 17874 : if ( ALLOWED_SIDE_LOBE_FLUCTUATION * powerSpectrum[upperIdx] < powerSpectrum[upperIdx + 1] )
626 : {
627 9051 : break;
628 : }
629 : /* Check for further decrease after a side lobe increase */
630 11124 : for ( j = upperIdx + 1; j < nSamples - 1; j++ )
631 : {
632 11124 : if ( powerSpectrum[j] < ALLOWED_SIDE_LOBE_FLUCTUATION * powerSpectrum[j + 1] )
633 : {
634 8823 : break;
635 : }
636 : }
637 : /* Side lobe increase must be 2 times smaller than the decrease to the foot */
638 : /* Eq. to 2.0f*powerSpectrum[lowerIdx-1]/powerSpectrum[lowerIdx] > powerSpectrum[lowerIdx]/powerSpectrum[j] */
639 8823 : if ( 2.0f * powerSpectrum[upperIdx + 1] * powerSpectrum[j] > powerSpectrum[upperIdx] * powerSpectrum[upperIdx] )
640 : {
641 7191 : break;
642 : }
643 1632 : upperIdx = j - 1; /* Reinitialize pointers due to the change of upperIdx */
644 : }
645 : }
646 : /* left foot */
647 115335 : for ( lowerIdx = k - 1; lowerIdx > 0; lowerIdx-- )
648 : {
649 114993 : if ( powerSpectrum[lowerIdx] < powerSpectrum[lowerIdx - 1] )
650 : {
651 : /* Side lobes may increase for certain ammount */
652 17328 : if ( ALLOWED_SIDE_LOBE_FLUCTUATION * powerSpectrum[lowerIdx] < powerSpectrum[lowerIdx - 1] )
653 : {
654 9045 : break;
655 : }
656 : /* Check for further decrease after a side lobe increase */
657 10212 : for ( j = lowerIdx - 1; j > 0; j-- )
658 : {
659 10212 : if ( powerSpectrum[j] < ALLOWED_SIDE_LOBE_FLUCTUATION * powerSpectrum[j - 1] )
660 : {
661 8283 : break;
662 : }
663 : }
664 : /* Side lobe increase must be 2 times smaller than the decrease to the foot */
665 : /* Eq. to 2.0f*powerSpectrum[lowerIdx-1]/powerSpectrum[lowerIdx] > powerSpectrum[lowerIdx]/powerSpectrum[j] */
666 8283 : if ( 2.0f * powerSpectrum[lowerIdx - 1] * powerSpectrum[j] > powerSpectrum[lowerIdx] * powerSpectrum[lowerIdx] )
667 : {
668 7023 : break;
669 : }
670 1260 : lowerIdx = j + 1; /* Reinitialize pointers due to the change of lowerIdx */
671 : }
672 : }
673 :
674 : /* Check if there is a bigger peak up to the next peak foot */
675 302796 : for ( j = max( GROUP_LENGTH / 2, lowerIdx ); j <= min( upperIdx, nSamples - ( GROUP_LENGTH - GROUP_LENGTH / 2 ) ); j++ )
676 : {
677 286386 : if ( powerSpectrum[j] > powerSpectrum[k] )
678 : {
679 153 : k = j; /* PTR_INIT for powerSpectrum[k] */
680 : }
681 : }
682 :
683 : /* Modify thresholds for the following frame */
684 65640 : for ( j = k - 1; j < k + 2; j++ )
685 : {
686 49230 : if ( smoothedSpectrum[j] > envelope[j] )
687 : {
688 32034 : thresholdModificationNew[j] = SMALL_THRESHOLD;
689 : }
690 : else
691 : {
692 17196 : thresholdModificationNew[j] = BIG_THRESHOLD;
693 : }
694 : }
695 : /* Jump to the next foot of the peak. */
696 16410 : k = upperIdx; /* Reinitialize pointers due to the change of k */
697 : }
698 : }
699 : }
700 :
701 6066 : return;
702 : }
703 :
704 6066 : static void RefineThresholdsUsingPitch(
705 : const uint16_t nSamples,
706 : const uint16_t nSamplesCore,
707 : const float *powerSpectrum,
708 : const float lastPitchLag,
709 : const float currentPitchLag,
710 : float *pF0,
711 : float *thresholdModification )
712 : {
713 : int16_t pitchIsStable;
714 : float origF0;
715 :
716 6066 : pitchIsStable = ( fabs( lastPitchLag - currentPitchLag ) < 0.25f );
717 6066 : if ( pitchIsStable )
718 : {
719 6066 : GetF0( nSamples, nSamplesCore, powerSpectrum, lastPitchLag, &origF0, pF0 );
720 :
721 6066 : modifyThresholds( *pF0, origF0, thresholdModification );
722 : }
723 : else
724 : {
725 0 : *pF0 = 0;
726 : }
727 :
728 6066 : return;
729 : }
730 :
731 6066 : static void findTonalComponents(
732 : uint16_t *indexOfTonalPeak, /* OUT */
733 : uint16_t *lowerIndex, /* OUT */
734 : uint16_t *upperIndex, /* OUT */
735 : uint16_t *numIndexes, /* OUT */
736 : const uint16_t nSamples, /* IN */
737 : const float *powerSpectrum, /* IN */
738 : const float F0, /* IN */
739 : const float *thresholdModification /* IN */
740 : )
741 : {
742 : float envelope[L_FRAME_MAX];
743 : float smoothedSpectrum[L_FRAME_MAX];
744 :
745 : uint16_t nrOfFIS;
746 : uint16_t upperIdx, lowerIdx, lowerBound;
747 : int16_t k, j;
748 :
749 6066 : getEnvelope( nSamples, powerSpectrum, F0, envelope, smoothedSpectrum );
750 6066 : nrOfFIS = 0;
751 6066 : lowerBound = 0;
752 2872104 : for ( k = GROUP_LENGTH / 2; k <= nSamples - ( GROUP_LENGTH - GROUP_LENGTH / 2 ); k++ )
753 : {
754 2866038 : if ( smoothedSpectrum[k] > envelope[k] * thresholdModification[k] )
755 : {
756 : /* The check that bin at k is bigger than bins at k-1 and k+1 is needed to avoid deadlocks when the thresholds are low. */
757 : /* It removes some true peaks, especially if non weighted sum is used for the smoothed spectrum. */
758 : float biggerNeighbor;
759 18033 : biggerNeighbor = max( powerSpectrum[k - 1], powerSpectrum[k + 1] );
760 18033 : if ( powerSpectrum[k] >= biggerNeighbor )
761 : {
762 : /* Find the right foot */
763 121812 : for ( upperIdx = k + 1; upperIdx < nSamples - 1; upperIdx++ )
764 : {
765 121716 : if ( powerSpectrum[upperIdx] < powerSpectrum[upperIdx + 1] )
766 : {
767 : /* Side lobes may increase for certain ammount */
768 13728 : if ( ALLOWED_SIDE_LOBE_FLUCTUATION * powerSpectrum[upperIdx] < powerSpectrum[upperIdx + 1] )
769 : {
770 5655 : break;
771 : }
772 : /* Check for further decrease after a side lobe increase */
773 10215 : for ( j = upperIdx + 1; j < nSamples - 1; j++ )
774 : {
775 10215 : if ( powerSpectrum[j] < ALLOWED_SIDE_LOBE_FLUCTUATION * powerSpectrum[j + 1] )
776 : {
777 8073 : break;
778 : }
779 : }
780 : /* Side lobe increase must be 2 times smaller than the decrease to the foot */
781 : /* Eq. to 2.0f*powerSpectrum[lowerIdx-1]/powerSpectrum[lowerIdx] > powerSpectrum[lowerIdx]/powerSpectrum[j] */
782 8073 : if ( 2.0f * powerSpectrum[upperIdx + 1] * powerSpectrum[j] > powerSpectrum[upperIdx] * powerSpectrum[upperIdx] )
783 : {
784 6813 : break;
785 : }
786 1260 : upperIdx = j - 1; /* Reinitialize pointers due to the change of upperIdx */
787 : }
788 : }
789 : /* left foot */
790 53967 : for ( lowerIdx = k - 1; lowerIdx > lowerBound; lowerIdx-- )
791 : {
792 52149 : if ( powerSpectrum[lowerIdx] < powerSpectrum[lowerIdx - 1] )
793 : {
794 : /* Side lobes may increase for certain ammount */
795 12030 : if ( ALLOWED_SIDE_LOBE_FLUCTUATION * powerSpectrum[lowerIdx] < powerSpectrum[lowerIdx - 1] )
796 : {
797 4620 : break;
798 : }
799 : /* Check for further decrease after a side lobe increase */
800 9531 : for ( j = lowerIdx - 1; j > 0; j-- )
801 : {
802 9531 : if ( powerSpectrum[j] < ALLOWED_SIDE_LOBE_FLUCTUATION * powerSpectrum[j - 1] )
803 : {
804 7410 : break;
805 : }
806 : }
807 : /* Side lobe increase must be 2 times smaller than the decrease to the foot */
808 : /* Eq. to 2.0f*powerSpectrum[lowerIdx-1]/powerSpectrum[lowerIdx] > powerSpectrum[lowerIdx]/powerSpectrum[j] */
809 7410 : if ( 2.0f * powerSpectrum[lowerIdx - 1] * powerSpectrum[j] > powerSpectrum[lowerIdx] * powerSpectrum[lowerIdx] )
810 : {
811 6126 : break;
812 : }
813 1284 : lowerIdx = j + 1; /* Reinitialize pointers due to the change of lowerIdx */
814 : }
815 : }
816 12564 : lowerBound = upperIdx;
817 :
818 : /* Check if there is a bigger peak up to the next peak foot */
819 203193 : for ( j = max( GROUP_LENGTH / 2, lowerIdx ); j <= min( upperIdx, nSamples - ( GROUP_LENGTH - GROUP_LENGTH / 2 ) ); j++ )
820 : {
821 190629 : if ( powerSpectrum[j] > powerSpectrum[k] )
822 : {
823 9 : k = j; /* PTR_INIT for powerSpectrum[k] */
824 : }
825 : }
826 :
827 12564 : assert( ( nrOfFIS == 0 ) || ( indexOfTonalPeak[nrOfFIS - 1] < k ) );
828 12564 : lowerIndex[nrOfFIS] = k - GROUP_LENGTH / 2;
829 12564 : upperIndex[nrOfFIS] = k + ( GROUP_LENGTH - GROUP_LENGTH / 2 - 1 );
830 12564 : if ( ( nrOfFIS > 0 ) && ( lowerIndex[nrOfFIS] <= upperIndex[nrOfFIS - 1] ) )
831 : {
832 84 : int16_t m = ( k + indexOfTonalPeak[nrOfFIS - 1] ) / 2;
833 84 : upperIndex[nrOfFIS - 1] = m;
834 84 : lowerIndex[nrOfFIS] = m + 1;
835 : }
836 12564 : indexOfTonalPeak[nrOfFIS++] = k;
837 12564 : if ( nrOfFIS == MAX_NUMBER_OF_IDX )
838 : {
839 0 : break;
840 : }
841 : /* Jump to the next foot of the peak. */
842 12564 : k = upperIdx; /* Reinitialize pointers due to the change of k */
843 : }
844 : }
845 : }
846 6066 : *numIndexes = nrOfFIS;
847 :
848 6066 : return;
849 : }
|
