1 /* Ltas.cpp
2 *
3 * Copyright (C) 1992-2012,2015,2016,2017 Paul Boersma
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.
13 * See the GNU 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 * a selection of changes:
21 * pb 2005/11/26 pitch-corrected Ltas
22 */
23
24 #include "Ltas.h"
25 #include "Sound_and_Spectrum.h"
26 #include "Sound_to_PointProcess.h"
27
28 Thing_implement (Ltas, Vector, 2);
29
v_info()30 void structLtas :: v_info () {
31 double meanPowerDensity;
32 structDaata :: v_info ();
33 MelderInfo_writeLine (U"Frequency domain:");
34 MelderInfo_writeLine (U" Lowest frequency: ", xmin, U" Hz");
35 MelderInfo_writeLine (U" Highest frequency: ", xmax, U" Hz");
36 MelderInfo_writeLine (U" Total frequency domain: ", xmax - xmin, U" Hz");
37 MelderInfo_writeLine (U"Frequency sampling:");
38 MelderInfo_writeLine (U" Number of frequency bands: ", nx);
39 MelderInfo_writeLine (U" Width of each band: ", dx, U" Hz");
40 MelderInfo_writeLine (U" First band centred at: ", x1, U" Hz");
41 meanPowerDensity = Sampled_getMean (this, xmin, xmax, 0, 1, false);
42 MelderInfo_writeLine (U"Total SPL: ", Melder_single (10.0 * log10 (meanPowerDensity * (xmax - xmin))), U" dB");
43 }
44
v_convertStandardToSpecialUnit(double value,integer,int unit)45 double structLtas :: v_convertStandardToSpecialUnit (double value, integer /* level */, int unit) {
46 if (unit == 1) {
47 return pow (10.0, 0.1 * value); // energy
48 } else if (unit == 2) {
49 return pow (2.0, 0.1 * value); // sones
50 }
51 return value;
52 }
53
v_convertSpecialToStandardUnit(double value,integer,int unit)54 double structLtas :: v_convertSpecialToStandardUnit (double value, integer /* level */, int unit) {
55 return
56 unit == 1 ?
57 10.0 * log10 (value) : // value = energy
58 unit == 2 ?
59 10.0 * NUMlog2 (value) : // value = sones
60 value; // value = dB
61 }
62
Ltas_create(integer nx,double dx)63 autoLtas Ltas_create (integer nx, double dx) {
64 try {
65 autoLtas me = Thing_new (Ltas);
66 Matrix_init (me.get(), 0.0, nx * dx, nx, dx, 0.5 * dx, 1.0, 1.0, 1, 1.0, 1.0);
67 return me;
68 } catch (MelderError) {
69 Melder_throw (U"Ltas not created.");
70 }
71 }
72
Ltas_draw(Ltas me,Graphics g,double fmin,double fmax,double minimum,double maximum,bool garnish,conststring32 method)73 void Ltas_draw (Ltas me, Graphics g, double fmin, double fmax, double minimum, double maximum, bool garnish, conststring32 method) {
74 Vector_draw (me, g, & fmin, & fmax, & minimum, & maximum, 1.0, method);
75 if (garnish) {
76 Graphics_drawInnerBox (g);
77 Graphics_textBottom (g, true, U"Frequency (Hz)");
78 Graphics_marksBottom (g, 2, true, true, false);
79 Graphics_textLeft (g, true, U"Sound pressure level (dB/Hz)");
80 Graphics_marksLeft (g, 2, true, true, false);
81 }
82 }
83
Ltas_getSlope(Ltas me,double f1min,double f1max,double f2min,double f2max,int averagingUnits)84 double Ltas_getSlope (Ltas me, double f1min, double f1max, double f2min, double f2max, int averagingUnits) {
85 double low = Sampled_getMean (me, f1min, f1max, 0, averagingUnits, false);
86 double high = Sampled_getMean (me, f2min, f2max, 0, averagingUnits, false);
87 if (isundef (low) || isundef (high)) return undefined;
88 return averagingUnits == 3 ? high - low : Function_convertSpecialToStandardUnit (me, high / low, 0, averagingUnits);
89 }
90
Ltas_getLocalPeakHeight(Ltas me,double environmentMin,double environmentMax,double peakMin,double peakMax,int averagingUnits)91 double Ltas_getLocalPeakHeight (Ltas me, double environmentMin, double environmentMax, double peakMin, double peakMax, int averagingUnits) {
92 double environmentLow = Sampled_getMean (me, environmentMin, peakMin, 0, averagingUnits, false);
93 double environmentHigh = Sampled_getMean (me, peakMax, environmentMax, 0, averagingUnits, false);
94 double peak = Sampled_getMean (me, peakMin, peakMax, 0, averagingUnits, false);
95 if (isundef (environmentLow) || isundef (environmentHigh) || isundef (peak)) return undefined;
96 return averagingUnits == 3 ? peak - 0.5 * (environmentLow + environmentHigh) :
97 Function_convertSpecialToStandardUnit (me, peak / (0.5 * (environmentLow + environmentHigh)), 0, averagingUnits);
98 }
99
Ltas_to_Matrix(Ltas me)100 autoMatrix Ltas_to_Matrix (Ltas me) {
101 try {
102 autoMatrix thee = Thing_new (Matrix);
103 my structMatrix :: v_copy (thee.get());
104 return thee;
105 } catch (MelderError) {
106 Melder_throw (me, U": not converted to Matrix.");
107 }
108 }
109
Matrix_to_Ltas(Matrix me)110 autoLtas Matrix_to_Ltas (Matrix me) {
111 try {
112 autoLtas thee = Thing_new (Ltas);
113 my structMatrix :: v_copy (thee.get()); // because copying to a descendant of Matrix with additional members should not crash
114 return thee;
115 } catch (MelderError) {
116 Melder_throw (me, U": not converted to Ltas.");
117 }
118 }
119
Ltases_merge(LtasBag ltases)120 autoLtas Ltases_merge (LtasBag ltases) {
121 try {
122 if (ltases->size < 1)
123 Melder_throw (U"Cannot merge zero Ltas objects.");
124 Ltas me = ltases->at [1];
125 autoLtas thee = Data_copy (me);
126 /*
127 * Convert to energy.
128 */
129 for (integer iband = 1; iband <= thy nx; iband ++) {
130 thy z [1] [iband] = pow (10.0, thy z [1] [iband] / 10.0);
131 }
132 for (integer ispec = 2; ispec <= ltases->size; ispec ++) {
133 Ltas him = ltases->at [ispec];
134 if (his xmin != thy xmin || his xmax != thy xmax)
135 Melder_throw (U"Frequency domains do not match.");
136 if (his dx != thy dx)
137 Melder_throw (U"Bandwidths do not match.");
138 if (his nx != thy nx || his x1 != thy x1)
139 Melder_throw (U"Frequency bands do not match.");
140 /*
141 * Add band energies.
142 */
143 for (integer iband = 1; iband <= thy nx; iband ++) {
144 thy z [1] [iband] += pow (10.0, his z [1] [iband] / 10.0);
145 }
146 }
147 /*
148 * Convert back to dB.
149 */
150 for (integer iband = 1; iband <= thy nx; iband ++) {
151 thy z [1] [iband] = 10.0 * log10 (thy z [1] [iband]);
152 }
153 return thee;
154 } catch (MelderError) {
155 Melder_throw (U"Ltas objects not merged.");
156 }
157 }
158
159 Thing_implement (LtasBag, Collection, 0);
160
Ltases_average(LtasBag ltases)161 autoLtas Ltases_average (LtasBag ltases) {
162 try {
163 double factor = -10.0 * log10 (ltases->size);
164 autoLtas thee = Ltases_merge (ltases);
165 for (integer iband = 1; iband <= thy nx; iband ++) {
166 thy z [1] [iband] += factor;
167 }
168 return thee;
169 } catch (MelderError) {
170 Melder_throw (U"Ltas objects not averaged.");
171 }
172 }
173
Ltas_computeTrendLine(Ltas me,double fmin,double fmax)174 autoLtas Ltas_computeTrendLine (Ltas me, double fmin, double fmax) {
175 try {
176 /*
177 * Find the first and last bin.
178 */
179 integer imin, imax, n;
180 if ((n = Sampled_getWindowSamples (me, fmin, fmax, & imin, & imax)) < 2)
181 Melder_throw (U"Number of bins too low (", n, U"). Should be at least 2.");
182 autoLtas thee = Data_copy (me);
183 /*
184 * Compute average amplitude and frequency.
185 */
186 longdouble sum = 0.0, numerator = 0.0, denominator = 0.0;
187 for (integer i = imin; i <= imax; i ++) {
188 sum += thy z [1] [i];
189 }
190 double amean = double (sum / n);
191 double fmean = thy x1 + (0.5 * (imin + imax) - 1) * thy dx;
192 /*
193 * Compute slope.
194 */
195 for (integer i = imin; i <= imax; i ++) {
196 double da = thy z [1] [i] - amean, df = thy x1 + (i - 1) * thy dx - fmean;
197 numerator += da * df;
198 denominator += df * df;
199 }
200 double slope = double (numerator / denominator);
201 /*
202 * Modify bins.
203 */
204 for (integer i = 1; i <= thy nx; i ++) {
205 double df = thy x1 + (i - 1) * thy dx - fmean;
206 thy z [1] [i] = amean + slope * df;
207 }
208 return thee;
209 } catch (MelderError) {
210 Melder_throw (me, U": trend line not computed.");
211 }
212 }
213
Ltas_subtractTrendLine(Ltas me,double fmin,double fmax)214 autoLtas Ltas_subtractTrendLine (Ltas me, double fmin, double fmax) {
215 try {
216 /*
217 * Find the first and last bin.
218 */
219 integer imin, imax, n;
220 if ((n = Sampled_getWindowSamples (me, fmin, fmax, & imin, & imax)) < 2)
221 Melder_throw (U"Number of bins too low (", n, U"). Should be at least 2.");
222 autoLtas thee = Data_copy (me);
223 /*
224 * Compute average amplitude and frequency.
225 */
226 longdouble sum = 0.0;
227 for (integer i = imin; i <= imax; i ++) {
228 sum += thy z [1] [i];
229 }
230 double amean = (double) sum / n;
231 double fmean = thy x1 + (0.5 * (imin + imax) - 1) * thy dx;
232 /*
233 * Compute slope.
234 */
235 longdouble numerator = 0.0, denominator = 0.0;
236 for (integer i = imin; i <= imax; i ++) {
237 double da = thy z [1] [i] - amean, df = thy x1 + (i - 1) * thy dx - fmean;
238 numerator += da * df;
239 denominator += df * df;
240 }
241 double slope = (double) (numerator / denominator);
242 /*
243 * Modify bins.
244 */
245 for (integer i = 1; i < imin; i ++) {
246 thy z [1] [i] = 0.0;
247 }
248 for (integer i = imin; i <= imax; i ++) {
249 double df = thy x1 + (i - 1) * thy dx - fmean;
250 thy z [1] [i] -= amean + slope * df;
251 }
252 for (integer i = imax + 1; i <= thy nx; i ++) {
253 thy z [1] [i] = 0.0;
254 }
255 return thee;
256 } catch (MelderError) {
257 Melder_throw (me, U": trend line not subtracted.");
258 }
259 }
260
Spectrum_to_Ltas(Spectrum me,double bandWidth)261 autoLtas Spectrum_to_Ltas (Spectrum me, double bandWidth) {
262 try {
263 integer numberOfBands = Melder_iceiling ((my xmax - my xmin) / bandWidth);
264 if (bandWidth <= my dx)
265 Melder_throw (U"Bandwidth must be greater than ", my dx, U".");
266 autoLtas thee = Thing_new (Ltas);
267 Matrix_init (thee.get(), my xmin, my xmax, numberOfBands, bandWidth, my xmin + 0.5 * bandWidth, 1.0, 1.0, 1, 1.0, 1.0);
268 for (integer iband = 1; iband <= numberOfBands; iband ++) {
269 double fmin = thy xmin + (iband - 1) * bandWidth;
270 double meanEnergyDensity = Sampled_getMean (me, fmin, fmin + bandWidth, 0, 1, false);
271 double meanPowerDensity = meanEnergyDensity * my dx; // as an approximation for a division by the original duration
272 thy z [1] [iband] = meanPowerDensity == 0.0 ? -300.0 : 10.0 * log10 (meanPowerDensity / 4.0e-10);
273 }
274 return thee;
275 } catch (MelderError) {
276 Melder_throw (me, U": not converted to Ltas.");
277 }
278 }
279
Spectrum_to_Ltas_1to1(Spectrum me)280 autoLtas Spectrum_to_Ltas_1to1 (Spectrum me) {
281 try {
282 autoLtas thee = Thing_new (Ltas);
283 Matrix_init (thee.get(), my xmin, my xmax, my nx, my dx, my x1, 1.0, 1.0, 1, 1.0, 1.0);
284 for (integer iband = 1; iband <= my nx; iband ++) {
285 thy z [1] [iband] = Sampled_getValueAtSample (me, iband, 0, 2);
286 }
287 return thee;
288 } catch (MelderError) {
289 Melder_throw (me, U": not converted to Ltas.");
290 }
291 }
292
Sound_to_Ltas(Sound me,double bandwidth)293 autoLtas Sound_to_Ltas (Sound me, double bandwidth) {
294 try {
295 autoSpectrum thee = Sound_to_Spectrum (me, true);
296 autoLtas him = Spectrum_to_Ltas (thee.get(), bandwidth);
297 double correction = -10.0 * log10 (thy dx * my nx * my dx);
298 for (integer iband = 1; iband <= his nx; iband ++) {
299 his z [1] [iband] += correction;
300 }
301 return him;
302 } catch (MelderError) {
303 Melder_throw (me, U": LTAS analysis not performed.");
304 }
305 }
306
PointProcess_Sound_to_Ltas(PointProcess pulses,Sound sound,double maximumFrequency,double bandWidth,double shortestPeriod,double longestPeriod,double maximumPeriodFactor)307 autoLtas PointProcess_Sound_to_Ltas (PointProcess pulses, Sound sound,
308 double maximumFrequency, double bandWidth,
309 double shortestPeriod, double longestPeriod, double maximumPeriodFactor)
310 {
311 try {
312 integer numberOfPeriods = pulses -> nt - 2, totalNumberOfEnergies = 0;
313 autoLtas ltas = Ltas_create (Melder_ifloor (maximumFrequency / bandWidth), bandWidth);
314 ltas -> xmax = maximumFrequency;
315 autoLtas numbers = Data_copy (ltas.get());
316 if (numberOfPeriods < 1)
317 Melder_throw (U"Cannot compute an Ltas if there are no periods in the point process.");
318 autoMelderProgress progress (U"Ltas analysis...");
319 for (integer ipulse = 2; ipulse < pulses -> nt; ipulse ++) {
320 double leftInterval = pulses -> t [ipulse] - pulses -> t [ipulse - 1];
321 double rightInterval = pulses -> t [ipulse + 1] - pulses -> t [ipulse];
322 double intervalFactor = leftInterval > rightInterval ? leftInterval / rightInterval : rightInterval / leftInterval;
323 Melder_progress ((double) ipulse / pulses -> nt, U"Sound & PointProcess: To Ltas: pulse ", ipulse, U" out of ", pulses -> nt);
324 if (leftInterval >= shortestPeriod && leftInterval <= longestPeriod &&
325 rightInterval >= shortestPeriod && rightInterval <= longestPeriod &&
326 intervalFactor <= maximumPeriodFactor)
327 {
328 /*
329 * We have a period! Compute the spectrum.
330 */
331 autoSound period = Sound_extractPart (sound,
332 pulses -> t [ipulse] - 0.5 * leftInterval, pulses -> t [ipulse] + 0.5 * rightInterval,
333 kSound_windowShape::RECTANGULAR, 1.0, false);
334 autoSpectrum spectrum = Sound_to_Spectrum (period.get(), false);
335 for (integer ifreq = 1; ifreq <= spectrum -> nx; ifreq ++) {
336 double frequency = spectrum -> xmin + (ifreq - 1) * spectrum -> dx;
337 double realPart = spectrum -> z [1] [ifreq];
338 double imaginaryPart = spectrum -> z [2] [ifreq];
339 double energy = (realPart * realPart + imaginaryPart * imaginaryPart) * 2.0 * spectrum -> dx /* OLD: * sound -> nx */;
340 integer iband = Melder_iceiling (frequency / bandWidth);
341 if (iband >= 1 && iband <= ltas -> nx) {
342 ltas -> z [1] [iband] += energy;
343 numbers -> z [1] [iband] += 1;
344 totalNumberOfEnergies += 1;
345 }
346 }
347 } else {
348 numberOfPeriods -= 1;
349 }
350 }
351 if (numberOfPeriods < 1)
352 Melder_throw (U"There are no periods in the point process.");
353 for (integer iband = 1; iband <= ltas -> nx; iband ++) {
354 if (numbers -> z [1] [iband] == 0.0) {
355 ltas -> z [1] [iband] = undefined;
356 } else {
357 /*
358 * Each bin now contains a total energy in Pa2 sec.
359 * To convert this to power density, we
360 */
361 double totalEnergyInThisBand = ltas -> z [1] [iband];
362 if (false /* i.e. if you just want to have a spectrum of the voiced parts... */) {
363 double energyDensityInThisBand = totalEnergyInThisBand / ltas -> dx;
364 double powerDensityInThisBand = energyDensityInThisBand / (sound -> xmax - sound -> xmin);
365 ltas -> z [1] [iband] = 10.0 * log10 (powerDensityInThisBand / 4.0e-10);
366 } else {
367 /*
368 * And this is what we really want. The total energy has to be redistributed.
369 */
370 double meanEnergyInThisBand = totalEnergyInThisBand / numbers -> z [1] [iband];
371 double meanNumberOfEnergiesPerBand = (double) totalNumberOfEnergies / ltas -> nx;
372 double redistributedEnergyInThisBand = meanEnergyInThisBand * meanNumberOfEnergiesPerBand;
373 double redistributedEnergyDensityInThisBand = redistributedEnergyInThisBand / ltas -> dx;
374 double redistributedPowerDensityInThisBand = redistributedEnergyDensityInThisBand / (sound -> xmax - sound -> xmin);
375 ltas -> z [1] [iband] = 10.0 * log10 (redistributedPowerDensityInThisBand / 4.0e-10);
376 /* OLD: ltas -> z [1] [iband] = 10.0 * log10 (ltas -> z [1] [iband] / numbers -> z [1] [iband] * sound -> nx);*/
377 }
378 }
379 }
380 for (integer iband = 1; iband <= ltas -> nx; iband ++) {
381 if (isundef (ltas -> z [1] [iband])) {
382 integer ibandleft = iband - 1, ibandright = iband + 1;
383 while (ibandleft >= 1 && isundef (ltas -> z [1] [ibandleft])) ibandleft --;
384 while (ibandright <= ltas -> nx && isundef (ltas -> z [1] [ibandright])) ibandright ++;
385 if (ibandleft < 1 && ibandright > ltas -> nx)
386 Melder_throw (U"Cannot create an Ltas without energy in any bins.");
387 if (ibandleft < 1) {
388 ltas -> z [1] [iband] = ltas -> z [1] [ibandright];
389 } else if (ibandright > ltas -> nx) {
390 ltas -> z [1] [iband] = ltas -> z [1] [ibandleft];
391 } else {
392 double frequency = ltas -> x1 + (iband - 1) * ltas -> dx;
393 double fleft = ltas -> x1 + (ibandleft - 1) * ltas -> dx;
394 double fright = ltas -> x1 + (ibandright - 1) * ltas -> dx;
395 ltas -> z [1] [iband] = ((fright - frequency) * ltas -> z [1] [ibandleft]
396 + (frequency - fleft) * ltas -> z [1] [ibandright]) / (fright - fleft);
397 }
398 }
399 }
400 return ltas;
401 } catch (MelderError) {
402 Melder_throw (sound, U" & ", pulses, U": LTAS analysis not performed.");
403 }
404 }
405
Sound_to_Ltas_pitchCorrected(Sound sound,double minimumPitch,double maximumPitch,double maximumFrequency,double bandWidth,double shortestPeriod,double longestPeriod,double maximumPeriodFactor)406 autoLtas Sound_to_Ltas_pitchCorrected (Sound sound, double minimumPitch, double maximumPitch,
407 double maximumFrequency, double bandWidth,
408 double shortestPeriod, double longestPeriod, double maximumPeriodFactor)
409 {
410 try {
411 autoPointProcess pulses = Sound_to_PointProcess_periodic_cc (sound, minimumPitch, maximumPitch);
412 autoLtas ltas = PointProcess_Sound_to_Ltas (pulses.get(), sound, maximumFrequency, bandWidth,
413 shortestPeriod, longestPeriod, maximumPeriodFactor);
414 return ltas;
415 } catch (MelderError) {
416 Melder_throw (sound, U": pitch-corrected LTAS analysis not performed.");
417 }
418 }
419
PointProcess_Sound_to_Ltas_harmonics(PointProcess pulses,Sound sound,integer maximumHarmonic,double shortestPeriod,double longestPeriod,double maximumPeriodFactor)420 autoLtas PointProcess_Sound_to_Ltas_harmonics (PointProcess pulses, Sound sound,
421 integer maximumHarmonic,
422 double shortestPeriod, double longestPeriod, double maximumPeriodFactor)
423 {
424 try {
425 integer numberOfPeriods = pulses -> nt - 2;
426 autoLtas ltas = Ltas_create (maximumHarmonic, 1.0);
427 ltas -> xmax = maximumHarmonic;
428 if (numberOfPeriods < 1)
429 Melder_throw (U"There are no periods in the point process.");
430 autoMelderProgress progress (U"LTAS (harmonics) analysis...");
431 for (integer ipulse = 2; ipulse < pulses -> nt; ipulse ++) {
432 double leftInterval = pulses -> t [ipulse] - pulses -> t [ipulse - 1];
433 double rightInterval = pulses -> t [ipulse + 1] - pulses -> t [ipulse];
434 double intervalFactor = leftInterval > rightInterval ? leftInterval / rightInterval : rightInterval / leftInterval;
435 Melder_progress ((double) ipulse / pulses -> nt, U"Sound & PointProcess: To Ltas: pulse ", ipulse, U" out of ", pulses -> nt);
436 if (leftInterval >= shortestPeriod && leftInterval <= longestPeriod &&
437 rightInterval >= shortestPeriod && rightInterval <= longestPeriod &&
438 intervalFactor <= maximumPeriodFactor)
439 {
440 /*
441 * We have a period! Compute the spectrum.
442 */
443 integer localMaximumHarmonic;
444 autoSound period = Sound_extractPart (sound,
445 pulses -> t [ipulse] - 0.5 * leftInterval, pulses -> t [ipulse] + 0.5 * rightInterval,
446 kSound_windowShape::RECTANGULAR, 1.0, false);
447 autoSpectrum spectrum = Sound_to_Spectrum (period.get(), false);
448 localMaximumHarmonic = maximumHarmonic < spectrum -> nx ? maximumHarmonic : spectrum -> nx;
449 for (integer iharm = 1; iharm <= localMaximumHarmonic; iharm ++) {
450 double realPart = spectrum -> z [1] [iharm];
451 double imaginaryPart = spectrum -> z [2] [iharm];
452 double energy = (realPart * realPart + imaginaryPart * imaginaryPart) * 2.0 * spectrum -> dx;
453 ltas -> z [1] [iharm] += energy;
454 }
455 } else {
456 numberOfPeriods -= 1;
457 }
458 }
459 if (numberOfPeriods < 1)
460 Melder_throw (U"There are no periods in the point process.");
461 for (integer iharm = 1; iharm <= ltas -> nx; iharm ++) {
462 if (ltas -> z [1] [iharm] == 0.0) {
463 ltas -> z [1] [iharm] = -300.0;
464 } else {
465 double energyInThisBand = ltas -> z [1] [iharm];
466 double powerInThisBand = energyInThisBand / (sound -> xmax - sound -> xmin);
467 ltas -> z [1] [iharm] = 10.0 * log10 (powerInThisBand / 4.0e-10);
468 }
469 }
470 return ltas;
471 } catch (MelderError) {
472 Melder_throw (sound, U" & ", pulses, U": LTAS analysis (harmonics) not performed.");
473 }
474 }
475
476 /* End of file Ltas.cpp */
477