1 /* Sampled.cpp
2 *
3 * Copyright (C) 1992-2005,2007,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 #include "Sampled.h"
20
21 #include "oo_DESTROY.h"
22 #include "Sampled_def.h"
23 #include "oo_COPY.h"
24 #include "Sampled_def.h"
25 #include "oo_EQUAL.h"
26 #include "Sampled_def.h"
27 #include "oo_CAN_WRITE_AS_ENCODING.h"
28 #include "Sampled_def.h"
29 #include "oo_WRITE_TEXT.h"
30 #include "Sampled_def.h"
31 #include "oo_READ_TEXT.h"
32 #include "Sampled_def.h"
33 #include "oo_WRITE_BINARY.h"
34 #include "Sampled_def.h"
35 #include "oo_READ_BINARY.h"
36 #include "Sampled_def.h"
37 #include "oo_DESCRIPTION.h"
38 #include "Sampled_def.h"
39
40 Thing_implement (Sampled, Function, 0);
41
v_shiftX(double xfrom,double xto)42 void structSampled :: v_shiftX (double xfrom, double xto) {
43 Sampled_Parent :: v_shiftX (xfrom, xto);
44 NUMshift (& our x1, xfrom, xto);
45 }
46
v_scaleX(double xminfrom,double xmaxfrom,double xminto,double xmaxto)47 void structSampled :: v_scaleX (double xminfrom, double xmaxfrom, double xminto, double xmaxto) {
48 Sampled_Parent :: v_scaleX (xminfrom, xmaxfrom, xminto, xmaxto);
49 NUMscale (& our x1, xminfrom, xmaxfrom, xminto, xmaxto);
50 our dx *= (xmaxto - xminto) / (xmaxfrom - xminfrom);
51 }
52
Sampled_getWindowSamples(Sampled me,double xmin,double xmax,integer * ixmin,integer * ixmax)53 integer Sampled_getWindowSamples (Sampled me, double xmin, double xmax, integer *ixmin, integer *ixmax) {
54 const double ixmin_real = 1.0 + Melder_roundUp ((xmin - my x1) / my dx);
55 const double ixmax_real = 1.0 + Melder_roundDown ((xmax - my x1) / my dx); // could be above 32-bit LONG_MAX
56 *ixmin = ( ixmin_real < 1.0 ? 1 : (integer) ixmin_real );
57 *ixmax = ( ixmax_real > (double) my nx ? my nx : (integer) ixmax_real );
58 if (*ixmin > *ixmax)
59 return 0;
60 return *ixmax - *ixmin + 1;
61 }
62
Sampled_init(Sampled me,double xmin,double xmax,integer nx,double dx,double x1)63 void Sampled_init (Sampled me, double xmin, double xmax, integer nx, double dx, double x1) {
64 my xmin = xmin;
65 my xmax = xmax;
66 my nx = nx;
67 my dx = dx;
68 my x1 = x1;
69 }
70
Sampled_shortTermAnalysis(Sampled me,double windowDuration,double timeStep,integer * numberOfFrames,double * firstTime)71 void Sampled_shortTermAnalysis (Sampled me, double windowDuration, double timeStep, integer *numberOfFrames, double *firstTime) {
72 Melder_assert (windowDuration > 0.0);
73 Melder_assert (timeStep > 0.0);
74 volatile double myDuration = my dx * my nx;
75 if (windowDuration > myDuration)
76 Melder_throw (me, U": shorter than window length.");
77 *numberOfFrames = Melder_ifloor ((myDuration - windowDuration) / timeStep) + 1;
78 Melder_assert (*numberOfFrames >= 1);
79 const double ourMidTime = my x1 - 0.5 * my dx + 0.5 * myDuration;
80 const double thyDuration = *numberOfFrames * timeStep;
81 *firstTime = ourMidTime - 0.5 * thyDuration + 0.5 * timeStep;
82 }
83
Sampled_getValueAtSample(Sampled me,integer sampleNumber,integer levelNumber,int unit)84 double Sampled_getValueAtSample (Sampled me, integer sampleNumber, integer levelNumber, int unit) {
85 if (sampleNumber < 1 || sampleNumber > my nx)
86 return undefined;
87 return my v_getValueAtSample (sampleNumber, levelNumber, unit);
88 }
89
Sampled_listValuesOfAllSamples(Sampled me,integer levelNumber,int unit)90 autoVEC Sampled_listValuesOfAllSamples (Sampled me, integer levelNumber, int unit) {
91 autoVEC result = raw_VEC (my nx);
92 for (integer isamp = 1; isamp <= my nx; isamp ++)
93 result [isamp] = my v_getValueAtSample (isamp, levelNumber, unit);
94 return result;
95 }
96
Sampled_listValuesAtXes(Sampled me,constVECVU const & xes,integer levelNumber,int unit,bool interpolate)97 autoVEC Sampled_listValuesAtXes (Sampled me, constVECVU const& xes, integer levelNumber, int unit, bool interpolate) {
98 autoVEC result = raw_VEC (xes.size);
99 for (integer ix = 1; ix <= xes.size; ix ++)
100 result [ix] = Sampled_getValueAtX (me, xes [ix], levelNumber, unit, interpolate);
101 return result;
102 }
103
Sampled_getValueAtX(Sampled me,double x,integer levelNumber,int unit,bool interpolate)104 double Sampled_getValueAtX (Sampled me, double x, integer levelNumber, int unit, bool interpolate) {
105 if (x < my xmin || x > my xmax)
106 return undefined;
107 if (interpolate) {
108 const double index_real = Sampled_xToIndex (me, x);
109 const integer leftIndex = Melder_ifloor (index_real);
110 integer nearIndex, farIndex;
111 double phase = index_real - leftIndex;
112 if (phase < 0.5) {
113 nearIndex = leftIndex;
114 farIndex = leftIndex + 1;
115 } else {
116 farIndex = leftIndex;
117 nearIndex = leftIndex + 1;
118 phase = 1.0 - phase; // so that 0.0 < phase <= 0.5
119 }
120 if (nearIndex < 1 || nearIndex > my nx) // x out of range?
121 return undefined;
122 const double nearValue = my v_getValueAtSample (nearIndex, levelNumber, unit);
123 if (isundef (nearValue)) // function value not defined?
124 return undefined;
125 if (farIndex < 1 || farIndex > my nx) // at edge?
126 return nearValue; // extrapolate
127 const double farValue = my v_getValueAtSample (farIndex, levelNumber, unit);
128 if (isundef (farValue)) // neighbour undefined?
129 return nearValue; // extrapolate
130 return nearValue + phase * (farValue - nearValue); // interpolate
131 }
132 return Sampled_getValueAtSample (me, Sampled_xToNearestIndex (me, x), levelNumber, unit);
133 }
134
autoWindowDomainSamples(Sampled me,double * inout_xmin,double * inout_xmax,integer * out_imin,integer * out_imax)135 static integer autoWindowDomainSamples (Sampled me, double *inout_xmin, double *inout_xmax, integer *out_imin, integer *out_imax) {
136 Function_unidirectionalAutowindow (me, inout_xmin, inout_xmax);
137 if (! Function_intersectRangeWithDomain (me, inout_xmin, inout_xmax))
138 return 0;
139 return Sampled_getWindowSamples (me, *inout_xmin, *inout_xmax, out_imin, out_imax);
140 }
141
Sampled_countDefinedSamples(Sampled me,double xmin,double xmax,integer levelNumber,int unit)142 integer Sampled_countDefinedSamples (Sampled me, double xmin, double xmax, integer levelNumber, int unit) {
143 integer imin, imax;
144 integer numberOfFrames = autoWindowDomainSamples (me, & xmin, & xmax, & imin, & imax);
145 if (numberOfFrames < 1)
146 return 0;
147 integer numberOfDefinedSamples = 0;
148 for (integer isamp = imin; isamp <= imax; isamp ++) {
149 const double value = my v_getValueAtSample (isamp, levelNumber, unit);
150 if (isdefined (value))
151 numberOfDefinedSamples += 1;
152 }
153 return numberOfDefinedSamples;
154 }
155
Sampled_getSortedValues(Sampled me,double xmin,double xmax,integer levelNumber,int unit)156 autoVEC Sampled_getSortedValues (Sampled me, double xmin, double xmax, integer levelNumber, int unit) {
157 integer numberOfDefinedSamples = Sampled_countDefinedSamples (me, xmin, xmax, levelNumber, unit);
158 if (numberOfDefinedSamples == 0)
159 return autoVEC();
160 autoVEC definedValues = raw_VEC (numberOfDefinedSamples);
161 integer imin, imax;
162 autoWindowDomainSamples (me, & xmin, & xmax, & imin, & imax);
163 integer definedSampleNumber = 0;
164 for (integer isamp = imin; isamp <= imax; isamp ++) {
165 const double value = my v_getValueAtSample (isamp, levelNumber, unit);
166 if (isdefined (value))
167 definedValues [++ definedSampleNumber] = value;
168 }
169 Melder_assert (definedSampleNumber == numberOfDefinedSamples);
170 sort_VEC_inout (definedValues.get());
171 return definedValues;
172 }
173
Sampled_getQuantile(Sampled me,double xmin,double xmax,double quantile,integer levelNumber,int unit)174 double Sampled_getQuantile (Sampled me, double xmin, double xmax, double quantile, integer levelNumber, int unit) {
175 try {
176 autoVEC values = Sampled_getSortedValues (me, xmin, xmax, levelNumber, unit);
177 return NUMquantile (values.get(), quantile);
178 } catch (MelderError) {
179 Melder_throw (me, U": quantile not computed.");
180 }
181 }
182
Sampled_getSumAndDefinitionRange(Sampled me,double xmin,double xmax,integer levelNumber,int unit,bool interpolate,double * out_sum,double * out_definitionRange)183 static void Sampled_getSumAndDefinitionRange
184 (Sampled me, double xmin, double xmax, integer levelNumber, int unit, bool interpolate, double *out_sum, double *out_definitionRange)
185 {
186 /*
187 This function computes the area under the linearly interpolated curve between xmin and xmax.
188 Outside [x1-dx/2, xN+dx/2], the curve is undefined and neither times nor values are counted.
189 In [x1-dx/2,x1] and [xN,xN+dx/2], the curve is linearly extrapolated.
190 */
191 longdouble sum = 0.0, definitionRange = 0.0;
192 Function_unidirectionalAutowindow (me, & xmin, & xmax);
193 if (Function_intersectRangeWithDomain (me, & xmin, & xmax)) {
194 if (interpolate) {
195 integer imin, imax;
196 if (Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax) > 0) {
197 double leftEdge = my x1 - 0.5 * my dx, rightEdge = leftEdge + my nx * my dx;
198 for (integer isamp = imin; isamp <= imax; isamp ++) {
199 const double value = my v_getValueAtSample (isamp, levelNumber, unit); // a fast way to integrate a linearly interpolated curve; works everywhere except at the edges
200 if (isdefined (value)) {
201 definitionRange += 1.0;
202 sum += value;
203 }
204 }
205 /*
206 Corrections within the first and last sampling intervals.
207 */
208 if (xmin > leftEdge) { // otherwise, constant extrapolation over 0.5 sample is OK
209 double phase = (my x1 + (imin - 1) * my dx - xmin) / my dx; // this fraction of sampling interval is still to be determined
210 const double rightValue = Sampled_getValueAtSample (me, imin, levelNumber, unit);
211 const double leftValue = Sampled_getValueAtSample (me, imin - 1, levelNumber, unit);
212 if (isdefined (rightValue)) {
213 definitionRange -= 0.5; // delete constant extrapolation over 0.5 sample
214 sum -= 0.5 * rightValue;
215 if (isdefined (leftValue)) {
216 definitionRange += phase; // add current fraction
217 sum += phase * (rightValue + 0.5 * phase * (leftValue - rightValue)); // interpolate to outside sample
218 } else {
219 if (phase > 0.5)
220 phase = 0.5;
221 definitionRange += phase; // add current fraction, but never more than 0.5
222 sum += phase * rightValue;
223 }
224 } else if (isdefined (leftValue) && phase > 0.5) {
225 definitionRange += phase - 0.5;
226 sum += (phase - 0.5) * leftValue;
227 }
228 }
229 if (xmax < rightEdge) { // otherwise, constant extrapolation is OK
230 double phase = (xmax - (my x1 + (imax - 1) * my dx)) / my dx; // this fraction of sampling interval is still to be determined
231 const double leftValue = Sampled_getValueAtSample (me, imax, levelNumber, unit);
232 const double rightValue = Sampled_getValueAtSample (me, imax + 1, levelNumber, unit);
233 if (isdefined (leftValue)) {
234 definitionRange -= 0.5; // delete constant extrapolation over 0.5 sample
235 sum -= 0.5 * leftValue;
236 if (isdefined (rightValue)) {
237 definitionRange += phase; // add current fraction
238 sum += phase * (leftValue + 0.5 * phase * (rightValue - leftValue)); // interpolate to outside sample
239 } else {
240 if (phase > 0.5)
241 phase = 0.5;
242 definitionRange += phase; // add current fraction, but never more than 0.5
243 sum += phase * leftValue;
244 }
245 } else if (isdefined (rightValue) && phase > 0.5) {
246 definitionRange += phase - 0.5;
247 sum += (phase - 0.5) * rightValue;
248 }
249 }
250 } else { // no sample centres between xmin and xmax
251 /*
252 Try to return the mean of the interpolated values at these two points.
253 Thus, a small (xmin, xmax) range gives the same value as the (xmin+xmax)/2 point.
254 */
255 const double leftValue = Sampled_getValueAtSample (me, imax, levelNumber, unit);
256 const double rightValue = Sampled_getValueAtSample (me, imin, levelNumber, unit);
257 double phase1 = (xmin - (my x1 + (imax - 1) * my dx)) / my dx;
258 double phase2 = (xmax - (my x1 + (imax - 1) * my dx)) / my dx;
259 if (imin == imax + 1) { // not too far from sample definition region
260 if (isdefined (leftValue)) {
261 if (isdefined (rightValue)) {
262 definitionRange += phase2 - phase1;
263 sum += (phase2 - phase1) * (leftValue + 0.5 * (phase1 + phase2) * (rightValue - leftValue));
264 } else if (phase1 < 0.5) {
265 if (phase2 > 0.5)
266 phase2 = 0.5;
267 definitionRange += phase2 - phase1;
268 sum += (phase2 - phase1) * leftValue;
269 }
270 } else if (isdefined (rightValue) && phase2 > 0.5) {
271 if (phase1 < 0.5)
272 phase1 = 0.5;
273 definitionRange += phase2 - phase1;
274 sum += (phase2 - phase1) * rightValue;
275 }
276 }
277 }
278 } else { // no interpolation
279 const double rimin = Sampled_xToIndex (me, xmin), rimax = Sampled_xToIndex (me, xmax);
280 if (rimax >= 0.5 && rimin < my nx + 0.5) {
281 const integer imin = ( rimin < 0.5 ? 0 : Melder_iround (rimin) );
282 const integer imax = ( rimax >= my nx + 0.5 ? my nx + 1 : Melder_iround (rimax) );
283 for (integer isamp = imin + 1; isamp < imax; isamp ++) {
284 const double value = my v_getValueAtSample (isamp, levelNumber, unit);
285 if (isdefined (value)) {
286 definitionRange += 1.0;
287 sum += value;
288 }
289 }
290 if (imin == imax) {
291 const double value = my v_getValueAtSample (imin, levelNumber, unit);
292 if (isdefined (value)) {
293 const double phase = rimax - rimin;
294 definitionRange += phase;
295 sum += phase * value;
296 }
297 } else {
298 if (imin >= 1) {
299 const double value = my v_getValueAtSample (imin, levelNumber, unit);
300 if (isdefined (value)) {
301 const double phase = imin - rimin + 0.5;
302 definitionRange += phase;
303 sum += phase * value;
304 }
305 }
306 if (imax <= my nx) {
307 const double value = my v_getValueAtSample (imax, levelNumber, unit);
308 if (isdefined (value)) {
309 const double phase = rimax - imax + 0.5;
310 definitionRange += phase;
311 sum += phase * value;
312 }
313 }
314 }
315 }
316 }
317 }
318 if (out_sum)
319 *out_sum = double (sum);
320 if (out_definitionRange)
321 *out_definitionRange = double (definitionRange);
322 }
323
Sampled_getMean(Sampled me,double xmin,double xmax,integer levelNumber,int unit,bool interpolate)324 double Sampled_getMean (Sampled me, double xmin, double xmax, integer levelNumber, int unit, bool interpolate) {
325 double sum, definitionRange;
326 Sampled_getSumAndDefinitionRange (me, xmin, xmax, levelNumber, unit, interpolate, & sum, & definitionRange);
327 return definitionRange <= 0.0 ? undefined : sum / definitionRange;
328 }
329
Sampled_getMean_standardUnit(Sampled me,double xmin,double xmax,integer levelNumber,int averagingUnit,bool interpolate)330 double Sampled_getMean_standardUnit (Sampled me, double xmin, double xmax, integer levelNumber, int averagingUnit, bool interpolate) {
331 const double mean = Sampled_getMean (me, xmin, xmax, levelNumber, averagingUnit, interpolate);
332 return Function_convertSpecialToStandardUnit (me, mean, levelNumber, averagingUnit);
333 }
334
Sampled_getIntegral(Sampled me,double xmin,double xmax,integer levelNumber,int unit,bool interpolate)335 double Sampled_getIntegral (Sampled me, double xmin, double xmax, integer levelNumber, int unit, bool interpolate) {
336 double sum, definitionRange;
337 Sampled_getSumAndDefinitionRange (me, xmin, xmax, levelNumber, unit, interpolate, & sum, & definitionRange);
338 return sum * my dx;
339 }
340
Sampled_getIntegral_standardUnit(Sampled me,double xmin,double xmax,integer levelNumber,int averagingUnit,bool interpolate)341 double Sampled_getIntegral_standardUnit (Sampled me, double xmin, double xmax, integer levelNumber, int averagingUnit, bool interpolate) {
342 const double integral = Sampled_getIntegral (me, xmin, xmax, levelNumber, averagingUnit, interpolate);
343 return Function_convertSpecialToStandardUnit (me, integral, levelNumber, averagingUnit);
344 }
345
Sampled_getSum2AndDefinitionRange(Sampled me,double xmin,double xmax,integer levelNumber,int unit,double mean,bool interpolate,double * out_sum2,double * out_definitionRange)346 static void Sampled_getSum2AndDefinitionRange
347 (Sampled me, double xmin, double xmax, integer levelNumber, int unit, double mean, bool interpolate, double *out_sum2, double *out_definitionRange)
348 {
349 /*
350 This function computes the area under the linearly interpolated squared difference curve between xmin and xmax.
351 Outside [x1-dx/2, xN+dx/2], the curve is undefined and neither times nor values are counted.
352 In [x1-dx/2,x1] and [xN,xN+dx/2], the curve is linearly extrapolated.
353 */
354 longdouble sum2 = 0.0, definitionRange = 0.0;
355 Function_unidirectionalAutowindow (me, & xmin, & xmax);
356 if (Function_intersectRangeWithDomain (me, & xmin, & xmax)) {
357 if (interpolate) {
358 integer imin, imax;
359 if (Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax) > 0) {
360 const double leftEdge = my x1 - 0.5 * my dx, rightEdge = leftEdge + my nx * my dx;
361 for (integer isamp = imin; isamp <= imax; isamp ++) {
362 double value = my v_getValueAtSample (isamp, levelNumber, unit); // a fast way to integrate a linearly interpolated curve; works everywhere except at the edges
363 if (isdefined (value)) {
364 value -= mean;
365 value *= value;
366 definitionRange += 1.0;
367 sum2 += value;
368 }
369 }
370 /*
371 Corrections within the first and last sampling intervals.
372 */
373 if (xmin > leftEdge) { // otherwise, constant extrapolation over 0.5 sample is OK
374 double phase = (my x1 + (imin - 1) * my dx - xmin) / my dx; // this fraction of sampling interval is still to be determined
375 double rightValue = Sampled_getValueAtSample (me, imin, levelNumber, unit);
376 double leftValue = Sampled_getValueAtSample (me, imin - 1, levelNumber, unit);
377 if (isdefined (rightValue)) {
378 rightValue -= mean;
379 rightValue *= rightValue;
380 definitionRange -= 0.5; // delete constant extrapolation over 0.5 sample
381 sum2 -= 0.5 * rightValue;
382 if (isdefined (leftValue)) {
383 leftValue -= mean;
384 leftValue *= leftValue;
385 definitionRange += phase; // add current fraction
386 sum2 += phase * (rightValue + 0.5 * phase * (leftValue - rightValue)); // interpolate to outside sample
387 } else {
388 if (phase > 0.5)
389 phase = 0.5;
390 definitionRange += phase; // add current fraction, but never more than 0.5
391 sum2 += phase * rightValue;
392 }
393 } else if (isdefined (leftValue) && phase > 0.5) {
394 leftValue -= mean;
395 leftValue *= leftValue;
396 definitionRange += phase - 0.5;
397 sum2 += (phase - 0.5) * leftValue;
398 }
399 }
400 if (xmax < rightEdge) { // otherwise, constant extrapolation is OK
401 double phase = (xmax - (my x1 + (imax - 1) * my dx)) / my dx; // this fraction of sampling interval is still to be determined
402 double leftValue = Sampled_getValueAtSample (me, imax, levelNumber, unit);
403 double rightValue = Sampled_getValueAtSample (me, imax + 1, levelNumber, unit);
404 if (isdefined (leftValue)) {
405 leftValue -= mean;
406 leftValue *= leftValue;
407 definitionRange -= 0.5; // delete constant extrapolation over 0.5 sample
408 sum2 -= 0.5 * leftValue;
409 if (isdefined (rightValue)) {
410 rightValue -= mean;
411 rightValue *= rightValue;
412 definitionRange += phase; // add current fraction
413 sum2 += phase * (leftValue + 0.5 * phase * (rightValue - leftValue)); // interpolate to outside sample
414 } else {
415 if (phase > 0.5)
416 phase = 0.5;
417 definitionRange += phase; // add current fraction, but never more than 0.5
418 sum2 += phase * leftValue;
419 }
420 } else if (isdefined (rightValue) && phase > 0.5) {
421 rightValue -= mean;
422 rightValue *= rightValue;
423 definitionRange += phase - 0.5;
424 sum2 += (phase - 0.5) * rightValue;
425 }
426 }
427 } else { // no sample centres between xmin and xmax
428 /*
429 Try to return the mean of the interpolated values at these two points.
430 Thus, a small (xmin, xmax) range gives the same value as the (xmin+xmax)/2 point.
431 */
432 double leftValue = Sampled_getValueAtSample (me, imax, levelNumber, unit);
433 double rightValue = Sampled_getValueAtSample (me, imin, levelNumber, unit);
434 double phase1 = (xmin - (my x1 + (imax - 1) * my dx)) / my dx;
435 double phase2 = (xmax - (my x1 + (imax - 1) * my dx)) / my dx;
436 if (imin == imax + 1) { // not too far from sample definition region
437 if (isdefined (leftValue)) {
438 leftValue -= mean;
439 leftValue *= leftValue;
440 if (isdefined (rightValue)) {
441 rightValue -= mean;
442 rightValue *= rightValue;
443 definitionRange += phase2 - phase1;
444 sum2 += (phase2 - phase1) * (leftValue + 0.5 * (phase1 + phase2) * (rightValue - leftValue));
445 } else if (phase1 < 0.5) {
446 if (phase2 > 0.5)
447 phase2 = 0.5;
448 definitionRange += phase2 - phase1;
449 sum2 += (phase2 - phase1) * leftValue;
450 }
451 } else if (isdefined (rightValue) && phase2 > 0.5) {
452 rightValue -= mean;
453 rightValue *= rightValue;
454 if (phase1 < 0.5)
455 phase1 = 0.5;
456 definitionRange += phase2 - phase1;
457 sum2 += (phase2 - phase1) * rightValue;
458 }
459 }
460 }
461 } else { // no interpolation
462 const double rimin = Sampled_xToIndex (me, xmin), rimax = Sampled_xToIndex (me, xmax);
463 if (rimax >= 0.5 && rimin < my nx + 0.5) {
464 const integer imin = rimin < 0.5 ? 0 : Melder_iround (rimin);
465 const integer imax = rimax >= my nx + 0.5 ? my nx + 1 : Melder_iround (rimax);
466 for (integer isamp = imin + 1; isamp < imax; isamp ++) {
467 double value = my v_getValueAtSample (isamp, levelNumber, unit);
468 if (isdefined (value)) {
469 value -= mean;
470 value *= value;
471 definitionRange += 1.0;
472 sum2 += value;
473 }
474 }
475 if (imin == imax) {
476 double value = my v_getValueAtSample (imin, levelNumber, unit);
477 if (isdefined (value)) {
478 const double phase = rimax - rimin;
479 value -= mean;
480 value *= value;
481 definitionRange += phase;
482 sum2 += phase * value;
483 }
484 } else {
485 if (imin >= 1) {
486 double value = my v_getValueAtSample (imin, levelNumber, unit);
487 if (isdefined (value)) {
488 const double phase = imin - rimin + 0.5;
489 value -= mean;
490 value *= value;
491 definitionRange += phase;
492 sum2 += phase * value;
493 }
494 }
495 if (imax <= my nx) {
496 double value = my v_getValueAtSample (imax, levelNumber, unit);
497 if (isdefined (value)) {
498 const double phase = rimax - imax + 0.5;
499 value -= mean;
500 value *= value;
501 definitionRange += phase;
502 sum2 += phase * value;
503 }
504 }
505 }
506 }
507 }
508 }
509 if (out_sum2)
510 *out_sum2 = double (sum2);
511 if (out_definitionRange)
512 *out_definitionRange = double (definitionRange);
513 }
514
Sampled_getStandardDeviation(Sampled me,double xmin,double xmax,integer levelNumber,int unit,bool interpolate)515 double Sampled_getStandardDeviation (Sampled me, double xmin, double xmax, integer levelNumber, int unit, bool interpolate) {
516 double sum, sum2, definitionRange;
517 Sampled_getSumAndDefinitionRange (me, xmin, xmax, levelNumber, unit, interpolate, & sum, & definitionRange);
518 if (definitionRange < 2.0)
519 return undefined;
520 Sampled_getSum2AndDefinitionRange (me, xmin, xmax, levelNumber, unit, sum / definitionRange, interpolate, & sum2, & definitionRange);
521 return sqrt (sum2 / (definitionRange - 1.0));
522 }
523
Sampled_getStandardDeviation_standardUnit(Sampled me,double xmin,double xmax,integer levelNumber,int averagingUnit,bool interpolate)524 double Sampled_getStandardDeviation_standardUnit (Sampled me, double xmin, double xmax, integer levelNumber, int averagingUnit, bool interpolate) {
525 const double stdev = Sampled_getStandardDeviation (me, xmin, xmax, levelNumber, averagingUnit, interpolate);
526 return Function_convertSpecialToStandardUnit (me, stdev, levelNumber, averagingUnit);
527 }
528
Sampled_getMinimumAndX(Sampled me,double xmin,double xmax,integer levelNumber,int unit,bool interpolate,double * return_minimum,double * return_xOfMinimum)529 void Sampled_getMinimumAndX (Sampled me, double xmin, double xmax, integer levelNumber, int unit, bool interpolate,
530 double *return_minimum, double *return_xOfMinimum)
531 {
532 double minimum = 1e301, xOfMinimum = 0.0;
533 if (isundef (xmin) || isundef (xmax)) {
534 minimum = xOfMinimum = undefined;
535 goto end;
536 }
537 Function_unidirectionalAutowindow (me, & xmin, & xmax);
538 if (! Function_intersectRangeWithDomain (me, & xmin, & xmax)) {
539 minimum = xOfMinimum = undefined; // requested range and logical domain do not intersect
540 goto end;
541 }
542 integer imin, imax;
543 if (! Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax)) {
544 /*
545 No sample centres between xmin and xmax.
546 Try to return the lesser of the values at these two points.
547 */
548 const double fleft = Sampled_getValueAtX (me, xmin, levelNumber, unit, interpolate);
549 const double fright = Sampled_getValueAtX (me, xmax, levelNumber, unit, interpolate);
550 if (isdefined (fleft) && fleft < minimum) {
551 minimum = fleft;
552 xOfMinimum = xmin;
553 }
554 if (isdefined (fright) && fright < minimum) {
555 minimum = fright;
556 xOfMinimum = xmax;
557 }
558 } else {
559 for (integer i = imin; i <= imax; i ++) {
560 const double fmid = my v_getValueAtSample (i, levelNumber, unit);
561 if (isundef (fmid)) continue;
562 if (! interpolate) {
563 if (fmid < minimum) {
564 minimum = fmid;
565 xOfMinimum = i;
566 }
567 } else {
568 /*
569 Try an interpolation, possibly even taking into account a sample just outside the selection.
570 */
571 const double fleft = ( i <= 1 ? undefined : my v_getValueAtSample (i - 1, levelNumber, unit) );
572 const double fright = ( i >= my nx ? undefined : my v_getValueAtSample (i + 1, levelNumber, unit) );
573 if (isundef (fleft) || isundef (fright)) {
574 if (fmid < minimum) {
575 minimum = fmid;
576 xOfMinimum = i;
577 }
578 } else if (fmid < fleft && fmid <= fright) {
579 const double y [] = { fleft, fmid, fright };
580 double i_real;
581 const double localMinimum = NUMimproveMinimum (ARRAY_TO_VEC (y), 2, NUM_PEAK_INTERPOLATE_PARABOLIC, & i_real);
582 if (localMinimum < minimum) {
583 minimum = localMinimum;
584 xOfMinimum = i_real + i - 2;
585 }
586 }
587 }
588 }
589 xOfMinimum = my x1 + (xOfMinimum - 1) * my dx; // from index plus phase to time
590 /*
591 Check boundary values.
592 */
593 if (interpolate) {
594 const double fleft = Sampled_getValueAtX (me, xmin, levelNumber, unit, true);
595 const double fright = Sampled_getValueAtX (me, xmax, levelNumber, unit, true);
596 if (isdefined (fleft) && fleft < minimum) {
597 minimum = fleft;
598 xOfMinimum = xmin;
599 }
600 if (isdefined (fright) && fright < minimum) {
601 minimum = fright;
602 xOfMinimum = xmax;
603 }
604 }
605 Melder_clip (xmin, & xOfMinimum, xmax);
606 }
607 if (minimum == 1e301)
608 minimum = xOfMinimum = undefined;
609 end:
610 if (return_minimum)
611 *return_minimum = minimum;
612 if (return_xOfMinimum)
613 *return_xOfMinimum = xOfMinimum;
614 }
615
Sampled_getMinimum(Sampled me,double xmin,double xmax,integer levelNumber,int unit,bool interpolate)616 double Sampled_getMinimum (Sampled me, double xmin, double xmax, integer levelNumber, int unit, bool interpolate) {
617 double minimum;
618 Sampled_getMinimumAndX (me, xmin, xmax, levelNumber, unit, interpolate, & minimum, nullptr);
619 return minimum;
620 }
621
Sampled_getXOfMinimum(Sampled me,double xmin,double xmax,integer levelNumber,int unit,bool interpolate)622 double Sampled_getXOfMinimum (Sampled me, double xmin, double xmax, integer levelNumber, int unit, bool interpolate) {
623 double time;
624 Sampled_getMinimumAndX (me, xmin, xmax, levelNumber, unit, interpolate, nullptr, & time);
625 return time;
626 }
627
Sampled_getMaximumAndX(Sampled me,double xmin,double xmax,integer levelNumber,int unit,bool interpolate,double * return_maximum,double * return_xOfMaximum)628 void Sampled_getMaximumAndX (Sampled me, double xmin, double xmax, integer levelNumber, int unit, bool interpolate,
629 double *return_maximum, double *return_xOfMaximum)
630 {
631 double maximum = -1e301, xOfMaximum = 0.0;
632 if (isundef (xmin) || isundef (xmax)) {
633 maximum = xOfMaximum = undefined;
634 goto end;
635 }
636 Function_unidirectionalAutowindow (me, & xmin, & xmax);
637 if (! Function_intersectRangeWithDomain (me, & xmin, & xmax)) {
638 maximum = xOfMaximum = undefined; // requested range and logical domain do not intersect
639 goto end;
640 }
641 integer imin, imax;
642 if (! Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax)) {
643 /*
644 No sample centres between tmin and tmax.
645 Try to return the greater of the values at these two points.
646 */
647 const double fleft = Sampled_getValueAtX (me, xmin, levelNumber, unit, interpolate);
648 const double fright = Sampled_getValueAtX (me, xmax, levelNumber, unit, interpolate);
649 if (isdefined (fleft) && fleft > maximum) {
650 maximum = fleft;
651 xOfMaximum = xmin;
652 }
653 if (isdefined (fright) && fright > maximum) {
654 maximum = fright;
655 xOfMaximum = xmax;
656 }
657 } else {
658 for (integer i = imin; i <= imax; i ++) {
659 const double fmid = my v_getValueAtSample (i, levelNumber, unit);
660 if (isundef (fmid)) continue;
661 if (! interpolate) {
662 if (fmid > maximum) {
663 maximum = fmid;
664 xOfMaximum = i;
665 }
666 } else {
667 /*
668 Try an interpolation, possibly even taking into account a sample just outside the selection.
669 */
670 const double fleft =
671 ( i <= 1 ? undefined : my v_getValueAtSample (i - 1, levelNumber, unit) );
672 const double fright =
673 ( i >= my nx ? undefined : my v_getValueAtSample (i + 1, levelNumber, unit) );
674 if (isundef (fleft) || isundef (fright)) {
675 if (fmid > maximum) {
676 maximum = fmid;
677 xOfMaximum = i;
678 }
679 } else if (fmid > fleft && fmid >= fright) {
680 const double y [] = { fleft, fmid, fright };
681 double i_real;
682 const double localMaximum = NUMimproveMaximum (ARRAY_TO_VEC (y), 2, NUM_PEAK_INTERPOLATE_PARABOLIC, & i_real);
683 if (localMaximum > maximum) {
684 maximum = localMaximum;
685 xOfMaximum = i_real + i - 2;
686 }
687 }
688 }
689 }
690 xOfMaximum = my x1 + (xOfMaximum - 1) * my dx; // from index plus phase to time
691 /*
692 Check boundary values.
693 */
694 if (interpolate) {
695 const double fleft = Sampled_getValueAtX (me, xmin, levelNumber, unit, true);
696 const double fright = Sampled_getValueAtX (me, xmax, levelNumber, unit, true);
697 if (isdefined (fleft) && fleft > maximum) {
698 maximum = fleft;
699 xOfMaximum = xmin;
700 }
701 if (isdefined (fright) && fright > maximum) {
702 maximum = fright;
703 xOfMaximum = xmax;
704 }
705 }
706 Melder_clip (xmin, & xOfMaximum, xmax);
707 }
708 if (maximum == -1e301)
709 maximum = xOfMaximum = undefined;
710 end:
711 if (return_maximum)
712 *return_maximum = maximum;
713 if (return_xOfMaximum)
714 *return_xOfMaximum = xOfMaximum;
715 }
716
Sampled_getMaximum(Sampled me,double xmin,double xmax,integer levelNumber,int unit,bool interpolate)717 double Sampled_getMaximum (Sampled me, double xmin, double xmax, integer levelNumber, int unit, bool interpolate) {
718 double maximum;
719 Sampled_getMaximumAndX (me, xmin, xmax, levelNumber, unit, interpolate, & maximum, nullptr);
720 return maximum;
721 }
722
Sampled_getXOfMaximum(Sampled me,double xmin,double xmax,integer levelNumber,int unit,bool interpolate)723 double Sampled_getXOfMaximum (Sampled me, double xmin, double xmax, integer levelNumber, int unit, bool interpolate) {
724 double time;
725 Sampled_getMaximumAndX (me, xmin, xmax, levelNumber, unit, interpolate, nullptr, & time);
726 return time;
727 }
728
Sampled_speckleInside(Sampled me,Graphics g,double xmin,double xmax,double ymin,double ymax,integer levelNumber,int unit)729 static void Sampled_speckleInside (Sampled me, Graphics g, double xmin, double xmax, double ymin, double ymax,
730 integer levelNumber, int unit)
731 {
732 Function_unidirectionalAutowindow (me, & xmin, & xmax);
733 integer ixmin, ixmax;
734 integer numberOfSamples = Sampled_getWindowSamples (me, xmin, xmax, & ixmin, & ixmax);
735 if (numberOfSamples <= 0)
736 return;
737 if (Function_isUnitLogarithmic (me, levelNumber, unit)) {
738 ymin = Function_convertStandardToSpecialUnit (me, ymin, levelNumber, unit);
739 ymax = Function_convertStandardToSpecialUnit (me, ymax, levelNumber, unit);
740 }
741 if (ymax <= ymin)
742 return;
743 Graphics_setWindow (g, xmin, xmax, ymin, ymax);
744 for (integer ix = ixmin; ix <= ixmax; ix ++) {
745 const double value = Sampled_getValueAtSample (me, ix, levelNumber, unit);
746 if (isdefined (value)) {
747 const double x = Sampled_indexToX (me, ix);
748 if (value >= ymin && value <= ymax)
749 Graphics_speckle (g, x, value);
750 }
751 }
752 }
753
Sampled_drawInside(Sampled me,Graphics g,double xmin,double xmax,double ymin,double ymax,bool speckle,integer levelNumber,int unit)754 void Sampled_drawInside (Sampled me, Graphics g, double xmin, double xmax, double ymin, double ymax,
755 bool speckle, integer levelNumber, int unit)
756 {
757 try {
758 if (speckle) {
759 Sampled_speckleInside (me, g, xmin, xmax, ymin, ymax, levelNumber, unit);
760 return;
761 }
762 Function_unidirectionalAutowindow (me, & xmin, & xmax);
763 integer ixmin, ixmax, startOfDefinedStretch = -1;
764 integer numberOfSamples = Sampled_getWindowSamples (me, xmin, xmax, & ixmin, & ixmax);
765 if (numberOfSamples <= 0)
766 return;
767 if (Function_isUnitLogarithmic (me, levelNumber, unit)) {
768 ymin = Function_convertStandardToSpecialUnit (me, ymin, levelNumber, unit);
769 ymax = Function_convertStandardToSpecialUnit (me, ymax, levelNumber, unit);
770 }
771 if (ymax <= ymin)
772 return;
773 Graphics_setWindow (g, xmin, xmax, ymin, ymax);
774 const integer lowIndex = ixmin - 1, highIndex = ixmax + 1;
775 const integer nbuffer = highIndex - lowIndex + 1;
776 autoVEC xbuffer = zero_VEC (nbuffer), ybuffer = zero_VEC (nbuffer);
777 const integer bufferShift = 1 - lowIndex;
778 double *xarray = & xbuffer [bufferShift];
779 double *yarray = & ybuffer [bufferShift];
780 double previousValue = Sampled_getValueAtSample (me, lowIndex, levelNumber, unit);
781 if (isdefined (previousValue)) {
782 startOfDefinedStretch = lowIndex;
783 xarray [lowIndex] = Sampled_indexToX (me, lowIndex);
784 yarray [lowIndex] = previousValue;
785 }
786 for (integer ix = ixmin; ix <= ixmax; ix ++) {
787 const double x = Sampled_indexToX (me, ix);
788 const double value = Sampled_getValueAtSample (me, ix, levelNumber, unit);
789 if (isdefined (value)) {
790 if (isdefined (previousValue)) {
791 xarray [ix] = x;
792 yarray [ix] = value;
793 } else {
794 startOfDefinedStretch = ix - 1;
795 xarray [ix - 1] = x - 0.5 * my dx;
796 yarray [ix - 1] = value;
797 xarray [ix] = x;
798 yarray [ix] = value;
799 }
800 } else if (isdefined (previousValue)) {
801 Melder_assert (startOfDefinedStretch >= lowIndex);
802 if (ix > ixmin) {
803 xarray [ix] = x - 0.5 * my dx;
804 yarray [ix] = previousValue;
805 if (xarray [startOfDefinedStretch] < xmin) {
806 const double phase = (xmin - xarray [startOfDefinedStretch]) / my dx;
807 xarray [startOfDefinedStretch] = xmin;
808 yarray [startOfDefinedStretch] = phase * yarray [startOfDefinedStretch + 1] + (1.0 - phase) * yarray [startOfDefinedStretch];
809 }
810 Graphics_polyline (g, ix + 1 - startOfDefinedStretch, & xarray [startOfDefinedStretch], & yarray [startOfDefinedStretch]);
811 }
812 startOfDefinedStretch = -1;
813 }
814 previousValue = value;
815 }
816 if (startOfDefinedStretch > -1) {
817 const double x = Sampled_indexToX (me, highIndex);
818 const double value = Sampled_getValueAtSample (me, ixmax + 1, levelNumber, unit);
819 Melder_assert (isdefined (previousValue));
820 if (isdefined (value)) {
821 xarray [highIndex] = x;
822 yarray [highIndex] = value;
823 } else {
824 xarray [highIndex] = x - 0.5 * my dx;
825 yarray [highIndex] = previousValue;
826 }
827 if (xarray [startOfDefinedStretch] < xmin) {
828 const double phase = (xmin - xarray [startOfDefinedStretch]) / my dx;
829 xarray [startOfDefinedStretch] = xmin;
830 yarray [startOfDefinedStretch] = phase * yarray [startOfDefinedStretch + 1] + (1.0 - phase) * yarray [startOfDefinedStretch];
831 }
832 if (xarray [highIndex] > xmax) {
833 const double phase = (xarray [highIndex] - xmax) / my dx;
834 xarray [highIndex] = xmax;
835 yarray [highIndex] = phase * yarray [ixmax] + (1.0 - phase) * yarray [highIndex];
836 }
837 Graphics_polyline (g, ixmax + 2 - startOfDefinedStretch, & xarray [startOfDefinedStretch], & yarray [startOfDefinedStretch]);
838 }
839 } catch (MelderError) {
840 Melder_clearError ();
841 }
842 }
843
844 /* End of file Sampled.cpp */
845