Coverage for src/thunderfish/waveanalysis.py: 34%
860 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 20:05 +0000
1"""
2Analysis of wave-type EODs.
4## Analysis
6- `extract_wave()`: retrieve average EOD waveform via Fourier transform.
7- `condition_wave()`: subtract offset, flip, and shift wave-type EOD waveform.
8- `analyze_wave_properties()`: characterize basic properties of a wave-type EOD.
9- `analyze_wave_phases()`: characterize all phases of a wave-type EOD.
10- `analyze_wave_spectrum()`: analyze the spectrum of a wave-type EOD.
12### Complete analysis
14Calls all the functions listed above:
16- `analyze_wave()`: full analysis the EOD waveform of a wave fish.
18## Quality assessment
20- `wave_quality()`: asses quality of EOD waveform of a wave fish.
22## Visualization
24- `plot_wave_eod()`: plot and annotate a wave-type EOD waveform.
25- `plot_wave_spectrum()`: plot and annotate spectrum of wave-type EODs.
27## Storage
29- `save_wave_eodfs()`: save frequencies of wave EODs to file.
30- `load_wave_eodfs()`: load frequencies of wave EODs from file.
31- `save_wave_fish()`: save properties of wave EODs to file.
32- `load_wave_fish()`: load properties of wave EODs from file.
33- `save_wave_phases()`: save phase properties of wave-type EOD to file.
34- `load_wave_phases()`: load phase properties of wave-type EOD from file.
35- `save_wave_spectrum()`: save amplitude and phase spectrum of wave EOD to file.
36- `load_wave_spectrum()`: load amplitude and phase spectrum of wave EOD from file.
38## Configuration
40- `add_extract_wave_config()`: add parameters for `extract_wave()` to configuration.
41- `extract_wave_args()`: retrieve parameters for `extract_wave()` from configuration.
42- `add_analyze_wave_config()`: add parameters for `analyze_wave()` to configuration.
43- `analyze_wave_args()`: retrieve parameters for `analyze_wave()` from configuration.
44- `wave_quality_args()`: retrieve parameters for `wave_quality()` from configuration.
46"""
48import numpy as np
50try:
51 import matplotlib.pyplot as plt
52 from matplotlib.ticker import MultipleLocator, FuncFormatter
53 from matplotlib.artist import setp
54except ImportError:
55 pass
57from pathlib import Path
58from scipy.stats import linregress
59from numba import jit
60from thunderlab.eventdetection import detect_peaks
61from thunderlab.eventdetection import threshold_crossings, threshold_crossing_times, merge_events
62from thunderlab.fourier import fourier_coeffs, normalize_fourier_coeffs
63from thunderlab.fourier import fourier_synthesis
64from thunderlab.powerspectrum import decibel
65from thunderlab.tabledata import TableData
67from .files import open_from_zip, close_zip
68from .harmonics import fundamental_freqs_and_power
71def extract_wave(data, rate, freq, deltaf,
72 win_fac='auto', min_segments=6, periods=5,
73 frate=1e6, max_harmonics=20,
74 min_corr=0.99, min_ampl_frac=0.5,
75 verbose=0, plot_level=0):
76 """Retrieve average EOD waveform via Fourier decomposition.
78 Fourier series are extracted for \\(n\\) frequencies
79 from `freq` \\(\\pm\\) `deltaf` with `n = 1 + 4*win_fac`
80 in multiple data segments of `1/(deltaf*win_fac)` duration.
82 For each data segment, the frequency for which the corresponding
83 Fourier decomposition results in the largest peak-peak amplitude
84 is taken.
86 Then a correlation matrix between the waveform estimates of each
87 data segments is computed and the segments that are most similar
88 to each other are selected. From these the ones whose standard
89 deviation is smaller than `min_ampl_fac` of the maximum one are
90 also discarded. The returned waveform and Fourier series is the
91 average of the remaining data segments.
93 Parameters
94 ----------
95 data: 1-D array of float
96 The data to be analysed.
97 rate: float
98 Sampling rate of the data in Hertz.
99 freq: float
100 Estimated EOD frequency.
101 deltaf: float
102 Frequency resolution in Hertz. Ususally the resolution of the power spectrum
103 from which `freq` was deduced.
104 win_fac: int or 'auto'
105 Integer factor by which the time window and
106 thus the frequency resolution is reduced.
107 Reducing the time window increases the chance for
108 stable waveforms, but might miss the right frequency.
109 If 'auto' then it is set to `max(1, int(np.round(freq/400)))`,
110 i.e. for higher EOD frequencies thisfactor gets larger.
111 min_segments: int
112 Minimum required data segments that fit into the data with
113 an overlap of three quarters. If necessary, the segment sizes
114 and thus also the frequency resolution is reduced by
115 powers of two.
116 periods: float
117 The duration of the returned waveform estimate is the
118 EOD period (`1/freq`) times `periods`.
119 frate: float
120 Sampling rate used for the returned waveform estimates.
121 max_harmonics: int
122 Highest harmonics used for the Fourier decomposition.
123 min_corr: float
124 Minimum required correlation between two waveform estimates.
125 min_ampl_fac: float
126 Minimum required standard deviation of waveform estimate
127 relative to the largest one.
128 verbose: int
129 Verbosity level.
130 plot_level: int
131 When larger than zero, plot some debug information.
133 Returns
134 -------
135 mean_coeffs: 1-D array of complex
136 The averaged Fourier coefficients of the extracted EOD waveform.
137 They correspond to the returnd `mean_eod`.
138 mean_eod: 2-D array of float
139 Average of the EOD snippets. First column is time in seconds,
140 second column the mean eod, third column the standard error.
141 eod_freq: float
142 Refined EOD frequency.
143 times: 1-D array
144 Start times of data segments in which Fourier coefficients
145 have been extracted.
146 n_eods: int
147 Number of EODs that went into the estimate of the mean waveform.
148 skip_reason: str
149 An empty string if the waveform is good, otherwise a string
150 indicating the failure.
152 """
154 def fourier_freq_range(data, rate, frange, nh, n):
155 """ Extract Fourier coefficients for some frequencies and
156 return waveform with largest p-p amplitude.
157 """
158 wave = np.zeros(1)
159 freq = 0.0
160 twave = np.arange(n)/rate
161 for f in frange:
162 w = np.zeros(len(twave))
163 t = np.arange(len(data))/rate
164 for k in range(nh):
165 Xk = np.trapz(data*np.exp(-2j*np.pi*k*f*t), t)*2/t[-1]
166 w += np.real(Xk*np.exp(2j*np.pi*k*f*twave))
167 if np.max(w) - np.min(w) > np.max(wave) - np.min(wave):
168 wave = w
169 freq = f
170 return wave, freq
172 # reduce frequency resolution and time window for high frequency fish:
173 if win_fac == 'auto':
174 win_fac = max(1, int(np.round(freq/400)))
175 deltaf *= win_fac
176 t_segment = min(len(data)/rate, 1/deltaf)
177 nfreqs = 1 + 2*2*win_fac # twice the frequency resolution is necessary and sufficient!
178 step = max(16, int(t_segment*rate))
179 # extract Fourier series from data segments:
180 n = int(periods/freq*rate)
181 freqs = []
182 indices = np.arange(0, max(1, len(data) - step + 1), step//4)
183 while len(indices) < min_segments and step > 8:
184 t_segment /= 2
185 step = max(8, int(t_segment*rate))
186 indices = np.arange(0, max(1, len(data) - step + 1), step//4)
187 times = indices/rate
188 frange = np.linspace(freq - deltaf, freq + deltaf, nfreqs)
189 for i in indices:
190 w, f = fourier_freq_range(data[i:i + step], rate, frange, 6, n)
191 freqs.append(f)
192 freqs = np.array(freqs)
193 mean_coeffs = np.zeros(0, dtype=complex)
194 mean_eod = np.zeros((0, 3))
195 if len(freqs) == 0:
196 # TODO: Why??? How can indices be empty?
197 return mean_coeffs, mean_eod, freq, np.array([]), 0, f'no frequencies detected ({len(indicies)} indices, freqs={freqs})'
198 # refined Fourier series and waveforms:
199 n = int(periods/np.mean(freqs)*frate)
200 coeffs = np.zeros((len(indices), max_harmonics + 1), dtype=complex)
201 waves = np.zeros((len(indices), n))
202 for k in range(len(indices)):
203 i = indices[k]
204 c = fourier_coeffs(data[i:i + step], rate, freqs[k], max_harmonics)
205 c = normalize_fourier_coeffs(c)
206 w = fourier_synthesis(freqs[k], c, frate, n)
207 coeffs[k] = c
208 waves[k] = w
209 if plot_level > 0:
210 fig, axs = plt.subplots(2, 3, width_ratios=[8, 8, 1],
211 layout='constrained')
212 axs[1, 2].set_visible(False)
213 axs[0, 0].set_title(f'EODf={np.mean(freqs):.2f}Hz')
214 t = np.arange(waves.shape[1])*1000/frate
215 for w in waves:
216 axs[0, 0].plot(t, w)
217 axs[0, 0].set_xlabel('time [ms]')
218 # only snippets that are most similar:
219 if len(waves) <= 1:
220 eodf = np.mean(freqs) if len(freqs) > 0 else freq
221 print(f'extract {freq:7.2f}Hz wave fish: {len(waves)} segments, EODf={eodf:.2f}Hz')
222 # TODO: what to do with single segment?
223 else:
224 corr = np.corrcoef(waves)
225 np.fill_diagonal(corr, 0.0)
226 corr_vals = np.sort(corr[corr > min_corr])
227 if len(corr_vals) == 0:
228 if plot_level > 0:
229 plt.show()
230 return mean_coeffs, mean_eod, freq, times, 0, f'waveforms not stable (max_corr={np.max(corr):.4f} SMALLER than {min_corr:.4f})'
231 # higher correlation threshold:
232 min_c = corr_vals[len(corr_vals)//2]
233 # number of high correlations for each segment:
234 num_c = np.sum(corr > min_c, axis=1)
235 num_cmax = np.max(num_c)
236 num_cthresh = max(1, num_cmax//2)
237 # collect pairs with high correlations,
238 # to stay within a single cluster in the correlation matrix.
239 mask = []
240 sub_corr = np.array(corr)
241 # pair with maximum correlation:
242 while not np.all(np.isnan(sub_corr)):
243 i_max_corr = np.nanargmax(sub_corr)
244 max_corr_a = i_max_corr//len(corr)
245 max_corr_b = i_max_corr%len(corr)
246 if corr[max_corr_a, max_corr_b] < max(min_c, 1 - 0.25*(1 - min_corr)):
247 break
248 sub_corr[:, max_corr_a] = np.nan
249 sub_corr[max_corr_a, :] = np.nan
250 sub_corr[:, max_corr_b] = np.nan
251 sub_corr[max_corr_b, :] = np.nan
252 if num_c[max_corr_a] > num_cthresh:
253 mask.append(max_corr_a)
254 if num_c[max_corr_b] > num_cthresh:
255 mask.append(max_corr_b)
256 for k in range(len(corr)):
257 if k >= len(mask):
258 break
259 corr_col = corr[:, mask[k]]
260 idx = np.argsort(corr_col)[::-1]
261 for i in idx:
262 if corr_col[i] > min_c and num_c[i] > num_cthresh and \
263 i not in mask:
264 mask.append(i)
265 break
266 # if only one was selected, add pair:
267 if len(mask) == 1:
268 mask.append(np.argmax(corr[:, mask[0]]))
269 if len(mask) > 0:
270 break
272 waves = waves[mask]
273 coeffs = coeffs[mask]
274 freqs = freqs[mask]
275 times = times[mask]
276 if plot_level > 0:
277 axs[0, 1].set_title('Correlations')
278 m = axs[0, 1].pcolormesh(corr, vmin=min_corr, vmax=1)
279 axs[0, 1].plot(i_max_corr%len(corr) + 0.5,
280 i_max_corr//len(corr) + 0.5, 'ok')
281 axs[0, 1].plot(i_max_corr//len(corr) + 0.5,
282 i_max_corr%len(corr) + 0.5, 'ok')
283 axs[0, 1].set_aspect('equal')
284 axs[0, 1].set_xlabel('segment')
285 axs[0, 1].set_ylabel('segment')
286 plt.colorbar(m, axs[0, 2], label='correlation coefficient')
287 selected = np.zeros(len(corr))
288 selected[mask] = 1
289 axs[1, 1].axhline(num_cthresh + 0.5, color='k',
290 label='num_cthresh')
291 axs[1, 1].plot(np.arange(len(num_c)), num_c, '-o',
292 label='num_c')
293 axs[1, 1].plot(np.arange(len(num_c))[mask], num_c[mask], 'o',
294 label='selected')
295 axs[1, 1].plot(np.arange(len(selected)), selected, '-o',
296 label='mask')
297 axs[1, 1].set_ylim(bottom=0)
298 axs[1, 1].set_xlabel('segment')
299 axs[1, 1].set_ylabel('num_c')
300 if len(num_c) < 10:
301 axs[1, 1].xaxis.set_major_locator(MultipleLocator(1))
302 if np.max(num_c) < 10:
303 axs[1, 1].yaxis.set_major_locator(MultipleLocator(1))
304 axs[1, 1].legend()
305 if verbose > 0:
306 eodf = np.mean(freqs) if len(freqs) > 0 else np.nan
307 kwargs = dict(legacy='1.25')
308 while True:
309 try:
310 with np.printoptions(formatter={'float': lambda x: f'{x:.2f}'},
311 linewidth=10000, **kwargs):
312 print(f'extract {freq:7.2f}Hz wave fish: min_corr={min_c:.4f}, max_corr={corr_vals[-1]:.4f}, num_cmax={num_cmax}, segments={len(corr)}, num_selected={len(mask)}, selected={mask}, EODfs={freqs}, EODf={eodf:.2f}Hz')
313 break
314 except KeyError:
315 # old numpy versions do not have the legacy='1.25' option:
316 kwargs = dict()
317 if len(waves) == 0:
318 if plot_level > 0:
319 plt.show()
320 return mean_coeffs, mean_eod, freq, times, 0, f'waveforms not stable (min_corr={min_c:.4f}, max_corr={corr_vals[-1]:.4f}, num_cmax={num_cmax})'
321 # only the largest snippets:
322 ampls = np.std(waves, axis=1)
323 mask = ampls >= min_ampl_frac*np.max(ampls)
324 if verbose > 0 and np.sum(mask) < len(ampls):
325 print(f' removed {len(ampls) - np.sum(mask)} small amplitude segments')
326 waves = waves[mask]
327 coeffs = coeffs[mask]
328 freqs = freqs[mask]
329 times = times[mask]
330 if len(waves) == 0:
331 if plot_level > 0:
332 plt.show()
333 return mean_coeffs, mean_eod, freq, times, 0, 'no large waveform'
334 if plot_level > 0:
335 axs[1, 0].set_title(f'selected EODs: EODf={np.mean(freqs):.2f}Hz')
336 t = np.arange(waves.shape[1])*1000/frate
337 for w in waves:
338 axs[1, 0].plot(t, w)
339 axs[1, 0].set_xlabel('time [ms]')
340 plt.show()
342 eod_freq = np.mean(freqs)
343 mean_coeffs = np.mean(coeffs, 0)
344 mean_coeffs[0] = 0 # no offset
345 mean_eod = np.zeros((n, 3))
346 mean_eod[:, 0] = np.arange(len(mean_eod))/frate
347 mean_eod[:, 1] = np.mean(waves, axis=0)
348 mean_eod[:, 1] = fourier_synthesis(eod_freq, mean_coeffs, frate, n)
349 mean_eod[:, 2] = np.std(waves, axis=0)
350 t_covered = 0
351 t_end = -1
352 for t in times:
353 if t_end < t:
354 t_end = t
355 t_covered += t + t_segment - t_end
356 t_end = t + t_segment
357 n_eods = int(t_covered*eod_freq)
358 return mean_coeffs, mean_eod, eod_freq, times, n_eods, ''
361def condition_wave(eod, ratetime, freq, coeffs=None, flip_wave='none'):
362 """Subtract offset, flip, and shift wave-type EOD waveform.
364 Parameters
365 ----------
366 eod: 1-D or 2-D array
367 The eod waveform to be analyzed. If an 1-D array, then this
368 is the waveform and you need to also pass a time array or
369 sampling rate in `ratetime`. If a 2-D array, then first
370 column is time in seconds, second column the EOD
371 waveform, and the last column - if present - is the waveform
372 obtained from Fourier decomposition. All other columns are optional
373 and are not used.
374 ratetime: None or float or array of float
375 If a 1-D array is passed on to `eod` then either the sampling
376 rate in Hertz or the time array corresponding to `eod`.
377 freq: float
378 The frequency of the EOD.
379 coeffs: None or 1-D array of complex
380 The Fourier coefficients of an EOD waveform.
381 flip_wave: 'auto', 'none', 'flip'
382 - 'auto' flip waveform such that the larger extremum is positive.
383 - 'flip' flip waveform.
384 - 'none' do not flip waveform.
386 Returns
387 -------
388 eod: 1-D or 2-D array
389 Same shape as the input `eod`.
390 If no `coeffs` are provided, the mean over integer multiples of
391 the period of the waveform is subtracted.
392 If waveform is flipped, it is flipped in the second and last column.
393 If a last column is present, the waveform in the last column
394 is recomputed from the flipped Fourier coefficients.
395 In the first time column, time zero is set to the maximum
396 of the first Fourier component.
397 time: 1-D array
398 If `eod` is an 1-D array, then a time array is returned,
399 with zero set to the maximum of the first Fourier component.
400 coeffs: None or 1-D array of complex
401 The Fourier coefficients of an EOD waveform with \\(\\pi\\)
402 added to the phases if waveform has been flipped.
403 The Fourier coefficients are then normalized such
404 that the phase of the fundamental is zero.
405 This is only returned, if `coeffs` is not None.
406 flipped: bool
407 True if waveform was flipped.
408 """
409 if eod.ndim == 2:
410 time = eod[:, 0]
411 eodw = eod[:, -1] if eod.shape[1] > 2 else eod[:, 1]
412 else:
413 eodw = eod
414 if isinstance(ratetime, (list, tuple, np.ndarray)):
415 time = ratetime
416 else:
417 time = np.arange(len(eod))/ratetime
419 # subtract mean:
420 if coeffs is None:
421 period = 1/freq
422 p_time = np.floor((time[-1] - time[0])/period)*periods
423 mask = (time >= time[-1] - p_time)
424 if eod.ndim == 2:
425 eod[:, 1] -= np.mean(eod[mask, 1])
426 else:
427 eod -= np.mean(eod[mask])
429 # flip:
430 flipped = False
431 if flip_wave.lower() in ['flip', 'true', 'yes']:
432 flipped = True
433 elif flip_wave.lower() == 'auto':
434 if -np.min(eodw) > np.max(eodw):
435 flipped = True
436 if flipped:
437 if eod.ndim == 2:
438 eod[:, 1] *= -1
439 if eod.shape[1] > 2:
440 eod[:, -1] *= -1
441 else:
442 eod *= -1
443 if coeffs is not None:
444 coeffs *= np.exp(1j*np.pi)
446 # shift time:
447 if coeffs is None:
448 c = fourier_coeffs(eodw, time, freq, 1)
449 else:
450 c = coeffs
451 period = 1/freq
452 phi1 = np.angle(c[1])
453 t_zero = -phi1/(2*np.pi*freq)
454 min_t = period
455 if time[-1] - time[0] < 3*period:
456 min_t = 0.5*period
457 while t_zero < min_t:
458 t_zero += period
459 time -= time[0] + t_zero
460 if coeffs is not None:
461 for k in range(1, len(coeffs)):
462 coeffs[k] *= np.exp(-1j*k*phi1)
463 if eod.ndim == 2:
464 eod[:, 0] = time
466 # return:
467 r = [eod]
468 if eod.ndim == 1:
469 r.append(time)
470 if coeffs is not None:
471 r.append(coeffs)
472 r.append(flipped)
473 return r
476def analyze_wave_properties(eod, ratetime, freq):
477 """Characterize basic properties of a wave-type EOD.
479 Parameters
480 ----------
481 eod: 1-D or 2-D array
482 The eod waveform to be analyzed. If an 1-D array, then this
483 is the waveform and you need to also pass a time array or
484 sampling rate in `ratetime`. If a 2-D array, then first
485 column is time in seconds, second column the EOD
486 waveform, and the last column - if present - is the waveform
487 obtained from Fourier decomposition. Allother columns are optional
488 and are not used.
489 ratetime: None or float or array of float
490 If a 1-D array is passed on to `eod` then either the sampling
491 rate in Hertz or the time array corresponding to `eod`.
492 freq: float
493 The frequency of the EOD.
495 Returns
496 -------
497 pp_ampl: float
498 Peak-to-peak amplitude.
499 rel_max_ampl: float
500 Amplitude of minimum or maximum, whichever is larger, relative to p-p amplitude.
501 distance: float
502 Temporal distance between largest negative trough and positive peak.
503 min_distance: float
504 Temporal distance between largest negative trough and positive peak
505 normalized to first half of the EOD cycle.
506 That is, if `distance` is larger than half an EOD cycle, then the
507 duration of an EOD cycle minus `distance`.
508 rms_sem: None or float
509 Root-mean squared standard deviation of the extracted
510 EOD waveform relative to the p-p amplitude.
511 rms_error: float
512 Root-mean-square difference between Fourier decomposition and
513 EOD waveform relative to the p-p amplitude.
514 """
515 if eod.ndim == 2:
516 time = eod[:, 0]
517 eodw = eod[:, -1] if eod.shape[1] > 2 else eod[:, 1]
518 else:
519 eodw = eod
520 if isinstance(ratetime, (list, tuple, np.ndarray)):
521 time = ratetime
522 else:
523 time = np.arange(len(eod))/ratetime
525 # cut out one period:
526 period = 1/freq
527 mask = (time >= 0) & (time <= period)
528 eodp = eodw[mask]
529 timep = time[mask]
531 # amplitudes:
532 pos_idx = np.argmax(eodp)
533 pos_ampl = abs(eodp[pos_idx])
534 neg_idx = np.argmin(eodp)
535 neg_ampl = abs(eodp[neg_idx])
536 pp_ampl = pos_ampl + neg_ampl
537 max_ampl = max(pos_ampl, neg_ampl)
538 rel_max_ampl = max_ampl/pp_ampl
540 # timing:
541 distance = abs(timep[neg_idx] - timep[pos_idx])
542 min_distance = distance
543 if distance > period/2:
544 min_distance = period - distance
546 # variance and fit error:
547 rms_sem = None
548 if eod.ndim == 2 and eod.shape[1] > 2:
549 rms_sem = np.sqrt(np.mean(eod[mask, 2]**2.0))/pp_ampl
550 rms_error = np.sqrt(np.mean((eod[mask, 1] - eod[mask, -1])**2.0))/pp_ampl
552 return pp_ampl, rel_max_ampl, distance, min_distance, rms_sem, rms_error
555def analyze_wave_phases(eod, ratetime, freq, thresh_frac=0.05):
556 """Characterize all phases of a wave-type EOD.
558 Parameters
559 ----------
560 eod: 1-D or 2-D array
561 The eod waveform to be analyzed. If an 1-D array, then this
562 is the waveform and you need to also pass a time array or
563 sampling rate in `ratetime`. If a 2-D array, then first
564 column is time in seconds, second column the EOD
565 waveform, and the last column - if present - is the waveform
566 obtained from Fourier decomposition. Allother columns are optional
567 and are not used.
568 ratetime: None or float or array of float
569 If a 1-D array is passed on to `eod` then either the sampling
570 rate in Hertz or the time array corresponding to `eod`.
571 freq: float
572 The frequency of the EOD.
573 thresh_frac: float
574 Threshold for detecting peaks and troughs as a fraction of
575 the p-p amplitude.
577 Returns
578 -------
579 phases: dict
580 Dictionary with
582 - "indices": indices of each phase
583 (1 is P1, i.e. the largest positive peak)
584 - "times": times of each phase within an EOD cycle
585 - "amplitudes": amplitudes of each phase
586 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
587 - "widths": widths of each phase computed from zeros
588 - "zeros": time point where amplitude between this and the next phase is half their difference.
590 """
591 if eod.ndim == 2:
592 time = eod[:, 0]
593 eod = eod[:, -1] if eod.shape[1] > 2 else eod[:, 1]
594 elif isinstance(ratetime, (list, tuple, np.ndarray)):
595 time = ratetime
596 else:
597 time = np.arange(len(eod))/ratetime
598 dt = np.mean(np.diff(time))
599 period = 1/freq
601 # threshold:
602 mask = (time >= 0) & (time <= period)
603 pp_ampl = np.max(eod[mask]) - np.min(eod[mask])
604 thresh = thresh_frac*pp_ampl
606 # find peaks and troughs:
607 peak_idx, trough_idx = detect_peaks(eod, thresh)
608 pt_idx = np.sort(np.concatenate((peak_idx, trough_idx)))
610 # maximum peak in first period after zero:
611 pt_pidx = pt_idx[(time[pt_idx] >= -0.25*period) & (time[pt_idx] < 0.75*period)]
612 pi = np.argmax(eod[pt_pidx])
613 max_inx = np.nonzero(pt_idx == pt_pidx[pi])[0][0]
614 max_time = time[pt_idx[max_inx]]
616 # analyse phases:
617 times = []
618 amplitudes = []
619 widths = []
620 zero_times = []
621 p_time = 0
622 sign_fac = 1
623 i = 0
624 for k in range(max_inx - 1, len(pt_idx)):
625 idx = pt_idx[k]
626 if time[idx] - max_time >= period - dt:
627 break
628 n_idx = pt_idx[k + 1]
629 th = 0.5*(eod[idx] + eod[n_idx])
630 snippet = eod[idx:n_idx] - th
631 stimes = time[idx:n_idx]
632 zidx = np.nonzero(snippet[:-1]*snippet[1:] < 0)[0]
633 if len(zidx) == 0:
634 zero_times,append(np.nan)
635 else:
636 zidx = zidx[len(zidx)//2] # reduce to single zero crossing
637 snippet = snippet[zidx:zidx + 2]
638 stimes = stimes[zidx:zidx + 2]
639 if sign_fac > 0:
640 z_time = np.interp(0, snippet[::-1], stimes[::-1])
641 else:
642 z_time = np.interp(0, snippet, stimes)
643 if i > 0:
644 times.append(time[idx])
645 amplitudes.append(eod[idx])
646 widths.append(z_time - p_time)
647 zero_times.append(z_time)
648 p_time = z_time
649 i += 1
650 sign_fac *= -1
651 amplitudes = np.array(amplitudes)
653 # store phase properties:
654 phases = dict(indices=np.arange(len(times)) + 1,
655 times=np.array(times),
656 amplitudes=amplitudes,
657 relamplitudes=amplitudes/amplitudes[0],
658 widths=np.array(widths),
659 zeros=np.array(zero_times))
660 return phases
663def analyze_wave_spectrum(freq, coeffs, n_phase_harmonics=8):
664 """Analyze the spectrum of a wave-type EOD.
666 Parameters
667 ----------
668 freq: float or 2-D array
669 The frequency of the EOD or the list of harmonics (rows) with
670 frequency and peak height (columns) as returned from
671 `harmonics.harmonic_groups()`.
672 coeffs: None or 1-D array of complex
673 The Fourier coefficients of an EOD waveform, first one is the offset.
674 n_phase_harmonics: int
675 Number of harmonics over which to compute the slope of the phases.
677 Returns
678 -------
679 spec: 2-D array of float
680 Amplitudes and phases for each harmonic.
681 Rows are the harmonics, first row is the fundamental frequency
682 with multiplier 1, relative amplitude of one, relative power of 0dB,
683 and phase shift of zero.
684 First six columns are from the spectrum of the extracted
685 waveform, optional column 6 is from the spectrum of the recording:
687 - column 0: multiplier of the harmonics (fundamental is one)
688 - column 1: frequency in Hertz
689 - column 2: amplitude
690 - column 3: amplitude relative to the one of the fundamental
691 - column 4: power of harmonics relative to fundamental in decibel
692 - column 5: phase shift relative to the fundamental
693 - column 6: if `freq` is a list of harmonics, the powers of
694 the harmonics from the power spectrum of the raw data.
696 power: float
697 Total power, i.e. sum of squared Fourier amplitudes.
698 data_power: float or None
699 Total power (sum of data powers) in the data, if available.
700 Only sum over as many harmonics as there are Fourier coefficients.
701 thd: float
702 Total harmonic distortion of the amplitudes \\(a_i\\) of the harmonics
703 is power in the higher harmonics relative to the fundamental:
704 \\[ \\text{thd} = \\sqrt{\\sum_{i=2}^n \\left(\\frac{a_i}{a_1}\\right)^2} \\]
705 max_harmonics: int
706 Harmonics with the largest power. Fundamental is one.
707 db_diff: float
708 Standard deviation of the differences of the decibel powers.
709 A measure of smoothness of the spectrum.
710 phase_slope: float
711 Slope of a linear regression between phases and multipliers of
712 the first `n_phase_harmonics` harmonics.
713 """
714 if hasattr(freq, 'shape'):
715 freq1 = freq[0][0]
716 n = len(coeffs) - 1
717 n += np.sum(freq[:, 0] > (len(coeffs) - 0.5)*freq1)
718 spec = np.zeros((n, 7))
719 spec[:, :] = np.nan
720 k = 0
721 for i in range(len(coeffs) - 1):
722 while k < len(freq) and freq[k, 0] < (i + 0.5)*freq1:
723 k += 1
724 if k >= len(freq):
725 break
726 if freq[k, 0] < (i + 1.5)*freq1:
727 spec[i, 6] = freq[k, 1]
728 k += 1
729 for i in range(len(coeffs) - 1, n):
730 if k >= len(freq):
731 break
732 spec[i, :2] = [np.round(freq[k, 0]/freq1), freq[k, 0]]
733 spec[i, 6] = freq[k, 1]
734 k += 1
735 else:
736 freq1 = freq
737 spec = np.zeros((len(coeffs) - 1, 6))
738 ampl1 = np.abs(coeffs[1])
739 for i in range(1, min(len(coeffs), len(spec) + 1)):
740 ampl = np.abs(coeffs[i])
741 phase = np.angle(coeffs[i])
742 spec[i - 1, :6] = [i, i*freq1, ampl, ampl/ampl1,
743 decibel((ampl/ampl1)**2.0), phase]
745 # harmonics, frequency, amplitude, relative amplitude, power, phase, data power
747 # total power:
748 pnh = len(coeffs) - 1
749 power = decibel(np.sum(spec[:pnh, 2]**2))
750 data_power = None
751 if spec.shape[1] > 6:
752 data_power = decibel(np.nansum(spec[:pnh, 6]))
754 # total harmonic distortion:
755 thd = np.sqrt(np.nansum(spec[1:, 3])**2)
757 # harmonic with maximum power:
758 db_powers = spec[:len(coeffs) - 1, 4]
759 max_harmonics = np.argmax(db_powers) + 1
761 # smoothness of power spectrum:
762 db_diff = np.std(np.diff(db_powers))
764 # slope of unwraped phases:
765 phases = spec[:n_phase_harmonics, 5]
766 phases = np.unwrap(phases)
767 r = linregress(np.arange(len(phases)), phases)
768 phase_slope = r.slope
770 return spec, power, data_power, thd, max_harmonics, db_diff, phase_slope
773def analyze_wave(eod, ratetime, freq, coeffs=None,
774 max_harmonics=20, flip_wave='none', thresh_frac=0.05,
775 n_phase_harmonics=8):
776 """Full analysis of the EOD waveform of a wave fish.
778 Parameters
779 ----------
780 eod: 1-D or 2-D array
781 The eod waveform to be analyzed.
782 If an 1-D array, then this is the waveform and you
783 need to also pass a time array or sampling rate in `ratetime`.
784 If a 2-D array, then first column is time in seconds, second
785 column the EOD waveform, third column, if present, is the
786 standard error of the EOD waveform. Further columns are
787 optional and are not used.
788 ratetime: None or float or array of float
789 If a 1-D array is passed on to `eod` then either the sampling
790 rate in Hertz or the time array corresponding to `eod`.
791 freq: float or 2-D array
792 The frequency of the EOD or the list of harmonics (rows) with
793 frequency and peak height (columns) as returned from
794 `harmonics.harmonic_groups()`.
795 coeffs: None or 1-D array of complex
796 The Fourier coefficients of an EOD waveform.
797 If provided, they are taken for the spectrum and the waveform
798 is updated from them.
799 max_harmonics: int
800 Highest harmonics used for the Fourier decomposition.
801 flip_wave: 'auto', 'none', 'flip'
802 - 'auto' flip waveform such that the larger extremum is positive.
803 - 'flip' flip waveform.
804 - 'none' do not flip waveform.
805 thresh_frac: float
806 Threshold for detecting peaks and troughs as a fraction of
807 the p-p amplitude.
808 n_phase_harmonics: int
809 Number of harmonics over which to compute the slope of the phases.
811 Returns
812 -------
813 meod: 2-D array of floats
814 The eod waveform. First column is time in seconds, second
815 column the eod waveform. Further columns are kept from the
816 input `eod`. And a column is added with the Fourier series.
817 props: dict
818 A dictionary with properties of the analyzed EOD waveform.
820 - type: set to 'wave'.
821 - flipped: true if the waveform was flipped.
822 - EODf: EOD fundamental frequency.
823 - period: period of the EOD, i.e. 1/EODf.
824 - ppampl: peak-to-peak amplitude of the Fourier decomposed EOD waveform.
825 - relmaxampl: amplitude of peak or trough, whichever is larger,
826 relative to p-p amplitude.
827 - power: summed power of all harmonics of the extracted
828 EOD waveform in decibel relative to one.
829 - datapower: summed power of all harmonics of the original
830 data in decibel relative to one.
831 Only if `freq` is a list of harmonics.
832 - noise: root-mean squared standard deviation of the extracted
833 EOD waveform relative to the p-p amplitude.
834 - rmserror: root-mean-square difference between
835 Fourier-decomposition and EOD waveform relative
836 to the p-p amplitude.
837 - ppdist: time between peak and trough relative to EOD period.
838 - minppdist: p-p-distance or EOD period minus p-p-distance,
839 whichever is smaller, relative to EOD period.
840 - peakwidth: width of the maximum peak phase relative
841 to EOD period.
842 - troughwidth: width of minimum trough phase relative
843 to EOD period.
844 - minwidth: peakwidth or troughwidth, whichever is smaller.
845 - phasethresh: threshold used for detecting phases as a
846 factor of the p-p-amplitude.
847 - nphases: number of phases of one EOD cycle.
848 - thd: total harmonic distortion, i.e. square root of sum
849 of squared relative amplitudes of higher harmonics.
850 - maxharmonics: harmonics with maximum amplitude.
851 - power2: power of the second harmonics.
852 - dbdiff: smoothness of power spectrum as standard deviation of
853 differences in decibel power.
854 - phaseslope: slope of a linear regression between phases and
855 multipliers of the harmonics.
857 phases: dict
858 Dictionary with
860 - "indices": indices of each phase
861 (1 is P1, i.e. the largest positive peak)
862 - "times": times of each phase within an EOD cycle
863 - "amplitudes": amplitudes of each phase
864 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
865 - "widths": widths of each phase computed from zeros
866 - "zeros": time point where amplitude between this and the next phase
867 is half their difference.
869 spec: 2-D array of floats
870 spec: 2-D array of float
871 Amplitudes and phases for each harmonic.
872 Rows are the harmonics, first row is the fundamental frequency
873 with multiplier 1, relative amplitude of one, relative power of 0dB,
874 and phase shift of zero.
875 First six columns are from the spectrum of the extracted
876 waveform, optional column 6 is from the spectrum of the recording:
878 - column 0: multiplier of the harmonics (fundamental is one)
879 - column 1: frequency in Hertz
880 - column 2: amplitude
881 - column 3: amplitude relative to the one of the fundamental
882 - column 4: power of harmonics relative to fundamental in decibel
883 - column 5: phase shift relative to the fundamental
884 - column 6: if `freq` is a list of harmonics, the powers of
885 the harmonics from the power spectrum of the raw data.
886 """
887 if eod.ndim == 2:
888 if eod.shape[1] > 2:
889 eeod = eod
890 else:
891 eeod = np.column_stack((eod, np.zeros(len(eod))))
892 rate = 1/np.mean(np.diff(eeod[:, 0]))
893 else:
894 if isinstance(ratetime, (list, tuple, np.ndarray)):
895 time = ratetime
896 rate = 1/np.mean(np.diff(time))
897 else:
898 rate = ratetime
899 time = np.arange(len(eod))/rate
900 eeod = np.zeros((len(eod), 3))
901 eeod[:, 0] = time
902 eeod[:, 1] = eod
903 # storage:
904 meod = np.zeros((eeod.shape[0], eeod.shape[1] + 1))
905 meod[:, :eeod.shape[1]] = eeod
906 meod[:, -1] = np.nan
908 freq1 = freq
909 if hasattr(freq, 'shape'):
910 freq1 = freq[0][0]
911 period = 1/freq1
913 # spectrum:
914 has_spec = coeffs is not None
915 if not has_spec:
916 coeffs = fourier_coeffs(meod[:, 1], meod[:, 0], freq1, max_harmonics)
917 phase1 = np.angle(coeffs[1])
918 deltat = phase1/(2*np.pi*freq1)
919 meod[:, 0] -= deltat # TODO: test direction of shift
921 # update waveform:
922 coeffs = normalize_fourier_coeffs(coeffs)
923 eodw = fourier_synthesis(freq1, coeffs, meod[:, 0])
924 if has_spec:
925 meod[:, 1] = eodw
926 meod[:, -1] = eodw
928 # subtract mean and flip:
929 meod, coeffs, flipped = condition_wave(meod, None, freq1, coeffs, flip_wave)
930 if has_spec:
931 meod[:, 1] = meod[:, -1]
933 # waveform properties:
934 pp_ampl, rel_max_ampl, distance, min_distance, rms_sem, rms_error = \
935 analyze_wave_properties(meod, None, freq1)
937 # phases:
938 phases = analyze_wave_phases(meod, None, freq1, thresh_frac=thresh_frac)
939 p_inx = np.argmax(phases['amplitudes'])
940 peak_width = phases['widths'][p_inx]
941 t_inx = np.argmin(phases['amplitudes'])
942 trough_width = phases['widths'][t_inx]
944 # spectral analysis:
945 spec, power, data_power, thd, max_harmonics, db_diff, phase_slope = \
946 analyze_wave_spectrum(freq, coeffs,
947 n_phase_harmonics=n_phase_harmonics)
949 # store results:
950 props = {}
951 props['type'] = 'wave'
952 props['flipped'] = flipped
953 props['EODf'] = freq1
954 props['period'] = period
955 props['ppampl'] = pp_ampl
956 props['relmaxampl'] = rel_max_ampl
957 props['power'] = power
958 if data_power is not None:
959 props['datapower'] = data_power
960 if rms_sem:
961 props['noise'] = rms_sem
962 props['rmserror'] = rms_error
963 props['ppdist'] = distance/period
964 props['minppdist'] = min_distance/period
965 props['peakwidth'] = peak_width/period
966 props['troughwidth'] = trough_width/period
967 props['minwidth'] = min(peak_width, trough_width)/period
968 props['phasethresh'] = thresh_frac
969 props['nphases'] = len(phases['times'])
970 props['thd'] = thd
971 props['maxharmonics'] = max_harmonics
972 props['power2'] = spec[1, 4]
973 props['dbdiff'] = db_diff
974 props['phaseslope'] = phase_slope
976 return meod, props, phases, spec
979def wave_quality(props, harm_relampl=None, min_freq=0.0,
980 max_freq=2000.0, max_clipped_frac=0.1,
981 max_phases=4, max_rms_sem=0.0, max_rms_error=0.05,
982 min_power=-100.0, max_thd=0.0, max_db_diff=20.0,
983 max_relampl_harm2=0.0, max_relampl_harm3=0.0,
984 max_relampl_harm4=0.0):
985 """Assess the quality of an EOD waveform of a wave fish.
987 Parameters
988 ----------
989 props: dict
990 A dictionary with properties of the analyzed EOD waveform
991 as returned by `analyze_wave()`.
992 harm_relampl: 1-D array of floats or None
993 Relative amplitude of at least the first 3 harmonics without
994 the fundamental.
995 min_freq: float
996 Minimum EOD frequency (`props['EODf']`).
997 max_freq: float
998 Maximum EOD frequency (`props['EODf']`).
999 max_clipped_frac: float
1000 If larger than zero, maximum allowed fraction of clipped data
1001 (`props['clipped']`).
1002 max_phases: int
1003 If larger than zero, maximum number of phases per EOD period
1004 (`props['nphases']`).
1005 max_rms_sem: float
1006 If larger than zero, maximum allowed standard deviation of the
1007 data relative to p-p amplitude (`props['noise']`).
1008 max_rms_error: float
1009 If larger than zero, maximum allowed root-mean-square difference
1010 between EOD waveform and Fourier series relative to p-p amplitude
1011 (`props['rmserror']`).
1012 min_power: float
1013 Minimum power of the EOD in dB (`props['power']`).
1014 max_thd: float
1015 If larger than zero, then maximum total harmonic distortion
1016 (`props['thd']`).
1017 max_db_diff: float
1018 If larger than zero, maximum standard deviation of differences between
1019 logarithmic powers of harmonics in decibel (`props['dbdiff']`).
1020 Low values enforce smoother power spectra.
1021 max_relampl_harm2: float
1022 If larger than zero, maximum allowed amplitude of second harmonic
1023 relative to fundamental (=first harmonics).
1024 max_relampl_harm3: float
1025 If larger than zero, maximum allowed amplitude of third harmonic
1026 relative to fundamental (=first harmonics).
1027 max_relampl_harm4: float
1028 If larger than zero, maximum allowed amplitude of fourth harmonic
1029 relative to fundamental (=first harmonics).
1031 Returns
1032 -------
1033 remove: bool
1034 If True then this is most likely not an electric fish. Remove
1035 it from both the waveform properties and the list of EOD
1036 frequencies. If False, keep it in the list of EOD
1037 frequencies, but remove it from the waveform properties if
1038 `skip_reason` is not empty.
1039 skip_reason: str
1040 An empty string if the waveform is good, otherwise a string
1041 indicating the failure.
1042 msg: str
1043 A textual representation of the values tested.
1044 """
1045 remove = False
1046 msg = []
1047 skip_reason = []
1048 # EOD frequency:
1049 if 'EODf' in props:
1050 eodf = props['EODf']
1051 msg += [f'EODf={eodf:6.1f}Hz']
1052 if eodf < min_freq or eodf > max_freq:
1053 remove = True
1054 skip_reason += [f'invalid EODf={eodf:6.1f}Hz (minimumFrequency={min_freq:6.1f}Hz, maximumFrequency={max_freq:6.1f}Hz)']
1055 # clipped fraction:
1056 if 'clipped' in props:
1057 clipped_frac = props['clipped']
1058 msg += [f'clipped={100*clipped_frac:3.0f}%']
1059 if max_clipped_frac > 0 and clipped_frac >= max_clipped_frac:
1060 skip_reason += [f'clipped={100*clipped_frac:3.0f}% (maximumClippedFraction={100*max_clipped_frac:3.0f}%)']
1061 # too many zero crossings:
1062 if 'nphases' in props:
1063 nphases = props['nphases']
1064 msg += [f'phases={nphases}']
1065 if max_phases > 0 and nphases > max_phases:
1066 skip_reason += [f'too many phases={nphases} (maximumPhases={max_phases})']
1067 # noise:
1068 rms_sem = None
1069 if 'rmssem' in props:
1070 rms_sem = props['rmssem']
1071 if 'noise' in props:
1072 rms_sem = props['noise']
1073 if rms_sem is not None:
1074 msg += [f'rms sem waveform={100*rms_sem:6.2f}%']
1075 if max_rms_sem > 0.0 and rms_sem >= max_rms_sem:
1076 skip_reason += [f'noisy waveform s.e.m.={100*rms_sem:6.2f}% (max {100*max_rms_sem:6.2f}%)']
1077 # fit error:
1078 if 'rmserror' in props:
1079 rms_error = props['rmserror']
1080 msg += [f'rmserror={100*rms_error:6.2f}%']
1081 if max_rms_error > 0.0 and rms_error >= max_rms_error:
1082 skip_reason += [f'noisy rmserror={100*rms_error:6.2f}% (maximumVariance={100*max_rms_error:6.2f}%)']
1083 # wave power:
1084 if 'power' in props:
1085 power = props['power']
1086 msg += [f'power={power:6.1f}dB']
1087 if power < min_power:
1088 skip_reason += [f'small power={power:6.1f}dB (minimumPower={min_power:6.1f}dB)']
1089 # total harmonic distortion:
1090 if 'thd' in props:
1091 thd = props['thd']
1092 msg += [f'thd={100*thd:5.1f}%']
1093 if max_thd > 0.0 and thd > max_thd:
1094 skip_reason += [f'large THD={100*thd:5.1f}% (maxximumTotalHarmonicDistortion={100*max_thd:5.1f}%)']
1095 # smoothness of spectrum:
1096 if 'dbdiff' in props:
1097 db_diff = props['dbdiff']
1098 msg += [f'dBdiff={db_diff:5.1f}dB']
1099 if max_db_diff > 0.0 and db_diff > max_db_diff:
1100 remove = True
1101 skip_reason += [f'not smooth s.d. diff={db_diff:5.1f}dB (maxximumPowerDifference={max_db_diff:5.1f}dB)']
1102 # relative amplitude of harmonics:
1103 if harm_relampl is not None:
1104 for k, max_relampl in enumerate([max_relampl_harm2, max_relampl_harm3, max_relampl_harm4]):
1105 if k >= len(harm_relampl):
1106 break
1107 msg += [f'ampl{k + 2}={100*harm_relampl[k]:5.1f}%']
1108 if max_relampl > 0.0 and k < len(harm_relampl) and harm_relampl[k] >= max_relampl:
1109 num_str = ['Second', 'Third', 'Fourth']
1110 skip_reason += [f'distorted ampl{k + 2}={100*harm_relampl[k]:5.1f}% (maximum{num_str[k]}HarmonicAmplitude={100*max_relampl:5.1f}%)']
1111 return remove, ', '.join(skip_reason), ', '.join(msg)
1114def plot_wave_eod(ax, eod_waveform, props, phases=None,
1115 unit=None, wave_periods=2, rel_width=True,
1116 wave_style=dict(lw=1.5, color='tab:red'),
1117 sem_style=dict(color='0.8'),
1118 phase_style=dict(zorder=0, ls='', marker='o', color='tab:red',
1119 markersize=5, mec='white', mew=1),
1120 zerox_style=dict(zorder=50, ls='', marker='o', color='black',
1121 markersize=5, mec='white', mew=1),
1122 zero_style=dict(lw=0.5, color='0.7'),
1123 fontsize='medium'):
1124 """Plot and annotate a wave-type EOD waveform.
1126 Parameters
1127 ----------
1128 ax: matplotlib axes
1129 Axes used for plotting.
1130 eod_waveform: 2-D array
1131 EOD waveform. First column is time in seconds, second column
1132 the (mean) eod waveform. The optional third column is the
1133 standard error. Further columns are ignored.
1134 props: dict
1135 A dictionary with properties of the analyzed EOD waveform as
1136 returned by `analyze_wave()`.
1137 phases: dict
1138 Dictionary with phase properties as returned by
1139 `analyze_wave_phases()`, `analyze_wave()`, and
1140 `load_wave_phases()`.
1141 unit: str
1142 Optional unit of the data used for y-label.
1143 wave_periods: float
1144 How many periods of a wave EOD are shown.
1145 rel_width: bool
1146 If True annotate width in percent relative to the EOD period,
1147 otherwise report them in ms.
1148 wave_style: dict
1149 Arguments passed on to the plot command for the EOD waveform.
1150 sem_style: dict
1151 Arguments passed on to the fill_between command for the
1152 standard error of the EOD.
1153 phase_style: dict
1154 Arguments passed on to the plot command for marking EOD phases.
1155 zerox_style: dict
1156 Arguments passed on to the plot command for marking zero crossings.
1157 zero_style: dict
1158 Arguments passed on to the plot command for the zero line.
1159 fontsize: str or float or int
1160 Fontsize for annotation text.
1162 """
1163 time = 1000*eod_waveform[:, 0]
1164 eod = eod_waveform[:, 1]
1165 # time axis:
1166 period = 1000/props['EODf']
1167 xlim_l = -0.5*wave_periods*period
1168 xlim_r = +0.5*wave_periods*period + period
1169 xlim = xlim_r - xlim_l
1170 ax.set_xlim(xlim_l, xlim_r)
1171 if xlim < 2:
1172 ax.xaxis.set_major_locator(MultipleLocator(0.5))
1173 elif xlim < 4:
1174 ax.xaxis.set_major_locator(MultipleLocator(1))
1175 elif xlim < 8:
1176 ax.xaxis.set_major_locator(MultipleLocator(2))
1177 ax.set_xlabel('Time [msec]')
1178 # amplitude axis:
1179 ylim = np.max(np.abs(eod[(time >= xlim_l) & (time <= xlim_r)]))
1180 ax.set_ylim(-1.2*ylim, +1.2*ylim)
1181 if unit:
1182 ax.set_ylabel(f'Amplitude [{unit}]')
1183 else:
1184 ax.set_ylabel('Amplitude')
1185 # ax height dimensions:
1186 t = ax.text(0, 0, 'test', fontsize=fontsize)
1187 fs = t.get_fontsize()
1188 t.remove()
1189 pixelx = np.abs(np.diff(ax.get_window_extent().get_points()[:, 0]))[0]
1190 dxu = 2*xlim/pixelx
1191 xfs = fs*dxu
1192 pixely = np.abs(np.diff(ax.get_window_extent().get_points()[:, 1]))[0]
1193 dyu = 2*ylim/pixely
1194 yfs = fs*dyu
1195 # plot zero line:
1196 ax.axhline(0.0, zorder=10, **zero_style)
1197 # plot waveform:
1198 ax.plot(time, eod, zorder=45, **wave_style)
1199 # plot standard error:
1200 if eod_waveform.shape[1] > 2:
1201 std_eod = eod_waveform[:, 2]
1202 ax.fill_between(time, eod + std_eod, eod - std_eod,
1203 zorder=20, **sem_style)
1204 # plot and annotate phases:
1205 if phases is not None and len(phases) > 0:
1206 # mark zero crossings:
1207 zeros = 1000*phases['zeros']
1208 zeros_ampls = np.interp(zeros, time, eod)
1209 ax.plot(zeros, zeros_ampls, **zerox_style)
1210 ax.plot(zeros[-1] - period, zeros_ampls[-1], **zerox_style)
1211 # mark phase peak/trough:
1212 times = 1000*phases['times']
1213 amplitudes = phases['amplitudes']
1214 ax.plot(times, amplitudes, **phase_style)
1215 # annotate period:
1216 ax.plot(times[0] + period, amplitudes[0], **phase_style)
1217 ax.text(times[0] + period + xfs, amplitudes[0], f'T={period:.3g}ms',
1218 va='center', zorder=100, fontsize=fontsize)
1219 # phase peaks and troughs:
1220 widths = phases['widths']
1221 for i in range(len(times)):
1222 index = phases['indices'][i]
1223 ptime = times[i]
1224 if ptime < xlim_l or ptime > xlim_r:
1225 continue
1226 pi = np.argmin(np.abs(time - ptime))
1227 pampl = amplitudes[i]
1228 relampl = phases['relamplitudes'][i]
1229 # text for phase label:
1230 label = f'P{index:.0f}'
1231 if np.abs(ptime) < 1:
1232 ts = f'{1000*ptime:.0f}\u00b5s'
1233 elif np.abs(ptime) < 10:
1234 ts = f'{ptime:.2f}ms'
1235 else:
1236 ts = f'{ptime:.3g}ms'
1237 if index == 1:
1238 label += f'(@ {ts})'
1239 else:
1240 if np.abs(relampl) < 0.05:
1241 ps = f'{100*relampl:.1f}%'
1242 else:
1243 ps = f'{100*relampl:.0f}%'
1244 label += f'({ps} @ {ts})'
1245 # position of phase label:
1246 nampl = amplitudes[i - 1 if i > 0 else i + 1]
1247 local_min = pampl < nampl
1248 if i == 0:
1249 ax.text(ptime + xfs, pampl, label,
1250 ha='left', va='center', zorder=100, fontsize=fontsize)
1251 else:
1252 dy = 1.5*yfs
1253 valign = 'bottom'
1254 if local_min:
1255 dy = -dy
1256 valign = 'top'
1257 ax.text(ptime, pampl + dy, label,
1258 ha='center', va=valign, zorder=100, fontsize=fontsize)
1259 # text for width label:
1260 pwidth = 1000*widths[i]
1261 if rel_width:
1262 ws = f'{100*pwidth/period:.2g}%'
1263 elif pwidth < 1:
1264 ws = f'{1000*pwidth:.0f}\u00b5s'
1265 elif pwidth < 10:
1266 ws = f'{pwidth:.2f}ms'
1267 else:
1268 ws = f'{pwidth:.3g}ms'
1269 if i > 0:
1270 x = 0.5*(zeros[i - 1] + zeros[i])
1271 else:
1272 x = 0.5*(zeros[-1] - period + zeros[i])
1273 y = 0.5*(zeros_ampls[i - 1] + zeros_ampls[i])
1274 if local_min:
1275 if y > amplitudes[i - 1] - yfs:
1276 y = amplitudes[i - 1] - yfs
1277 if y > amplitudes[(i + 1)%len(amplitudes)] - yfs:
1278 y = amplitudes[(i + 1)%len(amplitudes)] - yfs
1279 else:
1280 if y < amplitudes[i - 1] + yfs:
1281 y = amplitudes[i - 1] + yfs
1282 if y < amplitudes[(i + 1)%len(amplitudes)] + yfs:
1283 y = amplitudes[(i + 1)%len(amplitudes)] + yfs
1284 if pwidth < 0.25*period:
1285 ax.text(x, y, ws, ha='center', va='center', rotation='vertical', zorder=100, fontsize=fontsize)
1286 else:
1287 ax.text(x, y, ws, ha='center', va='center', zorder=100, fontsize=fontsize)
1288 # annotate plot:
1289 if props is not None:
1290 label = ''
1291 if 'n' in props:
1292 eods = 'EODs' if props['n'] > 1 else 'EOD'
1293 segs = ''
1294 if 'nsegments' in props:
1295 n_segs = props['nsegments']
1296 segs = f' (in {n_segs} segment{"s" if n_segs > 1 else ""})'
1297 label += f'n = {props["n"]} {eods}{segs}'
1298 if 'flipped' in props and props['flipped']:
1299 if len(label) > 0:
1300 label += '\n'
1301 label += 'flipped'
1302 ax.text(0.03, 1.03, label, transform=ax.transAxes,
1303 va='top', ha='left', zorder=100)
1306def plot_wave_spectrum(axa, axp, spec, props, unit=None,
1307 ampl_style=dict(ls='', marker='o', color='tab:blue', markersize=6),
1308 ampl_stem_style=dict(color='tab:blue', alpha=0.5, lw=2),
1309 phase_style=dict(ls='', marker='p', color='tab:blue', markersize=6),
1310 phase_stem_style=dict(color='tab:blue', alpha=0.5, lw=2),
1311 phase_twopi_style=dict(color='k', ls='--', lw=0.5)):
1312 """Plot and annotate spectrum of wave EOD.
1314 Parameters
1315 ----------
1316 axa: matplotlib axes
1317 Axes for amplitude plot.
1318 axp: matplotlib axes
1319 Axes for phase plot.
1320 spec: 2-D array
1321 The spectrum of the wave-type EOD as returned by
1322 `analyze_wave()`. First column is the harmonics,
1323 second column its frequency, third column its amplitude,
1324 fourth column its amplitude relative to the fundamental, fifth
1325 column is power of harmonics relative to fundamental in
1326 decibel, and sixth column the phase shift relative to the
1327 fundamental.
1328 props: dict
1329 A dictionary with properties of the analyzed EOD waveform as
1330 returned by `analyze_wave()`.
1331 unit: string
1332 Optional unit of the data used for y-label.
1333 ampl_style: dict
1334 Properties of the markers of the amplitude plot.
1335 ampl_stem_style: dict
1336 Properties of the stems of the amplitude plot.
1337 phase_style: dict
1338 Properties of the markers of the phase plot.
1339 phase_stem_style: dict
1340 Properties of the stems of the phase plot.
1341 phase_twopi_style: dict
1342 Properties of the horizontal lines marking \\(2\\pi\\)
1343 in the phase plot.
1344 """
1345 def pi_fmt(x, pos):
1346 i = int(np.round(x/np.pi))
1347 if i == 1:
1348 return '\u03c0'
1349 elif i == 0:
1350 return '0'
1351 elif i == -1:
1352 return '-\u03c0'
1353 else:
1354 return f'{i}\u03c0'
1356 n = min(9, np.sum(np.isfinite(spec[:, 2])))
1357 # amplitudes:
1358 markers, stemlines, _ = axa.stem(spec[:n, 0], spec[:n, 2],
1359 basefmt='none')
1360 setp(markers, clip_on=False, **ampl_style)
1361 setp(stemlines, **ampl_stem_style)
1362 if props and 'thd' in props:
1363 axa.text(1, 1, f'thd={100*props["thd"]:.0f}%',
1364 ha='right', va='top', transform=axa.transAxes)
1365 axa.set_xlim(0.5, n + 0.5)
1366 axa.set_ylim(bottom=0)
1367 axa.xaxis.set_major_locator(MultipleLocator(1))
1368 axa.tick_params('x', direction='out')
1369 if unit:
1370 axa.set_ylabel(f'Amplitude [{unit}]')
1371 else:
1372 axa.set_ylabel('Amplitude')
1373 # phases:
1374 phases = spec[:n, 5]
1375 phases = np.unwrap(phases)
1376 min_p = np.floor(np.min(phases)/np.pi)*np.pi
1377 max_p = np.ceil(np.max(phases)/np.pi)*np.pi
1378 if max_p - min_p < 1.5*np.pi:
1379 max_p += np.pi
1380 for p in np.arange(min_p + np.pi, max_p, np.pi):
1381 if np.abs(p - np.round(p/2/np.pi)*2*np.pi) < 1e-8:
1382 axp.axhline(p, **phase_twopi_style)
1383 markers, stemlines, _ = axp.stem(spec[:n, 0], phases[:n],
1384 basefmt='none')
1385 setp(markers, clip_on=False, **phase_style)
1386 setp(stemlines, **phase_stem_style)
1387 axp.set_xlim(0.5, n + 0.5)
1388 axp.xaxis.set_major_locator(MultipleLocator(1))
1389 axp.tick_params('x', direction='out')
1390 if props and 'phaseslope' in props:
1391 axp.text(0.03, 1, f'slope={props["phaseslope"]:.2g}',
1392 va='top', transform=axp.transAxes)
1393 axp.set_ylim(min_p, max_p)
1394 if max_p - min_p < 3.5*np.pi:
1395 axp.yaxis.set_major_locator(MultipleLocator(1*np.pi))
1396 else:
1397 axp.yaxis.set_major_locator(MultipleLocator(2*np.pi))
1398 axp.yaxis.set_major_formatter(FuncFormatter(pi_fmt))
1399 axp.set_xlabel('Harmonics')
1400 axp.set_ylabel('Phase')
1403def save_wave_eodfs(wave_eodfs, wave_indices, basename, **kwargs):
1404 """Save frequencies of wave EODs to file.
1406 Parameters
1407 ----------
1408 wave_eodfs: list of 2D arrays
1409 Each item is a matrix with the frequencies and powers
1410 (columns) of the fundamental and harmonics (rows) as returned
1411 by `harmonics.harmonic_groups()`.
1412 wave_indices: array
1413 Indices identifying each fish or NaN.
1414 If None no index column is inserted.
1415 basename: string or stream
1416 If string, path and basename of file.
1417 If `basename` does not have an extension,
1418 '-waveeodfs' and a file extension according to `kwargs` are appended.
1419 If stream, write EOD frequencies data into this stream.
1420 kwargs:
1421 Arguments passed on to `TableData.write()`.
1423 Returns
1424 -------
1425 filename: Path
1426 Path and full name of the written file in case of `basename`
1427 being a string. Otherwise, the file name and extension that
1428 would have been appended to a basename.
1430 See Also
1431 --------
1432 load_wave_eodfs()
1434 """
1435 eodfs = fundamental_freqs_and_power(wave_eodfs)
1436 td = TableData()
1437 if wave_indices is not None:
1438 td.append('index', '', '%d',
1439 value=[wi if wi >= 0 else np.nan
1440 for wi in wave_indices])
1441 td.append('EODf', 'Hz', '%7.2f', value=eodfs[:, 0])
1442 td.append('datapower', 'dB', '%7.2f', value=eodfs[:, 1])
1443 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
1444 fp = '-waveeodfs' if not ext else ''
1445 return td.write_file_stream(basename, fp, **kwargs)
1448def load_wave_eodfs(file_handle):
1449 """Load frequencies of wave EODs from file.
1451 Parameters
1452 ----------
1453 file_handle: str, Path, or stream
1454 The file to be loaded.
1455 If str or Path and file does not exist, try to open
1456 possible corresponding zip file and read file from there.
1458 Returns
1459 -------
1460 eodfs: 2D array of floats
1461 EODfs and power of wave type fish.
1462 indices: array of ints
1463 Corresponding indices of fish, can contain negative numbers to
1464 indicate frequencies without fish.
1466 Raises
1467 ------
1468 FileNotFoundError:
1469 If `file_path` does not exist.
1471 See Also
1472 --------
1473 save_wave_eodfs()
1474 """
1475 file_handle, zf = open_from_zip(file_handle)
1476 data = TableData(file_handle)
1477 close_zip(file_handle, zf)
1478 eodfs = data.array()
1479 if 'index' in data:
1480 indices = data[:, 'index']
1481 indices[~np.isfinite(indices)] = -1
1482 indices = np.array(indices, dtype=int)
1483 eodfs = eodfs[:, 1:]
1484 else:
1485 indices = np.zeros(data.rows(), dtype=int) - 1
1486 return eodfs, indices
1489def save_wave_fish(eod_props, unit, basename, **kwargs):
1490 """Save properties of wave EODs to file.
1492 Parameters
1493 ----------
1494 eod_props: list of dict
1495 Properties of EODs as returned by `analyze_wave()` and
1496 `analyze_pulse()`. Only properties of wave fish are saved.
1497 unit: string
1498 Unit of the waveform data.
1499 basename: string or stream
1500 If string, path and basename of file.
1501 If `basename` does not have an extension,
1502 '-wavefish' and a file extension are appended.
1503 If stream, write wave fish properties into this stream.
1504 kwargs:
1505 Arguments passed on to `TableData.write()`.
1507 Returns
1508 -------
1509 filename: Path or None
1510 Path and full name of the written file in case of `basename`
1511 being a string. Otherwise, the file name and extension that
1512 would have been appended to a basename.
1513 None if no wave fish are contained in eod_props and
1514 consequently no file was written.
1516 See Also
1517 --------
1518 load_wave_fish()
1519 """
1520 wave_props = [p for p in eod_props if p['type'] == 'wave']
1521 if len(wave_props) == 0:
1522 return None
1523 td = TableData()
1524 if 'twin' in wave_props[0] or 'samplerate' in wave_props[0] or \
1525 'nfft' in wave_props[0]:
1526 td.append_section('recording')
1527 if 'twin' in wave_props[0]:
1528 td.append('twin', 's', '%7.2f', value=wave_props)
1529 td.append('window', 's', '%7.2f', value=wave_props)
1530 td.append('channel', '', '%d', value=wave_props)
1531 td.append('winclipped', '%', '%.2f',
1532 value=wave_props, fac=100.0)
1533 if 'nsegments' in wave_props[0]:
1534 td.append('nsegments', '', '%d', value=wave_props)
1535 if 'samplerate' in wave_props[0]:
1536 td.append('samplerate', 'kHz', '%.3f',
1537 value=wave_props, fac=0.001)
1538 if 'nfft' in wave_props[0]:
1539 td.append('nfft', '', '%d', value=wave_props)
1540 td.append('dfreq', 'Hz', '%.2f', value=wave_props)
1541 td.append_section('waveform')
1542 td.append('index', '', '%d', value=wave_props)
1543 td.append('n', '', '%5d', value=wave_props)
1544 td.append('flipped', '', '%d', value=wave_props)
1545 td.append('EODf', 'Hz', '%7.2f', value=wave_props)
1546 td.append('period', 'ms', '%7.3f', value=wave_props, fac=1000)
1547 td.append('ppampl', unit, '%.5g', value=wave_props)
1548 td.append('relmaxampl', '%', '%.2f', value=wave_props, fac=100)
1549 td.append('power', 'dB', '%7.2f', value=wave_props)
1550 if 'datapower' in wave_props[0]:
1551 td.append('datapower', 'dB', '%7.2f', value=wave_props)
1552 if 'noise' in wave_props[0]:
1553 td.append('noise', '%', '%.1f', value=wave_props, fac=100)
1554 td.append('rmserror', '%', '%.2f', value=wave_props, fac=100)
1555 if 'clipped' in wave_props[0]:
1556 td.append('clipped', '%', '%.1f', value=wave_props, fac=100)
1557 td.append_section('timing')
1558 td.append('ppdist', '%', '%.2f', value=wave_props, fac=100)
1559 td.append('minppdist', '%', '%.2f', value=wave_props, fac=100)
1560 td.append('peakwidth', '%', '%.2f', value=wave_props, fac=100)
1561 td.append('troughwidth', '%', '%.2f', value=wave_props, fac=100)
1562 td.append('minwidth', '%', '%.2f', value=wave_props, fac=100)
1563 td.append_section('phases')
1564 td.append('phasethresh', '%', '%.1f', value=wave_props, fac=100)
1565 td.append('nphases', '', '%d', value=wave_props)
1566 td.append_section('spectrum')
1567 td.append('thd', '%', '%.2f', value=wave_props, fac=100)
1568 td.append('maxharmonics', '', '%d', value=wave_props)
1569 td.append('power2', 'dB', '%7.2f', value=wave_props)
1570 td.append('dbdiff', 'dB', '%7.2f', value=wave_props)
1571 td.append('phaseslope', '', '%.3f', value=wave_props)
1572 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
1573 fp = '-wavefish' if not ext else ''
1574 return td.write_file_stream(basename, fp, **kwargs)
1577def load_wave_fish(file_handle):
1578 """Load properties of wave EODs from file.
1580 All times are scaled to seconds, all frequencies to Hertz and all
1581 percentages to fractions.
1583 Parameters
1584 ----------
1585 file_handle: str, Path, or stream
1586 The file to be loaded.
1587 If str or Path and file does not exist, try to open
1588 possible corresponding zip file and read file from there.
1590 Returns
1591 -------
1592 eod_props: list of dict
1593 Properties of EODs.
1595 Raises
1596 ------
1597 FileNotFoundError:
1598 If `file_path` does not exist.
1600 See Also
1601 --------
1602 save_wave_fish()
1604 """
1605 file_handle, zf = open_from_zip(file_handle)
1606 data = TableData(file_handle)
1607 close_zip(file_handle, zf)
1608 eod_props = data.dicts()
1609 for props in eod_props:
1610 if 'winclipped' in props:
1611 props['winclipped'] /= 100
1612 if 'nsegments' in props:
1613 props['nsegments'] = int(props['nsegments'])
1614 if 'samplerate' in props:
1615 props['samplerate'] *= 1000
1616 if 'nfft' in props:
1617 props['nfft'] = int(props['nfft'])
1618 props['type'] = 'wave'
1619 props['channel'] = int(props['channel'])
1620 props['index'] = int(props['index'])
1621 props['n'] = int(props['n'])
1622 props['period'] /= 1000
1623 props['relmaxampl'] /= 100
1624 props['ppdist'] /= 100
1625 props['minppdist'] /= 100
1626 props['phasethresh'] /= 100
1627 props['nphases'] = int(props['nphases'])
1628 props['thd'] /= 100
1629 props['noise'] /= 100
1630 props['rmserror'] /= 100
1631 if 'clipped' in props:
1632 props['clipped'] /= 100
1633 props['peakwidth'] /= 100
1634 props['troughwidth'] /= 100
1635 props['minwidth'] /= 100
1636 return eod_props
1639def save_wave_phases(phases, unit, idx, basename, **kwargs):
1640 """Save phase properties of wave-type EOD to file.
1642 Parameters
1643 ----------
1644 phases: dict
1645 Dictionary with
1647 - "indices": indices of each phase
1648 (1 is P1, i.e. the largest positive peak)
1649 - "times": times of each phase within an EOD cycle
1650 - "amplitudes": amplitudes of each phase
1651 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
1652 - "widths": widths of each phase computed from zeros
1653 - "zeros": time point where amplitude between this and the next phase is half the difference.
1655 as returned by `analyze_wave_phases()` and `analyze_wave()`.
1656 unit: string
1657 Unit of the waveform data.
1658 idx: int or None
1659 Index of fish.
1660 basename: string or stream
1661 If string, path and basename of file.
1662 If `basename` does not have an extension,
1663 '-wavephases', the fish index, and a file extension are appended.
1664 If stream, write wave phases into this stream.
1665 kwargs:
1666 Arguments passed on to `TableData.write()`.
1668 Returns
1669 -------
1670 filename: Path
1671 Path and full name of the written file in case of `basename`
1672 being a string. Otherwise, the file name and extension that
1673 would have been appended to a basename.
1675 See Also
1676 --------
1677 load_wave_phases()
1678 """
1679 if not phases or len(phases['times']) == 0:
1680 return None
1681 td = TableData()
1682 td.append('index', '', '%.0f', value=phases['indices'])
1683 td.append('time', 'ms', '%.4f', value=phases['times'], fac=1000)
1684 td.append('amplitude', unit, '%.5f', value=phases['amplitudes'])
1685 td.append('relampl', '%', '%.2f', value=phases['relamplitudes'], fac=100)
1686 td.append('width', 'ms', '%.4f', value=phases['widths'], fac=1000)
1687 td.append('zeros', 'ms', '%.4f', value=phases['zeros'], fac=1000)
1688 fp = ''
1689 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
1690 if not ext:
1691 fp = '-wavephases'
1692 if idx is not None:
1693 fp += f'-{idx}'
1694 return td.write_file_stream(basename, fp, **kwargs)
1697def load_wave_phases(file_handle):
1698 """Load phase properties of wave-type EOD from file.
1700 Parameters
1701 ----------
1702 file_handle: str, Path, or stream
1703 The file to be loaded.
1704 If str or Path and file does not exist, try to open
1705 possible corresponding zip file and read file from there.
1707 Returns
1708 -------
1709 phases: dict
1710 Dictionary with
1712 - "indices": indices of each phase
1713 (1 is P1, i.e. the largest positive peak)
1714 - "times": times of each phase within an EOD cycle
1715 - "amplitudes": amplitudes of each phase
1716 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
1717 - "widths": widths of each phase computed from zeros
1718 - "zeros": time point where amplitude between this and the next phase is half the difference.
1720 unit: string
1721 Unit of phase amplitudes.
1723 Raises
1724 ------
1725 FileNotFoundError:
1726 If `file_path` does not exist.
1728 See Also
1729 --------
1730 save_wave_phases()
1731 """
1732 file_handle, zf = open_from_zip(file_handle)
1733 data = TableData(file_handle)
1734 close_zip(file_handle, zf)
1735 phases = dict(indices=data['index'].astype(int),
1736 times=data['time']*0.001,
1737 amplitudes=data['amplitude'],
1738 relamplitudes=data['relampl']*0.01,
1739 widths=data['width']*0.001,
1740 zeros=data['zeros']*0.001)
1741 return phases, data.unit('amplitude')
1744def save_wave_spectrum(spec_data, unit, idx, basename, **kwargs):
1745 """Save amplitude and phase spectrum of wave EOD to file.
1747 Parameters
1748 ----------
1749 spec_data: 2D array of floats
1750 Amplitude and phase spectrum of wave EOD as returned by
1751 `analyze_wave()`.
1752 unit: string
1753 Unit of the waveform data.
1754 idx: int or None
1755 Index of fish.
1756 basename: string or stream
1757 If string, path and basename of file.
1758 If `basename` does not have an extension,
1759 '-wavespectrum', the fish index, and a file extension are appended.
1760 If stream, write wave spectrum into this stream.
1761 kwargs:
1762 Arguments passed on to `TableData.write()`.
1764 Returns
1765 -------
1766 filename: Path
1767 Path and full name of the written file in case of `basename`
1768 being a string. Otherwise, the file name and extension that
1769 would have been appended to a basename.
1771 See Also
1772 --------
1773 load_wave_spectrum()
1775 """
1776 td = TableData(spec_data[:, :6]*[1, 1, 1, 100, 1, 1],
1777 ['harmonics', 'frequency', 'amplitude',
1778 'relampl', 'relpower', 'phase'],
1779 ['', 'Hz', unit, '%', 'dB', 'rad'],
1780 ['%.0f', '%.2f', '%.6f', '%10.2f', '%6.2f', '%8.4f'])
1781 if spec_data.shape[1] > 6:
1782 td.append('datapower', f'{unit}^2/Hz', '%11.4e',
1783 value=spec_data[:, 6])
1784 fp = ''
1785 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
1786 if not ext:
1787 fp = '-wavespectrum'
1788 if idx is not None:
1789 fp += f'-{idx}'
1790 return td.write_file_stream(basename, fp, **kwargs)
1793def load_wave_spectrum(file_handle):
1794 """Load amplitude and phase spectrum of wave EOD from file.
1796 Parameters
1797 ----------
1798 file_handle: str, Path, or stream
1799 The file to be loaded.
1800 If str or Path and file does not exist, try to open
1801 possible corresponding zip file and read file from there.
1803 Returns
1804 -------
1805 spec: 2D array of floats
1806 Amplitude and phase spectrum of wave EOD:
1807 harmonics, frequency, amplitude, relative amplitude in dB,
1808 relative power in dB, phase, data power in unit squared.
1809 Can contain NaNs.
1810 unit: string
1811 Unit of amplitudes.
1813 Raises
1814 ------
1815 FileNotFoundError:
1816 If `file_path` does not exist.
1818 See Also
1819 --------
1820 save_wave_spectrum()
1821 """
1822 file_handle, zf = open_from_zip(file_handle)
1823 data = TableData(file_handle)
1824 close_zip(file_handle, zf)
1825 spec = data.array()
1826 spec[:, 3] *= 0.01
1827 return spec, data.unit('amplitude')
1830def add_extract_wave_config(cfg, win_fac='auto', min_segments=6,
1831 periods=5, frate=1e6, max_harmonics=20,
1832 min_corr=0.99, min_ampl_frac=0.5):
1833 """Add all parameters needed for `extract_wave()` as a new
1834 section to a configuration.
1836 Parameters
1837 ----------
1838 cfg: ConfigFile
1839 The configuration.
1841 See `extract_wave()` for details on the remaining arguments.
1842 """
1843 cfg.add_section('Wave-type EOD extraction:')
1844 cfg.add('dataSegmentDivider', win_fac, '', 'Reduce duration of data segments for EOD waveform estimation by this factor.')
1845 cfg.add('minimumDataSegments', min_segments, '', 'Minimum number of data segments required for EOD waveform estimation.')
1846 cfg.add('periodsWaveEOD', periods, '', 'Number of periods computed for an EOD estimate.')
1847 cfg.add('samplingRateWaveEOD', 0.001*frate, 'kHz', 'Sampling rate of EOD waveform estimates.')
1848 cfg.add('eodHarmonics', max_harmonics, '', 'Highest harmonics of the Fourier decomposition.')
1849 cfg.add('minimumEODCorrelations', min_corr, '', 'Minimum correlation between good EOD waveform estimates.')
1850 cfg.add('minimumEODAmplitude', 100*min_ampl_frac, '%', 'Minimum amplitude of an EOD estimate relative to the largest one.')
1853def extract_wave_args(cfg):
1854 """Retrieve parameters for `extract_wave()` from configuration.
1856 The return value can then be passed as key-word arguments to this
1857 function.
1859 Parameters
1860 ----------
1861 cfg: ConfigFile
1862 The configuration.
1864 Returns
1865 -------
1866 a: dict
1867 Dictionary with names of arguments of the `extract_wave()` function
1868 and their values as supplied by `cfg`.
1869 """
1870 a = cfg.map(win_fac='dataSegmentDivider',
1871 min_segments='minimumDataSegments',
1872 periods='periodsWaveEOD',
1873 frate='samplingRateWaveEOD',
1874 max_harmonics='eodHarmonics',
1875 min_corr='minimumEODCorrelations',
1876 min_ampl_frac='minimumEODAmplitude')
1877 if a['win_fac'] != 'auto':
1878 a['win_fac'] = int(a['win_fac'])
1879 a['frate'] *= 1000
1880 a['min_ampl_frac'] *= 0.01
1881 return a
1884def add_analyze_wave_config(cfg, max_harmonics=20, flip_wave='none',
1885 thresh_frac=0.05, n_phase_harmonics=8):
1886 """Add all parameters needed for `analyze_wave()` as a new
1887 section to a configuration.
1889 Parameters
1890 ----------
1891 cfg: ConfigFile
1892 The configuration.
1894 See `analyze_wave()` for details on the remaining arguments.
1895 """
1896 cfg.add_section('Wave-type EOD analysis:')
1897 cfg.add('eodHarmonics', max_harmonics, '', 'Highest harmonics of the Fourier decomposition.')
1898 cfg.add('flipWaveEOD', flip_wave, '', 'Flip EOD of wave fish to make largest extremum positive (flip, none, or auto).')
1899 cfg.add('waveEODThresholdFraction', 100*thresh_frac, '%', 'Threshold for detecting peaks and troughs in wave-type EODs as a fraction of the p-p amplitude.')
1900 cfg.add('waveEODPhaseHarmonics', n_phase_harmonics, '', 'Number of harmonics over which to compute the slope of the phases..')
1903def analyze_wave_args(cfg):
1904 """Retrieve parameters for `analyze_wave()` from configuration.
1906 The return value can then be passed as key-word arguments to this
1907 function.
1909 Parameters
1910 ----------
1911 cfg: ConfigFile
1912 The configuration.
1914 Returns
1915 -------
1916 a: dict
1917 Dictionary with names of arguments of the `analyze_wave()` function
1918 and their values as supplied by `cfg`.
1919 """
1920 a = cfg.map(max_harmonics='eodHarmonics',
1921 flip_wave='flipWaveEOD',
1922 thresh_frac='waveEODThresholdFraction',
1923 n_phase_harmonics='waveEODPhaseHarmonics')
1924 a['thresh_frac'] *= 0.01
1925 return a
1928def wave_quality_args(cfg):
1929 """Retrieve parameters for `wave_quality()` from configuration.
1931 The return value can then be passed as key-word arguments to this
1932 function.
1934 Parameters
1935 ----------
1936 cfg: ConfigFile
1937 The configuration.
1939 Returns
1940 -------
1941 a: dict
1942 Dictionary with names of arguments of the `wave_quality()` function
1943 and their values as supplied by `cfg`.
1945 See Also
1946 --------
1947 eodanalysis.add_eod_quality_config()
1949 """
1950 a = cfg.map(min_freq='minimumFrequency',
1951 max_freq='maximumFrequency',
1952 max_clipped_frac='maximumClippedFraction',
1953 max_phases='maximumPhases',
1954 max_rms_sem='maximumVariance',
1955 max_rms_error='maximumRMSError',
1956 min_power='minimumPower',
1957 max_thd='maximumTotalHarmonicDistortion',
1958 max_db_diff='maximumPowerDifferences',
1959 max_relampl_harm2='maximumSecondHarmonicAmplitude',
1960 max_relampl_harm3='maximumThirdHarmonicAmplitude',
1961 max_relampl_harm4='maximumFourthHarmonicAmplitude')
1962 a['max_clipped_frac'] *= 0.01
1963 return a
1966def main():
1967 import matplotlib.pyplot as plt
1968 from thunderlab.eventdetection import snippets
1969 from .fakefish import wavefish_eods, export_wavefish
1971 print('Analysis of wave-type EOD waveforms.')
1973 # data:
1974 eodf = 456 # Hz
1975 rate = 96_000
1976 tmax = 5 # s
1977 data = wavefish_eods('Eigenmannia', eodf, rate, tmax, noise_std=0.02)
1978 unit = 'mV'
1979 eod_times = np.arange(0.5/eodf, tmax - 1/eodf, 1/eodf)
1980 eod_idx = (eod_times*rate).astype(int)
1982 # mean EOD:
1983 snips = snippets(data, eod_idx, start=-300, stop=300)
1984 mean_eod = np.mean(snips, 0)
1986 # analyse EOD:
1987 mean_eod, props, power, error_str = \
1988 analyze_wave(mean_eod, rate, eodf)
1990 # write out as python code:
1991 export_wavefish(power, '', 'Eigenmannia spec')
1993 # plot:
1994 fig, axs = plt.subplot_mosaic('wa\nwp', layout='constrained')
1995 plot_wave_eod(axs['w'], mean_eod, props, unit=unit)
1996 axs['w'].set_title(f'wave fish: EODf = {props["EODf"]:.1f} Hz')
1997 plot_wave_spectrum(axs['a'], axs['p'], power, props)
1998 plt.show()
2001if __name__ == '__main__':
2002 main()