1 /* SPINET_to_Pitch.cpp
2 *
3 * Copyright (C) 1993-2019 David Weenink
4 *
5 * This code is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or (at
8 * your option) any later version.
9 *
10 * This code is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this work. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 /*
20 djmw 19970408
21 djmw 20020813 GPL header
22 */
23
24 #include "SPINET_to_Pitch.h"
25 #include "Pitch_extensions.h"
26 #include "NUM2.h"
27
28 /*
29 from Erb-scale to log2-scale
30 */
31
SPINET_to_Pitch(SPINET me,double harmonicFallOffSlope,double ceiling,integer maxnCandidates)32 autoPitch SPINET_to_Pitch (SPINET me, double harmonicFallOffSlope, double ceiling, integer maxnCandidates) {
33 try {
34 constexpr integer nPointsPerOctave = 48;
35 const double fmin = NUMerbToHertz (SampledXY_indexToY (me, 1L));
36 const double fmax = NUMerbToHertz (SampledXY_indexToY (me, my ny));
37 const double fminl2 = NUMlog2 (fmin), fmaxl2 = NUMlog2 (fmax);
38 const double points = (fmaxl2 - fminl2) * nPointsPerOctave;
39 const double dfl2 = (fmaxl2 - fminl2) / (points - 1);
40 const integer numberOfFrequencyPoints = Melder_ifloor (points);
41 const integer maxHarmonic = Melder_ifloor (fmax / fmin);
42 const double unvoicedCriterium = 0.45;
43
44 Melder_require (numberOfFrequencyPoints > 1,
45 U"Frequency range too small.");
46 Melder_require (fmin < ceiling,
47 U"The centre frequency of the lowest filter should be smaller than the ceiling.");
48
49 autoPitch thee = Pitch_create (my xmin, my xmax, my nx, my dx, my x1, ceiling, maxnCandidates);
50 autoVEC power = raw_VEC (my nx);
51 autoVEC pitch = raw_VEC (numberOfFrequencyPoints);
52 autoVEC sumspec = raw_VEC (numberOfFrequencyPoints);
53 autoVEC y = raw_VEC (my ny);
54 autoVEC yv2 = raw_VEC (my ny);
55 autoVEC fl2 = raw_VEC (my ny);
56 /*
57 From ERB's to log (f)
58 */
59 for (integer i = 1; i <= my ny; i ++) {
60 const double f = NUMerbToHertz (my y1 + (i - 1) * my dy);
61 fl2 [i] = NUMlog2 (f);
62 }
63 /*
64 Determine global maximum power in frame
65 */
66 double maxStrength = 0.0, maxPower = 0.0;
67 for (integer j = 1; j <= my nx; j ++) {
68 const double p = NUMsum (my s.column (j));
69 if (p > maxPower)
70 maxPower = p;
71 power [j] = p;
72 }
73 Melder_require (maxPower != 0.0,
74 U"The sound should not have all amplitudes equal to zero.");
75
76 for (integer j = 1; j <= my nx; j ++) {
77 const Pitch_Frame pitchFrame = & thy frames [j];
78
79 pitchFrame -> intensity = power [j] / maxPower;
80 y.all() <<= my s.column (j);
81
82 NUMcubicSplineInterpolation_getSecondDerivatives (yv2.get(), fl2.get(), y.get(), 1e30, 1e30);
83 for (integer k = 1; k <= numberOfFrequencyPoints; k ++) {
84 const double f = fminl2 + (k - 1) * dfl2;
85 pitch [k] = NUMcubicSplineInterpolation (fl2.get(), y.get(), yv2.get(), f);
86 sumspec [k] = 0.0;
87 }
88 /*
89 Formula (8): weighted harmonic summation.
90 */
91 for (integer m = 1; m <= maxHarmonic; m ++) {
92 const double hm = 1 - harmonicFallOffSlope * NUMlog2 (m);
93 const integer kb = 1 + Melder_ifloor (nPointsPerOctave * NUMlog2 (m)); // TODO: what is kb?
94 for (integer k = kb; k <= numberOfFrequencyPoints; k ++)
95 if (pitch [k] > 0.0)
96 sumspec [k - kb + 1] += pitch [k] * hm;
97 }
98 /*
99 Into Pitch object
100 */
101 Pitch_Frame_init (pitchFrame, maxnCandidates);
102 pitchFrame -> candidates. resize (pitchFrame -> nCandidates = 0); // !!!!!
103 Pitch_Frame_addPitch (pitchFrame, 0, 0, maxnCandidates); // unvoiced
104
105 for (integer k = 2; k <= numberOfFrequencyPoints - 1; k ++) {
106 const double y1 = sumspec [k - 1], y2 = sumspec [k], y3 = sumspec [k + 1];
107 if (y2 > y1 && y2 >= y3) {
108 const double denum = y1 - 2.0 * y2 + y3, tmp = y3 - 4.0 * y2;
109 const double x = dfl2 * (y1 - y3) / (2 * denum);
110 const double f = pow (2.0, fminl2 + (k - 1) * dfl2 + x);
111 const double strength = (2.0 * y1 * (4.0 * y2 + y3) - y1 * y1 - tmp * tmp) / (8.0 * denum);
112 if (strength > maxStrength)
113 maxStrength = strength;
114 Pitch_Frame_addPitch (pitchFrame, f, strength, maxnCandidates);
115 }
116 }
117 }
118
119 // Scale the pitch strengths
120
121 for (integer j = 1; j <= my nx; j ++) {
122 double f0, localStrength;
123 Pitch_Frame_getPitch (& thy frames [j], & f0, & localStrength);
124 Pitch_Frame_resizeStrengths (& thy frames [j], localStrength / maxStrength, unvoicedCriterium);
125 }
126 return thee;
127 } catch (MelderError) {
128 Melder_throw (me, U": no Pitch created.");
129 }
130 }
131
132 /* End of file SPINET_to_Pitch.cpp */
133