Coverage for src/thunderfish/bestwindow.py: 83%
292 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 20:05 +0000
1"""
2Select the best region within a recording with the most stable signal of largest amplitude that is not clipped.
4## Main functions
5- `clip_amplitudes()`: estimated clipping amplitudes from the data.
6- `best_window_indices()`: select start- and end-indices of the best window
7- `best_window_times()`: select start end end-time of the best window
8- `best_window()`: return data of the best window
9- `analysis_window()`: set clipping amplitudes and find analysis window.
11## Configuration
12- `add_clip_config()`: add parameters for `clip_amplitudes()` to configuration.
13- `clip_args()`: retrieve parameters for `clip_amplitudes()` from configuration.
14- `add_best_window_config()`: add parameters for `best_window()` to configuration.
15- `best_window_args()`: retrieve parameters for `best_window*()` from configuration.
17## Visualization
18- `plot_clipping()`: visualization of the algorithm for detecting clipped amplitudes in `clip_amplitudes()`.
19- `plot_best_window()`: visualization of the algorithm used in `best_window_indices()`.
20- `plot_data_window()`: plot the data and the selected analysis window.
21"""
23import numpy as np
24try:
25 import matplotlib.pyplot as plt
26 import matplotlib.ticker as ticker
27except ImportError:
28 pass
30from thunderlab.eventdetection import percentile_threshold, detect_peaks
31from thunderlab.eventdetection import trim_to_peak
32from audioio import unwrap
35def clip_amplitudes(data, win_indices=10000,
36 min_fac=2.0, nbins=20,
37 min_ampl=None, max_ampl=1.0,
38 min_clip=None, max_clip=0,
39 plot_hist_func=None, **kwargs):
40 """Find the amplitudes where the signal clips.
42 Histograms in data segements of win_indices length are analyzed.
43 If the bins at the edges are more than min_fac times as large as
44 the neighboring bins, clipping at the bin's amplitude is assumed.
46 Parameters
47 ----------
48 data: 1-D array
49 The data.
50 win_indices: int
51 Size of the analysis window in indices.
52 min_fac: float
53 If the first or the second bin is at least `min_fac` times
54 as large as the third bin, their upper bin edge is set as min_clip.
55 Likewise for the last and next-to last bin.
56 nbins: int
57 Number of bins used for computing a histogram within `min_ampl` and `max_ampl`.
58 min_ampl: float or None
59 Minimum to be expected amplitude of the data.
60 If set to None, set it to the negative of max_ampl.
61 max_ampl: float
62 Maximum to be expected amplitude of the data
63 min_clip: float or None
64 Minimum amplitude that is not clipped. If zero, estimate from data.
65 If positive, set to this factor of the minimum amplitude of the data.
66 If set to None, set it to the negative of max_clip.
67 max_clip: float
68 Maximum amplitude that is not clipped. If zero, estimate from data.
69 If negativ, set to this factor of the maximum amplitude of the data.
70 plot_hist_func: function
71 Function for visualizing the histograms, is called for every window.
72 `plot_clipping()` is a simple function that can be passed as `plot_hist_func`
73 to quickly visualize what is going on in `clip_amplitudes()`.
75 Signature:
77 `plot_hist_func(data, winx0, winx1, bins, h, min_clip, max_clip,
78 min_ampl, max_ampl, kwargs)`
80 with the arguments:
82 - `data` (array): the full data array.
83 - `winx0` (int): the start index of the current window.
84 - `winx1` (int): the end index of the current window.
85 - `bins` (array): the bin edges of the histogram.
86 - `h` (array): the histogram, plot it with
87 ```
88 plt.bar(bins[:-1], h, width=np.mean(np.diff(bins)))
89 ```
90 - `min_clip` (float): the current value of the minimum clip amplitude.
91 - `max_clip` (float): the current value of the minimum clip amplitude.
92 - `min_ampl` (float): the minimum amplitude of the data.
93 - `max_ampl` (float): the maximum amplitude of the data.
94 - `kwargs` (dict): further user supplied key-word arguments.
96 Returns
97 -------
98 min_clip: float
99 Minimum amplitude that is not clipped.
100 max_clip: float
101 Maximum amplitude that is not clipped.
102 """
103 if min_ampl is None:
104 min_ampl = -max_ampl
105 if min_clip is None:
106 min_clip = -max_clip
107 if min_clip > 0:
108 min_clip = min_clip*min_ampl
109 if max_clip < 0:
110 max_clip = abs(max_clip)*max_ampl
111 if min_clip > 0 and max_clip > 0:
112 return min_clip, max_clip
113 min_clip = min_ampl
114 max_clip = max_ampl
115 bins = np.linspace(min_ampl, max_ampl, nbins, endpoint=True)
116 win_tinxs = np.arange(0, len(data) - win_indices, win_indices)
117 for wtinx in win_tinxs:
118 h, b = np.histogram(data[wtinx:wtinx + win_indices], bins)
119 if h[0] > min_fac * h[2] and b[0] < 0.4*min_ampl:
120 if h[1] > min_fac * h[2] and b[2] > min_clip:
121 min_clip = b[2]
122 elif b[1] > min_clip:
123 min_clip = b[1]
124 if h[-1] > min_fac * h[-3] and b[-1] > 0.4*max_ampl:
125 if h[-2] > min_fac * h[-3] and b[-3] < max_clip:
126 max_clip = b[-3]
127 elif b[-2] < max_clip:
128 max_clip = b[-2]
129 if plot_hist_func:
130 plot_hist_func(data, wtinx, wtinx + win_indices,
131 b, h, min_clip, max_clip,
132 min_ampl, max_ampl, **kwargs)
133 return min_clip, max_clip
136def plot_clipping(data, winx0, winx1, bins,
137 h, min_clip, max_clip, min_ampl, max_ampl):
138 """Visualize the data histograms and the detected clipping amplitudes.
140 Pass this function as the `plot_hist_func` argument to `clip_amplitudes()`.
141 """
142 fig, (ax1, ax2) = plt.subplots(2, 1)
143 ax1.plot(data[winx0:winx1], 'b')
144 ax1.axhline(min_clip, color='r')
145 ax1.axhline(max_clip, color='r')
146 ax1.set_ylim(-1.0, 1.0)
147 ax2.bar(bins[:-1], h, width=np.mean(np.diff(bins)))
148 ax2.axvline(min_clip, color='r')
149 ax2.axvline(max_clip, color='r')
150 ax2.set_xlim(-1.0, 1.0)
151 plt.show()
154def add_clip_config(cfg, min_clip=0.0, max_clip=0.0,
155 window=1.0, min_fac=2.0, nbins=20,
156 min_ampl=-1.0, max_ampl=1.0):
157 """Add parameter needed for `clip_amplitudes()` as a new section to a configuration.
159 Parameters
160 ----------
161 cfg: ConfigFile
162 The configuration.
163 min_clip: float
164 Default minimum clip amplitude.
165 max_clip: float
166 Default maximum clip amplitude.
168 See `clip_amplitudes()` for details on the remaining arguments.
169 """
170 cfg.add_section('Clipping amplitudes:')
171 cfg.add('minClipAmplitude', min_clip, '', 'Minimum amplitude that is not clipped. If zero estimate from data. If positive, set to this factor of the minimum amplitude of the data.')
172 cfg.add('maxClipAmplitude', max_clip, '', 'Maximum amplitude that is not clipped. If zero estimate from data. If negativ, set to this factor of the maximum amplitude of the data.')
173 cfg.add('clipWindow', window, 's', 'Window size for estimating clip amplitudes.')
174 cfg.add('clipBins', nbins, '', 'Number of bins used for constructing histograms of signal amplitudes.')
175 cfg.add('minClipFactor', min_fac, '',
176 'Edge bins of the histogram of clipped signals have to be larger then their neighbors by this factor.')
179def clip_args(cfg, rate):
180 """Translates a configuration to the respective parameter names of `clip_amplitudes()`.
182 The return value can then be passed as key-word arguments to this
183 function.
185 Parameters
186 ----------
187 cfg: ConfigFile
188 The configuration.
189 rate: float
190 Sampling rate of the data.
192 Returns
193 -------
194 a: dict
195 Dictionary with names of arguments of the `clip_amplitudes()` function
196 and their values as supplied by `cfg`.
197 """
198 a = cfg.map(min_fac='minClipFactor',
199 nbins='clipBins',
200 min_clip='minClipAmplitude',
201 max_clip='maxClipAmplitude')
202 if 'clipWindow' in cfg:
203 a['win_indices'] = int(cfg.value('clipWindow')*rate)
204 return a
207def best_window_indices(data, rate, expand=False, win_size=1.,
208 win_shift=0.5, thresh_fac=0.8, percentile=0.1,
209 min_clip=-np.inf, max_clip=np.inf,
210 w_cv_interv=1.0, w_ampl=1.0, w_cv_ampl=1.0,
211 tolerance=0.2, plot_data_func=None, **kwargs):
212 """Find the window within data most suitable for subsequent analysis.
214 First, large peaks and troughs of the data are detected. Peaks and
215 troughs have to be separated in amplitude by at least the value of a
216 dynamic threshold. The threshold is computed in `win_shift` wide
217 windows as `thresh_fac` times the interpercentile range at
218 the `percentile`-th and 100.0-`percentile`-th percentile of the data
219 using the `thunderlab.eventdetection.percentile_threshold()` function.
221 Second, criteria for selecting the best window are computed for each
222 window of width `win_size` shifted by `win_shift` trough the data. The
223 three criteria are:
225 - the mean peak-to-trough amplitude multiplied with the fraction of
226 non clipped peak and trough amplitudes.
227 - the coefficient of variation of the peak-to-trough amplitude.
228 - the coefficient of variation of the inter-peak and inter-trough
229 intervals.
231 Third, a cost function is computed as a weighted sum of the three
232 criteria (the mean amplitude is taken negatively). The respective
233 weights are given by `w_ampl`, `w_cv_ampl`, and `w_cv_interv`.
235 Finally, a threshold is set to the minimum value of the cost
236 function plus tolerance. Then the largest region with the cost
237 function below this threshold is selected as the best window. If
238 `expand` is `False`, then only the single window with smallest cost
239 within the selected largest region is returned.
241 The data used by best window algorithm can be visualized by supplying the
242 function `plot_data_func`. Additional arguments for this function can
243 be supplied via key-word arguments `kwargs`.
245 Parameters
246 ----------
247 data: 1-D array
248 The data to be analyzed.
249 rate: float
250 Sampling rate of the data in Hertz.
251 expand: boolean
252 If `False` return only the single window with the smallest cost.
253 If `True` return the largest window with the cost below the minimum cost
254 plus tolerance.
255 win_size: float
256 Minimum size of the desired best window in seconds.
257 Choose it large enough for the subsequent analysis.
258 win_shift: float
259 Time shift in seconds between windows. Should be smaller or equal to `win_size`.
260 percentile: float
261 `percentile` parameter for the `thunderlab.eventdetection.percentile_threshold()` function
262 used to estimate thresholds for detecting peaks in the data.
263 thresh_fac: float
264 `thresh_fac` parameter for the `thunderlab.eventdetection.percentile_threshold()` function
265 used to estimate thresholds for detecting peaks in the data.
266 min_clip: float
267 Minimum amplitude below which data are clipped.
268 max_clip: float
269 Maximum amplitude above which data are clipped.
270 w_cv_interv: float
271 Weight for the coefficient of variation of the intervals between detected
272 peaks and throughs.
273 w_ampl: float
274 Weight for the mean peak-to-trough amplitude.
275 w_cv_ampl: float
276 Weight for the coefficient of variation of the amplitudes.
277 tolerance: float
278 Added to the minimum cost for expanding the region of the best window.
279 plot_data_func: function
280 Function for plotting the raw data, detected peaks and troughs, the criteria,
281 the cost function and the selected best window.
282 `plot_best_window()` is a simple function that can be passed as the `plot_data_func`
283 parameter to quickly visualize what is going on in selecting the best window.
285 Signature:
287 ````
288 plot_data_func(data, rate, peak_thresh, peak_idx, trough_idx,
289 idx0, idx1, win_start_times, cv_interv, mean_ampl,
290 cv_ampl, clipped_frac, cost_thresh, thresh,
291 valid_wins, **kwargs)
292 ```
294 with the arguments:
296 - `data` (array): raw data.
297 - `rate` (float): sampling rate of the data.
298 - `peak_thresh` (array): thresholds used for detecting peaks and troughs in each data window.
299 - `peak_idx` (array): indices into raw data indicating detected peaks.
300 - `trough_idx` (array): indices into raw data indicating detected troughs.
301 - `idx0` (int): index of the start of the best window.
302 - `idx1` (int): index of the end of the best window.
303 - `win_start_times` (array): times of the analysis windows.
304 - `cv_interv` (array): coefficients of variation of the inter-peak and -trough
305 intervals.
306 - `mean_ampl` (array): mean peak-to-trough amplitudes.
307 - `cv_ampl` (array): coefficients of variation of the peak-to-trough amplitudes.
308 - `clipped_frac` (array): fraction of clipped peaks or troughs.
309 - `cost` (array): cost function.
310 - `cost_thresh` (float): threshold for the cost function.
311 - `valid_wins` (array): boolean array indicating the windows which fulfill
312 all three criteria.
313 - `**kwargs` (dict): further user supplied key-word arguments.
314 kwargs: dict
315 Keyword arguments passed to `plot_data_func`.
317 Returns
318 -------
319 start_index: int
320 Index of the start of the best window.
321 end_index: int
322 Index of the end of the best window.
323 clipped: float.
324 The fraction of clipped peaks or troughs.
326 Raises
327 ------
328 UserWarning
329 - Not enough data for requested `win_size`.
330 - No peaks detected.
331 - No finite amplitudes detected.
332 - No valid interval CV detected.
333 - No valid amplitude CV detected.
334 """
335 # too little data:
336 if len(data)/rate < win_size:
337 raise UserWarning(f'not enough data (data={len(data) / rate:g}s, win={win_size:g}s)')
339 # threshold for peak detection:
340 threshold = percentile_threshold(data, int(win_shift*rate),
341 thresh_fac=thresh_fac,
342 percentile=percentile)
344 # detect large peaks and troughs:
345 peak_idx, trough_idx = detect_peaks(data, threshold)
346 if len(peak_idx) == 0 or len(trough_idx) == 0:
347 raise UserWarning('no peaks or troughs detected')
349 # compute cv of intervals, mean peak amplitude and its cv:
350 invalid_cv = 1000.0
351 win_size_indices = int(win_size*rate)
352 win_start_inxs = np.arange(0, len(data) - win_size_indices,
353 int(0.5*win_shift*rate))
354 if len(win_start_inxs) == 0:
355 win_start_inxs = [0]
356 cv_interv = np.zeros(len(win_start_inxs))
357 mean_ampl = np.zeros(len(win_start_inxs))
358 cv_ampl = np.zeros(len(win_start_inxs))
359 clipped_frac = np.zeros(len(win_start_inxs))
360 for i, wtinx in enumerate(win_start_inxs):
361 # indices of peaks and troughs inside analysis window:
362 pinx = (peak_idx >= wtinx) & (peak_idx <= wtinx + win_size_indices)
363 tinx = (trough_idx >= wtinx) & (trough_idx <= wtinx + win_size_indices)
364 p_idx, t_idx = trim_to_peak(peak_idx[pinx], trough_idx[tinx])
365 # interval statistics:
366 ipis = np.diff(p_idx)
367 itis = np.diff(t_idx)
368 if len(ipis) > 2:
369 cv_interv[i] = 0.5 * (np.std(ipis) / np.mean(ipis) + np.std(itis) / np.mean(itis))
370 # penalize regions without detected peaks:
371 mean_interv = np.mean(ipis)
372 if p_idx[0] - wtinx > mean_interv:
373 cv_interv[i] *= (p_idx[0] - wtinx) / mean_interv
374 if wtinx + win_size_indices - p_idx[-1] > mean_interv:
375 cv_interv[i] *= (wtinx + win_size_indices - p_idx[-1]) / mean_interv
376 else:
377 cv_interv[i] = invalid_cv
378 # statistics of peak-to-trough amplitude:
379 p2t_ampl = data[p_idx] - data[t_idx]
380 if len(p2t_ampl) > 2:
381 mean_ampl[i] = np.mean(p2t_ampl)
382 cv_ampl[i] = np.std(p2t_ampl) / mean_ampl[i]
383 # penalize for clipped peaks:
384 clipped_frac[i] = float(np.sum(data[p_idx] > max_clip) +
385 np.sum(data[t_idx] < min_clip)) / 2.0 / len(p2t_ampl)
386 mean_ampl[i] *= (1.0 - clipped_frac[i]) ** 2.0
387 else:
388 mean_ampl[i] = 0.0
389 cv_ampl[i] = invalid_cv
391 # check:
392 if len(mean_ampl[mean_ampl >= 0.0]) < 0:
393 raise UserWarning('no finite amplitudes detected')
394 if len(cv_interv[cv_interv < invalid_cv]) <= 0:
395 raise UserWarning('no valid interval cv detected')
396 if len(cv_ampl[cv_ampl < invalid_cv]) <= 0:
397 raise UserWarning('no valid amplitude cv detected')
399 # cost function:
400 cost = w_cv_interv * cv_interv + w_cv_ampl * cv_ampl - w_ampl * mean_ampl
401 thresh = np.min(cost) + tolerance
403 # find largest region with low costs:
404 valid_win_idx = np.nonzero(cost <= thresh)[0]
405 cidx0 = valid_win_idx[0] # start of current window
406 cidx1 = cidx0 + 1 # end of current window
407 win_idx0 = cidx0 # start of largest window
408 win_idx1 = cidx1 # end of largest window
409 i = 1
410 while i < len(valid_win_idx): # loop through all valid window positions
411 if valid_win_idx[i] == valid_win_idx[i - 1] + 1:
412 cidx1 = valid_win_idx[i] + 1
413 else:
414 cidx0 = valid_win_idx[i]
415 if cidx1 - cidx0 > win_idx1 - win_idx0: # current window is largest
416 win_idx0 = cidx0
417 win_idx1 = cidx1
418 i += 1
420 # find single best window within the largest region:
421 if not expand:
422 win_idx0 += np.argmin(cost[win_idx0:win_idx1])
423 win_idx1 = win_idx0 + 1
425 # retrive indices of best window for data:
426 idx0 = win_start_inxs[win_idx0]
427 idx1 = win_start_inxs[win_idx1 - 1] + win_size_indices
429 # clipped data?
430 clipped = np.mean(clipped_frac[win_idx0:win_idx1])
432 if plot_data_func:
433 plot_data_func(data, rate, threshold, peak_idx, trough_idx,
434 idx0, idx1, win_start_inxs/rate, cv_interv,
435 mean_ampl, cv_ampl, clipped_frac, cost, thresh,
436 win_idx0, win_idx1, **kwargs)
438 return idx0, idx1, clipped
441def best_window_times(data, rate, expand=False, win_size=1.,
442 win_shift=0.5, thresh_fac=0.8, percentile=0.1,
443 min_clip=-np.inf, max_clip=np.inf,
444 w_cv_interv=1.0, w_ampl=1.0, w_cv_ampl=1.0,
445 tolerance=0.2, plot_data_func=None, **kwargs):
446 """Find the window within data most suitable for subsequent analysis.
448 See `best_window_indices()` for details.
450 Returns
451 -------
452 start_time: float
453 Time of the start of the best window.
454 end_time: float
455 Time of the end of the best window.
456 clipped: float
457 The fraction of clipped peaks or troughs.
458 """
459 start_inx, end_inx, clipped = best_window_indices(data, rate,
460 expand,
461 win_size,
462 win_shift,
463 thresh_fac,
464 percentile,
465 min_clip,
466 max_clip,
467 w_cv_interv,
468 w_ampl,
469 w_cv_ampl,
470 tolerance,
471 plot_data_func,
472 **kwargs)
473 return start_inx/rate, end_inx/rate, clipped
476def best_window(data, rate, expand=False, win_size=1., win_shift=0.5,
477 thresh_fac=0.8, percentile=0.1, min_clip=-np.inf, max_clip=np.inf,
478 w_cv_interv=1.0, w_ampl=1.0, w_cv_ampl=1.0, tolerance=0.2,
479 plot_data_func=None, **kwargs):
480 """Find the window within data most suitable for subsequent analysis.
482 See `best_window_indices()` for details.
484 Returns
485 -------
486 data: array
487 The data of the best window.
488 clipped: float
489 The fraction of clipped peaks or troughs.
490 """
491 start_inx, end_inx, clipped = best_window_indices(data, rate,
492 expand,
493 win_size,
494 win_shift,
495 thresh_fac,
496 percentile,
497 min_clip,
498 max_clip,
499 w_cv_interv,
500 w_ampl,
501 w_cv_ampl,
502 tolerance,
503 plot_data_func,
504 **kwargs)
505 return data[start_inx:end_inx], clipped
508def plot_best_window(data, rate, threshold,
509 peak_idx, trough_idx, idx0, idx1,
510 win_times, cv_interv, mean_ampl, cv_ampl, clipped_frac,
511 cost, thresh, win_idx0, win_idx1, axs):
512 """Visualize the cost function of used for finding the best window for analysis.
514 Pass this function as the `plot_data_func` to the `best_window_*` functions.
516 Parameters
517 ----------
518 See documentation of the `best_window_indices()` functions.
519 """
520 # raw data:
521 time = np.arange(0.0, len(data))/rate
522 axs[0].plot(time, data, 'b', lw=3)
523 if np.mean(clipped_frac[win_idx0:win_idx1]) > 0.01:
524 axs[0].plot(time[idx0:idx1], data[idx0:idx1], color='magenta', lw=3)
525 else:
526 axs[0].plot(time[idx0:idx1], data[idx0:idx1], color='grey', lw=3)
527 axs[0].plot(time[peak_idx], data[peak_idx], 'o', mfc='red', ms=6)
528 axs[0].plot(time[trough_idx], data[trough_idx], 'o', mfc='green', ms=6)
529 axs[0].plot(time, threshold, '#CCCCCC', lw=2)
530 up_lim = np.max(data) * 1.05
531 down_lim = np.min(data) * .95
532 axs[0].set_ylim((down_lim, up_lim))
533 axs[0].set_ylabel('Amplitude [a.u]')
535 # cv of inter-peak intervals:
536 axs[1].plot(win_times[cv_interv < 1000.0], cv_interv[cv_interv < 1000.0],
537 'o', ms=10, color='grey', mew=2., mec='black', alpha=0.6)
538 axs[1].plot(win_times[win_idx0:win_idx1], cv_interv[win_idx0:win_idx1],
539 'o', ms=10, color='red', mew=2., mec='black', alpha=0.6)
540 axs[1].set_ylabel('CV intervals')
541 axs[1].set_ylim(bottom=0.0)
543 # mean amplitude:
544 axs[2].plot(win_times[mean_ampl > 0.0], mean_ampl[mean_ampl > 0.0],
545 'o', ms=10, color='grey', mew=2., mec='black', alpha=0.6)
546 axs[2].plot(win_times[win_idx0:win_idx1], mean_ampl[win_idx0:win_idx1],
547 'o', ms=10, color='red', mew=2., mec='black', alpha=0.6)
548 axs[2].set_ylabel('Mean amplitude [a.u]')
549 axs[2].set_ylim(bottom=0.0)
551 # cv:
552 axs[3].plot(win_times[cv_ampl < 1000.0], cv_ampl[cv_ampl < 1000.0],
553 'o', ms=10, color='grey', mew=2., mec='black', alpha=0.6)
554 axs[3].plot(win_times[win_idx0:win_idx1], cv_ampl[win_idx0:win_idx1],
555 'o', ms=10, color='red', mew=2., mec='black', alpha=0.6)
556 axs[3].set_ylabel('CV amplitude')
557 axs[3].set_ylim(bottom=0.0)
559 # cost:
560 axs[4].plot(win_times[cost < thresh + 10], cost[cost < thresh + 10],
561 'o', ms=10, color='grey', mew=2., mec='black', alpha=0.6)
562 axs[4].plot(win_times[win_idx0:win_idx1], cost[win_idx0:win_idx1],
563 'o', ms=10, color='red', mew=2., mec='black', alpha=0.6)
564 axs[4].axhline(thresh, color='k')
565 axs[4].set_ylabel('Cost')
566 axs[4].set_xlabel('Time [sec]')
569def plot_data_window(ax, data, rate, unit, idx0, idx1, clipped,
570 dstyle=dict(color='tab:blue', lw=1),
571 wstyle=dict(color='tab:red', lw=1)):
572 """Plot the data and mark the analysis window.
574 Parameters
575 ----------
576 ax: matplotlib axes
577 Axes used for plotting.
578 data: 1-D array
579 The full data trace.
580 rate: float
581 Sampling rate of the data in Hertz.
582 unit: string
583 The unit of the data.
584 idx0: int
585 Start index of the best window.
586 idx1: int
587 Stop index of the best window.
588 clipped: float
589 Fraction of clipped peaks.
590 data_color:
591 Color used for plotting the data trace.
592 window_color:
593 Color used for plotting the selected best window.
594 """
595 time = np.arange(len(data))/rate
596 ax.plot(time[:idx0], data[:idx0], **dstyle)
597 ax.plot(time[idx1:], data[idx1:], **dstyle)
598 if idx1 > idx0:
599 ax.plot(time[idx0:idx1], data[idx0:idx1], **wstyle)
600 label = 'analysis\nwindow'
601 if clipped > 0.0:
602 label += f'\n{100.0*clipped:.0f}% clipped'
603 ax.text(time[(idx0+idx1)//2], 0.0, label, ha='center', va='center')
604 ax.set_xlim(time[0], time[-1])
605 ax.set_xlabel('Time [sec]')
606 if len(unit) == 0 or unit == 'a.u.':
607 ax.set_ylabel('Amplitude')
608 else:
609 ax.set_ylabel(f'Amplitude [{unit}]')
610 ax.yaxis.set_major_locator(ticker.MaxNLocator(3))
613def add_best_window_config(cfg, win_pos='best', win_size=1., win_shift=0.5,
614 thresh_fac=0.8, percentile=0.1,
615 min_clip=-np.inf, max_clip=np.inf,
616 w_cv_interv=1.0, w_ampl=1.0, w_cv_ampl=1.0,
617 tolerance=0.2, expand=False):
618 """ Add parameter needed for the `best_window()` functions as a new section to a configuration.
620 Parameters
621 ----------
622 cfg: ConfigFile
623 The configuration.
625 See `best_window_indices()` for details on the remaining arguments.
626 """
627 cfg.add_section('Analysis window:')
628 cfg.add('windowPosition', win_pos, '', 'Position of the analysis window: "beginning", "center", "end", "best", or a time in seconds.')
629 cfg.add('windowSize', win_size, 's', 'Size of the analysis window. This should be much larger than the expected period of the signal. If 0 select the whole time series.')
630 cfg.add('bestWindowShift', win_shift, 's',
631 'For finding the best analysis window, increment for shifting a windows trough the data. Should be larger than the expected period of the signal.')
632 cfg.add('bestWindowThresholdPercentile', percentile, '%',
633 'Percentile for estimating interpercentile range. Should be smaller than the duty cycle of the periodic signal.')
634 cfg.add('bestWindowThresholdFactor', thresh_fac, '',
635 'Threshold for detecting peaks is interperecntile range of the data times this factor.')
636 cfg.add('weightCVInterval', w_cv_interv, '',
637 'Weight factor for the coefficient of variation of the inter-peak and inter-trough intervals.')
638 cfg.add('weightAmplitude', w_ampl, '',
639 'Weight factor for the mean peak-to-trough amplitudes.')
640 cfg.add('weightCVAmplitude', w_cv_ampl, '',
641 'Weight factor for the coefficient of variation of the peak-to-trough amplitude.')
642 cfg.add('bestWindowTolerance', tolerance, '',
643 'Add this to the minimum value of the cost function to get a threshold for selecting the largest best window.')
644 cfg.add('expandBestWindow', expand, '',
645 'Return the largest valid best window. If False return the sole best window. ')
648def best_window_args(cfg):
649 """Translates a configuration to the respective parameter names of the functions `best_window*()`.
651 The return value can then be passed as key-word arguments to these
652 functions.
654 Parameters
655 ----------
656 cfg: ConfigFile
657 The configuration.
659 Returns
660 -------
661 a: dict
662 Dictionary with names of arguments of the `best_window*()` functions
663 and their values as supplied by `cfg`.
664 """
665 return cfg.map(win_size='windowSize',
666 win_shift='bestWindowShift',
667 percentile='bestWindowThresholdPercentile',
668 thresh_fac='bestWindowThresholdFactor',
669 w_cv_interv='weightCVInterval',
670 w_ampl='weightAmplitude',
671 w_cv_ampl='weightCVAmplitude',
672 tolerance='bestWindowTolerance',
673 expand='expandBestWindow')
676def analysis_window(data, rate, ampl_max, win_pos,
677 cfg, plot_level=0):
678 """Set clipping amplitudes and find analysis window.
680 Parameters
681 ----------
682 data: 1-D array
683 The data to be analyzed.
684 rate: float
685 Sampling rate of the data in Hertz.
686 ampl_max: float
687 Maximum value of input range.
688 win_pos: string or float
689 Position of the analysis window: "beginning", "center", "end" or "best".
690 Alternatively the beginning of the analysis window in seconds.
691 cfg: ConfigFile
692 Configuration for clipping and best window.
693 plot_level: int
694 If larger than zero show a plot with the best window cost functions.
696 Returns
697 -------
698 data: 1-D array
699 The data array of the best window
700 idx0: int
701 The start index of the best window in the original data.
702 idx1: int
703 The end index of the best window in the original data.
704 clipped: float
705 Fraction of clipped amplitudes.
706 min_clip: float
707 Minimum amplitude that is not clipped.
708 max_clip: float
709 Maximum amplitude that is not clipped.
710 """
711 found_bestwindow = True
712 clipped = 0
713 min_clip, max_clip = clip_amplitudes(data, max_ampl=ampl_max,
714 **clip_args(cfg, rate))
715 if 'unwrapData' in cfg and cfg.value('unwrapData'):
716 unwrap(data, 1.5, ampl_max)
717 min_clip *= 2
718 max_clip *= 2
719 # window size parameter:
720 bwa = best_window_args(cfg)
721 if 'win_size' in bwa:
722 del bwa['win_size']
723 window_size = cfg.value('windowSize') if 'windowSize' in cfg else 8.0
724 if window_size <= 0.0:
725 window_size = (len(data)-1)/rate
726 # plot cost function:
727 if win_pos == 'best' and plot_level > 0:
728 fig, axs = plt.subplots(5, 1, sharex=True, figsize=(14., 10.),
729 layout='constrained')
730 axs[0].set_title('Best window for analysis')
731 try:
732 idx0, idx1, clipped = \
733 best_window_indices(data, rate,
734 min_clip=min_clip, max_clip=max_clip,
735 win_size=window_size,
736 plot_data_func=plot_best_window,
737 axs=axs, **bwa)
738 plt.show()
739 except UserWarning as e:
740 found_bestwindow = False
741 else:
742 # too little data:
743 n_win = int(window_size*rate)
744 if len(data) < n_win:
745 return data, 0, 0, False, min_clip, max_clip
746 if win_pos == 'best':
747 try:
748 idx0, idx1, clipped = best_window_indices(data, rate,
749 min_clip=min_clip,
750 max_clip=max_clip,
751 win_size=window_size,
752 **bwa)
753 except UserWarning as e:
754 found_bestwindow = False
755 else:
756 if win_pos[:5] == 'begin':
757 idx0 = 0
758 elif win_pos == 'end':
759 idx0 = len(data) - n_win
760 elif win_pos == 'center':
761 idx0 = (len(data) - n_win)//2
762 else:
763 try:
764 if win_pos[-1] == 's':
765 win_pos = win_pos[:-1]
766 t0 = float(win_pos)
767 except ValueError:
768 found_bestwindow = False
769 t0 = 0.0
770 idx0 = int(t0*rate)
771 idx1 = idx0 + n_win
772 if not found_bestwindow or idx1 > len(data):
773 return data, 0, 0, False, min_clip, max_clip
774 data_seg = data[idx0:idx1]
775 # check for clipping:
776 win_shift = cfg.value('bestWindowShift')
777 thresh_fac = cfg.value('bestWindowThresholdFactor')
778 percentile = cfg.value('bestWindowThresholdPercentile')
779 threshold = percentile_threshold(data_seg,
780 int(win_shift*rate),
781 thresh_fac=thresh_fac,
782 percentile=percentile)
783 peak_idx, trough_idx = detect_peaks(data_seg, threshold)
784 p_idx, t_idx = trim_to_peak(peak_idx, trough_idx)
785 if len(p_idx) > 0:
786 p2t_ampl = data_seg[p_idx] - data_seg[t_idx]
787 clipped = float(np.sum(data_seg[p_idx] > max_clip) +
788 np.sum(data_seg[t_idx] < min_clip))/2/len(p2t_ampl)
789 if found_bestwindow:
790 return data[idx0:idx1], idx0, idx1, clipped, min_clip, max_clip
791 else:
792 return data, 0, 0, False, min_clip, max_clip
795def main(data_file=None):
796 title = 'bestwindow'
797 if data_file is None:
798 # generate data:
799 print('generate waveform...')
800 rate = 100000.0
801 time = np.arange(0.0, 1.0, 1.0 / rate)
802 f = 600.0
803 snippets = []
804 amf = 20.0
805 for ampl in [0.2, 0.5, 0.8]:
806 for am_ampl in [0.0, 0.3, 0.9]:
807 data = ampl * np.sin(2.0 * np.pi * f * time) * (1.0 + am_ampl * np.sin(2.0 * np.pi * amf * time))
808 data[data > 1.3] = 1.3
809 data[data < -1.3] = -1.3
810 snippets.extend(data)
811 data = np.asarray(snippets)
812 title = 'test sines'
813 data += 0.01 * np.random.randn(len(data))
814 else:
815 from thunderlab.dataloader import load_data
816 print(f'load {data_file} ...')
817 data, rate, unit, amax = load_data(data_file)
818 data = data[:,0]
819 title = data_file
821 # determine clipping amplitudes:
822 clip_win_size = 0.5
823 min_clip_fac = 2.0
824 min_clip, max_clip = clip_amplitudes(data, int(clip_win_size * rate),
825 min_fac=min_clip_fac)
826 # min_clip, max_clip = clip_amplitudes(data, int(clip_win_size*rate),
827 # min_fac=min_clip_fac,
828 # plot_hist_func=plot_clipping)
830 # setup plots:
831 fig, axs = plt.subplots(5, 1, sharex=True, figsize=(20, 12),
832 layout='constrained')
833 fig.canvas.manager.set_window_title(title)
835 # compute best window:
836 best_window_indices(data, rate, expand=False, win_size=4.0,
837 win_shift=0.5, thresh_fac=0.8, percentile=0.1,
838 min_clip=min_clip, max_clip=max_clip,
839 w_cv_ampl=10.0, tolerance=0.5,
840 plot_data_func=plot_best_window, axs=axs)
842 plt.tight_layout()
843 plt.show()
846if __name__ == '__main__':
847 import sys
848 data_file = sys.argv[1] if len(sys.argv) > 1 else None
849 main(data_file)