1 /* Spectrum.cpp
2 *
3 * Copyright (C) 1992-2008,2011,2012,2014-2021 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 2001/11/30 Spectral moments
22 * pb 2002/05/22 changed sign of imaginary part
23 * pb 2006/02/06 better cepstral smoothing
24 */
25
26 #include "Sound_and_Spectrum.h"
27 #include "SpectrumTier.h"
28
29 #include "oo_DESTROY.h"
30 #include "Spectrum_def.h"
31 #include "oo_COPY.h"
32 #include "Spectrum_def.h"
33 #include "oo_EQUAL.h"
34 #include "Spectrum_def.h"
35 #include "oo_CAN_WRITE_AS_ENCODING.h"
36 #include "Spectrum_def.h"
37 #include "oo_WRITE_TEXT.h"
38 #include "Spectrum_def.h"
39 #include "oo_READ_TEXT.h"
40 #include "Spectrum_def.h"
41 #include "oo_WRITE_BINARY.h"
42 #include "Spectrum_def.h"
43 #include "oo_READ_BINARY.h"
44 #include "Spectrum_def.h"
45 #include "oo_DESCRIPTION.h"
46 #include "Spectrum_def.h"
47
48 Thing_implement (Spectrum, Matrix, 2);
49
v_info()50 void structSpectrum :: v_info () {
51 structDaata :: v_info ();
52 MelderInfo_writeLine (U"Frequency domain:");
53 MelderInfo_writeLine (U" Lowest frequency: ", xmin, U" Hz");
54 MelderInfo_writeLine (U" Highest frequency: ", xmax, U" Hz");
55 MelderInfo_writeLine (U" Total bandwidth: ", xmax - xmin, U" Hz");
56 MelderInfo_writeLine (U"Frequency sampling:");
57 MelderInfo_writeLine (U" Number of frequency bands (bins): ", nx);
58 MelderInfo_writeLine (U" Frequency step (bin width): ", dx, U" Hz");
59 MelderInfo_writeLine (U" First frequency band around (bin centre at): ", x1, U" Hz");
60 MelderInfo_writeLine (U"Total energy: ", Melder_single (Spectrum_getBandEnergy (this, 0.0, 0.0)), U" Pa2 sec");
61 }
62
v_getValueAtSample(integer isamp,integer which,int units)63 double structSpectrum :: v_getValueAtSample (integer isamp, integer which, int units) {
64 if (units == 0) {
65 return which == 1 ? z [1] [isamp] : which == 2 ? z [2] [isamp] : undefined;
66 } else {
67 /*
68 The energy in a bin is 2 * (re^2 + im^2) times the bin width.
69 The factor of 2 derives from the assumption that the spectrum contains positive-frequency values only,
70 and that the negative-frequency values have the same norm, since they are the complex conjugates
71 of the positive-frequency values.
72 The bin width is usually my dx, although it is 0.5 * my dx for the first bin,
73 and if the last bin is centred on my xmax [the Nyquist] then it is also 0.5 * my dx for the last bin.
74 */
75 const double energyDensity = 2.0 * (sqr (z [1] [isamp]) + sqr (z [2] [isamp]));
76 /* Pa2/Hz2; sum of positive and negative frequencies */
77 if (units == 1) {
78 return energyDensity;
79 } else {
80 const double powerDensity = energyDensity * our dx; // Pa^2 Hz-2 s-1, after division by approximate duration
81 if (units == 2)
82 /* "dB/Hz" */
83 return powerDensity == 0.0 ? -300.0 : 10.0 * log10 (powerDensity / 4.0e-10);
84 }
85 }
86 return undefined;
87 }
88
Spectrum_create(double fmax,integer nf)89 autoSpectrum Spectrum_create (double fmax, integer nf) {
90 try {
91 autoSpectrum me = Thing_new (Spectrum);
92 Matrix_init (me.get(), 0.0, fmax, nf, fmax / (nf - 1), 0.0, 1.0, 2.0, 2, 1.0, 1.0);
93 return me;
94 } catch (MelderError) {
95 Melder_throw (U"Spectrum not created.");
96 }
97 }
98
Spectrum_getPowerDensityRange(Spectrum me,double * minimum,double * maximum)99 int Spectrum_getPowerDensityRange (Spectrum me, double *minimum, double *maximum) {
100 *minimum = 1e308;
101 *maximum = 0.0;
102 for (integer ifreq = 1; ifreq <= my nx; ifreq ++) {
103 const double oneSidedPowerSpectralDensity = // Pa2 Hz-2 s-1
104 2.0 * (sqr (my z [1] [ifreq]) + sqr (my z [2] [ifreq])) * my dx;
105 if (oneSidedPowerSpectralDensity < *minimum)
106 *minimum = oneSidedPowerSpectralDensity;
107 if (oneSidedPowerSpectralDensity > *maximum)
108 *maximum = oneSidedPowerSpectralDensity;
109 }
110 if (*maximum == 0.0)
111 return 0;
112 *minimum = 10.0 * log10 (*minimum / 4.0e-10);
113 *maximum = 10.0 * log10 (*maximum / 4.0e-10);
114 return 1;
115 }
116
Spectrum_drawInside(Spectrum me,Graphics g,double fmin,double fmax,double minimum,double maximum)117 void Spectrum_drawInside (Spectrum me, Graphics g, double fmin, double fmax, double minimum, double maximum) {
118 const bool autoscaling = ( minimum >= maximum );
119 if (fmax <= fmin) {
120 fmin = my xmin;
121 fmax = my xmax;
122 }
123 integer ifmin, ifmax;
124 const integer nf = Matrix_getWindowSamplesX (me, fmin, fmax, & ifmin, & ifmax);
125 if (nf == 0)
126 return;
127 autoVEC ybuffer = zero_VEC (nf);
128 double *yWC = & ybuffer [1 - ifmin];
129
130 /*
131 First pass: compute power density.
132 */
133 if (autoscaling)
134 maximum = -1e308;
135 for (integer ifreq = ifmin; ifreq <= ifmax; ifreq ++) {
136 const double y = my v_getValueAtSample (ifreq, 0, 2);
137 if (autoscaling && y > maximum)
138 maximum = y;
139 yWC [ifreq] = y;
140 }
141 if (autoscaling) {
142 constexpr double defaultDynamicRange_dB = 60.0;
143 minimum = maximum - defaultDynamicRange_dB;
144 if (minimum == maximum) { // because infinity minus something is still infinity
145 Graphics_setWindow (g, 0.0, 1.0, 0.0, 1.0);
146 Graphics_setTextAlignment (g, kGraphics_horizontalAlignment::CENTRE, Graphics_HALF);
147 Graphics_text (g, 0.5, 0.5, U"(undefined spectrum values cannot be drawn)");
148 return;
149 }
150 }
151
152 /*
153 Second pass: clip.
154 */
155 for (integer ifreq = ifmin; ifreq <= ifmax; ifreq ++)
156 Melder_clip (minimum, & yWC [ifreq], maximum);
157
158 Graphics_setWindow (g, fmin, fmax, minimum, maximum);
159 Graphics_function (g, yWC, ifmin, ifmax, Matrix_columnToX (me, ifmin), Matrix_columnToX (me, ifmax));
160 }
161
Spectrum_draw(Spectrum me,Graphics g,double fmin,double fmax,double minimum,double maximum,bool garnish)162 void Spectrum_draw (Spectrum me, Graphics g, double fmin, double fmax, double minimum, double maximum, bool garnish) {
163 Graphics_setInner (g);
164 Spectrum_drawInside (me, g, fmin, fmax, minimum, maximum);
165 Graphics_unsetInner (g);
166 if (garnish) {
167 Graphics_drawInnerBox (g);
168 Graphics_textBottom (g, true, U"Frequency (Hz)");
169 Graphics_marksBottom (g, 2, true, true, false);
170 Graphics_textLeft (g, true, U"Sound pressure level (dB/Hz)");
171 Graphics_marksLeftEvery (g, 1.0, 20.0, true, true, false);
172 }
173 }
174
Spectrum_drawLogFreq(Spectrum me,Graphics g,double fmin,double fmax,double minimum,double maximum,bool garnish)175 void Spectrum_drawLogFreq (Spectrum me, Graphics g, double fmin, double fmax, double minimum, double maximum, bool garnish) {
176 const bool autoscaling = ( minimum >= maximum );
177 if (fmax <= fmin) {
178 fmin = my xmin;
179 fmax = my xmax;
180 }
181 integer ifmin, ifmax;
182 const integer nf = Matrix_getWindowSamplesX (me, fmin, fmax, & ifmin, & ifmax);
183 if (nf == 0)
184 return;
185 if(ifmin==1)ifmin=2; /* BUG */
186 auto xbuffer = zero_VEC (nf), ybuffer = zero_VEC (nf);
187 double *xWC = & xbuffer [1 - ifmin], *yWC = & ybuffer [1 - ifmin];
188
189 /*
190 First pass: compute power density.
191 */
192 if (autoscaling)
193 maximum = -1e6;
194 for (integer ifreq = ifmin; ifreq <= ifmax; ifreq ++) {
195 xWC [ifreq] = log10 (my x1 + (ifreq - 1) * my dx);
196 yWC [ifreq] = my v_getValueAtSample (ifreq, 0, 2);
197 if (autoscaling && yWC [ifreq] > maximum)
198 maximum = yWC [ifreq];
199 }
200 if (autoscaling)
201 minimum = maximum - 60.0; // default dynamic range is 60 dB
202
203 /*
204 Second pass: clip.
205 */
206 for (integer ifreq = ifmin; ifreq <= ifmax; ifreq ++)
207 Melder_clip (minimum, & yWC [ifreq], maximum);
208
209 Graphics_setInner (g);
210 Graphics_setWindow (g, log10 (fmin), log10 (fmax), minimum, maximum);
211 Graphics_polyline (g, ifmax - ifmin + 1, & xWC [ifmin], & yWC [ifmin]);
212 Graphics_unsetInner (g);
213 if (garnish) {
214 Graphics_drawInnerBox (g);
215 Graphics_textBottom (g, true, U"Frequency (Hz)");
216 Graphics_marksBottomLogarithmic (g, 3, true, true, false);
217 Graphics_textLeft (g, true, U"Sound pressure level (dB/Hz)");
218 Graphics_marksLeftEvery (g, 1.0, 20.0, true, true, false);
219 }
220 }
221
Spectrum_tabulate(Spectrum me,bool includeBinNumbers,bool includeFrequency,bool includeRealPart,bool includeImaginaryPart,bool includeEnergyDensity,bool includePowerDensity)222 autoTable Spectrum_tabulate (Spectrum me, bool includeBinNumbers, bool includeFrequency,
223 bool includeRealPart, bool includeImaginaryPart, bool includeEnergyDensity, bool includePowerDensity)
224 {
225 try {
226 autoTable thee = Table_createWithoutColumnNames (my nx,
227 includeBinNumbers + includeFrequency + includeRealPart + includeImaginaryPart + includeEnergyDensity + includePowerDensity);
228 integer icol = 0;
229 if (includeBinNumbers)
230 Table_setColumnLabel (thee.get(), ++ icol, U"bin");
231 if (includeFrequency)
232 Table_setColumnLabel (thee.get(), ++ icol, U"freq(Hz)");
233 if (includeRealPart)
234 Table_setColumnLabel (thee.get(), ++ icol, U"re(Pa/Hz)");
235 if (includeImaginaryPart)
236 Table_setColumnLabel (thee.get(), ++ icol, U"im(Pa/Hz)");
237 if (includeEnergyDensity)
238 Table_setColumnLabel (thee.get(), ++ icol, U"energy(Pa^2/Hz^2)");
239 if (includePowerDensity)
240 Table_setColumnLabel (thee.get(), ++ icol, U"pow(dB/Hz)");
241 for (integer ibin = 1; ibin <= my nx; ibin ++) {
242 icol = 0;
243 if (includeBinNumbers)
244 Table_setNumericValue (thee.get(), ibin, ++ icol, ibin);
245 if (includeFrequency)
246 Table_setNumericValue (thee.get(), ibin, ++ icol, my x1 + (ibin - 1) * my dx);
247 if (includeRealPart)
248 Table_setNumericValue (thee.get(), ibin, ++ icol, my z [1] [ibin]);
249 if (includeImaginaryPart)
250 Table_setNumericValue (thee.get(), ibin, ++ icol, my z [2] [ibin]);
251 if (includeEnergyDensity)
252 Table_setNumericValue (thee.get(), ibin, ++ icol, Sampled_getValueAtSample (me, ibin, 0, 1));
253 if (includePowerDensity)
254 Table_setNumericValue (thee.get(), ibin, ++ icol, Sampled_getValueAtSample (me, ibin, 0, 2));
255 }
256 return thee;
257 } catch (MelderError) {
258 Melder_throw (me, U": not converted to Table.");
259 }
260 }
261
Spectrum_tabulate_verbose(Spectrum me)262 autoTable Spectrum_tabulate_verbose (Spectrum me) {
263 try {
264 static const conststring32 columnNames [] = { U"bin", U"frequency(Hz)", U"re(Pa/Hz)", U"im(Pa/Hz)",
265 U"energySpectralDensity(Pa^2s/Hz)", U"startOfBinWithinDomain(Hz)", U"endOfBinWithinDomain(Hz)", U"binWidthWithinDomain(Hz)", U"binEnergy(Pa^2s)",
266 U"powerSpectralDensity(Pa^2/Hz)", U"powerSpectralDensityInAir(W/m^2/Hz)", U"auditorySpectralDensityLevel(dB/Hz)",
267 U"binPower(Pa^2)", U"binPowerInAir(W/m^2)"
268 };
269 autoTable thee = Table_createWithColumnNames (my nx, ARRAY_TO_STRVEC (columnNames));
270 for (integer ibin = 1; ibin <= my nx; ibin ++) {
271 Table_setNumericValue (thee.get(), ibin, 1, ibin); // column "bin"
272 const double frequency = my x1 + (ibin - 1) * my dx;
273 Table_setNumericValue (thee.get(), ibin, 2, frequency); // column "frequency(Hz)"
274 const double re = my z [1] [ibin], im = my z [2] [ibin];
275 Table_setNumericValue (thee.get(), ibin, 3, re); // column "re(Pa/Hz)"
276 Table_setNumericValue (thee.get(), ibin, 4, im); // column "im(Pa/Hz)"
277 const double energySpectralDensity = Sampled_getValueAtSample (me, ibin, 0, 1);
278 Table_setNumericValue (thee.get(), ibin, 5, energySpectralDensity); // column "energySpectralDensity(Pa^2s/Hz)"
279 /*
280 The energy in a bin is integrated over only that part of the bin that falls within the frequency domain.
281 Practically speaking, the integration time (bin width) is my dx for most bins,
282 but almost always 0.5 * my dx for the first bin,
283 and also 0.5 * my dx for the last bin if the Nyquist falls within the bin
284 (i.e. if the original number of samples was even).
285 */
286 Melder_assert (my xmax >= my xmin); // required for Melder_clipped
287 const double startOfBinWithinDomain = Melder_clipped (my xmin, frequency - 0.5 * my dx, my xmax);
288 Table_setNumericValue (thee.get(), ibin, 6, startOfBinWithinDomain); // column "startOfBinWithinDomain(Hz)"
289 const double endOfBinWithinDomain = Melder_clipped (my xmin, frequency + 0.5 * my dx, my xmax);
290 Table_setNumericValue (thee.get(), ibin, 7, endOfBinWithinDomain); // column "endOfBinWithinDomain(Hz)"
291 const double binWidthWithinDomain = endOfBinWithinDomain - startOfBinWithinDomain;
292 Table_setNumericValue (thee.get(), ibin, 8, binWidthWithinDomain); // column "binWidthWithinDomain(Hz)"
293 Table_setNumericValue (thee.get(), ibin, 9, energySpectralDensity * binWidthWithinDomain); // column "binEnergy(Pa^2s)"
294 const double inverseOfEstimatedOriginalDuration = my dx;
295 const double powerSpectralDensity = energySpectralDensity * inverseOfEstimatedOriginalDuration;
296 Table_setNumericValue (thee.get(), ibin, 10, powerSpectralDensity); // column "powerSpectralDensity(Pa^2/Hz)"
297 constexpr double RHO_C = 400.0;
298 const double powerSpectralDensityInAir = powerSpectralDensity / RHO_C;
299 Table_setNumericValue (thee.get(), ibin, 11, powerSpectralDensityInAir); // column "powerSpectralDensityInAir(W/m^2/Hz)"
300 const double auditorySpectralDensityLevel = Sampled_getValueAtSample (me, ibin, 0, 2);
301 Table_setNumericValue (thee.get(), ibin, 12, auditorySpectralDensityLevel); // column "auditorySpectralDensityLevel(dB/Hz)"
302 Table_setNumericValue (thee.get(), ibin, 13, powerSpectralDensity * binWidthWithinDomain); // column "binPower(Pa^2)"
303 Table_setNumericValue (thee.get(), ibin, 14, powerSpectralDensity * binWidthWithinDomain / RHO_C); // column "binPowerInAir(W/m^2)"
304 }
305 return thee;
306 } catch (MelderError) {
307 Melder_throw (me, U": not converted to Table.");
308 }
309 }
310
Matrix_to_Spectrum(Matrix me)311 autoSpectrum Matrix_to_Spectrum (Matrix me) {
312 try {
313 if (my ny != 2)
314 Melder_throw (U"The Matrix should have exactly 2 rows.");
315 autoSpectrum thee = Thing_new (Spectrum);
316 my structMatrix :: v_copy (thee.get());
317 return thee;
318 } catch (MelderError) {
319 Melder_throw (me, U": not converted to Spectrum.");
320 }
321 }
322
Spectrum_to_Matrix(Spectrum me)323 autoMatrix Spectrum_to_Matrix (Spectrum me) {
324 try {
325 autoMatrix thee = Thing_new (Matrix);
326 my structMatrix :: v_copy (thee.get());
327 return thee;
328 } catch (MelderError) {
329 Melder_throw (me, U": not converted to Matrix.");
330 }
331 }
332
Spectrum_cepstralSmoothing(Spectrum me,double bandWidth)333 autoSpectrum Spectrum_cepstralSmoothing (Spectrum me, double bandWidth) {
334 try {
335 /*
336 dB-spectrum is log (power).
337 */
338 autoSpectrum dBspectrum = Data_copy (me);
339 VEC re = dBspectrum -> z.row (1), im = dBspectrum -> z.row (2);
340 for (integer i = 1; i <= dBspectrum -> nx; i ++) {
341 re [i] = log (sqr (re [i]) + sqr (im [i]) + 1e-308);
342 im [i] = 0.0;
343 }
344
345 /*
346 Cepstrum is Fourier transform of dB-spectrum.
347 */
348 autoSound cepstrum = Spectrum_to_Sound (dBspectrum.get());
349
350 /*
351 Multiply cepstrum by a Gaussian.
352 */
353 const double factor = - bandWidth * bandWidth;
354 for (integer i = 1; i <= cepstrum -> nx; i ++) {
355 const double t = (i - 1) * cepstrum -> dx;
356 cepstrum -> z [1] [i] *= exp (factor * sqr (t)) * ( i == 1 ? 1.0 : 2.0 );
357 }
358
359 /*
360 Smoothed power spectrum is original power spectrum convolved with a Gaussian.
361 */
362 autoSpectrum thee = Sound_to_Spectrum (cepstrum.get(), true);
363
364 /*
365 Convert power spectrum back into a "complex" spectrum without phase information.
366 */
367 re = thy z.row (1);
368 im = thy z.row (2);
369 for (integer i = 1; i <= thy nx; i ++) {
370 re [i] = exp (0.5 * re [i]); // i.e., sqrt (exp (re [i]))
371 im [i] = 0.0;
372 }
373 return thee;
374 } catch (MelderError) {
375 Melder_throw (me, U": cepstral smoothing not computed.");
376 }
377 }
378
Spectrum_passHannBand(Spectrum me,double fmin,double fmax0,double smooth)379 void Spectrum_passHannBand (Spectrum me, double fmin, double fmax0, double smooth) {
380 const double fmax = ( fmax0 == 0.0 ? my xmax : fmax0 );
381 const double f1 = fmin - smooth, f2 = fmin + smooth, f3 = fmax - smooth, f4 = fmax + smooth;
382 const double halfpibysmooth = ( smooth != 0.0 ? NUMpi / (2.0 * smooth) : 0.0 );
383 const VEC re = my z.row (1), im = my z.row (2);
384 for (integer i = 1; i <= my nx; i ++) {
385 const double frequency = my x1 + (i - 1) * my dx;
386 if (frequency < f1 || frequency > f4)
387 re [i] = im [i] = 0.0;
388 if (frequency < f2 && fmin > 0.0) {
389 const double factor = 0.5 - 0.5 * cos (halfpibysmooth * (frequency - f1));
390 re [i] *= factor;
391 im [i] *= factor;
392 } else if (frequency > f3 && fmax < my xmax) {
393 const double factor = 0.5 + 0.5 * cos (halfpibysmooth * (frequency - f3));
394 re [i] *= factor;
395 im [i] *= factor;
396 }
397 }
398 }
399
Spectrum_stopHannBand(Spectrum me,double fmin,double fmax0,double smooth)400 void Spectrum_stopHannBand (Spectrum me, double fmin, double fmax0, double smooth) {
401 const double fmax = ( fmax0 == 0.0 ? my xmax : fmax0 );
402 const double f1 = fmin - smooth, f2 = fmin + smooth, f3 = fmax - smooth, f4 = fmax + smooth;
403 const double halfpibysmooth = ( smooth != 0.0 ? NUMpi / (2.0 * smooth) : 0.0 );
404 const VEC re = my z.row (1), im = my z.row (2);
405 for (integer i = 1; i <= my nx; i ++) {
406 const double frequency = my x1 + (i - 1) * my dx;
407 if (frequency < f1 || frequency > f4)
408 continue;
409 if (frequency < f2 && fmin > 0.0) {
410 const double factor = 0.5 + 0.5 * cos (halfpibysmooth * (frequency - f1));
411 re [i] *= factor;
412 im [i] *= factor;
413 } else if (frequency > f3 && fmax < my xmax) {
414 const double factor = 0.5 - 0.5 * cos (halfpibysmooth * (frequency - f3));
415 re [i] *= factor;
416 im [i] *= factor;
417 } else
418 re [i] = im [i] = 0.0;
419 }
420 }
421
Spectrum_getBandEnergy(Spectrum me,double fmin,double fmax)422 double Spectrum_getBandEnergy (Spectrum me, double fmin, double fmax) {
423 /*
424 The computation requires that the negative-frequency values are the complex conjugates
425 of the positive-frequency values, and that only the positive-frequency values are included
426 in the spectrum.
427 */
428 if (my xmin < 0.0)
429 return undefined;
430 /*
431 Any energy outside [my xmin, my xmax] is ignored.
432 This is very important, since my xmin and my xmax determine the meaning of the first and last bins; see below.
433
434 The width of most bins is my dx, but the first and last bins can be truncated by fmin and fmax.
435
436 This truncation even applies in the case of autowindowing,
437 i.e. if the total energy in the spectrum is computed.
438 In that case, the first and last bins can be truncated by my xmin and my xmax.
439 This happens almost always for the first bin,
440 which is usually centred at f=0, hence has a width of 0.5 * my dx,
441 and quite often for the last bin as well (namely if the original sound had an even number of samples),
442 which is then centred at f=nyquist, hence has a width of 0.5 * my dx.
443
444 All this truncation is automatically performed by Sampled_getIntegral ().
445 */
446 return Sampled_getIntegral (me, fmin, fmax, 0, 1, false);
447 }
448
Spectrum_getBandDensity(Spectrum me,double fmin,double fmax)449 double Spectrum_getBandDensity (Spectrum me, double fmin, double fmax) {
450 if (my xmin < 0.0)
451 return undefined; // no negative frequencies allowed in one-sided spectral density
452 return Sampled_getMean (me, fmin, fmax, 0, 1, false);
453 }
454
Spectrum_getBandDensityDifference(Spectrum me,double lowBandMin,double lowBandMax,double highBandMin,double highBandMax)455 double Spectrum_getBandDensityDifference (Spectrum me, double lowBandMin, double lowBandMax, double highBandMin, double highBandMax) {
456 const double lowBandDensity = Spectrum_getBandDensity (me, lowBandMin, lowBandMax);
457 const double highBandDensity = Spectrum_getBandDensity (me, highBandMin, highBandMax);
458 if (isundef (lowBandDensity) || isundef (highBandDensity))
459 return undefined;
460 if (lowBandDensity == 0.0 || highBandDensity == 0.0)
461 return undefined;
462 return 10.0 * log10 (highBandDensity / lowBandDensity);
463 }
464
Spectrum_getBandEnergyDifference(Spectrum me,double lowBandMin,double lowBandMax,double highBandMin,double highBandMax)465 double Spectrum_getBandEnergyDifference (Spectrum me, double lowBandMin, double lowBandMax, double highBandMin, double highBandMax) {
466 const double lowBandEnergy = Spectrum_getBandEnergy (me, lowBandMin, lowBandMax);
467 const double highBandEnergy = Spectrum_getBandEnergy (me, highBandMin, highBandMax);
468 if (isundef (lowBandEnergy) || isundef (highBandEnergy))
469 return undefined;
470 if (lowBandEnergy == 0.0 || highBandEnergy == 0.0)
471 return undefined;
472 return 10.0 * log10 (highBandEnergy / lowBandEnergy);
473 }
474
Spectrum_getCentreOfGravity(Spectrum me,double power)475 double Spectrum_getCentreOfGravity (Spectrum me, double power) {
476 const double halfpower = 0.5 * power;
477 longdouble sumenergy = 0.0, sumfenergy = 0.0;
478 for (integer i = 1; i <= my nx; i ++) {
479 const double re = my z [1] [i], im = my z [2] [i];
480 double energy = sqr (re) + sqr (im);
481 const double f = my x1 + (i - 1) * my dx;
482 if (halfpower != 1.0)
483 energy = pow (energy, halfpower);
484 sumenergy += energy;
485 sumfenergy += f * energy;
486 }
487 return sumenergy == 0.0 ? undefined : double (sumfenergy / sumenergy);
488 }
489
Spectrum_getCentralMoment(Spectrum me,double moment,double power)490 double Spectrum_getCentralMoment (Spectrum me, double moment, double power) {
491 const double fmean = Spectrum_getCentreOfGravity (me, power);
492 if (isundef (fmean))
493 return undefined;
494 const double halfpower = 0.5 * power;
495 longdouble sumenergy = 0.0, sumfenergy = 0.0;
496 for (integer i = 1; i <= my nx; i ++) {
497 const double re = my z [1] [i], im = my z [2] [i];
498 double energy = sqr (re) + sqr (im);
499 const double f = my x1 + (i - 1) * my dx;
500 if (halfpower != 1.0)
501 energy = pow (energy, halfpower);
502 sumenergy += energy;
503 sumfenergy += pow (f - fmean, moment) * energy;
504 }
505 return double (sumfenergy / sumenergy);
506 }
507
Spectrum_getStandardDeviation(Spectrum me,double power)508 double Spectrum_getStandardDeviation (Spectrum me, double power) {
509 return sqrt (Spectrum_getCentralMoment (me, 2.0, power));
510 }
511
Spectrum_getSkewness(Spectrum me,double power)512 double Spectrum_getSkewness (Spectrum me, double power) {
513 const double m2 = Spectrum_getCentralMoment (me, 2.0, power);
514 const double m3 = Spectrum_getCentralMoment (me, 3.0, power);
515 if (isundef (m2) || isundef (m3) || m2 == 0.0)
516 return undefined;
517 return m3 / (m2 * sqrt (m2));
518 }
519
Spectrum_getKurtosis(Spectrum me,double power)520 double Spectrum_getKurtosis (Spectrum me, double power) {
521 const double m2 = Spectrum_getCentralMoment (me, 2.0, power);
522 const double m4 = Spectrum_getCentralMoment (me, 4.0, power);
523 if (isundef (m2) || isundef (m4) || m2 == 0.0)
524 return undefined;
525 return m4 / (m2 * m2) - 3.0;
526 }
527
Spectrum_getNearestMaximum(Spectrum me,double frequency)528 MelderPoint Spectrum_getNearestMaximum (Spectrum me, double frequency) {
529 try {
530 autoSpectrumTier thee = Spectrum_to_SpectrumTier_peaks (me);
531 const integer index = AnyTier_timeToNearestIndex (thee.get()->asAnyTier(), frequency);
532 if (index == 0)
533 Melder_throw (U"No peak.");
534 RealPoint point = thy points.at [index];
535 MelderPoint result;
536 result. x = point -> number;
537 result. y = point -> value;
538 return result;
539 } catch (MelderError) {
540 Melder_throw (me, U": no nearest maximum found.");
541 }
542 }
543
544 /* End of file Spectrum.cpp */
545