1 /* LPC_and_Tube.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 20020612 GPL header
21 djmw 20041020 struct Tube_Frame -> struct structTube_Frame; struct LPC_Frame -> struct structLPC_Frame;
22 struct Formant_Frame->struct structFormant_Frame
23 djmw 20051005 Always make a VocalTract with length 0.01 m when isundef(wakita_length).
24 */
25
26 #include "LPC_and_Tube.h"
27 #include "LPC_and_Formant.h"
28 #include "LPC_to_Spectrum.h"
29 #include "SpectrumTier.h"
30 #include "VocalTract_to_Spectrum.h"
31 #include "NUM2.h"
32
33 // IEEE: Programs fo digital signal processing section 4.3 LPTRN
34
LPC_Frame_into_Tube_Frame_rc(LPC_Frame me,Tube_Frame thee)35 void LPC_Frame_into_Tube_Frame_rc (LPC_Frame me, Tube_Frame thee) {
36 Melder_assert (my nCoefficients == my a.size); // check invariant
37 thy c.resize (my nCoefficients);
38 thy numberOfSegments = thy c.size; // maintain invariant
39 VECrc_from_lpc (thy c.get(), my a.get());
40 }
41
LPC_Frame_into_Tube_Frame_area(LPC_Frame me,Tube_Frame thee)42 void LPC_Frame_into_Tube_Frame_area (LPC_Frame me, Tube_Frame thee) {
43 VECarea_from_lpc (thy c.part (1, my nCoefficients), my a.part (1, my nCoefficients));
44 }
45
VocalTract_LPC_Frame_getMatchingLength(VocalTract me,LPC_Frame thee,double glottalDamping,bool radiationDamping,bool internalDamping)46 double VocalTract_LPC_Frame_getMatchingLength (VocalTract me, LPC_Frame thee, double glottalDamping, bool radiationDamping, bool internalDamping) {
47 try {
48 /*
49 Match the average distance between the first two formants in the VocaTract and the LPC spectrum
50 */
51 constexpr integer numberOfFrequencies = 1000;
52 const double maximumFrequency = 5000.0;
53 autoSpectrum vts = VocalTract_to_Spectrum (me, numberOfFrequencies, maximumFrequency, glottalDamping, radiationDamping, internalDamping);
54 const double samplingFrequency = 1000.0 * my nx;
55 autoSpectrum lps = Spectrum_create (0.5 * samplingFrequency, numberOfFrequencies);
56 LPC_Frame_into_Spectrum (thee, lps.get(), 0, 50);
57 autoSpectrumTier vtst = Spectrum_to_SpectrumTier_peaks (vts.get());
58 autoSpectrumTier lpst = Spectrum_to_SpectrumTier_peaks (lps.get());
59 const double vt_f1 = vtst -> points.at [1] -> number, vt_f2 = vtst -> points.at [2] -> number;
60 const double lp_f1 = lpst -> points.at [1] -> number, lp_f2 = lpst -> points.at [2] -> number;
61 const double df1 = lp_f1 - vt_f1, df2 = lp_f2 - vt_f2, df = 0.5 * (df1 + df2);
62 const double dl = - df / lp_f2;
63 return my dx * my nx * (1 + dl);
64 } catch (MelderError) {
65 Melder_throw (U"Length could not be determined from VocalTract and LPC_Frame.");
66 }
67 }
68
LPC_Frame_getVTL_wakita(LPC_Frame me,double samplingPeriod,double refLength)69 double LPC_Frame_getVTL_wakita (LPC_Frame me, double samplingPeriod, double refLength) {
70 struct structLPC_Frame lpc_struct;
71 const LPC_Frame lpc = & lpc_struct;
72 struct structFormant_Frame f_struct;
73 const Formant_Frame f = & f_struct;
74 struct structTube_Frame rc_struct, af_struct;
75 const Tube_Frame rc = & rc_struct, af = & af_struct;
76 try {
77 const integer numberOfCoefficients = my nCoefficients;
78 const double dlength = 0.001;
79 double length, wakita_length = undefined;
80 double varMin = 1e308;
81
82 memset (& lpc_struct, 0, sizeof (lpc_struct));
83 memset (& f_struct, 0, sizeof (f_struct));
84 memset (& rc_struct, 0, sizeof (rc_struct));
85 memset (& af_struct, 0, sizeof (af_struct));
86
87 LPC_Frame_init (lpc, numberOfCoefficients);
88 Tube_Frame_init (rc, numberOfCoefficients, refLength);
89 Tube_Frame_init (af, numberOfCoefficients, refLength);
90 /*
91 Step 2
92 */
93 LPC_Frame_into_Formant_Frame (me, f, samplingPeriod, 0.0);
94 /*
95 LPC_Frame_into_Formant_Frame performs the Formant_Frame_init !!
96 */
97 Melder_require (f -> numberOfFormants > 0,
98 U"Not enough formants.");
99 VEC area = af -> c.get(); // TODO
100 double lmin = length = 0.10;
101 double plength = refLength;
102 while (length <= 0.25) {
103 /*
104 Step 3
105 */
106 const double fscale = plength / length;
107 for (integer i = 1; i <= f -> numberOfFormants; i ++) {
108 f -> formant [i]. frequency *= fscale;
109 f -> formant [i]. bandwidth *= fscale;
110 }
111 /*
112 20000125: Bovenstaande schaling van f1/b1 kan ook gedaan worden door
113 MGfb_to_a (f, b, nf, samplingFrequency*length/refLength, a1)
114 De berekening is extreem gevoelig voor de samplefrequentie: een zelfde
115 stel f,b waardes geven andere lengtes afhankelijk van Fs. Ook het
116 weglaten van een hogere formant heeft consekwenties.
117 De refLength zou eigenlijk vast moeten liggen op
118 refLength=c*aantalFormanten/Fs waarbij c=340 m/s (geluidssnelheid).
119 Bij Fs=10000 zou aantalFormanten=5 zijn en refLength -> 0.17 m
120 */
121 /*
122 Step 4
123 */
124 Formant_Frame_into_LPC_Frame (f, lpc, samplingPeriod);
125 /*
126 Step 5
127 */
128 rc -> length = length;
129 LPC_Frame_into_Tube_Frame_rc (lpc, rc);
130 /*
131 Step 6.1
132 */
133 Tube_Frames_rc_into_area (rc, af);
134 /*
135 Step 6.2 Log(areas)
136 */
137 double logSum = 0.0;
138 for (integer i = 1; i <= af -> numberOfSegments; i ++) {
139 area [i] = log (area [i]);
140 logSum += area [i];
141 }
142 /*
143 Steps 6.3 and 7
144 */
145 double var = 0.0;
146 for (integer i = 1; i <= af -> numberOfSegments; i ++) {
147 const double delta = area [i] - logSum / af -> numberOfSegments;
148 var += delta * delta;
149 }
150 if (var < varMin) {
151 lmin = length;
152 varMin = var;
153 }
154 plength = length;
155 length += dlength;
156 }
157
158 wakita_length = lmin;
159 f -> destroy ();
160 lpc -> destroy ();
161 rc -> destroy ();
162 af -> destroy ();
163 return wakita_length;
164 } catch (MelderError) {
165 f -> destroy ();
166 lpc -> destroy ();
167 rc -> destroy ();
168 af -> destroy ();
169 return undefined;
170 }
171 }
172
Tube_Frame_into_LPC_Frame_area(Tube_Frame me,LPC_Frame thee)173 int Tube_Frame_into_LPC_Frame_area (Tube_Frame me, LPC_Frame thee) {
174 (void) me;
175 (void) thee;
176 return 0;
177 }
178
Tube_Frame_into_LPC_Frame_rc(Tube_Frame me,LPC_Frame thee)179 int Tube_Frame_into_LPC_Frame_rc (Tube_Frame me, LPC_Frame thee) {
180 (void) me;
181 (void) thee;
182 return 0;
183 }
184
VocalTract_setLength(VocalTract me,double newLength)185 void VocalTract_setLength (VocalTract me, double newLength) {
186 my xmax = newLength;
187 my dx = newLength / my nx;
188 my x1 = 0.5 * my dx;
189 }
190
LPC_to_VocalTract_slice_special(LPC me,double time,double glottalDamping,bool radiationDamping,bool internalDamping)191 autoVocalTract LPC_to_VocalTract_slice_special (LPC me, double time, double glottalDamping, bool radiationDamping, bool internalDamping) {
192 try {
193 integer frameNumber = Sampled_xToNearestIndex (me, time);
194 Melder_clip (1_integer, & frameNumber, my nx); // constant extrapolation
195 LPC_Frame lpc = & my d_frames [frameNumber];
196 autoVocalTract thee = LPC_Frame_to_VocalTract (lpc, 0.17);
197 const double length = VocalTract_LPC_Frame_getMatchingLength (thee.get(), lpc, glottalDamping, radiationDamping, internalDamping);
198 VocalTract_setLength (thee.get(), length);
199 return thee;
200 } catch (MelderError) {
201 Melder_throw (me, U": no VocalTract created.");
202 }
203 }
204
LPC_Frame_to_VocalTract(LPC_Frame me,double length)205 autoVocalTract LPC_Frame_to_VocalTract (LPC_Frame me, double length) {
206 try {
207 const integer m = my nCoefficients;
208 autoVocalTract thee = VocalTract_create (m, length / m);
209 VECarea_from_lpc (thy z.row (1), my a.part (1, m));
210 // area [lips..glottis] (m^2) to VocalTract [glottis..lips] (m^2)
211
212 for (integer i = 1; i <= m / 2; i ++)
213 std::swap (thy z [1] [i], thy z [1] [m + 1 - i]);
214 return thee;
215 } catch (MelderError) {
216 Melder_throw (U"No VocalTract created from LPC_Frame.");
217 }
218 }
219
LPC_to_VocalTract_slice(LPC me,double time,double length)220 autoVocalTract LPC_to_VocalTract_slice (LPC me, double time, double length) {
221 try {
222 integer frameNumber = Sampled_xToNearestIndex (me, time);
223 Melder_clip (1_integer, & frameNumber, my nx); // constant extrapolation
224 const LPC_Frame lpc = & my d_frames [frameNumber];
225 autoVocalTract thee = LPC_Frame_to_VocalTract (lpc, length);
226 return thee;
227 } catch (MelderError) {
228 Melder_throw (me, U": no VocalTract created.");
229 }
230 }
231
232 /* End of file LPC_and_Tube.cpp */
233