Coverage for src / thunderfish / eodanalysis.py: 14%
930 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-02-09 14:25 +0000
1"""
2Analysis of EOD waveforms.
4## EOD waveform analysis
6- `eod_waveform()`: compute an averaged EOD waveform.
7- `adjust_eodf()`: adjust EOD frequencies to a standard temperature.
9## Similarity of EOD waveforms
11- `wave_similarity()`: root-mean squared difference between two wave fish EODs.
12- `pulse_similarity()`: root-mean squared difference between two pulse fish EODs.
13- `load_species_waveforms()`: load template EOD waveforms for species matching.
15## Quality assessment
17- `clipped_fraction()`: compute fraction of clipped EOD waveform snippets.
18- `wave_quality()`: asses quality of EOD waveform of a wave fish.
19- `pulse_quality()`: asses quality of EOD waveform of a pulse fish.
21## Visualization
23- `plot_eod_recording()`: plot a zoomed in range of the recorded trace.
24- `plot_eod_snippets()`: plot a few EOD waveform snippets.
25- `plot_eod_waveform()`: plot and annotate the averaged EOD-waveform with standard error.
27## Storage
29- `save_eod_waveform()`: save mean EOD waveform to file.
30- `load_eod_waveform()`: load EOD waveform from file.
31- `parse_filename()`: parse components of an EOD analysis file name.
32- `save_analysis(): save EOD analysis results to files.
33- `load_analysis()`: load EOD analysis files.
34- `load_recording()`: load recording.
36## Filter functions
38- `unfilter()`: apply inverse low-pass filter on data.
40## Configuration
42- `add_eod_analysis_config()`: add parameters for EOD analysis functions to configuration.
43- `eod_waveform_args()`: retrieve parameters for `eod_waveform()` from configuration.
44- `analyze_wave_args()`: retrieve parameters for `analyze_wave()` from configuration.
45- `analyze_pulse_args()`: retrieve parameters for `analyze_pulse()` from configuration.
46- `add_species_config()`: add parameters needed for assigning EOD waveforms to species.
47- `add_eod_quality_config()`: add parameters needed for assesing the quality of an EOD waveform.
48- `wave_quality_args()`: retrieve parameters for `wave_quality()` from configuration.
49- `pulse_quality_args()`: retrieve parameters for `pulse_quality()` from configuration.
50"""
52import os
53import io
54import zipfile
55import numpy as np
57try:
58 from matplotlib.ticker import MultipleLocator
59except ImportError:
60 pass
62from pathlib import Path
63from audioio import get_str
64from thunderlab.eventdetection import detect_peaks, snippets
65from thunderlab.powerspectrum import decibel
66from thunderlab.tabledata import TableData
67from thunderlab.dataloader import DataLoader
69from .fakefish import normalize_pulsefish, export_pulsefish
70from .fakefish import normalize_wavefish, export_wavefish
71from .waveanalysis import save_wave_eodfs, load_wave_eodfs
72from .waveanalysis import save_wave_fish, load_wave_fish
73from .waveanalysis import save_wave_spectrum, load_wave_spectrum
74from .pulseanalysis import save_pulse_fish, load_pulse_fish
75from .pulseanalysis import save_pulse_spectrum, load_pulse_spectrum
76from .pulseanalysis import save_pulse_phases, load_pulse_phases
77from .pulseanalysis import save_pulse_gaussians, load_pulse_gaussians
78from .pulseanalysis import save_pulse_times, load_pulse_times
81def eod_waveform(data, rate, eod_times, win_fac=2.0, min_win=0.01,
82 min_sem=False, max_eods=None):
83 """Extract data snippets around each EOD, and compute a mean waveform with standard error.
85 Retrieving the EOD waveform of a wave fish works under the following
86 conditions: (i) at a signal-to-noise ratio \\(SNR = P_s/P_n\\),
87 i.e. the power \\(P_s\\) of the EOD of interest relative to the
88 largest other EOD \\(P_n\\), we need to average over at least \\(n >
89 (SNR \\cdot c_s^2)^{-1}\\) snippets to bring the standard error of the
90 averaged EOD waveform down to \\(c_s\\) relative to its
91 amplitude. For a s.e.m. less than 5% ( \\(c_s=0.05\\) ) and an SNR of
92 -10dB (the signal is 10 times smaller than the noise, \\(SNR=0.1\\) ) we
93 get \\(n > 0.00025^{-1} = 4000\\) data snippets - a recording a
94 couple of seconds long. (ii) Very important for wave fish is that
95 they keep their frequency constant. Slight changes in the EOD
96 frequency will corrupt the average waveform. If the period of the
97 waveform changes by \\(c_f=\\Delta T/T\\), then after \\(n =
98 1/c_f\\) periods moved the modified waveform through a whole period.
99 This is in the range of hundreds or thousands waveforms.
101 NOTE: we need to take into account a possible error in the estimate
102 of the EOD period. This will limit the maximum number of snippets to
103 be averaged.
105 If `min_sem` is set, the algorithm checks for a global minimum of
106 the s.e.m. as a function of snippet number. If there is one then
107 the average is computed for this number of snippets, otherwise all
108 snippets are taken from the provided data segment. Note that this
109 check only works for the strongest EOD in a recording. For weaker
110 EOD the s.e.m. always decays with snippet number (empirical
111 observation).
113 TODO: use power spectra to check for changes in EOD frequency!
115 Parameters
116 ----------
117 data: 1-D array of float
118 The data to be analysed.
119 rate: float
120 Sampling rate of the data in Hertz.
121 eod_times: 1-D array of float
122 Array of EOD times in seconds over which the waveform should be
123 averaged.
124 WARNING: The first data point must be at time zero!
125 win_fac: float
126 The snippet size is the EOD period times `win_fac`. The EOD period
127 is determined as the minimum interval between EOD times.
128 min_win: float
129 The minimum size of the snippets in seconds.
130 min_sem: bool
131 If set, check for minimum in s.e.m. to set the maximum numbers
132 of EODs to be used for computing the average waveform.
133 max_eods: int or None
134 Maximum number of EODs to be used for averaging.
135 unfilter_cutoff: float
136 If not zero, the cutoff frequency for an inverse high-pass filter
137 applied to the mean EOD waveform.
139 Returns
140 -------
141 mean_eod: 2-D array
142 Average of the EOD snippets. First column is time in seconds,
143 second column the mean eod, third column the standard error.
144 eod_times: 1-D array
145 Times of EOD peaks in seconds that have been actually used to calculate the
146 averaged EOD waveform.
147 """
148 # indices of EOD times:
149 eod_idx = np.round(eod_times*rate).astype(int)
151 # window size:
152 period = np.min(np.diff(eod_times))
153 win = 0.5*win_fac*period
154 if 2*win < min_win:
155 win = 0.5*min_win
156 win_inx = int(win*rate)
158 # extract snippets:
159 eod_times = eod_times[(eod_idx >= win_inx) & (eod_idx < len(data)-win_inx)]
160 eod_idx = eod_idx[(eod_idx >= win_inx) & (eod_idx < len(data)-win_inx)]
161 if max_eods and max_eods > 0 and len(eod_idx) > max_eods:
162 dn = (len(eod_idx) - max_eods)//2
163 eod_times = eod_times[dn:dn+max_eods]
164 eod_idx = eod_idx[dn:dn+max_eods]
165 eod_snippets = snippets(data, eod_idx, -win_inx, win_inx)
166 if len(eod_snippets) == 0:
167 return np.zeros((0, 3)), eod_times
169 # optimal number of snippets:
170 step = 10
171 if min_sem and len(eod_snippets) > step:
172 sems = [np.mean(np.std(eod_snippets[:k], axis=0, ddof=1)/np.sqrt(k))
173 for k in range(step, len(eod_snippets), step)]
174 idx = np.argmin(sems)
175 # there is a local minimum:
176 if idx > 0 and idx < len(sems)-1:
177 maxn = step*(idx+1)
178 eod_snippets = eod_snippets[:maxn]
179 eod_times = eod_times[:maxn]
181 # mean and std of snippets:
182 mean_eod = np.zeros((len(eod_snippets[0]), 3))
183 mean_eod[:, 1] = np.mean(eod_snippets, axis=0)
184 if len(eod_snippets) > 1:
185 mean_eod[:, 2] = np.std(eod_snippets, axis=0, ddof=1)/np.sqrt(len(eod_snippets))
187 # time axis:
188 mean_eod[:, 0] = (np.arange(len(mean_eod)) - win_inx) / rate
190 return mean_eod, eod_times
193def adjust_eodf(eodf, temp, temp_adjust=25.0, q10=1.62):
194 """Adjust EOD frequencies to a standard temperature using Q10.
196 Parameters
197 ----------
198 eodf: float or ndarray
199 EOD frequencies.
200 temp: float
201 Temperature in degree celsisus at which EOD frequencies in
202 `eodf` were measured.
203 temp_adjust: float
204 Standard temperature in degree celsisus to which EOD
205 frequencies are adjusted.
206 q10: float
207 Q10 value describing temperature dependence of EOD
208 frequencies. The default of 1.62 is from Dunlap, Smith, Yetka
209 (2000) Brain Behav Evol, measured for Apteronotus
210 lepthorhynchus in the lab.
212 Returns
213 -------
214 eodf_corrected: float or array
215 EOD frequencies adjusted to `temp_adjust` using `q10`.
216 """
217 return eodf*q10**((temp_adjust - temp) / 10.0)
220def unfilter(data, rate, cutoff):
221 """Apply inverse high-pass filter on data.
223 Assumes high-pass filter
224 \\[ \\tau \\dot y = -y + \\tau \\dot x \\]
225 has been applied on the original data \\(x\\), where
226 \\(\\tau=(2\\pi f_{cutoff})^{-1}\\) is the time constant of the
227 filter. To recover \\(x\\) the ODE
228 \\[ \\tau \\dot x = y + \\tau \\dot y \\]
229 is applied on the filtered data \\(y\\).
231 Parameters
232 ----------
233 data: ndarray
234 High-pass filtered original data.
235 rate: float
236 Sampling rate of `data` in Hertz.
237 cutoff: float
238 Cutoff frequency \\(f_{cutoff}\\) of the high-pass filter in Hertz.
240 Returns
241 -------
242 data: ndarray
243 Recovered original data.
244 """
245 tau = 0.5/np.pi/cutoff
246 fac = tau*rate
247 data -= np.mean(data)
248 d0 = data[0]
249 x = d0
250 for k in range(len(data)):
251 d1 = data[k]
252 x += (d1 - d0) + d0/fac
253 data[k] = x
254 d0 = d1
255 return data
258def load_species_waveforms(species_file='none'):
259 """Load template EOD waveforms for species matching.
261 Parameters
262 ----------
263 species_file: str
264 Name of file containing species definitions. The content of
265 this file is as follows:
267 - Empty lines and line starting with a hash ('#') are skipped.
268 - A line with the key-word 'wavefish' marks the beginning of the
269 table for wave fish.
270 - A line with the key-word 'pulsefish' marks the beginning of the
271 table for pulse fish.
272 - Each line in a species table has three fields,
273 separated by colons (':'):
275 1. First field is an abbreviation of the species name.
276 2. Second field is the filename of the recording containing the
277 EOD waveform.
278 3. The optional third field is the EOD frequency of the EOD waveform.
280 The EOD frequency is used to normalize the time axis of a
281 wave fish EOD to one EOD period. If it is not specified in
282 the third field, it is taken from the corresponding
283 *-wavespectrum-* file, if present. Otherwise the species is
284 excluded and a warning is issued.
286 Example file content:
287 ``` plain
288 Wavefish
289 Aptero : F_91009L20-eodwaveform-0.csv : 823Hz
290 Eigen : C_91008L01-eodwaveform-0.csv
292 Pulsefish
293 Gymnotus : pulsefish/gymnotus.csv
294 Brachy : H_91009L11-eodwaveform-0.csv
295 ```
297 Returns
298 -------
299 wave_names: list of str
300 List of species names of wave-type fish.
301 wave_eods: list of 2-D arrays
302 List of EOD waveforms of wave-type fish corresponding to
303 `wave_names`. First column of a waveform is time in seconds,
304 second column is the EOD waveform. The waveforms contain
305 exactly one EOD period.
306 pulse_names: list of str
307 List of species names of pulse-type fish.
308 pulse_eods: list of 2-D arrays
309 List of EOD waveforms of pulse-type fish corresponding to
310 `pulse_names`. First column of a waveform is time in seconds,
311 second column is the EOD waveform.
312 """
313 if not Path(species_file).is_file():
314 return [], [], [], []
315 wave_names = []
316 wave_eods = []
317 pulse_names = []
318 pulse_eods = []
319 fish_type = 'wave'
320 with open(species_file, 'r') as sf:
321 for line in sf:
322 line = line.strip()
323 if len(line) == 0 or line[0] == '#':
324 continue
325 if line.lower() == 'wavefish':
326 fish_type = 'wave'
327 elif line.lower() == 'pulsefish':
328 fish_type = 'pulse'
329 else:
330 ls = [s.strip() for s in line.split(':')]
331 if len(ls) < 2:
332 continue
333 name = ls[0]
334 waveform_file = ls[1]
335 eod = TableData(waveform_file).array()
336 eod[:, 0] *= 0.001
337 if fish_type == 'wave':
338 eodf = None
339 if len(ls) > 2:
340 eodf = float(ls[2].replace('Hz', '').strip())
341 else:
342 spectrum_file = waveform_file.replace('eodwaveform', 'wavespectrum')
343 try:
344 wave_spec = TableData(spectrum_file)
345 eodf = wave_spec[0, 1]
346 except FileNotFoundError:
347 pass
348 if eodf is None:
349 print('warning: unknown EOD frequency of %s. Skip.' % name)
350 continue
351 eod[:, 0] *= eodf
352 wave_names.append(name)
353 wave_eods.append(eod[:, :2])
354 elif fish_type == 'pulse':
355 pulse_names.append(name)
356 pulse_eods.append(eod[:, :2])
357 return wave_names, wave_eods, pulse_names, pulse_eods
360def wave_similarity(eod1, eod2, eod1f=1.0, eod2f=1.0):
361 """Root-mean squared difference between two wave fish EODs.
363 Compute the root-mean squared difference between two wave fish
364 EODs over one period. The better sampled signal is down-sampled to
365 the worse sampled one. Amplitude are normalized to peak-to-peak
366 amplitude before computing rms difference. Also compute the rms
367 difference between the one EOD and the other one inverted in
368 amplitude. The smaller of the two rms values is returned.
370 Parameters
371 ----------
372 eod1: 2-D array
373 Time and amplitude of reference EOD.
374 eod2: 2-D array
375 Time and amplitude of EOD that is to be compared to `eod1`.
376 eod1f: float
377 EOD frequency of `eod1` used to transform the time axis of `eod1`
378 to multiples of the EOD period. If already normalized to EOD period,
379 as for example by the `load_species_waveforms()` function, then
380 set the EOD frequency to one (default).
381 eod2f: float
382 EOD frequency of `eod2` used to transform the time axis of `eod2`
383 to multiples of the EOD period. If already normalized to EOD period,
384 as for example by the `load_species_waveforms()` function, then
385 set the EOD frequency to one (default).
387 Returns
388 -------
389 rmse: float
390 Root-mean-squared difference between the two EOD waveforms relative to
391 their standard deviation over one period.
392 """
393 # copy:
394 eod1 = np.array(eod1[:, :2])
395 eod2 = np.array(eod2[:, :2])
396 # scale to multiples of EOD period:
397 eod1[:, 0] *= eod1f
398 eod2[:, 0] *= eod2f
399 # make eod1 the waveform with less samples per period:
400 n1 = int(1.0/(eod1[1,0]-eod1[0,0]))
401 n2 = int(1.0/(eod2[1,0]-eod2[0,0]))
402 if n1 > n2:
403 eod1, eod2 = eod2, eod1
404 n1, n2 = n2, n1
405 # one period around time zero:
406 i0 = np.argmin(np.abs(eod1[:, 0]))
407 k0 = i0-n1//2
408 if k0 < 0:
409 k0 = 0
410 k1 = k0 + n1 + 1
411 if k1 >= len(eod1):
412 k1 = len(eod1)
413 # cut out one period from the reference EOD around maximum:
414 i = k0 + np.argmax(eod1[k0:k1,1])
415 k0 = i-n1//2
416 if k0 < 0:
417 k0 = 0
418 k1 = k0 + n1 + 1
419 if k1 >= len(eod1):
420 k1 = len(eod1)
421 eod1 = eod1[k0:k1,:]
422 # normalize amplitudes of first EOD:
423 eod1[:, 1] -= np.min(eod1[:, 1])
424 eod1[:, 1] /= np.max(eod1[:, 1])
425 sigma = np.std(eod1[:, 1])
426 # set time zero to maximum of second EOD:
427 i0 = np.argmin(np.abs(eod2[:, 0]))
428 k0 = i0-n2//2
429 if k0 < 0:
430 k0 = 0
431 k1 = k0 + n2 + 1
432 if k1 >= len(eod2):
433 k1 = len(eod2)
434 i = k0 + np.argmax(eod2[k0:k1,1])
435 eod2[:, 0] -= eod2[i,0]
436 # interpolate eod2 to the time base of eod1:
437 eod2w = np.interp(eod1[:, 0], eod2[:, 0], eod2[:, 1])
438 # normalize amplitudes of second EOD:
439 eod2w -= np.min(eod2w)
440 eod2w /= np.max(eod2w)
441 # root-mean-square difference:
442 rmse1 = np.sqrt(np.mean((eod1[:, 1] - eod2w)**2))/sigma
443 # root-mean-square difference of the flipped signal:
444 i = k0 + np.argmin(eod2[k0:k1,1])
445 eod2[:, 0] -= eod2[i,0]
446 eod2w = np.interp(eod1[:, 0], eod2[:, 0], -eod2[:, 1])
447 eod2w -= np.min(eod2w)
448 eod2w /= np.max(eod2w)
449 rmse2 = np.sqrt(np.mean((eod1[:, 1] - eod2w)**2))/sigma
450 # take the smaller value:
451 rmse = min(rmse1, rmse2)
452 return rmse
455def pulse_similarity(eod1, eod2, wfac=10.0):
456 """Root-mean squared difference between two pulse fish EODs.
458 Compute the root-mean squared difference between two pulse fish
459 EODs over `wfac` times the distance between minimum and maximum of
460 the wider EOD. The waveforms are normalized to their maxima prior
461 to computing the rms difference. Also compute the rms difference
462 between the one EOD and the other one inverted in amplitude. The
463 smaller of the two rms values is returned.
465 Parameters
466 ----------
467 eod1: 2-D array
468 Time and amplitude of reference EOD.
469 eod2: 2-D array
470 Time and amplitude of EOD that is to be compared to `eod1`.
471 wfac: float
472 Multiply the distance between minimum and maximum by this factor
473 to get the window size over which to compute the rms difference.
475 Returns
476 -------
477 rmse: float
478 Root-mean-squared difference between the two EOD waveforms
479 relative to their standard deviation over the analysis window.
480 """
481 # copy:
482 eod1 = np.array(eod1[:, :2])
483 eod2 = np.array(eod2[:, :2])
484 # width of the pulses:
485 imin1 = np.argmin(eod1[:, 1])
486 imax1 = np.argmax(eod1[:, 1])
487 w1 = np.abs(eod1[imax1,0]-eod1[imin1,0])
488 imin2 = np.argmin(eod2[:, 1])
489 imax2 = np.argmax(eod2[:, 1])
490 w2 = np.abs(eod2[imax2,0]-eod2[imin2,0])
491 w = wfac*max(w1, w2)
492 # cut out signal from the reference EOD:
493 n = int(w/(eod1[1,0]-eod1[0,0]))
494 i0 = imax1-n//2
495 if i0 < 0:
496 i0 = 0
497 i1 = imax1+n//2+1
498 if i1 >= len(eod1):
499 i1 = len(eod1)
500 eod1[:, 0] -= eod1[imax1,0]
501 eod1 = eod1[i0:i1,:]
502 # normalize amplitude of first EOD:
503 eod1[:, 1] /= np.max(eod1[:, 1])
504 sigma = np.std(eod1[:, 1])
505 # interpolate eod2 to the time base of eod1:
506 eod2[:, 0] -= eod2[imax2,0]
507 eod2w = np.interp(eod1[:, 0], eod2[:, 0], eod2[:, 1])
508 # normalize amplitude of second EOD:
509 eod2w /= np.max(eod2w)
510 # root-mean-square difference:
511 rmse1 = np.sqrt(np.mean((eod1[:, 1] - eod2w)**2))/sigma
512 # root-mean-square difference of the flipped signal:
513 eod2[:, 0] -= eod2[imin2,0]
514 eod2w = np.interp(eod1[:, 0], eod2[:, 0], -eod2[:, 1])
515 eod2w /= np.max(eod2w)
516 rmse2 = np.sqrt(np.mean((eod1[:, 1] - eod2w)**2))/sigma
517 # take the smaller value:
518 rmse = min(rmse1, rmse2)
519 return rmse
522def clipped_fraction(data, rate, eod_times, mean_eod,
523 min_clip=-np.inf, max_clip=np.inf):
524 """Compute fraction of clipped EOD waveform snippets.
526 Cut out snippets at each `eod_times` based on time axis of
527 `mean_eod`. Check which fraction of snippets exceeds clipping
528 amplitude `min_clip` and `max_clip`.
530 Parameters
531 ----------
532 data: 1-D array of float
533 The data to be analysed.
534 rate: float
535 Sampling rate of the data in Hertz.
536 eod_times: 1-D array of float
537 Array of EOD times in seconds.
538 mean_eod: 2-D array with time, mean, sem, and fit.
539 Averaged EOD waveform of wave fish. Only the time axis is used
540 to set width of snippets.
541 min_clip: float
542 Minimum amplitude that is not clipped.
543 max_clip: float
544 Maximum amplitude that is not clipped.
546 Returns
547 -------
548 clipped_frac: float
549 Fraction of snippets that are clipped.
550 """
551 # snippets:
552 idx0 = np.argmin(np.abs(mean_eod[:, 0])) # index of time zero
553 w0 = -idx0
554 w1 = len(mean_eod[:, 0]) - idx0
555 eod_idx = np.round(eod_times*rate).astype(int)
556 eod_snippets = snippets(data, eod_idx, w0, w1)
557 # fraction of clipped snippets:
558 clipped_frac = np.sum(np.any((eod_snippets > max_clip) |
559 (eod_snippets < min_clip), axis=1))\
560 / len(eod_snippets)
561 return clipped_frac
564def wave_quality(props, harm_relampl=None, min_freq=0.0,
565 max_freq=2000.0, max_clipped_frac=0.1,
566 max_crossings=4, max_rms_sem=0.0, max_rms_error=0.05,
567 min_power=-100.0, max_thd=0.0, max_db_diff=20.0,
568 max_harmonics_db=-5.0, max_relampl_harm1=0.0,
569 max_relampl_harm2=0.0, max_relampl_harm3=0.0):
570 """Assess the quality of an EOD waveform of a wave fish.
572 Parameters
573 ----------
574 props: dict
575 A dictionary with properties of the analyzed EOD waveform
576 as returned by `analyze_wave()`.
577 harm_relampl: 1-D array of floats or None
578 Relative amplitude of at least the first 3 harmonics without
579 the fundamental.
580 min_freq: float
581 Minimum EOD frequency (`props['EODf']`).
582 max_freq: float
583 Maximum EOD frequency (`props['EODf']`).
584 max_clipped_frac: float
585 If larger than zero, maximum allowed fraction of clipped data
586 (`props['clipped']`).
587 max_crossings: int
588 If larger than zero, maximum number of zero crossings per EOD period
589 (`props['ncrossings']`).
590 max_rms_sem: float
591 If larger than zero, maximum allowed standard error of the
592 data relative to p-p amplitude (`props['noise']`).
593 max_rms_error: float
594 If larger than zero, maximum allowed root-mean-square error
595 between EOD waveform and Fourier fit relative to p-p amplitude
596 (`props['rmserror']`).
597 min_power: float
598 Minimum power of the EOD in dB (`props['power']`).
599 max_thd: float
600 If larger than zero, then maximum total harmonic distortion
601 (`props['thd']`).
602 max_db_diff: float
603 If larger than zero, maximum standard deviation of differences between
604 logarithmic powers of harmonics in decibel (`props['dbdiff']`).
605 Low values enforce smoother power spectra.
606 max_harmonics_db:
607 Maximum power of higher harmonics relative to peak power in
608 decibel (`props['maxdb']`).
609 max_relampl_harm1: float
610 If larger than zero, maximum allowed amplitude of first harmonic
611 relative to fundamental.
612 max_relampl_harm2: float
613 If larger than zero, maximum allowed amplitude of second harmonic
614 relative to fundamental.
615 max_relampl_harm3: float
616 If larger than zero, maximum allowed amplitude of third harmonic
617 relative to fundamental.
619 Returns
620 -------
621 remove: bool
622 If True then this is most likely not an electric fish. Remove
623 it from both the waveform properties and the list of EOD
624 frequencies. If False, keep it in the list of EOD
625 frequencies, but remove it from the waveform properties if
626 `skip_reason` is not empty.
627 skip_reason: str
628 An empty string if the waveform is good, otherwise a string
629 indicating the failure.
630 msg: str
631 A textual representation of the values tested.
632 """
633 remove = False
634 msg = []
635 skip_reason = []
636 # EOD frequency:
637 if 'EODf' in props:
638 eodf = props['EODf']
639 msg += ['EODf=%5.1fHz' % eodf]
640 if eodf < min_freq or eodf > max_freq:
641 remove = True
642 skip_reason += ['invalid EODf=%5.1fHz (minimumFrequency=%5.1fHz, maximumFrequency=%5.1f))' %
643 (eodf, min_freq, max_freq)]
644 # clipped fraction:
645 if 'clipped' in props:
646 clipped_frac = props['clipped']
647 msg += ['clipped=%3.0f%%' % (100.0*clipped_frac)]
648 if max_clipped_frac > 0 and clipped_frac >= max_clipped_frac:
649 skip_reason += ['clipped=%3.0f%% (maximumClippedFraction=%3.0f%%)' %
650 (100.0*clipped_frac, 100.0*max_clipped_frac)]
651 # too many zero crossings:
652 if 'ncrossings' in props:
653 ncrossings = props['ncrossings']
654 msg += ['zero crossings=%d' % ncrossings]
655 if max_crossings > 0 and ncrossings > max_crossings:
656 skip_reason += ['too many zero crossings=%d (maximumCrossings=%d)' %
657 (ncrossings, max_crossings)]
658 # noise:
659 rms_sem = None
660 if 'rmssem' in props:
661 rms_sem = props['rmssem']
662 if 'noise' in props:
663 rms_sem = props['noise']
664 if rms_sem is not None:
665 msg += ['rms sem waveform=%6.2f%%' % (100.0*rms_sem)]
666 if max_rms_sem > 0.0 and rms_sem >= max_rms_sem:
667 skip_reason += ['noisy waveform s.e.m.=%6.2f%% (max %6.2f%%)' %
668 (100.0*rms_sem, 100.0*max_rms_sem)]
669 # fit error:
670 if 'rmserror' in props:
671 rms_error = props['rmserror']
672 msg += ['rmserror=%6.2f%%' % (100.0*rms_error)]
673 if max_rms_error > 0.0 and rms_error >= max_rms_error:
674 skip_reason += ['noisy rmserror=%6.2f%% (maximumVariance=%6.2f%%)' %
675 (100.0*rms_error, 100.0*max_rms_error)]
676 # wave power:
677 if 'power' in props:
678 power = props['power']
679 msg += ['power=%6.1fdB' % power]
680 if power < min_power:
681 skip_reason += ['small power=%6.1fdB (minimumPower=%6.1fdB)' %
682 (power, min_power)]
683 # total harmonic distortion:
684 if 'thd' in props:
685 thd = props['thd']
686 msg += ['thd=%5.1f%%' % (100.0*thd)]
687 if max_thd > 0.0 and thd > max_thd:
688 skip_reason += ['large THD=%5.1f%% (maxximumTotalHarmonicDistortion=%5.1f%%)' %
689 (100.0*thd, 100.0*max_thd)]
690 # smoothness of spectrum:
691 if 'dbdiff' in props:
692 db_diff = props['dbdiff']
693 msg += ['dBdiff=%5.1fdB' % db_diff]
694 if max_db_diff > 0.0 and db_diff > max_db_diff:
695 remove = True
696 skip_reason += ['not smooth s.d. diff=%5.1fdB (maxximumPowerDifference=%5.1fdB)' %
697 (db_diff, max_db_diff)]
698 # maximum power of higher harmonics:
699 if 'maxdb' in props:
700 max_harmonics = props['maxdb']
701 msg += ['max harmonics=%5.1fdB' % max_harmonics]
702 if max_harmonics > max_harmonics_db:
703 remove = True
704 skip_reason += ['maximum harmonics=%5.1fdB too strong (maximumHarmonicsPower=%5.1fdB)' %
705 (max_harmonics, max_harmonics_db)]
706 # relative amplitude of harmonics:
707 if harm_relampl is not None:
708 for k, max_relampl in enumerate([max_relampl_harm1, max_relampl_harm2, max_relampl_harm3]):
709 if k >= len(harm_relampl):
710 break
711 msg += ['ampl%d=%5.1f%%' % (k+1, 100.0*harm_relampl[k])]
712 if max_relampl > 0.0 and k < len(harm_relampl) and harm_relampl[k] >= max_relampl:
713 num_str = ['First', 'Second', 'Third']
714 skip_reason += ['distorted ampl%d=%5.1f%% (maximum%sHarmonicAmplitude=%5.1f%%)' %
715 (k+1, 100.0*harm_relampl[k], num_str[k], 100.0*max_relampl)]
716 return remove, ', '.join(skip_reason), ', '.join(msg)
719def pulse_quality(props, max_clipped_frac=0.1, max_rms_sem=0.0):
720 """Assess the quality of an EOD waveform of a pulse fish.
722 Parameters
723 ----------
724 props: dict
725 A dictionary with properties of the analyzed EOD waveform
726 as returned by `analyze_pulse()`.
727 max_clipped_frac: float
728 Maximum allowed fraction of clipped data.
729 max_rms_sem: float
730 If not zero, maximum allowed standard error of the data
731 relative to p-p amplitude.
733 Returns
734 -------
735 skip_reason: str
736 An empty string if the waveform is good, otherwise a string
737 indicating the failure.
738 msg: str
739 A textual representation of the values tested.
740 skipped_clipped: bool
741 True if waveform was skipped because of clipping.
742 """
743 msg = []
744 skip_reason = []
745 skipped_clipped = False
746 # clipped fraction:
747 if 'clipped' in props:
748 clipped_frac = props['clipped']
749 msg += ['clipped=%3.0f%%' % (100.0*clipped_frac)]
750 if clipped_frac >= max_clipped_frac:
751 skip_reason += ['clipped=%3.0f%% (maximumClippedFraction=%3.0f%%)' %
752 (100.0*clipped_frac, 100.0*max_clipped_frac)]
753 skipped_clipped = True
754 # noise:
755 rms_sem = None
756 if 'rmssem' in props:
757 rms_sem = props['rmssem']
758 if 'noise' in props:
759 rms_sem = props['noise']
760 if rms_sem is not None:
761 msg += ['rms sem waveform=%6.2f%%' % (100.0*rms_sem)]
762 if max_rms_sem > 0.0 and rms_sem >= max_rms_sem:
763 skip_reason += ['noisy waveform s.e.m.=%6.2f%% (maximumRMSNoise=%6.2f%%)' %
764 (100.0*rms_sem, 100.0*max_rms_sem)]
765 return ', '.join(skip_reason), ', '.join(msg), skipped_clipped
768def plot_eod_recording(ax, data, rate, unit=None, width=0.1,
769 toffs=0.0, rec_style=dict(lw=2, color='tab:red')):
770 """Plot a zoomed in range of the recorded trace.
772 Parameters
773 ----------
774 ax: matplotlib axes
775 Axes used for plotting.
776 data: 1D ndarray
777 Recorded data to be plotted.
778 rate: float
779 Sampling rate of the data in Hertz.
780 unit: str
781 Optional unit of the data used for y-label.
782 width: float
783 Width of data segment to be plotted in seconds.
784 toffs: float
785 Time of first data value in seconds.
786 rec_style: dict
787 Arguments passed on to the plot command for the recorded trace.
788 """
789 widx2 = int(width*rate)//2
790 i0 = len(data)//2 - widx2
791 i0 = (i0//widx2)*widx2
792 i1 = i0 + 2*widx2
793 if i0 < 0:
794 i0 = 0
795 if i1 >= len(data):
796 i1 = len(data)
797 time = np.arange(len(data))/rate + toffs
798 tunit = 'sec'
799 if np.abs(time[i0]) < 1.0 and np.abs(time[i1]) < 1.0:
800 time *= 1000.0
801 tunit = 'ms'
802 ax.plot(time, data, **rec_style)
803 ax.set_xlim(time[i0], time[i1])
805 ax.set_xlabel('Time [%s]' % tunit)
806 ymin = np.min(data[i0:i1])
807 ymax = np.max(data[i0:i1])
808 dy = ymax - ymin
809 ax.set_ylim(ymin-0.05*dy, ymax+0.05*dy)
810 if len(unit) == 0 or unit == 'a.u.':
811 ax.set_ylabel('Amplitude')
812 else:
813 ax.set_ylabel('Amplitude [%s]' % unit)
816def plot_eod_snippets(ax, data, rate, tmin, tmax, eod_times,
817 n_snippets=10, flip=False, aoffs=0,
818 snippet_style=dict(scaley=False,
819 lw=0.5, color='0.6')):
820 """Plot a few EOD waveform snippets.
822 Parameters
823 ----------
824 ax: matplotlib axes
825 Axes used for plotting.
826 data: 1D ndarray
827 Recorded data from which the snippets are taken.
828 rate: float
829 Sampling rate of the data in Hertz.
830 tmin: float
831 Start time of each snippet.
832 tmax: float
833 End time of each snippet.
834 eod_times: 1-D array
835 EOD peak times from which a few are selected to be plotted.
836 n_snippets: int
837 Number of snippets to be plotted. If zero do not plot anything.
838 flip: bool
839 If True flip the snippets upside down.
840 aoffs: float
841 Offset that was subtracted from the average EOD waveform.
842 snippet_style: dict
843 Arguments passed on to the plot command for plotting the snippets.
844 """
845 if data is None or n_snippets <= 0:
846 return
847 i0 = int(tmin*rate)
848 i1 = int(tmax*rate)
849 time = 1000.0*np.arange(i0, i1)/rate
850 step = len(eod_times)//n_snippets
851 if step < 1:
852 step = 1
853 for t in eod_times[n_snippets//2::step]:
854 idx = int(np.round(t*rate))
855 if idx + i0 < 0 or idx + i1 >= len(data):
856 continue
857 snippet = data[idx + i0:idx + i1] - aoffs
858 if flip:
859 snippet *= -1
860 ax.plot(time, snippet - np.mean(snippet[:len(snippet)//4]),
861 zorder=-5, **snippet_style)
864def plot_eod_waveform(ax, eod_waveform, props, phases=None,
865 unit=None, wave_periods=2,
866 magnification_factor=20,
867 wave_style=dict(lw=1.5, color='tab:red'),
868 magnified_style=dict(lw=0.8, color='tab:red'),
869 positive_style=dict(facecolor='tab:green', alpha=0.2,
870 edgecolor='none'),
871 negative_style=dict(facecolor='tab:blue', alpha=0.2,
872 edgecolor='none'),
873 sem_style=dict(color='0.8'),
874 fit_style=dict(lw=1.5, color='tab:blue'),
875 phase_style=dict(zorder=0, ls='', marker='o', color='tab:red',
876 markersize=5, mec='white', mew=1),
877 zerox_style=dict(zorder=50, ls='', marker='o', color='black',
878 markersize=5, mec='white', mew=1),
879 zero_style=dict(lw=0.5, color='0.7'),
880 fontsize='medium'):
881 """Plot mean EOD, its standard error, and an optional fit to the EOD.
883 Parameters
884 ----------
885 ax: matplotlib axes
886 Axes used for plotting.
887 eod_waveform: 2-D array
888 EOD waveform. First column is time in seconds, second column
889 the (mean) eod waveform. The optional third column is the
890 standard error, the optional fourth column is a fit of the
891 whole waveform, and the optional fourth column is a fit of
892 the tails of a pulse waveform.
893 props: dict
894 A dictionary with properties of the analyzed EOD waveform as
895 returned by `analyze_wave()` and `analyze_pulse()`.
896 phases: dict
897 Dictionary with phase properties as returned by
898 `analyze_pulse_phases()`, `analyze_pulse()`, and
899 `load_pulse_phases()`.
900 unit: str
901 Optional unit of the data used for y-label.
902 wave_periods: float
903 How many periods of a wave EOD are shown.
904 magnification_factor: float
905 If larger than one, plot a magnified version of the EOD
906 waveform magnified by this factor.
907 wave_style: dict
908 Arguments passed on to the plot command for the EOD waveform.
909 magnified_style: dict
910 Arguments passed on to the plot command for the magnified EOD waveform.
911 positive_style: dict
912 Arguments passed on to the fill_between command for coloring
913 positive phases.
914 negative_style: dict
915 Arguments passed on to the fill_between command for coloring
916 negative phases.
917 sem_style: dict
918 Arguments passed on to the fill_between command for the
919 standard error of the EOD.
920 fit_style: dict
921 Arguments passed on to the plot command for the fitted EOD.
922 phase_style: dict
923 Arguments passed on to the plot command for marking EOD phases.
924 zerox_style: dict
925 Arguments passed on to the plot command for marking zero crossings.
926 zero_style: dict
927 Arguments passed on to the plot command for the zero line.
928 fontsize: str or float or int
929 Fontsize for annotation text.
931 """
932 time = 1000*eod_waveform[:, 0]
933 eod = eod_waveform[:, 1]
934 # time axis:
935 if props is not None and props['type'] == 'wave':
936 period = 1000.0/props['EODf']
937 xlim = 0.5*wave_periods*period
938 xlim_l = -xlim
939 xlim_r = +xlim
940 elif props is not None and props['type'] == 'pulse':
941 # width of maximum peak:
942 meod = np.abs(eod)
943 ip = np.argmax(meod)
944 thresh = 0.5*meod[ip]
945 i0 = ip - np.argmax(meod[ip::-1] < thresh)
946 i1 = ip + np.argmax(meod[ip:] < thresh)
947 w = 4*(time[i1] - time[i0])
948 w = np.ceil(w/0.5)*0.5
949 # make sure tstart, tend, and time constant are included:
950 if props is not None:
951 if 'tstart' in props and 1000*props['tstart'] < -w:
952 w = np.ceil(abs(1000*props['tstart'])/0.5)*0.5
953 if 'tend' in props and 1000*props['tend'] > 2*w:
954 w = np.ceil(0.5*abs(1000*props['tend'])/0.5)*0.5
955 if 'taustart' in props and props['taustart'] is not None and \
956 1100*props['taustart'] > 2*w:
957 w = np.ceil(0.5*abs(1100*props['taustart'])/0.5)*0.5
958 # set xaxis limits:
959 xlim_l = -w
960 xlim_r = 2*w
961 xlim = (xlim_r - xlim_l)/2
962 else:
963 w = (time[-1] - time[0])/2
964 w = np.floor(w/0.5)*0.5
965 xlim_l = -w
966 xlim_r = +w
967 xlim = w
968 ax.set_xlim(xlim_l, xlim_r)
969 if xlim < 2:
970 ax.xaxis.set_major_locator(MultipleLocator(0.5))
971 elif xlim < 4:
972 ax.xaxis.set_major_locator(MultipleLocator(1))
973 elif xlim < 8:
974 ax.xaxis.set_major_locator(MultipleLocator(2))
975 ax.set_xlabel('Time [msec]')
976 # amplitude axis:
977 ylim = np.max(np.abs(eod[(time >= xlim_l) & (time <= xlim_r)]))
978 ax.set_ylim(-1.15*ylim, +1.15*ylim)
979 if unit:
980 ax.set_ylabel(f'Amplitude [{unit}]')
981 else:
982 ax.set_ylabel('Amplitude')
983 # ax height dimensions:
984 t = ax.text(0, 0, 'test', fontsize=fontsize)
985 fs = t.get_fontsize()
986 t.remove()
987 pixelx = np.abs(np.diff(ax.get_window_extent().get_points()[:, 0]))[0]
988 dxu = 2*xlim/pixelx
989 xfs = fs*dxu
990 pixely = np.abs(np.diff(ax.get_window_extent().get_points()[:, 1]))[0]
991 dyu = 2*ylim/pixely
992 yfs = fs*dyu
993 texts = []
994 quadrants = np.zeros((2, 2), dtype=int)
995 # magnification threshold:
996 if magnification_factor > 1:
997 mag_thresh = 0.95*np.max(np.abs(eod))/magnification_factor
998 else:
999 mag_thresh = 0
1000 # plot zero line:
1001 ax.axhline(0.0, zorder=10, **zero_style)
1002 # plot areas:
1003 if phases is not None and len(phases) > 0:
1004 if positive_style is not None and len(positive_style) > 0:
1005 ax.fill_between(time, eod, 0, eod >= 0, zorder=0,
1006 **positive_style)
1007 if negative_style is not None and len(negative_style) > 0:
1008 ax.fill_between(time, eod, 0, eod <= 0, zorder=0,
1009 **negative_style)
1010 # plot Fourier/Gaussian fit:
1011 if eod_waveform.shape[1] > 3 and np.all(np.isfinite(eod_waveform[:, 3])):
1012 ax.plot(time, eod_waveform[:, 3], zorder=30, **fit_style)
1013 # plot time constant fit:
1014 tau_magnified = False
1015 if eod_waveform.shape[1] > 4:
1016 if np.nanmax(np.abs(eod_waveform[:, 4])) < mag_thresh:
1017 tau_magnified = True
1018 ax.plot(time, magnification_factor*eod_waveform[:, 4],
1019 zorder=35, **fit_style)
1020 else:
1021 fs = dict(**fit_style)
1022 if 'lw' in fs:
1023 fs['lw'] *= 2
1024 ax.plot(time, eod_waveform[:, 4], zorder=35, **fs)
1025 # plot waveform:
1026 ax.plot(time, eod, zorder=45, **wave_style)
1027 # plot standard error:
1028 if eod_waveform.shape[1] > 2:
1029 std_eod = eod_waveform[:, 2]
1030 ax.fill_between(time, eod + std_eod, eod - std_eod,
1031 zorder=20, **sem_style)
1032 # plot magnified pulse waveform:
1033 magnification_mask = np.zeros(len(time), dtype=bool)
1034 if magnification_factor > 1 and phases is not None and len(phases) > 0:
1035 i0 = np.argmax(np.abs(eod) > mag_thresh)
1036 if i0 > 0:
1037 left_eod = magnification_factor*eod[:i0]
1038 magnification_mask[:i0] = True
1039 ax.plot(time[:i0], left_eod, zorder=40, **magnified_style)
1040 ta = None
1041 if left_eod[-1] > 0:
1042 it = np.argmax(left_eod > 0.95*np.max(eod))
1043 if it < len(left_eod)//2:
1044 it = len(left_eod) - 1
1045 if time[it] - 3*xfs > xlim_l:
1046 ty = left_eod[it] if left_eod[it] < np.max(eod) else np.max(eod)
1047 ta = ax.text(time[it], ty, f'x{magnification_factor:.0f} ',
1048 ha='right', va='top', zorder=100,
1049 fontsize=fontsize)
1050 if ty > 0.5*ylim:
1051 quadrants[0, 0] += 1
1052 else:
1053 it = np.argmax(left_eod < 0.95*np.min(eod))
1054 if it < len(left_eod)//2:
1055 it = len(left_eod) - 1
1056 if time[it] - 3*xfs > xlim_l:
1057 ty = left_eod[it] if left_eod[it] > np.min(eod) else np.min(eod)
1058 ta = ax.text(time[it], ty, f'x{magnification_factor:.0f} ',
1059 ha='right', va='bottom', zorder=100,
1060 fontsize=fontsize)
1061 if ty < -0.5*ylim:
1062 quadrants[1, 0] += 1
1063 if ta is not None:
1064 texts.append(ta)
1065 left_snip = left_eod[time[:i0] < 0.1*xlim_l]
1066 if len(left_snip) > 0:
1067 if quadrants[0, 0] == 0:
1068 quadrants[0, 0] += np.max(left_snip) > 0.5*ylim
1069 if quadrants[1, 0] == 0:
1070 quadrants[1, 0] += np.min(left_snip) < -0.5*ylim
1071 i1 = len(eod) - np.argmax(np.abs(eod[::-1]) > mag_thresh)
1072 right_eod = magnification_factor*eod[i1:]
1073 magnification_mask[i1:] = True
1074 ax.plot(time[i1:], right_eod, zorder=40, **magnified_style)
1075 right_snip = right_eod[time[i1:] > 0.4*xlim_r]
1076 if len(right_snip) > 0:
1077 quadrants[0, 1] += np.max(right_snip) > 0.5*ylim
1078 quadrants[1, 1] += np.min(right_snip) < -0.5*ylim
1079 # annotate time constant fit:
1080 tau = None if props is None else props.get('tau', None)
1081 if tau is not None and eod_waveform.shape[1] > 4:
1082 if tau < 0.001:
1083 label = f'\u03c4={1.e6*tau:.0f}\u00b5s'
1084 else:
1085 label = f'\u03c4={1.e3*tau:.2f}ms'
1086 inx = np.argmin(np.isnan(eod_waveform[:, 4]))
1087 x0 = time[inx]
1088 x = x0 + 1000*np.log(2)*tau
1089 if x + 4*xfs >= xlim_r:
1090 if xlim_r - x0 >= 4*xfs:
1091 x = xlim_r - 8*xfs
1092 else:
1093 x = x0
1094 elif x + 8*xfs > xlim_r:
1095 x = xlim_r - 8*xfs
1096 if x < x0:
1097 x = x0
1098 y = eod_waveform[np.argmin(np.abs(time - x)), 4]
1099 if tau_magnified:
1100 y *= magnification_factor
1101 va = 'bottom' if eod[inx] > 0 else 'top'
1102 if eod[inx] < 0:
1103 y -= 0.5*yfs
1104 ta = ax.text(x + xfs, y, label, ha='left', va=va,
1105 zorder=100, fontsize=fontsize)
1106 texts.append(ta)
1107 if x + xfs > 0.4*xlim_r:
1108 if y > 0.5*ylim:
1109 quadrants[0, 1] += 1
1110 elif y < -0.5*ylim:
1111 quadrants[1, 1] += 1
1112 if props is not None:
1113 # mark start and end:
1114 if 'tstart' in props:
1115 ax.axvline(1000*props['tstart'], 0.45, 0.55,
1116 color='k', lw=0.5, zorder=25)
1117 if 'tend' in props:
1118 ax.axvline(1000*props['tend'], 0.45, 0.55,
1119 color='k', lw=0.5, zorder=25)
1120 # mark cumulative:
1121 if 'median' in props:
1122 y = -1.07*ylim
1123 m = 1000*props['median']
1124 q1 = 1000*props['quartile1']
1125 q3 = 1000*props['quartile3']
1126 w = q3 - q1
1127 ax.plot([q1, q3], [y, y], 'gray', lw=4, zorder=25)
1128 ax.plot(m, y, 'o', color='white', ms=3, zorder=26)
1129 label = f'{w:.2f}ms' if w >= 1 else f'{1000*w:.0f}\u00b5s'
1130 ax.text(q3 + xfs, y, label,
1131 va='center', zorder=100, fontsize=fontsize)
1132 # plot and annotate phases:
1133 if phases is not None and len(phases) > 0:
1134 upper_area_text = False
1135 lower_area_text = False
1136 # mark zero crossings:
1137 zeros = 1000*phases['zeros']
1138 ax.plot(zeros, np.zeros(len(zeros)), **zerox_style)
1139 # phase peaks and troughs:
1140 max_peak_idx = np.argmax(phases['amplitudes'])
1141 min_trough_idx = np.argmin(phases['amplitudes'])
1142 for i in range(len(phases['times'])):
1143 index = phases['indices'][i]
1144 ptime = 1000*phases['times'][i]
1145 if ptime < xlim_l or ptime > xlim_r:
1146 continue
1147 pi = np.argmin(np.abs(time - ptime))
1148 mfac = magnification_factor if magnification_mask[pi] else 1
1149 pampl = mfac*phases['amplitudes'][i]
1150 relampl = phases['relamplitudes'][i]
1151 relarea = phases['relareas'][i]
1152 # classify phase:
1153 ampl_phase = phases['amplitudes'][i]
1154 ampl_left = phases['amplitudes'][i - 1] if i > 0 else 0
1155 ampl_right = phases['amplitudes'][i + 1] if i + 1 < len(phases['amplitudes']) else 0
1156 local_maximum = ampl_phase > ampl_left and ampl_phase > ampl_right
1157 if local_maximum:
1158 right_phase = (i >= max_peak_idx)
1159 min_max_phase = (i == max_peak_idx)
1160 local_phase = (ampl_phase < 0)
1161 else:
1162 right_phase = i >= min_trough_idx
1163 min_max_phase = (i == min_trough_idx)
1164 local_phase = (ampl_phase > 0)
1165 sign = np.sign(pampl)
1166 # mark phase peak/trough:
1167 ax.plot(ptime, pampl, **phase_style)
1168 # text for phase label:
1169 label = f'P{index:.0f}'
1170 if index != 1 and not local_phase:
1171 if np.abs(ptime) < 1:
1172 ts = f'{1000*ptime:.0f}\u00b5s'
1173 elif np.abs(ptime) < 10:
1174 ts = f'{ptime:.2f}ms'
1175 else:
1176 ts = f'{ptime:.3g}ms'
1177 if np.abs(relampl) < 0.05:
1178 ps = f'{100*relampl:.1f}%'
1179 else:
1180 ps = f'{100*relampl:.0f}%'
1181 label += f'({ps} @ {ts})'
1182 # position of phase label:
1183 ltime = ptime
1184 lampl = pampl
1185 valign = 'top' if sign < 0 else 'baseline'
1186 add = True
1187 if local_phase or (min_max_phase and abs(pampl)/ylim < 0.8):
1188 halign = 'center'
1189 dx = 0
1190 dy = 0.6*yfs
1191 if local_phase:
1192 add = False
1193 elif min_max_phase:
1194 halign = 'left' if right_phase else 'right'
1195 dx = xfs if right_phase else -xfs
1196 dy = 0
1197 if abs(relampl) > 0.85:
1198 dx *= 2
1199 dy = -1.5*yfs
1200 else:
1201 dx = 0
1202 dy = 0.8*yfs
1203 if right_phase:
1204 halign = 'left'
1205 if i > 0 and np.isfinite(phases['zeros'][i - 1]):
1206 ltime = 1000*phases['zeros'][i - 1]
1207 else:
1208 dx = -2*xfs
1209 #np.sum(phases['amplitudes'][i + 1:]*pampl > 0)
1210 else:
1211 halign = 'right'
1212 if np.isfinite(phases['zeros'][i]):
1213 ltime = 1000*phases['zeros'][i]
1214 else:
1215 dx = 2*xfs
1216 if sign < 0:
1217 dy = -dy
1218 ta = ax.text(ltime + dx, lampl + dy, label,
1219 ha=halign, va=valign, zorder=100, fontsize=fontsize)
1220 if add:
1221 texts.append(ta)
1222 # area:
1223 if np.abs(relarea) < 0.01:
1224 continue
1225 elif np.abs(relarea) < 0.05:
1226 label = f'{100*relarea:.1f}%'
1227 else:
1228 label = f'{100*relarea:.0f}%'
1229 x = ptime
1230 if i > 0 and i < len(phases['times']) - 1:
1231 xl = 1000*phases['times'][i - 1]
1232 xr = 1000*phases['times'][i + 1]
1233 tsnippet = time[(time > xl) & (time < xr)]
1234 snippet = eod[(time > xl) & (time < xr)]
1235 tsnippet = tsnippet[np.sign(pampl)*snippet > 0]
1236 snippet = snippet[np.sign(pampl)*snippet > 0]
1237 x = np.sum(tsnippet*snippet)/np.sum(snippet)
1238 if abs(relampl) > 0.5:
1239 ax.text(x, sign*0.6*yfs, label,
1240 rotation='vertical',
1241 va='top' if sign < 0 else 'bottom',
1242 ha='center', zorder=35, fontsize=fontsize)
1243 elif abs(relampl) > 0.25 and abs(relarea) > 0.19:
1244 ax.text(x, sign*0.6*yfs, label,
1245 va='top' if sign < 0 else 'baseline',
1246 ha='center', zorder=35, fontsize=fontsize)
1247 else:
1248 dx = 0.5*xfs if right_phase else -0.5*xfs
1249 ta = ax.text(ltime + dx, -sign*0.4*yfs, label,
1250 va='baseline' if sign < 0 else 'top',
1251 ha=halign, zorder=100, fontsize=fontsize)
1252 if -sign > 0 and not upper_area_text:
1253 texts.append(ta)
1254 upper_area_text = True
1255 if -sign < 0 and not lower_area_text:
1256 texts.append(ta)
1257 lower_area_text = True
1258 # arrange text vertically to avoid overlaps:
1259 ul_texts = []
1260 ur_texts = []
1261 ll_texts = []
1262 lr_texts = []
1263 for t in texts:
1264 x, y = t.get_position()
1265 if y > 0:
1266 if x >= phases['times'][max_peak_idx]:
1267 ur_texts.append(t)
1268 else:
1269 ul_texts.append(t)
1270 else:
1271 if x >= phases['times'][min_trough_idx]:
1272 lr_texts.append(t)
1273 else:
1274 ll_texts.append(t)
1275 for ts, (j, k) in zip([ul_texts, ur_texts, ll_texts, lr_texts],
1276 [(0, 0), (0, 1), (1, 0), (1, 1)]):
1277 if len(ts) > 1:
1278 ys = []
1279 for t in ts:
1280 # alternative:
1281 #renderer = ax.get_fig().canvas.renderer
1282 #bbox = t.get_window_extent(renderer).transformed(ax.transData.inverted())
1283 x, y = t.get_position()
1284 ys.append(abs(y))
1285 idx = np.argsort(ys)
1286 x, y = ts[idx[0]].get_position()
1287 yp = abs(y)
1288 for i in idx[1:]:
1289 t = ts[i]
1290 x, y = t.get_position()
1291 s = t.get_text()
1292 if abs(y) < abs(yp) + 2*yfs and \
1293 len(s) > 4 and s[:2] != '\u03c4=':
1294 y = np.sign(y)*(abs(yp) + 2*yfs)
1295 t.set_y(y)
1296 if len(s) >= 4 and abs(y) > 0.5*ylim:
1297 quadrants[j, k] += 1
1298 yp = y
1299 # annotate plot:
1300 if unit is None or len(unit) == 0 or unit == 'a.u.':
1301 unit = ''
1302 if props is not None:
1303 label = '' # f'p-p amplitude = {props["p-p-amplitude"]:.3g} {unit}\n'
1304 if 'n' in props:
1305 eods = 'EODs' if props['n'] > 1 else 'EOD'
1306 label += f'n = {props["n"]} {eods}\n'
1307 if 'flipped' in props and props['flipped']:
1308 label += 'flipped\n'
1309 if 'polaritybalance' in props:
1310 label += f'PB={100*props["polaritybalance"]:.0f} %\n'
1311 # weigh left quadrants less:
1312 quadrants *= 2
1313 quadrants[quadrants[:, 1] > 0, 1] -= 1
1314 # find free quadrant:
1315 q_row, q_col = np.unravel_index(np.argmin(quadrants), quadrants.shape)
1316 # place text in quadrant:
1317 y = 1 if q_row == 0 else 0
1318 va = 'top' if q_row == 0 else 'bottom'
1319 x = 0.03 if q_col == 0 else 0.97
1320 ha = 'left' if q_col == 0 else 'right'
1321 ax.text(x, y, label, transform=ax.transAxes,
1322 va=va, ha=ha, zorder=100)
1325def save_eod_waveform(mean_eod, unit, idx, basename, **kwargs):
1326 """Save mean EOD waveform to file.
1328 Parameters
1329 ----------
1330 mean_eod: 2D array of floats
1331 Averaged EOD waveform as returned by `eod_waveform()`,
1332 `analyze_wave()`, and `analyze_pulse()`.
1333 unit: str
1334 Unit of the waveform data.
1335 idx: int or None
1336 Index of fish.
1337 basename: str or Path or stream
1338 If string, path and basename of file.
1339 If `basename` does not have an extension,
1340 '-eodwaveform', the fish index, and a file extension are appended.
1341 If stream, write EOD waveform data into this stream.
1342 kwargs:
1343 Arguments passed on to `TableData.write()`.
1345 Returns
1346 -------
1347 filename: Path
1348 Path and full name of the written file in case of `basename`
1349 being a string. Otherwise, the file name and extension that
1350 would have been appended to a basename.
1352 See Also
1353 --------
1354 load_eod_waveform()
1355 """
1356 td = TableData(mean_eod[:, :3]*[1000.0, 1.0, 1.0],
1357 ['time', 'mean', 'sem'],
1358 ['ms', unit, unit],
1359 ['%.3f', '%.6g', '%.6g'])
1360 if mean_eod.shape[1] > 3:
1361 td.append('fit', unit, '%.5f', value=mean_eod[:, 3])
1362 if mean_eod.shape[1] > 4:
1363 td.append('tailfit', unit, '%.5f', value=mean_eod[:, 4])
1364 fp = ''
1365 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
1366 if not ext:
1367 fp = '-eodwaveform'
1368 if idx is not None:
1369 fp += f'-{idx}'
1370 return td.write_file_stream(basename, fp, **kwargs)
1373def load_eod_waveform(file_path):
1374 """Load EOD waveform from file.
1376 Parameters
1377 ----------
1378 file_path: str or Path
1379 Path of the file to be loaded.
1381 Returns
1382 -------
1383 mean_eod: 2D array of floats
1384 Averaged EOD waveform: time in seconds, mean, standard deviation, fit.
1385 unit: str
1386 Unit of EOD waveform.
1388 Raises
1389 ------
1390 FileNotFoundError:
1391 If `file_path` does not exist.
1393 See Also
1394 --------
1395 save_eod_waveform()
1396 """
1397 data = TableData(file_path)
1398 mean_eod = data.array()
1399 mean_eod[:, 0] *= 0.001
1400 return mean_eod, data.unit('mean')
1403file_types = ['waveeodfs', 'wavefish', 'pulsefish', 'eodwaveform',
1404 'wavespectrum', 'pulsephases', 'pulsegaussians', 'pulsespectrum', 'pulsetimes']
1405"""List of all file types generated and supported by the `save_*` and `load_*` functions."""
1408def parse_filename(file_path):
1409 """Parse components of an EOD analysis file name.
1411 Analysis files generated by the `eodanalysis` module are named
1412 according to
1413 ```plain
1414 PATH/RECORDING-CHANNEL-TIME-FTYPE-N.EXT
1415 ```
1417 Parameters
1418 ----------
1419 file_path: str or Path
1420 Path of the file to be parsed.
1422 Returns
1423 -------
1424 recording: str
1425 Path and basename of the recording, i.e. 'PATH/RECORDING'.
1426 A leading './' is removed.
1427 base_path: Path
1428 Path and basename of the analysis results,
1429 i.e. 'PATH/RECORDING-CHANNEL-TIME'.
1430 channel: int
1431 Channel of the recording
1432 ('CHANNEL' component of the file name if present).
1433 -1 if not present in `file_path`.
1434 time: float
1435 Start time of analysis window in seconds
1436 ('TIME' component of the file name if present).
1437 `None` if not present in `file_path`.
1438 ftype: str
1439 Type of analysis file (e.g. 'wavespectrum', 'pulsephases', etc.),
1440 ('FTYPE' component of the file name if present).
1441 See `file_types` for a list of all supported file types.
1442 Empty string if not present in `file_path`.
1443 index: int
1444 Index of the EOD.
1445 ('N' component of the file name if present).
1446 -1 if not present in `file_path`.
1447 ext: str
1448 File extension *without* leading period
1449 ('EXT' component of the file name).
1451 """
1452 file_path = Path(file_path)
1453 ext = file_path.suffix
1454 ext = ext[1:]
1455 parts = file_path.stem.split('-')
1456 index = -1
1457 if len(parts) > 0 and parts[-1].isdigit():
1458 index = int(parts[-1])
1459 parts = parts[:-1]
1460 ftype = ''
1461 if len(parts) > 0:
1462 ftype = parts[-1]
1463 parts = parts[:-1]
1464 base_path = file_path.parent / '-'.join(parts)
1465 time = None
1466 if len(parts) > 0 and len(parts[-1]) > 0 and \
1467 parts[-1][0] == 't' and parts[-1][-1] == 's' and \
1468 parts[-1][1:-1].isdigit():
1469 time = float(parts[-1][1:-1])
1470 parts = parts[:-1]
1471 channel = -1
1472 if len(parts) > 0 and len(parts[-1]) > 0 and \
1473 parts[-1][0] == 'c' and parts[-1][1:].isdigit():
1474 channel = int(parts[-1][1:])
1475 parts = parts[:-1]
1476 recording = '-'.join(parts)
1477 return recording, base_path, channel, time, ftype, index, ext
1480def save_analysis(output_basename, zip_file, eod_props, mean_eods, spec_data,
1481 phase_data, pulse_data, wave_eodfs, wave_indices, unit,
1482 verbose, **kwargs):
1483 """Save EOD analysis results to files.
1485 Parameters
1486 ----------
1487 output_basename: str or Path
1488 Path and basename of files to be saved.
1489 zip_file: bool
1490 If `True`, write all analysis results into a zip archive.
1491 eod_props: list of dict
1492 Properties of EODs as returned by `analyze_wave()` and
1493 `analyze_pulse()`.
1494 mean_eods: list of 2D array of floats
1495 Averaged EOD waveforms as returned by `eod_waveform()`,
1496 `analyze_wave()`, and `analyze_pulse()`.
1497 spec_data: list of 2D array of floats
1498 Energy spectra of single pulses as returned by
1499 `analyze_pulse()`.
1500 phase_data: list of dict
1501 Properties of phases of pulse EODs as returned by
1502 `analyze_pulse()` and `analyze_pulse_phases()`.
1503 pulse_data: list of dict
1504 For each pulse fish a dictionary with phase times, amplitudes and standard
1505 deviations of Gaussians fitted to the pulse waveform.
1506 wave_eodfs: list of 2D array of float
1507 Each item is a matrix with the frequencies and powers
1508 (columns) of the fundamental and harmonics (rows) as returned
1509 by `harmonics.harmonic_groups()`.
1510 wave_indices: array of int
1511 Indices identifying each fish in `wave_eodfs` or NaN.
1512 unit: str
1513 Unit of the waveform data.
1514 verbose: int
1515 Verbosity level.
1516 kwargs:
1517 Arguments passed on to `TableData.write()`.
1518 """
1519 def write_file_zip(zf, save_func, output, *args, **kwargs):
1520 if zf is None:
1521 fp = save_func(*args, basename=output, **kwargs)
1522 if verbose > 0 and fp is not None:
1523 print('wrote file', fp)
1524 else:
1525 with io.StringIO() as df:
1526 fp = save_func(*args, basename=df, **kwargs)
1527 if fp is not None:
1528 fp = Path(output.stem + str(fp))
1529 zf.writestr(fp.name, df.getvalue())
1530 if verbose > 0:
1531 print('zipped file', fp.name)
1534 output_basename = Path(output_basename)
1535 if 'table_format' in kwargs and kwargs['table_format'] == 'py':
1536 with open(output_basename.with_suffix('.py'), 'w') as f:
1537 name = output_basename.stem
1538 for k in range(len((spec_data))):
1539 species = eod_props[k].get('species', '')
1540 if len(pulse_data[k]) > 0:
1541 fish = normalize_pulsefish(pulse_data[k])
1542 export_pulsefish(fish, f'{name}-{k}_phases',
1543 species, f)
1544 f.write('\n')
1545 else:
1546 sdata = spec_data[k]
1547 if len(sdata) > 0 and sdata.shape[1] > 2:
1548 fish = dict(amplitudes=sdata[:, 3], phases=sdata[:, 5])
1549 fish = normalize_wavefish(fish)
1550 export_wavefish(fish, f'{name}-{k}_harmonics',
1551 species, f)
1552 f.write('\n')
1553 else:
1554 zf = None
1555 if zip_file:
1556 zf = zipfile.ZipFile(output_basename.with_suffix('.zip'), 'w')
1557 # all wave fish in wave_eodfs:
1558 if len(wave_eodfs) > 0:
1559 write_file_zip(zf, save_wave_eodfs, output_basename,
1560 wave_eodfs, wave_indices, **kwargs)
1561 # all wave and pulse fish:
1562 for i, (mean_eod, sdata, pdata, pulse, props) in enumerate(zip(mean_eods, spec_data, phase_data,
1563 pulse_data, eod_props)):
1564 write_file_zip(zf, save_eod_waveform, output_basename,
1565 mean_eod, unit, i, **kwargs)
1566 # spectrum:
1567 if len(sdata)>0:
1568 if sdata.shape[1] <= 3:
1569 write_file_zip(zf, save_pulse_spectrum, output_basename,
1570 sdata, unit, i, **kwargs)
1571 else:
1572 write_file_zip(zf, save_wave_spectrum, output_basename,
1573 sdata, unit, i, **kwargs)
1574 # phases:
1575 write_file_zip(zf, save_pulse_phases, output_basename,
1576 pdata, unit, i, **kwargs)
1577 # pulses:
1578 write_file_zip(zf, save_pulse_gaussians, output_basename,
1579 pulse, unit, i, **kwargs)
1580 # times:
1581 write_file_zip(zf, save_pulse_times, output_basename,
1582 props, i, **kwargs)
1583 # wave fish properties:
1584 write_file_zip(zf, save_wave_fish, output_basename,
1585 eod_props, unit, **kwargs)
1586 # pulse fish properties:
1587 write_file_zip(zf, save_pulse_fish, output_basename,
1588 eod_props, unit, **kwargs)
1591def load_analysis(file_pathes):
1592 """Load all EOD analysis files.
1594 Parameters
1595 ----------
1596 file_pathes: list of str or Path
1597 Pathes of the analysis files of a single recording to be loaded.
1599 Returns
1600 -------
1601 mean_eods: list of 2D array of floats
1602 Averaged EOD waveforms: time in seconds, mean, standard deviation, fit.
1603 wave_eodfs: 2D array of floats
1604 EODfs and power of wave type fish.
1605 wave_indices: array of ints
1606 Corresponding indices of fish, can contain negative numbers to
1607 indicate frequencies without fish.
1608 eod_props: list of dict
1609 Properties of EODs. The 'index' property is an index into the
1610 reurned lists.
1611 spec_data: list of 2D array of floats
1612 Amplitude and phase spectrum of wave-type EODs with columns
1613 harmonics, frequency, amplitude, relative amplitude in dB,
1614 relative power in dB, phase, data power in unit squared.
1615 Energy spectrum of single pulse-type EODs with columns
1616 frequency and energy.
1617 phase_data: list of dict
1618 Properties of phases of pulse-type EODs with keys
1619 indices, times, amplitudes, relamplitudes, widths, areas, relareas, zeros
1620 pulse_data: list of dict
1621 For each pulse fish a dictionary with phase times, amplitudes and standard
1622 deviations of Gaussians fitted to the pulse waveform. Use the
1623 functions provided in thunderfish.fakefish to simulate pulse
1624 fish EODs from this data.
1625 recording: str
1626 Path and base name of the recording file.
1627 channel: int
1628 Analysed channel of the recording.
1629 unit: str
1630 Unit of EOD waveform.
1631 """
1632 recording = None
1633 channel = -1
1634 eod_props = []
1635 zf = None
1636 if len(file_pathes) == 1 and Path(file_pathes[0]).suffix[1:] == 'zip':
1637 zf = zipfile.ZipFile(file_pathes[0])
1638 file_pathes = sorted(zf.namelist())
1639 # read wave- and pulse-fish summaries:
1640 pulse_fish = False
1641 wave_fish = False
1642 for f in file_pathes:
1643 recording, _, channel, _, ftype, _, _ = parse_filename(f)
1644 if zf is not None:
1645 f = io.TextIOWrapper(zf.open(f, 'r'))
1646 if ftype == 'wavefish':
1647 eod_props.extend(load_wave_fish(f))
1648 wave_fish = True
1649 elif ftype == 'pulsefish':
1650 eod_props.extend(load_pulse_fish(f))
1651 pulse_fish = True
1652 idx_offs = 0
1653 if wave_fish and not pulse_fish:
1654 idx_offs = sorted([ep['index'] for ep in eod_props])[0]
1655 # load all other files:
1656 neods = len(eod_props)
1657 if neods < 1:
1658 neods = 1
1659 eod_props = [None]
1660 wave_eodfs = np.array([])
1661 wave_indices = np.array([])
1662 mean_eods = [None]*neods
1663 spec_data = [None]*neods
1664 phase_data = [None]*neods
1665 pulse_data = [None]*neods
1666 unit = None
1667 for f in file_pathes:
1668 recording, _, channel, _, ftype, idx, _ = parse_filename(f)
1669 if neods == 1 and idx > 0:
1670 idx = 0
1671 idx -= idx_offs
1672 if zf is not None:
1673 f = io.TextIOWrapper(zf.open(f, 'r'))
1674 if ftype == 'waveeodfs':
1675 wave_eodfs, wave_indices = load_wave_eodfs(f)
1676 elif ftype == 'eodwaveform':
1677 mean_eods[idx], unit = load_eod_waveform(f)
1678 elif ftype == 'wavespectrum':
1679 spec_data[idx], unit = load_wave_spectrum(f)
1680 elif ftype == 'pulsephases':
1681 phase_data[idx], unit = load_pulse_phases(f)
1682 elif ftype == 'pulsegaussians':
1683 pulse_data[idx], unit = load_pulse_gaussians(f)
1684 elif ftype == 'pulsetimes':
1685 pulse_times = load_pulse_times(f)
1686 eod_props[idx]['times'] = pulse_times
1687 eod_props[idx]['peaktimes'] = pulse_times
1688 elif ftype == 'pulsespectrum':
1689 spec_data[idx] = load_pulse_spectrum(f)
1690 # fix wave spectra:
1691 wave_eodfs = [fish.reshape(1, 2) if len(fish)>0 else fish
1692 for fish in wave_eodfs]
1693 if len(wave_eodfs) > 0 and len(spec_data) > 0:
1694 eodfs = []
1695 for idx, fish in zip(wave_indices, wave_eodfs):
1696 if idx >= 0:
1697 spec = spec_data[idx]
1698 specd = np.zeros((np.sum(np.isfinite(spec[:, -1])),
1699 2))
1700 specd[:, 0] = spec[np.isfinite(spec[:, -1]),1]
1701 specd[:, 1] = spec[np.isfinite(spec[:, -1]),-1]
1702 eodfs.append(specd)
1703 else:
1704 specd = np.zeros((10, 2))
1705 specd[:, 0] = np.arange(len(specd))*fish[0,0]
1706 specd[:, 1] = np.nan
1707 eodfs.append(specd)
1708 wave_eodfs = eodfs
1709 return mean_eods, wave_eodfs, wave_indices, eod_props, spec_data, \
1710 phase_data, pulse_data, recording, channel, unit
1713def load_recording(file_path, channel=0, load_kwargs={},
1714 eod_props=None, verbose=0):
1715 """Load recording.
1717 Parameters
1718 ----------
1719 file_path: str or Path
1720 Full path of the file with the recorded data.
1721 Extension is optional. If absent, look for the first file
1722 with a reasonable extension.
1723 channel: int
1724 Channel of the recording to be returned.
1725 load_kwargs: dict
1726 Keyword arguments that are passed on to the
1727 format specific loading functions.
1728 eod_props: list of dict or None
1729 List of EOD properties from which start and end times of
1730 analysis window are extracted.
1731 verbose: int
1732 Verbosity level passed on to load function.
1734 Returns
1735 -------
1736 data: array of float
1737 Data of the requested `channel`.
1738 rate: float
1739 Sampling rate in Hertz.
1740 idx0: int
1741 Start index of the analysis window.
1742 idx1: int
1743 End index of the analysis window.
1744 info_dict: dict
1745 Dictionary with path, name, species, channel, chanstr, time.
1746 """
1747 data = None
1748 rate = 0.0
1749 idx0 = 0
1750 idx1 = 0
1751 info_dict = dict(path='',
1752 name='',
1753 species='',
1754 channel=0,
1755 chanstr='',
1756 time='')
1757 for k in range(1, 10):
1758 info_dict[f'part{k}'] = ''
1759 data_file = Path()
1760 file_path = Path(file_path)
1761 if len(file_path.suffix) > 1:
1762 data_file = file_path
1763 else:
1764 data_files = file_path.parent.glob(file_path.stem + '*')
1765 for dfile in data_files:
1766 if not dfile.suffix[1:] in ['zip'] + list(TableData.ext_formats.values()):
1767 data_file = dfile
1768 break
1769 if data_file.is_file():
1770 all_data = DataLoader(data_file, verbose=verbose, **load_kwargs)
1771 rate = all_data.rate
1772 unit = all_data.unit
1773 ampl_max = all_data.ampl_max
1774 data = all_data[:, channel]
1775 species = get_str(all_data.metadata(), ['species'], default='')
1776 if len(species) > 0:
1777 species += ' '
1778 info_dict.update(path=os.fsdecode(all_data.filepath),
1779 name=all_data.basename(),
1780 species=species,
1781 channel=channel)
1782 offs = 1
1783 for k, part in enumerate(all_data.filepath.parts):
1784 if k == 0 and part == all_data.filepath.anchor:
1785 offs = 0
1786 continue
1787 if part == all_data.filepath.name:
1788 break
1789 info_dict[f'part{k + offs}'] = part
1790 if all_data.channels > 1:
1791 if all_data.channels > 100:
1792 info_dict['chanstr'] = f'-c{channel:03d}'
1793 elif all_data.channels > 10:
1794 info_dict['chanstr'] = f'-c{channel:02d}'
1795 else:
1796 info_dict['chanstr'] = f'-c{channel:d}'
1797 else:
1798 info_dict['chanstr'] = ''
1799 idx0 = 0
1800 idx1 = len(data)
1801 if eod_props is not None and len(eod_props) > 0 and 'twin' in eod_props[0]:
1802 idx0 = int(eod_props[0]['twin']*rate)
1803 if len(eod_props) > 0 and 'window' in eod_props[0]:
1804 idx1 = idx0 + int(eod_props[0]['window']*rate)
1805 info_dict['time'] = f'-t{idx0/rate:.0f}s'
1806 all_data.close()
1808 return data, rate, idx0, idx1, info_dict
1811def add_eod_analysis_config(cfg, win_fac=2.0, min_win=0.01, max_eods=None,
1812 min_sem=False, unfilter_cutoff=0.0,
1813 flip_wave='none', flip_pulse='none',
1814 n_harm=10, min_pulse_win=0.001,
1815 start_end_thresh_fac=0.01, peak_thresh_fac=0.002,
1816 min_dist=50.0e-6, width_frac = 0.5,
1817 fit_frac = 0.5, fit_gaussians=True,
1818 freq_resolution=1.0, fade_frac=0.0,
1819 ipi_cv_thresh=0.5, ipi_percentile=30.0):
1820 """Add all parameters needed for the eod analysis functions as a new
1821 section to a configuration.
1823 Parameters
1824 ----------
1825 cfg: ConfigFile
1826 The configuration.
1828 See `eod_waveform()`, `analyze_wave()`, and `analyze_pulse()` for
1829 details on the remaining arguments.
1830 """
1831 cfg.add_section('EOD analysis:')
1832 cfg.add('eodSnippetFac', win_fac, '', 'The duration of EOD snippets is the EOD period times this factor.')
1833 cfg.add('eodMinSnippet', min_win, 's', 'Minimum duration of cut out EOD snippets.')
1834 cfg.add('eodMaxEODs', max_eods or 0, '', 'The maximum number of EODs used to compute the average EOD. If 0 use all EODs.')
1835 cfg.add('eodMinSem', min_sem, '', 'Use minimum of s.e.m. to set maximum number of EODs used to compute the average EOD.')
1836 cfg.add('unfilterCutoff', unfilter_cutoff, 'Hz', 'If non-zero remove effect of high-pass filter with this cut-off frequency.')
1837 cfg.add('flipWaveEOD', flip_wave, '', 'Flip EOD of wave fish to make largest extremum positive (flip, none, or auto).')
1838 cfg.add('flipPulseEOD', flip_pulse, '', 'Flip EOD of pulse fish to make the first large peak positive (flip, none, or auto).')
1839 cfg.add('eodHarmonics', n_harm, '', 'Number of harmonics fitted to the EOD waveform.')
1840 cfg.add('eodMinPulseSnippet', min_pulse_win, 's', 'Minimum duration of cut out EOD snippets for a pulse fish.')
1841 cfg.add('eodPeakThresholdFactor', peak_thresh_fac, '', 'Threshold for detection of peaks and troughs in pulse EODs as a fraction of the p-p amplitude.')
1842 cfg.add('eodStartEndThresholdFactor', start_end_thresh_fac, '', 'Threshold for for start and end time of pulse EODs as a fraction of the p-p amplitude.')
1843 cfg.add('eodMinimumDistance', min_dist, 's', 'Minimum distance between peaks and troughs in a EOD pulse.')
1844 cfg.add('eodPulseWidthFraction', 100*width_frac, '%', 'The width of a pulse is measured at this fraction of the pulse height.')
1845 cfg.add('eodExponentialFitFraction', 100*fit_frac, '%', 'An exponential function is fitted on the tail of a pulse starting at this fraction of the height of the last peak.')
1846 cfg.add('eodFitGaussians', fit_gaussians, '', 'Fit sum of Gaussians to pulse-type EOD waveform.')
1847 cfg.add('eodPulseFrequencyResolution', freq_resolution, 'Hz', 'Frequency resolution of single pulse spectrum.')
1848 cfg.add('eodPulseFadeFraction', 100*fade_frac, '%', 'Fraction of time of the EOD waveform snippet that is used to fade in and out to zero baseline.')
1849 cfg.add('ipiCVThresh', ipi_cv_thresh, '', 'If coefficient of variation of interpulse intervals is smaller than this threshold, then use all intervals for computing EOD frequency.')
1850 cfg.add('ipiPercentile', ipi_percentile, '%', 'Use only interpulse intervals shorter than this percentile to compute EOD frequency.')
1853def eod_waveform_args(cfg):
1854 """Translates a configuration to the respective parameter names of
1855 the function `eod_waveform()`.
1857 The return value can then be passed as key-word arguments to this
1858 function.
1860 Parameters
1861 ----------
1862 cfg: ConfigFile
1863 The configuration.
1865 Returns
1866 -------
1867 a: dict
1868 Dictionary with names of arguments of the `eod_waveform()` function
1869 and their values as supplied by `cfg`.
1870 """
1871 a = cfg.map({'win_fac': 'eodSnippetFac',
1872 'min_win': 'eodMinSnippet',
1873 'max_eods': 'eodMaxEODs',
1874 'min_sem': 'eodMinSem'})
1875 return a
1878def analyze_wave_args(cfg):
1879 """Translates a configuration to the respective parameter names of
1880 the function `analyze_wave()`.
1882 The return value can then be passed as key-word arguments to this
1883 function.
1885 Parameters
1886 ----------
1887 cfg: ConfigFile
1888 The configuration.
1890 Returns
1891 -------
1892 a: dict
1893 Dictionary with names of arguments of the `analyze_wave()` function
1894 and their values as supplied by `cfg`.
1895 """
1896 a = cfg.map({'n_harm': 'eodHarmonics',
1897 'power_n_harmonics': 'powerNHarmonics',
1898 'flip_wave': 'flipWaveEOD'})
1899 return a
1902def analyze_pulse_args(cfg):
1903 """Translates a configuration to the respective parameter names of
1904 the function `analyze_pulse()`.
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_pulse()` function
1918 and their values as supplied by `cfg`.
1919 """
1920 a = cfg.map({'min_pulse_win': 'eodMinPulseSnippet',
1921 'start_end_thresh_fac': 'eodStartEndThresholdFactor',
1922 'peak_thresh_fac': 'eodPeakThresholdFactor',
1923 'min_dist': 'eodMinimumDistance',
1924 'width_frac': 'eodPulseWidthFraction',
1925 'fit_frac': 'eodExponentialFitFraction',
1926 'flip_pulse': 'flipPulseEOD',
1927 'fit_gaussians': 'eodFitGaussians',
1928 'freq_resolution': 'eodPulseFrequencyResolution',
1929 'fade_frac': 'eodPulseFadeFraction',
1930 'ipi_cv_thresh': 'ipiCVThresh',
1931 'ipi_percentile': 'ipiPercentile'})
1932 a['width_frac'] *= 0.01
1933 a['fit_frac'] *= 0.01
1934 a['fade_frac'] *= 0.01
1935 return a
1938def add_species_config(cfg, species_file='none', wave_max_rms=0.2,
1939 pulse_max_rms=0.2):
1940 """Add parameters needed for assigning EOD waveforms to species.
1942 Parameters
1943 ----------
1944 cfg: ConfigFile
1945 The configuration.
1946 species_file: str
1947 File path to a file containing species names and corresponding
1948 file names of EOD waveform templates. If 'none', no species
1949 assignemnt is performed.
1950 wave_max_rms: float
1951 Maximum allowed rms difference (relative to standard deviation
1952 of EOD waveform) to an EOD waveform template for assignment to
1953 a wave fish species.
1954 pulse_max_rms: float
1955 Maximum allowed rms difference (relative to standard deviation
1956 of EOD waveform) to an EOD waveform template for assignment to
1957 a pulse fish species.
1958 """
1959 cfg.add_section('Species assignment:')
1960 cfg.add('speciesFile', species_file, '', 'File path to a file containing species names and corresponding file names of EOD waveform templates.')
1961 cfg.add('maximumWaveSpeciesRMS', wave_max_rms, '', 'Maximum allowed rms difference (relative to standard deviation of EOD waveform) to an EOD waveform template for assignment to a wave fish species.')
1962 cfg.add('maximumPulseSpeciesRMS', pulse_max_rms, '', 'Maximum allowed rms difference (relative to standard deviation of EOD waveform) to an EOD waveform template for assignment to a pulse fish species.')
1965def add_eod_quality_config(cfg, max_clipped_frac=0.1, max_variance=0.0,
1966 max_rms_error=0.05, min_power=-100.0, max_thd=0.0,
1967 max_crossings=4, max_relampl_harm1=0.0,
1968 max_relampl_harm2=0.0, max_relampl_harm3=0.0):
1969 """Add parameters needed for assesing the quality of an EOD waveform.
1971 Parameters
1972 ----------
1973 cfg: ConfigFile
1974 The configuration.
1976 See `wave_quality()` and `pulse_quality()` for details on
1977 the remaining arguments.
1978 """
1979 cfg.add_section('Waveform selection:')
1980 cfg.add('maximumClippedFraction', 100*max_clipped_frac, '%', 'Take waveform of the fish with the highest power only if the fraction of clipped signals is below this value.')
1981 cfg.add('maximumVariance', max_variance, '', 'Skip waveform of fish if the standard error of the EOD waveform relative to the peak-to-peak amplitude is larger than this number. A value of zero allows any variance.')
1982 cfg.add('maximumRMSError', max_rms_error, '', 'Skip waveform of wave fish if the root-mean-squared error of the fit relative to the peak-to-peak amplitude is larger than this number.')
1983 cfg.add('minimumPower', min_power, 'dB', 'Skip waveform of wave fish if its power is smaller than this value.')
1984 cfg.add('maximumTotalHarmonicDistortion', max_thd, '', 'Skip waveform of wave fish if its total harmonic distortion is larger than this value. If set to zero do not check.')
1985 cfg.add('maximumCrossings', max_crossings, '', 'Maximum number of zero crossings per EOD period.')
1986 cfg.add('maximumFirstHarmonicAmplitude', max_relampl_harm1, '', 'Skip waveform of wave fish if the amplitude of the first harmonic is higher than this factor times the amplitude of the fundamental. If set to zero do not check.')
1987 cfg.add('maximumSecondHarmonicAmplitude', max_relampl_harm2, '', 'Skip waveform of wave fish if the ampltude of the second harmonic is higher than this factor times the amplitude of the fundamental. That is, the waveform appears to have twice the frequency than the fundamental. If set to zero do not check.')
1988 cfg.add('maximumThirdHarmonicAmplitude', max_relampl_harm3, '', 'Skip waveform of wave fish if the ampltude of the third harmonic is higher than this factor times the amplitude of the fundamental. If set to zero do not check.')
1991def wave_quality_args(cfg):
1992 """Translates a configuration to the respective parameter names of
1993 the function `wave_quality()`.
1995 The return value can then be passed as key-word arguments to this
1996 function.
1998 Parameters
1999 ----------
2000 cfg: ConfigFile
2001 The configuration.
2003 Returns
2004 -------
2005 a: dict
2006 Dictionary with names of arguments of the `wave_quality()` function
2007 and their values as supplied by `cfg`.
2008 """
2009 a = cfg.map({'max_clipped_frac': 'maximumClippedFraction',
2010 'max_rms_sem': 'maximumVariance',
2011 'max_rms_error': 'maximumRMSError',
2012 'min_power': 'minimumPower',
2013 'max_crossings': 'maximumCrossings',
2014 'min_freq': 'minimumFrequency',
2015 'max_freq': 'maximumFrequency',
2016 'max_thd': 'maximumTotalHarmonicDistortion',
2017 'max_db_diff': 'maximumPowerDifference',
2018 'max_harmonics_db': 'maximumHarmonicsPower',
2019 'max_relampl_harm1': 'maximumFirstHarmonicAmplitude',
2020 'max_relampl_harm2': 'maximumSecondHarmonicAmplitude',
2021 'max_relampl_harm3': 'maximumThirdHarmonicAmplitude'})
2022 a['max_clipped_frac'] *= 0.01
2023 return a
2026def pulse_quality_args(cfg):
2027 """Translates a configuration to the respective parameter names of
2028 the function `pulse_quality()`.
2030 The return value can then be passed as key-word arguments to this
2031 function.
2033 Parameters
2034 ----------
2035 cfg: ConfigFile
2036 The configuration.
2038 Returns
2039 -------
2040 a: dict
2041 Dictionary with names of arguments of the `pulse_quality()` function
2042 and their values as supplied by `cfg`.
2043 """
2044 a = cfg.map({'max_clipped_frac': 'maximumClippedFraction',
2045 'max_rms_sem': 'maximumRMSNoise'})
2046 a['max_clipped_frac'] *= 0.01
2047 return a
2050def main():
2051 import matplotlib.pyplot as plt
2052 from .fakefish import pulsefish_eods
2053 from .pulseanalysis import analyze_pulse
2055 print('Analysis of EOD waveforms.')
2057 # data:
2058 rate = 96_000
2059 data = pulsefish_eods('Triphasic', 83.0, rate, 5.0, noise_std=0.02)
2060 unit = 'mV'
2061 eod_idx, _ = detect_peaks(data, 1.0)
2062 eod_times = eod_idx/rate
2064 # analyse EOD:
2065 mean_eod, eod_times = eod_waveform(data, rate, eod_times)
2066 mean_eod, props, peaks, pulses, energy = \
2067 analyze_pulse(mean_eod, None, eod_times)
2069 # plot:
2070 fig, axs = plt.subplots(1, 2)
2071 plot_eod_waveform(axs[0], mean_eod, props, peaks, unit=unit)
2072 axs[0].set_title(f'{props["type"]} fish: EODf = {props["EODf"]:.1f} Hz')
2073 plot_pulse_spectrum(axs[1], energy, props)
2074 plt.show()
2077if __name__ == '__main__':
2078 main()