Coverage for src/thunderfish/pulseanalysis.py: 43%
1311 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 20:05 +0000
1"""
2Analysis of pulse-type EOD waveforms.
4## Analysis of pulse-type EODs
6### Analysis
8- `condition_pulse()`: subtract offset, flip, shift, and cut out pulse EOD waveform.
9- `analyze_pulse_properties()`: characterize basic properties of a pulse-type EOD.
10- `analyze_pulse_phases()`: characterize all phases of a pulse-type EOD.
11- `decompose_pulse()`: decompose single pulse waveform into sum of Gaussians.
12- `analyze_pulse_tail()`: fit exponential to last peak/trough of pulse EOD.
13- `pulse_spectrum()`: spectrum of a single pulse-type EOD.
14- `pulsetrain(): synthesize train of pulse-type EODs.
15- `pulsetrain_spectrum()`: power spectrum of train of pulse-type EODs.
16- `analyze_pulse_spectrum()`: analyze the spectrum of a pulse-type EOD.
17- `analyze_pulse_intervals()`: basic statistics of interpulse intervals.
19### Complete analysis
21Calls all the functions listed above:
23- `analyze_pulse()`: analyze the EOD waveform of a pulse fish.
25## Quality assessment
27- `pulse_quality()`: asses quality of EOD waveform of a pulse fish.
29## Visualization
31- `plot_pulse_eodtimes()`: mark pulse EODs in a plot of an EOD recording.
32- `plot_pulse_rate()`: plot pulse rates.
33- `plot_pulse_eod()`: plot and annotate a pulse-type EOD waveform.
34- `plot_pulse_spectrum()`: plot and annotate spectrum of single pulse EOD.
35- `plot_pulse_distribution()`: plot EOD waveform, amplitude distribution and cumulative.
37## Storage
39- `save_pulse_fish()`: save properties of pulse EODs to file.
40- `load_pulse_fish()`: load properties of pulse EODs from file.
41- `save_pulse_spectrum()`: save spectrum of pulse EOD to file.
42- `load_pulse_spectrum()`: load spectrum of pulse EOD from file.
43- `save_pulse_phases()`: save phase properties of pulse EOD to file.
44- `load_pulse_phases()`: load phase properties of pulse EOD from file.
45- `save_pulse_gaussians()`: save Gaussian phase properties of pulse EOD to file.
46- `load_pulse_gaussians()`: load Gaussian phase properties of pulse EOD from file.
47- `save_pulse_times()`: save times of pulse EOD to file.
48- `load_pulse_times()`: load times of pulse EOD from file.
50## Fit functions
52- `gaussian_sum()`: sum of Gaussians.
53- `gaussian_sum_spectrum()`: energy spectrum of sum of Gaussians.
54- `gaussian_sum_costs()`: cost function for fitting sum of Gaussians.
55- `exp_decay()`: exponential decay.
57## Configuration
59- `add_analyze_pulse_config()`: add parameters for `analyze_pulse()` to configuration.
60- `analyze_pulse_args()`: retrieve parameters for `analyze_pulse()` from configuration.
61- `add_pulse_quality_config()`: add parameters for `pulse_quality()` to configuration.
62- `pulse_quality_args()`: retrieve parameters for `pulse_quality()` from configuration.
64"""
66import numpy as np
68try:
69 import matplotlib.pyplot as plt
70 from matplotlib.ticker import MultipleLocator
71except ImportError:
72 pass
74from pathlib import Path
75from scipy.optimize import curve_fit, minimize
76from scipy.signal import savgol_filter, find_peaks
77from thunderlab.eventdetection import detect_peaks, peak_width
78from thunderlab.powerspectrum import next_power_of_two, nfft, psd, decibel
79from thunderlab.tabledata import TableData
81from .files import open_from_zip, close_zip
82from .fakefish import pulsefish_waveform, pulsefish_spectrum
85def condition_pulse(eod, ratetime=None, sem=None, flip='none',
86 baseline_mode='mode', baseline_frac=0.05,
87 large_phase_frac=0.2, min_pulse_win=0.001):
88 """Subtract offset, flip, shift, and cut out pulse EOD waveform.
90 Parameters
91 ----------
92 eod: 1-D or 2-D array of float
93 The eod waveform to be conditioned.
94 If an 1-D array, then this is the waveform and you
95 need to also pass a sampling rate in `rate`.
96 If a 2-D array, then first column is time in seconds and second
97 column is the eod waveform. If a third column is present,
98 it is interpreted as the standard error of the mean
99 corresponding to the averaged waveform.
100 Further columns are ignored.
101 ratetime: None or float or array of float
102 If a 1-D array is passed on to `eod` then either the sampling
103 rate in Hertz or the time array corresponding to `eod`.
104 sem: None or float or array of float
105 If a 1-D array is passed on to `eod`, then the optional
106 standard error of the averaged waveform. Either as a single value
107 or as an 1-D array for the whole waveform.
108 flip: 'auto', 'none', 'flip'
109 - 'auto' flip waveform such that the first large extremum is positive.
110 - 'flip' flip waveform.
111 - 'none' do not flip waveform.
112 baseline_mode: "left", "right", "mean", "mode"
113 Method to compute baseline voltage:
114 "left": mean of leftmost `baseline_frac` data points.
115 "right": mean of rightmost `baseline_frac` data points.
116 "mean": mean of both leftmost and rightmost `baseline_frac` data points.
117 "mode": mode of the distribution of the data values,
118 i.e. most common value.
119 baseline_frac: float
120 Fraction of data points from which the amplitude offset is computed.
121 large_phase_frac: float
122 Minimum amplitude of a large phase as a fraction of the largest one.
123 min_pulse_win: float
124 The minimum size of the cut-out EOD waveform.
126 Returns
127 -------
128 eod: 1-D or 2-D array of float
129 The input `eod` in the format of the input `eod`
130 with the conditioned waveform and optional time.
131 time: 1-D array of float
132 In case `eod` is a 1-D array, then the time array is also returned.
133 toffs: float
134 Time that was subtracted from the time axis,
135 such that the maximum of the EOD waveform was shifted to zero.
136 aoffs: float
137 Offset that was subtracted from the EOD waveform.
138 This is the average over the first `baseline_frac` data points
139 of the EOD waveform.
140 flipped: bool
141 True if waveform was flipped.
142 noise_thresh: float
143 Minimum threshold that is just higher than the noise level
144 within the first `baseline_frac` data points of the EOD waveform.
146 """
147 def mode(data):
148 n = max(10, len(data)//100)
149 order = 1
150 sdata = np.sort(data)
151 dsd = sdata[-1] - sdata[0]
152 mask = (sdata > sdata[0] + 0.01*dsd) & (sdata < sdata[0] + 0.99*dsd)
153 filtered = savgol_filter(sdata, n, order, 1)
154 return sdata[mask][np.argmin(filtered[mask])]
156 if eod.ndim == 2:
157 time = eod[:, 0]
158 meod = eod[:, 1]
159 if eod.shape[1] > 2:
160 sem = eod[:, 2]
161 else:
162 meod = eod
163 if isinstance(ratetime, (list, tuple, np.ndarray)):
164 time = ratetime
165 else:
166 time = np.arange(len(meod))/ratetime
167 if np.isscalar(sem):
168 sem = np.ones(len(meod))*sem
169 dt = time[1] - time[0]
171 # check for additional high peaks:
172 l_idx = 0
173 r_idx = len(meod)
174 data = meod - mode(meod)
175 max_ampl = np.max(data)
176 min_ampl = np.min(data)
177 if max_ampl >= min_ampl:
178 peak_indices, _ = find_peaks(data, height=0.7*max_ampl,
179 distance=int(np.ceil(0.001/dt)))
180 else:
181 peak_indices, _ = find_peaks(-data, height=0.7*min_ampl,
182 distance=int(np.ceil(0.001/dt)))
183 main_idx = np.argmin(np.abs(peak_indices - len(data)//3))
184 if main_idx > 0:
185 l_idx = peak_indices[main_idx]
186 l_idx -= (peak_indices[main_idx] - peak_indices[main_idx - 1])//3
187 if main_idx < len(peak_indices) - 1:
188 r_idx = peak_indices[main_idx]
189 r_idx += 2*(peak_indices[main_idx + 1] - peak_indices[main_idx])//3
190 # cut out relevant signal:
191 if time[r_idx - 1] - time[l_idx] < min_pulse_win:
192 tmax = time[main_idx]
193 mask = (time >= tmax - min_pulse_win/3) & (time <= tmax + 2*min_pulse_win/3)
194 meod = meod[mask]
195 time = time[mask]
196 if eod.ndim == 2:
197 eod = eod[mask, :]
198 else:
199 meod = meod[l_idx:r_idx]
200 time = time[l_idx:r_idx]
201 if eod.ndim == 2:
202 eod = eod[l_idx:r_idx, :]
204 # plot_pulse_distribution(eod)
206 n_base = max(5, int(baseline_frac*len(meod)))
207 # subtract baseline:
208 if baseline_mode.lower() == 'left':
209 aoffs = np.mean(meod[:n_base])
210 elif baseline_mode.lower() == 'right':
211 aoffs = np.mean(meod[-n_base:])
212 elif baseline_mode.lower() == 'mean':
213 aoffs = 0.5*(np.mean(meod[:n_base]) + np.mean(meod[-n_base:]))
214 elif baseline_mode.lower() == 'mode':
215 aoffs = mode(meod)
216 else:
217 raise ValueError(f'ERROR in condition_pulse(): invalid baseline_mode "{baseline_mode}"! must be one of "left", "right", "mean" or "mode".')
218 meod -= aoffs
220 # flip waveform:
221 max_idx = np.argmax(meod)
222 min_idx = np.argmin(meod)
223 flipped = False
224 if flip.lower() in ['flip', 'true', 'yes']:
225 flipped = True
226 elif flip.lower() == 'auto':
227 max_ampl = abs(meod[max_idx])
228 min_ampl = abs(meod[min_idx])
229 amplitude = max(max_ampl, min_ampl)
230 if max_ampl > large_phase_frac*amplitude and \
231 min_ampl > large_phase_frac*amplitude:
232 # two major peaks:
233 if min_idx < max_idx:
234 flipped = True
235 elif min_ampl > large_phase_frac*amplitude:
236 flipped = True
237 if flipped:
238 meod = -meod
239 idx = min_idx
240 min_idx = max_idx
241 max_idx = idx
242 max_ampl = abs(meod[max_idx])
243 min_ampl = abs(meod[min_idx])
245 # shift maximum to zero:
246 toffs = time[max_idx]
247 time -= toffs
249 # threshold from baseline maximum and minimum:
250 th_max = np.max(meod[:n_base])
251 th_min = np.min(meod[:n_base])
252 range_thresh = 2*(th_max - th_min)
253 # threshold from standard error:
254 if sem is not None:
255 msem = np.mean(sem[:n_base])
256 sem_thresh = 8*msem
257 # noise threshold:
258 noise_thresh = max(range_thresh, sem_thresh)
259 if noise_thresh > 0.5*max_ampl:
260 noise_thresh = 0.5*max_ampl
262 # return offset, flipped, shifted, and cut out waveform:
263 if eod.ndim == 2:
264 eod[:, 0] = time
265 eod[:, 1] = meod
266 return eod, toffs, aoffs, flipped, noise_thresh
267 else:
268 return meod, time, toffs, aoffs, flipped, noise_thresh
271def analyze_pulse_properties(noise_thresh, eod, ratetime=None):
272 """Characterize basic properties of a pulse-type EOD.
274 Parameters
275 ----------
276 noise_thresh: float
277 Minimum threshold that is just higher than the baseline noise level.
278 As returned by `condition_pulse()`.
279 eod: 1-D or 2-D array
280 The eod waveform to be analyzed.
281 If an 1-D array, then this is the waveform and you
282 need to also pass a sampling rate in `rate`.
283 If a 2-D array, then first column is time in seconds and second
284 column is the eod waveform. Further columns are ignored.
285 ratetime: None or float or array of float
286 If a 1-D array is passed on to `eod` then either the sampling
287 rate in Hertz or the time array corresponding to `eod`.
289 Returns
290 -------
291 pp_ampl: float
292 Peak-to-peak amplitude.
293 rel_max_ampl: float
294 Amplitude of minimum or maximum, whichever is larger, relative to p-p amplitude.
295 dist: float
296 Temporal distance between largest negative trough and positive peak.
297 total_area: float
298 Integral under absolute values of EOD waveform, i.e. sum of areas
299 under positive and negative phases.
300 polarity_balance: float
301 Contrast between positive and negative areas of EOD waveform, i.e.
302 (pos_area - neg_area)/(pos_area + neg_area).
303 center: float
304 Center of mass (first moment when treating the absolute value of
305 the waveform as a distribution).
306 stdev: float
307 Standard deviation of mass (square root of second central moment
308 when treating the absolute value of the waveform as a distribution).
309 median: float
310 Median of the distribution of the absolute EOD waveform.
311 quartile1: float
312 First quartile of the distribution of the absolute EOD waveform.
313 quartile3: float
314 Third quartile of the distribution of the absolute EOD waveform.
315 """
316 if eod.ndim == 2:
317 time = eod[:, 0]
318 eod = eod[:, 1]
319 elif isinstance(ratetime, (list, tuple, np.ndarray)):
320 time = ratetime
321 else:
322 time = np.arange(len(eod))/rate
323 dt = time[1] - time[0]
325 # amplitudes:
326 pos_idx = np.argmax(eod)
327 pos_ampl = abs(eod[pos_idx])
328 neg_idx = np.argmin(eod)
329 neg_ampl = abs(eod[neg_idx])
330 if neg_ampl < noise_thresh:
331 neg_ampl = 0
332 pp_ampl = pos_ampl + neg_ampl
333 max_ampl = max(pos_ampl, neg_ampl)
334 rel_max_ampl = max_ampl/pp_ampl if pp_ampl > 0 else 1
336 # timing:
337 dist = time[neg_idx] - time[pos_idx]
338 if pos_ampl < noise_thresh or neg_ampl < noise_thresh:
339 dist = np.inf
341 # integrals (areas) and polarity balance:
342 pos_area = abs(np.sum(eod[eod >= 0]))*dt
343 neg_area = abs(np.sum(eod[eod < 0]))*dt
344 total_area = np.sum(np.abs(eod))*dt
345 integral_area = np.sum(eod)*dt # = pos_area - neg_area
346 polarity_balance = integral_area/total_area
348 # cumulative:
349 cumul = np.cumsum(np.abs(eod))/np.sum(np.abs(eod))
350 median = time[np.argmax(cumul > 0.5)]
351 quartile1 = time[np.argmax(cumul > 0.25)]
352 quartile3 = time[np.argmax(cumul > 0.75)]
354 return pp_ampl, rel_max_ampl, dist, \
355 total_area, polarity_balance, \
356 median, quartile1, quartile3
359def analyze_pulse_phases(peak_thresh, startend_thresh,
360 eod, ratetime=None,
361 min_dist=50.0e-6, width_frac=0.5):
362 """Characterize all phases of a pulse-type EOD.
364 Parameters
365 ----------
366 peak_thresh: float
367 Threshold for detecting peaks and troughs.
368 startend_thresh: float or None
369 Threshold for detecting start and end time of EOD.
370 If None, use `peak_thresh`.
371 eod: 1-D or 2-D array
372 The eod waveform to be analyzed.
373 If an 1-D array, then this is the waveform and you
374 need to also pass a sampling rate in `rate`.
375 If a 2-D array, then first column is time in seconds and second
376 column is the eod waveform. Further columns are ignored.
377 ratetime: None or float or array of float
378 If a 1-D array is passed on to `eod` then either the sampling
379 rate in Hertz or the time array corresponding to `eod`.
380 min_dist: float
381 Minimum distance between peak and troughs of the pulse.
382 width_frac: float
383 The width of an EOD phase is measured at this fraction of a peak's
384 height (0-1).
386 Returns
387 -------
388 tstart: float
389 Start time of EOD waveform, i.e. the first time it crosses `threshold`.
390 tend: float
391 End time of EOD waveform, i.e. the last time it falls under `threshold`.
392 phases: dict
393 Dictionary with
395 - "indices": indices of each phase
396 (1 is P1, i.e. the largest positive peak)
397 - "times": times of each phase relative to P1 in seconds
398 - "amplitudes": amplitudes of each phase
399 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
400 - "widths": widths of each phase at `width_frac` height
401 - "areas": areas of each phase
402 - "relareas": areas of the phases relative to total area
403 - "zeros": time of zero crossing towards next phase in seconds
405 Empty dictionary if waveform is not a pulse EOD.
407 """
408 if eod.ndim == 2:
409 time = eod[:, 0]
410 eod = eod[:, 1]
411 elif isinstance(ratetime, (list, tuple, np.ndarray)):
412 time = ratetime
413 else:
414 time = np.arange(len(eod))/rate
415 dt = time[1] - time[0]
417 # start and end time:
418 if startend_thresh is None:
419 startend_thresh = peak_thresh
420 l_idx = np.argmax(np.abs(eod) > startend_thresh)
421 r_idx = len(eod) - 1 - np.argmax(np.abs(eod[::-1]) > startend_thresh)
422 tstart = time[l_idx]
423 tend = time[r_idx]
425 # find peaks and troughs:
426 peak_idx, trough_idx = detect_peaks(eod, peak_thresh)
427 if len(peak_idx) == 0:
428 return tstart, tend, {}
429 # and their width:
430 peak_widths = peak_width(time, eod, peak_idx, trough_idx,
431 peak_frac=width_frac, base='max')
432 trough_widths = peak_width(time, -eod, trough_idx, peak_idx,
433 peak_frac=width_frac, base='max')
434 # combine peaks and troughs:
435 pt_idx = np.concatenate((peak_idx, trough_idx))
436 pt_widths = np.concatenate((peak_widths, trough_widths))
437 pts_idx = np.argsort(pt_idx)
438 phase_list = pt_idx[pts_idx]
439 width_list = pt_widths[pts_idx]
440 # remove phases at the start and end of the signal:
441 n = len(eod)//20
442 if n < 5:
443 n = 5
444 mask = (phase_list > n) & (phase_list < len(eod) - 20)
445 phase_list = phase_list[mask]
446 width_list = width_list[mask]
447 if len(phase_list) == 0:
448 return tstart, tend, {}
449 # remove multiple peaks that are too close:
450 # TODO: XXX replace by Dexters function that keeps the maximum peak
451 rmidx = [(k, k+1) for k in np.where(np.diff(time[phase_list]) < min_dist)[0]]
452 # flatten and keep maximum and minimum phase:
453 max_idx = np.argmax(eod)
454 min_idx = np.argmin(eod)
455 rmidx = np.unique([k for kk in rmidx for k in kk
456 if phase_list[k] != max_idx and phase_list[k] != min_idx])
457 # delete:
458 if len(rmidx) > 0:
459 phase_list = np.delete(phase_list, rmidx)
460 width_list = np.delete(width_list, rmidx)
461 if len(phase_list) == 0:
462 return tstart, tend, {}
463 # find P1:
464 p1i = np.argmax(phase_list == max_idx)
465 # amplitudes:
466 amplitudes = eod[phase_list]
467 max_ampl = np.abs(amplitudes[p1i])
468 # phase areas and zeros:
469 phase_areas = np.zeros(len(phase_list))
470 zero_times = np.zeros(len(phase_list))
471 for i in range(len(phase_list)):
472 sign_fac = np.sign(eod[phase_list[i]])
473 i0 = phase_list[i - 1] if i > 0 else 0
474 i1 = phase_list[i + 1] if i + 1 < len(phase_list) else len(eod)
475 if i0 > 0 and sign_fac*eod[i0] > 0 and \
476 i1 < len(eod) and sign_fac*eod[i1] > 0:
477 phase_areas[i] = 0
478 else:
479 snippet = sign_fac*eod[i0:i1]
480 phase_areas[i] = np.sum(snippet[snippet > 0])*dt
481 phase_areas[i] *= sign_fac
482 if i < len(phase_list) - 1:
483 i0 = phase_list[i]
484 snippet = eod[i0:i1]
485 stimes = time[i0:i1]
486 zidx = np.nonzero(snippet[:-1]*snippet[1:] < 0)[0]
487 if len(zidx) == 0:
488 zero_times[i] = np.nan
489 else:
490 zidx = zidx[len(zidx)//2] # reduce to single zero crossing
491 snippet = snippet[zidx:zidx + 2]
492 stimes = stimes[zidx:zidx + 2]
493 if sign_fac > 0:
494 zero_times[i] = np.interp(0, snippet[::-1], stimes[::-1])
495 else:
496 zero_times[i] = np.interp(0, snippet, stimes)
497 else:
498 zero_times[i] = np.nan
499 total_area = np.sum(np.abs(phase_areas))
500 # store phase properties:
501 phases = dict(indices=np.arange(len(phase_list)) + 1 - p1i,
502 times=time[phase_list],
503 amplitudes=amplitudes,
504 relamplitudes=amplitudes/max_ampl,
505 widths=width_list,
506 areas=phase_areas,
507 relareas=phase_areas/total_area,
508 zeros=zero_times)
509 return tstart, tend, phases
512def gaussian_sum(x, *pas):
513 """Sum of Gaussians.
515 Parameters
516 ----------
517 x: array of float
518 The x array over which the sum of Gaussians is evaluated.
519 *pas: list of floats
520 The parameters of the Gaussians in a flat list.
521 Position, amplitude, and standard deviation of first Gaussian,
522 position, amplitude, and standard deviation of second Gaussian,
523 and so on.
525 Returns
526 -------
527 sg: array of float
528 The sum of Gaussians for the times given in `t`.
529 """
530 sg = np.zeros(len(x))
531 for pos, ampl, std in zip(pas[0:-2:3], pas[1:-1:3], pas[2::3]):
532 sg += ampl*np.exp(-0.5*((x - pos)/std)**2)
533 return sg
536def gaussian_sum_spectrum(f, *pas):
537 """Energy spectrum of sum of Gaussians.
539 Parameters
540 ----------
541 f: 1-D array of float
542 The frequencies at which to evaluate the spectrum.
543 *pas: list of floats
544 The parameters of the Gaussians in a flat list.
545 Position, amplitude, and standard deviation of first Gaussian,
546 position, amplitude, and standard deviation of second Gaussian,
547 and so on.
549 Returns
550 -------
551 energy: 1-D array of float
552 The one-sided energy spectrum of the sum of Gaussians.
553 """
554 spec = np.zeros(len(f), dtype=complex)
555 for dt, a, s in zip(pas[0:-2:3], pas[1:-1:3], pas[2::3]):
556 gauss = a*np.sqrt(2*np.pi)*s*np.exp(-0.5*(2*np.pi*s*f)**2)
557 shift = np.exp(-2j*np.pi*f*dt)
558 spec += gauss*shift
559 spec *= np.sqrt(2) # because of one-sided spectrum
560 return np.abs(spec)**2
563def gaussian_sum_costs(pas, time, eod, freqs, energy):
564 """ Cost function for fitting sum of Gaussian to pulse EOD.
566 Parameters
567 ----------
568 pas: list of floats
569 The pulse parameters in a flat list.
570 Position, amplitude, and standard deviation of first phase,
571 position, amplitude, and standard deviation of second phase,
572 and so on.
573 time: 1-D array of float
574 Time points of the EOD waveform.
575 eod: 1-D array of float
576 The real EOD waveform.
577 freqs: 1-D array of float
578 The frequency components of the spectrum.
579 energy: 1-D array of float
580 The energy spectrum of the real pulse.
582 Returns
583 -------
584 costs: float
585 Weighted sum of rms waveform difference and rms spectral difference.
586 """
587 eod_fit = gaussian_sum(time, *pas)
588 eod_rms = np.sqrt(np.mean((eod_fit - eod)**2))/np.max(np.abs(eod))
589 level = decibel(energy)
590 level_range = 30
591 n = np.argmax(level)
592 energy_fit = gaussian_sum_spectrum(freqs, *pas)
593 level_fit = decibel(energy_fit)
594 weight = np.ones(n)
595 weight[freqs[:n] < 10] = 100
596 weight /= np.sum(weight)
597 spec_rms = np.sqrt(np.mean(weight*(level_fit[:n] - level[:n])**2))/level_range
598 costs = eod_rms + 5*spec_rms
599 #print(f'{costs:.4f}, {eod_rms:.4f}, {spec_rms:.4f}')
600 return costs
603def decompose_pulse(eod, freqs, energy, phases, width_frac=0.5, verbose=0):
604 """Decompose single pulse waveform into sum of Gaussians.
606 Use the output to simulate pulse-type EODs using the functions
607 provided in the thunderfish.fakefish module.
609 Parameters
610 ----------
611 eod: 2-D array of float
612 The eod waveform. First column is time in seconds and second
613 column is the eod waveform. Further columns are ignored.
614 freqs: 1-D array of float
615 The frequency components of the spectrum.
616 energy: 1-D array of float
617 The energy spectrum of the real pulse.
618 phases: dict
619 Properties of the EOD phases as returned by analyze_pulse_phases().
620 width_frac: float
621 The width of a peak is measured at this fraction of a peak's
622 height (0-1).
623 verbose: int
624 Verbosity level passed for error and info messages.
626 Returns
627 -------
628 pulse: dict
629 Dictionary with
631 - "times": phase times in seconds,
632 - "amplitudes": amplitudes, and
633 - "stdevs": standard deviations in seconds
635 of Gaussians fitted to the pulse waveform. Use the functions
636 provided in thunderfish.fakefish to simulate pulse fish EODs
637 from this data.
639 """
640 pulse = {}
641 if len(phases) == 0:
642 return pulse
643 if eod.ndim == 2:
644 time = eod[:, 0]
645 eod = eod[:, 1]
646 elif isinstance(ratetime, (list, tuple, np.ndarray)):
647 time = ratetime
648 else:
649 time = np.arange(len(eod))/rate
650 # convert half width to standard deviation:
651 fac = 0.5/np.sqrt(-2*np.log(width_frac))
652 # fit parameter as single list:
653 tas = []
654 for t, a, s in zip(phases['times'], phases['amplitudes'],
655 phases['widths']*fac):
656 tas.extend((t, a, s))
657 tas = np.asarray(tas)
658 # fit EOD waveform:
659 try:
660 tas, _ = curve_fit(gaussian_sum, time, eod, tas)
661 except RuntimeError as e:
662 if verbose > 0:
663 print('Fit gaussian_sum failed in decompose_pulse():', e)
664 return pulse
665 # fit EOD waveform and spectrum:
666 bnds = [(1e-5, None) if k%3 == 2 else (None, None)
667 for k in range(len(tas))]
668 res = minimize(gaussian_sum_costs, tas,
669 args=(time, eod, freqs, energy), bounds=bnds)
670 if not res.success and verbose > 0:
671 print('warning: optimization gaussian_sum_costs failed in decompose_pulse():',
672 res.message)
673 else:
674 tas = res.x
675 # add another Gaussian:
676 rms_norm = np.max(np.abs(eod))
677 rms_old = np.sqrt(np.mean((gaussian_sum(time, *tas) - eod)**2))/rms_norm
678 eod_diff = np.abs(gaussian_sum(time, *tas) - eod)/rms_norm
679 if np.max(eod_diff) > 0.1:
680 if verbose > 1:
681 print(f'decompose_pulse(): added Gaussian because maximum rms error was {100*np.max(eod_diff):.0f}%')
682 ntas = np.concatenate((tas, (time[np.argmax(eod_diff)], np.max(eod_diff),
683 np.mean(tas[2::3]))))
684 bnds = [(1e-5, None) if k%3 == 2 else (None, None)
685 for k in range(len(ntas))]
686 res = minimize(gaussian_sum_costs, ntas,
687 args=(time, eod, freqs, energy), bounds=bnds)
688 if res.success:
689 rms_new = np.sqrt(np.mean((gaussian_sum(time, *res.x) - eod)**2))/rms_norm
690 if rms_new < 0.8*rms_old:
691 tas = res.x
692 elif verbose > 1:
693 print('decompose_pulse(): removed added Gaussian because it did not improve the fit')
694 elif verbose > 0:
695 print('warnong: optimization gaussian_sum_costs for additional Gaussian failed in decompose_pulse():',
696 res.message)
697 times = np.asarray(tas[0::3])
698 ampls = np.asarray(tas[1::3])
699 stdevs = np.asarray(tas[2::3])
700 pulse = dict(times=times, amplitudes=ampls, stdevs=stdevs)
701 return pulse
704def exp_decay(t, tau, ampl, offs):
705 """Exponential decay function.
707 x(t) = ampl*exp(-t/tau) + offs
709 Parameters
710 ----------
711 t: float or array
712 Time.
713 tau: float
714 Time constant of exponential decay.
715 ampl: float
716 Amplitude of exponential decay, i.e. initial value minus
717 steady-state value.
718 offs: float
719 Steady-state value.
721 Returns
722 -------
723 x: float or array
724 The exponential decay evaluated at times `t`.
726 """
727 return offs + ampl*np.exp(-t/tau)
730def analyze_pulse_tail(peak_index, eod, ratetime=None,
731 threshold=0.0, fit_frac=0.5, verbose=0):
732 """ Fit exponential to last peak/trough of pulse EOD.
734 Parameters
735 ----------
736 peak_index: int
737 Index of last peak in `eod`.
738 eod: 1-D or 2-D array
739 The eod waveform to be analyzed.
740 If an 1-D array, then this is the waveform and you
741 need to also pass a sampling rate in `rate`.
742 If a 2-D array, then first column is time in seconds and second
743 column is the eod waveform. Further columns are ignored.
744 ratetime: None or float or array of float
745 If a 1-D array is passed on to `eod` then either the sampling
746 rate in Hertz or the time array corresponding to `eod`.
747 threshold: float
748 Maximum noise level of the pulse waveform.
749 fit_frac: float or None
750 An exponential is fitted to the tail of the last peak/trough
751 starting where the waveform falls below this fraction of the
752 peak's height (0-1).
753 If None, do not attempt to fit.
754 verbose: int
755 Verbosity level passed for error and info messages.
757 Returns
758 -------
759 tau: float or np.nan
760 Time constant of pulse tail in seconds.
761 tstart: float or np.nan
762 Time where fit started in seconds.
763 fit: 1-D array of float or np.nan
764 Time trace of the fit corresponding to `eod`.
765 """
766 if fit_frac is None:
767 return np.nan, np.nan, None
768 if eod.ndim == 2:
769 time = eod[:, 0]
770 eod = eod[:, 1]
771 elif isinstance(ratetime, (list, tuple, np.ndarray)):
772 time = ratetime
773 else:
774 time = np.arange(len(eod))/rate
775 dt = np.mean(np.diff(time))
776 pi = peak_index
777 # positive or negative decay:
778 sign = 1
779 eodpp = eod[pi:] - 0.5*threshold
780 eodpn = -eod[pi:] - 0.5*threshold
781 if np.sum(eodpn[eodpn > 0]) > np.sum(eodpp[eodpp > 0]):
782 sign = -1
783 if sign*eod[pi] < 0:
784 pi += np.argmax(sign*eod[pi:])
785 pi_ampl = np.abs(eod[pi])
786 # no sufficiently large initial value:
787 sampl = sign*eod[pi]*fit_frac
788 if sampl <= threshold:
789 if verbose > 0:
790 print(f'exponential fit to tail of pulse failed: initial amplitude {sampl:.5f} smaller than threshold {threshold:.5f}')
791 return np.nan, np.nan, None
792 # no sufficiently long decay:
793 sinx = pi + np.argmax(sign*eod[pi:] < sampl)
794 n = 2*len(eod[sinx:])//3
795 if n < 10:
796 if verbose > 0:
797 print(f'exponential fit to tail of pulse failed: less than 10 samples {n}')
798 return np.nan, np.nan, None
799 # not decaying towards zero:
800 max_line = sampl - (sampl - threshold)*np.arange(n)/n + 1e-8
801 above_frac = np.sum(sign*eod[sinx:sinx + n] > max_line)/n
802 if above_frac > 0.05:
803 if verbose > 0:
804 print(f'exponential fit to tail of pulse failed: not decaying towards zero {100*above_frac:.1f}% > 5%')
805 return np.nan, np.nan, None
806 # estimate tau:
807 thresh = eod[sinx]*np.exp(-1.0)
808 tau_inx = np.argmax(sign*eod[sinx:] < sign*thresh)
809 if tau_inx < 2:
810 tau_inx = 2
811 rinx = sinx + 6*tau_inx
812 if rinx > len(eod) - 1:
813 rinx = len(eod) - 1
814 if rinx - sinx < 2*tau_inx:
815 if verbose > 0:
816 print(f'exponential fit to tail of pulse failed: less samples {rinx - sinx} than two time constants 3*{tau_inx}')
817 return np.nan, np.nan, None
818 tau = time[sinx + tau_inx] - time[sinx]
819 params = [tau, eod[sinx] - eod[rinx], eod[rinx]]
820 try:
821 popt, pcov = curve_fit(exp_decay, time[sinx:rinx] - time[sinx],
822 eod[sinx:rinx], params,
823 bounds=([0.0, -np.inf, -np.inf], np.inf))
824 except TypeError:
825 popt, pcov = curve_fit(exp_decay, time[sinx:rinx] - time[sinx],
826 eod[sinx:rinx], params)
827 if popt[0] > 1.2*tau:
828 tau_inx = int(np.round(popt[0]/dt))
829 rinx = sinx + 6*tau_inx
830 if rinx > len(eod) - 1:
831 rinx = len(eod) - 1
832 try:
833 popt, pcov = curve_fit(exp_decay, time[sinx:rinx] - time[sinx],
834 eod[sinx:rinx], popt,
835 bounds=([0.0, -np.inf, -np.inf], np.inf))
836 except TypeError:
837 popt, pcov = curve_fit(exp_decay, time[sinx:rinx] - time[sinx],
838 eod[sinx:rinx], popt)
839 tau = popt[0]
840 tstart = time[sinx]
841 fit = np.zeros(len(eod))
842 fit[:] = np.nan
843 fit[sinx:rinx] = exp_decay(time[sinx:rinx] - time[sinx], *popt)
844 if verbose > 0:
845 print(f'exponential fit to tail of pulse: got time constant {1000*tau:.3f}ms')
846 return tau, tstart, fit
849def pulse_spectrum(eod, ratetime=None, freq_resolution=1.0, fade_frac=0.0):
850 """Spectrum of a single pulse-type EOD.
852 Parameters
853 ----------
854 eod: 1-D or 2-D array
855 The eod waveform of which the spectrum is computed.
856 If an 1-D array, then this is the waveform and you
857 need to also pass a sampling rate in `rate`.
858 If a 2-D array, then first column is time in seconds and second
859 column is the eod waveform. Further columns are ignored.
860 ratetime: None or float or array of float
861 If a 1-D array is passed on to `eod` then either the sampling
862 rate in Hertz or the time array corresponding to `eod`.
863 freq_resolution: float
864 The frequency resolution of the spectrum.
865 fade_frac: float
866 Fraction of time of the EOD waveform that is used to fade in
867 and out to zero baseline.
869 Returns
870 -------
871 freqs: 1-D array of float
872 The frequency components of the energy spectrum.
873 energy: 1-D array of float
874 The energy spectrum of the single pulse EOD
875 with unit (x s)^2 = x^2 s/Hz.
876 The integral over the energy spectrum `np.sum(energy)*freqs[1]`
877 equals the integral over the squared eod, `np.sum(eod**2)/rate`.
878 That is, by making the energy spectrum a power sepctrum
879 (dividing the energy by the FFT window duration), the integral
880 over the power spectrum equals the mean-squared signal
881 (variance). But the single-pulse spectrum is not a power-spectrum.
882 because in the limit to infinitely long window, the power vanishes!
884 See Also
885 --------
886 thunderfish.fakefish.pulsefish_spectrum(): analytically computed spectra for simulated pulse EODs.
887 pulsetrain_spectrum(): power spectrum of train of pulse-type EODs.
889 """
890 if eod.ndim == 2:
891 rate = 1.0/(eod[1, 0] - eod[0, 0])
892 eod = eod[:, 1]
893 elif isinstance(ratetime, (list, tuple, np.ndarray)):
894 rate = 1.0/(ratetime[1] - ratetime[0])
895 else:
896 rate = ratetime
897 n_fft = nfft(rate, freq_resolution)
898 # subtract mean computed from the ends of the EOD snippet:
899 n = len(eod)//20 if len(eod) >= 20 else 1
900 eod = eod - 0.5*(np.mean(eod[:n]) + np.mean(eod[-n:]))
901 # zero padding:
902 max_idx = np.argmax(eod)
903 n0 = max_idx
904 n1 = len(eod) - max_idx
905 n = 2*max(n0, n1)
906 if n_fft < n:
907 n_fft = next_power_of_two(n)
908 data = np.zeros(n_fft)
909 data[n_fft//2 - n0:n_fft//2 + n1] = eod
910 # fade in and out:
911 if fade_frac > 0:
912 fn = int(fade_frac*len(eod))
913 data[n_fft//2 - n0:n_fft//2 - n0 + fn] *= np.arange(fn)/fn
914 data[n_fft//2 + n1 - fn:n_fft//2 + n1] *= np.arange(fn)[::-1]/fn
915 # spectrum:
916 dt = 1/rate
917 freqs = np.fft.rfftfreq(n_fft, dt)
918 fourier = np.fft.rfft(data)*dt
919 energy = 2*np.abs(fourier)**2 # one-sided spectrum!
920 return freqs, energy
923def pulsetrain(eod_times, eod, ratetime=None,
924 dwin=None, drate=None, fade_frac=0.0):
925 """Synthesize train of pulse-type EODs.
927 Place mean EOD waveform at `eod_times`.
929 Parameters
930 ----------
931 eod_times: 1-D array or None
932 List of times of detected EODs.
933 eod: 1-D or 2-D array
934 The EOD waveform used to synthezie the pulse train.
935 If an 1-D array, then this is the waveform and you
936 need to also pass a sampling rate in `rate`.
937 If a 2-D array, then first column is time in seconds and second
938 column is the EOD waveform. Further columns are ignored.
939 ratetime: None or float or array of float
940 If a 1-D array is passed on to `eod` then either the sampling
941 rate in Hertz or the time array corresponding to `eod`.
942 dwin: float or None
943 Duration of data window from which `eod_times` where extracted.
944 If None estimate from `eod_times`.
945 drate: float or None
946 Sampling rate of the generated data. If None set to sampling rate
947 of EOD waveform.
948 fade_frac: float
949 Fraction of time of the EOD waveform that is used to fade in
950 and out to zero baseline.
952 Returns
953 -------
954 data: 1-D array of float
955 Synthesized pulse train with samplingrate `drate`.
956 """
957 if eod.ndim == 2:
958 rate = 1.0/(eod[1, 0] - eod[0, 0])
959 time = eod[:, 0]
960 eod = eod[:, 1]
961 elif isinstance(ratetime, (list, tuple, np.ndarray)):
962 time = ratetime
963 rate = 1.0/(time[1] - time[0])
964 else:
965 time = np.arange(len(eod))/ratetime
967 # interpolate:
968 if drate is None:
969 drate = rate
970 dtime = time
971 else:
972 dtime = np.arange(time[0], time[-1] + 0.5/drate, 1/drate)
973 eod = np.interp(dtime, time, eod)
975 ipi = np.median(np.diff(eod_times))
976 n_ipi = int(ipi*drate)
977 if dwin is None:
978 n = int(eod_times[-1]*drate) + n_ipi//2
979 else:
980 n = int(dwin*drate)
981 data = np.zeros(n)
982 i0 = np.argmin(np.abs(dtime))
983 i1 = len(eod) - i0
984 fn = int(fade_frac*len(eod))
985 # place eod waveforms at eod_times:
986 prev_idx = -n_ipi
987 for t in eod_times:
988 idx = int(np.round(t*drate))
989 if idx >= n:
990 break
991 ii0 = i0 if idx - i0 >= 0 else idx
992 if idx - ii0 < (idx + prev_idx)//2:
993 ii0 = idx - (idx + prev_idx)//2
994 if ii0 < 0:
995 ii0 = 0
996 ii1 = i1 if idx + i1 < len(data) else len(data) - 1 - idx
997 data[idx - ii0:idx + ii1] = eod[i0 - ii0:i0 + ii1]
998 # fade in and out:
999 if fn > 0:
1000 data[idx - ii0:idx - ii0 + fn] *= np.arange(fn)/fn
1001 data[idx + ii1 - fn:idx + ii1] *= np.arange(fn)[::-1]/fn
1002 prev_idx = idx
1003 return data
1006def pulsetrain_spectrum(eod_times, eod, ratetime=None,
1007 dwin=None, drate=None,
1008 fade_frac=0.0, freq_resolution=1.0, **kwargs):
1009 """Power spectrum of train of pulse-type EODs.
1011 Place mean EOD waveform at `eod_times` and compute from that the
1012 power spectrum.
1014 Parameters
1015 ----------
1016 eod_times: 1-D array or None
1017 List of times of detected EODs.
1018 eod: 1-D or 2-D array
1019 The EOD waveform of which the spectrum is computed.
1020 If an 1-D array, then this is the waveform and you
1021 need to also pass a sampling rate in `rate`.
1022 If a 2-D array, then first column is time in seconds and second
1023 column is the EOD waveform. Further columns are ignored.
1024 ratetime: None or float or array of float
1025 If a 1-D array is passed on to `eod` then either the sampling
1026 rate in Hertz or the time array corresponding to `eod`.
1027 dwin: float or None
1028 Duration of data window from which `eod_times` where extracted.
1029 If None estimate from `eod_times`.
1030 drate: float or None
1031 Sampling rate of the generated data. If None set to sampling rate
1032 of EOD waveform.
1033 fade_frac: float
1034 Fraction of time of the EOD waveform that is used to fade in
1035 and out to zero baseline.
1036 freq_resolution: float
1037 The frequency resolution of the spectrum.
1038 kwargs: dict
1039 Further key-word arguments for `thunderlab.powerspectrum.psd()`.
1041 Returns
1042 -------
1043 freq: 1-D array of float
1044 Frequencies corresponding to power array.
1045 power: 1-D array of float
1046 Power spectral density in [eod]^2/Hz.
1048 """
1049 data = pulsetrain(eod_times, eod, ratetime=ratetime,
1050 dwin=dwin, drate=drate, fade_frac=fade_frac)
1051 freq, power = psd(data - np.mean(data), drate, freq_resolution, **kwargs)
1052 return freq, power
1055def analyze_pulse_spectrum(freqs, energy, eodf=None):
1056 """Analyze the spectrum of a pulse-type EOD.
1058 Parameters
1059 ----------
1060 freqs: 1-D array of float
1061 The frequency components of the energy spectrum.
1062 energy: 1-D array of float
1063 The energy spectrum of the single pulse-type EOD.
1064 eodf: float or None
1065 EOD frequency.
1067 Returns
1068 -------
1069 peak_freq: float
1070 Frequency at peak energy of the spectrum in Hertz.
1071 peak_energy: float
1072 Peak energy of the pulse spectrum in x^2 s/Hz.
1073 trough_freq: float
1074 Frequency at trough before peak in Hertz.
1075 trough_energy: float
1076 Energy of trough before peak in x^2 s/Hz.
1077 eodf_energy: float or None
1078 Energy at the EOD frequency if `eodf` is provided.
1079 att5: float
1080 Attenuation of average energy below 5 Hz relative to
1081 peak energy in decibel.
1082 att50: float
1083 Attenuation of average energy below 50 Hz relative to
1084 peak energy in decibel.
1085 low_cutoff: float
1086 Frequency at which the energy reached half of the
1087 peak energy relative to the DC energy in Hertz.
1088 high_cutoff: float
1089 3dB roll-off frequency in Hertz.
1090 """
1091 ip = np.argmax(energy)
1092 peak_freq = freqs[ip]
1093 peak_energy = energy[ip]
1094 it = np.argmin(energy[:ip]) if ip > 0 else 0
1095 trough_freq = freqs[it]
1096 trough_energy = energy[it]
1097 eodf_energy = None
1098 if eodf is not None:
1099 ie = np.argmin(np.abs(freqs - eodf))
1100 eodf_energy = energy[ie]
1101 att5 = decibel(np.mean(energy[freqs<5.0]), peak_energy)
1102 att50 = decibel(np.mean(energy[freqs<50.0]), peak_energy)
1103 low_cutoff = freqs[decibel(energy, peak_energy) > 0.5*att5][0]
1104 high_cutoff = freqs[decibel(energy, peak_energy) > -3.0][-1]
1105 return peak_freq, peak_energy, trough_freq, trough_energy, \
1106 eodf_energy, att5, att50, low_cutoff, high_cutoff
1109def analyze_pulse_intervals(eod_times, ipi_cv_thresh=0.5,
1110 ipi_percentile=30.0):
1111 """ Basic statistics of interpulse intervals.
1113 Parameters
1114 ----------
1115 eod_times: 1-D array or None
1116 List of times of detected EODs.
1117 ipi_cv_thresh: float
1118 If the coefficient of variation of the interpulse intervals
1119 (IPIs) is smaller than this threshold, then the statistics of
1120 IPIs is estimated from all IPIs. Otherwise only intervals
1121 smaller than a certain percentile are used.
1122 ipi_percentile: float
1123 When computing the statistics of IPIs from a subset of the
1124 IPIs, only intervals smaller than this percentile (between 0
1125 and 100) are used.
1127 Returns
1128 -------
1129 ipi_median: float
1130 Median inter-pulse interval.
1131 ipi_mean: float
1132 Mean inter-pulse interval.
1133 ipi_std: float
1134 Standard deviation of inter-pulse intervals.
1136 """
1137 if eod_times is None:
1138 return None, None, None
1139 inter_pulse_intervals = np.diff(eod_times)
1140 ipi_cv = np.std(inter_pulse_intervals)/np.mean(inter_pulse_intervals)
1141 if ipi_cv < ipi_cv_thresh:
1142 ipi_median = np.median(inter_pulse_intervals)
1143 ipi_mean = np.mean(inter_pulse_intervals)
1144 ipi_std = np.std(inter_pulse_intervals)
1145 else:
1146 intervals = inter_pulse_intervals[inter_pulse_intervals <
1147 np.percentile(inter_pulse_intervals, ipi_percentile)]
1148 ipi_median = np.median(intervals)
1149 ipi_mean = np.mean(intervals)
1150 ipi_std = np.std(intervals)
1151 return ipi_median, ipi_mean, ipi_std
1154def analyze_pulse(eod, ratetime=None, eod_times=None,
1155 baseline_mode='mode', min_pulse_win=0.001,
1156 start_end_thresh_fac=0.01, peak_thresh_frac=0.002,
1157 min_dist=50.0e-6, width_frac=0.5, fit_frac=0.5,
1158 freq_resolution=1.0, fade_frac=0.0,
1159 flip_pulse='none', fit_gaussians=True, ipi_cv_thresh=0.5,
1160 ipi_percentile=30.0, verbose=0):
1161 """Analyze the EOD waveform of a pulse fish.
1163 Parameters
1164 ----------
1165 eod: 1-D or 2-D array
1166 The eod waveform to be analyzed.
1167 If an 1-D array, then this is the waveform and you
1168 need to also pass a time array or sampling rate in `ratetime`.
1169 If a 2-D array, then first column is time in seconds, second
1170 column the EOD waveform, third column, if present, is the
1171 standard error of the EOD waveform. Further columns are
1172 optional and are not used.
1173 ratetime: None or float or array of float
1174 If a 1-D array is passed on to `eod` then either the sampling
1175 rate in Hertz or the time array corresponding to `eod`.
1176 eod_times: 1-D array or None
1177 List of times of detected EOD peaks.
1178 baseline_mode: "left", "right", "mean", "mode"
1179 Method to compute baseline voltage:
1180 "left": mean of leftmost `baseline_frac` data points.
1181 "right": mean of rightmost `baseline_frac` data points.
1182 "mean": mean of both leftmost and rightmost `baseline_frac` data points.
1183 "mode": mode of the distribution of the data values,
1184 i.e. most common value.
1185 min_pulse_win: float
1186 The minimum size of cut-out EOD waveform.
1187 start_end_thresh_fac: float
1188 Set the threshold for the start and end time to the p-p amplitude
1189 times this factor.
1190 peak_thresh_frac: float
1191 Set the threshold for peak and trough detection as a fraction
1192 of the p-p amplitude.
1193 min_dist: float
1194 Minimum distance between peak and troughs of the pulse.
1195 width_frac: float
1196 The width of a peak is measured at this fraction of a peak's
1197 height (0-1).
1198 fit_frac: float or None
1199 An exponential is fitted to the tail of the last peak/trough
1200 starting where the waveform falls below this fraction of the
1201 peak's height (0-1).
1202 freq_resolution: float
1203 The frequency resolution of the spectrum of the single pulse.
1204 fade_frac: float
1205 Fraction of time of the EOD waveform that is fade in and out
1206 to zero baseline.
1207 flip_pulse: 'auto', 'none', 'flip'
1208 - 'auto' flip waveform such that the first large extremum is positive.
1209 - 'flip' flip waveform.
1210 - 'none' do not flip waveform.
1211 fit_gaussians: bool
1212 If True, fit sum of Gaussians to pulse waveform.
1213 Otherwise, do not fit, and return empty dictionary for `pulses`.
1214 ipi_cv_thresh: float
1215 If the coefficient of variation of the interpulse intervals
1216 (IPIs) is smaller than this threshold, then the statistics of
1217 IPIs is estimated from all IPIs. Otherwise only intervals
1218 smaller than a certain percentile are used.
1219 ipi_percentile: float
1220 When computing the statistics of IPIs from a subset of the
1221 IPIs, only intervals smaller than this percentile (between 0
1222 and 100) are used.
1223 verbose: int
1224 Verbosity level passed for error and info messages.
1226 Returns
1227 -------
1228 meod: 2-D array of floats
1229 The eod waveform. First column is time in seconds,
1230 second column the eod waveform.
1231 Further columns are kept from the input `eod`.
1232 As the two last columns the waveform resulting from the
1233 decomposition into Gaussians and the fit to the tail of the
1234 last peak are appended.
1235 props: dict
1236 A dictionary with properties of the analyzed EOD waveform.
1238 - type: set to 'pulse'.
1239 - flipped: true if the waveform was flipped.
1240 - n: number of pulses analyzed (i.e. `len(eod_times)`),
1241 if provided.
1242 - times: the times of the detected EOD pulses
1243 (i.e. `eod_times`), if provided.
1244 - EODf: the inverse of the median interval between `eod_times`,
1245 if provided.
1246 - period: median interval between `eod_times`, if provided.
1247 - IPI-mean: mean interval between `eod_times`, if provided.
1248 - IPI-std: standard deviation of the intervals between
1249 `eod_times`, if provided.
1250 - IPI-CV: coefficient of variation of the intervals between
1251 `eod_times`, if provided.
1252 - aoffs: Offset that was subtracted from the average
1253 EOD waveform.
1254 - ppampl: peak-to-peak amplitude of the EOD waveform.
1255 - relmaxampl: amplitude of largest peak or trough,
1256 whichever is larger, relative to p-p amplitude.
1257 - noise: average standard error mean of the averaged
1258 EOD waveform relative to the p-p amplitude.
1259 - rmserror: root-mean-square error between fit with
1260 sum of Gaussians and EOD waveform relative to the
1261 p-p amplitude. Infinity if fit failed.
1262 - ppdist: distance between minimum and maximum phase
1263 in seconds.
1264 - startendthresh: threshold for start and end time is
1265 this factor times p-p-amplitude.
1266 - tstart: time in seconds where the pulse starts,
1267 i.e. crosses the threshold for the first time.
1268 - tend: time in seconds where the pulse ends,
1269 i.e. crosses the threshold for the last time.
1270 - width: total width of the pulse in seconds (tend-tstart).
1271 - totalarea: sum of areas under positive and negative peaks.
1272 - polaritybalance: contrast between areas under positive and
1273 negative phases.
1274 - median: median of the distribution of the absolute
1275 EOD waveform.
1276 - quartile1: first quartile of the distribution of the
1277 absolute EOD waveform.
1278 - quartile3: third quartile of the distribution of the
1279 absolute EOD waveform.
1280 - iq-range: inter-quartile range of the distribution of
1281 the absolute EOD waveform.
1282 - phasethresh: threshold for detecting peaks and troughs
1283 is this factor times p-p-amplitude.
1284 - nphases: number of phases of the single pulse.
1285 - firstphase: index of the first phase in the pulse
1286 (i.e. -1 for P-1)
1287 - lastphase: index of the last phase in the pulse
1288 (i.e. 3 for P3)
1289 - tau: time constant of exponential decay of pulse tail
1290 in seconds.
1291 - taustart: time where exponential fit started in seconds.
1292 - peakfreq: frequency at peak energy of the single pulse
1293 spectrum in Hertz.
1294 - peakenergy: peak energy of the single pulse spectrum.
1295 - troughfreq: frequency at trough before peak in Hertz.
1296 - troughenergy: energy of trough before peak in x^2 s/Hz.
1297 - eodfenergy: energy at EOD frequency in x^2 s/Hz.
1298 - energyatt5: attenuation of average energy below 5 Hz
1299 relative to peak energy in decibel.
1300 - energyatt50: attenuation of average energy below 50 Hz
1301 relative to peak energy in decibel.
1302 - lowcutoff: frequency at which the energy reached half of the
1303 peak energy relative to the DC energy in Hertz.
1304 - highcutoff: 3dB roll-off frequency of spectrum in Hertz.
1306 phases: dict
1307 Dictionary with
1309 - "indices": indices of each phase
1310 (1 is P1, i.e. the largest positive peak)
1311 - "times": times of each phase relative to P1 in seconds
1312 - "amplitudes": amplitudes of each phase
1313 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
1314 - "widths": widths of each phase at `width_frac` height
1315 - "areas": areas of each phase
1316 - "relareas": areas of the phases relative to total area
1317 - "zeros": time of zero crossing towards next phase in seconds
1319 Empty dictionary if waveform is not a pulse EOD.
1320 pulse: dict
1321 Dictionary with
1323 - "times": phase times in seconds,
1324 - "amplitudes": amplitudes, and
1325 - "stdevs": standard deviations in seconds
1327 of Gaussians fitted to the pulse waveform. Use the functions
1328 provided in thunderfish.fakefish to simulate pulse fish EODs
1329 from this data. Empty dictionary if `fit_gaussians` is False.
1330 energy: 2-D array
1331 The energy spectrum of a single pulse. First column are the
1332 frequencies, second column the energy in x^2 s/Hz.
1333 Empty if waveform is not a pulse EOD.
1335 """
1336 if eod.ndim == 2:
1337 if eod.shape[1] > 2:
1338 eeod = eod
1339 else:
1340 eeod = np.column_stack((eod, np.zeros(len(eod))))
1341 else:
1342 if isinstance(ratetime, (list, tuple, np.ndarray)):
1343 time = ratetime
1344 else:
1345 time = np.arange(len(eod))/ratetime
1346 eeod = np.zeros((len(eod), 3))
1347 eeod[:, 0] = time
1348 eeod[:, 1] = eod
1349 # storage:
1350 meod = np.zeros((eeod.shape[0], eeod.shape[1] + 2))
1351 meod[:, :eeod.shape[1]] = eeod
1352 meod[:, -2] = np.nan
1353 meod[:, -1] = np.nan
1355 # conditioning of the waveform:
1356 meod, toffs, aoffs, flipped, noise_thresh = \
1357 condition_pulse(meod, flip=flip_pulse, baseline_mode=baseline_mode,
1358 baseline_frac=0.05, large_phase_frac=0.2,
1359 min_pulse_win=min_pulse_win)
1361 # analysis of pulse waveform:
1362 pp_ampl, rel_max_ampl, dist, total_area, \
1363 polarity_balance, median, quartile1, quartile3 = \
1364 analyze_pulse_properties(noise_thresh, meod)
1366 # threshold for start and end time:
1367 start_end_thresh = pp_ampl*start_end_thresh_fac
1368 if start_end_thresh < 2*noise_thresh:
1369 start_end_thresh = 2*noise_thresh
1370 start_end_thresh_fac = start_end_thresh/pp_ampl if pp_ampl > 0 else 1
1372 # threshold for peak detection:
1373 peak_thresh = pp_ampl*peak_thresh_frac
1374 if peak_thresh < noise_thresh:
1375 peak_thresh = noise_thresh
1376 peak_thresh_frac = peak_thresh/pp_ampl if pp_ampl > 0 else 1
1378 # characterize EOD phases:
1379 tstart, tend, phases = analyze_pulse_phases(peak_thresh,
1380 start_end_thresh, meod,
1381 min_dist=min_dist,
1382 width_frac=width_frac)
1384 # fit exponential to last phase:
1385 tau = np.nan
1386 tau_start = np.nan
1387 if len(phases) > 0 and len(phases['times']) > 1:
1388 pi = np.argmin(np.abs(meod[:, 0] - phases['times'][-1]))
1389 tau, tau_start, fit = analyze_pulse_tail(pi, meod, None,
1390 threshold=noise_thresh,
1391 fit_frac=fit_frac,
1392 verbose=verbose)
1393 if fit is not None:
1394 meod[:, -1] = fit
1396 # analyze pulse intervals:
1397 ipi_median, ipi_mean, ipi_std = \
1398 analyze_pulse_intervals(eod_times, ipi_cv_thresh, ipi_percentile)
1399 eodf = 1/ipi_median if ipi_median is not None else None
1401 # energy spectrum of single EOD pulse:
1402 freqs, energy = pulse_spectrum(meod, None, freq_resolution, fade_frac)
1403 # store spectrum:
1404 eenergy = np.zeros((len(energy), 2))
1405 eenergy[:, 0] = freqs
1406 eenergy[:, 1] = energy
1407 # analyse spectrum:
1408 peakfreq, peakenergy, troughfreq, troughenergy, \
1409 eodfenergy, att5, att50, lowcutoff, highcutoff = \
1410 analyze_pulse_spectrum(freqs, energy, eodf)
1412 # decompose EOD waveform:
1413 rmserror = np.inf
1414 pulse = {}
1415 if fit_gaussians:
1416 pulse = decompose_pulse(meod, freqs, energy, phases,
1417 width_frac, verbose=verbose)
1418 if len(pulse) > 0:
1419 eod_fit = pulsefish_waveform(pulse, meod[:, 0])
1420 meod[:, -2] = eod_fit
1421 rmserror = np.sqrt(np.mean((meod[:, 1] - meod[:, -2])**2))/pp_ampl
1422 spec = pulsefish_spectrum(pulse, freqs)
1423 spec = np.abs(spec)**2
1424 eenergy = np.hstack((eenergy, spec.reshape((-1, 1))))
1426 # store properties:
1427 props = {}
1428 props['type'] = 'pulse'
1429 props['flipped'] = flipped
1430 if eod_times is not None:
1431 props['n'] = len(eod_times)
1432 props['times'] = eod_times + toffs
1433 props['EODf'] = eodf
1434 props['period'] = ipi_median
1435 props['IPI-mean'] = ipi_mean
1436 props['IPI-std'] = ipi_std
1437 props['IPI-CV'] = ipi_std/ipi_mean
1438 props['aoffs'] = aoffs
1439 props['ppampl'] = pp_ampl
1440 props['relmaxampl'] = rel_max_ampl
1441 if meod.shape[1] > 2:
1442 props['noise'] = np.mean(meod[:, 2])/pp_ampl if pp_ampl > 0 else 1
1443 props['rmserror'] = rmserror
1444 props['ppdist'] = dist
1445 props['startendthresh'] = start_end_thresh_fac
1446 props['tstart'] = tstart
1447 props['tend'] = tend
1448 props['width'] = tend - tstart
1449 props['totalarea'] = total_area
1450 props['polaritybalance'] = polarity_balance
1451 props['median'] = median
1452 props['quartile1'] = quartile1
1453 props['quartile3'] = quartile3
1454 props['iqrange'] = quartile3 - quartile1
1455 props['phasethresh'] = peak_thresh_frac
1456 props['nphases'] = len(phases['times']) if len(phases) > 0 else 0
1457 props['firstphase'] = phases['indices'][0] if len(phases) > 0 else 1
1458 props['lastphase'] = phases['indices'][-1] if len(phases) > 0 else 1
1459 props['tau'] = tau
1460 props['taustart'] = tau_start
1461 props['peakfreq'] = peakfreq
1462 props['peakenergy'] = peakenergy
1463 props['troughfreq'] = troughfreq
1464 props['troughenergy'] = troughenergy
1465 props['eodfenergy'] = eodfenergy
1466 props['energyatt5'] = att5
1467 props['energyatt50'] = att50
1468 props['lowcutoff'] = lowcutoff
1469 props['highcutoff'] = highcutoff
1471 return meod, props, phases, pulse, eenergy
1474def pulse_quality(props, max_clipped_frac=0.1, max_rms_sem=0.0):
1475 """Assess the quality of an EOD waveform of a pulse fish.
1477 Parameters
1478 ----------
1479 props: dict
1480 A dictionary with properties of the analyzed EOD waveform
1481 as returned by `analyze_pulse()`.
1482 max_clipped_frac: float
1483 Maximum allowed fraction of clipped data.
1484 max_rms_sem: float
1485 If not zero, maximum allowed standard error of the data
1486 relative to p-p amplitude.
1488 Returns
1489 -------
1490 skip_reason: str
1491 An empty string if the waveform is good, otherwise a string
1492 indicating the failure.
1493 msg: str
1494 A textual representation of the values tested.
1495 skipped_clipped: bool
1496 True if waveform was skipped because of clipping.
1497 """
1498 msg = []
1499 skip_reason = []
1500 skipped_clipped = False
1501 # clipped fraction:
1502 if 'clipped' in props:
1503 clipped_frac = props['clipped']
1504 msg += [f'clipped={100*clipped_frac:3.0f}%']
1505 if clipped_frac >= max_clipped_frac:
1506 skip_reason += [f'clipped={100*clipped_frac:3.0f}% (maximumClippedFraction={100*max_clipped_frac:3.0f}%)']
1507 skipped_clipped = True
1508 # noise:
1509 rms_sem = None
1510 if 'rmssem' in props:
1511 rms_sem = props['rmssem']
1512 if 'noise' in props:
1513 rms_sem = props['noise']
1514 if rms_sem is not None:
1515 msg += [f'rms sem waveform={100*rms_sem:6.2f}%']
1516 if max_rms_sem > 0.0 and rms_sem >= max_rms_sem:
1517 skip_reason += [f'noisy waveform s.e.m.={100*rms_sem:6.2f}% (maximumRMSNoise={100*max_rms_sem:6.2f}%)']
1518 return ', '.join(skip_reason), ', '.join(msg), skipped_clipped
1521def plot_pulse_eodtimes(ax, data, rate, width, eod_props,
1522 toffs=0, colors=None, markers=None,
1523 marker_size=10, legend_rows=8, **kwargs):
1524 """Mark pulse EODs in a plot of an EOD recording.
1526 Parameters
1527 ----------
1528 ax: matplotlib axes
1529 Axes used for plotting.
1530 data: 1D ndarray
1531 Recorded data (these are not plotted!).
1532 rate: float
1533 Sampling rate of the data in Hertz.
1534 width: float
1535 Minimum width of the data to be plotted in seconds.
1536 eod_props: list of dictionaries
1537 Lists of EOD properties as returned by `analyze_pulse()`
1538 and `analyze_wave()`. From the entries with 'type' ==
1539 'pulse' the properties 'EODf', 'peaktimes', and 'ppampl' are used.
1540 'EODf' is the averaged EOD frequency, 'peaktimes' is a list of
1541 detected EOD pulse times, and 'ppampl' is the peak-to-peak amplitude
1542 of the EOD waveform.
1543 Three new elements with keys 'color', 'marker', 'markersize' are added.
1544 These are the color, marker, andrelative size actually used to mark
1545 pulse times.
1546 toffs: float
1547 Time of first data value in seconds that will be added
1548 to the pulse times in `eod_props`.
1549 colors: list of colors or None
1550 If not None list of colors for plotting each cluster.
1551 markers: list of markers or None
1552 If not None list of markers for plotting each cluster.
1553 marker_size: float
1554 Size of markers used to mark the pulses.
1555 legend_rows: int
1556 Maximum number of rows to be used for the legend.
1557 kwargs:
1558 Key word arguments for the legend of the plot.
1560 Returns
1561 -------
1562 fish_dict: dict
1563 Dictionary that maps artist handles marking the pusle fish times
1564 with the "index" from `eod_props`.
1565 """
1566 amax = 0
1567 for eod in eod_props:
1568 if eod['type'] != 'pulse' or 'peaktimes' not in eod:
1569 continue
1570 if 'ppampl' in eod and amax < eod['ppampl']:
1571 amax = eod['ppampl']
1572 fish_indices = []
1573 fish_dict = {}
1574 k = 0
1575 for eod in eod_props:
1576 if eod['type'] != 'pulse':
1577 continue
1578 if 'peaktimes' not in eod:
1579 continue
1580 if 'ppampl' not in eod:
1581 continue
1582 x = eod['peaktimes'] + toffs
1583 pidx = np.round(eod['peaktimes']*rate).astype(int)
1584 mask = (pidx > 0) & (pidx < len(data))
1585 y = data[pidx[mask]]
1586 x = x[mask]
1587 style = dict(mec=None, mew=0)
1588 if colors is not None:
1589 style['color'] = colors[k%len(colors)]
1590 if markers is None:
1591 style['marker'] = 'o'
1592 else:
1593 style['marker'] = markers[k%len(markers)]
1594 style['ms'] = marker_size*(0.3 + 0.7*(min(1, eod['ppampl']/amax))**0.25)
1595 eod['color'] = style['color']
1596 eod['marker'] = style['marker']
1597 eod['markersize'] = style['ms']/marker_size
1598 label = f'{eod["EODf"]:6.1f} Hz'
1599 aa = ax.plot(x, y, linestyle='none', label=label,
1600 zorder=-1, **style)
1601 for a in aa:
1602 a.set_picker(8)
1603 fish_dict[a] = eod['index']
1604 fish_indices.append(eod['index'])
1605 k += 1
1606 # legend:
1607 if k > 0:
1608 if not 'ncol' in kwargs and legend_rows > 0:
1609 ncol = (k - 1) // legend_rows + 1
1610 kwargs['ncol'] = ncol
1611 if not 'numpoints' in kwargs:
1612 kwargs['numpoints'] = 1
1613 handles = ax.legend(**kwargs)
1614 lines = handles.get_lines()
1615 for j in range(len(lines)):
1616 a = lines[j]
1617 a.set_picker(8)
1618 fish_dict[a] = fish_indices[j]
1619 return fish_dict
1622def plot_pulse_rate(ax, eod_props, toffs=0, colors=None,
1623 markers=None, marker_size=10,
1624 legend_rows=8, **kwargs):
1625 """Plot pulse rates.
1627 Parameters
1628 ----------
1629 ax: matplotlib axes
1630 Axes used for plotting.
1631 eod_props: list of dictionaries
1632 Lists of EOD properties as returned by `analyze_pulse()`
1633 and `analyze_wave()`. From the entries with 'type' ==
1634 'pulse' the properties 'EODf' and 'peaktimes' are used. 'EODf'
1635 is the averaged EOD frequency, and 'peaktimes' is a list of
1636 detected EOD pulse times.
1637 toffs: float
1638 Time of first data value in seconds that will be added
1639 to the pulse times in `eod_props`.
1640 colors: list of colors or None
1641 If not None list of colors for plotting the rate of each cluster.
1642 markers: list of markers or None
1643 If not None list of markers for plotting each cluster.
1644 marker_size: float
1645 Size of markers used to mark the pulses.
1646 legend_rows: int
1647 Maximum number of rows to be used for the legend.
1648 kwargs:
1649 Key word arguments for the legend of the plot.
1651 Returns
1652 -------
1653 fish_dict: dict
1654 Dictionary that maps artist handles marking the pusle fish times
1655 with the "index" from `eod_props`.
1656 """
1657 fish_indices = []
1658 fish_dict = {}
1659 k = 0
1660 for props in eod_props:
1661 if props['type'] != 'pulse':
1662 continue
1663 if 'peaktimes' not in props:
1664 continue
1665 times = props['peaktimes'] + toffs
1666 rate = 1/np.diff(times)
1667 mask = np.append(True, np.abs(np.diff(rate))/rate[:-1] < 0.1)
1668 style = dict(mec=None, mew=0)
1669 if colors is not None:
1670 style['color'] = colors[k%len(colors)]
1671 if markers is None:
1672 style['marker'] = 'o'
1673 else:
1674 style['marker'] = markers[k%len(markers)]
1675 if marker_size is not None:
1676 style['ms'] = marker_size
1677 label = f'{props["EODf"]:6.1f} Hz'
1678 if 'IPI-CV' in props:
1679 cv = 100*props["IPI-CV"]
1680 label += f' CV={cv:.1f}%' if cv > 0.5 else f' CV={cv:.2f}%'
1681 aa = ax.plot(times[:-1][mask], rate[mask], label=label,
1682 linestyle='-', **style)
1683 ab = ax.plot(times[:-1][~mask], rate[~mask],
1684 linestyle='none', **style)
1685 for a in aa + ab:
1686 a.set_picker(8)
1687 fish_dict[a] = props['index']
1688 fish_indices.append(props['index'])
1689 k += 1
1690 ax.set_xlabel('Time [s]')
1691 ax.set_ylim(bottom=0)
1692 ax.set_ylabel('Rate [Hz]')
1693 # legend:
1694 if k > 0:
1695 if not 'ncol' in kwargs and legend_rows > 0:
1696 ncol = (k - 1) // legend_rows + 1
1697 kwargs['ncol'] = ncol
1698 if not 'numpoints' in kwargs:
1699 kwargs['numpoints'] = 1
1700 handles = ax.legend(**kwargs)
1701 lines = handles.get_lines()
1702 for j in range(len(lines)):
1703 a = lines[j]
1704 a.set_picker(8)
1705 fish_dict[a] = fish_indices[j]
1706 return fish_dict
1709def plot_pulse_eod(ax, eod_waveform, props, phases=None,
1710 unit=None, wave_periods=2,
1711 magnification_factor=20,
1712 wave_style=dict(lw=1.5, color='tab:red'),
1713 magnified_style=dict(lw=0.8, color='tab:red'),
1714 positive_style=dict(facecolor='tab:green', alpha=0.2,
1715 edgecolor='none'),
1716 negative_style=dict(facecolor='tab:blue', alpha=0.2,
1717 edgecolor='none'),
1718 sem_style=dict(color='0.8'),
1719 fit_style=dict(lw=1.5, color='tab:blue'),
1720 phase_style=dict(zorder=0, ls='', marker='o', color='tab:red',
1721 markersize=5, mec='white', mew=1),
1722 zerox_style=dict(zorder=50, ls='', marker='o', color='black',
1723 markersize=5, mec='white', mew=1),
1724 zero_style=dict(lw=0.5, color='0.7'),
1725 fontsize='medium'):
1726 """Plot and annotate a pulse-type EOD waveform.
1728 Parameters
1729 ----------
1730 ax: matplotlib axes
1731 Axes used for plotting.
1732 eod_waveform: 2-D array
1733 EOD waveform. First column is time in seconds, second column
1734 the (mean) eod waveform. The optional third column is the
1735 standard error, the optional fourth column is a fit of the
1736 whole waveform, and the optional fourth column is a fit of
1737 the tails of a pulse waveform.
1738 props: dict
1739 A dictionary with properties of the analyzed EOD waveform as
1740 returned by `analyze_wave()` and `analyze_pulse()`.
1741 phases: dict
1742 Dictionary with phase properties as returned by
1743 `analyze_pulse_phases()`, `analyze_pulse()`, and
1744 `load_pulse_phases()`.
1745 unit: str
1746 Optional unit of the data used for y-label.
1747 wave_periods: float
1748 How many periods of a wave EOD are shown.
1749 magnification_factor: float
1750 If larger than one, plot a magnified version of the EOD
1751 waveform magnified by this factor.
1752 wave_style: dict
1753 Arguments passed on to the plot command for the EOD waveform.
1754 magnified_style: dict
1755 Arguments passed on to the plot command for the magnified EOD waveform.
1756 positive_style: dict
1757 Arguments passed on to the fill_between command for coloring
1758 positive phases.
1759 negative_style: dict
1760 Arguments passed on to the fill_between command for coloring
1761 negative phases.
1762 sem_style: dict
1763 Arguments passed on to the fill_between command for the
1764 standard error of the EOD.
1765 fit_style: dict
1766 Arguments passed on to the plot command for the fitted EOD.
1767 phase_style: dict
1768 Arguments passed on to the plot command for marking EOD phases.
1769 zerox_style: dict
1770 Arguments passed on to the plot command for marking zero crossings.
1771 zero_style: dict
1772 Arguments passed on to the plot command for the zero line.
1773 fontsize: str or float or int
1774 Fontsize for annotation text.
1776 """
1777 time = 1000*eod_waveform[:, 0]
1778 eod = eod_waveform[:, 1]
1779 # time axis:
1780 # width of maximum peak:
1781 meod = np.abs(eod)
1782 ip = np.argmax(meod)
1783 thresh = 0.5*meod[ip]
1784 i0 = ip - np.argmax(meod[ip::-1] < thresh)
1785 i1 = ip + np.argmax(meod[ip:] < thresh)
1786 w = 4*(time[i1] - time[i0])
1787 w = np.ceil(w/0.5)*0.5
1788 # make sure tstart, tend, and time constant are included:
1789 if props is not None:
1790 if 'tstart' in props and 1000*props['tstart'] < -w:
1791 w = np.ceil(abs(1000*props['tstart'])/0.5)*0.5
1792 if 'tend' in props and 1000*props['tend'] > 2*w:
1793 w = np.ceil(0.5*abs(1000*props['tend'])/0.5)*0.5
1794 if 'taustart' in props and props['taustart'] is not None and \
1795 1100*props['taustart'] > 2*w:
1796 w = np.ceil(0.5*abs(1100*props['taustart'])/0.5)*0.5
1797 # set xaxis limits:
1798 xlim_l = -w
1799 xlim_r = 2*w
1800 xlim = (xlim_r - xlim_l)/2
1801 ax.set_xlim(xlim_l, xlim_r)
1802 if xlim < 2:
1803 ax.xaxis.set_major_locator(MultipleLocator(0.5))
1804 elif xlim < 4:
1805 ax.xaxis.set_major_locator(MultipleLocator(1))
1806 elif xlim < 8:
1807 ax.xaxis.set_major_locator(MultipleLocator(2))
1808 ax.set_xlabel('Time [msec]')
1809 # amplitude axis:
1810 ylim = np.max(np.abs(eod[(time >= xlim_l) & (time <= xlim_r)]))
1811 ax.set_ylim(-1.15*ylim, +1.15*ylim)
1812 if unit:
1813 ax.set_ylabel(f'Amplitude [{unit}]')
1814 else:
1815 ax.set_ylabel('Amplitude')
1816 # ax height dimensions:
1817 t = ax.text(0, 0, 'test', fontsize=fontsize)
1818 fs = t.get_fontsize()
1819 t.remove()
1820 pixelx = np.abs(np.diff(ax.get_window_extent().get_points()[:, 0]))[0]
1821 dxu = 2*xlim/pixelx
1822 xfs = fs*dxu
1823 pixely = np.abs(np.diff(ax.get_window_extent().get_points()[:, 1]))[0]
1824 dyu = 2*ylim/pixely
1825 yfs = fs*dyu
1826 texts = []
1827 quadrants = np.zeros((2, 2), dtype=int)
1828 # magnification threshold:
1829 if magnification_factor > 1:
1830 mag_thresh = 0.95*np.max(np.abs(eod))/magnification_factor
1831 else:
1832 mag_thresh = 0
1833 # plot zero line:
1834 ax.axhline(0.0, zorder=10, **zero_style)
1835 # plot areas:
1836 if phases is not None and len(phases) > 0:
1837 if positive_style is not None and len(positive_style) > 0:
1838 ax.fill_between(time, eod, 0, eod >= 0, zorder=0,
1839 **positive_style)
1840 if negative_style is not None and len(negative_style) > 0:
1841 ax.fill_between(time, eod, 0, eod <= 0, zorder=0,
1842 **negative_style)
1843 # plot Fourier/Gaussian fit:
1844 if eod_waveform.shape[1] > 3 and np.all(np.isfinite(eod_waveform[:, 3])):
1845 ax.plot(time, eod_waveform[:, 3], zorder=30, **fit_style)
1846 # plot time constant fit:
1847 tau_magnified = False
1848 if eod_waveform.shape[1] > 4:
1849 if np.nanmax(np.abs(eod_waveform[:, 4])) < mag_thresh:
1850 tau_magnified = True
1851 ax.plot(time, magnification_factor*eod_waveform[:, 4],
1852 zorder=35, **fit_style)
1853 else:
1854 fs = dict(**fit_style)
1855 if 'lw' in fs:
1856 fs['lw'] *= 2
1857 ax.plot(time, eod_waveform[:, 4], zorder=35, **fs)
1858 # plot waveform:
1859 ax.plot(time, eod, zorder=45, **wave_style)
1860 # plot standard error:
1861 if eod_waveform.shape[1] > 2:
1862 std_eod = eod_waveform[:, 2]
1863 ax.fill_between(time, eod + std_eod, eod - std_eod,
1864 zorder=20, **sem_style)
1865 # plot magnified pulse waveform:
1866 magnification_mask = np.zeros(len(time), dtype=bool)
1867 if magnification_factor > 1 and phases is not None and len(phases) > 0:
1868 i0 = np.argmax(np.abs(eod) > mag_thresh)
1869 if i0 > 0:
1870 left_eod = magnification_factor*eod[:i0]
1871 magnification_mask[:i0] = True
1872 ax.plot(time[:i0], left_eod, zorder=40, **magnified_style)
1873 ta = None
1874 if left_eod[-1] > 0:
1875 it = np.argmax(left_eod > 0.95*np.max(eod))
1876 if it < len(left_eod)//2:
1877 it = len(left_eod) - 1
1878 if time[it] - 3*xfs > xlim_l:
1879 ty = left_eod[it] if left_eod[it] < np.max(eod) else np.max(eod)
1880 ta = ax.text(time[it], ty, f'x{magnification_factor:.0f} ',
1881 ha='right', va='top', zorder=100,
1882 fontsize=fontsize)
1883 if ty > 0.5*ylim:
1884 quadrants[0, 0] += 1
1885 else:
1886 it = np.argmax(left_eod < 0.95*np.min(eod))
1887 if it < len(left_eod)//2:
1888 it = len(left_eod) - 1
1889 if time[it] - 3*xfs > xlim_l:
1890 ty = left_eod[it] if left_eod[it] > np.min(eod) else np.min(eod)
1891 ta = ax.text(time[it], ty, f'x{magnification_factor:.0f} ',
1892 ha='right', va='bottom', zorder=100,
1893 fontsize=fontsize)
1894 if ty < -0.5*ylim:
1895 quadrants[1, 0] += 1
1896 if ta is not None:
1897 texts.append(ta)
1898 left_snip = left_eod[time[:i0] < 0.1*xlim_l]
1899 if len(left_snip) > 0:
1900 if quadrants[0, 0] == 0:
1901 quadrants[0, 0] += np.max(left_snip) > 0.5*ylim
1902 if quadrants[1, 0] == 0:
1903 quadrants[1, 0] += np.min(left_snip) < -0.5*ylim
1904 i1 = len(eod) - np.argmax(np.abs(eod[::-1]) > mag_thresh)
1905 right_eod = magnification_factor*eod[i1:]
1906 magnification_mask[i1:] = True
1907 ax.plot(time[i1:], right_eod, zorder=40, **magnified_style)
1908 right_snip = right_eod[time[i1:] > 0.4*xlim_r]
1909 if len(right_snip) > 0:
1910 quadrants[0, 1] += np.max(right_snip) > 0.5*ylim
1911 quadrants[1, 1] += np.min(right_snip) < -0.5*ylim
1912 # annotate time constant fit:
1913 tau = None if props is None else props.get('tau', None)
1914 if tau is not None and eod_waveform.shape[1] > 4:
1915 if tau < 0.001:
1916 label = f'\u03c4={1.e6*tau:.0f}\u00b5s'
1917 else:
1918 label = f'\u03c4={1.e3*tau:.2f}ms'
1919 inx = np.argmin(np.isnan(eod_waveform[:, 4]))
1920 x0 = time[inx]
1921 x = x0 + 1000*np.log(2)*tau
1922 if x + 4*xfs >= xlim_r:
1923 if xlim_r - x0 >= 4*xfs:
1924 x = xlim_r - 8*xfs
1925 else:
1926 x = x0
1927 elif x + 8*xfs > xlim_r:
1928 x = xlim_r - 8*xfs
1929 if x < x0:
1930 x = x0
1931 y = eod_waveform[np.argmin(np.abs(time - x)), 4]
1932 if tau_magnified:
1933 y *= magnification_factor
1934 va = 'bottom' if eod[inx] > 0 else 'top'
1935 if eod[inx] < 0:
1936 y -= 0.5*yfs
1937 ta = ax.text(x + xfs, y, label, ha='left', va=va,
1938 zorder=100, fontsize=fontsize)
1939 texts.append(ta)
1940 if x + xfs > 0.4*xlim_r:
1941 if y > 0.5*ylim:
1942 quadrants[0, 1] += 1
1943 elif y < -0.5*ylim:
1944 quadrants[1, 1] += 1
1945 if props is not None:
1946 # mark start and end:
1947 if 'tstart' in props:
1948 ax.axvline(1000*props['tstart'], 0.45, 0.55,
1949 color='k', lw=0.5, zorder=25)
1950 if 'tend' in props:
1951 ax.axvline(1000*props['tend'], 0.45, 0.55,
1952 color='k', lw=0.5, zorder=25)
1953 # mark cumulative:
1954 if 'median' in props:
1955 y = -1.07*ylim
1956 m = 1000*props['median']
1957 q1 = 1000*props['quartile1']
1958 q3 = 1000*props['quartile3']
1959 w = q3 - q1
1960 ax.plot([q1, q3], [y, y], 'gray', lw=4, zorder=25)
1961 ax.plot(m, y, 'o', color='white', ms=3, zorder=26)
1962 label = f'{w:.2f}ms' if w >= 1 else f'{1000*w:.0f}\u00b5s'
1963 ax.text(q3 + xfs, y, label,
1964 va='center', zorder=100, fontsize=fontsize)
1965 # plot and annotate phases:
1966 if phases is not None and len(phases) > 0:
1967 upper_area_text = False
1968 lower_area_text = False
1969 # mark zero crossings:
1970 zeros = 1000*phases['zeros']
1971 ax.plot(zeros, np.zeros(len(zeros)), **zerox_style)
1972 # phase peaks and troughs:
1973 max_peak_idx = np.argmax(phases['amplitudes'])
1974 min_trough_idx = np.argmin(phases['amplitudes'])
1975 for i in range(len(phases['times'])):
1976 index = phases['indices'][i]
1977 ptime = 1000*phases['times'][i]
1978 if ptime < xlim_l or ptime > xlim_r:
1979 continue
1980 pi = np.argmin(np.abs(time - ptime))
1981 mfac = magnification_factor if magnification_mask[pi] else 1
1982 pampl = mfac*phases['amplitudes'][i]
1983 relampl = phases['relamplitudes'][i]
1984 relarea = phases['relareas'][i]
1985 # classify phase:
1986 ampl_phase = phases['amplitudes'][i]
1987 ampl_left = phases['amplitudes'][i - 1] if i > 0 else 0
1988 ampl_right = phases['amplitudes'][i + 1] if i + 1 < len(phases['amplitudes']) else 0
1989 local_maximum = ampl_phase > ampl_left and ampl_phase > ampl_right
1990 if local_maximum:
1991 right_phase = (i >= max_peak_idx)
1992 min_max_phase = (i == max_peak_idx)
1993 local_phase = (ampl_phase < 0)
1994 else:
1995 right_phase = i >= min_trough_idx
1996 min_max_phase = (i == min_trough_idx)
1997 local_phase = (ampl_phase > 0)
1998 sign = np.sign(pampl)
1999 # mark phase peak/trough:
2000 ax.plot(ptime, pampl, **phase_style)
2001 # text for phase label:
2002 label = f'P{index:.0f}'
2003 if index != 1 and not local_phase:
2004 if np.abs(ptime) < 1:
2005 ts = f'{1000*ptime:.0f}\u00b5s'
2006 elif np.abs(ptime) < 10:
2007 ts = f'{ptime:.2f}ms'
2008 else:
2009 ts = f'{ptime:.3g}ms'
2010 if np.abs(relampl) < 0.05:
2011 ps = f'{100*relampl:.1f}%'
2012 else:
2013 ps = f'{100*relampl:.0f}%'
2014 label += f'({ps} @ {ts})'
2015 # position of phase label:
2016 ltime = ptime
2017 lampl = pampl
2018 valign = 'top' if sign < 0 else 'baseline'
2019 add = True
2020 if local_phase or (min_max_phase and abs(pampl)/ylim < 0.8):
2021 halign = 'center'
2022 dx = 0
2023 dy = 0.6*yfs
2024 if local_phase:
2025 add = False
2026 elif min_max_phase:
2027 halign = 'left' if right_phase else 'right'
2028 dx = xfs if right_phase else -xfs
2029 dy = 0
2030 if abs(relampl) > 0.85:
2031 dx *= 2
2032 dy = -1.5*yfs
2033 else:
2034 dx = 0
2035 dy = 0.8*yfs
2036 if right_phase:
2037 halign = 'left'
2038 if i > 0 and np.isfinite(phases['zeros'][i - 1]):
2039 ltime = 1000*phases['zeros'][i - 1]
2040 else:
2041 dx = -2*xfs
2042 #np.sum(phases['amplitudes'][i + 1:]*pampl > 0)
2043 else:
2044 halign = 'right'
2045 if np.isfinite(phases['zeros'][i]):
2046 ltime = 1000*phases['zeros'][i]
2047 else:
2048 dx = 2*xfs
2049 if sign < 0:
2050 dy = -dy
2051 ta = ax.text(ltime + dx, lampl + dy, label,
2052 ha=halign, va=valign, zorder=100, fontsize=fontsize)
2053 if add:
2054 texts.append(ta)
2055 # area:
2056 if np.abs(relarea) < 0.01:
2057 continue
2058 elif np.abs(relarea) < 0.05:
2059 label = f'{100*relarea:.1f}%'
2060 else:
2061 label = f'{100*relarea:.0f}%'
2062 x = ptime
2063 if i > 0 and i < len(phases['times']) - 1:
2064 xl = 1000*phases['times'][i - 1]
2065 xr = 1000*phases['times'][i + 1]
2066 tsnippet = time[(time > xl) & (time < xr)]
2067 snippet = eod[(time > xl) & (time < xr)]
2068 tsnippet = tsnippet[np.sign(pampl)*snippet > 0]
2069 snippet = snippet[np.sign(pampl)*snippet > 0]
2070 x = np.sum(tsnippet*snippet)/np.sum(snippet)
2071 if abs(relampl) > 0.5:
2072 ax.text(x, sign*0.6*yfs, label,
2073 rotation='vertical',
2074 va='top' if sign < 0 else 'bottom',
2075 ha='center', zorder=35, fontsize=fontsize)
2076 elif abs(relampl) > 0.25 and abs(relarea) > 0.19:
2077 ax.text(x, sign*0.6*yfs, label,
2078 va='top' if sign < 0 else 'baseline',
2079 ha='center', zorder=35, fontsize=fontsize)
2080 else:
2081 dx = 0.5*xfs if right_phase else -0.5*xfs
2082 ta = ax.text(ltime + dx, -sign*0.4*yfs, label,
2083 va='baseline' if sign < 0 else 'top',
2084 ha=halign, zorder=100, fontsize=fontsize)
2085 if -sign > 0 and not upper_area_text:
2086 texts.append(ta)
2087 upper_area_text = True
2088 if -sign < 0 and not lower_area_text:
2089 texts.append(ta)
2090 lower_area_text = True
2091 # arrange text vertically to avoid overlaps:
2092 ul_texts = []
2093 ur_texts = []
2094 ll_texts = []
2095 lr_texts = []
2096 for t in texts:
2097 x, y = t.get_position()
2098 if y > 0:
2099 if x >= phases['times'][max_peak_idx]:
2100 ur_texts.append(t)
2101 else:
2102 ul_texts.append(t)
2103 else:
2104 if x >= phases['times'][min_trough_idx]:
2105 lr_texts.append(t)
2106 else:
2107 ll_texts.append(t)
2108 for ts, (j, k) in zip([ul_texts, ur_texts, ll_texts, lr_texts],
2109 [(0, 0), (0, 1), (1, 0), (1, 1)]):
2110 if len(ts) > 1:
2111 ys = []
2112 for t in ts:
2113 # alternative:
2114 #renderer = ax.get_fig().canvas.renderer
2115 #bbox = t.get_window_extent(renderer).transformed(ax.transData.inverted())
2116 x, y = t.get_position()
2117 ys.append(abs(y))
2118 idx = np.argsort(ys)
2119 x, y = ts[idx[0]].get_position()
2120 yp = abs(y)
2121 for i in idx[1:]:
2122 t = ts[i]
2123 x, y = t.get_position()
2124 s = t.get_text()
2125 if abs(y) < abs(yp) + 2*yfs and \
2126 len(s) > 4 and s[:2] != '\u03c4=':
2127 y = np.sign(y)*(abs(yp) + 2*yfs)
2128 t.set_y(y)
2129 if len(s) >= 4 and abs(y) > 0.5*ylim:
2130 quadrants[j, k] += 1
2131 yp = y
2132 # annotate plot:
2133 if unit is None or len(unit) == 0 or unit == 'a.u.':
2134 unit = ''
2135 if props is not None:
2136 label = '' # f'p-p amplitude = {props["p-p-amplitude"]:.3g} {unit}\n'
2137 if 'n' in props:
2138 eods = 'EODs' if props['n'] > 1 else 'EOD'
2139 segs = ''
2140 if 'nsegments' in props:
2141 n_segs = props['nsegments']
2142 segs = f' (in {n_segs} segment{"s" if n_segs > 1 else ""})'
2143 label += f'n = {props["n"]} {eods}{segs}\n'
2144 if 'flipped' in props and props['flipped']:
2145 label += 'flipped\n'
2146 if 'polaritybalance' in props:
2147 label += f'PB={100*props["polaritybalance"]:.0f} %\n'
2148 # weigh left quadrants less:
2149 quadrants *= 2
2150 quadrants[quadrants[:, 1] > 0, 1] -= 1
2151 # find free quadrant:
2152 q_row, q_col = np.unravel_index(np.argmin(quadrants), quadrants.shape)
2153 # place text in quadrant:
2154 y = 1.05 if q_row == 0 else 0
2155 va = 'top' if q_row == 0 else 'bottom'
2156 x = 0.03 if q_col == 0 else 0.97
2157 ha = 'left' if q_col == 0 else 'right'
2158 ax.text(x, y, label, transform=ax.transAxes,
2159 va=va, ha=ha, zorder=100)
2162def plot_pulse_spectrum(ax, energy, props, min_freq=1.0, max_freq=10000.0,
2163 min_db = -60, ref_energy=None,
2164 spec_style=dict(lw=3, color='tab:blue'),
2165 analytic_style=dict(lw=4, color='tab:cyan'),
2166 peak_style=dict(ls='', marker='o', markersize=6,
2167 color='tab:blue', mec='none', mew=0,
2168 alpha=0.4),
2169 cutoff_style=dict(lw=1, ls='-', color='0.5'),
2170 att5_color='0.8', att50_color='0.9',
2171 fontsize='medium'):
2172 """Plot and annotate spectrum of single pulse EOD.
2174 Parameters
2175 ----------
2176 ax: matplotlib axes
2177 Axes used for plotting.
2178 energy: 2-D array
2179 The energy spectrum of a single pulse as returned by `analyze_pulse()`.
2180 First column are the frequencies, second column the energy.
2181 An optional third column is an analytically computed spectrum.
2182 props: dict
2183 A dictionary with properties of the analyzed EOD waveform as
2184 returned by `analyze_pulse()`.
2185 min_freq: float
2186 Minimun frequency of the spectrum to be plotted (logscale!).
2187 max_freq: float
2188 Maximun frequency of the spectrum to be plotted (logscale!).
2189 min_db: float
2190 Minimum decibel level shown.
2191 ref_energy: float or None
2192 Reference energy for computing decibel.
2193 If None take the maximum energy.
2194 spec_style: dict
2195 Arguments passed on to the plot command for the energy spectrum
2196 computed from the data.
2197 analytic_style: dict
2198 Arguments passed on to the plot command for the energy spectrum
2199 that was analytically computed from the Gaussian fits
2200 (optional third column in `energy`).
2201 peak_style: dict
2202 Arguments passed on to the plot commands for marking the peak
2203 and trough frequency.
2204 cutoff_style: dict
2205 Arguments passed on to the plot command for the lines marking
2206 cutoff frequencies.
2207 att5_color: matplotlib color specification
2208 Color for the rectangular patch marking the first 5 Hz.
2209 att50_color: matplotlib color specification
2210 Color for the rectangular patch marking the first 50 Hz.
2211 fontsize: str or float or int
2212 Fontsize for annotation text.
2214 Returns
2215 -------
2216 ref_energy: float
2217 The actually used reference energy for computing decibels.
2218 """
2219 ax.axvspan(1, 50, color=att50_color, zorder=-20)
2220 att = props['energyatt50']
2221 if att < -10:
2222 y = att + 1
2223 if y < min_db:
2224 y = min_db + 2
2225 ax.text(10, y, f'{att:.0f}dB',
2226 ha='left', va='bottom', zorder=100, fontsize=fontsize)
2227 else:
2228 ax.text(10, att - 1, f'{att:.0f}dB',
2229 ha='left', va='top', zorder=100, fontsize=fontsize)
2230 ax.axvspan(1, 5, color=att5_color, zorder=-10)
2231 att = props['energyatt5']
2232 if att < -10:
2233 y = att + 1
2234 if y < min_db:
2235 y = min_db + 2
2236 ax.text(4, y, f'{att:.0f}dB',
2237 ha='right', va='bottom', zorder=100, fontsize=fontsize)
2238 else:
2239 ax.text(4, att - 1, f'{att:.0f}dB',
2240 ha='right', va='top', zorder=100, fontsize=fontsize)
2241 lowcutoff = props['lowcutoff']
2242 if lowcutoff >= min_freq:
2243 ax.plot([lowcutoff, lowcutoff, 1], [min_db, 0.5*att, 0.5*att],
2244 zorder=30, **cutoff_style)
2245 ax.text(1.2*lowcutoff, 0.5*att - 1, f'{lowcutoff:.0f}Hz',
2246 ha='left', va='top', zorder=100, fontsize=fontsize)
2247 highcutoff = props['highcutoff']
2248 ax.plot([highcutoff, highcutoff], [min_db, -3], zorder=30, **cutoff_style)
2249 ax.text(1.2*highcutoff, -3, f'{highcutoff:.0f}Hz',
2250 ha='left', va='center', zorder=100, fontsize=fontsize)
2251 if ref_energy is None:
2252 ref_energy = np.max(energy[:, 1])
2253 if 'EODf' in props and 'eodfenergy' in props and \
2254 props['eodfenergy'] is not None:
2255 eodf = props['EODf']
2256 db = decibel(props['eodfenergy'], ref_energy)
2257 ax.plot([eodf, eodf], [min_db, db], zorder=30, **cutoff_style)
2258 ax.plot(eodf, db, zorder=60, **peak_style)
2259 if lowcutoff < eodf or lowcutoff > 2*eodf:
2260 ax.text(1.1*eodf, 0.5*(min_db + db),
2261 f'EODf={eodf:.0f}Hz: {db:.0f}dB',
2262 rotation='vertical', va='center', fontsize=fontsize)
2263 else:
2264 ax.text(0.95*eodf, min_db + 2, f'EODf={eodf:.0f}Hz: {db:.0f}dB',
2265 rotation='vertical', ha='right', va='bottom',
2266 fontsize=fontsize)
2267 if energy.shape[1] > 2 and np.all(np.isfinite(energy[:, 2])) and len(analytic_style) > 0:
2268 db = decibel(energy[:, 2], ref_energy)
2269 ax.plot(energy[:, 0], db, zorder=45, **analytic_style)
2270 db = decibel(energy[:, 1], ref_energy)
2271 ax.plot(energy[:, 0], db, zorder=50, **spec_style)
2272 peakfreq = props['peakfreq']
2273 if peakfreq >= min_freq:
2274 ax.plot(peakfreq, 0, zorder=60, clip_on=False, **peak_style)
2275 ax.text(peakfreq*1.2, 1, f'{peakfreq:.0f}Hz',
2276 va='bottom', zorder=100, fontsize=fontsize)
2277 troughfreq = props['troughfreq']
2278 if troughfreq >= min_freq:
2279 troughenergy = decibel(props['troughenergy'], ref_energy)
2280 ax.plot(troughfreq, troughenergy, zorder=60, **peak_style)
2281 ax.text(troughfreq, troughenergy - 3,
2282 f'{troughenergy:.0f}dB @ {troughfreq:.0f}Hz',
2283 ha='center', va='top', zorder=100, fontsize=fontsize)
2284 ax.set_xlim(min_freq, max_freq)
2285 ax.set_xscale('log')
2286 ax.set_ylim(min_db, 2)
2287 ax.set_xlabel('Frequency [Hz]')
2288 ax.set_ylabel('Energy [dB]')
2289 return ref_energy
2292def plot_pulse_distribution(eod_waveform):
2293 """ Plot EOD waveform, amplitude distribution and cumulative.
2295 This is for debug. It shows a new figure.
2297 Parameters
2298 ----------
2299 eod_waveform: 2-D array
2300 EOD waveform. First column is time in seconds, second column
2301 the (mean) eod waveform.
2302 """
2303 time = 1000*eod_waveform[:, 0]
2304 eod = eod_waveform[:, 1]
2305 sdata = np.sort(eod)
2306 icuml = np.arange(len(sdata))/len(sdata)
2307 n = max(10, len(sdata)//100)
2308 order = 1
2309 filtered = savgol_filter(sdata, n, order, 0)
2310 derived = savgol_filter(sdata, n, order, 1)
2311 dsd = sdata[-1] - sdata[0]
2312 mask = (sdata > sdata[0] + 0.01*dsd) & (sdata < sdata[0] + 0.99*dsd)
2313 baseline = sdata[mask][np.argmin(derived[mask])]
2314 fig, axs = plt.subplots(1, 4, sharey=True)
2315 axs[0].plot(time, eod)
2316 axs[0].axhline(baseline, color='k')
2317 axs[0].set_xlabel('Time [ms]')
2318 axs[0].set_ylabel('Amplitude')
2319 axs[1].hist(eod, bins=len(eod)//10, orientation='horizontal')
2320 axs[1].axhline(baseline, color='k')
2321 axs[1].set_xlabel('Histogram')
2322 axs[2].plot(icuml, sdata, '-o', ms=4)
2323 #axs[2].plot(icuml, iampl)
2324 axs[2].plot(icuml, filtered)
2325 axs[2].axhline(baseline, color='k')
2326 axs[2].set_xlabel('Cumulative')
2327 axs[3].plot(derived, sdata)
2328 axs[3].axhline(baseline, color='k')
2329 axs[3].set_xlabel('Derivative')
2330 plt.show()
2333def save_pulse_fish(eod_props, unit, basename, **kwargs):
2334 """Save properties of pulse EODs to file.
2336 Parameters
2337 ----------
2338 eod_props: list of dict
2339 Properties of EODs as returned by `analyze_wave()` and
2340 `analyze_pulse()`. Only properties of pulse fish are saved.
2341 unit: string
2342 Unit of the waveform data.
2343 basename: string or stream
2344 If string, path and basename of file.
2345 If `basename` does not have an extension,
2346 '-pulsefish' and a file extension are appended.
2347 If stream, write pulse fish properties into this stream.
2348 kwargs:
2349 Arguments passed on to `TableData.write()`.
2351 Returns
2352 -------
2353 filename: Path or None
2354 Path and full name of the written file in case of `basename`
2355 being a string. Otherwise, the file name and extension that
2356 would have been appended to a basename.
2357 None if no pulse fish are contained in eod_props and
2358 consequently no file was written.
2360 See Also
2361 --------
2362 load_pulse_fish()
2363 """
2364 pulse_props = [p for p in eod_props if p['type'] == 'pulse']
2365 if len(pulse_props) == 0:
2366 return None
2367 td = TableData()
2368 if 'twin' in pulse_props[0] or 'samplerate' in pulse_props[0] or \
2369 'nfft' in pulse_props[0]:
2370 td.append_section('recording')
2371 if 'twin' in pulse_props[0]:
2372 td.append('twin', 's', '%7.2f', value=pulse_props)
2373 td.append('window', 's', '%7.2f', value=pulse_props)
2374 td.append('channel', '', '%d', value=pulse_props)
2375 td.append('winclipped', '%', '%.2f',
2376 value=pulse_props, fac=100)
2377 if 'samplerate' in pulse_props[0]:
2378 td.append('samplerate', 'kHz', '%.3f', value=pulse_props,
2379 fac=0.001)
2380 if 'nfft' in pulse_props[0]:
2381 td.append('nfft', '', '%d', value=pulse_props)
2382 td.append('dfreq', 'Hz', '%.2f', value=pulse_props)
2383 td.append_section('waveform')
2384 td.append('index', '', '%d', value=pulse_props)
2385 td.append('n', '', '%d', value=pulse_props)
2386 td.append('flipped', '', '%d', value=pulse_props)
2387 td.append('EODf', 'Hz', '%7.2f', value=pulse_props)
2388 td.append('period', 'ms', '%7.2f', value=pulse_props, fac=1000)
2389 td.append('aoffs', unit, '%.5g', value=pulse_props)
2390 td.append('ppampl', unit, '%.5g', value=pulse_props)
2391 td.append('relmaxampl', '%', '%.2f', value=pulse_props, fac=100)
2392 if 'noise' in pulse_props[0]:
2393 td.append('noise', '%', '%.2f', value=pulse_props, fac=100)
2394 td.append('rmserror', '%', '%.2f', value=pulse_props, fac=100)
2395 if 'clipped' in pulse_props[0]:
2396 td.append('clipped', '%', '%.2f', value=pulse_props, fac=100)
2397 td.append_section('timing')
2398 td.append('ppdist', 'ms', '%.3f', value=pulse_props, fac=1000)
2399 td.append('startendthresh', '%', '%.2f', value=pulse_props, fac=100)
2400 td.append('tstart', 'ms', '%.3f', value=pulse_props, fac=1000)
2401 td.append('tend', 'ms', '%.3f', value=pulse_props, fac=1000)
2402 td.append('width', 'ms', '%.3f', value=pulse_props, fac=1000)
2403 td.append_section('areas')
2404 td.append('totalarea', f'{unit}*ms', '%.4f', value=pulse_props, fac=1000)
2405 td.append('polaritybalance', '%', '%.2f', value=pulse_props, fac=100)
2406 td.append('median', 'ms', '%.3f', value=pulse_props, fac=1000)
2407 td.append('quartile1', 'ms', '%.3f', value=pulse_props, fac=1000)
2408 td.append('quartile3', 'ms', '%.3f', value=pulse_props, fac=1000)
2409 td.append('iqrange', 'ms', '%.3f', value=pulse_props, fac=1000)
2410 td.append_section('phases')
2411 td.append('phasethresh', '%', '%.2f', value=pulse_props, fac=100)
2412 td.append('nphases', '', '%d', value=pulse_props)
2413 td.append('firstphase', '', '%d', value=pulse_props)
2414 td.append('lastphase', '', '%d', value=pulse_props)
2415 td.append_section('tail')
2416 td.append('tau', 'ms', '%.3f', value=pulse_props, fac=1000)
2417 td.append('taustart', 'ms', '%.3f', value=pulse_props, fac=1000)
2418 td.append_section('spectrum')
2419 td.append('peakfreq', 'Hz', '%.2f', value=pulse_props)
2420 td.append('peakenergy', f'{unit}^2s/Hz', '%.3g', value=pulse_props)
2421 td.append('troughfreq', 'Hz', '%.2f', value=pulse_props)
2422 td.append('troughenergy', f'{unit}^2s/Hz', '%.3g', value=pulse_props)
2423 td.append('eodfenergy', f'{unit}^2s/Hz', '%.3g', value=pulse_props)
2424 td.append('energyatt5', 'dB', '%.2f', value=pulse_props)
2425 td.append('energyatt50', 'dB', '%.2f', value=pulse_props)
2426 td.append('lowcutoff', 'Hz', '%.2f', value=pulse_props)
2427 td.append('highcutoff', 'Hz', '%.2f', value=pulse_props)
2428 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
2429 fp = '-pulsefish' if not ext else ''
2430 return td.write_file_stream(basename, fp, **kwargs)
2433def load_pulse_fish(file_handle):
2434 """Load properties of pulse EODs from file.
2436 All times are scaled to seconds, all frequencies to Hertz, and all
2437 percentages to fractions.
2439 Parameters
2440 ----------
2441 file_handle: str, Path, or stream
2442 The file to be loaded.
2443 If str or Path and file does not exist, try to open
2444 possible corresponding zip file and read file from there.
2446 Returns
2447 -------
2448 eod_props: list of dict
2449 Properties of EODs.
2451 Raises
2452 ------
2453 FileNotFoundError:
2454 If `file_path` does not exist.
2456 See Also
2457 --------
2458 save_pulse_fish()
2460 """
2461 file_handle, zf = open_from_zip(file_handle)
2462 data = TableData(file_handle)
2463 close_zip(file_handle, zf)
2464 eod_props = data.dicts()
2465 for props in eod_props:
2466 if 'winclipped' in props:
2467 props['winclipped'] /= 100
2468 if 'samplerate' in props:
2469 props['samplerate'] *= 1000
2470 if 'nfft' in props:
2471 props['nfft'] = int(props['nfft'])
2472 props['type'] = 'pulse'
2473 props['index'] = int(props['index'])
2474 props['channel'] = int(props['channel'])
2475 props['n'] = int(props['n'])
2476 props['period'] /= 1000
2477 props['relmaxampl'] /= 100
2478 props['ppdist'] /= 1000
2479 props['rmserror'] /= 100
2480 props['noise'] /= 100
2481 if 'clipped' in props:
2482 props['clipped'] /= 100
2483 props['startendthresh'] /= 100
2484 props['phasethresh'] /= 100
2485 props['nphases'] = int(props['nphases'])
2486 props['tstart'] /= 1000
2487 props['tend'] /= 1000
2488 props['width'] /= 1000
2489 props['totalarea'] /= 1000
2490 props['polaritybalance'] /= 100
2491 props['median'] /= 1000
2492 props['quartile1'] /= 1000
2493 props['quartile3'] /= 1000
2494 props['iqrange'] /= 1000
2495 props['tau'] /= 1000
2496 props['taustart'] /= 1000
2497 props['firstphase'] = int(props['firstphase'])
2498 props['lastphase'] = int(props['lastphase'])
2499 return eod_props
2502def save_pulse_spectrum(spec_data, unit, idx, basename, **kwargs):
2503 """Save energy spectrum of pulse EOD to file.
2505 Parameters
2506 ----------
2507 spec_data: 2D array of floats
2508 Energy spectrum of single pulse as returned by `analyze_pulse()`.
2509 unit: string
2510 Unit of the waveform data.
2511 idx: int or None
2512 Index of fish.
2513 basename: string or stream
2514 If string, path and basename of file.
2515 If `basename` does not have an extension,
2516 '-pulsespectrum', the fish index, and a file extension are appended.
2517 If stream, write pulse spectrum into this stream.
2518 kwargs:
2519 Arguments passed on to `TableData.write()`.
2521 Returns
2522 -------
2523 filename: Path
2524 Path and full name of the written file in case of `basename`
2525 being a string. Otherwise, the file name and extension that
2526 would have been appended to a basename.
2528 See Also
2529 --------
2530 load_pulse_spectrum()
2531 """
2532 td = TableData(spec_data[:, :2], ['frequency', 'energy'],
2533 ['Hz', f'{unit}^2s/Hz'], ['%.2f', '%.4e'])
2534 if spec_data.shape[1] > 2:
2535 td.append('gaussian', f'{unit}^2s/Hz', '%.4e', value=spec_data[:, 2])
2536 fp = ''
2537 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
2538 if not ext:
2539 fp = '-pulsespectrum'
2540 if idx is not None:
2541 fp += f'-{idx}'
2542 return td.write_file_stream(basename, fp, **kwargs)
2545def load_pulse_spectrum(file_handle):
2546 """Load energy spectrum of pulse EOD from file.
2548 Parameters
2549 ----------
2550 file_handle: str, Path, or stream
2551 The file to be loaded.
2552 If str or Path and file does not exist, try to open
2553 possible corresponding zip file and read file from there.
2555 Returns
2556 -------
2557 spec: 2D array of floats
2558 Energy spectrum of single pulse: frequency, energy
2560 Raises
2561 ------
2562 FileNotFoundError:
2563 If `file_path` does not exist.
2565 See Also
2566 --------
2567 save_pulse_spectrum()
2568 """
2569 file_handle, zf = open_from_zip(file_handle)
2570 data = TableData(file_handle)
2571 close_zip(file_handle, zf)
2572 spec = data.array()
2573 return spec
2576def save_pulse_phases(phases, unit, idx, basename, **kwargs):
2577 """Save phase properties of pulse EOD to file.
2579 Parameters
2580 ----------
2581 phases: dict
2582 Dictionary with
2584 - "indices": indices of each phase
2585 (1 is P1, i.e. the largest positive peak)
2586 - "times": times of each phase relative to P1 in seconds
2587 - "amplitudes": amplitudes of each phase
2588 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
2589 - "widths": widths of each phase at `width_frac` height
2590 - "areas": areas of each phase
2591 - "relareas": areas of the phases relative to total area
2592 - "zeros": time of zero crossing towards next phase in seconds
2594 as returned by `analyze_pulse_phases()` and `analyze_pulse()`.
2595 unit: string
2596 Unit of the waveform data.
2597 idx: int or None
2598 Index of fish.
2599 basename: string or stream
2600 If string, path and basename of file.
2601 If `basename` does not have an extension,
2602 '-pulsephases', the fish index, and a file extension are appended.
2603 If stream, write pulse phases into this stream.
2604 kwargs:
2605 Arguments passed on to `TableData.write()`.
2607 Returns
2608 -------
2609 filename: Path
2610 Path and full name of the written file in case of `basename`
2611 being a string. Otherwise, the file name and extension that
2612 would have been appended to a basename.
2614 See Also
2615 --------
2616 load_pulse_phases()
2617 """
2618 if not phases or len(phases['times']) == 0:
2619 return None
2620 td = TableData()
2621 td.append('type', '', '%s', value=['P']*len(phases['indices']))
2622 td.append('index', '', '%.0f', value=phases['indices'])
2623 td.append('time', 'ms', '%.4f', value=phases['times'], fac=1000)
2624 td.append('amplitude', unit, '%.5f', value=phases['amplitudes'])
2625 td.append('relampl', '%', '%.2f', value=phases['relamplitudes'], fac=100)
2626 td.append('width', 'ms', '%.4f', value=phases['widths'], fac=1000)
2627 td.append('area', f'{unit}*ms', '%.4f', value=phases['areas'], fac=1000)
2628 td.append('relarea', '%', '%.2f', value=phases['relareas'], fac=100)
2629 td.append('zeros', 'ms', '%.4f', value=phases['zeros'], fac=1000)
2630 fp = ''
2631 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
2632 if not ext:
2633 fp = '-pulsephases'
2634 if idx is not None:
2635 fp += f'-{idx}'
2636 return td.write_file_stream(basename, fp, **kwargs)
2639def load_pulse_phases(file_handle):
2640 """Load phase properties of pulse EOD from file.
2642 Parameters
2643 ----------
2644 file_handle: str, Path, or stream
2645 The file to be loaded.
2646 If str or Path and file does not exist, try to open
2647 possible corresponding zip file and read file from there.
2649 Returns
2650 -------
2651 phases: dict
2652 Dictionary with
2654 - "indices": indices of each phase
2655 (1 is P1, i.e. the largest positive peak)
2656 - "times": times of each phase relative to P1 in seconds
2657 - "amplitudes": amplitudes of each phase
2658 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
2659 - "widths": widths of each phase at `width_frac` height
2660 - "areas": areas of each phase
2661 - "relareas": areas of the phases relative to total area
2662 - "zeros": time of zero crossing towards next phase in seconds
2664 unit: string
2665 Unit of phase amplitudes.
2667 Raises
2668 ------
2669 FileNotFoundError:
2670 If `file_path` does not exist.
2672 See Also
2673 --------
2674 save_pulse_phases()
2675 """
2676 file_handle, zf = open_from_zip(file_handle)
2677 data = TableData(file_handle)
2678 close_zip(file_handle, zf)
2679 phases = dict(indices=data['index'].astype(int),
2680 times=data['time']*0.001,
2681 amplitudes=data['amplitude'],
2682 relamplitudes=data['relampl']*0.01,
2683 widths=data['width']*0.001,
2684 areas=data['area']*0.001,
2685 relareas=data['relarea']*0.01,
2686 zeros=data['zeros']*0.001)
2687 return phases, data.unit('amplitude')
2690def save_pulse_gaussians(pulse, unit, idx, basename, **kwargs):
2691 """Save Gaussian phase properties of pulse EOD to file.
2693 Parameters
2694 ----------
2695 pulse: dict
2696 Dictionary with
2698 - "times": phase times in seconds,
2699 - "amplitudes": amplitudes, and
2700 - "stdevs": standard deviations in seconds
2702 of Gaussians fitted to the pulse waveform as returned by
2703 `decompose_pulse()` and `analyze_pulse()`.
2704 unit: string
2705 Unit of the waveform data.
2706 idx: int or None
2707 Index of fish.
2708 basename: string or stream
2709 If string, path and basename of file.
2710 If `basename` does not have an extension,
2711 '-pulsegaussians', the fish index, and a file extension are appended.
2712 If stream, write pulse phases into this stream.
2713 kwargs:
2714 Arguments passed on to `TableData.write()`.
2716 Returns
2717 -------
2718 filename: Path
2719 Path and full name of the written file in case of `basename`
2720 being a string. Otherwise, the file name and extension that
2721 would have been appended to a basename.
2723 See Also
2724 --------
2725 load_pulse_gaussians()
2727 """
2728 if len(pulse) == 0:
2729 return None
2730 td = TableData(pulse,
2731 units=['ms', unit, 'ms'],
2732 formats=['%.3f', '%.5f', '%.3f'])
2733 td['times'] *= 1000
2734 td['stdevs'] *= 1000
2735 fp = ''
2736 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
2737 if not ext:
2738 fp = '-pulsegaussians'
2739 if idx is not None:
2740 fp += f'-{idx}'
2741 return td.write_file_stream(basename, fp, **kwargs)
2744def load_pulse_gaussians(file_handle):
2745 """Load Gaussian phase properties of pulse EOD from file.
2747 Parameters
2748 ----------
2749 file_handle: str, Path, or stream
2750 The file to be loaded.
2751 If str or Path and file does not exist, try to open
2752 possible corresponding zip file and read file from there.
2754 Returns
2755 -------
2756 pulse: dict
2757 Dictionary with
2759 - "times": phase times in seconds,
2760 - "amplitudes": amplitudes, and
2761 - "stdevs": standard deviations in seconds
2763 of Gaussians fitted to the pulse waveform.
2764 Use the functions provided in thunderfish.fakefish to simulate
2765 pulse fish EODs from this data.
2766 unit: string
2767 Unit of Gaussian amplitudes.
2769 Raises
2770 ------
2771 FileNotFoundError:
2772 If `file_path` does not exist.
2774 See Also
2775 --------
2776 save_pulse_gaussians()
2778 """
2779 file_handle, zf = open_from_zip(file_handle)
2780 data = TableData(file_handle)
2781 close_zip(file_handle, zf)
2782 pulses = data.dict()
2783 pulses['times'] = 0.001*np.array(data['times'])
2784 pulses['amplitudes'] = np.array(data['amplitudes'])
2785 pulses['stdevs'] = 0.001*np.array(data['stdevs'])
2786 return pulses, data.unit('amplitudes')
2789def save_pulse_times(pulse_times, idx, basename, **kwargs):
2790 """Save times of pulse EOD to file.
2792 Parameters
2793 ----------
2794 pulse_times: dict or array of floats
2795 Times of EOD pulses. Either as array of times or
2796 `props['peaktimes']` or `props['times']` as returned by
2797 `analyze_pulse()`.
2798 idx: int or None
2799 Index of fish.
2800 basename: string or stream
2801 If string, path and basename of file.
2802 If `basename` does not have an extension,
2803 '-pulsetimes', the fish index, and a file extension are appended.
2804 If stream, write pulse times into this stream.
2805 kwargs:
2806 Arguments passed on to `TableData.write()`.
2808 Returns
2809 -------
2810 filename: Path
2811 Path and full name of the written file in case of `basename`
2812 being a string. Otherwise, the file name and extension that
2813 would have been appended to a basename.
2815 See Also
2816 --------
2817 load_pulse_times()
2818 """
2819 if isinstance(pulse_times, dict):
2820 props = pulse_times
2821 pulse_times = props.get('times', [])
2822 pulse_times = props.get('peaktimes', pulse_times)
2823 if len(pulse_times) == 0:
2824 return None
2825 td = TableData()
2826 td.append('time', 's', '%.4f', value=pulse_times)
2827 fp = ''
2828 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
2829 if not ext:
2830 fp = '-pulsetimes'
2831 if idx is not None:
2832 fp += f'-{idx}'
2833 return td.write_file_stream(basename, fp, **kwargs)
2836def load_pulse_times(file_handle):
2837 """Load times of pulse EOD from file.
2839 Parameters
2840 ----------
2841 file_handle: str, Path, or stream
2842 The file to be loaded.
2843 If str or Path and file does not exist, try to open
2844 possible corresponding zip file and read file from there.
2846 Returns
2847 -------
2848 pulse_times: array of floats
2849 Times of pulse EODs in seconds.
2851 Raises
2852 ------
2853 FileNotFoundError:
2854 If `file_path` does not exist.
2856 See Also
2857 --------
2858 save_pulse_times()
2859 """
2860 file_handle, zf = open_from_zip(file_handle)
2861 data = TableData(file_handle)
2862 close_zip(file_handle, zf)
2863 pulse_times = data.array()[:, 0]
2864 return pulse_times
2867def add_analyze_pulse_config(cfg, baseline_mode='mode', min_pulse_win=0.001,
2868 start_end_thresh_fac=0.01, peak_thresh_frac=0.002,
2869 min_dist=50.0e-6, width_frac=0.5, fit_frac=0.5,
2870 freq_resolution=1.0, fade_frac=0.0,
2871 flip_pulse='none', fit_gaussians=False,
2872 ipi_cv_thresh=0.5, ipi_percentile=30.0):
2873 """Add all parameters needed for the eod analysis functions as a new
2874 section to a configuration.
2876 Parameters
2877 ----------
2878 cfg: ConfigFile
2879 The configuration.
2881 See `analyze_pulse()` for details on the remaining arguments.
2882 """
2883 cfg.add_section('Pulse-type EOD analysis:')
2884 cfg.add('eodPulseBaseline', baseline_mode, '', 'Method to compute baseline voltage. One of "left". "right", "mean", or "mode".')
2885 cfg.add('eodMinPulseSnippet', min_pulse_win, 's', 'Minimum duration of cut out EOD snippets for a pulse fish.')
2886 cfg.add('eodStartEndThresholdFactor', start_end_thresh_fac, '', 'Threshold for for start and end time of pulse EODs as a fraction of the p-p amplitude.')
2887 cfg.add('eodPeakThresholdFactor', 100*peak_thresh_frac, '%', 'Threshold for detection of peaks and troughs in pulse EODs as a fraction of the p-p amplitude.')
2888 cfg.add('eodMinimumDistance', min_dist, 's', 'Minimum distance between peaks and troughs in a EOD pulse.')
2889 cfg.add('eodPulseWidthFraction', 100*width_frac, '%', 'The width of a pulse is measured at this fraction of the pulse height.')
2890 cfg.add('eodExponentialFitFraction', 100*fit_frac, '%', 'An exponential function is fitted on the tail of a pulse starting at this fraction of the height of the last peak.')
2891 cfg.add('eodPulseFrequencyResolution', freq_resolution, 'Hz', 'Frequency resolution of single pulse spectrum.')
2892 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.')
2893 cfg.add('flipPulseEOD', flip_pulse, '', 'Flip EOD of pulse fish to make the first large peak positive (flip, none, or auto).')
2894 cfg.add('eodFitGaussians', fit_gaussians, '', 'Fit sum of Gaussians to pulse-type EOD waveform.')
2895 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.')
2896 cfg.add('ipiPercentile', ipi_percentile, '%', 'Use only interpulse intervals shorter than this percentile to compute EOD frequency.')
2899def analyze_pulse_args(cfg):
2900 """Retrieve parameters for `analyze_pulse()` from configuration.
2902 The return value can then be passed as key-word arguments to this
2903 function.
2905 Parameters
2906 ----------
2907 cfg: ConfigFile
2908 The configuration.
2910 Returns
2911 -------
2912 a: dict
2913 Dictionary with names of arguments of the `analyze_pulse()` function
2914 and their values as supplied by `cfg`.
2915 """
2916 a = cfg.map(baseline_mode='eodPulseBaseline',
2917 min_pulse_win='eodMinPulseSnippet',
2918 start_end_thresh_fac='eodStartEndThresholdFactor',
2919 peak_thresh_frac='eodPeakThresholdFactor',
2920 min_dist='eodMinimumDistance',
2921 width_frac='eodPulseWidthFraction',
2922 fit_frac='eodExponentialFitFraction',
2923 freq_resolution='eodPulseFrequencyResolution',
2924 fade_frac='eodPulseFadeFraction',
2925 flip_pulse='flipPulseEOD',
2926 fit_gaussians='eodFitGaussians',
2927 ipi_cv_thresh='ipiCVThresh',
2928 ipi_percentile='ipiPercentile')
2929 a['peak_thresh_frac'] *= 0.01
2930 a['width_frac'] *= 0.01
2931 a['fit_frac'] *= 0.01
2932 a['fade_frac'] *= 0.01
2933 return a
2936def pulse_quality_args(cfg):
2937 """Retrieve parameters for `pulse_quality()` from configuration.
2939 The return value can then be passed as key-word arguments to this
2940 function.
2942 Parameters
2943 ----------
2944 cfg: ConfigFile
2945 The configuration.
2947 Returns
2948 -------
2949 a: dict
2950 Dictionary with names of arguments of the `pulse_quality()` function
2951 and their values as supplied by `cfg`.
2953 See Also
2954 --------
2955 eodanalysis.add_eod_quality_config()
2956 """
2957 a = cfg.map(max_clipped_frac='maximumClippedFraction',
2958 max_rms_sem='maximumVariance')
2959 a['max_clipped_frac'] *= 0.01
2960 return a
2963def main():
2964 import matplotlib.pyplot as plt
2965 from thunderlab.eventdetection import snippets
2966 from .fakefish import pulsefish_eods, export_pulsefish
2968 print('Analysis of pulse-type EOD waveforms.')
2970 deltaf = 2 # frequency resolution
2971 eodf = 80 # EOD frequency
2973 # data:
2974 rate = 100_000
2975 data = pulsefish_eods('Triphasic', eodf, rate, 5.0,
2976 jitter_cv=0.01, noise_std=0.01)
2977 unit = 'mV'
2978 eod_idx, _ = detect_peaks(data, 1.0)
2979 eod_times = eod_idx/rate
2981 # mean EOD:
2982 ns = int(0.003*rate)
2983 snips = snippets(data, eod_idx, start=-ns, stop=ns)
2984 mean_eod = np.mean(snips, 0)
2986 # analyse EOD:
2987 mean_eod, props, peaks, pulses, energy = \
2988 analyze_pulse(mean_eod, rate, eod_times,
2989 freq_resolution=deltaf, fit_gaussians=False)
2991 # write out as python code:
2992 export_pulsefish(pulses, '', 'Triphasic spec')
2994 # full spectrum:
2995 freq, power = pulsetrain_spectrum(eod_times, mean_eod, None,
2996 freq_resolution=deltaf, fade_frac=0.05)
2997 dfreq, dpower = psd(data, rate, freq_resolution=deltaf)
2999 # plot:
3000 fig, axs = plt.subplots(1, 2, layout='constrained')
3001 plot_pulse_eod(axs[0], mean_eod, props, peaks, unit=unit)
3002 axs[0].set_title(f'pulse fish: EODf = {props["EODf"]:.1f} Hz')
3003 ref = plot_pulse_spectrum(axs[1], energy, props)
3004 fac = 6300/deltaf # conversion power to energy spectrum (independent of recording duration and sampling rate, dependent on EOD frequency)
3005 print(f'FFT window={1/freq[1]:.2f}s, number of pulses={len(eod_times)}')
3006 axs[1].plot(freq, decibel(power, ref*fac), 'k')
3007 axs[1].plot(dfreq, decibel(dpower, ref*fac), 'gray')
3008 plt.show()
3011if __name__ == '__main__':
3012 main()