Coverage for src / thunderfish / eodanalysis.py: 34%
1687 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-15 17:50 +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## Analysis of wave-type EODs
11- `waveeod_waveform()`: retrieve average EOD waveform via Fourier transform.
12- `analyze_wave()`: analyze the EOD waveform of a wave fish.
14## Analysis of pulse-type EODs
16- `pulse_spectrum()`: compute the spectrum of a single pulse-type EOD.
17- `analyze_pulse_spectrum()`: analyze the spectrum of a pulse-type EOD.
18- `analyze_pulse_phases()`: characterize all phases of a pulse-type EOD.
19- `decompose_pulse()`: decompose single pulse waveform into sum of Gaussians.
20- `analyze_pulse_tail()`: fit exponential to last peak/trough of pulse EOD.
21- `analyze_pulse_intervals()`: basic statistics of interpulse intervals.
22- `analyze_pulse()`: analyze the EOD waveform of a pulse fish.
24## Similarity of EOD waveforms
26- `wave_similarity()`: root-mean squared difference between two wave fish EODs.
27- `pulse_similarity()`: root-mean squared difference between two pulse fish EODs.
28- `load_species_waveforms()`: load template EOD waveforms for species matching.
30## Quality assessment
32- `clipped_fraction()`: compute fraction of clipped EOD waveform snippets.
33- `wave_quality()`: asses quality of EOD waveform of a wave fish.
34- `pulse_quality()`: asses quality of EOD waveform of a pulse fish.
36## Visualization
38- `plot_eod_recording()`: plot a zoomed in range of the recorded trace.
39- `plot_pulse_eods()`: mark pulse EODs in a plot of an EOD recording.
40- `plot_eod_snippets()`: plot a few EOD waveform snippets.
41- `plot_eod_waveform()`: plot and annotate the averaged EOD-waveform with standard error.
42- `plot_wave_spectrum()`: plot and annotate spectrum of wave EODs.
43- `plot_pulse_spectrum()`: plot and annotate spectrum of single pulse EOD.
45## Storage
47- `save_eod_waveform()`: save mean EOD waveform to file.
48- `load_eod_waveform()`: load EOD waveform from file.
49- `save_wave_eodfs()`: save frequencies of wave EODs to file.
50- `load_wave_eodfs()`: load frequencies of wave EODs from file.
51- `save_wave_fish()`: save properties of wave EODs to file.
52- `load_wave_fish()`: load properties of wave EODs from file.
53- `save_pulse_fish()`: save properties of pulse EODs to file.
54- `load_pulse_fish()`: load properties of pulse EODs from file.
55- `save_wave_spectrum()`: save amplitude and phase spectrum of wave EOD to file.
56- `load_wave_spectrum()`: load amplitude and phase spectrum of wave EOD from file.
57- `save_pulse_spectrum()`: save spectrum of pulse EOD to file.
58- `load_pulse_spectrum()`: load spectrum of pulse EOD from file.
59- `save_pulse_phases()`: save phase properties of pulse EOD to file.
60- `load_pulse_phases()`: load phase properties of pulse EOD from file.
61- `save_pulse_gaussians()`: save Gaussian phase properties of pulse EOD to file.
62- `load_pulse_gaussians()`: load Gaussian phase properties of pulse EOD from file.
63- `save_pulse_times()`: save times of pulse EOD to file.
64- `load_pulse_times()`: load times of pulse EOD from file.
65- `parse_filename()`: parse components of an EOD analysis file name.
66- `save_analysis(): save EOD analysis results to files.
67- `load_analysis()`: load EOD analysis files.
68- `load_recording()`: load recording.
70## Fit functions
72- `fourier_series()`: Fourier series of sine waves with amplitudes and phases.
73- `exp_decay()`: exponential decay.
74- `gaussian_sum()`: sum of Gaussians.
75- `gaussian_sum_spectrum()`: energy spectrum of sum of Gaussians.
76- `gaussian_sum_costs()`: cost function for fitting sum of Gaussians.
78## Filter functions
80- `unfilter()`: apply inverse low-pass filter on data.
82## Configuration
84- `add_eod_analysis_config()`: add parameters for EOD analysis functions to configuration.
85- `eod_waveform_args()`: retrieve parameters for `eod_waveform()` from configuration.
86- `analyze_wave_args()`: retrieve parameters for `analyze_wave()` from configuration.
87- `analyze_pulse_args()`: retrieve parameters for `analyze_pulse()` from configuration.
88- `add_species_config()`: add parameters needed for assigning EOD waveforms to species.
89- `add_eod_quality_config()`: add parameters needed for assesing the quality of an EOD waveform.
90- `wave_quality_args()`: retrieve parameters for `wave_quality()` from configuration.
91- `pulse_quality_args()`: retrieve parameters for `pulse_quality()` from configuration.
92"""
94import io
95import glob
96import zipfile
97import numpy as np
98try:
99 import matplotlib.pyplot as plt
100except ImportError:
101 pass
103from pathlib import Path
104from scipy.optimize import curve_fit, minimize
105from numba import jit
106from thunderlab.eventdetection import percentile_threshold, detect_peaks, snippets, peak_width
107from thunderlab.eventdetection import threshold_crossings, threshold_crossing_times, merge_events
108from thunderlab.powerspectrum import next_power_of_two, nfft, decibel
109from thunderlab.tabledata import TableData
110from thunderlab.dataloader import DataLoader
112from .harmonics import fundamental_freqs_and_power
113from .fakefish import pulsefish_waveform
114from .fakefish import normalize_pulsefish, export_pulsefish
115from .fakefish import normalize_wavefish, export_wavefish
118def eod_waveform(data, rate, eod_times, win_fac=2.0, min_win=0.01,
119 min_sem=False, max_eods=None, unfilter_cutoff=0.0):
120 """Extract data snippets around each EOD, and compute a mean waveform with standard error.
122 Retrieving the EOD waveform of a wave fish works under the following
123 conditions: (i) at a signal-to-noise ratio \\(SNR = P_s/P_n\\),
124 i.e. the power \\(P_s\\) of the EOD of interest relative to the
125 largest other EOD \\(P_n\\), we need to average over at least \\(n >
126 (SNR \\cdot c_s^2)^{-1}\\) snippets to bring the standard error of the
127 averaged EOD waveform down to \\(c_s\\) relative to its
128 amplitude. For a s.e.m. less than 5% ( \\(c_s=0.05\\) ) and an SNR of
129 -10dB (the signal is 10 times smaller than the noise, \\(SNR=0.1\\) ) we
130 get \\(n > 0.00025^{-1} = 4000\\) data snippets - a recording a
131 couple of seconds long. (ii) Very important for wave fish is that
132 they keep their frequency constant. Slight changes in the EOD
133 frequency will corrupt the average waveform. If the period of the
134 waveform changes by \\(c_f=\\Delta T/T\\), then after \\(n =
135 1/c_f\\) periods moved the modified waveform through a whole period.
136 This is in the range of hundreds or thousands waveforms.
138 NOTE: we need to take into account a possible error in the estimate
139 of the EOD period. This will limit the maximum number of snippets to
140 be averaged.
142 If `min_sem` is set, the algorithm checks for a global minimum of
143 the s.e.m. as a function of snippet number. If there is one then
144 the average is computed for this number of snippets, otherwise all
145 snippets are taken from the provided data segment. Note that this
146 check only works for the strongest EOD in a recording. For weaker
147 EOD the s.e.m. always decays with snippet number (empirical
148 observation).
150 TODO: use power spectra to check for changes in EOD frequency!
152 Parameters
153 ----------
154 data: 1-D array of float
155 The data to be analysed.
156 rate: float
157 Sampling rate of the data in Hertz.
158 eod_times: 1-D array of float
159 Array of EOD times in seconds over which the waveform should be
160 averaged.
161 WARNING: The first data point must be at time zero!
162 win_fac: float
163 The snippet size is the EOD period times `win_fac`. The EOD period
164 is determined as the minimum interval between EOD times.
165 min_win: float
166 The minimum size of the snippets in seconds.
167 min_sem: bool
168 If set, check for minimum in s.e.m. to set the maximum numbers
169 of EODs to be used for computing the average waveform.
170 max_eods: int or None
171 Maximum number of EODs to be used for averaging.
172 unfilter_cutoff: float
173 If not zero, the cutoff frequency for an inverse high-pass filter
174 applied to the mean EOD waveform.
176 Returns
177 -------
178 mean_eod: 2-D array
179 Average of the EOD snippets. First column is time in seconds,
180 second column the mean eod, third column the standard error.
181 eod_times: 1-D array
182 Times of EOD peaks in seconds that have been actually used to calculate the
183 averaged EOD waveform.
184 """
185 # indices of EOD times:
186 eod_idx = np.round(eod_times*rate).astype(int)
188 # window size:
189 period = np.min(np.diff(eod_times))
190 win = 0.5*win_fac*period
191 if 2*win < min_win:
192 win = 0.5*min_win
193 win_inx = int(win*rate)
195 # extract snippets:
196 eod_times = eod_times[(eod_idx >= win_inx) & (eod_idx < len(data)-win_inx)]
197 eod_idx = eod_idx[(eod_idx >= win_inx) & (eod_idx < len(data)-win_inx)]
198 if max_eods and max_eods > 0 and len(eod_idx) > max_eods:
199 dn = (len(eod_idx) - max_eods)//2
200 eod_times = eod_times[dn:dn+max_eods]
201 eod_idx = eod_idx[dn:dn+max_eods]
202 eod_snippets = snippets(data, eod_idx, -win_inx, win_inx)
203 if len(eod_snippets) == 0:
204 return np.zeros((0, 3)), eod_times
206 # optimal number of snippets:
207 step = 10
208 if min_sem and len(eod_snippets) > step:
209 sems = [np.mean(np.std(eod_snippets[:k], axis=0, ddof=1)/np.sqrt(k))
210 for k in range(step, len(eod_snippets), step)]
211 idx = np.argmin(sems)
212 # there is a local minimum:
213 if idx > 0 and idx < len(sems)-1:
214 maxn = step*(idx+1)
215 eod_snippets = eod_snippets[:maxn]
216 eod_times = eod_times[:maxn]
218 # mean and std of snippets:
219 mean_eod = np.zeros((len(eod_snippets[0]), 3))
220 mean_eod[:, 1] = np.mean(eod_snippets, axis=0)
221 if len(eod_snippets) > 1:
222 mean_eod[:, 2] = np.std(eod_snippets, axis=0, ddof=1)/np.sqrt(len(eod_snippets))
224 # apply inverse filter:
225 if unfilter_cutoff and unfilter_cutoff > 0.0:
226 unfilter(mean_eod[:, 1], rate, unfilter_cutoff)
228 # time axis:
229 mean_eod[:, 0] = (np.arange(len(mean_eod)) - win_inx) / rate
231 return mean_eod, eod_times
234def adjust_eodf(eodf, temp, temp_adjust=25.0, q10=1.62):
235 """Adjust EOD frequencies to a standard temperature using Q10.
237 Parameters
238 ----------
239 eodf: float or ndarray
240 EOD frequencies.
241 temp: float
242 Temperature in degree celsisus at which EOD frequencies in
243 `eodf` were measured.
244 temp_adjust: float
245 Standard temperature in degree celsisus to which EOD
246 frequencies are adjusted.
247 q10: float
248 Q10 value describing temperature dependence of EOD
249 frequencies. The default of 1.62 is from Dunlap, Smith, Yetka
250 (2000) Brain Behav Evol, measured for Apteronotus
251 lepthorhynchus in the lab.
253 Returns
254 -------
255 eodf_corrected: float or array
256 EOD frequencies adjusted to `temp_adjust` using `q10`.
257 """
258 return eodf*q10**((temp_adjust - temp) / 10.0)
261def waveeod_waveform(data, rate, freq, win_fac=2.0, unfilter_cutoff=0.0):
262 """Retrieve average EOD waveform via Fourier transform.
264 TODO: use power spectra to check minimum data segment needed and
265 check for changes in frequency over several segments!
267 TODO: return waveform with higher samplign rate? (important for
268 2kHz waves on 24kHz sampling). But this seems to render some EODs
269 inacceptable in the further thunderfish processing pipeline.
271 Parameters
272 ----------
273 data: 1-D array of float
274 The data to be analysed.
275 rate: float
276 Sampling rate of the data in Hertz.
277 freq: float
278 EOD frequency.
279 win_fac: float
280 The snippet size is the EOD period times `win_fac`. The EOD period
281 is determined as the minimum interval between EOD times.
282 unfilter_cutoff: float
283 If not zero, the cutoff frequency for an inverse high-pass filter
284 applied to the mean EOD waveform.
286 Returns
287 -------
288 mean_eod: 2-D array
289 Average of the EOD snippets. First column is time in seconds,
290 second column the mean eod, third column the standard error.
291 eod_times: 1-D array
292 Times of EOD peaks in seconds that have been actually used to
293 calculate the averaged EOD waveform.
295 """
297 @jit(nopython=True)
298 def fourier_wave(data, rate, freq):
299 """
300 extracting wave via fourier coefficients
301 """
302 twave = np.arange(0, (1+win_fac)/freq, 1/rate)
303 wave = np.zeros(len(twave))
304 t = np.arange(len(data))/rate
305 for k in range(0, 31):
306 Xk = np.trapz(data*np.exp(-1j*2*np.pi*k*freq*t), t)*2/t[-1]
307 wave += np.real(Xk*np.exp(1j*2*np.pi*k*freq*twave))
308 return wave
310 @jit(nopython=True)
311 def fourier_range(data, rate, f0, f1, df):
312 wave = np.zeros(1)
313 freq = f0
314 for f in np.arange(f0, f1, df):
315 w = fourier_wave(data, rate, f)
316 if np.max(w) - np.min(w) > np.max(wave) - np.min(wave):
317 wave = w
318 freq = f
319 return wave, freq
321 # TODO: parameterize!
322 tsnippet = 2
323 min_corr = 0.98
324 min_ampl_frac = 0.5
325 frange = 0.1
326 fstep = 0.1
327 waves = []
328 freqs = []
329 times = []
330 step = int(tsnippet*rate)
331 for i in range(0, len(data) - step//2, step//2):
332 w, f = fourier_range(data[i:i + step], rate, freq - frange,
333 freq + frange + fstep/2, fstep)
334 waves.append(w)
335 freqs.append(f)
336 """
337 waves.append(np.zeros(1))
338 freqs.append(freq)
339 for f in np.arange(freq - frange, freq + frange + fstep/2, fstep):
340 w = fourier_wave(data[i:i + step], rate, f)
341 if np.max(w) - np.min(w) > np.max(waves[-1]) - np.min(waves[-1]):
342 waves[-1] = w
343 freqs[-1] = f
344 """
345 times.append(np.arange(i/rate, (i + step)/rate, 1/freqs[-1]))
346 eod_freq = np.mean(freqs)
347 mean_eod = np.zeros((0, 3))
348 eod_times = np.zeros((0))
349 if len(waves) == 0:
350 return mean_eod, eod_times
351 for k in range(len(waves)):
352 period = int(np.ceil(rate/freqs[k]))
353 i = np.argmax(waves[k][:period])
354 waves[k] = waves[k][i:]
355 n = np.min([len(w) for w in waves])
356 waves = np.array([w[:n] for w in waves])
357 # only snippets that are similar:
358 if len(waves) > 1:
359 corr = np.corrcoef(waves)
360 nmax = np.argmax(np.sum(corr > min_corr, axis=1))
361 if nmax <= 1:
362 nmax = 2
363 select = np.sum(corr > min_corr, axis=1) >= nmax
364 waves = waves[select]
365 times = [times[k] for k in range(len(times)) if select[k]]
366 if len(waves) == 0:
367 return mean_eod, eod_times
368 # only the largest snippets:
369 ampls = np.std(waves, axis=1)
370 select = ampls >= min_ampl_frac*np.max(ampls)
371 waves = waves[select]
372 times = [times[k] for k in range(len(times)) if select[k]]
373 if len(waves) == 0:
374 return mean_eod, eod_times
375 """
376 #plt.plot(freqs)
377 plt.plot(waves.T)
378 plt.show()
379 """
380 mean_eod = np.zeros((n, 3))
381 mean_eod[:, 0] = np.arange(len(mean_eod))/rate
382 mean_eod[:, 1] = np.mean(waves, axis=0)
383 mean_eod[:, 2] = np.std(waves, axis=0)
384 eod_times = np.concatenate(times)
386 # apply inverse filter:
387 if unfilter_cutoff and unfilter_cutoff > 0.0:
388 unfilter(mean_eod[:, 1], rate, unfilter_cutoff)
390 return mean_eod, eod_times
393def unfilter(data, rate, cutoff):
394 """Apply inverse high-pass filter on data.
396 Assumes high-pass filter
397 \\[ \\tau \\dot y = -y + \\tau \\dot x \\]
398 has been applied on the original data \\(x\\), where
399 \\(\\tau=(2\\pi f_{cutoff})^{-1}\\) is the time constant of the
400 filter. To recover \\(x\\) the ODE
401 \\[ \\tau \\dot x = y + \\tau \\dot y \\]
402 is applied on the filtered data \\(y\\).
404 Parameters
405 ----------
406 data: ndarray
407 High-pass filtered original data.
408 rate: float
409 Sampling rate of `data` in Hertz.
410 cutoff: float
411 Cutoff frequency \\(f_{cutoff}\\) of the high-pass filter in Hertz.
413 Returns
414 -------
415 data: ndarray
416 Recovered original data.
417 """
418 tau = 0.5/np.pi/cutoff
419 fac = tau*rate
420 data -= np.mean(data)
421 d0 = data[0]
422 x = d0
423 for k in range(len(data)):
424 d1 = data[k]
425 x += (d1 - d0) + d0/fac
426 data[k] = x
427 d0 = d1
428 return data
431def fourier_series(t, freq, *ap):
432 """Fourier series of sine waves with amplitudes and phases.
434 x(t) = sum_{i=0}^n ap[2*i]*sin(2 pi (i+1) freq t + ap[2*i+1])
436 Parameters
437 ----------
438 t: float or array
439 Time.
440 freq: float
441 Fundamental frequency.
442 *ap: list of floats
443 The amplitudes and phases (in rad) of the fundamental and harmonics.
445 Returns
446 -------
447 x: float or array
448 The Fourier series evaluated at times `t`.
449 """
450 omega = 2.0*np.pi*freq
451 x = 0.0
452 for i, (a, p) in enumerate(zip(ap[0:-1:2], ap[1::2])):
453 x += a*np.sin((i+1)*omega*t+p)
454 return x
457def analyze_wave(eod, freq, n_harm=10, power_n_harmonics=0,
458 n_harmonics=3, flip_wave='none'):
459 """Analyze the EOD waveform of a wave fish.
461 Parameters
462 ----------
463 eod: 2-D array
464 The eod waveform. First column is time in seconds, second
465 column the EOD waveform, third column, if present, is the
466 standard error of the EOD waveform, Further columns are
467 optional but not used.
468 freq: float or 2-D array
469 The frequency of the EOD or the list of harmonics (rows) with
470 frequency and peak height (columns) as returned from
471 `harmonics.harmonic_groups()`.
472 n_harm: int
473 Maximum number of harmonics used for the Fourier decomposition.
474 power_n_harmonics: int
475 Sum over the first `power_n_harmonics` harmonics for computing
476 the total power. If 0 sum over all harmonics.
477 n_harmonics: int
478 The maximum power of higher harmonics is computed from
479 harmonics higher than the maximum harmonics within the first
480 three harmonics plus `n_harmonics`.
481 flip_wave: 'auto', 'none', 'flip'
482 - 'auto' flip waveform such that the larger extremum is positive.
483 - 'flip' flip waveform.
484 - 'none' do not flip waveform.
486 Returns
487 -------
488 meod: 2-D array of floats
489 The eod waveform. First column is time in seconds, second
490 column the eod waveform. Further columns are kept from the
491 input `eod`. And a column is added with the fit of the fourier
492 series to the waveform.
493 props: dict
494 A dictionary with properties of the analyzed EOD waveform.
496 - type: set to 'wave'.
497 - EODf: is set to the EOD fundamental frequency.
498 - p-p-amplitude: peak-to-peak amplitude of the Fourier fit.
499 - flipped: True if the waveform was flipped.
500 - amplitude: amplitude factor of the Fourier fit.
501 - noise: root-mean squared standard error mean of the averaged
502 EOD waveform relative to the p-p amplitude.
503 - rmserror: root-mean-square error between Fourier-fit and
504 EOD waveform relative to the p-p amplitude. If larger than
505 about 0.05 the data are bad.
506 - ncrossings: number of zero crossings per period
507 - peakwidth: width of the peak at the averaged amplitude relative
508 to EOD period.
509 - troughwidth: width of the trough at the averaged amplitude
510 relative to EOD period.
511 - minwidth: peakwidth or troughwidth, whichever is smaller.
512 - leftpeak: time from positive zero crossing to peak relative
513 to EOD period.
514 - rightpeak: time from peak to negative zero crossing relative to
515 EOD period.
516 - lefttrough: time from negative zero crossing to trough relative to
517 EOD period.
518 - righttrough: time from trough to positive zero crossing relative to
519 EOD period.
520 - p-p-distance: time between peak and trough relative to EOD period.
521 - min-p-p-distance: p-p-distance or EOD period minus p-p-distance,
522 whichever is smaller, relative to EOD period.
523 - relpeakampl: amplitude of peak or trough, whichever is larger,
524 relative to p-p amplitude.
525 - power: summed power of all harmonics of the extracted EOD waveform
526 in decibel relative to one.
527 - datapower: summed power of all harmonics of the original data in
528 decibel relative to one. Only if `freq` is a list of harmonics.
529 - thd: total harmonic distortion, i.e. square root of sum of
530 amplitudes squared of harmonics relative to amplitude
531 of fundamental.
532 - dbdiff: smoothness of power spectrum as standard deviation of
533 differences in decibel power.
534 - maxdb: maximum power of higher harmonics relative to peak power
535 in decibel.
537 spec_data: 2-D array of floats
538 First six columns are from the spectrum of the extracted
539 waveform. First column is the index of the harmonics, second
540 column its frequency, third column its amplitude, fourth
541 column its amplitude relative to the fundamental, fifth column
542 is power of harmonics relative to fundamental in decibel, and
543 sixth column the phase shift relative to the fundamental.
544 If `freq` is a list of harmonics, a seventh column is added to
545 `spec_data` that contains the powers of the harmonics from the
546 original power spectrum of the raw data. Rows are the
547 harmonics, first row is the fundamental frequency with index
548 0, relative amplitude of one, relative power of 0dB, and phase
549 shift of zero.
550 error_str: string
551 If fitting of the fourier series failed,
552 this is reported in this string.
554 Raises
555 ------
556 IndexError:
557 EOD data is less than one period long.
558 """
559 error_str = ''
561 freq0 = freq
562 if hasattr(freq, 'shape'):
563 freq0 = freq[0][0]
565 # storage:
566 meod = np.zeros((eod.shape[0], eod.shape[1]+1))
567 meod[:, :-1] = eod
569 # subtract mean and flip:
570 period = 1.0/freq0
571 pinx = int(np.ceil(period/(meod[1,0]-meod[0,0])))
572 maxn = (len(meod)//pinx)*pinx
573 if maxn < pinx: maxn = len(meod)
574 offs = (len(meod) - maxn)//2
575 meod[:, 1] -= np.mean(meod[offs:offs+pinx,1])
576 flipped = False
577 if 'flip' in flip_wave or ('auto' in flip_wave and -np.min(meod[:, 1]) > np.max(meod[:, 1])):
578 meod[:, 1] = -meod[:, 1]
579 flipped = True
581 # move peak of waveform to zero:
582 offs = len(meod)//4
583 maxinx = offs+np.argmax(meod[offs:3*offs,1])
584 meod[:, 0] -= meod[maxinx,0]
586 # indices of exactly one or two periods around peak:
587 if len(meod) < pinx:
588 raise IndexError('data need to contain at least one EOD period')
589 if len(meod) >= 2*pinx:
590 i0 = maxinx - pinx if maxinx >= pinx else 0
591 i1 = i0 + 2*pinx
592 if i1 > len(meod):
593 i1 = len(meod)
594 i0 = i1 - 2*pinx
595 else:
596 i0 = maxinx - pinx//2 if maxinx >= pinx//2 else 0
597 i1 = i0 + pinx
599 # subtract mean:
600 meod[:, 1] -= np.mean(meod[i0:i1,1])
602 # zero crossings:
603 ui, di = threshold_crossings(meod[:, 1], 0.0)
604 ut, dt = threshold_crossing_times(meod[:, 0], meod[:, 1], 0.0, ui, di)
605 ut, dt = merge_events(ut, dt, 0.02/freq0)
606 ncrossings = int(np.round((len(ut) + len(dt))/(meod[-1,0]-meod[0,0])/freq0))
607 if np.any(ut<0.0):
608 up_time = ut[ut<0.0][-1]
609 else:
610 up_time = 0.0
611 error_str += '%.1f Hz wave fish: no upward zero crossing. ' % freq0
612 if np.any(dt>0.0):
613 down_time = dt[dt>0.0][0]
614 else:
615 down_time = 0.0
616 error_str += '%.1f Hz wave fish: no downward zero crossing. ' % freq0
617 peak_width = down_time - up_time
618 trough_width = period - peak_width
619 peak_time = 0.0
620 trough_time = meod[maxinx+np.argmin(meod[maxinx:maxinx+pinx,1]),0]
621 phase1 = peak_time - up_time
622 phase2 = down_time - peak_time
623 phase3 = trough_time - down_time
624 phase4 = up_time + period - trough_time
625 distance = trough_time - peak_time
626 min_distance = distance
627 if distance > period/2:
628 min_distance = period - distance
630 # fit fourier series:
631 ampl = 0.5*(np.max(meod[:, 1])-np.min(meod[:, 1]))
632 while n_harm > 1:
633 params = [freq0]
634 for i in range(1, n_harm+1):
635 params.extend([ampl/i, 0.0])
636 try:
637 popt, pcov = curve_fit(fourier_series, meod[i0:i1,0],
638 meod[i0:i1,1], params, maxfev=2000)
639 break
640 except (RuntimeError, TypeError):
641 error_str += '%.1f Hz wave fish: fit of fourier series failed for %d harmonics. ' % (freq0, n_harm)
642 n_harm //= 2
643 # store fourier fit:
644 meod[:, -1] = fourier_series(meod[:, 0], *popt)
645 # make all amplitudes positive:
646 for i in range(n_harm):
647 if popt[i*2+1] < 0.0:
648 popt[i*2+1] *= -1.0
649 popt[i*2+2] += np.pi
650 # phases relative to fundamental:
651 # phi0 = 2*pi*f0*dt <=> dt = phi0/(2*pi*f0)
652 # phik = 2*pi*i*f0*dt = i*phi0
653 phi0 = popt[2]
654 for i in range(n_harm):
655 popt[i*2+2] -= (i + 1)*phi0
656 # all phases in the range -pi to pi:
657 popt[i*2+2] %= 2*np.pi
658 if popt[i*2+2] > np.pi:
659 popt[i*2+2] -= 2*np.pi
660 # store fourier spectrum:
661 if hasattr(freq, 'shape'):
662 n = n_harm
663 n += np.sum(freq[:, 0] > (n_harm+0.5)*freq[0,0])
664 spec_data = np.zeros((n, 7))
665 spec_data[:, :] = np.nan
666 k = 0
667 for i in range(n_harm):
668 while k < len(freq) and freq[k,0] < (i+0.5)*freq0:
669 k += 1
670 if k >= len(freq):
671 break
672 if freq[k,0] < (i+1.5)*freq0:
673 spec_data[i,6] = freq[k,1]
674 k += 1
675 for i in range(n_harm, n):
676 if k >= len(freq):
677 break
678 spec_data[i,:2] = [np.round(freq[k,0]/freq0)-1, freq[k,0]]
679 spec_data[i,6] = freq[k,1]
680 k += 1
681 else:
682 spec_data = np.zeros((n_harm, 6))
683 for i in range(n_harm):
684 spec_data[i,:6] = [i, (i+1)*freq0, popt[i*2+1], popt[i*2+1]/popt[1],
685 decibel((popt[i*2+1]/popt[1])**2.0), popt[i*2+2]]
686 # smoothness of power spectrum:
687 db_powers = decibel(spec_data[:n_harm,2]**2)
688 db_diff = np.std(np.diff(db_powers))
689 # maximum relative power of higher harmonics:
690 p_max = np.argmax(db_powers[:3])
691 db_powers -= db_powers[p_max]
692 if len(db_powers[p_max+n_harmonics:]) == 0:
693 max_harmonics_power = -100.0
694 else:
695 max_harmonics_power = np.max(db_powers[p_max+n_harmonics:])
696 # total harmonic distortion:
697 thd = np.sqrt(np.nansum(spec_data[1:, 3]))
699 # peak-to-peak and trough amplitudes:
700 ppampl = np.max(meod[i0:i1,1]) - np.min(meod[i0:i1,1])
701 relpeakampl = max(np.max(meod[i0:i1,1]), np.abs(np.min(meod[i0:i1,1])))/ppampl
703 # variance and fit error:
704 rmssem = np.sqrt(np.mean(meod[i0:i1,2]**2.0))/ppampl if eod.shape[1] > 2 else None
705 rmserror = np.sqrt(np.mean((meod[i0:i1,1] - meod[i0:i1,-1])**2.0))/ppampl
707 # store results:
708 props = {}
709 props['type'] = 'wave'
710 props['EODf'] = freq0
711 props['p-p-amplitude'] = ppampl
712 props['flipped'] = flipped
713 props['amplitude'] = 0.5*ppampl # remove it
714 props['rmserror'] = rmserror
715 if rmssem:
716 props['noise'] = rmssem
717 props['ncrossings'] = ncrossings
718 props['peakwidth'] = peak_width/period
719 props['troughwidth'] = trough_width/period
720 props['minwidth'] = min(peak_width, trough_width)/period
721 props['leftpeak'] = phase1/period
722 props['rightpeak'] = phase2/period
723 props['lefttrough'] = phase3/period
724 props['righttrough'] = phase4/period
725 props['p-p-distance'] = distance/period
726 props['min-p-p-distance'] = min_distance/period
727 props['relpeakampl'] = relpeakampl
728 pnh = power_n_harmonics if power_n_harmonics > 0 else n_harm
729 pnh = min(n_harm, pnh)
730 props['power'] = decibel(np.sum(spec_data[:pnh,2]**2.0))
731 if hasattr(freq, 'shape'):
732 props['datapower'] = decibel(np.sum(freq[:pnh,1]))
733 props['thd'] = thd
734 props['dbdiff'] = db_diff
735 props['maxdb'] = max_harmonics_power
737 return meod, props, spec_data, error_str
740def exp_decay(t, tau, ampl, offs):
741 """Exponential decay function.
743 x(t) = ampl*exp(-t/tau) + offs
745 Parameters
746 ----------
747 t: float or array
748 Time.
749 tau: float
750 Time constant of exponential decay.
751 ampl: float
752 Amplitude of exponential decay, i.e. initial value minus
753 steady-state value.
754 offs: float
755 Steady-state value.
757 Returns
758 -------
759 x: float or array
760 The exponential decay evaluated at times `t`.
762 """
763 return offs + ampl*np.exp(-t/tau)
766def pulse_spectrum(eod, ratetime=None, freq_resolution=1.0, fade_frac=0.0):
767 """Compute the spectrum of a single pulse-type EOD.
769 Parameters
770 ----------
771 eod: 1-D or 2-D array
772 The eod waveform of which the spectrum is computed.
773 If an 1-D array, then this is the waveform and you
774 need to also pass a sampling rate in `rate`.
775 If a 2-D array, then first column is time in seconds and second
776 column is the eod waveform. Further columns are ignored.
777 ratetime: None or float or array of float
778 If a 1-D array is passed on to `eod` then either the sampling
779 rate in Hertz or the time array corresponding to `eod`.
780 freq_resolution: float
781 The frequency resolution of the spectrum.
782 fade_frac: float
783 Fraction of time of the EOD waveform that is used to fade in
784 and out to zero baseline.
786 Returns
787 -------
788 freqs: 1-D array of float
789 The frequency components of the energy spectrum.
790 energy: 1-D array of float
791 The energy spectrum of the single pulse EOD
792 with unit (x s)^2 = x^2 s/Hz.
793 The integral over the energy spectrum `np.sum(energy)*freqs[1]`
794 equals the integral over the squared eod, `np.sum(eod**2)/rate`.
795 That is, by making the energy spectrum a power sepctrum
796 (dividing the energy by the FFT window duration), the integral
797 over the power spectrum equals the mean-squared signal
798 (variance). But the single-pulse spectrum is not a power-spectrum.
799 because in the limit to infinitely long window, the power vanishes!
801 See Also
802 --------
803 thunderfish.fakefish.pulsefish_spectrum(): analytically computed spectra for simulated pulse EODs.
805 """
806 if eod.ndim == 2:
807 rate = 1.0/(eod[1, 0] - eod[0, 0])
808 eod = eod[:, 1]
809 elif isinstance(ratetime, (list, tuple, np.ndarray)):
810 rate = 1.0/(ratetime[1] - ratetime[0])
811 else:
812 rate = ratetime
813 n_fft = nfft(rate, freq_resolution)
814 # subtract mean computed from the ends of the EOD snippet:
815 n = len(eod)//20 if len(eod) >= 20 else 1
816 eod = eod - 0.5*(np.mean(eod[:n]) + np.mean(eod[-n:]))
817 # zero padding:
818 max_idx = np.argmax(eod)
819 n0 = max_idx
820 n1 = len(eod) - max_idx
821 n = 2*max(n0, n1)
822 if n_fft < n:
823 n_fft = next_power_of_two(n)
824 data = np.zeros(n_fft)
825 data[n_fft//2 - n0:n_fft//2 + n1] = eod
826 # fade in and out:
827 if fade_frac > 0:
828 fn = int(fade_frac*len(eod))
829 data[n_fft//2 - n0:n_fft//2 - n0 + fn] *= np.arange(fn)/fn
830 data[n_fft//2 + n1 - fn:n_fft//2 + n1] *= np.arange(fn)[::-1]/fn
831 # spectrum:
832 dt = 1/rate
833 freqs = np.fft.rfftfreq(n_fft, dt)
834 fourier = np.fft.rfft(data)*dt
835 energy = 2*np.abs(fourier)**2 # one-sided spectrum!
836 return freqs, energy
839def analyze_pulse_spectrum(freqs, energy):
840 """Analyze the spectrum of a pulse-type EOD.
842 Parameters
843 ----------
844 freqs: 1-D array of float
845 The frequency components of the energy spectrum.
846 energy: 1-D array of float
847 The energy spectrum of the single pulse-type EOD.
849 Returns
850 -------
851 peak_freq: float
852 Frequency at peak energy of the spectrum in Hertz.
853 peak_energy: float
854 Peak energy of the pulse spectrum in x^2 s/Hz.
855 trough_freq: float
856 Frequency at trough before peak in Hertz.
857 trough_energy: float
858 Energy of trough before peak in x^2 s/Hz.
859 att5: float
860 Attenuation of average energy below 5 Hz relative to
861 peak energy in decibel.
862 att50: float
863 Attenuation of average energy below 50 Hz relative to
864 peak energy in decibel.
865 low_cutoff: float
866 Frequency at which the energy reached half of the
867 peak energy relative to the DC energy in Hertz.
868 high_cutoff: float
869 3dB roll-off frequency in Hertz.
870 """
871 ip = np.argmax(energy)
872 peak_freq = freqs[ip]
873 peak_energy = energy[ip]
874 it = np.argmin(energy[:ip]) if ip > 0 else 0
875 trough_freq = freqs[it]
876 trough_energy = energy[it]
877 att5 = decibel(np.mean(energy[freqs<5.0]), peak_energy)
878 att50 = decibel(np.mean(energy[freqs<50.0]), peak_energy)
879 low_cutoff = freqs[decibel(energy, peak_energy) > 0.5*att5][0]
880 high_cutoff = freqs[decibel(energy, peak_energy) > -3.0][-1]
881 return peak_freq, peak_energy, trough_freq, trough_energy, \
882 att5, att50, low_cutoff, high_cutoff
885def analyze_pulse_phases(threshold, eod, ratetime=None,
886 min_dist=50.0e-6, width_frac=0.5):
887 """ Characterize all phases of a pulse-type EOD.
889 Parameters
890 ----------
891 threshold: float
892 Threshold for detecting peaks and troughs.
893 eod: 1-D or 2-D array
894 The eod waveform to be analyzed.
895 If an 1-D array, then this is the waveform and you
896 need to also pass a sampling rate in `rate`.
897 If a 2-D array, then first column is time in seconds and second
898 column is the eod waveform. Further columns are ignored.
899 ratetime: None or float or array of float
900 If a 1-D array is passed on to `eod` then either the sampling
901 rate in Hertz or the time array corresponding to `eod`.
902 min_dist: float
903 Minimum distance between peak and troughs of the pulse.
904 width_frac: float
905 The width of an EOD phase is measured at this fraction of a peak's
906 height (0-1).
908 Returns
909 -------
910 phases: 2-D array
911 For each peak and trough (rows) of the EOD waveform
912 7 columns: the peak index (1 is P1, i.e. the largest positive peak),
913 time relative to largest positive peak, amplitude,
914 amplitude normalized to P1 phase,
915 width of peak/trough at `width_frac` height,
916 area under the peak, and area under the peak relative to total area.
917 Empty if waveform is not a pulse EOD.
918 """
919 if eod.ndim == 2:
920 time = eod[:, 0]
921 eod = eod[:, 1]
922 elif isinstance(ratetime, (list, tuple, np.ndarray)):
923 time = ratetime
924 else:
925 time = np.arange(len(eod))/rate
926 dt = time[1] - time[0]
927 # find peaks and troughs:
928 peak_idx, trough_idx = detect_peaks(eod, threshold)
929 if len(peak_idx) == 0:
930 return np.zeros((0, 7))
931 # and their width:
932 peak_widths = peak_width(time, eod, peak_idx, trough_idx,
933 peak_frac=width_frac, base='max')
934 trough_widths = peak_width(time, -eod, trough_idx, peak_idx,
935 peak_frac=width_frac, base='max')
936 # combine peaks and troughs:
937 pt_idx = np.concatenate((peak_idx, trough_idx))
938 pt_widths = np.concatenate((peak_widths, trough_widths))
939 pts_idx = np.argsort(pt_idx)
940 phase_list = pt_idx[pts_idx]
941 width_list = pt_widths[pts_idx]
942 # remove multiple peaks that are too close:
943 # TODO: XXX replace by Dexters function that keeps the maximum peak
944 rmidx = [(k, k+1) for k in np.where(np.diff(time[phase_list]) < min_dist)[0]]
945 # flatten and keep maximum and minimum phase:
946 max_idx = np.argmax(eod)
947 min_idx = np.argmin(eod)
948 rmidx = np.unique([k for kk in rmidx for k in kk
949 if phase_list[k] != max_idx and phase_list[k] != min_idx])
950 # delete:
951 if len(rmidx) > 0:
952 phase_list = np.delete(phase_list, rmidx)
953 width_list = np.delete(width_list, rmidx)
954 if len(phase_list) == 0:
955 return np.zeros((0, 7))
956 # find P1:
957 p1i = np.argmax(phase_list == max_idx)
958 ## truncate peaks to the left: # TODO: REALLY? WHY?
959 #offs = 0 if p1i <= 2 else p1i - 2
960 #phase_list = phase_list[offs:]
961 #width_list = width_list[offs:]
962 #p1i -= offs
963 # amplitudes:
964 amplitudes = eod[phase_list]
965 max_ampl = np.abs(amplitudes[p1i])
966 # phase areas:
967 phase_areas = np.zeros(len(phase_list))
968 for i in range(len(phase_list)):
969 sign_fac = np.sign(eod[phase_list[i]])
970 i0 = phase_list[i - 1] if i > 0 else 0
971 i1 = phase_list[i + 1] if i + 1 < len(phase_list) else len(eod)
972 if i0 > 0 and sign_fac*eod[i0] > 0 and \
973 i1 < len(eod) and sign_fac*eod[i1] > 0:
974 phase_areas[i] = 0
975 else:
976 snippet = sign_fac*eod[i0:i1]
977 phase_areas[i] = np.sum(snippet[snippet > 0])*dt
978 phase_areas[i] *= sign_fac
979 total_area = np.sum(np.abs(phase_areas))
980 # store phase properties:
981 phases = np.zeros((len(phase_list), 7))
982 for i, pi in enumerate(phase_list):
983 phases[i, :] = [i + 1 - p1i, time[pi], amplitudes[i],
984 amplitudes[i]/max_ampl, width_list[i],
985 phase_areas[i], phase_areas[i]/total_area]
986 return phases
989def gaussian_sum(x, *pas):
990 """Sum of Gaussians.
992 Parameters
993 ----------
994 x: array of float
995 The x array over which the sum of Gaussians is evaluated.
996 *pas: list of floats
997 The parameters of the Gaussians in a flat list.
998 Position, amplitude, and standard deviation of first Gaussian,
999 position, amplitude, and standard deviation of second Gaussian,
1000 and so on.
1002 Returns
1003 -------
1004 sg: array of float
1005 The sum of Gaussians for the times given in `t`.
1006 """
1007 sg = np.zeros(len(x))
1008 for pos, ampl, std in zip(pas[0:-2:3], pas[1:-1:3], pas[2::3]):
1009 sg += ampl*np.exp(-0.5*((x - pos)/std)**2)
1010 return sg
1013def gaussian_sum_spectrum(f, *pas):
1014 """Energy spectrum of sum of Gaussians.
1016 Parameters
1017 ----------
1018 f: 1-D array of float
1019 The frequencies at which to evaluate the spectrum.
1020 *pas: list of floats
1021 The parameters of the Gaussians in a flat list.
1022 Position, amplitude, and standard deviation of first Gaussian,
1023 position, amplitude, and standard deviation of second Gaussian,
1024 and so on.
1026 Returns
1027 -------
1028 energy: 1-D array of float
1029 The one-sided energy spectrum of the sum of Gaussians.
1030 """
1031 spec = np.zeros(len(f), dtype=complex)
1032 for dt, a, s in zip(pas[0:-2:3], pas[1:-1:3], pas[2::3]):
1033 gauss = a*np.sqrt(2*np.pi)*s*np.exp(-0.5*(2*np.pi*s*f)**2)
1034 shift = np.exp(-2j*np.pi*f*dt)
1035 spec += gauss*shift
1036 spec *= np.sqrt(2) # because of one-sided spectrum
1037 return np.abs(spec)**2
1040def gaussian_sum_costs(pas, time, eod, freqs, energy):
1041 """ Cost function for fitting sum of Gaussian to pulse EOD.
1043 Parameters
1044 ----------
1045 pas: list of floats
1046 The pulse parameters in a flat list.
1047 Position, amplitude, and standard deviation of first phase,
1048 position, amplitude, and standard deviation of second phase,
1049 and so on.
1050 time: 1-D array of float
1051 Time points of the EOD waveform.
1052 eod: 1-D array of float
1053 The real EOD waveform.
1054 freqs: 1-D array of float
1055 The frequency components of the spectrum.
1056 energy: 1-D array of float
1057 The energy spectrum of the real pulse.
1059 Returns
1060 -------
1061 costs: float
1062 Weighted sum of rms waveform difference and rms spectral difference.
1063 """
1064 eod_fit = gaussian_sum(time, *pas)
1065 eod_rms = np.sqrt(np.mean((eod_fit - eod)**2))/np.max(np.abs(eod))
1066 level = decibel(energy)
1067 level_range = 30
1068 n = np.argmax(level)
1069 energy_fit = gaussian_sum_spectrum(freqs, *pas)
1070 level_fit = decibel(energy_fit)
1071 weight = np.ones(n)
1072 weight[freqs[:n] < 10] = 100
1073 weight /= np.sum(weight)
1074 spec_rms = np.sqrt(np.mean(weight*(level_fit[:n] - level[:n])**2))/level_range
1075 costs = eod_rms + 5*spec_rms
1076 #print(f'{costs:.4f}, {eod_rms:.4f}, {spec_rms:.4f}')
1077 return costs
1080def decompose_pulse(eod, freqs, energy, phases, width_frac=0.5):
1081 """Decompose single pulse waveform into sum of Gaussians.
1083 Use the output to simulate pulse-type EODs using the functions
1084 provided in the thunderfish.fakefish module.
1086 Parameters
1087 ----------
1088 eod: 2-D array of float
1089 The eod waveform. First column is time in seconds and second
1090 column is the eod waveform. Further columns are ignored.
1091 freqs: 1-D array of float
1092 The frequency components of the spectrum.
1093 energy: 1-D array of float
1094 The energy spectrum of the real pulse.
1095 phases: 2-D array of float
1096 Properties of the EOD phases as returned by analyze_pulse_phases().
1097 width_frac: float
1098 The width of a peak is measured at this fraction of a peak's
1099 height (0-1).
1101 Returns
1102 -------
1103 pulse: dict
1104 Dictionary with phase times, amplitudes and standard
1105 deviations of Gaussians fitted to the pulse waveform. Use the
1106 functions provided in thunderfish.fakefish to simulate pulse
1107 fish EODs from this data.
1109 """
1110 pulse = dict(times=np.zeros(0), amplitudes=np.zeros(0),
1111 stdevs=np.zeros(0))
1112 if len(phases) < 1:
1113 return pulse
1114 if eod.ndim == 2:
1115 time = eod[:, 0]
1116 eod = eod[:, 1]
1117 elif isinstance(ratetime, (list, tuple, np.ndarray)):
1118 time = ratetime
1119 else:
1120 time = np.arange(len(eod))/rate
1121 # convert half width to standard deviation:
1122 fac = 0.5/np.sqrt(-2*np.log(width_frac))
1123 # fit parameter as single list:
1124 tas = []
1125 for t, a, s in zip(phases[:, 1], phases[:, 2], phases[:, 4]*fac):
1126 tas.extend((t, a, s))
1127 tas = np.asarray(tas)
1128 # fit EOD waveform:
1129 try:
1130 tas, _ = curve_fit(gaussian_sum, time, eod, tas)
1131 except RuntimeError as e:
1132 print('Fit gaussian_sum failed in decompose_pulse():', e)
1133 return pulse
1134 # fit EOD waveform and spectrum:
1135 bnds = [(1e-5, None) if k%3 == 2 else (None, None)
1136 for k in range(len(tas))]
1137 res = minimize(gaussian_sum_costs, tas,
1138 args=(time, eod, freqs, energy), bounds=bnds)
1139 if not res.success:
1140 print('Optimization gaussian_sum_costs failed in decompose_pulse():',
1141 res.message)
1142 else:
1143 tas = res.x
1144 # add another Gaussian:
1145 rms_norm = np.max(np.abs(eod))
1146 rms_old = np.sqrt(np.mean((gaussian_sum(time, *tas) - eod)**2))/rms_norm
1147 eod_diff = np.abs(gaussian_sum(time, *tas) - eod)/rms_norm
1148 if np.max(eod_diff) > 0.1:
1149 print('add Gaussian', np.max(eod_diff))
1150 ntas = np.concatenate((tas, (time[np.argmax(eod_diff)], np.max(eod_diff),
1151 np.mean(tas[2::3]))))
1152 bnds = [(1e-5, None) if k%3 == 2 else (None, None)
1153 for k in range(len(ntas))]
1154 res = minimize(gaussian_sum_costs, ntas,
1155 args=(time, eod, freqs, energy), bounds=bnds)
1156 if res.success:
1157 rms_new = np.sqrt(np.mean((gaussian_sum(time, *res.x) - eod)**2))/rms_norm
1158 if rms_new < 0.8*rms_old:
1159 tas = res.x
1160 else:
1161 print('removed new Gaussian')
1162 else:
1163 print('Optimization gaussian_sum_costs for additional Gaussian failed in decompose_pulse():',
1164 res.message)
1165 times = np.asarray(tas[0::3])
1166 ampls = np.asarray(tas[1::3])
1167 stdevs = np.asarray(tas[2::3])
1168 pulse = dict(times=times, amplitudes=ampls, stdevs=stdevs)
1169 return pulse
1172def analyze_pulse_tail(peak_index, eod, ratetime=None,
1173 threshold=0.0, fit_frac=0.5):
1174 """ Fit exponential to last peak/trough of pulse EOD.
1176 Parameters
1177 ----------
1178 peak_index: int
1179 Index of last peak in `eod`.
1180 eod: 1-D or 2-D array
1181 The eod waveform to be analyzed.
1182 If an 1-D array, then this is the waveform and you
1183 need to also pass a sampling rate in `rate`.
1184 If a 2-D array, then first column is time in seconds and second
1185 column is the eod waveform. Further columns are ignored.
1186 ratetime: None or float or array of float
1187 If a 1-D array is passed on to `eod` then either the sampling
1188 rate in Hertz or the time array corresponding to `eod`.
1189 threshold: float
1190 Minimum size of peaks of the pulse waveform.
1191 fit_frac: float or None
1192 An exponential is fitted to the tail of the last peak/trough
1193 starting where the waveform falls below this fraction of the
1194 peak's height (0-1).
1195 If None, do not attempt to fit.
1197 Returns
1198 -------
1199 tau: float
1200 Time constant of pulse tail in seconds.
1201 fit: 1-D array of float
1202 Time trace of the fit corresponding to `eod`.
1203 """
1204 if fit_frac is None:
1205 return None, None
1206 if eod.ndim == 2:
1207 time = eod[:, 0]
1208 eod = eod[:, 1]
1209 elif isinstance(ratetime, (list, tuple, np.ndarray)):
1210 time = ratetime
1211 else:
1212 time = np.arange(len(eod))/rate
1213 dt = np.mean(np.diff(time))
1214 pi = peak_index
1215 # positive or negative decay:
1216 sign = 1
1217 eodpp = eod[pi:] - 0.5*threshold
1218 eodpn = -eod[pi:] - 0.5*threshold
1219 if np.sum(eodpn[eodpn > 0]) > np.sum(eodpp[eodpp > 0]):
1220 sign = -1
1221 if sign*eod[pi] < 0:
1222 pi += np.argmax(sign*eod[pi:])
1223 pi_ampl = np.abs(eod[pi])
1224 n = len(eod[pi:])
1225 # no sufficiently large initial value:
1226 if pi_ampl*fit_frac <= 0.5*threshold:
1227 fit_frac = False
1228 # no sufficiently long decay:
1229 if n < 10:
1230 fit_frac = False
1231 # not decaying towards zero:
1232 max_line = pi_ampl - (pi_ampl - threshold)*np.arange(n)/n + 1e-8
1233 if np.sum(np.abs(eod[pi + 2:]) > max_line[2:])/n > 0.05:
1234 fit_frac = False
1235 if not fit_frac:
1236 return None, None
1237 thresh = eod[pi]*fit_frac
1238 inx = pi + np.argmax(sign*eod[pi:] < sign*thresh)
1239 thresh = eod[inx]*np.exp(-1.0)
1240 tau_inx = np.argmax(sign*eod[inx:] < sign*thresh)
1241 if tau_inx < 2:
1242 tau_inx = 2
1243 rridx = inx + 6*tau_inx
1244 if rridx > len(eod) - 1:
1245 rridx = len(eod) - 1
1246 if rridx - inx < 3*tau_inx:
1247 return None, None
1248 tau = time[inx + tau_inx] - time[inx]
1249 params = [tau, eod[inx] - eod[rridx], eod[rridx]]
1250 try:
1251 popt, pcov = curve_fit(exp_decay, time[inx:rridx] - time[inx],
1252 eod[inx:rridx], params,
1253 bounds=([0.0, -np.inf, -np.inf], np.inf))
1254 except TypeError:
1255 popt, pcov = curve_fit(exp_decay, time[inx:rridx] - time[inx],
1256 eod[inx:rridx], params)
1257 if popt[0] > 1.2*tau:
1258 tau_inx = int(np.round(popt[0]/dt))
1259 rridx = inx + 6*tau_inx
1260 if rridx > len(eod) - 1:
1261 rridx = len(eod) - 1
1262 try:
1263 popt, pcov = curve_fit(exp_decay, time[inx:rridx] - time[inx],
1264 eod[inx:rridx], popt,
1265 bounds=([0.0, -np.inf, -np.inf], np.inf))
1266 except TypeError:
1267 popt, pcov = curve_fit(exp_decay, time[inx:rridx] - time[inx],
1268 eod[inx:rridx], popt)
1269 tau = popt[0]
1270 fit = np.zeros(len(eod))
1271 fit[:] = np.nan
1272 fit[inx:rridx] = exp_decay(time[inx:rridx] - time[inx], *popt)
1273 return tau, fit
1276def analyze_pulse_intervals(eod_times, ipi_cv_thresh=0.5,
1277 ipi_percentile=30.0):
1278 """ Basic statistics of interpulse intervals.
1280 Parameters
1281 ----------
1282 eod_times: 1-D array or None
1283 List of times of detected EODs.
1284 ipi_cv_thresh: float
1285 If the coefficient of variation of the interpulse intervals
1286 (IPIs) is smaller than this threshold, then the statistics of
1287 IPIs is estimated from all IPIs. Otherwise only intervals
1288 smaller than a certain percentile are used.
1289 ipi_percentile: float
1290 When computing the statistics of IPIs from a subset of the
1291 IPIs, only intervals smaller than this percentile (between 0
1292 and 100) are used.
1294 Returns
1295 -------
1296 ipi_median: float
1297 Median inter-pulse interval.
1298 ipi_mean: float
1299 Mean inter-pulse interval.
1300 ipi_std: float
1301 Standard deviation of inter-pulse intervals.
1303 """
1304 if eod_times is None:
1305 return None, None, None
1306 inter_pulse_intervals = np.diff(eod_times)
1307 ipi_cv = np.std(inter_pulse_intervals)/np.mean(inter_pulse_intervals)
1308 if ipi_cv < ipi_cv_thresh:
1309 ipi_median = np.median(inter_pulse_intervals)
1310 ipi_mean = np.mean(inter_pulse_intervals)
1311 ipi_std = np.std(inter_pulse_intervals)
1312 else:
1313 intervals = inter_pulse_intervals[inter_pulse_intervals <
1314 np.percentile(inter_pulse_intervals, ipi_percentile)]
1315 ipi_median = np.median(intervals)
1316 ipi_mean = np.mean(intervals)
1317 ipi_std = np.std(intervals)
1318 return ipi_median, ipi_mean, ipi_std
1321def analyze_pulse(eod, eod_times=None, min_pulse_win=0.001,
1322 peak_thresh_fac=0.01, min_dist=50.0e-6,
1323 width_frac=0.5, fit_frac=0.5,
1324 freq_resolution=1.0, fade_frac=0.0,
1325 flip_pulse='none', ipi_cv_thresh=0.5,
1326 ipi_percentile=30.0):
1327 """Analyze the EOD waveform of a pulse fish.
1329 Parameters
1330 ----------
1331 eod: 2-D array
1332 The eod waveform. First column is time in seconds, second
1333 column the EOD waveform, third column, if present, is the
1334 standard error of the EOD waveform, Further columns are
1335 optional but not used.
1336 eod_times: 1-D array or None
1337 List of times of detected EOD peaks.
1338 min_pulse_win: float
1339 The minimum size of cut-out EOD waveform.
1340 peak_thresh_fac: float
1341 Set the threshold for peak detection to the maximum pulse
1342 amplitude times this factor.
1343 min_dist: float
1344 Minimum distance between peak and troughs of the pulse.
1345 width_frac: float
1346 The width of a peak is measured at this fraction of a peak's
1347 height (0-1).
1348 fit_frac: float or None
1349 An exponential is fitted to the tail of the last peak/trough
1350 starting where the waveform falls below this fraction of the
1351 peak's height (0-1).
1352 freq_resolution: float
1353 The frequency resolution of the spectrum of the single pulse.
1354 fade_frac: float
1355 Fraction of time of the EOD waveform that is fade in and out
1356 to zero baseline.
1357 flip_pulse: 'auto', 'none', 'flip'
1358 - 'auto' flip waveform such that the first large extremum is positive.
1359 - 'flip' flip waveform.
1360 - 'none' do not flip waveform.
1361 ipi_cv_thresh: float
1362 If the coefficient of variation of the interpulse intervals
1363 (IPIs) is smaller than this threshold, then the statistics of
1364 IPIs is estimated from all IPIs. Otherwise only intervals
1365 smaller than a certain percentile are used.
1366 ipi_percentile: float
1367 When computing the statistics of IPIs from a subset of the
1368 IPIs, only intervals smaller than this percentile (between 0
1369 and 100) are used.
1371 Returns
1372 -------
1373 meod: 2-D array of floats
1374 The eod waveform. First column is time in seconds,
1375 second column the eod waveform.
1376 Further columns are kept from the input `eod`.
1377 As the two last columns the waveform resulting from the
1378 decomposition into Gaussians and the fit to the tail of the
1379 last peak is appended.
1380 props: dict
1381 A dictionary with properties of the analyzed EOD waveform.
1383 - type: set to 'pulse'.
1384 - EODf: the inverse of the median interval between `eod_times`,
1385 if provided.
1386 - period: the median interval between `eod_times`, if provided.
1387 - IPI-mean: the mean interval between `eod_times`, if provided.
1388 - IPI-std: the standard deviation of the intervals between
1389 `eod_times`, if provided.
1390 - max-ampl: the amplitude of the largest positive peak (P1).
1391 - min-ampl: the amplitude of the largest negative peak (P2).
1392 - p-p-amplitude: peak-to-peak amplitude of the EOD waveform.
1393 - noise: root-mean squared standard error mean of the averaged
1394 EOD waveform relative to the p-p amplitude.
1395 - tstart: time in seconds where the pulse starts,
1396 i.e. crosses the threshold for the first time.
1397 - tend: time in seconds where the pulse ends,
1398 i.e. crosses the threshold for the last time.
1399 - width: total width of the pulse in seconds (tend-tstart).
1400 - peak-dist: distance between minimum and maximum phase in seconds.
1401 - tau: time constant of exponential decay of pulse tail in seconds.
1402 - firstphase: index of the first phase in the pulse (i.e. -1 for P-1)
1403 - lastphase: index of the last phase in the pulse (i.e. 3 for P3)
1404 - totalarea: sum of areas under positive and negative peaks.
1405 - positivearea: area under positive phases relative to total area.
1406 - negativearea: area under negative phases relative to total area.
1407 - polaritybalance: contrast between areas under positive and
1408 negative phases.
1409 - peakfreq: frequency at peak energy of the single pulse spectrum
1410 in Hertz.
1411 - peakenergy: peak energy of the single pulse spectrum.
1412 - troughfreq: frequency at trough before peak in Hertz.
1413 - troughenergy: energy of trough before peak in x^2 s/Hz.
1414 - energyatt5: attenuation of average energy below 5 Hz relative to
1415 peak energy in decibel.
1416 - energyatt50: attenuation of average energy below 50 Hz relative to
1417 peak energy in decibel.
1418 - lowcutoff: frequency at which the energy reached half of the
1419 peak energy relative to the DC energy in Hertz.
1420 - highcutoff: 3dB roll-off frequency of spectrum in Hertz.
1421 - flipped: True if the waveform was flipped.
1422 - n: number of pulses analyzed (i.e. `len(eod_times)`), if provided.
1423 - times: the times of the detected EOD pulses (i.e. `eod_times`),
1424 if provided.
1426 Empty if waveform is not a pulse EOD.
1427 phases: 2-D array
1428 For each phase (rows) of the EOD waveform
1429 7 columns: the phase index (1 is P1, i.e. the largest positive phase),
1430 time relative to largest positive phase, amplitude,
1431 amplitude normalized to P1 phase,
1432 width of phase at half height,
1433 area under the phase, and area under the phase relative to total area.
1434 Empty if waveform is not a pulse EOD.
1435 pulse: dict
1436 Dictionary with phase times, amplitudes and standard
1437 deviations of Gaussians fitted to the pulse waveform. Use the
1438 functions provided in thunderfish.fakefish to simulate pulse
1439 fish EODs from this data.
1440 energy: 2-D array
1441 The energy spectrum of a single pulse. First column are the
1442 frequencies, second column the energy in x^2 s/Hz.
1443 Empty if waveform is not a pulse EOD.
1445 """
1446 # storage:
1447 meod = np.zeros((eod.shape[0], eod.shape[1] + 2))
1448 meod[:, :eod.shape[1]] = eod
1449 meod[:, -2] = np.nan
1450 meod[:, -1] = np.nan
1451 toffs = 0
1453 # cut out stable estimate if standard deviation is available:
1454 if eod.shape[1] > 2 and np.max(meod[:, 2]) > 3*np.min(meod[:, 2]):
1455 n = len(meod)
1456 idx0 = np.argmax(np.abs(meod[n//10:9*n//10,1])) + n//10
1457 toffs += meod[idx0,0]
1458 meod[:, 0] -= meod[idx0,0]
1459 # minimum in standard deviation:
1460 lstd_idx = np.argmin(meod[:idx0-5,2])
1461 rstd_idx = np.argmin(meod[idx0+5:, 2]) + idx0
1462 # central, left, and right maximum of standard deviation:
1463 max_std = np.max(meod[lstd_idx:rstd_idx,2])
1464 l_std = np.max(meod[:len(meod)//4,2])
1465 r_std = np.max(meod[len(meod)*3//4:, 2])
1466 lidx = 0
1467 ridx = len(meod)
1468 if l_std > max_std and lstd_idx > lidx:
1469 lidx = lstd_idx - np.argmax(meod[lstd_idx:0:-1,2] >= 0.25*l_std + 0.75*meod[lstd_idx,2])
1470 if r_std > max_std and rstd_idx < ridx:
1471 ridx = rstd_idx + np.argmax(meod[rstd_idx:, 2] >= 0.25*r_std + 0.75*meod[rstd_idx,2])
1472 #plt.plot(meod[:, 0], meod[:, 1])
1473 #plt.plot(meod[:, 0], meod[:, 2], '-r')
1474 #plt.plot([meod[lidx,0], meod[lidx,0]], [-0.1, 0.1], '-k')
1475 #plt.plot([meod[ridx-1,0], meod[ridx-1,0]], [-0.1, 0.1], '-b')
1476 #plt.show()
1477 meod = meod[lidx:ridx,:]
1479 # subtract mean computed from the ends of the snippet:
1480 n = len(meod)//20 if len(meod) >= 20 else 1
1481 meod[:, 1] -= 0.5*(np.mean(meod[:n,1]) + np.mean(meod[-n:, 1]))
1483 # largest positive and negative peak and flip: # TODO move after cutting!
1484 flipped = False
1485 max_idx = np.argmax(meod[:, 1])
1486 max_ampl = np.abs(meod[max_idx,1])
1487 min_idx = np.argmin(meod[:, 1])
1488 min_ampl = np.abs(meod[min_idx,1])
1489 amplitude = np.max((max_ampl, min_ampl))
1490 if max_ampl > 0.2*amplitude and min_ampl > 0.2*amplitude:
1491 # two major peaks:
1492 if 'flip' in flip_pulse or ('auto' in flip_pulse and min_idx < max_idx):
1493 # flip:
1494 meod[:, 1] = -meod[:, 1]
1495 peak_idx = min_idx
1496 min_idx = max_idx
1497 max_idx = peak_idx
1498 flipped = True
1499 elif 'flip' in flip_pulse or ('auto' in flip_pulse and min_ampl > 0.2*amplitude):
1500 # flip:
1501 meod[:, 1] = -meod[:, 1]
1502 peak_idx = min_idx
1503 min_idx = max_idx
1504 max_idx = peak_idx
1505 flipped = True
1506 max_ampl = np.abs(meod[max_idx, 1])
1507 min_ampl = np.abs(meod[min_idx, 1])
1509 # minimum threshold for peak detection: # return from cutting function
1510 n = len(meod[:, 1])//10 if len(meod) >= 20 else 2
1511 thl_max = np.max(meod[:n, 1])
1512 thl_min = np.min(meod[:n, 1])
1513 thr_max = np.max(meod[-n:, 1])
1514 thr_min = np.min(meod[-n:, 1])
1515 min_thresh = 2*np.max([thl_max, thr_max]) - np.min([thl_min, thr_min])
1516 if min_thresh > 0.5*(max_ampl + min_ampl):
1517 min_thresh = 0.5*(max_ampl + min_ampl)
1518 fit_frac = None
1519 # threshold for peak detection:
1520 threshold = max_ampl*peak_thresh_fac
1521 if threshold < min_thresh:
1522 threshold = min_thresh
1523 if threshold <= 0.0:
1524 return meod, {}, [], []
1526 # cut out relevant signal:
1527 # TODO: this needs to be improved!
1528 # more to the right in case of a shallow decay to baseline!
1529 lidx = np.argmax(np.abs(meod[:, 1]) > threshold)
1530 ridx = len(meod) - 1 - np.argmax(np.abs(meod[::-1,1]) > threshold)
1531 t0 = meod[lidx,0]
1532 t1 = meod[ridx,0]
1533 width = t1 - t0
1534 if width < min_pulse_win:
1535 width = min_pulse_win
1536 dt = meod[1, 0] - meod[0, 0]
1537 width_idx = int(np.round(width/dt))
1538 # expand width:
1539 leidx = lidx - width_idx//2
1540 if leidx < 0:
1541 leidx = 0
1542 reidx = ridx + width_idx//2
1543 if reidx >= len(meod):
1544 reidx = len(meod)
1545 meod = meod[leidx:reidx,:]
1546 max_idx -= leidx
1547 min_idx -= leidx
1549 # move peak of waveform to zero:
1550 toffs += meod[max_idx, 0]
1551 meod[:, 0] -= meod[max_idx, 0]
1553 tau = None
1554 dist = 0.0
1555 peaks = []
1557 # amplitude and variance:
1558 ppampl = max_ampl + min_ampl
1559 dist = meod[min_idx, 0] - meod[max_idx, 0] # TODO: move to analyze_pulse_phases() and handle non existent minimum
1560 rmssem = np.sqrt(np.mean(meod[:, 2]**2.0))/ppampl if eod.shape[1] > 2 else None
1562 # integrals and polarity balance:
1563 pos_area = np.sum(meod[meod[:, 1] >= 0, 1])*dt
1564 neg_area = np.sum(meod[meod[:, 1] < 0, 1])*dt
1565 total_area = np.sum(np.abs(meod[:, 1]))*dt
1566 integral_area = np.sum(meod[:, 1])*dt # = pos_area - neg_area
1567 polarity_balance = integral_area/total_area
1569 # characterize EOD phases:
1570 phases = analyze_pulse_phases(threshold, meod, None,
1571 min_dist=min_dist, width_frac=width_frac)
1573 # fit exponential to last peak/trough:
1574 if len(phases) > 1:
1575 pi = np.argmin(np.abs(meod[:, 0] - phases[-1, 1]))
1576 tau, fit = analyze_pulse_tail(pi, meod, None,
1577 threshold, fit_frac)
1578 if fit is not None:
1579 meod[:, -1] = fit
1581 # energy spectrum of single pulse:
1582 freqs, energy = pulse_spectrum(meod, None, freq_resolution, fade_frac)
1583 # store spectrum:
1584 eenergy = np.zeros((len(energy), 2))
1585 eenergy[:, 0] = freqs
1586 eenergy[:, 1] = energy
1587 # analyse spectrum:
1588 peakfreq, peakenergy, troughfreq, troughenergy, \
1589 att5, att50, lowcutoff, highcutoff = \
1590 analyze_pulse_spectrum(freqs, energy)
1592 # decompose EOD waveform:
1593 pulse = decompose_pulse(meod, freqs, energy, phases, width_frac)
1594 if len(pulse['amplitudes']) > 0:
1595 eod_fit = pulsefish_waveform(pulse, meod[:, 0])
1596 meod[:, -2] = eod_fit
1598 # analyze pulse intervals:
1599 ipi_median, ipi_mean, ipi_std = analyze_pulse_intervals(eod_times,
1600 ipi_cv_thresh,
1601 ipi_percentile)
1603 # store properties:
1604 props = {}
1605 props['type'] = 'pulse'
1606 if eod_times is not None:
1607 props['EODf'] = 1.0/ipi_median
1608 props['period'] = ipi_median
1609 props['IPI-mean'] = ipi_mean
1610 props['IPI-std'] = ipi_std
1611 props['max-ampl'] = max_ampl
1612 props['min-ampl'] = min_ampl
1613 props['p-p-amplitude'] = ppampl
1614 if rmssem:
1615 props['noise'] = rmssem
1616 props['tstart'] = t0
1617 props['tend'] = t1
1618 props['width'] = t1 - t0
1619 props['peak-dist'] = dist
1620 if tau:
1621 props['tau'] = tau
1622 props['firstphase'] = phases[0, 0] if len(phases) > 0 else 1
1623 props['lastphase'] = phases[-1, 0] if len(phases) > 0 else 1
1624 props['totalarea'] = total_area
1625 props['positivearea'] = pos_area/total_area
1626 props['negativearea'] = neg_area/total_area
1627 props['polaritybalance'] = polarity_balance
1628 props['peakfreq'] = peakfreq
1629 props['peakenergy'] = peakenergy
1630 props['troughfreq'] = troughfreq
1631 props['troughenergy'] = troughenergy
1632 props['energyatt5'] = att5
1633 props['energyatt50'] = att50
1634 props['lowcutoff'] = lowcutoff
1635 props['highcutoff'] = highcutoff
1636 props['flipped'] = flipped
1637 if eod_times is not None:
1638 props['n'] = len(eod_times)
1639 props['times'] = eod_times + toffs
1641 return meod, props, phases, pulse, eenergy
1644def load_species_waveforms(species_file='none'):
1645 """Load template EOD waveforms for species matching.
1647 Parameters
1648 ----------
1649 species_file: string
1650 Name of file containing species definitions. The content of
1651 this file is as follows:
1653 - Empty lines and line starting with a hash ('#') are skipped.
1654 - A line with the key-word 'wavefish' marks the beginning of the
1655 table for wave fish.
1656 - A line with the key-word 'pulsefish' marks the beginning of the
1657 table for pulse fish.
1658 - Each line in a species table has three fields,
1659 separated by colons (':'):
1661 1. First field is an abbreviation of the species name.
1662 2. Second field is the filename of the recording containing the
1663 EOD waveform.
1664 3. The optional third field is the EOD frequency of the EOD waveform.
1666 The EOD frequency is used to normalize the time axis of a
1667 wave fish EOD to one EOD period. If it is not specified in
1668 the third field, it is taken from the corresponding
1669 *-wavespectrum-* file, if present. Otherwise the species is
1670 excluded and a warning is issued.
1672 Example file content:
1673 ``` plain
1674 Wavefish
1675 Aptero : F_91009L20-eodwaveform-0.csv : 823Hz
1676 Eigen : C_91008L01-eodwaveform-0.csv
1678 Pulsefish
1679 Gymnotus : pulsefish/gymnotus.csv
1680 Brachy : H_91009L11-eodwaveform-0.csv
1681 ```
1683 Returns
1684 -------
1685 wave_names: list of strings
1686 List of species names of wave-type fish.
1687 wave_eods: list of 2-D arrays
1688 List of EOD waveforms of wave-type fish corresponding to
1689 `wave_names`. First column of a waveform is time in seconds,
1690 second column is the EOD waveform. The waveforms contain
1691 exactly one EOD period.
1692 pulse_names: list of strings
1693 List of species names of pulse-type fish.
1694 pulse_eods: list of 2-D arrays
1695 List of EOD waveforms of pulse-type fish corresponding to
1696 `pulse_names`. First column of a waveform is time in seconds,
1697 second column is the EOD waveform.
1698 """
1699 if not Path(species_file).is_file():
1700 return [], [], [], []
1701 wave_names = []
1702 wave_eods = []
1703 pulse_names = []
1704 pulse_eods = []
1705 fish_type = 'wave'
1706 with open(species_file, 'r') as sf:
1707 for line in sf:
1708 line = line.strip()
1709 if len(line) == 0 or line[0] == '#':
1710 continue
1711 if line.lower() == 'wavefish':
1712 fish_type = 'wave'
1713 elif line.lower() == 'pulsefish':
1714 fish_type = 'pulse'
1715 else:
1716 ls = [s.strip() for s in line.split(':')]
1717 if len(ls) < 2:
1718 continue
1719 name = ls[0]
1720 waveform_file = ls[1]
1721 eod = TableData(waveform_file).array()
1722 eod[:, 0] *= 0.001
1723 if fish_type == 'wave':
1724 eodf = None
1725 if len(ls) > 2:
1726 eodf = float(ls[2].replace('Hz', '').strip())
1727 else:
1728 spectrum_file = waveform_file.replace('eodwaveform', 'wavespectrum')
1729 try:
1730 wave_spec = TableData(spectrum_file)
1731 eodf = wave_spec[0, 1]
1732 except FileNotFoundError:
1733 pass
1734 if eodf is None:
1735 print('warning: unknown EOD frequency of %s. Skip.' % name)
1736 continue
1737 eod[:, 0] *= eodf
1738 wave_names.append(name)
1739 wave_eods.append(eod[:, :2])
1740 elif fish_type == 'pulse':
1741 pulse_names.append(name)
1742 pulse_eods.append(eod[:, :2])
1743 return wave_names, wave_eods, pulse_names, pulse_eods
1746def wave_similarity(eod1, eod2, eod1f=1.0, eod2f=1.0):
1747 """Root-mean squared difference between two wave fish EODs.
1749 Compute the root-mean squared difference between two wave fish
1750 EODs over one period. The better sampled signal is down-sampled to
1751 the worse sampled one. Amplitude are normalized to peak-to-peak
1752 amplitude before computing rms difference. Also compute the rms
1753 difference between the one EOD and the other one inverted in
1754 amplitude. The smaller of the two rms values is returned.
1756 Parameters
1757 ----------
1758 eod1: 2-D array
1759 Time and amplitude of reference EOD.
1760 eod2: 2-D array
1761 Time and amplitude of EOD that is to be compared to `eod1`.
1762 eod1f: float
1763 EOD frequency of `eod1` used to transform the time axis of `eod1`
1764 to multiples of the EOD period. If already normalized to EOD period,
1765 as for example by the `load_species_waveforms()` function, then
1766 set the EOD frequency to one (default).
1767 eod2f: float
1768 EOD frequency of `eod2` used to transform the time axis of `eod2`
1769 to multiples of the EOD period. If already normalized to EOD period,
1770 as for example by the `load_species_waveforms()` function, then
1771 set the EOD frequency to one (default).
1773 Returns
1774 -------
1775 rmse: float
1776 Root-mean-squared difference between the two EOD waveforms relative to
1777 their standard deviation over one period.
1778 """
1779 # copy:
1780 eod1 = np.array(eod1[:, :2])
1781 eod2 = np.array(eod2[:, :2])
1782 # scale to multiples of EOD period:
1783 eod1[:, 0] *= eod1f
1784 eod2[:, 0] *= eod2f
1785 # make eod1 the waveform with less samples per period:
1786 n1 = int(1.0/(eod1[1,0]-eod1[0,0]))
1787 n2 = int(1.0/(eod2[1,0]-eod2[0,0]))
1788 if n1 > n2:
1789 eod1, eod2 = eod2, eod1
1790 n1, n2 = n2, n1
1791 # one period around time zero:
1792 i0 = np.argmin(np.abs(eod1[:, 0]))
1793 k0 = i0-n1//2
1794 if k0 < 0:
1795 k0 = 0
1796 k1 = k0 + n1 + 1
1797 if k1 >= len(eod1):
1798 k1 = len(eod1)
1799 # cut out one period from the reference EOD around maximum:
1800 i = k0 + np.argmax(eod1[k0:k1,1])
1801 k0 = i-n1//2
1802 if k0 < 0:
1803 k0 = 0
1804 k1 = k0 + n1 + 1
1805 if k1 >= len(eod1):
1806 k1 = len(eod1)
1807 eod1 = eod1[k0:k1,:]
1808 # normalize amplitudes of first EOD:
1809 eod1[:, 1] -= np.min(eod1[:, 1])
1810 eod1[:, 1] /= np.max(eod1[:, 1])
1811 sigma = np.std(eod1[:, 1])
1812 # set time zero to maximum of second EOD:
1813 i0 = np.argmin(np.abs(eod2[:, 0]))
1814 k0 = i0-n2//2
1815 if k0 < 0:
1816 k0 = 0
1817 k1 = k0 + n2 + 1
1818 if k1 >= len(eod2):
1819 k1 = len(eod2)
1820 i = k0 + np.argmax(eod2[k0:k1,1])
1821 eod2[:, 0] -= eod2[i,0]
1822 # interpolate eod2 to the time base of eod1:
1823 eod2w = np.interp(eod1[:, 0], eod2[:, 0], eod2[:, 1])
1824 # normalize amplitudes of second EOD:
1825 eod2w -= np.min(eod2w)
1826 eod2w /= np.max(eod2w)
1827 # root-mean-square difference:
1828 rmse1 = np.sqrt(np.mean((eod1[:, 1] - eod2w)**2))/sigma
1829 # root-mean-square difference of the flipped signal:
1830 i = k0 + np.argmin(eod2[k0:k1,1])
1831 eod2[:, 0] -= eod2[i,0]
1832 eod2w = np.interp(eod1[:, 0], eod2[:, 0], -eod2[:, 1])
1833 eod2w -= np.min(eod2w)
1834 eod2w /= np.max(eod2w)
1835 rmse2 = np.sqrt(np.mean((eod1[:, 1] - eod2w)**2))/sigma
1836 # take the smaller value:
1837 rmse = min(rmse1, rmse2)
1838 return rmse
1841def pulse_similarity(eod1, eod2, wfac=10.0):
1842 """Root-mean squared difference between two pulse fish EODs.
1844 Compute the root-mean squared difference between two pulse fish
1845 EODs over `wfac` times the distance between minimum and maximum of
1846 the wider EOD. The waveforms are normalized to their maxima prior
1847 to computing the rms difference. Also compute the rms difference
1848 between the one EOD and the other one inverted in amplitude. The
1849 smaller of the two rms values is returned.
1851 Parameters
1852 ----------
1853 eod1: 2-D array
1854 Time and amplitude of reference EOD.
1855 eod2: 2-D array
1856 Time and amplitude of EOD that is to be compared to `eod1`.
1857 wfac: float
1858 Multiply the distance between minimum and maximum by this factor
1859 to get the window size over which to compute the rms difference.
1861 Returns
1862 -------
1863 rmse: float
1864 Root-mean-squared difference between the two EOD waveforms
1865 relative to their standard deviation over the analysis window.
1866 """
1867 # copy:
1868 eod1 = np.array(eod1[:, :2])
1869 eod2 = np.array(eod2[:, :2])
1870 # width of the pulses:
1871 imin1 = np.argmin(eod1[:, 1])
1872 imax1 = np.argmax(eod1[:, 1])
1873 w1 = np.abs(eod1[imax1,0]-eod1[imin1,0])
1874 imin2 = np.argmin(eod2[:, 1])
1875 imax2 = np.argmax(eod2[:, 1])
1876 w2 = np.abs(eod2[imax2,0]-eod2[imin2,0])
1877 w = wfac*max(w1, w2)
1878 # cut out signal from the reference EOD:
1879 n = int(w/(eod1[1,0]-eod1[0,0]))
1880 i0 = imax1-n//2
1881 if i0 < 0:
1882 i0 = 0
1883 i1 = imax1+n//2+1
1884 if i1 >= len(eod1):
1885 i1 = len(eod1)
1886 eod1[:, 0] -= eod1[imax1,0]
1887 eod1 = eod1[i0:i1,:]
1888 # normalize amplitude of first EOD:
1889 eod1[:, 1] /= np.max(eod1[:, 1])
1890 sigma = np.std(eod1[:, 1])
1891 # interpolate eod2 to the time base of eod1:
1892 eod2[:, 0] -= eod2[imax2,0]
1893 eod2w = np.interp(eod1[:, 0], eod2[:, 0], eod2[:, 1])
1894 # normalize amplitude of second EOD:
1895 eod2w /= np.max(eod2w)
1896 # root-mean-square difference:
1897 rmse1 = np.sqrt(np.mean((eod1[:, 1] - eod2w)**2))/sigma
1898 # root-mean-square difference of the flipped signal:
1899 eod2[:, 0] -= eod2[imin2,0]
1900 eod2w = np.interp(eod1[:, 0], eod2[:, 0], -eod2[:, 1])
1901 eod2w /= np.max(eod2w)
1902 rmse2 = np.sqrt(np.mean((eod1[:, 1] - eod2w)**2))/sigma
1903 # take the smaller value:
1904 rmse = min(rmse1, rmse2)
1905 return rmse
1908def clipped_fraction(data, rate, eod_times, mean_eod,
1909 min_clip=-np.inf, max_clip=np.inf):
1910 """Compute fraction of clipped EOD waveform snippets.
1912 Cut out snippets at each `eod_times` based on time axis of
1913 `mean_eod`. Check which fraction of snippets exceeds clipping
1914 amplitude `min_clip` and `max_clip`.
1916 Parameters
1917 ----------
1918 data: 1-D array of float
1919 The data to be analysed.
1920 rate: float
1921 Sampling rate of the data in Hertz.
1922 eod_times: 1-D array of float
1923 Array of EOD times in seconds.
1924 mean_eod: 2-D array with time, mean, sem, and fit.
1925 Averaged EOD waveform of wave fish. Only the time axis is used
1926 to set width of snippets.
1927 min_clip: float
1928 Minimum amplitude that is not clipped.
1929 max_clip: float
1930 Maximum amplitude that is not clipped.
1932 Returns
1933 -------
1934 clipped_frac: float
1935 Fraction of snippets that are clipped.
1936 """
1937 # snippets:
1938 idx0 = np.argmin(np.abs(mean_eod[:, 0])) # index of time zero
1939 w0 = -idx0
1940 w1 = len(mean_eod[:, 0]) - idx0
1941 eod_idx = np.round(eod_times*rate).astype(int)
1942 eod_snippets = snippets(data, eod_idx, w0, w1)
1943 # fraction of clipped snippets:
1944 clipped_frac = np.sum(np.any((eod_snippets > max_clip) |
1945 (eod_snippets < min_clip), axis=1))\
1946 / len(eod_snippets)
1947 return clipped_frac
1950def wave_quality(props, harm_relampl=None, min_freq=0.0,
1951 max_freq=2000.0, max_clipped_frac=0.1,
1952 max_crossings=4, max_rms_sem=0.0, max_rms_error=0.05,
1953 min_power=-100.0, max_thd=0.0, max_db_diff=20.0,
1954 max_harmonics_db=-5.0, max_relampl_harm1=0.0,
1955 max_relampl_harm2=0.0, max_relampl_harm3=0.0):
1956 """Assess the quality of an EOD waveform of a wave fish.
1958 Parameters
1959 ----------
1960 props: dict
1961 A dictionary with properties of the analyzed EOD waveform
1962 as returned by `analyze_wave()`.
1963 harm_relampl: 1-D array of floats or None
1964 Relative amplitude of at least the first 3 harmonics without
1965 the fundamental.
1966 min_freq: float
1967 Minimum EOD frequency (`props['EODf']`).
1968 max_freq: float
1969 Maximum EOD frequency (`props['EODf']`).
1970 max_clipped_frac: float
1971 If larger than zero, maximum allowed fraction of clipped data
1972 (`props['clipped']`).
1973 max_crossings: int
1974 If larger than zero, maximum number of zero crossings per EOD period
1975 (`props['ncrossings']`).
1976 max_rms_sem: float
1977 If larger than zero, maximum allowed standard error of the
1978 data relative to p-p amplitude (`props['noise']`).
1979 max_rms_error: float
1980 If larger than zero, maximum allowed root-mean-square error
1981 between EOD waveform and Fourier fit relative to p-p amplitude
1982 (`props['rmserror']`).
1983 min_power: float
1984 Minimum power of the EOD in dB (`props['power']`).
1985 max_thd: float
1986 If larger than zero, then maximum total harmonic distortion
1987 (`props['thd']`).
1988 max_db_diff: float
1989 If larger than zero, maximum standard deviation of differences between
1990 logarithmic powers of harmonics in decibel (`props['dbdiff']`).
1991 Low values enforce smoother power spectra.
1992 max_harmonics_db:
1993 Maximum power of higher harmonics relative to peak power in
1994 decibel (`props['maxdb']`).
1995 max_relampl_harm1: float
1996 If larger than zero, maximum allowed amplitude of first harmonic
1997 relative to fundamental.
1998 max_relampl_harm2: float
1999 If larger than zero, maximum allowed amplitude of second harmonic
2000 relative to fundamental.
2001 max_relampl_harm3: float
2002 If larger than zero, maximum allowed amplitude of third harmonic
2003 relative to fundamental.
2005 Returns
2006 -------
2007 remove: bool
2008 If True then this is most likely not an electric fish. Remove
2009 it from both the waveform properties and the list of EOD
2010 frequencies. If False, keep it in the list of EOD
2011 frequencies, but remove it from the waveform properties if
2012 `skip_reason` is not empty.
2013 skip_reason: string
2014 An empty string if the waveform is good, otherwise a string
2015 indicating the failure.
2016 msg: string
2017 A textual representation of the values tested.
2018 """
2019 remove = False
2020 msg = []
2021 skip_reason = []
2022 # EOD frequency:
2023 if 'EODf' in props:
2024 eodf = props['EODf']
2025 msg += ['EODf=%5.1fHz' % eodf]
2026 if eodf < min_freq or eodf > max_freq:
2027 remove = True
2028 skip_reason += ['invalid EODf=%5.1fHz (minimumFrequency=%5.1fHz, maximumFrequency=%5.1f))' %
2029 (eodf, min_freq, max_freq)]
2030 # clipped fraction:
2031 if 'clipped' in props:
2032 clipped_frac = props['clipped']
2033 msg += ['clipped=%3.0f%%' % (100.0*clipped_frac)]
2034 if max_clipped_frac > 0 and clipped_frac >= max_clipped_frac:
2035 skip_reason += ['clipped=%3.0f%% (maximumClippedFraction=%3.0f%%)' %
2036 (100.0*clipped_frac, 100.0*max_clipped_frac)]
2037 # too many zero crossings:
2038 if 'ncrossings' in props:
2039 ncrossings = props['ncrossings']
2040 msg += ['zero crossings=%d' % ncrossings]
2041 if max_crossings > 0 and ncrossings > max_crossings:
2042 skip_reason += ['too many zero crossings=%d (maximumCrossings=%d)' %
2043 (ncrossings, max_crossings)]
2044 # noise:
2045 rms_sem = None
2046 if 'rmssem' in props:
2047 rms_sem = props['rmssem']
2048 if 'noise' in props:
2049 rms_sem = props['noise']
2050 if rms_sem is not None:
2051 msg += ['rms sem waveform=%6.2f%%' % (100.0*rms_sem)]
2052 if max_rms_sem > 0.0 and rms_sem >= max_rms_sem:
2053 skip_reason += ['noisy waveform s.e.m.=%6.2f%% (max %6.2f%%)' %
2054 (100.0*rms_sem, 100.0*max_rms_sem)]
2055 # fit error:
2056 if 'rmserror' in props:
2057 rms_error = props['rmserror']
2058 msg += ['rmserror=%6.2f%%' % (100.0*rms_error)]
2059 if max_rms_error > 0.0 and rms_error >= max_rms_error:
2060 skip_reason += ['noisy rmserror=%6.2f%% (maximumVariance=%6.2f%%)' %
2061 (100.0*rms_error, 100.0*max_rms_error)]
2062 # wave power:
2063 if 'power' in props:
2064 power = props['power']
2065 msg += ['power=%6.1fdB' % power]
2066 if power < min_power:
2067 skip_reason += ['small power=%6.1fdB (minimumPower=%6.1fdB)' %
2068 (power, min_power)]
2069 # total harmonic distortion:
2070 if 'thd' in props:
2071 thd = props['thd']
2072 msg += ['thd=%5.1f%%' % (100.0*thd)]
2073 if max_thd > 0.0 and thd > max_thd:
2074 skip_reason += ['large THD=%5.1f%% (maxximumTotalHarmonicDistortion=%5.1f%%)' %
2075 (100.0*thd, 100.0*max_thd)]
2076 # smoothness of spectrum:
2077 if 'dbdiff' in props:
2078 db_diff = props['dbdiff']
2079 msg += ['dBdiff=%5.1fdB' % db_diff]
2080 if max_db_diff > 0.0 and db_diff > max_db_diff:
2081 remove = True
2082 skip_reason += ['not smooth s.d. diff=%5.1fdB (maxximumPowerDifference=%5.1fdB)' %
2083 (db_diff, max_db_diff)]
2084 # maximum power of higher harmonics:
2085 if 'maxdb' in props:
2086 max_harmonics = props['maxdb']
2087 msg += ['max harmonics=%5.1fdB' % max_harmonics]
2088 if max_harmonics > max_harmonics_db:
2089 remove = True
2090 skip_reason += ['maximum harmonics=%5.1fdB too strong (maximumHarmonicsPower=%5.1fdB)' %
2091 (max_harmonics, max_harmonics_db)]
2092 # relative amplitude of harmonics:
2093 if harm_relampl is not None:
2094 for k, max_relampl in enumerate([max_relampl_harm1, max_relampl_harm2, max_relampl_harm3]):
2095 if k >= len(harm_relampl):
2096 break
2097 msg += ['ampl%d=%5.1f%%' % (k+1, 100.0*harm_relampl[k])]
2098 if max_relampl > 0.0 and k < len(harm_relampl) and harm_relampl[k] >= max_relampl:
2099 num_str = ['First', 'Second', 'Third']
2100 skip_reason += ['distorted ampl%d=%5.1f%% (maximum%sHarmonicAmplitude=%5.1f%%)' %
2101 (k+1, 100.0*harm_relampl[k], num_str[k], 100.0*max_relampl)]
2102 return remove, ', '.join(skip_reason), ', '.join(msg)
2105def pulse_quality(props, max_clipped_frac=0.1, max_rms_sem=0.0):
2106 """Assess the quality of an EOD waveform of a pulse fish.
2108 Parameters
2109 ----------
2110 props: dict
2111 A dictionary with properties of the analyzed EOD waveform
2112 as returned by `analyze_pulse()`.
2113 max_clipped_frac: float
2114 Maximum allowed fraction of clipped data.
2115 max_rms_sem: float
2116 If not zero, maximum allowed standard error of the data
2117 relative to p-p amplitude.
2119 Returns
2120 -------
2121 skip_reason: string
2122 An empty string if the waveform is good, otherwise a string
2123 indicating the failure.
2124 msg: string
2125 A textual representation of the values tested.
2126 skipped_clipped: bool
2127 True if waveform was skipped because of clipping.
2128 """
2129 msg = []
2130 skip_reason = []
2131 skipped_clipped = False
2132 # clipped fraction:
2133 if 'clipped' in props:
2134 clipped_frac = props['clipped']
2135 msg += ['clipped=%3.0f%%' % (100.0*clipped_frac)]
2136 if clipped_frac >= max_clipped_frac:
2137 skip_reason += ['clipped=%3.0f%% (maximumClippedFraction=%3.0f%%)' %
2138 (100.0*clipped_frac, 100.0*max_clipped_frac)]
2139 skipped_clipped = True
2140 # noise:
2141 rms_sem = None
2142 if 'rmssem' in props:
2143 rms_sem = props['rmssem']
2144 if 'noise' in props:
2145 rms_sem = props['noise']
2146 if rms_sem is not None:
2147 msg += ['rms sem waveform=%6.2f%%' % (100.0*rms_sem)]
2148 if max_rms_sem > 0.0 and rms_sem >= max_rms_sem:
2149 skip_reason += ['noisy waveform s.e.m.=%6.2f%% (maximumRMSNoise=%6.2f%%)' %
2150 (100.0*rms_sem, 100.0*max_rms_sem)]
2151 return ', '.join(skip_reason), ', '.join(msg), skipped_clipped
2154def plot_eod_recording(ax, data, rate, unit=None, width=0.1,
2155 toffs=0.0, pstyle=dict(lw=2, color='tab:red')):
2156 """Plot a zoomed in range of the recorded trace.
2158 Parameters
2159 ----------
2160 ax: matplotlib axes
2161 Axes used for plotting.
2162 data: 1D ndarray
2163 Recorded data to be plotted.
2164 rate: float
2165 Sampling rate of the data in Hertz.
2166 unit: string
2167 Optional unit of the data used for y-label.
2168 width: float
2169 Width of data segment to be plotted in seconds.
2170 toffs: float
2171 Time of first data value in seconds.
2172 pstyle: dict
2173 Arguments passed on to the plot command for the recorded trace.
2174 """
2175 widx2 = int(width*rate)//2
2176 i0 = len(data)//2 - widx2
2177 i0 = (i0//widx2)*widx2
2178 i1 = i0 + 2*widx2
2179 if i0 < 0:
2180 i0 = 0
2181 if i1 >= len(data):
2182 i1 = len(data)
2183 time = np.arange(len(data))/rate + toffs
2184 tunit = 'sec'
2185 if np.abs(time[i0]) < 1.0 and np.abs(time[i1]) < 1.0:
2186 time *= 1000.0
2187 tunit = 'ms'
2188 ax.plot(time, data, **pstyle)
2189 ax.set_xlim(time[i0], time[i1])
2191 ax.set_xlabel('Time [%s]' % tunit)
2192 ymin = np.min(data[i0:i1])
2193 ymax = np.max(data[i0:i1])
2194 dy = ymax - ymin
2195 ax.set_ylim(ymin-0.05*dy, ymax+0.05*dy)
2196 if len(unit) == 0 or unit == 'a.u.':
2197 ax.set_ylabel('Amplitude')
2198 else:
2199 ax.set_ylabel('Amplitude [%s]' % unit)
2202def plot_pulse_eods(ax, data, rate, zoom_window, width, eod_props,
2203 toffs=0.0, colors=None, markers=None, marker_size=10,
2204 legend_rows=8, **kwargs):
2205 """Mark pulse EODs in a plot of an EOD recording.
2207 Parameters
2208 ----------
2209 ax: matplotlib axes
2210 Axes used for plotting.
2211 data: 1D ndarray
2212 Recorded data (these are not plotted!).
2213 rate: float
2214 Sampling rate of the data in Hertz.
2215 zoom_window: tuple of floats
2216 Start and end time of the data to be plotted in seconds.
2217 width: float
2218 Minimum width of the data to be plotted in seconds.
2219 eod_props: list of dictionaries
2220 Lists of EOD properties as returned by `analyze_pulse()`
2221 and `analyze_wave()`. From the entries with 'type' ==
2222 'pulse' the properties 'EODf' and 'times' are used. 'EODf'
2223 is the averaged EOD frequency, and 'times' is a list of
2224 detected EOD pulse times.
2225 toffs: float
2226 Time of first data value in seconds that will be added
2227 to the pulse times in `eod_props`.
2228 colors: list of colors or None
2229 If not None list of colors for plotting each cluster
2230 markers: list of markers or None
2231 If not None list of markers for plotting each cluster
2232 marker_size: float
2233 Size of markers used to mark the pulses.
2234 legend_rows: int
2235 Maximum number of rows to be used for the legend.
2236 kwargs:
2237 Key word arguments for the legend of the plot.
2238 """
2239 k = 0
2240 for eod in eod_props:
2241 if eod['type'] != 'pulse':
2242 continue
2243 if 'times' not in eod:
2244 continue
2246 width = np.min([width, np.diff(zoom_window)[0]])
2247 while len(eod['peaktimes'][(eod['peaktimes']>(zoom_window[1]-width)) & (eod['peaktimes']<(zoom_window[1]))]) == 0:
2248 width *= 2
2249 if zoom_window[1] - width < 0:
2250 width = width/2
2251 break
2253 x = eod['peaktimes'] + toffs
2254 pidx = np.round(eod['peaktimes']*rate).astype(int)
2255 y = data[pidx[(pidx>0)&(pidx<len(data))]]
2256 x = x[(pidx>0)&(pidx<len(data))]
2257 color_kwargs = {}
2258 #if colors is not None:
2259 # color_kwargs['color'] = colors[k%len(colors)]
2260 if marker_size is not None:
2261 color_kwargs['ms'] = marker_size
2262 label = '%6.1f Hz' % eod['EODf']
2263 if legend_rows > 5 and k >= legend_rows:
2264 label = None
2265 if markers is None:
2266 ax.plot(x, y, 'o', label=label, zorder=-1, **color_kwargs)
2267 else:
2268 ax.plot(x, y, linestyle='none', label=label,
2269 marker=markers[k%len(markers)], mec=None, mew=0.0,
2270 zorder=-1, **color_kwargs)
2271 k += 1
2273 # legend:
2274 if k > 1:
2275 if legend_rows > 0:
2276 if legend_rows > 5:
2277 ncol = 1
2278 else:
2279 ncol = (len(idx)-1) // legend_rows + 1
2280 ax.legend(numpoints=1, ncol=ncol, **kwargs)
2281 else:
2282 ax.legend(numpoints=1, **kwargs)
2284 # reset window so at least one EOD of each cluster is visible
2285 if len(zoom_window)>0:
2286 ax.set_xlim(max(toffs,toffs+zoom_window[1]-width), toffs+zoom_window[1])
2288 i0 = max(0,int((zoom_window[1]-width)*rate))
2289 i1 = int(zoom_window[1]*rate)
2291 ymin = np.min(data[i0:i1])
2292 ymax = np.max(data[i0:i1])
2293 dy = ymax - ymin
2294 ax.set_ylim(ymin-0.05*dy, ymax+0.05*dy)
2297def plot_eod_snippets(ax, data, rate, tmin, tmax, eod_times,
2298 n_snippets=10, flip=False,
2299 sstyle=dict(scaley=False,
2300 lw=0.5, color='0.6')):
2301 """Plot a few EOD waveform snippets.
2303 Parameters
2304 ----------
2305 ax: matplotlib axes
2306 Axes used for plotting.
2307 data: 1D ndarray
2308 Recorded data from which the snippets are taken.
2309 rate: float
2310 Sampling rate of the data in Hertz.
2311 tmin: float
2312 Start time of each snippet.
2313 tmax: float
2314 End time of each snippet.
2315 eod_times: 1-D array
2316 EOD peak times from which a few are selected to be plotted.
2317 n_snippets: int
2318 Number of snippets to be plotted. If zero do not plot anything.
2319 flip: bool
2320 If True flip the snippets upside down.
2321 sstyle: dict
2322 Arguments passed on to the plot command for plotting the snippets.
2323 """
2324 if data is None or n_snippets <= 0:
2325 return
2326 i0 = int(tmin*rate)
2327 i1 = int(tmax*rate)
2328 time = 1000.0*np.arange(i0, i1)/rate
2329 step = len(eod_times)//n_snippets
2330 if step < 1:
2331 step = 1
2332 for t in eod_times[n_snippets//2::step]:
2333 idx = int(np.round(t*rate))
2334 if idx + i0 < 0 or idx + i1 >= len(data):
2335 continue
2336 snippet = data[idx + i0:idx + i1]
2337 if flip:
2338 snippet *= -1
2339 ax.plot(time, snippet - np.mean(snippet[:len(snippet)//4]),
2340 zorder=-5, **sstyle)
2343def plot_eod_waveform(ax, eod_waveform, props, phases=None,
2344 unit=None, tfac=1,
2345 mstyle=dict(lw=2, color='tab:red'),
2346 pstyle=dict(facecolor='tab:green', alpha=0.2,
2347 edgecolor='none'),
2348 nstyle=dict(facecolor='tab:blue', alpha=0.2,
2349 edgecolor='none'),
2350 sstyle=dict(color='0.8'),
2351 fstyle=dict(lw=3, color='tab:blue'),
2352 zstyle=dict(lw=1, color='0.7')):
2353 """Plot mean EOD, its standard error, and an optional fit to the EOD.
2355 Parameters
2356 ----------
2357 ax: matplotlib axes
2358 Axes used for plotting.
2359 eod_waveform: 2-D array
2360 EOD waveform. First column is time in seconds, second column
2361 the (mean) eod waveform. The optional third column is the
2362 standard error, the optional fourth column is a fit of the
2363 whole waveform, and the optional fourth column is a fit of
2364 the tails of a pulse waveform.
2365 props: dict
2366 A dictionary with properties of the analyzed EOD waveform as
2367 returned by `analyze_wave()` and `analyze_pulse()`.
2368 phases: 2_D arrays or None
2369 List of phase properties (index, time, amplitude,
2370 relative amplitude, width, area) of a EOD pulse
2371 as returned by `analyze_pulse()`.
2372 unit: string
2373 Optional unit of the data used for y-label.
2374 tfac: float
2375 Factor scaling the time axis limits.
2376 mstyle: dict
2377 Arguments passed on to the plot command for the mean EOD.
2378 pstyle: dict
2379 Arguments passed on to the fill_between command for coloring
2380 positive phases.
2381 nstyle: dict
2382 Arguments passed on to the fill_between command for coloring
2383 negative phases.
2384 sstyle: dict
2385 Arguments passed on to the fill_between command for the
2386 standard error of the EOD.
2387 fstyle: dict
2388 Arguments passed on to the plot command for the fitted EOD.
2389 zstyle: dict
2390 Arguments passed on to the plot command for the zero line.
2391 """
2392 ax.autoscale(True)
2393 time = 1000 * eod_waveform[:, 0]
2394 mean_eod = eod_waveform[:, 1]
2395 # plot zero line:
2396 ax.axhline(0.0, zorder=-5, **zstyle)
2397 # plot areas:
2398 if phases is not None and len(phases) > 0:
2399 if pstyle is not None and len(pstyle) > 0:
2400 ax.fill_between(time, mean_eod, 0, mean_eod >= 0, zorder=4,
2401 **pstyle)
2402 if nstyle is not None and len(nstyle) > 0:
2403 ax.fill_between(time, mean_eod, 0, mean_eod <= 0, zorder=4,
2404 **nstyle)
2405 # plot fits:
2406 if eod_waveform.shape[1] > 3 and np.all(np.isfinite(eod_waveform[:, 3])):
2407 ax.plot(time, eod_waveform[:, 3], zorder=4, **fstyle)
2408 if eod_waveform.shape[1] > 4:
2409 fs = dict(**fstyle)
2410 if 'lw' in fs:
2411 fs['lw'] *= 2
2412 ax.plot(time, eod_waveform[:, 4], zorder=5, **fs)
2413 # plot waveform:
2414 ax.plot(time, mean_eod, zorder=10, **mstyle)
2415 # plot standard error:
2416 if eod_waveform.shape[1] > 2:
2417 std_eod = eod_waveform[:, 2]
2418 if np.mean(std_eod)/(np.max(mean_eod) - np.min(mean_eod)) > 0.1:
2419 ax.autoscale_view(False)
2420 ax.autoscale(False)
2421 ax.fill_between(time, mean_eod + std_eod, mean_eod - std_eod,
2422 zorder=-10, **sstyle)
2423 # ax height dimensions:
2424 pixelx = np.abs(np.diff(ax.get_window_extent().get_points()[:, 0]))[0]
2425 dxu = (time[-1] - time[0])/pixelx
2426 xfs = plt.rcParams['font.size']*dxu
2427 pixely = np.abs(np.diff(ax.get_window_extent().get_points()[:, 1]))[0]
2428 ymin, ymax = ax.get_ylim()
2429 unity = ymax - ymin
2430 dyu = np.abs(unity)/pixely
2431 yfs = plt.rcParams['font.size']*dyu
2432 # annotate fit:
2433 tau = None if props is None else props.get('tau', None)
2434 ty = 0.0
2435 if tau is not None and eod_waveform.shape[1] > 4:
2436 if tau < 0.001:
2437 label = f'\u03c4={1.e6*tau:.0f}\u00b5s'
2438 else:
2439 label = f'\u03c4={1.e3*tau:.2f}ms'
2440 inx = np.argmin(np.isnan(eod_waveform[:, 4]))
2441 x = eod_waveform[inx, 0] + 1.5*tau
2442 ty = 0.7*eod_waveform[inx, 4]
2443 if np.abs(ty) < 0.5*yfs:
2444 ty = 0.5*yfs*np.sign(ty)
2445 va = 'bottom' if ty > 0.0 else 'top'
2446 ax.text(1000*x, ty, label, ha='left', va=va, zorder=20)
2447 # annotate phases:
2448 if phases is not None and len(phases) > 0:
2449 for i, p in enumerate(phases):
2450 ax.plot(1000*p[1], p[2], 'o', clip_on=False, zorder=0,
2451 alpha=0.4, color=mstyle['color'], ms=12,
2452 mec='none', mew=0)
2453 label = f'P{p[0]:.0f}'
2454 if p[0] != 1:
2455 if np.abs(p[1]) < 0.001:
2456 ts = f'{1.0e6*p[1]:.0f}\u00b5s'
2457 elif np.abs(p[1]) < 0.01:
2458 ts = f'{1.0e3*p[1]:.2f}ms'
2459 else:
2460 ts = f'{1.0e3*p[1]:.3g}ms'
2461 if np.abs(p[3]) < 0.05:
2462 ps = f'{100*p[3]:.1f}%'
2463 else:
2464 ps = f'{100*p[3]:.0f}%'
2465 label += f'({ps} @ {ts})'
2466 left_text = (p[2] > 0 and p[1] >= 0.0) or \
2467 (p[2] < 0 and p[1] >= phases[np.argmin(phases[:, 2]), 1])
2468 va = 'baseline'
2469 dy = 0.6*yfs
2470 sign = np.sign(p[2])
2471 if i > 0 and i < len(phases) - 1:
2472 if phases[i - 1][2] > p[2] and phases[i + 1][2] > p[2]:
2473 va = 'top'
2474 dy = -dy
2475 elif sign < 0:
2476 va = 'top'
2477 dy = -dy
2478 if p[0] == 1:
2479 dy = 0.0
2480 if p[2] + dy < ymin + 1.3*yfs:
2481 dy = ymin + 1.3*yfs - p[2]
2482 if p[0] == np.max(phases[:, 0]) and ty*p[2] > 0.0 and \
2483 sign*p[2] + dy < sign*ty + 1.2*yfs:
2484 dy = ty + sign*1.2*yfs - p[2]
2485 dx = 0.5*xfs
2486 if left_text:
2487 ax.text(1000*p[1] + dx, p[2] + dy, label,
2488 ha='left', va=va, zorder=20)
2489 else:
2490 ax.text(1000*p[1] - dx, p[2] + dy, label,
2491 ha='right', va=va, zorder=20)
2492 # area:
2493 if len(p) > 6:
2494 if np.abs(p[6]) < 1e-8:
2495 continue
2496 elif np.abs(p[6]) < 0.05:
2497 label = f'{100*p[6]:.1f}%'
2498 else:
2499 label = f'{100*p[6]:.0f}%'
2500 x = 1000*p[1]
2501 if i > 0 and i < len(phases) - 1:
2502 xl = 1000*phases[i - 1][1]
2503 xr = 1000*phases[i + 1][1]
2504 tsnippet = time[(time > xl) & (time < xr)]
2505 snippet = mean_eod[(time > xl) & (time < xr)]
2506 tsnippet = tsnippet[np.sign(p[2])*snippet > 0]
2507 snippet = snippet[np.sign(p[2])*snippet > 0]
2508 xc = np.sum(tsnippet*snippet)/np.sum(snippet)
2509 x = xc
2510 if abs(p[3]) > 0.5:
2511 ax.text(x, sign*0.6*yfs, label,
2512 rotation='vertical',
2513 va='top' if sign < 0 else 'bottom',
2514 ha='center', zorder=20)
2515 elif abs(p[3]) > 0.25:
2516 ax.text(x, sign*0.6*yfs, label,
2517 va='top' if sign < 0 else 'baseline',
2518 ha='center', zorder=20)
2519 else:
2520 ax.text(x, -sign*0.4*yfs, label,
2521 va='baseline' if sign < 0 else 'top',
2522 ha='center', zorder=20)
2523 # annotate plot:
2524 if unit is None or len(unit) == 0 or unit == 'a.u.':
2525 unit = ''
2526 if props is not None:
2527 label = f'p-p amplitude = {props["p-p-amplitude"]:.3g} {unit}\n'
2528 if 'n' in props:
2529 eods = 'EODs' if props['n'] > 1 else 'EOD'
2530 label += f'n = {props["n"]} {eods}\n'
2531 if 'flipped' in props and props['flipped']:
2532 label += 'flipped\n'
2533 if 'polaritybalance' in props:
2534 label += f'PB={100*props["polaritybalance"]:.0f} %\n'
2535 if -eod_waveform[0,0] < 0.6*eod_waveform[-1,0]:
2536 ax.text(0.97, 0.97, label, transform=ax.transAxes,
2537 va='top', ha='right', zorder=20)
2538 else:
2539 ax.text(0.03, 0.97, label, transform=ax.transAxes,
2540 va='top', zorder=20)
2541 # axis:
2542 if props is not None and props['type'] == 'wave':
2543 lim = tfac*1000.0/props['EODf']
2544 ax.set_xlim([-lim, +lim])
2545 else:
2546 ax.set_xlim(tfac*time[0], tfac*time[-1])
2547 ax.set_xlabel('Time [msec]')
2548 if unit:
2549 ax.set_ylabel(f'Amplitude [{unit}]')
2550 else:
2551 ax.set_ylabel('Amplitude')
2554def plot_wave_spectrum(axa, axp, spec, props, unit=None,
2555 mstyle=dict(color='tab:blue', markersize=10),
2556 sstyle=dict(color='tab:blue', alpha=0.5, lw=2)):
2557 """Plot and annotate spectrum of wave EOD.
2559 Parameters
2560 ----------
2561 axa: matplotlib axes
2562 Axes for amplitude plot.
2563 axp: matplotlib axes
2564 Axes for phase plot.
2565 spec: 2-D array
2566 The amplitude spectrum of a single pulse as returned by
2567 `analyze_wave()`. First column is the index of the harmonics,
2568 second column its frequency, third column its amplitude,
2569 fourth column its amplitude relative to the fundamental, fifth
2570 column is power of harmonics relative to fundamental in
2571 decibel, and sixth column the phase shift relative to the
2572 fundamental.
2573 props: dict
2574 A dictionary with properties of the analyzed EOD waveform as
2575 returned by `analyze_wave()`.
2576 unit: string
2577 Optional unit of the data used for y-label.
2578 mstyle: dict
2579 Arguments passed on to the stem plot command for the markers.
2580 sstyle: dict
2581 Arguments passed on to the stem plot command for the stem lines.
2582 """
2583 n = min(9, np.sum(np.isfinite(spec[:, 2])))
2584 # amplitudes:
2585 markers, stemlines, _ = axa.stem(spec[:n,0]+1, spec[:n,2], basefmt='none')
2586 plt.setp(markers, clip_on=False, **mstyle)
2587 plt.setp(stemlines, **sstyle)
2588 axa.set_xlim(0.5, n+0.5)
2589 axa.set_ylim(bottom=0)
2590 axa.xaxis.set_major_locator(plt.MultipleLocator(1))
2591 axa.tick_params('x', direction='out')
2592 if unit:
2593 axa.set_ylabel(f'Amplitude [{unit}]')
2594 else:
2595 axa.set_ylabel('Amplitude')
2596 # phases:
2597 phases = spec[:n,5]
2598 phases[phases<0.0] = phases[phases<0.0] + 2.0*np.pi
2599 markers, stemlines, _ = axp.stem(spec[:n,0]+1, phases[:n], basefmt='none')
2600 plt.setp(markers, clip_on=False, **mstyle)
2601 plt.setp(stemlines, **sstyle)
2602 axp.set_xlim(0.5, n+0.5)
2603 axp.xaxis.set_major_locator(plt.MultipleLocator(1))
2604 axp.tick_params('x', direction='out')
2605 axp.set_ylim(0, 2.0*np.pi)
2606 axp.set_yticks([0, np.pi, 2.0*np.pi])
2607 axp.set_yticklabels(['0', '\u03c0', '2\u03c0'])
2608 axp.set_xlabel('Harmonics')
2609 axp.set_ylabel('Phase')
2612def plot_pulse_spectrum(ax, energy, props, min_freq=1.0, max_freq=10000.0,
2613 sstyle=dict(lw=3, color='tab:blue'),
2614 astyle=dict(lw=4, color='tab:cyan'),
2615 pstyle=dict(ls='', marker='o', markersize=10,
2616 color='tab:blue', mec='none', mew=0,
2617 alpha=0.4),
2618 cstyle=dict(lw=1, ls='-', color='0.5'),
2619 att5_color='0.8', att50_color='0.9'):
2620 """Plot and annotate spectrum of single pulse EOD.
2622 Parameters
2623 ----------
2624 ax: matplotlib axes
2625 Axes used for plotting.
2626 energy: 2-D array
2627 The energy spectrum of a single pulse as returned by `analyze_pulse()`.
2628 First column are the frequencies, second column the energy.
2629 An optional third column is an analytically computed spectrum.
2630 props: dict
2631 A dictionary with properties of the analyzed EOD waveform as
2632 returned by `analyze_pulse()`.
2633 min_freq: float
2634 Minimun frequency of the spectrum to be plotted (logscale!).
2635 max_freq: float
2636 Maximun frequency of the spectrum to be plotted (logscale!).
2637 sstyle: dict
2638 Arguments passed on to the plot command for the energy spectrum
2639 computed from the data.
2640 astyle: dict
2641 Arguments passed on to the plot command for the energy spectrum
2642 that was analytically computed from the Gaussian fits
2643 (optional third column in `energy`).
2644 pstyle: dict
2645 Arguments passed on to the plot command for marking the peak frequency.
2646 cstyle: dict
2647 Arguments passed on to the plot command for the line marking
2648 the lower cutoff frequency.
2649 att5_color: matplotlib color specification
2650 Color for the rectangular patch marking the first 5 Hz.
2651 att50_color: matplotlib color specification
2652 Color for the rectangular patch marking the first 50 Hz.
2653 """
2654 ax.axvspan(1, 50, color=att50_color, zorder=10)
2655 att = props['energyatt50']
2656 if att < -10:
2657 ax.text(10, att + 1, f'{att:.0f}dB',
2658 ha='left', va='bottom', zorder=100)
2659 else:
2660 ax.text(10, att - 1, f'{att:.0f}dB',
2661 ha='left', va='top', zorder=100)
2662 ax.axvspan(1, 5, color=att5_color, zorder=20)
2663 att = props['energyatt5']
2664 if att < -10:
2665 ax.text(4, att + 1, f'{att:.0f}dB',
2666 ha='right', va='bottom', zorder=100)
2667 else:
2668 ax.text(4, att - 1, f'{att:.0f}dB',
2669 ha='right', va='top', zorder=100)
2670 lowcutoff = props['lowcutoff']
2671 if lowcutoff >= min_freq:
2672 ax.plot([lowcutoff, lowcutoff, 1], [-60, 0.5*att, 0.5*att],
2673 zorder=30, **cstyle)
2674 ax.text(1.2*lowcutoff, 0.5*att - 1, f'{lowcutoff:.0f}Hz',
2675 ha='left', va='top', zorder=100)
2676 highcutoff = props['highcutoff']
2677 ax.plot([highcutoff, highcutoff], [-60, -3], zorder=30, **cstyle)
2678 ax.text(1.2*highcutoff, -3, f'{highcutoff:.0f}Hz',
2679 ha='left', va='center', zorder=100)
2680 ref_energy = np.max(energy[:, 1])
2681 if energy.shape[1] > 2 and np.all(np.isfinite(energy[:, 2])) and len(astyle) > 0:
2682 db = decibel(energy[:, 2], ref_energy)
2683 ax.plot(energy[:, 0], db, zorder=45, **astyle)
2684 db = decibel(energy[:, 1], ref_energy)
2685 ax.plot(energy[:, 0], db, zorder=50, **sstyle)
2686 peakfreq = props['peakfreq']
2687 if peakfreq >= min_freq:
2688 ax.plot(peakfreq, 0, zorder=60, clip_on=False, **pstyle)
2689 ax.text(peakfreq*1.2, 1, f'{peakfreq:.0f}Hz',
2690 va='bottom', zorder=100)
2691 troughfreq = props['troughfreq']
2692 if troughfreq >= min_freq:
2693 troughenergy = decibel(props['troughenergy'], ref_energy)
2694 ax.plot(troughfreq, troughenergy, zorder=60, **pstyle)
2695 ax.text(troughfreq, troughenergy - 3,
2696 f'{troughenergy:.0f}dB @ {troughfreq:.0f}Hz',
2697 ha='center', va='top', zorder=100)
2698 ax.set_xlim(min_freq, max_freq)
2699 ax.set_xscale('log')
2700 ax.set_ylim(-60, 2)
2701 ax.set_xlabel('Frequency [Hz]')
2702 ax.set_ylabel('Energy [dB]')
2705def save_eod_waveform(mean_eod, unit, idx, basename, **kwargs):
2706 """Save mean EOD waveform to file.
2708 Parameters
2709 ----------
2710 mean_eod: 2D array of floats
2711 Averaged EOD waveform as returned by `eod_waveform()`,
2712 `analyze_wave()`, and `analyze_pulse()`.
2713 unit: string
2714 Unit of the waveform data.
2715 idx: int or None
2716 Index of fish.
2717 basename: string or stream
2718 If string, path and basename of file.
2719 If `basename` does not have an extension,
2720 '-eodwaveform', the fish index, and a file extension are appended.
2721 If stream, write EOD waveform data into this stream.
2722 kwargs:
2723 Arguments passed on to `TableData.write()`.
2725 Returns
2726 -------
2727 filename: Path
2728 Path and full name of the written file in case of `basename`
2729 being a string. Otherwise, the file name and extension that
2730 would have been appended to a basename.
2732 See Also
2733 --------
2734 load_eod_waveform()
2735 """
2736 td = TableData(mean_eod[:, :3]*[1000.0, 1.0, 1.0],
2737 ['time', 'mean', 'sem'],
2738 ['ms', unit, unit],
2739 ['%.3f', '%.6g', '%.6g'])
2740 if mean_eod.shape[1] > 3:
2741 td.append('fit', unit, '%.5f', value=mean_eod[:, 3])
2742 if mean_eod.shape[1] > 4:
2743 td.append('tailfit', unit, '%.5f', value=mean_eod[:, 4])
2744 fp = ''
2745 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
2746 if not ext:
2747 fp = '-eodwaveform'
2748 if idx is not None:
2749 fp += f'-{idx}'
2750 return td.write_file_stream(basename, fp, **kwargs)
2753def load_eod_waveform(file_path):
2754 """Load EOD waveform from file.
2756 Parameters
2757 ----------
2758 file_path: string
2759 Path of the file to be loaded.
2761 Returns
2762 -------
2763 mean_eod: 2D array of floats
2764 Averaged EOD waveform: time in seconds, mean, standard deviation, fit.
2765 unit: string
2766 Unit of EOD waveform.
2768 Raises
2769 ------
2770 FileNotFoundError:
2771 If `file_path` does not exist.
2773 See Also
2774 --------
2775 save_eod_waveform()
2776 """
2777 data = TableData(file_path)
2778 mean_eod = data.array()
2779 mean_eod[:, 0] *= 0.001
2780 return mean_eod, data.unit('mean')
2783def save_wave_eodfs(wave_eodfs, wave_indices, basename, **kwargs):
2784 """Save frequencies of wave EODs to file.
2786 Parameters
2787 ----------
2788 wave_eodfs: list of 2D arrays
2789 Each item is a matrix with the frequencies and powers
2790 (columns) of the fundamental and harmonics (rows) as returned
2791 by `harmonics.harmonic_groups()`.
2792 wave_indices: array
2793 Indices identifying each fish or NaN.
2794 If None no index column is inserted.
2795 basename: string or stream
2796 If string, path and basename of file.
2797 If `basename` does not have an extension,
2798 '-waveeodfs' and a file extension according to `kwargs` are appended.
2799 If stream, write EOD frequencies data into this stream.
2800 kwargs:
2801 Arguments passed on to `TableData.write()`.
2803 Returns
2804 -------
2805 filename: Path
2806 Path and full name of the written file in case of `basename`
2807 being a string. Otherwise, the file name and extension that
2808 would have been appended to a basename.
2810 See Also
2811 --------
2812 load_wave_eodfs()
2814 """
2815 eodfs = fundamental_freqs_and_power(wave_eodfs)
2816 td = TableData()
2817 if wave_indices is not None:
2818 td.append('index', '', '%d',
2819 value=[wi if wi >= 0 else np.nan
2820 for wi in wave_indices])
2821 td.append('EODf', 'Hz', '%7.2f', value=eodfs[:, 0])
2822 td.append('datapower', 'dB', '%7.2f', value=eodfs[:, 1])
2823 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
2824 fp = '-waveeodfs' if not ext else ''
2825 return td.write_file_stream(basename, fp, **kwargs)
2828def load_wave_eodfs(file_path):
2829 """Load frequencies of wave EODs from file.
2831 Parameters
2832 ----------
2833 file_path: string
2834 Path of the file to be loaded.
2836 Returns
2837 -------
2838 eodfs: 2D array of floats
2839 EODfs and power of wave type fish.
2840 indices: array of ints
2841 Corresponding indices of fish, can contain negative numbers to
2842 indicate frequencies without fish.
2844 Raises
2845 ------
2846 FileNotFoundError:
2847 If `file_path` does not exist.
2849 See Also
2850 --------
2851 save_wave_eodfs()
2852 """
2853 data = TableData(file_path)
2854 eodfs = data.array()
2855 if 'index' in data:
2856 indices = data[:, 'index']
2857 indices[~np.isfinite(indices)] = -1
2858 indices = np.array(indices, dtype=int)
2859 eodfs = eodfs[:, 1:]
2860 else:
2861 indices = np.zeros(data.rows(), dtype=int) - 1
2862 return eodfs, indices
2865def save_wave_fish(eod_props, unit, basename, **kwargs):
2866 """Save properties of wave EODs to file.
2868 Parameters
2869 ----------
2870 eod_props: list of dict
2871 Properties of EODs as returned by `analyze_wave()` and
2872 `analyze_pulse()`. Only properties of wave fish are saved.
2873 unit: string
2874 Unit of the waveform data.
2875 basename: string or stream
2876 If string, path and basename of file.
2877 If `basename` does not have an extension,
2878 '-wavefish' and a file extension are appended.
2879 If stream, write wave fish properties into this stream.
2880 kwargs:
2881 Arguments passed on to `TableData.write()`.
2883 Returns
2884 -------
2885 filename: Path or None
2886 Path and full name of the written file in case of `basename`
2887 being a string. Otherwise, the file name and extension that
2888 would have been appended to a basename.
2889 None if no wave fish are contained in eod_props and
2890 consequently no file was written.
2892 See Also
2893 --------
2894 load_wave_fish()
2895 """
2896 wave_props = [p for p in eod_props if p['type'] == 'wave']
2897 if len(wave_props) == 0:
2898 return None
2899 td = TableData()
2900 if 'twin' in wave_props[0] or 'rate' in wave_props[0] or \
2901 'nfft' in wave_props[0]:
2902 td.append_section('recording')
2903 if 'twin' in wave_props[0]:
2904 td.append('twin', 's', '%7.2f', value=wave_props)
2905 td.append('window', 's', '%7.2f', value=wave_props)
2906 td.append('winclipped', '%', '%.2f',
2907 value=wave_props, fac=100.0)
2908 if 'samplerate' in wave_props[0]:
2909 td.append('samplerate', 'kHz', '%.3f',
2910 value=wave_props, fac=0.001)
2911 if 'nfft' in wave_props[0]:
2912 td.append('nfft', '', '%d', value=wave_props)
2913 td.append('dfreq', 'Hz', '%.2f', value=wave_props)
2914 td.append_section('waveform')
2915 td.append('index', '', '%d', value=wave_props)
2916 td.append('EODf', 'Hz', '%7.2f', value=wave_props)
2917 td.append('p-p-amplitude', unit, '%.5f', value=wave_props)
2918 td.append('power', 'dB', '%7.2f', value=wave_props)
2919 if 'datapower' in wave_props[0]:
2920 td.append('datapower', 'dB', '%7.2f', value=wave_props)
2921 td.append('thd', '%', '%.2f', value=wave_props, fac=100.0)
2922 td.append('dbdiff', 'dB', '%7.2f', value=wave_props)
2923 td.append('maxdb', 'dB', '%7.2f', value=wave_props)
2924 if 'noise' in wave_props[0]:
2925 td.append('noise', '%', '%.1f', value=wave_props, fac=100.0)
2926 td.append('rmserror', '%', '%.2f', value=wave_props, fac=100.0)
2927 if 'clipped' in wave_props[0]:
2928 td.append('clipped', '%', '%.1f', value=wave_props, fac=100.0)
2929 td.append('flipped', '', '%d', value=wave_props)
2930 td.append('n', '', '%5d', value=wave_props)
2931 td.append_section('timing')
2932 td.append('ncrossings', '', '%d', value=wave_props)
2933 td.append('peakwidth', '%', '%.2f', value=wave_props, fac=100.0)
2934 td.append('troughwidth', '%', '%.2f', value=wave_props, fac=100.0)
2935 td.append('minwidth', '%', '%.2f', value=wave_props, fac=100.0)
2936 td.append('leftpeak', '%', '%.2f', value=wave_props, fac=100.0)
2937 td.append('rightpeak', '%', '%.2f', value=wave_props, fac=100.0)
2938 td.append('lefttrough', '%', '%.2f', value=wave_props, fac=100.0)
2939 td.append('righttrough', '%', '%.2f', value=wave_props, fac=100.0)
2940 td.append('p-p-distance', '%', '%.2f', value=wave_props, fac=100.0)
2941 td.append('min-p-p-distance', '%', '%.2f',
2942 value=wave_props, fac=100.0)
2943 td.append('relpeakampl', '%', '%.2f', value=wave_props, fac=100.0)
2944 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
2945 fp = '-wavefish' if not ext else ''
2946 return td.write_file_stream(basename, fp, **kwargs)
2949def load_wave_fish(file_path):
2950 """Load properties of wave EODs from file.
2952 All times are scaled to seconds, all frequencies to Hertz and all
2953 percentages to fractions.
2955 Parameters
2956 ----------
2957 file_path: string
2958 Path of the file to be loaded.
2960 Returns
2961 -------
2962 eod_props: list of dict
2963 Properties of EODs.
2965 Raises
2966 ------
2967 FileNotFoundError:
2968 If `file_path` does not exist.
2970 See Also
2971 --------
2972 save_wave_fish()
2974 """
2975 data = TableData(file_path)
2976 eod_props = data.dicts()
2977 for props in eod_props:
2978 if 'winclipped' in props:
2979 props['winclipped'] /= 100
2980 if 'samplerate' in props:
2981 props['samplerate'] *= 1000
2982 if 'nfft' in props:
2983 props['nfft'] = int(props['nfft'])
2984 props['index'] = int(props['index'])
2985 props['n'] = int(props['n'])
2986 props['type'] = 'wave'
2987 props['thd'] /= 100
2988 props['noise'] /= 100
2989 props['rmserror'] /= 100
2990 if 'clipped' in props:
2991 props['clipped'] /= 100
2992 props['ncrossings'] = int(props['ncrossings'])
2993 props['peakwidth'] /= 100
2994 props['troughwidth'] /= 100
2995 props['minwidth'] /= 100
2996 props['leftpeak'] /= 100
2997 props['rightpeak'] /= 100
2998 props['lefttrough'] /= 100
2999 props['righttrough'] /= 100
3000 props['p-p-distance'] /= 100
3001 props['min-p-p-distance'] /= 100
3002 props['relpeakampl'] /= 100
3003 return eod_props
3006def save_pulse_fish(eod_props, unit, basename, **kwargs):
3007 """Save properties of pulse EODs to file.
3009 Parameters
3010 ----------
3011 eod_props: list of dict
3012 Properties of EODs as returned by `analyze_wave()` and
3013 `analyze_pulse()`. Only properties of pulse fish are saved.
3014 unit: string
3015 Unit of the waveform data.
3016 basename: string or stream
3017 If string, path and basename of file.
3018 If `basename` does not have an extension,
3019 '-pulsefish' and a file extension are appended.
3020 If stream, write pulse fish properties into this stream.
3021 kwargs:
3022 Arguments passed on to `TableData.write()`.
3024 Returns
3025 -------
3026 filename: Path or None
3027 Path and full name of the written file in case of `basename`
3028 being a string. Otherwise, the file name and extension that
3029 would have been appended to a basename.
3030 None if no pulse fish are contained in eod_props and
3031 consequently no file was written.
3033 See Also
3034 --------
3035 load_pulse_fish()
3036 """
3037 pulse_props = [p for p in eod_props if p['type'] == 'pulse']
3038 if len(pulse_props) == 0:
3039 return None
3040 td = TableData()
3041 if 'twin' in pulse_props[0] or 'samplerate' in pulse_props[0] or \
3042 'nfft' in pulse_props[0]:
3043 td.append_section('recording')
3044 if 'twin' in pulse_props[0]:
3045 td.append('twin', 's', '%7.2f', value=pulse_props)
3046 td.append('window', 's', '%7.2f', value=pulse_props)
3047 td.append('winclipped', '%', '%.2f',
3048 value=pulse_props, fac=100.0)
3049 if 'samplerate' in pulse_props[0]:
3050 td.append('samplerate', 'kHz', '%.3f', value=pulse_props,
3051 fac=0.001)
3052 if 'nfft' in pulse_props[0]:
3053 td.append('nfft', '', '%d', value=pulse_props)
3054 td.append('dfreq', 'Hz', '%.2f', value=pulse_props)
3055 td.append_section('waveform')
3056 td.append('index', '', '%d', value=pulse_props)
3057 td.append('EODf', 'Hz', '%7.2f', value=pulse_props)
3058 td.append('period', 'ms', '%7.2f', value=pulse_props, fac=1000.0)
3059 td.append('max-ampl', unit, '%.5f', value=pulse_props)
3060 td.append('min-ampl', unit, '%.5f', value=pulse_props)
3061 td.append('p-p-amplitude', unit, '%.5f', value=pulse_props)
3062 if 'noise' in pulse_props[0]:
3063 td.append('noise', '%', '%.2f', value=pulse_props, fac=100.0)
3064 if 'clipped' in pulse_props[0]:
3065 td.append('clipped', '%', '%.2f', value=pulse_props, fac=100.0)
3066 td.append('flipped', '', '%d', value=pulse_props)
3067 td.append('tstart', 'ms', '%.3f', value=pulse_props, fac=1000.0)
3068 td.append('tend', 'ms', '%.3f', value=pulse_props, fac=1000.0)
3069 td.append('width', 'ms', '%.3f', value=pulse_props, fac=1000.0)
3070 td.append('peak-dist', 'ms', '%.3f', value=pulse_props, fac=1000.0)
3071 td.append('tau', 'ms', '%.3f', value=pulse_props, fac=1000.0)
3072 td.append('firstphase', '', '%d', value=pulse_props)
3073 td.append('lastphase', '', '%d', value=pulse_props)
3074 td.append('totalarea', f'{unit}*ms', '%.4f',
3075 value=pulse_props, fac=1000.0)
3076 td.append('positivearea', '%', '%.2f',
3077 value=pulse_props, fac=100.0)
3078 td.append('negativearea', '%', '%.2f',
3079 value=pulse_props, fac=100.0)
3080 td.append('polaritybalance', '%', '%.2f',
3081 value=pulse_props, fac=100.0)
3082 td.append('n', '', '%d', value=pulse_props)
3083 td.append_section('spectrum')
3084 td.append('peakfreq', 'Hz', '%.2f', value=pulse_props)
3085 td.append('peakenergy', f'{unit}^2s/Hz', '%.3g', value=pulse_props)
3086 td.append('troughfreq', 'Hz', '%.2f', value=pulse_props)
3087 td.append('troughenergy', f'{unit}^2s/Hz', '%.3g', value=pulse_props)
3088 td.append('energyatt5', 'dB', '%.2f', value=pulse_props)
3089 td.append('energyatt50', 'dB', '%.2f', value=pulse_props)
3090 td.append('lowcutoff', 'Hz', '%.2f', value=pulse_props)
3091 td.append('highcutoff', 'Hz', '%.2f', value=pulse_props)
3092 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
3093 fp = '-pulsefish' if not ext else ''
3094 return td.write_file_stream(basename, fp, **kwargs)
3097def load_pulse_fish(file_path):
3098 """Load properties of pulse EODs from file.
3100 All times are scaled to seconds, all frequencies to Hertz, and all
3101 percentages to fractions.
3103 Parameters
3104 ----------
3105 file_path: string
3106 Path of the file to be loaded.
3108 Returns
3109 -------
3110 eod_props: list of dict
3111 Properties of EODs.
3113 Raises
3114 ------
3115 FileNotFoundError:
3116 If `file_path` does not exist.
3118 See Also
3119 --------
3120 save_pulse_fish()
3122 """
3123 data = TableData(file_path)
3124 eod_props = data.dicts()
3125 for props in eod_props:
3126 if 'winclipped' in props:
3127 props['winclipped'] /= 100
3128 if 'samplerate' in props:
3129 props['samplerate'] *= 1000
3130 if 'nfft' in props:
3131 props['nfft'] = int(props['nfft'])
3132 props['index'] = int(props['index'])
3133 props['n'] = int(props['n'])
3134 props['firstphase'] = int(props['firstphase'])
3135 props['lastphase'] = int(props['lastphase'])
3136 if 'totalarea' in props:
3137 props['totalarea'] /= 1000
3138 if 'positivearea' in props:
3139 props['positivearea'] /= 100
3140 if 'negativearea' in props:
3141 props['negativearea'] /= 100
3142 if 'polaritybalance' in props:
3143 props['polaritybalance'] /= 100
3144 props['type'] = 'pulse'
3145 if 'clipped' in props:
3146 props['clipped'] /= 100
3147 props['period'] /= 1000
3148 props['noise'] /= 100
3149 props['tstart'] /= 1000
3150 props['tend'] /= 1000
3151 props['width'] /= 1000
3152 props['peak-dist'] /= 1000
3153 props['tau'] /= 1000
3154 return eod_props
3157def save_wave_spectrum(spec_data, unit, idx, basename, **kwargs):
3158 """Save amplitude and phase spectrum of wave EOD to file.
3160 Parameters
3161 ----------
3162 spec_data: 2D array of floats
3163 Amplitude and phase spectrum of wave EOD as returned by
3164 `analyze_wave()`.
3165 unit: string
3166 Unit of the waveform data.
3167 idx: int or None
3168 Index of fish.
3169 basename: string or stream
3170 If string, path and basename of file.
3171 If `basename` does not have an extension,
3172 '-wavespectrum', the fish index, and a file extension are appended.
3173 If stream, write wave spectrum into this stream.
3174 kwargs:
3175 Arguments passed on to `TableData.write()`.
3177 Returns
3178 -------
3179 filename: Path
3180 Path and full name of the written file in case of `basename`
3181 being a string. Otherwise, the file name and extension that
3182 would have been appended to a basename.
3184 See Also
3185 --------
3186 load_wave_spectrum()
3188 """
3189 td = TableData(spec_data[:, :6]*[1.0, 1.0, 1.0, 100.0, 1.0, 1.0],
3190 ['harmonics', 'frequency', 'amplitude', 'relampl', 'relpower', 'phase'],
3191 ['', 'Hz', unit, '%', 'dB', 'rad'],
3192 ['%.0f', '%.2f', '%.6f', '%10.2f', '%6.2f', '%8.4f'])
3193 if spec_data.shape[1] > 6:
3194 td.append('datapower', '%s^2/Hz' % unit, '%11.4e',
3195 value=spec_data[:, 6])
3196 fp = ''
3197 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
3198 if not ext:
3199 fp = '-wavespectrum'
3200 if idx is not None:
3201 fp += f'-{idx}'
3202 return td.write_file_stream(basename, fp, **kwargs)
3205def load_wave_spectrum(file_path):
3206 """Load amplitude and phase spectrum of wave EOD from file.
3208 Parameters
3209 ----------
3210 file_path: string
3211 Path of the file to be loaded.
3213 Returns
3214 -------
3215 spec: 2D array of floats
3216 Amplitude and phase spectrum of wave EOD:
3217 harmonics, frequency, amplitude, relative amplitude in dB,
3218 relative power in dB, phase, data power in unit squared.
3219 Can contain NaNs.
3220 unit: string
3221 Unit of amplitudes.
3223 Raises
3224 ------
3225 FileNotFoundError:
3226 If `file_path` does not exist.
3228 See Also
3229 --------
3230 save_wave_spectrum()
3231 """
3232 data = TableData(file_path)
3233 spec = data.array()
3234 spec[:, 3] *= 0.01
3235 return spec, data.unit('amplitude')
3238def save_pulse_spectrum(spec_data, unit, idx, basename, **kwargs):
3239 """Save energy spectrum of pulse EOD to file.
3241 Parameters
3242 ----------
3243 spec_data: 2D array of floats
3244 Energy spectrum of single pulse as returned by `analyze_pulse()`.
3245 unit: string
3246 Unit of the waveform data.
3247 idx: int or None
3248 Index of fish.
3249 basename: string or stream
3250 If string, path and basename of file.
3251 If `basename` does not have an extension,
3252 '-pulsespectrum', the fish index, and a file extension are appended.
3253 If stream, write pulse spectrum into this stream.
3254 kwargs:
3255 Arguments passed on to `TableData.write()`.
3257 Returns
3258 -------
3259 filename: Path
3260 Path and full name of the written file in case of `basename`
3261 being a string. Otherwise, the file name and extension that
3262 would have been appended to a basename.
3264 See Also
3265 --------
3266 load_pulse_spectrum()
3267 """
3268 td = TableData(spec_data[:, :2], ['frequency', 'energy'],
3269 ['Hz', '%s^2s/Hz' % unit], ['%.2f', '%.4e'])
3270 fp = ''
3271 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
3272 if not ext:
3273 fp = '-pulsespectrum'
3274 if idx is not None:
3275 fp += f'-{idx}'
3276 return td.write_file_stream(basename, fp, **kwargs)
3279def load_pulse_spectrum(file_path):
3280 """Load energy spectrum of pulse EOD from file.
3282 Parameters
3283 ----------
3284 file_path: string
3285 Path of the file to be loaded.
3287 Returns
3288 -------
3289 spec: 2D array of floats
3290 Energy spectrum of single pulse: frequency, energy
3292 Raises
3293 ------
3294 FileNotFoundError:
3295 If `file_path` does not exist.
3297 See Also
3298 --------
3299 save_pulse_spectrum()
3300 """
3301 data = TableData(file_path)
3302 spec = data.array()
3303 return spec
3306def save_pulse_phases(phase_data, unit, idx, basename, **kwargs):
3307 """Save phase properties of pulse EOD to file.
3309 Parameters
3310 ----------
3311 phase_data: 2D array of floats
3312 Properties of phases of pulse EOD as returned by
3313 `analyze_pulse()`.
3314 unit: string
3315 Unit of the waveform data.
3316 idx: int or None
3317 Index of fish.
3318 basename: string or stream
3319 If string, path and basename of file.
3320 If `basename` does not have an extension,
3321 '-pulsephases', the fish index, and a file extension are appended.
3322 If stream, write pulse phases into this stream.
3323 kwargs:
3324 Arguments passed on to `TableData.write()`.
3326 Returns
3327 -------
3328 filename: Path
3329 Path and full name of the written file in case of `basename`
3330 being a string. Otherwise, the file name and extension that
3331 would have been appended to a basename.
3333 See Also
3334 --------
3335 load_pulse_phases()
3336 """
3337 if len(phase_data) == 0:
3338 return None
3339 td = TableData(phase_data[:, :7]*[1, 1000, 1, 100, 1000, 1000, 100],
3340 ['P', 'time', 'amplitude', 'relampl', 'width', 'area', 'relarea'],
3341 ['', 'ms', unit, '%', 'ms', f'{unit}*ms', '%'],
3342 ['%.0f', '%.3f', '%.5f', '%.2f', '%.3f', '%.4f', '%.2f'])
3343 fp = ''
3344 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
3345 if not ext:
3346 fp = '-pulsephases'
3347 if idx is not None:
3348 fp += f'-{idx}'
3349 return td.write_file_stream(basename, fp, **kwargs)
3352def load_pulse_phases(file_path):
3353 """Load phase properties of pulse EOD from file.
3355 Parameters
3356 ----------
3357 file_path: string
3358 Path of the file to be loaded.
3360 Returns
3361 -------
3362 phase_data: 2D array of floats
3363 Properties of phases of pulse EOD:
3364 P, time, amplitude, relampl, width
3365 unit: string
3366 Unit of phase amplitudes.
3368 Raises
3369 ------
3370 FileNotFoundError:
3371 If `file_path` does not exist.
3373 See Also
3374 --------
3375 save_pulse_phases()
3376 """
3377 data = TableData(file_path)
3378 phases = data.array()
3379 phases[:, 1] *= 0.001
3380 phases[:, 3] *= 0.01
3381 phases[:, 4] *= 0.001
3382 phases[:, 5] *= 0.001
3383 phases[:, 6] *= 0.01
3384 return phases, data.unit('amplitude')
3387def save_pulse_gaussians(pulse_data, unit, idx, basename, **kwargs):
3388 """Save Gaussian phase properties of pulse EOD to file.
3390 Parameters
3391 ----------
3392 pulse: dict
3393 Dictionary with phase times, amplitudes and standard
3394 deviations of Gaussians fitted to the pulse waveform
3395 as returned by `decompose_pulse()` and `analyze_pulse()`.
3396 unit: string
3397 Unit of the waveform data.
3398 idx: int or None
3399 Index of fish.
3400 basename: string or stream
3401 If string, path and basename of file.
3402 If `basename` does not have an extension,
3403 '-pulsegaussians', the fish index, and a file extension are appended.
3404 If stream, write pulse phases into this stream.
3405 kwargs:
3406 Arguments passed on to `TableData.write()`.
3408 Returns
3409 -------
3410 filename: Path
3411 Path and full name of the written file in case of `basename`
3412 being a string. Otherwise, the file name and extension that
3413 would have been appended to a basename.
3415 See Also
3416 --------
3417 load_pulse_gaussians()
3418 """
3419 if len(pulse_data) == 0:
3420 return None
3421 td = TableData(pulse_data,
3422 units=['ms', unit, 'ms'],
3423 formats=['%.3f', '%.5f', '%.3f'])
3424 td[:, 'times'] *= 1000
3425 td[:, 'stdevs'] *= 1000
3426 fp = ''
3427 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
3428 if not ext:
3429 fp = '-pulsegaussians'
3430 if idx is not None:
3431 fp += f'-{idx}'
3432 return td.write_file_stream(basename, fp, **kwargs)
3435def load_pulse_gaussians(file_path):
3436 """Load Gaussian phase properties of pulse EOD from file.
3438 Parameters
3439 ----------
3440 file_path: string
3441 Path of the file to be loaded.
3443 Returns
3444 -------
3445 pulse: dict
3446 Dictionary with phase times, amplitudes and standard
3447 deviations of Gaussians fitted to the pulse waveform.
3448 Use the functions provided in thunderfish.fakefish to simulate
3449 pulse fish EODs from this data.
3450 unit: string
3451 Unit of Gaussian amplitudes.
3453 Raises
3454 ------
3455 FileNotFoundError:
3456 If `file_path` does not exist.
3458 See Also
3459 --------
3460 save_pulse_gaussians()
3462 """
3463 data = TableData(file_path)
3464 pulses = data.dict()
3465 pulses['times'] = np.asarray(pulses['times'])
3466 pulses['amplitudes'] = np.asarray(pulses['amplitudes'])
3467 pulses['stdevs'] = np.asarray(pulses['stdevs'])
3468 pulses['times'] *= 0.001
3469 pulses['stdevs'] *= 0.001
3470 return pulses, data.unit('amplitudes')
3473def save_pulse_times(pulse_times, idx, basename, **kwargs):
3474 """Save times of pulse EOD to file.
3476 Parameters
3477 ----------
3478 pulse_times: dict or array of floats
3479 Times of EOD pulses. Either as array of times or
3480 `props['peaktimes']` or `props['times']` as returned by
3481 `analyze_pulse()`.
3482 idx: int or None
3483 Index of fish.
3484 basename: string or stream
3485 If string, path and basename of file.
3486 If `basename` does not have an extension,
3487 '-pulsetimes', the fish index, and a file extension are appended.
3488 If stream, write pulse times into this stream.
3489 kwargs:
3490 Arguments passed on to `TableData.write()`.
3492 Returns
3493 -------
3494 filename: Path
3495 Path and full name of the written file in case of `basename`
3496 being a string. Otherwise, the file name and extension that
3497 would have been appended to a basename.
3499 See Also
3500 --------
3501 load_pulse_times()
3502 """
3503 if isinstance(pulse_times, dict):
3504 props = pulse_times
3505 pulse_times = props.get('times', [])
3506 pulse_times = props.get('peaktimes', pulse_times)
3507 if len(pulse_times) == 0:
3508 return None
3509 td = TableData()
3510 td.append('time', 's', '%.4f', value=pulse_times)
3511 fp = ''
3512 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
3513 if not ext:
3514 fp = '-pulsetimes'
3515 if idx is not None:
3516 fp += f'-{idx}'
3517 return td.write_file_stream(basename, fp, **kwargs)
3520def load_pulse_times(file_path):
3521 """Load times of pulse EOD from file.
3523 Parameters
3524 ----------
3525 file_path: string
3526 Path of the file to be loaded.
3528 Returns
3529 -------
3530 pulse_times: array of floats
3531 Times of pulse EODs in seconds.
3533 Raises
3534 ------
3535 FileNotFoundError:
3536 If `file_path` does not exist.
3538 See Also
3539 --------
3540 save_pulse_times()
3541 """
3542 data = TableData(file_path)
3543 pulse_times = data.array()[:, 0]
3544 return pulse_times
3547file_types = ['waveeodfs', 'wavefish', 'pulsefish', 'eodwaveform',
3548 'wavespectrum', 'pulsephases', 'pulsegaussians', 'pulsespectrum', 'pulsetimes']
3549"""List of all file types generated and supported by the `save_*` and `load_*` functions."""
3552def parse_filename(file_path):
3553 """Parse components of an EOD analysis file name.
3555 Analysis files generated by the `eodanalysis` module are named
3556 according to
3557 ```plain
3558 PATH/RECORDING-CHANNEL-TIME-FTYPE-N.EXT
3559 ```
3561 Parameters
3562 ----------
3563 file_path: string
3564 Path of the file to be parsed.
3566 Returns
3567 -------
3568 recording: string
3569 Path and basename of the recording, i.e. 'PATH/RECORDING'.
3570 A leading './' is removed.
3571 base_path: string
3572 Path and basename of the analysis results,
3573 i.e. 'PATH/RECORDING-CHANNEL-TIME'. A leading './' is removed.
3574 channel: int
3575 Channel of the recording
3576 ('CHANNEL' component of the file name if present).
3577 -1 if not present in `file_path`.
3578 time: float
3579 Start time of analysis window in seconds
3580 ('TIME' component of the file name if present).
3581 `None` if not present in `file_path`.
3582 ftype: string
3583 Type of analysis file (e.g. 'wavespectrum', 'pulsephases', etc.),
3584 ('FTYPE' component of the file name if present).
3585 See `file_types` for a list of all supported file types.
3586 Empty string if not present in `file_path`.
3587 index: int
3588 Index of the EOD.
3589 ('N' component of the file name if present).
3590 -1 if not present in `file_path`.
3591 ext: string
3592 File extension *without* leading period
3593 ('EXT' component of the file name).
3595 """
3596 file_path = Path(file_path)
3597 ext = file_path.suffix
3598 ext = ext[1:]
3599 parts = file_path.stem.split('-')
3600 index = -1
3601 if len(parts) > 0 and parts[-1].isdigit():
3602 index = int(parts[-1])
3603 parts = parts[:-1]
3604 ftype = ''
3605 if len(parts) > 0:
3606 ftype = parts[-1]
3607 parts = parts[:-1]
3608 base_path = file_path.parent / '-'.join(parts)
3609 time = None
3610 if len(parts) > 0 and len(parts[-1]) > 0 and \
3611 parts[-1][0] == 't' and parts[-1][-1] == 's' and \
3612 parts[-1][1:-1].isdigit():
3613 time = float(parts[-1][1:-1])
3614 parts = parts[:-1]
3615 channel = -1
3616 if len(parts) > 0 and len(parts[-1]) > 0 and \
3617 parts[-1][0] == 'c' and parts[-1][1:].isdigit():
3618 channel = int(parts[-1][1:])
3619 parts = parts[:-1]
3620 recording = '-'.join(parts)
3621 return recording, base_path, channel, time, ftype, index, ext
3624def save_analysis(output_basename, zip_file, eod_props, mean_eods, spec_data,
3625 phase_data, pulse_data, wave_eodfs, wave_indices, unit,
3626 verbose, **kwargs):
3627 """Save EOD analysis results to files.
3629 Parameters
3630 ----------
3631 output_basename: string
3632 Path and basename of files to be saved.
3633 zip_file: bool
3634 If `True`, write all analysis results into a zip archive.
3635 eod_props: list of dict
3636 Properties of EODs as returned by `analyze_wave()` and
3637 `analyze_pulse()`.
3638 mean_eods: list of 2D array of floats
3639 Averaged EOD waveforms as returned by `eod_waveform()`,
3640 `analyze_wave()`, and `analyze_pulse()`.
3641 spec_data: list of 2D array of floats
3642 Energy spectra of single pulses as returned by
3643 `analyze_pulse()`.
3644 phase_data: list of 2D array of floats
3645 Properties of phases of pulse EODs as returned by
3646 `analyze_pulse()`.
3647 pulse_data: list of dict
3648 For each pulse fish a dictionary with phase times, amplitudes and standard
3649 deviations of Gaussians fitted to the pulse waveform.
3650 wave_eodfs: list of 2D array of float
3651 Each item is a matrix with the frequencies and powers
3652 (columns) of the fundamental and harmonics (rows) as returned
3653 by `harmonics.harmonic_groups()`.
3654 wave_indices: array of int
3655 Indices identifying each fish in `wave_eodfs` or NaN.
3656 unit: string
3657 Unit of the waveform data.
3658 verbose: int
3659 Verbosity level.
3660 kwargs:
3661 Arguments passed on to `TableData.write()`.
3662 """
3663 def write_file_zip(zf, save_func, output, *args, **kwargs):
3664 if zf is None:
3665 fp = save_func(*args, basename=output, **kwargs)
3666 if verbose > 0 and fp is not None:
3667 print('wrote file', fp)
3668 else:
3669 with io.StringIO() as df:
3670 fp = save_func(*args, basename=df, **kwargs)
3671 if fp is not None:
3672 fp = Path(output + str(fp))
3673 zf.writestr(fp.name, df.getvalue())
3674 if verbose > 0:
3675 print('zipped file', fp.name)
3678 if 'table_format' in kwargs and kwargs['table_format'] == 'py':
3679 with open(output_basename + '.py', 'w') as f:
3680 name = Path(output_basename).stem
3681 for k in range(len((spec_data))):
3682 if len(pulse_data[k]) > 0:
3683 fish = normalize_pulsefish(pulse_data[k])
3684 export_pulsefish(fish, f'{name}-{k}_phases', f)
3685 else:
3686 sdata = spec_data[k]
3687 if len(sdata) > 0 and sdata.shape[1] > 2:
3688 fish = dict(amplitudes=sdata[:, 3], phases=sdata[:, 5])
3689 fish = normalize_wavefish(fish)
3690 export_wavefish(fish, f'{name}-{k}_harmonics', f)
3691 else:
3692 zf = None
3693 if zip_file:
3694 zf = zipfile.ZipFile(output_basename + '.zip', 'w')
3695 # all wave fish in wave_eodfs:
3696 if len(wave_eodfs) > 0:
3697 write_file_zip(zf, save_wave_eodfs, output_basename,
3698 wave_eodfs, wave_indices, **kwargs)
3699 # all wave and pulse fish:
3700 for i, (mean_eod, sdata, pdata, pulse, props) in enumerate(zip(mean_eods, spec_data, phase_data,
3701 pulse_data, eod_props)):
3702 write_file_zip(zf, save_eod_waveform, output_basename,
3703 mean_eod, unit, i, **kwargs)
3704 # spectrum:
3705 if len(sdata)>0:
3706 if sdata.shape[1] == 2:
3707 write_file_zip(zf, save_pulse_spectrum, output_basename,
3708 sdata, unit, i, **kwargs)
3709 else:
3710 write_file_zip(zf, save_wave_spectrum, output_basename,
3711 sdata, unit, i, **kwargs)
3712 # phases:
3713 write_file_zip(zf, save_pulse_phases, output_basename,
3714 pdata, unit, i, **kwargs)
3715 # pulses:
3716 write_file_zip(zf, save_pulse_gaussians, output_basename,
3717 pulse, unit, i, **kwargs)
3718 # times:
3719 write_file_zip(zf, save_pulse_times, output_basename,
3720 props, i, **kwargs)
3721 # wave fish properties:
3722 write_file_zip(zf, save_wave_fish, output_basename,
3723 eod_props, unit, **kwargs)
3724 # pulse fish properties:
3725 write_file_zip(zf, save_pulse_fish, output_basename,
3726 eod_props, unit, **kwargs)
3729def load_analysis(file_pathes):
3730 """Load all EOD analysis files.
3732 Parameters
3733 ----------
3734 file_pathes: list of string
3735 Pathes of the analysis files of a single recording to be loaded.
3737 Returns
3738 -------
3739 mean_eods: list of 2D array of floats
3740 Averaged EOD waveforms: time in seconds, mean, standard deviation, fit.
3741 wave_eodfs: 2D array of floats
3742 EODfs and power of wave type fish.
3743 wave_indices: array of ints
3744 Corresponding indices of fish, can contain negative numbers to
3745 indicate frequencies without fish.
3746 eod_props: list of dict
3747 Properties of EODs. The 'index' property is an index into the
3748 reurned lists.
3749 spec_data: list of 2D array of floats
3750 Amplitude and phase spectrum of wave-type EODs with columns
3751 harmonics, frequency, amplitude, relative amplitude in dB,
3752 relative power in dB, phase, data power in unit squared.
3753 Energy spectrum of single pulse-type EODs with columns
3754 frequency and energy.
3755 phase_data: list of 2D array of floats
3756 Properties of phases of pulse-type EODs with columns
3757 P, time, amplitude, relampl, width, area, relarea
3758 pulse_data: list of dict
3759 For each pulse fish a dictionary with phase times, amplitudes and standard
3760 deviations of Gaussians fitted to the pulse waveform. Use the
3761 functions provided in thunderfish.fakefish to simulate pulse
3762 fish EODs from this data.
3763 recording: string
3764 Path and base name of the recording file.
3765 channel: int
3766 Analysed channel of the recording.
3767 unit: string
3768 Unit of EOD waveform.
3769 """
3770 recording = None
3771 channel = -1
3772 eod_props = []
3773 zf = None
3774 if len(file_pathes) == 1 and Path(file_pathes[0]).suffix[1:] == 'zip':
3775 zf = zipfile.ZipFile(file_pathes[0])
3776 file_pathes = sorted(zf.namelist())
3777 # read wave- and pulse-fish summaries:
3778 pulse_fish = False
3779 wave_fish = False
3780 for f in file_pathes:
3781 recording, _, channel, _, ftype, _, _ = parse_filename(f)
3782 if zf is not None:
3783 f = io.TextIOWrapper(zf.open(f, 'r'))
3784 if ftype == 'wavefish':
3785 eod_props.extend(load_wave_fish(f))
3786 wave_fish = True
3787 elif ftype == 'pulsefish':
3788 eod_props.extend(load_pulse_fish(f))
3789 pulse_fish = True
3790 idx_offs = 0
3791 if wave_fish and not pulse_fish:
3792 idx_offs = sorted([ep['index'] for ep in eod_props])[0]
3793 # load all other files:
3794 neods = len(eod_props)
3795 if neods < 1:
3796 neods = 1
3797 eod_props = [None]
3798 wave_eodfs = np.array([])
3799 wave_indices = np.array([])
3800 mean_eods = [None]*neods
3801 spec_data = [None]*neods
3802 phase_data = [None]*neods
3803 pulse_data = [None]*neods
3804 unit = None
3805 for f in file_pathes:
3806 recording, _, channel, _, ftype, idx, _ = parse_filename(f)
3807 if neods == 1 and idx > 0:
3808 idx = 0
3809 idx -= idx_offs
3810 if zf is not None:
3811 f = io.TextIOWrapper(zf.open(f, 'r'))
3812 if ftype == 'waveeodfs':
3813 wave_eodfs, wave_indices = load_wave_eodfs(f)
3814 elif ftype == 'eodwaveform':
3815 mean_eods[idx], unit = load_eod_waveform(f)
3816 elif ftype == 'wavespectrum':
3817 spec_data[idx], unit = load_wave_spectrum(f)
3818 elif ftype == 'pulsephases':
3819 phase_data[idx], unit = load_pulse_phases(f)
3820 elif ftype == 'pulsegaussians':
3821 pulse_data[idx], unit = load_pulse_gaussians(f)
3822 elif ftype == 'pulsetimes':
3823 pulse_times = load_pulse_times(f)
3824 eod_props[idx]['times'] = pulse_times
3825 eod_props[idx]['peaktimes'] = pulse_times
3826 elif ftype == 'pulsespectrum':
3827 spec_data[idx] = load_pulse_spectrum(f)
3828 # fix wave spectra:
3829 wave_eodfs = [fish.reshape(1, 2) if len(fish)>0 else fish
3830 for fish in wave_eodfs]
3831 if len(wave_eodfs) > 0 and len(spec_data) > 0:
3832 eodfs = []
3833 for idx, fish in zip(wave_indices, wave_eodfs):
3834 if idx >= 0:
3835 spec = spec_data[idx]
3836 specd = np.zeros((np.sum(np.isfinite(spec[:, -1])),
3837 2))
3838 specd[:, 0] = spec[np.isfinite(spec[:, -1]),1]
3839 specd[:, 1] = spec[np.isfinite(spec[:, -1]),-1]
3840 eodfs.append(specd)
3841 else:
3842 specd = np.zeros((10, 2))
3843 specd[:, 0] = np.arange(len(specd))*fish[0,0]
3844 specd[:, 1] = np.nan
3845 eodfs.append(specd)
3846 wave_eodfs = eodfs
3847 return mean_eods, wave_eodfs, wave_indices, eod_props, spec_data, \
3848 phase_data, pulse_data, recording, channel, unit
3851def load_recording(file_path, channel=0, load_kwargs={},
3852 eod_props=None, verbose=0):
3853 """Load recording.
3855 Parameters
3856 ----------
3857 file_path: string or Path
3858 Full path of the file with the recorded data.
3859 Extension is optional. If absent, look for the first file
3860 with a reasonable extension.
3861 channel: int
3862 Channel of the recording to be returned.
3863 load_kwargs: dict
3864 Keyword arguments that are passed on to the
3865 format specific loading functions.
3866 eod_props: list of dict or None
3867 List of EOD properties from which start and end times of
3868 analysis window are extracted.
3869 verbose: int
3870 Verbosity level passed on to load function.
3872 Returns
3873 -------
3874 data: array of float
3875 Data of the requested `channel`.
3876 rate: float
3877 Sampling rate in Hertz.
3878 idx0: int
3879 Start index of the analysis window.
3880 idx1: int
3881 End index of the analysis window.
3882 info_dict: dict
3883 Dictionary with path, basename, species, channel, chanstr, time.
3884 """
3885 data = None
3886 rate = 0.0
3887 idx0 = 0
3888 idx1 = 0
3889 data_file = ''
3890 info_dict = {}
3891 file_path = Path(file_path)
3892 if len(file_path.suffix) > 1:
3893 data_file = file_path
3894 else:
3895 data_files = file_path.parent.glob(file_path.stem + '*')
3896 for dfile in data_files:
3897 if not dfile.suffix[1:] in ['zip'] + list(TableData.ext_formats.values()):
3898 data_file = dfile
3899 break
3900 if data_file.is_file():
3901 all_data = DataLoader(data_file, verbose=verbose, **load_kwargs)
3902 rate = all_data.rate
3903 unit = all_data.unit
3904 ampl_max = all_data.ampl_max
3905 data = all_data[:, channel]
3906 species = get_str(all_data.metadata(), ['species'], default='')
3907 if len(species) > 0:
3908 species += ' '
3909 info_dict = dict(path=all_data.filepath,
3910 basename=all_data.basename(),
3911 species=species,
3912 channel=channel)
3913 if all_data.channels > 1:
3914 if all_data.channels > 100:
3915 info_dict['chanstr'] = f'-c{channel:03d}'
3916 elif all_data.channels > 10:
3917 info_dict['chanstr'] = f'-c{channel:02d}'
3918 else:
3919 info_dict['chanstr'] = f'-c{channel:d}'
3920 else:
3921 info_dict['chanstr'] = ''
3922 idx0 = 0
3923 idx1 = len(data)
3924 if eod_props is not None and len(eod_props) > 0 and 'twin' in eod_props[0]:
3925 idx0 = int(eod_props[0]['twin']*rate)
3926 if len(eod_props) > 0 and 'window' in eod_props[0]:
3927 idx1 = idx0 + int(eod_props[0]['window']*rate)
3928 info_dict['time'] = '-t{idx0/rate:.0f}s'
3929 all_data.close()
3931 return data, rate, idx0, idx1, info_dict
3934def add_eod_analysis_config(cfg, thresh_fac=0.8, percentile=0.1,
3935 win_fac=2.0, min_win=0.01, max_eods=None,
3936 min_sem=False, unfilter_cutoff=0.0,
3937 flip_wave='none', flip_pulse='none',
3938 n_harm=10, min_pulse_win=0.001,
3939 peak_thresh_fac=0.01, min_dist=50.0e-6,
3940 width_frac = 0.5, fit_frac = 0.5,
3941 freq_resolution=1.0, fade_frac=0.0,
3942 ipi_cv_thresh=0.5, ipi_percentile=30.0):
3943 """Add all parameters needed for the eod analysis functions as a new
3944 section to a configuration.
3946 Parameters
3947 ----------
3948 cfg: ConfigFile
3949 The configuration.
3951 See `eod_waveform()`, `analyze_wave()`, and `analyze_pulse()` for
3952 details on the remaining arguments.
3953 """
3954 cfg.add_section('EOD analysis:')
3955 cfg.add('eodSnippetFac', win_fac, '', 'The duration of EOD snippets is the EOD period times this factor.')
3956 cfg.add('eodMinSnippet', min_win, 's', 'Minimum duration of cut out EOD snippets.')
3957 cfg.add('eodMaxEODs', max_eods or 0, '', 'The maximum number of EODs used to compute the average EOD. If 0 use all EODs.')
3958 cfg.add('eodMinSem', min_sem, '', 'Use minimum of s.e.m. to set maximum number of EODs used to compute the average EOD.')
3959 cfg.add('unfilterCutoff', unfilter_cutoff, 'Hz', 'If non-zero remove effect of high-pass filter with this cut-off frequency.')
3960 cfg.add('flipWaveEOD', flip_wave, '', 'Flip EOD of wave fish to make largest extremum positive (flip, none, or auto).')
3961 cfg.add('flipPulseEOD', flip_pulse, '', 'Flip EOD of pulse fish to make the first large peak positive (flip, none, or auto).')
3962 cfg.add('eodHarmonics', n_harm, '', 'Number of harmonics fitted to the EOD waveform.')
3963 cfg.add('eodMinPulseSnippet', min_pulse_win, 's', 'Minimum duration of cut out EOD snippets for a pulse fish.')
3964 cfg.add('eodPeakThresholdFactor', peak_thresh_fac, '', 'Threshold for detection of peaks and troughs in pulse EODs as a fraction of the pulse amplitude.')
3965 cfg.add('eodMinimumDistance', min_dist, 's', 'Minimum distance between peaks and troughs in a EOD pulse.')
3966 cfg.add('eodPulseWidthFraction', width_frac, '', 'The width of a pulse is measured at this fraction of the pulse height.')
3967 cfg.add('eodExponentialFitFraction', fit_frac, '', 'An exponential function is fitted on the tail of a pulse starting at this fraction of the height of the last peak.')
3968 cfg.add('eodPulseFrequencyResolution', freq_resolution, 'Hz', 'Frequency resolution of single pulse spectrum.')
3969 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.')
3970 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.')
3971 cfg.add('ipiPercentile', ipi_percentile, '%', 'Use only interpulse intervals shorter than this percentile to compute EOD frequency.')
3974def eod_waveform_args(cfg):
3975 """Translates a configuration to the respective parameter names of
3976 the function `eod_waveform()`.
3978 The return value can then be passed as key-word arguments to this
3979 function.
3981 Parameters
3982 ----------
3983 cfg: ConfigFile
3984 The configuration.
3986 Returns
3987 -------
3988 a: dict
3989 Dictionary with names of arguments of the `eod_waveform()` function
3990 and their values as supplied by `cfg`.
3991 """
3992 a = cfg.map({'win_fac': 'eodSnippetFac',
3993 'min_win': 'eodMinSnippet',
3994 'max_eods': 'eodMaxEODs',
3995 'min_sem': 'eodMinSem',
3996 'unfilter_cutoff': 'unfilterCutoff'})
3997 return a
4000def analyze_wave_args(cfg):
4001 """Translates a configuration to the respective parameter names of
4002 the function `analyze_wave()`.
4004 The return value can then be passed as key-word arguments to this
4005 function.
4007 Parameters
4008 ----------
4009 cfg: ConfigFile
4010 The configuration.
4012 Returns
4013 -------
4014 a: dict
4015 Dictionary with names of arguments of the `analyze_wave()` function
4016 and their values as supplied by `cfg`.
4017 """
4018 a = cfg.map({'n_harm': 'eodHarmonics',
4019 'power_n_harmonics': 'powerNHarmonics',
4020 'flip_wave': 'flipWaveEOD'})
4021 return a
4024def analyze_pulse_args(cfg):
4025 """Translates a configuration to the respective parameter names of
4026 the function `analyze_pulse()`.
4028 The return value can then be passed as key-word arguments to this
4029 function.
4031 Parameters
4032 ----------
4033 cfg: ConfigFile
4034 The configuration.
4036 Returns
4037 -------
4038 a: dict
4039 Dictionary with names of arguments of the `analyze_pulse()` function
4040 and their values as supplied by `cfg`.
4041 """
4042 a = cfg.map({'min_pulse_win': 'eodMinPulseSnippet',
4043 'peak_thresh_fac': 'eodPeakThresholdFactor',
4044 'min_dist': 'eodMinimumDistance',
4045 'width_frac': 'eodPulseWidthFraction',
4046 'fit_frac': 'eodExponentialFitFraction',
4047 'flip_pulse': 'flipPulseEOD',
4048 'freq_resolution': 'eodPulseFrequencyResolution',
4049 'fade_frac': 'eodPulseFadeFraction',
4050 'ipi_cv_thresh': 'ipiCVThresh',
4051 'ipi_percentile': 'ipiPercentile'})
4052 a['fade_frac'] *= 0.01
4053 return a
4056def add_species_config(cfg, species_file='none', wave_max_rms=0.2,
4057 pulse_max_rms=0.2):
4058 """Add parameters needed for assigning EOD waveforms to species.
4060 Parameters
4061 ----------
4062 cfg: ConfigFile
4063 The configuration.
4064 species_file: string
4065 File path to a file containing species names and corresponding
4066 file names of EOD waveform templates. If 'none', no species
4067 assignemnt is performed.
4068 wave_max_rms: float
4069 Maximum allowed rms difference (relative to standard deviation
4070 of EOD waveform) to an EOD waveform template for assignment to
4071 a wave fish species.
4072 pulse_max_rms: float
4073 Maximum allowed rms difference (relative to standard deviation
4074 of EOD waveform) to an EOD waveform template for assignment to
4075 a pulse fish species.
4076 """
4077 cfg.add_section('Species assignment:')
4078 cfg.add('speciesFile', species_file, '', 'File path to a file containing species names and corresponding file names of EOD waveform templates.')
4079 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.')
4080 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.')
4083def add_eod_quality_config(cfg, max_clipped_frac=0.1, max_variance=0.0,
4084 max_rms_error=0.05, min_power=-100.0, max_thd=0.0,
4085 max_crossings=4, max_relampl_harm1=0.0,
4086 max_relampl_harm2=0.0, max_relampl_harm3=0.0):
4087 """Add parameters needed for assesing the quality of an EOD waveform.
4089 Parameters
4090 ----------
4091 cfg: ConfigFile
4092 The configuration.
4094 See `wave_quality()` and `pulse_quality()` for details on
4095 the remaining arguments.
4096 """
4097 cfg.add_section('Waveform selection:')
4098 cfg.add('maximumClippedFraction', max_clipped_frac, '', 'Take waveform of the fish with the highest power only if the fraction of clipped signals is below this value.')
4099 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.')
4100 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.')
4101 cfg.add('minimumPower', min_power, 'dB', 'Skip waveform of wave fish if its power is smaller than this value.')
4102 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.')
4103 cfg.add('maximumCrossings', max_crossings, '', 'Maximum number of zero crossings per EOD period.')
4104 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.')
4105 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.')
4106 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.')
4109def wave_quality_args(cfg):
4110 """Translates a configuration to the respective parameter names of
4111 the function `wave_quality()`.
4113 The return value can then be passed as key-word arguments to this
4114 function.
4116 Parameters
4117 ----------
4118 cfg: ConfigFile
4119 The configuration.
4121 Returns
4122 -------
4123 a: dict
4124 Dictionary with names of arguments of the `wave_quality()` function
4125 and their values as supplied by `cfg`.
4126 """
4127 a = cfg.map({'max_clipped_frac': 'maximumClippedFraction',
4128 'max_rms_sem': 'maximumVariance',
4129 'max_rms_error': 'maximumRMSError',
4130 'min_power': 'minimumPower',
4131 'max_crossings': 'maximumCrossings',
4132 'min_freq': 'minimumFrequency',
4133 'max_freq': 'maximumFrequency',
4134 'max_thd': 'maximumTotalHarmonicDistortion',
4135 'max_db_diff': 'maximumPowerDifference',
4136 'max_harmonics_db': 'maximumHarmonicsPower',
4137 'max_relampl_harm1': 'maximumFirstHarmonicAmplitude',
4138 'max_relampl_harm2': 'maximumSecondHarmonicAmplitude',
4139 'max_relampl_harm3': 'maximumThirdHarmonicAmplitude'})
4140 return a
4143def pulse_quality_args(cfg):
4144 """Translates a configuration to the respective parameter names of
4145 the function `pulse_quality()`.
4147 The return value can then be passed as key-word arguments to this
4148 function.
4150 Parameters
4151 ----------
4152 cfg: ConfigFile
4153 The configuration.
4155 Returns
4156 -------
4157 a: dict
4158 Dictionary with names of arguments of the `pulse_quality()` function
4159 and their values as supplied by `cfg`.
4160 """
4161 a = cfg.map({'max_clipped_frac': 'maximumClippedFraction',
4162 'max_rms_sem': 'maximumRMSNoise'})
4163 return a
4166def main():
4167 import matplotlib.pyplot as plt
4168 from .fakefish import pulsefish_eods
4170 print('Analysis of EOD waveforms.')
4172 # data:
4173 rate = 96_000
4174 data = pulsefish_eods('Triphasic', 83.0, rate, 5.0, noise_std=0.02)
4175 unit = 'mV'
4176 eod_idx, _ = detect_peaks(data, 1.0)
4177 eod_times = eod_idx/rate
4179 # analyse EOD:
4180 mean_eod, eod_times = eod_waveform(data, rate, eod_times)
4181 mean_eod, props, peaks, pulses, energy = analyze_pulse(mean_eod, eod_times)
4183 # plot:
4184 fig, axs = plt.subplots(1, 2)
4185 plot_eod_waveform(axs[0], mean_eod, props, peaks, unit=unit)
4186 axs[0].set_title(f'{props["type"]} fish: EODf = {props["EODf"]:.1f} Hz')
4187 plot_pulse_spectrum(axs[1], energy, props)
4188 plt.show()
4191if __name__ == '__main__':
4192 main()