1 /* LPC_and_Formant.cpp
2 *
3 * Copyright (C) 1994-2020 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 20030616 Formant_Frame_into_LPC_Frame: remove formant with f >= Nyquist +
21 change lpc indexing from -1..m
22 djmw 20080122 float -> double
23 */
24
25 #include "LPC_and_Formant.h"
26 #include "LPC_and_Polynomial.h"
27 #include "NUM2.h"
28 #include <thread>
29 #include <atomic>
30 #include <vector>
31
Formant_Frame_init(Formant_Frame me,integer numberOfFormants)32 void Formant_Frame_init (Formant_Frame me, integer numberOfFormants) {
33 if (numberOfFormants > 0)
34 my formant = newvectorzero <structFormant_Formant> (numberOfFormants);
35 my numberOfFormants = my formant.size; // maintain invariant
36 }
37
Formant_Frame_scale(Formant_Frame me,double scale)38 void Formant_Frame_scale (Formant_Frame me, double scale) {
39 for (integer iformant = 1; iformant <= my numberOfFormants; iformant ++) {
40 my formant [iformant]. frequency *= scale;
41 my formant [iformant]. bandwidth *= scale;
42 }
43 }
44
Roots_into_Formant_Frame(Roots me,Formant_Frame thee,double samplingFrequency,double margin)45 void Roots_into_Formant_Frame (Roots me, Formant_Frame thee, double samplingFrequency, double margin) {
46 /*
47 Determine the formants and bandwidths
48 */
49 Melder_assert (my numberOfRoots == my roots.size); // check invariant
50 thy formant.resize (0);
51 const double nyquistFrequency = 0.5 * samplingFrequency;
52 const double fLow = margin, fHigh = nyquistFrequency - margin;
53 for (integer iroot = 1; iroot <= my numberOfRoots; iroot ++) {
54 if (my roots [iroot].imag() < 0.0)
55 continue;
56 const double frequency = fabs (atan2 (my roots [iroot].imag(), my roots [iroot].real())) * nyquistFrequency / NUMpi;
57 if (frequency >= fLow && frequency <= fHigh) {
58 const double bandwidth = - log (norm (my roots [iroot])) * nyquistFrequency / NUMpi;
59 Formant_Formant newff = thy formant . append ();
60 newff -> frequency = frequency;
61 newff -> bandwidth = bandwidth;
62 }
63 }
64 thy numberOfFormants = thy formant.size; // maintain invariant
65 }
66
LPC_Frame_into_Formant_Frame(LPC_Frame me,Formant_Frame thee,double samplingPeriod,double margin)67 void LPC_Frame_into_Formant_Frame (LPC_Frame me, Formant_Frame thee, double samplingPeriod, double margin) {
68 Melder_assert (my nCoefficients == my a.size); // check invariant
69 thy intensity = my gain;
70 if (my nCoefficients == 0) {
71 thy formant.resize (0);
72 thy numberOfFormants = thy formant.size; // maintain invariant
73 return;
74 }
75 autoPolynomial p = LPC_Frame_to_Polynomial (me);
76 autoRoots r = Polynomial_to_Roots (p.get());
77 Roots_fixIntoUnitCircle (r.get());
78 Roots_into_Formant_Frame (r.get(), thee, 1.0 / samplingPeriod, margin);
79 }
80
LPC_Frame_into_Formant_Frame_mt(LPC_Frame me,Formant_Frame thee,double samplingPeriod,double margin,Polynomial p,Roots r,VEC const & workspace)81 void LPC_Frame_into_Formant_Frame_mt (LPC_Frame me, Formant_Frame thee, double samplingPeriod, double margin, Polynomial p, Roots r, VEC const& workspace) {
82 Melder_assert (my nCoefficients == my a.size); // check invariant
83 thy intensity = my gain;
84 if (my nCoefficients == 0) {
85 thy formant.resize (0);
86 thy numberOfFormants = thy formant.size; // maintain invariant
87 return;
88 }
89 LPC_Frame_into_Polynomial (me, p);
90 Polynomial_into_Roots (p, r, workspace);
91 Roots_fixIntoUnitCircle (r);
92 Roots_into_Formant_Frame (r, thee, 1.0 / samplingPeriod, margin);
93 }
94
LPC_to_Formant_noThreads(LPC me,double margin)95 autoFormant LPC_to_Formant_noThreads (LPC me, double margin) {
96 try {
97 const double samplingFrequency = 1.0 / my samplingPeriod;
98 /*
99 In very extreme case all roots of the lpc polynomial might be real.
100 A real root gives either a frequency at 0 Hz or at the Nyquist frequency.
101 If margin > 0 these frequencies are filtered out and the number of formants can never exceed
102 my maxnCoefficients / 2.
103 */
104 const integer maximumNumberOfFormants = ( margin == 0.0 ? my maxnCoefficients : (my maxnCoefficients + 1) / 2 );
105 const integer maximumNumberOfPolynomialCoefficients = my maxnCoefficients + 1;
106 integer numberOfSuspectFrames = 0;
107 const integer interval = ( my maxnCoefficients > 20 ? 1 : 10 );
108 Melder_require (my maxnCoefficients < 100,
109 U"We cannot find the roots of a polynomial of order > 99.");
110 Melder_require (margin < samplingFrequency / 4.0,
111 U"Margin should be smaller than ", samplingFrequency / 4.0, U".");
112 autoFormant thee = Formant_create (my xmin, my xmax, my nx, my dx, my x1, maximumNumberOfFormants);
113 autoPolynomial polynomial = Polynomial_create (-1.0, 1.0, my maxnCoefficients);
114 autoRoots roots = Roots_create (my maxnCoefficients);
115 autoVEC workspace = raw_VEC (maximumNumberOfPolynomialCoefficients * (maximumNumberOfPolynomialCoefficients + 9));
116 autoMelderProgress progress (U"LPC to Formant");
117 for (integer iframe = 1; iframe <= my nx; iframe ++) {
118 const LPC_Frame lpcFrame = & my d_frames [iframe];
119 const Formant_Frame formantFrame = & thy frames [iframe];
120 Formant_Frame_init (formantFrame, maximumNumberOfFormants);
121 try {
122 LPC_Frame_into_Formant_Frame_mt (lpcFrame, formantFrame, my samplingPeriod, margin, polynomial.get(), roots.get(), workspace.get());
123 //LPC_Frame_into_Formant_Frame (lpcFrame, formant, my samplingPeriod, margin);
124 } catch (MelderError) {
125 Melder_clearError();
126 numberOfSuspectFrames ++;
127 }
128 if (interval == 1 || iframe % interval == 1)
129 Melder_progress ((double) iframe / my nx, U"LPC to Formant: frame ", iframe, U" out of ", my nx, U".");
130 }
131 Formant_sort (thee.get());
132 if (numberOfSuspectFrames > 0)
133 Melder_warning (numberOfSuspectFrames, U" formant frames out of ", my nx, U" are suspect.");
134 return thee;
135 } catch (MelderError) {
136 Melder_throw (me, U": no Formant created.");
137 }
138 }
139
LPC_to_Formant(LPC me,double margin)140 autoFormant LPC_to_Formant (LPC me, double margin) {
141 try {
142 const integer numberOfProcessors = std::thread::hardware_concurrency ();
143 if (numberOfProcessors <= 1) {
144 /*
145 We cannot use multithreading.
146 */
147 return LPC_to_Formant_noThreads (me, margin);
148 }
149 const double samplingFrequency = 1.0 / my samplingPeriod;
150 Melder_require (my maxnCoefficients < 100,
151 U"We cannot find the roots of a polynomial of order > 99.");
152 Melder_require (margin < samplingFrequency / 4.0,
153 U"Margin should be smaller than ", samplingFrequency / 4.0, U".");
154 /*
155 In a very extreme case all roots of the lpc polynomial might be real.
156 A real root gives either a frequency at 0 Hz or at the Nyquist frequency.
157 To play it safely, the maximum number of formants is than equal to the number of coefficients.
158 However, to get the acoustically "real" formant frequencies we always neglect frequencies at 0 Hz
159 and at the Nyquist.
160 If margin > 0 these frequencies are filtered out and the number of formants can never exceed
161 int ( my maxnCoefficients / 2 ) because if maxnCoefficients is uneven than at least one of the roots is real.
162 */
163 const integer maximumNumberOfFormants = ( margin == 0.0 ? my maxnCoefficients : (my maxnCoefficients + 1) / 2 );
164 const integer maximumNumberOfPolynomialCoefficients = my maxnCoefficients + 1;
165 const integer numberOfFrames = my nx;
166 autoFormant thee = Formant_create (my xmin, my xmax, numberOfFrames, my dx, my x1, maximumNumberOfFormants);
167 for (integer iframe = 1; iframe <= numberOfFrames; iframe ++) {
168 const Formant_Frame formantFrame = & thy frames [iframe];
169 Formant_Frame_init (formantFrame, maximumNumberOfFormants);
170 }
171
172 constexpr integer maximumNumberOfThreads = 16;
173 integer numberOfThreads, numberOfFramesPerThread = 25;
174 NUMgetThreadingInfo (numberOfFrames, std::min (numberOfProcessors, maximumNumberOfThreads), & numberOfFramesPerThread, & numberOfThreads);
175 /*
176 Reserve working memory for each thread
177 */
178 autoPolynomial polynomials [maximumNumberOfThreads + 1];
179 autoRoots roots [maximumNumberOfThreads + 1];
180 for (integer ithread = 1; ithread <= numberOfThreads; ithread ++) {
181 polynomials [ithread] = Polynomial_create (-1.0, 1.0, my maxnCoefficients);
182 roots [ithread] = Roots_create (my maxnCoefficients);
183 }
184 autoMAT workspaces = raw_MAT (numberOfThreads, maximumNumberOfPolynomialCoefficients * (maximumNumberOfPolynomialCoefficients + 9));
185 std::vector <std::thread> thread (numberOfThreads);
186 std::atomic<integer> numberOfSuspectFrames (0);
187
188 try {
189 for (integer ithread = 1; ithread <= numberOfThreads; ithread ++) {
190 Polynomial p = polynomials [ithread].get();
191 Roots r = roots [ithread].get();
192 Formant formant = thee.get();
193 LPC lpc = me;
194 VEC workspace = workspaces.row (ithread);
195 const integer firstFrame = 1 + (ithread - 1) * numberOfFramesPerThread;
196 const integer lastFrame = ( ithread == numberOfThreads ? numberOfFrames : firstFrame + numberOfFramesPerThread - 1 );
197
198 thread [ithread - 1] = std::thread ([=, & numberOfSuspectFrames] () {
199 for (integer iframe = firstFrame; iframe <= lastFrame; iframe ++) {
200 const LPC_Frame lpcFrame = & lpc -> d_frames [iframe];
201 const Formant_Frame formantFrame = & formant -> frames [iframe];
202 try {
203 LPC_Frame_into_Formant_Frame_mt (lpcFrame, formantFrame, my samplingPeriod, margin, p, r, workspace);
204 } catch (MelderError) {
205 numberOfSuspectFrames ++;
206 }
207 }
208 });
209 }
210 } catch (MelderError) {
211 for (integer ithread = 1; ithread <= numberOfThreads; ithread ++) {
212 if (thread [ithread - 1]. joinable ())
213 thread [ithread - 1]. join ();
214 }
215 throw;
216 }
217 for (integer ithread = 1; ithread <= numberOfThreads; ithread ++)
218 thread [ithread - 1]. join ();
219
220
221 Formant_sort (thee.get());
222 if (numberOfSuspectFrames > 0)
223 Melder_warning ((integer) numberOfSuspectFrames, U" formant frames out of ", numberOfFrames, U" are suspect.");
224 return thee;
225 } catch (MelderError) {
226 Melder_throw (me, U": no Formant created.");
227 }
228 }
229
Formant_Frame_into_LPC_Frame(Formant_Frame me,LPC_Frame thee,double samplingPeriod)230 void Formant_Frame_into_LPC_Frame (Formant_Frame me, LPC_Frame thee, double samplingPeriod) {
231 if (my numberOfFormants < 1)
232 return;
233 const double nyquistFrequency = 0.5 / samplingPeriod;
234 integer numberOfPoles = 2 * my numberOfFormants;
235 autoVEC lpc = zero_VEC (numberOfPoles + 2); // all odd coefficients have to be initialized to zero
236 lpc [2] = 1.0;
237 integer m = 2;
238 for (integer iformant = 1; iformant <= my numberOfFormants; iformant ++) {
239 const double formantFrequency = my formant [iformant]. frequency;
240 if (formantFrequency > nyquistFrequency)
241 continue;
242 /*
243 D(z): 1 + p z^-1 + q z^-2
244 */
245 const double r = exp (- NUMpi * my formant [iformant]. bandwidth * samplingPeriod);
246 const double p = - 2.0 * r * cos (NUM2pi * formantFrequency * samplingPeriod);
247 const double q = r * r;
248 /*
249 By setting the two extra elements (0, 1) in the lpc vector we can avoid boundary testing;
250 lpc [3..n+2] come to contain the coefficients.
251 */
252 for (integer j = m + 2; j > 2; j --)
253 lpc [j] += p * lpc [j - 1] + q * lpc [j - 2];
254 m += 2;
255 }
256 if (thy nCoefficients < numberOfPoles)
257 numberOfPoles = thy nCoefficients;
258 for (integer i = 1; i <= numberOfPoles; i ++)
259 thy a [i] = lpc [i + 2];
260 thy gain = my intensity;
261 }
262
Formant_to_LPC(Formant me,double samplingPeriod)263 autoLPC Formant_to_LPC (Formant me, double samplingPeriod) {
264 try {
265 autoLPC thee = LPC_create (my xmin, my xmax, my nx, my dx, my x1, 2 * my maxnFormants, samplingPeriod);
266
267 for (integer i = 1; i <= my nx; i ++) {
268 const Formant_Frame f = & my frames [i];
269 const LPC_Frame lpc = & thy d_frames [i];
270 const integer numberOfCoefficients = 2 * f -> numberOfFormants;
271 LPC_Frame_init (lpc, numberOfCoefficients);
272 Formant_Frame_into_LPC_Frame (f, lpc, samplingPeriod);
273 }
274 return thee;
275 } catch (MelderError) {
276 Melder_throw (me, U": no LPC created.");
277 }
278 }
279
280 /* End of file LPC_and_Formant.cpp */
281