Coverage for src / thunderfish / pulseanalysis.py: 53%
813 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-02-09 14:25 +0000
1"""
2Analysis of pulse-type EOD waveforms.
4## Analysis of pulse-type EODs
6### Analysis of various pulse EOD aspects
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_spectrum()`: power spectrum of train of pulse-type EODs.
15- `analyze_pulse_spectrum()`: analyze the spectrum of a pulse-type EOD.
16- `analyze_pulse_intervals()`: basic statistics of interpulse intervals.
18### Complete analysis
20Calls all the functions listed above:
22- `analyze_pulse()`: analyze the EOD waveform of a pulse fish.
24## Visualization
26- `plot_pulse_eods()`: mark pulse EODs in a plot of an EOD recording.
27- `plot_pulse_spectrum()`: plot and annotate spectrum of single pulse EOD.
29## Storage
31- `save_pulse_fish()`: save properties of pulse EODs to file.
32- `load_pulse_fish()`: load properties of pulse EODs from file.
33- `save_pulse_spectrum()`: save spectrum of pulse EOD to file.
34- `load_pulse_spectrum()`: load spectrum of pulse EOD from file.
35- `save_pulse_phases()`: save phase properties of pulse EOD to file.
36- `load_pulse_phases()`: load phase properties of pulse EOD from file.
37- `save_pulse_gaussians()`: save Gaussian phase properties of pulse EOD to file.
38- `load_pulse_gaussians()`: load Gaussian phase properties of pulse EOD from file.
39- `save_pulse_times()`: save times of pulse EOD to file.
40- `load_pulse_times()`: load times of pulse EOD from file.
42## Fit functions
44- `gaussian_sum()`: sum of Gaussians.
45- `gaussian_sum_spectrum()`: energy spectrum of sum of Gaussians.
46- `gaussian_sum_costs()`: cost function for fitting sum of Gaussians.
47- `exp_decay()`: exponential decay.
49"""
51import numpy as np
53from pathlib import Path
54from scipy.optimize import curve_fit, minimize
55from thunderlab.eventdetection import detect_peaks, peak_width
56from thunderlab.powerspectrum import next_power_of_two, nfft, psd, decibel
57from thunderlab.tabledata import TableData
59from .fakefish import pulsefish_waveform, pulsefish_spectrum
62def condition_pulse(eod, ratetime=None, sem=None, flip='none',
63 baseline_frac=0.05, large_phase_frac=0.2,
64 min_pulse_win=0.001):
65 """Subtract offset, flip, shift, and cut out pulse EOD waveform.
67 Parameters
68 ----------
69 eod: 1-D or 2-D array of float
70 The eod waveform of which the spectrum is computed.
71 If an 1-D array, then this is the waveform and you
72 need to also pass a sampling rate in `rate`.
73 If a 2-D array, then first column is time in seconds and second
74 column is the eod waveform. If a third column is present,
75 it is interpreted as the standard error of the mean
76 corresponding to the averaged waveform.
77 Further columns are ignored.
78 ratetime: None or float or array of float
79 If a 1-D array is passed on to `eod` then either the sampling
80 rate in Hertz or the time array corresponding to `eod`.
81 sem: None or float or array of float
82 If a 1-D array is passed on to `eod`, then the optional
83 standard error of the averaged waveform. Either as a single value
84 or as an 1-D array for the whole waveform.
85 flip: 'auto', 'none', 'flip'
86 - 'auto' flip waveform such that the first large extremum is positive.
87 - 'flip' flip waveform.
88 - 'none' do not flip waveform.
89 baseline_frac: float
90 Fraction of data points from which the amplitude offset is computed.
91 large_phase_frac: float
92 Minimum amplitude of a large phase as a fraction of the largest one.
93 min_pulse_win: float
94 The minimum size of the cut-out EOD waveform.
96 Returns
97 -------
98 eod: 1-D or 2-D array of float
99 The input `eod` in the format of the input `eod`
100 with the conditioned waveform and optional time.
101 time: 1-D array of float
102 In case `eod` is a 1-D array, then the time array is also returned.
103 toffs: float
104 Time that was subtracted from the time axis,
105 such that the maximum of the EOD waveform was shifted to zero.
106 aoffs: float
107 Offset that was subtracted from the EOD waveform.
108 This is the average over the first `baseline_frac` data points
109 of the EOD waveform.
110 flipped: bool
111 True if waveform was flipped.
112 noise_thresh: float
113 Minimum threshold that is just higher than the noise level
114 within the first `baseline_frac` data points of the EOD waveform.
116 """
117 if eod.ndim == 2:
118 time = eod[:, 0]
119 meod = eod[:, 1]
120 if eod.shape[1] > 2:
121 sem = eod[:, 2]
122 else:
123 meod = eod
124 if isinstance(ratetime, (list, tuple, np.ndarray)):
125 time = ratetime
126 else:
127 time = np.arange(len(meod))/rate
128 if np.isscalar(sem):
129 sem = np.ones(len(meod))*sem
131 # subtract mean computed from the left end:
132 n_base = int(baseline_frac*len(meod))
133 if n_base < 5:
134 n_base = 5
135 aoffs = np.mean(meod[:n_base]) # baseline
136 meod -= aoffs
138 # flip waveform:
139 max_idx = np.argmax(meod)
140 min_idx = np.argmin(meod)
141 flipped = 'flip' in flip.lower()
142 if 'auto' in flip.lower():
143 max_ampl = abs(meod[max_idx])
144 min_ampl = abs(meod[min_idx])
145 amplitude = max(max_ampl, min_ampl)
146 if max_ampl > large_phase_frac*amplitude and \
147 min_ampl > large_phase_frac*amplitude:
148 # two major peaks:
149 if min_idx < max_idx:
150 flipped = True
151 elif min_ampl > large_phase_frac*amplitude:
152 flipped = True
153 if flipped:
154 meod = -meod
155 idx = min_idx
156 min_idx = max_idx
157 max_idx = idx
158 max_ampl = abs(meod[max_idx])
159 min_ampl = abs(meod[min_idx])
161 # shift maximum to zero:
162 toffs = time[max_idx]
163 time -= toffs
165 # threshold from baseline maximum and minimum:
166 th_max = np.max(meod[:n_base])
167 th_min = np.min(meod[:n_base])
168 range_thresh = 2*(th_max - th_min)
169 # threshold from standard error:
170 if sem is not None:
171 msem = np.mean(sem[:n_base])
172 sem_thresh = 8*msem
173 # noise threshold:
174 noise_thresh = max(range_thresh, sem_thresh)
176 # generous left edge of waveform:
177 l1_idx = np.argmax(np.abs(meod) > noise_thresh)
178 l2_idx = np.argmax(np.abs(meod) > 2*noise_thresh)
179 w = 2*(l2_idx - l1_idx)
180 if w < n_base:
181 w = n_base
182 l_idx = l1_idx - w
183 if l_idx < 0:
184 l_idx = 0
185 # generous right edge of waveform:
186 r1_idx = len(meod) - 1 - np.argmax(np.abs(meod[::-1]) > noise_thresh)
187 r2_idx = len(meod) - 1 - np.argmax(np.abs(meod[::-1]) > 2*noise_thresh)
188 w = 2*(r1_idx - r2_idx)
189 if w < n_base:
190 w = n_base
191 r_idx = r1_idx + w
192 if r_idx >= len(meod):
193 r_idx = len(meod)
194 # cut out relevant signal:
195 if time[r_idx - 1] - time[l_idx] < min_pulse_win:
196 ct = time[(l_idx + r_idx)//2]
197 mask = (time >= ct - min_pulse_win/2) & (time <= ct + min_pulse_win/2)
198 meod = meod[mask]
199 time = time[mask]
200 if eod.ndim == 2:
201 eod = eod[mask, :]
202 else:
203 meod = meod[l_idx:r_idx]
204 time = time[l_idx:r_idx]
205 if eod.ndim == 2:
206 eod = eod[l_idx:r_idx, :]
208 # return offset, flipped and, shifted waveform:
209 if eod.ndim == 2:
210 eod[:, 0] = time
211 eod[:, 1] = meod
212 return eod, toffs, aoffs, flipped, noise_thresh
213 else:
214 return meod, time, toffs, aoffs, flipped, noise_thresh
217def analyze_pulse_properties(noise_thresh, eod, ratetime=None):
218 """Characterize basic properties of a pulse-type EOD.
220 Parameters
221 ----------
222 noise_thresh: float
223 Minimum threshold that is just higher than the baseline noise level.
224 As returned by `condition_pulse()`.
225 eod: 1-D or 2-D array
226 The eod waveform to be analyzed.
227 If an 1-D array, then this is the waveform and you
228 need to also pass a sampling rate in `rate`.
229 If a 2-D array, then first column is time in seconds and second
230 column is the eod waveform. Further columns are ignored.
231 ratetime: None or float or array of float
232 If a 1-D array is passed on to `eod` then either the sampling
233 rate in Hertz or the time array corresponding to `eod`.
235 Returns
236 -------
237 pos_ampl: float
238 Amplitude of largest positive peak.
239 neg_ampl: float
240 Amplitude of largest negative trough (absolute value).
241 dist: float
242 Temporal distance between largest negative trough and positive peak.
243 pos_area: float
244 Integral under all positive values of EOD waveform.
245 neg_area: float
246 Integral under all negative values of EOD waveform (absolute value).
247 polarity_balance: float
248 Contrast between positive and negative areas of EOD waveform, i.e.
249 (pos_area - neg_area)/(pos_area + neg_area).
250 center: float
251 Center of mass (first moment when treating the absolute value of
252 the waveform as a distribution).
253 stdev: float
254 Standard deviation of mass (square root of second central moment
255 when treating the absolute value of the waveform as a distribution).
256 median: float
257 Median of the distribution of the absolute EOD waveform.
258 quartile1: float
259 First quartile of the distribution of the absolute EOD waveform.
260 quartile3: float
261 Third quartile of the distribution of the absolute EOD waveform.
262 """
263 if eod.ndim == 2:
264 time = eod[:, 0]
265 eod = eod[:, 1]
266 elif isinstance(ratetime, (list, tuple, np.ndarray)):
267 time = ratetime
268 else:
269 time = np.arange(len(eod))/rate
270 dt = time[1] - time[0]
272 # amplitudes:
273 pos_idx = np.argmax(eod)
274 pos_ampl = abs(eod[pos_idx])
275 if pos_ampl < noise_thresh:
276 pos_ampl = 0
277 neg_idx = np.argmin(eod)
278 neg_ampl = abs(eod[neg_idx])
279 if neg_ampl < noise_thresh:
280 neg_ampl = 0
281 dist = time[neg_idx] - time[pos_idx]
282 if pos_ampl < noise_thresh or neg_ampl < noise_thresh:
283 dist = np.inf
285 # integrals (areas) and polarity balance:
286 pos_area = abs(np.sum(eod[eod >= 0]))*dt
287 neg_area = abs(np.sum(eod[eod < 0]))*dt
288 total_area = np.sum(np.abs(eod))*dt
289 integral_area = np.sum(eod)*dt # = pos_area - neg_area
290 polarity_balance = integral_area/total_area
292 # moments (EOD waveforms are not Gaussian!):
293 #center = np.sum(time*np.abs(eod))/np.sum(np.abs(eod))
294 #var = np.sum((time - center)**2*np.abs(eod))/np.sum(np.abs(eod))
295 #stdev = np.sqrt(var)
297 # cumulative:
298 cumul = np.cumsum(np.abs(eod))/np.sum(np.abs(eod))
299 median = time[np.argmax(cumul > 0.5)]
300 quartile1 = time[np.argmax(cumul > 0.25)]
301 quartile3 = time[np.argmax(cumul > 0.75)]
303 return pos_ampl, neg_ampl, dist, \
304 pos_area, neg_area, polarity_balance, \
305 median, quartile1, quartile3
308def analyze_pulse_phases(peak_thresh, startend_thresh,
309 eod, ratetime=None,
310 min_dist=50.0e-6, width_frac=0.5):
311 """Characterize all phases of a pulse-type EOD.
313 Parameters
314 ----------
315 peak_thresh: float
316 Threshold for detecting peaks and troughs.
317 startend_thresh: float or None
318 Threshold for detecting start and end time of EOD.
319 If None, use `peak_thresh`.
320 eod: 1-D or 2-D array
321 The eod waveform to be analyzed.
322 If an 1-D array, then this is the waveform and you
323 need to also pass a sampling rate in `rate`.
324 If a 2-D array, then first column is time in seconds and second
325 column is the eod waveform. Further columns are ignored.
326 ratetime: None or float or array of float
327 If a 1-D array is passed on to `eod` then either the sampling
328 rate in Hertz or the time array corresponding to `eod`.
329 min_dist: float
330 Minimum distance between peak and troughs of the pulse.
331 width_frac: float
332 The width of an EOD phase is measured at this fraction of a peak's
333 height (0-1).
335 Returns
336 -------
337 tstart: float
338 Start time of EOD waveform, i.e. the first time it crosses `threshold`.
339 tend: float
340 End time of EOD waveform, i.e. the last time it falls under `threshold`.
341 phases: dict
342 Dictionary with
344 - "indices": indices of each phase
345 (1 is P1, i.e. the largest positive peak)
346 - "times": times of each phase relative to P1 in seconds
347 - "amplitudes": amplitudes of each phase
348 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
349 - "widths": widths of each phase at `width_frac` height
350 - "areas": areas of each phase
351 - "relareas": areas of the phases relative to total area
352 - "zeros": time of zero crossing towards next phase in seconds
354 Empty dictionary if waveform is not a pulse EOD.
356 """
357 if eod.ndim == 2:
358 time = eod[:, 0]
359 eod = eod[:, 1]
360 elif isinstance(ratetime, (list, tuple, np.ndarray)):
361 time = ratetime
362 else:
363 time = np.arange(len(eod))/rate
364 dt = time[1] - time[0]
366 # start and end time:
367 if startend_thresh is None:
368 startend_thresh = peak_thresh
369 l_idx = np.argmax(np.abs(eod) > startend_thresh)
370 r_idx = len(eod) - 1 - np.argmax(np.abs(eod[::-1]) > startend_thresh)
371 tstart = time[l_idx]
372 tend = time[r_idx]
374 # find peaks and troughs:
375 peak_idx, trough_idx = detect_peaks(eod, peak_thresh)
376 if len(peak_idx) == 0:
377 return tstart, tend, {}
378 # and their width:
379 peak_widths = peak_width(time, eod, peak_idx, trough_idx,
380 peak_frac=width_frac, base='max')
381 trough_widths = peak_width(time, -eod, trough_idx, peak_idx,
382 peak_frac=width_frac, base='max')
383 # combine peaks and troughs:
384 pt_idx = np.concatenate((peak_idx, trough_idx))
385 pt_widths = np.concatenate((peak_widths, trough_widths))
386 pts_idx = np.argsort(pt_idx)
387 phase_list = pt_idx[pts_idx]
388 width_list = pt_widths[pts_idx]
389 # remove phases at the start and end of the signal:
390 n = len(eod)//20
391 if n < 5:
392 n = 5
393 mask = (phase_list > n) & (phase_list < len(eod) - 20)
394 phase_list = phase_list[mask]
395 width_list = width_list[mask]
396 if len(phase_list) == 0:
397 return tstart, tend, {}
398 # remove multiple peaks that are too close:
399 # TODO: XXX replace by Dexters function that keeps the maximum peak
400 rmidx = [(k, k+1) for k in np.where(np.diff(time[phase_list]) < min_dist)[0]]
401 # flatten and keep maximum and minimum phase:
402 max_idx = np.argmax(eod)
403 min_idx = np.argmin(eod)
404 rmidx = np.unique([k for kk in rmidx for k in kk
405 if phase_list[k] != max_idx and phase_list[k] != min_idx])
406 # delete:
407 if len(rmidx) > 0:
408 phase_list = np.delete(phase_list, rmidx)
409 width_list = np.delete(width_list, rmidx)
410 if len(phase_list) == 0:
411 return tstart, tend, {}
412 # find P1:
413 p1i = np.argmax(phase_list == max_idx)
414 # amplitudes:
415 amplitudes = eod[phase_list]
416 max_ampl = np.abs(amplitudes[p1i])
417 # phase areas and zeros:
418 phase_areas = np.zeros(len(phase_list))
419 zero_times = np.zeros(len(phase_list))
420 for i in range(len(phase_list)):
421 sign_fac = np.sign(eod[phase_list[i]])
422 i0 = phase_list[i - 1] if i > 0 else 0
423 i1 = phase_list[i + 1] if i + 1 < len(phase_list) else len(eod)
424 if i0 > 0 and sign_fac*eod[i0] > 0 and \
425 i1 < len(eod) and sign_fac*eod[i1] > 0:
426 phase_areas[i] = 0
427 else:
428 snippet = sign_fac*eod[i0:i1]
429 phase_areas[i] = np.sum(snippet[snippet > 0])*dt
430 phase_areas[i] *= sign_fac
431 if i < len(phase_list) - 1:
432 i0 = phase_list[i]
433 snippet = eod[i0:i1]
434 stimes = time[i0:i1]
435 zidx = np.nonzero(snippet[:-1]*snippet[1:] < 0)[0]
436 if len(zidx) == 0:
437 zero_times[i] = np.nan
438 else:
439 zidx = zidx[len(zidx)//2] # reduce to single zero crossing
440 snippet = snippet[zidx:zidx + 2]
441 stimes = stimes[zidx:zidx + 2]
442 if sign_fac > 0:
443 zero_times[i] = np.interp(0, snippet[::-1], stimes[::-1])
444 else:
445 zero_times[i] = np.interp(0, snippet, stimes)
446 else:
447 zero_times[i] = np.nan
448 total_area = np.sum(np.abs(phase_areas))
449 # store phase properties:
450 phases = dict(indices=np.arange(len(phase_list)) + 1 - p1i,
451 times=time[phase_list],
452 amplitudes=amplitudes,
453 relamplitudes=amplitudes/max_ampl,
454 widths=width_list,
455 areas=phase_areas,
456 relareas=phase_areas/total_area,
457 zeros=zero_times)
458 return tstart, tend, phases
461def gaussian_sum(x, *pas):
462 """Sum of Gaussians.
464 Parameters
465 ----------
466 x: array of float
467 The x array over which the sum of Gaussians is evaluated.
468 *pas: list of floats
469 The parameters of the Gaussians in a flat list.
470 Position, amplitude, and standard deviation of first Gaussian,
471 position, amplitude, and standard deviation of second Gaussian,
472 and so on.
474 Returns
475 -------
476 sg: array of float
477 The sum of Gaussians for the times given in `t`.
478 """
479 sg = np.zeros(len(x))
480 for pos, ampl, std in zip(pas[0:-2:3], pas[1:-1:3], pas[2::3]):
481 sg += ampl*np.exp(-0.5*((x - pos)/std)**2)
482 return sg
485def gaussian_sum_spectrum(f, *pas):
486 """Energy spectrum of sum of Gaussians.
488 Parameters
489 ----------
490 f: 1-D array of float
491 The frequencies at which to evaluate the spectrum.
492 *pas: list of floats
493 The parameters of the Gaussians in a flat list.
494 Position, amplitude, and standard deviation of first Gaussian,
495 position, amplitude, and standard deviation of second Gaussian,
496 and so on.
498 Returns
499 -------
500 energy: 1-D array of float
501 The one-sided energy spectrum of the sum of Gaussians.
502 """
503 spec = np.zeros(len(f), dtype=complex)
504 for dt, a, s in zip(pas[0:-2:3], pas[1:-1:3], pas[2::3]):
505 gauss = a*np.sqrt(2*np.pi)*s*np.exp(-0.5*(2*np.pi*s*f)**2)
506 shift = np.exp(-2j*np.pi*f*dt)
507 spec += gauss*shift
508 spec *= np.sqrt(2) # because of one-sided spectrum
509 return np.abs(spec)**2
512def gaussian_sum_costs(pas, time, eod, freqs, energy):
513 """ Cost function for fitting sum of Gaussian to pulse EOD.
515 Parameters
516 ----------
517 pas: list of floats
518 The pulse parameters in a flat list.
519 Position, amplitude, and standard deviation of first phase,
520 position, amplitude, and standard deviation of second phase,
521 and so on.
522 time: 1-D array of float
523 Time points of the EOD waveform.
524 eod: 1-D array of float
525 The real EOD waveform.
526 freqs: 1-D array of float
527 The frequency components of the spectrum.
528 energy: 1-D array of float
529 The energy spectrum of the real pulse.
531 Returns
532 -------
533 costs: float
534 Weighted sum of rms waveform difference and rms spectral difference.
535 """
536 eod_fit = gaussian_sum(time, *pas)
537 eod_rms = np.sqrt(np.mean((eod_fit - eod)**2))/np.max(np.abs(eod))
538 level = decibel(energy)
539 level_range = 30
540 n = np.argmax(level)
541 energy_fit = gaussian_sum_spectrum(freqs, *pas)
542 level_fit = decibel(energy_fit)
543 weight = np.ones(n)
544 weight[freqs[:n] < 10] = 100
545 weight /= np.sum(weight)
546 spec_rms = np.sqrt(np.mean(weight*(level_fit[:n] - level[:n])**2))/level_range
547 costs = eod_rms + 5*spec_rms
548 #print(f'{costs:.4f}, {eod_rms:.4f}, {spec_rms:.4f}')
549 return costs
552def decompose_pulse(eod, freqs, energy, phases, width_frac=0.5, verbose=0):
553 """Decompose single pulse waveform into sum of Gaussians.
555 Use the output to simulate pulse-type EODs using the functions
556 provided in the thunderfish.fakefish module.
558 Parameters
559 ----------
560 eod: 2-D array of float
561 The eod waveform. First column is time in seconds and second
562 column is the eod waveform. Further columns are ignored.
563 freqs: 1-D array of float
564 The frequency components of the spectrum.
565 energy: 1-D array of float
566 The energy spectrum of the real pulse.
567 phases: dict
568 Properties of the EOD phases as returned by analyze_pulse_phases().
569 width_frac: float
570 The width of a peak is measured at this fraction of a peak's
571 height (0-1).
572 verbose: int
573 Verbosity level passed for error and info messages.
575 Returns
576 -------
577 pulse: dict
578 Dictionary with
580 - "times": phase times in seconds,
581 - "amplitudes": amplitudes, and
582 - "stdevs": standard deviations in seconds
584 of Gaussians fitted to the pulse waveform. Use the functions
585 provided in thunderfish.fakefish to simulate pulse fish EODs
586 from this data.
588 """
589 pulse = {}
590 if len(phases) == 0:
591 return pulse
592 if eod.ndim == 2:
593 time = eod[:, 0]
594 eod = eod[:, 1]
595 elif isinstance(ratetime, (list, tuple, np.ndarray)):
596 time = ratetime
597 else:
598 time = np.arange(len(eod))/rate
599 # convert half width to standard deviation:
600 fac = 0.5/np.sqrt(-2*np.log(width_frac))
601 # fit parameter as single list:
602 tas = []
603 for t, a, s in zip(phases['times'], phases['amplitudes'],
604 phases['widths']*fac):
605 tas.extend((t, a, s))
606 tas = np.asarray(tas)
607 # fit EOD waveform:
608 try:
609 tas, _ = curve_fit(gaussian_sum, time, eod, tas)
610 except RuntimeError as e:
611 if verbose > 0:
612 print('Fit gaussian_sum failed in decompose_pulse():', e)
613 return pulse
614 # fit EOD waveform and spectrum:
615 bnds = [(1e-5, None) if k%3 == 2 else (None, None)
616 for k in range(len(tas))]
617 res = minimize(gaussian_sum_costs, tas,
618 args=(time, eod, freqs, energy), bounds=bnds)
619 if not res.success and verbose > 0:
620 print('warning: optimization gaussian_sum_costs failed in decompose_pulse():',
621 res.message)
622 else:
623 tas = res.x
624 # add another Gaussian:
625 rms_norm = np.max(np.abs(eod))
626 rms_old = np.sqrt(np.mean((gaussian_sum(time, *tas) - eod)**2))/rms_norm
627 eod_diff = np.abs(gaussian_sum(time, *tas) - eod)/rms_norm
628 if np.max(eod_diff) > 0.1:
629 if verbose > 1:
630 print(f'decompose_pulse(): added Gaussian because maximum rms error was {100*np.max(eod_diff):.0f}%')
631 ntas = np.concatenate((tas, (time[np.argmax(eod_diff)], np.max(eod_diff),
632 np.mean(tas[2::3]))))
633 bnds = [(1e-5, None) if k%3 == 2 else (None, None)
634 for k in range(len(ntas))]
635 res = minimize(gaussian_sum_costs, ntas,
636 args=(time, eod, freqs, energy), bounds=bnds)
637 if res.success:
638 rms_new = np.sqrt(np.mean((gaussian_sum(time, *res.x) - eod)**2))/rms_norm
639 if rms_new < 0.8*rms_old:
640 tas = res.x
641 elif verbose > 1:
642 print('decompose_pulse(): removed added Gaussian because it did not improve the fit')
643 elif verbose > 0:
644 print('warnong: optimization gaussian_sum_costs for additional Gaussian failed in decompose_pulse():',
645 res.message)
646 times = np.asarray(tas[0::3])
647 ampls = np.asarray(tas[1::3])
648 stdevs = np.asarray(tas[2::3])
649 pulse = dict(times=times, amplitudes=ampls, stdevs=stdevs)
650 return pulse
653def exp_decay(t, tau, ampl, offs):
654 """Exponential decay function.
656 x(t) = ampl*exp(-t/tau) + offs
658 Parameters
659 ----------
660 t: float or array
661 Time.
662 tau: float
663 Time constant of exponential decay.
664 ampl: float
665 Amplitude of exponential decay, i.e. initial value minus
666 steady-state value.
667 offs: float
668 Steady-state value.
670 Returns
671 -------
672 x: float or array
673 The exponential decay evaluated at times `t`.
675 """
676 return offs + ampl*np.exp(-t/tau)
679def analyze_pulse_tail(peak_index, eod, ratetime=None,
680 threshold=0.0, fit_frac=0.5, verbose=0):
681 """ Fit exponential to last peak/trough of pulse EOD.
683 Parameters
684 ----------
685 peak_index: int
686 Index of last peak in `eod`.
687 eod: 1-D or 2-D array
688 The eod waveform to be analyzed.
689 If an 1-D array, then this is the waveform and you
690 need to also pass a sampling rate in `rate`.
691 If a 2-D array, then first column is time in seconds and second
692 column is the eod waveform. Further columns are ignored.
693 ratetime: None or float or array of float
694 If a 1-D array is passed on to `eod` then either the sampling
695 rate in Hertz or the time array corresponding to `eod`.
696 threshold: float
697 Maximum noise level of the pulse waveform.
698 fit_frac: float or None
699 An exponential is fitted to the tail of the last peak/trough
700 starting where the waveform falls below this fraction of the
701 peak's height (0-1).
702 If None, do not attempt to fit.
703 verbose: int
704 Verbosity level passed for error and info messages.
706 Returns
707 -------
708 tau: float or np.nan
709 Time constant of pulse tail in seconds.
710 tstart: float or np.nan
711 Time where fit started in seconds.
712 fit: 1-D array of float or np.nan
713 Time trace of the fit corresponding to `eod`.
714 """
715 if fit_frac is None:
716 return np.nan, np.nan, None
717 if eod.ndim == 2:
718 time = eod[:, 0]
719 eod = eod[:, 1]
720 elif isinstance(ratetime, (list, tuple, np.ndarray)):
721 time = ratetime
722 else:
723 time = np.arange(len(eod))/rate
724 dt = np.mean(np.diff(time))
725 pi = peak_index
726 # positive or negative decay:
727 sign = 1
728 eodpp = eod[pi:] - 0.5*threshold
729 eodpn = -eod[pi:] - 0.5*threshold
730 if np.sum(eodpn[eodpn > 0]) > np.sum(eodpp[eodpp > 0]):
731 sign = -1
732 if sign*eod[pi] < 0:
733 pi += np.argmax(sign*eod[pi:])
734 pi_ampl = np.abs(eod[pi])
735 # no sufficiently large initial value:
736 sampl = sign*eod[pi]*fit_frac
737 if sampl <= threshold:
738 if verbose > 0:
739 print(f'exponential fit to tail of pulse failed: initial amplitude {sampl:.5f} smaller than threshold {threshold:.5f}')
740 return np.nan, np.nan, None
741 # no sufficiently long decay:
742 sinx = pi + np.argmax(sign*eod[pi:] < sampl)
743 n = 2*len(eod[sinx:])//3
744 if n < 10:
745 if verbose > 0:
746 print(f'exponential fit to tail of pulse failed: less than 10 samples {n}')
747 return np.nan, np.nan, None
748 # not decaying towards zero:
749 max_line = sampl - (sampl - threshold)*np.arange(n)/n + 1e-8
750 above_frac = np.sum(sign*eod[sinx:sinx + n] > max_line)/n
751 if above_frac > 0.05:
752 if verbose > 0:
753 print(f'exponential fit to tail of pulse failed: not decaying towards zero {100*above_frac:.1f}% > 5%')
754 return np.nan, np.nan, None
755 # estimate tau:
756 thresh = eod[sinx]*np.exp(-1.0)
757 tau_inx = np.argmax(sign*eod[sinx:] < sign*thresh)
758 if tau_inx < 2:
759 tau_inx = 2
760 rinx = sinx + 6*tau_inx
761 if rinx > len(eod) - 1:
762 rinx = len(eod) - 1
763 if rinx - sinx < 2*tau_inx:
764 if verbose > 0:
765 print(f'exponential fit to tail of pulse failed: less samples {rinx - sinx} than two time constants 3*{tau_inx}')
766 return np.nan, np.nan, None
767 tau = time[sinx + tau_inx] - time[sinx]
768 params = [tau, eod[sinx] - eod[rinx], eod[rinx]]
769 try:
770 popt, pcov = curve_fit(exp_decay, time[sinx:rinx] - time[sinx],
771 eod[sinx:rinx], params,
772 bounds=([0.0, -np.inf, -np.inf], np.inf))
773 except TypeError:
774 popt, pcov = curve_fit(exp_decay, time[sinx:rinx] - time[sinx],
775 eod[sinx:rinx], params)
776 if popt[0] > 1.2*tau:
777 tau_inx = int(np.round(popt[0]/dt))
778 rinx = sinx + 6*tau_inx
779 if rinx > len(eod) - 1:
780 rinx = len(eod) - 1
781 try:
782 popt, pcov = curve_fit(exp_decay, time[sinx:rinx] - time[sinx],
783 eod[sinx:rinx], popt,
784 bounds=([0.0, -np.inf, -np.inf], np.inf))
785 except TypeError:
786 popt, pcov = curve_fit(exp_decay, time[sinx:rinx] - time[sinx],
787 eod[sinx:rinx], popt)
788 tau = popt[0]
789 tstart = time[sinx]
790 fit = np.zeros(len(eod))
791 fit[:] = np.nan
792 fit[sinx:rinx] = exp_decay(time[sinx:rinx] - time[sinx], *popt)
793 if verbose > 0:
794 print(f'exponential fit to tail of pulse: got time constant {1000*tau:.3f}ms')
795 return tau, tstart, fit
798def pulse_spectrum(eod, ratetime=None, freq_resolution=1.0, fade_frac=0.0):
799 """Spectrum of a single pulse-type EOD.
801 Parameters
802 ----------
803 eod: 1-D or 2-D array
804 The eod waveform of which the spectrum is computed.
805 If an 1-D array, then this is the waveform and you
806 need to also pass a sampling rate in `rate`.
807 If a 2-D array, then first column is time in seconds and second
808 column is the eod waveform. Further columns are ignored.
809 ratetime: None or float or array of float
810 If a 1-D array is passed on to `eod` then either the sampling
811 rate in Hertz or the time array corresponding to `eod`.
812 freq_resolution: float
813 The frequency resolution of the spectrum.
814 fade_frac: float
815 Fraction of time of the EOD waveform that is used to fade in
816 and out to zero baseline.
818 Returns
819 -------
820 freqs: 1-D array of float
821 The frequency components of the energy spectrum.
822 energy: 1-D array of float
823 The energy spectrum of the single pulse EOD
824 with unit (x s)^2 = x^2 s/Hz.
825 The integral over the energy spectrum `np.sum(energy)*freqs[1]`
826 equals the integral over the squared eod, `np.sum(eod**2)/rate`.
827 That is, by making the energy spectrum a power sepctrum
828 (dividing the energy by the FFT window duration), the integral
829 over the power spectrum equals the mean-squared signal
830 (variance). But the single-pulse spectrum is not a power-spectrum.
831 because in the limit to infinitely long window, the power vanishes!
833 See Also
834 --------
835 thunderfish.fakefish.pulsefish_spectrum(): analytically computed spectra for simulated pulse EODs.
836 pulsetrain_spectrum(): power spectrum of train of pulse-type EODs.
838 """
839 if eod.ndim == 2:
840 rate = 1.0/(eod[1, 0] - eod[0, 0])
841 eod = eod[:, 1]
842 elif isinstance(ratetime, (list, tuple, np.ndarray)):
843 rate = 1.0/(ratetime[1] - ratetime[0])
844 else:
845 rate = ratetime
846 n_fft = nfft(rate, freq_resolution)
847 # subtract mean computed from the ends of the EOD snippet:
848 n = len(eod)//20 if len(eod) >= 20 else 1
849 eod = eod - 0.5*(np.mean(eod[:n]) + np.mean(eod[-n:]))
850 # zero padding:
851 max_idx = np.argmax(eod)
852 n0 = max_idx
853 n1 = len(eod) - max_idx
854 n = 2*max(n0, n1)
855 if n_fft < n:
856 n_fft = next_power_of_two(n)
857 data = np.zeros(n_fft)
858 data[n_fft//2 - n0:n_fft//2 + n1] = eod
859 # fade in and out:
860 if fade_frac > 0:
861 fn = int(fade_frac*len(eod))
862 data[n_fft//2 - n0:n_fft//2 - n0 + fn] *= np.arange(fn)/fn
863 data[n_fft//2 + n1 - fn:n_fft//2 + n1] *= np.arange(fn)[::-1]/fn
864 # spectrum:
865 dt = 1/rate
866 freqs = np.fft.rfftfreq(n_fft, dt)
867 fourier = np.fft.rfft(data)*dt
868 energy = 2*np.abs(fourier)**2 # one-sided spectrum!
869 return freqs, energy
872def pulsetrain_spectrum(eod_times, eod, ratetime=None,
873 freq_resolution=1.0, fade_frac=0.0):
874 """Power spectrum of train of pulse-type EODs.
876 Places mean eod waveform at `eod_times` and computes from that the
877 power spectrum.
879 Parameters
880 ----------
881 eod_times: 1-D array or None
882 List of times of detected EODs.
883 eod: 1-D or 2-D array
884 The eod waveform of which the spectrum is computed.
885 If an 1-D array, then this is the waveform and you
886 need to also pass a sampling rate in `rate`.
887 If a 2-D array, then first column is time in seconds and second
888 column is the eod waveform. Further columns are ignored.
889 ratetime: None or float or array of float
890 If a 1-D array is passed on to `eod` then either the sampling
891 rate in Hertz or the time array corresponding to `eod`.
892 freq_resolution: float
893 The frequency resolution of the spectrum.
894 fade_frac: float
895 Fraction of time of the EOD waveform that is used to fade in
896 and out to zero baseline.
898 Returns
899 -------
900 freq: 1-D array of float
901 Frequencies corresponding to power array.
902 power: 1-D array of float
903 Power spectral density in [eod]^2/Hz.
905 """
906 if eod.ndim == 2:
907 rate = 1.0/(eod[1, 0] - eod[0, 0])
908 time = eod[:, 0]
909 eod = eod[:, 1]
910 elif isinstance(ratetime, (list, tuple, np.ndarray)):
911 time = ratetime
912 rate = 1.0/(time[1] - time[0])
913 else:
914 time = np.arange(len(eod))/ratetime
916 ipi = np.median(np.diff(eod_times))
917 n_ipi = int(ipi*rate)
918 n = int(eod_times[-1]*rate) + n_ipi//2
919 data = np.zeros(n)
920 i0 = np.argmin(np.abs(time))
921 i1 = len(eod) - i0
922 fn = int(fade_frac*len(eod))
923 # place eod waveforms at eod_times:
924 prev_idx = -n_ipi
925 for t in eod_times:
926 idx = int(np.round(t*rate))
927 ii0 = i0 if idx - i0 >= 0 else idx
928 if idx - ii0 < prev_idx + n_ipi//2:
929 ii0 = idx - prev_idx - n_ipi//2
930 if ii0 < 0:
931 ii0 = 0
932 ii1 = i1 if idx + i1 < len(data) else len(data) - 1 - idx
933 data[idx - ii0:idx + ii1] = eod[i0 - ii0:i0 + ii1]
934 # fade in and out:
935 if fn > 0:
936 data[idx - ii0:idx - ii0 + fn] *= np.arange(fn)/fn
937 data[idx + ii1 - fn:idx + ii1] *= np.arange(fn)[::-1]/fn
938 prev_idx = idx
939 freq, power = psd(data - np.mean(data), rate, freq_resolution)
940 #power *= n/rate/ipi/len(eod_times)
941 return freq, power
944def analyze_pulse_spectrum(freqs, energy):
945 """Analyze the spectrum of a pulse-type EOD.
947 Parameters
948 ----------
949 freqs: 1-D array of float
950 The frequency components of the energy spectrum.
951 energy: 1-D array of float
952 The energy spectrum of the single pulse-type EOD.
954 Returns
955 -------
956 peak_freq: float
957 Frequency at peak energy of the spectrum in Hertz.
958 peak_energy: float
959 Peak energy of the pulse spectrum in x^2 s/Hz.
960 trough_freq: float
961 Frequency at trough before peak in Hertz.
962 trough_energy: float
963 Energy of trough before peak in x^2 s/Hz.
964 att5: float
965 Attenuation of average energy below 5 Hz relative to
966 peak energy in decibel.
967 att50: float
968 Attenuation of average energy below 50 Hz relative to
969 peak energy in decibel.
970 low_cutoff: float
971 Frequency at which the energy reached half of the
972 peak energy relative to the DC energy in Hertz.
973 high_cutoff: float
974 3dB roll-off frequency in Hertz.
975 """
976 ip = np.argmax(energy)
977 peak_freq = freqs[ip]
978 peak_energy = energy[ip]
979 it = np.argmin(energy[:ip]) if ip > 0 else 0
980 trough_freq = freqs[it]
981 trough_energy = energy[it]
982 att5 = decibel(np.mean(energy[freqs<5.0]), peak_energy)
983 att50 = decibel(np.mean(energy[freqs<50.0]), peak_energy)
984 low_cutoff = freqs[decibel(energy, peak_energy) > 0.5*att5][0]
985 high_cutoff = freqs[decibel(energy, peak_energy) > -3.0][-1]
986 return peak_freq, peak_energy, trough_freq, trough_energy, \
987 att5, att50, low_cutoff, high_cutoff
990def analyze_pulse_intervals(eod_times, ipi_cv_thresh=0.5,
991 ipi_percentile=30.0):
992 """ Basic statistics of interpulse intervals.
994 Parameters
995 ----------
996 eod_times: 1-D array or None
997 List of times of detected EODs.
998 ipi_cv_thresh: float
999 If the coefficient of variation of the interpulse intervals
1000 (IPIs) is smaller than this threshold, then the statistics of
1001 IPIs is estimated from all IPIs. Otherwise only intervals
1002 smaller than a certain percentile are used.
1003 ipi_percentile: float
1004 When computing the statistics of IPIs from a subset of the
1005 IPIs, only intervals smaller than this percentile (between 0
1006 and 100) are used.
1008 Returns
1009 -------
1010 ipi_median: float
1011 Median inter-pulse interval.
1012 ipi_mean: float
1013 Mean inter-pulse interval.
1014 ipi_std: float
1015 Standard deviation of inter-pulse intervals.
1017 """
1018 if eod_times is None:
1019 return None, None, None
1020 inter_pulse_intervals = np.diff(eod_times)
1021 ipi_cv = np.std(inter_pulse_intervals)/np.mean(inter_pulse_intervals)
1022 if ipi_cv < ipi_cv_thresh:
1023 ipi_median = np.median(inter_pulse_intervals)
1024 ipi_mean = np.mean(inter_pulse_intervals)
1025 ipi_std = np.std(inter_pulse_intervals)
1026 else:
1027 intervals = inter_pulse_intervals[inter_pulse_intervals <
1028 np.percentile(inter_pulse_intervals, ipi_percentile)]
1029 ipi_median = np.median(intervals)
1030 ipi_mean = np.mean(intervals)
1031 ipi_std = np.std(intervals)
1032 return ipi_median, ipi_mean, ipi_std
1035def analyze_pulse(eod, ratetime=None, eod_times=None,
1036 min_pulse_win=0.001,
1037 start_end_thresh_fac=0.01, peak_thresh_fac=0.002,
1038 min_dist=50.0e-6, width_frac=0.5, fit_frac=0.5,
1039 freq_resolution=1.0, fade_frac=0.0,
1040 flip_pulse='none', fit_gaussians=True, ipi_cv_thresh=0.5,
1041 ipi_percentile=30.0, verbose=0):
1042 """Analyze the EOD waveform of a pulse fish.
1044 Parameters
1045 ----------
1046 eod: 1-D or 2-D array
1047 The eod waveform to be analyzed.
1048 If an 1-D array, then this is the waveform and you
1049 need to also pass a time array or sampling rate in `ratetime`.
1050 If a 2-D array, then first column is time in seconds, second
1051 column the EOD waveform, third column, if present, is the
1052 standard error of the EOD waveform. Further columns are
1053 optional and are not used.
1054 ratetime: None or float or array of float
1055 If a 1-D array is passed on to `eod` then either the sampling
1056 rate in Hertz or the time array corresponding to `eod`.
1057 eod_times: 1-D array or None
1058 List of times of detected EOD peaks.
1059 min_pulse_win: float
1060 The minimum size of cut-out EOD waveform.
1061 start_end_thresh_fac: float
1062 Set the threshold for the start and end time to the p-p amplitude
1063 times this factor.
1064 peak_thresh_fac: float
1065 Set the threshold for peak and trough detection to the p-p amplitude
1066 times this factor.
1067 min_dist: float
1068 Minimum distance between peak and troughs of the pulse.
1069 width_frac: float
1070 The width of a peak is measured at this fraction of a peak's
1071 height (0-1).
1072 fit_frac: float or None
1073 An exponential is fitted to the tail of the last peak/trough
1074 starting where the waveform falls below this fraction of the
1075 peak's height (0-1).
1076 freq_resolution: float
1077 The frequency resolution of the spectrum of the single pulse.
1078 fade_frac: float
1079 Fraction of time of the EOD waveform that is fade in and out
1080 to zero baseline.
1081 flip_pulse: 'auto', 'none', 'flip'
1082 - 'auto' flip waveform such that the first large extremum is positive.
1083 - 'flip' flip waveform.
1084 - 'none' do not flip waveform.
1085 fit_gaussians: bool
1086 If True, fit sum of Gaussians to pulse waveform.
1087 Otherwise, do not fit, and return empty dictionary for `pulses`.
1088 ipi_cv_thresh: float
1089 If the coefficient of variation of the interpulse intervals
1090 (IPIs) is smaller than this threshold, then the statistics of
1091 IPIs is estimated from all IPIs. Otherwise only intervals
1092 smaller than a certain percentile are used.
1093 ipi_percentile: float
1094 When computing the statistics of IPIs from a subset of the
1095 IPIs, only intervals smaller than this percentile (between 0
1096 and 100) are used.
1097 verbose: int
1098 Verbosity level passed for error and info messages.
1100 Returns
1101 -------
1102 meod: 2-D array of floats
1103 The eod waveform. First column is time in seconds,
1104 second column the eod waveform.
1105 Further columns are kept from the input `eod`.
1106 As the two last columns the waveform resulting from the
1107 decomposition into Gaussians and the fit to the tail of the
1108 last peak are appended.
1109 props: dict
1110 A dictionary with properties of the analyzed EOD waveform.
1112 - type: set to 'pulse'.
1113 - flipped: True if the waveform was flipped.
1114 - n: number of pulses analyzed (i.e. `len(eod_times)`), if provided.
1115 - times: the times of the detected EOD pulses (i.e. `eod_times`),
1116 if provided.
1117 - EODf: the inverse of the median interval between `eod_times`,
1118 if provided.
1119 - period: median interval between `eod_times`, if provided.
1120 - IPI-mean: mean interval between `eod_times`, if provided.
1121 - IPI-std: standard deviation of the intervals between
1122 `eod_times`, if provided.
1123 - IPI-CV: coefficient of variation of the intervals between
1124 `eod_times`, if provided.
1125 - aoffs: Offset that was subtracted from the average EOD waveform.
1126 - pos-ampl: amplitude of the largest positive peak.
1127 - neg-ampl: amplitude of the largest negative trough.
1128 - max-ampl: maximum of largest peak or trough amplitude in the units of the input data.
1129 - p-p-amplitude: peak-to-peak amplitude of the EOD waveform.
1130 - p-p-dist: distance between minimum and maximum phase in seconds.
1131 - noise: average standard error mean of the averaged
1132 EOD waveform relative to the p-p amplitude.
1133 - noise: average standard error of the averaged EOD waveform relative to the peak-to_peak amplitude in percent.
1134 - rmserror: root-mean-square error between fit with sum of Gaussians and
1135 EOD waveform relative to the p-p amplitude. Infinity if fit failed.
1136 - peakthresh: Threshold for detecting peaks and troughs is at this factor times p-p-ampl.
1137 - startendthresh: Threshold for start and end time is at this factor times p-p-ampl.
1138 - tstart: time in seconds where the pulse starts,
1139 i.e. crosses the threshold for the first time.
1140 - tend: time in seconds where the pulse ends,
1141 i.e. crosses the threshold for the last time.
1142 - width: total width of the pulse in seconds (tend-tstart).
1143 - totalarea: sum of areas under positive and negative peaks.
1144 - pos-area: area under positive phases relative to total area.
1145 - neg-area: area under negative phases relative to total area.
1146 - polaritybalance: contrast between areas under positive and
1147 negative phases.
1148 - median: median of the distribution of the absolute EOD waveform.
1149 - quartile1: first quartile of the distribution of the absolute EOD waveform.
1150 - quartile3: third quartile of the distribution of the absolute EOD waveform.
1151 - iq-range: inter-quartile range of the distribution of the absolute EOD waveform.
1152 - tau: time constant of exponential decay of pulse tail in seconds.
1153 - firstphase: index of the first phase in the pulse (i.e. -1 for P-1)
1154 - lastphase: index of the last phase in the pulse (i.e. 3 for P3)
1155 - peakfreq: frequency at peak energy of the single pulse spectrum
1156 in Hertz.
1157 - peakenergy: peak energy of the single pulse spectrum.
1158 - troughfreq: frequency at trough before peak in Hertz.
1159 - troughenergy: energy of trough before peak in x^2 s/Hz.
1160 - energyatt5: attenuation of average energy below 5 Hz relative to
1161 peak energy in decibel.
1162 - energyatt50: attenuation of average energy below 50 Hz relative to
1163 peak energy in decibel.
1164 - lowcutoff: frequency at which the energy reached half of the
1165 peak energy relative to the DC energy in Hertz.
1166 - highcutoff: 3dB roll-off frequency of spectrum in Hertz.
1168 Empty if waveform is not a pulse EOD.
1169 phases: dict
1170 Dictionary with
1172 - "indices": indices of each phase
1173 (1 is P1, i.e. the largest positive peak)
1174 - "times": times of each phase relative to P1 in seconds
1175 - "amplitudes": amplitudes of each phase
1176 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
1177 - "widths": widths of each phase at `width_frac` height
1178 - "areas": areas of each phase
1179 - "relareas": areas of the phases relative to total area
1180 - "zeros": time of zero crossing towards next phase in seconds
1182 Empty dictionary if waveform is not a pulse EOD.
1183 pulse: dict
1184 Dictionary with
1186 - "times": phase times in seconds,
1187 - "amplitudes": amplitudes, and
1188 - "stdevs": standard deviations in seconds
1190 of Gaussians fitted to the pulse waveform. Use the functions
1191 provided in thunderfish.fakefish to simulate pulse fish EODs
1192 from this data. Empty dictionary if `fit_gaussians` is False.
1193 energy: 2-D array
1194 The energy spectrum of a single pulse. First column are the
1195 frequencies, second column the energy in x^2 s/Hz.
1196 Empty if waveform is not a pulse EOD.
1198 """
1199 if eod.ndim == 2:
1200 if eod.shape[1] > 2:
1201 eeod = eod
1202 else:
1203 eeod = np.column_stack((eod, np.zeros(len(eod))))
1204 else:
1205 if isinstance(ratetime, (list, tuple, np.ndarray)):
1206 time = ratetime
1207 else:
1208 time = np.arange(len(eod))/ratetime
1209 eeod = np.zeros((len(eod), 3))
1210 eeod[:, 0] = time
1211 eeod[:, 1] = eod
1212 # storage:
1213 meod = np.zeros((eeod.shape[0], eeod.shape[1] + 2))
1214 meod[:, :eeod.shape[1]] = eeod
1215 meod[:, -2] = np.nan
1216 meod[:, -1] = np.nan
1218 # conditioning of the waveform:
1219 meod, toffs, aoffs, flipped, noise_thresh = \
1220 condition_pulse(meod, flip=flip_pulse,
1221 baseline_frac=0.05, large_phase_frac=0.2,
1222 min_pulse_win=min_pulse_win)
1224 # analysis of pulse waveform:
1225 pos_ampl, neg_ampl, dist, pos_area, neg_area, \
1226 polarity_balance, median, quartile1, quartile3 = \
1227 analyze_pulse_properties(noise_thresh, meod)
1228 pp_ampl = pos_ampl + neg_ampl
1229 max_ampl = max(pos_ampl, neg_ampl)
1230 total_area = pos_area + neg_area
1232 # threshold for start and end time:
1233 start_end_thresh = pp_ampl*start_end_thresh_fac
1234 if start_end_thresh < 2*noise_thresh:
1235 start_end_thresh = 2*noise_thresh
1236 start_end_thresh_fac = start_end_thresh/pp_ampl if pp_ampl > 0 else 1
1238 # threshold for peak detection:
1239 peak_thresh = pp_ampl*peak_thresh_fac
1240 if peak_thresh < noise_thresh:
1241 peak_thresh = noise_thresh
1242 peak_thresh_fac = peak_thresh/pp_ampl if pp_ampl > 0 else 1
1244 # characterize EOD phases:
1245 tstart, tend, phases = analyze_pulse_phases(peak_thresh,
1246 start_end_thresh, meod,
1247 min_dist=min_dist,
1248 width_frac=width_frac)
1250 # fit exponential to last phase:
1251 tau = np.nan
1252 taustart = np.nan
1253 if len(phases) > 0 and len(phases['times']) > 1:
1254 pi = np.argmin(np.abs(meod[:, 0] - phases['times'][-1]))
1255 tau, taustart, fit = analyze_pulse_tail(pi, meod, None,
1256 threshold=noise_thresh,
1257 fit_frac=fit_frac,
1258 verbose=verbose)
1259 if fit is not None:
1260 meod[:, -1] = fit
1262 # energy spectrum of single EOD pulse:
1263 freqs, energy = pulse_spectrum(meod, None, freq_resolution, fade_frac)
1264 # store spectrum:
1265 eenergy = np.zeros((len(energy), 2))
1266 eenergy[:, 0] = freqs
1267 eenergy[:, 1] = energy
1268 # analyse spectrum:
1269 peakfreq, peakenergy, troughfreq, troughenergy, \
1270 att5, att50, lowcutoff, highcutoff = \
1271 analyze_pulse_spectrum(freqs, energy)
1273 # decompose EOD waveform:
1274 rmserror = np.inf
1275 pulse = {}
1276 if fit_gaussians:
1277 pulse = decompose_pulse(meod, freqs, energy, phases,
1278 width_frac, verbose=verbose)
1279 if len(pulse) > 0:
1280 eod_fit = pulsefish_waveform(pulse, meod[:, 0])
1281 meod[:, -2] = eod_fit
1282 rmserror = np.sqrt(np.mean((meod[:, 1] - meod[:, -2])**2))/pp_ampl
1283 spec = pulsefish_spectrum(pulse, freqs)
1284 spec = np.abs(spec)**2
1285 eenergy = np.hstack((eenergy, spec.reshape((-1, 1))))
1287 # analyze pulse intervals:
1288 ipi_median, ipi_mean, ipi_std = \
1289 analyze_pulse_intervals(eod_times, ipi_cv_thresh, ipi_percentile)
1291 # store properties:
1292 props = {}
1293 props['type'] = 'pulse'
1294 props['flipped'] = flipped
1295 if eod_times is not None:
1296 props['n'] = len(eod_times)
1297 props['times'] = eod_times + toffs
1298 props['EODf'] = 1.0/ipi_median
1299 props['period'] = ipi_median
1300 props['IPI-mean'] = ipi_mean
1301 props['IPI-std'] = ipi_std
1302 props['IPI-CV'] = ipi_std/ipi_mean
1303 props['aoffs'] = aoffs
1304 props['pos-ampl'] = pos_ampl
1305 props['neg-ampl'] = neg_ampl
1306 props['max-ampl'] = max_ampl
1307 props['p-p-amplitude'] = pp_ampl
1308 props['p-p-dist'] = dist
1309 if meod.shape[1] > 2:
1310 props['noise'] = np.mean(meod[:, 2])/pp_ampl if pp_ampl > 0 else 1
1311 props['rmserror'] = rmserror
1312 props['peakthresh'] = peak_thresh_fac
1313 props['startendthresh'] = start_end_thresh_fac
1314 props['tstart'] = tstart
1315 props['tend'] = tend
1316 props['width'] = tstart - tend
1317 props['totalarea'] = total_area
1318 props['pos-area'] = pos_area/total_area
1319 props['neg-area'] = neg_area/total_area
1320 props['polaritybalance'] = polarity_balance
1321 props['median'] = median
1322 props['quartile1'] = quartile1
1323 props['quartile3'] = quartile3
1324 props['iq-range'] = quartile3 - quartile1
1325 props['tau'] = tau
1326 props['taustart'] = taustart
1327 props['firstphase'] = phases['indices'][0] if len(phases) > 0 else 1
1328 props['lastphase'] = phases['indices'][-1] if len(phases) > 0 else 1
1329 props['peakfreq'] = peakfreq
1330 props['peakenergy'] = peakenergy
1331 props['troughfreq'] = troughfreq
1332 props['troughenergy'] = troughenergy
1333 props['energyatt5'] = att5
1334 props['energyatt50'] = att50
1335 props['lowcutoff'] = lowcutoff
1336 props['highcutoff'] = highcutoff
1338 return meod, props, phases, pulse, eenergy
1341def plot_pulse_eods(ax, data, rate, zoom_window, width, eod_props,
1342 toffs=0.0, colors=None, markers=None, marker_size=10,
1343 legend_rows=8, **kwargs):
1344 """Mark pulse EODs in a plot of an EOD recording.
1346 Parameters
1347 ----------
1348 ax: matplotlib axes
1349 Axes used for plotting.
1350 data: 1D ndarray
1351 Recorded data (these are not plotted!).
1352 rate: float
1353 Sampling rate of the data in Hertz.
1354 zoom_window: tuple of floats
1355 Start and end time of the data to be plotted in seconds.
1356 width: float
1357 Minimum width of the data to be plotted in seconds.
1358 eod_props: list of dictionaries
1359 Lists of EOD properties as returned by `analyze_pulse()`
1360 and `analyze_wave()`. From the entries with 'type' ==
1361 'pulse' the properties 'EODf' and 'times' are used. 'EODf'
1362 is the averaged EOD frequency, and 'times' is a list of
1363 detected EOD pulse times.
1364 toffs: float
1365 Time of first data value in seconds that will be added
1366 to the pulse times in `eod_props`.
1367 colors: list of colors or None
1368 If not None list of colors for plotting each cluster
1369 markers: list of markers or None
1370 If not None list of markers for plotting each cluster
1371 marker_size: float
1372 Size of markers used to mark the pulses.
1373 legend_rows: int
1374 Maximum number of rows to be used for the legend.
1375 kwargs:
1376 Key word arguments for the legend of the plot.
1377 """
1378 k = 0
1379 for eod in eod_props:
1380 if eod['type'] != 'pulse':
1381 continue
1382 if 'times' not in eod:
1383 continue
1385 width = np.min([width, np.diff(zoom_window)[0]])
1386 while len(eod['peaktimes'][(eod['peaktimes']>(zoom_window[1]-width)) & (eod['peaktimes']<(zoom_window[1]))]) == 0:
1387 width *= 2
1388 if zoom_window[1] - width < 0:
1389 width = width/2
1390 break
1392 x = eod['peaktimes'] + toffs
1393 pidx = np.round(eod['peaktimes']*rate).astype(int)
1394 y = data[pidx[(pidx>0)&(pidx<len(data))]]
1395 x = x[(pidx>0)&(pidx<len(data))]
1396 color_kwargs = {}
1397 #if colors is not None:
1398 # color_kwargs['color'] = colors[k%len(colors)]
1399 if marker_size is not None:
1400 color_kwargs['ms'] = marker_size
1401 label = '%6.1f Hz' % eod['EODf']
1402 if legend_rows > 5 and k >= legend_rows:
1403 label = None
1404 if markers is None:
1405 ax.plot(x, y, 'o', label=label, zorder=-1, **color_kwargs)
1406 else:
1407 ax.plot(x, y, linestyle='none', label=label,
1408 marker=markers[k%len(markers)], mec=None, mew=0.0,
1409 zorder=-1, **color_kwargs)
1410 k += 1
1412 # legend:
1413 if k > 1:
1414 if legend_rows > 0:
1415 if legend_rows > 5:
1416 ncol = 1
1417 else:
1418 ncol = (len(idx)-1) // legend_rows + 1
1419 ax.legend(numpoints=1, ncol=ncol, **kwargs)
1420 else:
1421 ax.legend(numpoints=1, **kwargs)
1423 # reset window so at least one EOD of each cluster is visible
1424 if len(zoom_window)>0:
1425 ax.set_xlim(max(toffs,toffs+zoom_window[1]-width), toffs+zoom_window[1])
1427 i0 = max(0,int((zoom_window[1]-width)*rate))
1428 i1 = int(zoom_window[1]*rate)
1430 ymin = np.min(data[i0:i1])
1431 ymax = np.max(data[i0:i1])
1432 dy = ymax - ymin
1433 ax.set_ylim(ymin-0.05*dy, ymax+0.05*dy)
1436def plot_pulse_spectrum(ax, energy, props, min_freq=1.0, max_freq=10000.0,
1437 min_db = -60, ref_energy=None,
1438 spec_style=dict(lw=3, color='tab:blue'),
1439 analytic_style=dict(lw=4, color='tab:cyan'),
1440 peak_style=dict(ls='', marker='o', markersize=6,
1441 color='tab:blue', mec='none', mew=0,
1442 alpha=0.4),
1443 cutoff_style=dict(lw=1, ls='-', color='0.5'),
1444 att5_color='0.8', att50_color='0.9',
1445 fontsize='medium'):
1446 """Plot and annotate spectrum of single pulse EOD.
1448 Parameters
1449 ----------
1450 ax: matplotlib axes
1451 Axes used for plotting.
1452 energy: 2-D array
1453 The energy spectrum of a single pulse as returned by `analyze_pulse()`.
1454 First column are the frequencies, second column the energy.
1455 An optional third column is an analytically computed spectrum.
1456 props: dict
1457 A dictionary with properties of the analyzed EOD waveform as
1458 returned by `analyze_pulse()`.
1459 min_freq: float
1460 Minimun frequency of the spectrum to be plotted (logscale!).
1461 max_freq: float
1462 Maximun frequency of the spectrum to be plotted (logscale!).
1463 min_db: float
1464 Minimum decibel level shown.
1465 ref_energy: float or None
1466 Reference energy for computing decibel.
1467 If None take the maximum energy.
1468 spec_style: dict
1469 Arguments passed on to the plot command for the energy spectrum
1470 computed from the data.
1471 analytic_style: dict
1472 Arguments passed on to the plot command for the energy spectrum
1473 that was analytically computed from the Gaussian fits
1474 (optional third column in `energy`).
1475 peak_style: dict
1476 Arguments passed on to the plot commands for marking the peak
1477 and trough frequency.
1478 cutoff_style: dict
1479 Arguments passed on to the plot command for the lines marking
1480 cutoff frequencies.
1481 att5_color: matplotlib color specification
1482 Color for the rectangular patch marking the first 5 Hz.
1483 att50_color: matplotlib color specification
1484 Color for the rectangular patch marking the first 50 Hz.
1485 fontsize: str or float or int
1486 Fontsize for annotation text.
1488 Returns
1489 -------
1490 ref_energy: float
1491 The actually used reference energy for computing decibels.
1492 """
1493 ax.axvspan(1, 50, color=att50_color, zorder=-20)
1494 att = props['energyatt50']
1495 if att < -10:
1496 y = att + 1
1497 if y < min_db:
1498 y = min_db + 2
1499 ax.text(10, y, f'{att:.0f}dB',
1500 ha='left', va='bottom', zorder=100, fontsize=fontsize)
1501 else:
1502 ax.text(10, att - 1, f'{att:.0f}dB',
1503 ha='left', va='top', zorder=100, fontsize=fontsize)
1504 ax.axvspan(1, 5, color=att5_color, zorder=-10)
1505 att = props['energyatt5']
1506 if att < -10:
1507 y = att + 1
1508 if y < min_db:
1509 y = min_db + 2
1510 ax.text(4, y, f'{att:.0f}dB',
1511 ha='right', va='bottom', zorder=100, fontsize=fontsize)
1512 else:
1513 ax.text(4, att - 1, f'{att:.0f}dB',
1514 ha='right', va='top', zorder=100, fontsize=fontsize)
1515 lowcutoff = props['lowcutoff']
1516 if lowcutoff >= min_freq:
1517 ax.plot([lowcutoff, lowcutoff, 1], [min_db, 0.5*att, 0.5*att],
1518 zorder=30, **cutoff_style)
1519 ax.text(1.2*lowcutoff, 0.5*att - 1, f'{lowcutoff:.0f}Hz',
1520 ha='left', va='top', zorder=100, fontsize=fontsize)
1521 highcutoff = props['highcutoff']
1522 ax.plot([highcutoff, highcutoff], [min_db, -3], zorder=30, **cutoff_style)
1523 ax.text(1.2*highcutoff, -3, f'{highcutoff:.0f}Hz',
1524 ha='left', va='center', zorder=100, fontsize=fontsize)
1525 if ref_energy is None:
1526 ref_energy = np.max(energy[:, 1])
1527 if energy.shape[1] > 2 and np.all(np.isfinite(energy[:, 2])) and len(analytic_style) > 0:
1528 db = decibel(energy[:, 2], ref_energy)
1529 ax.plot(energy[:, 0], db, zorder=45, **analytic_style)
1530 db = decibel(energy[:, 1], ref_energy)
1531 ax.plot(energy[:, 0], db, zorder=50, **spec_style)
1532 peakfreq = props['peakfreq']
1533 if peakfreq >= min_freq:
1534 ax.plot(peakfreq, 0, zorder=60, clip_on=False, **peak_style)
1535 ax.text(peakfreq*1.2, 1, f'{peakfreq:.0f}Hz',
1536 va='bottom', zorder=100, fontsize=fontsize)
1537 troughfreq = props['troughfreq']
1538 if troughfreq >= min_freq:
1539 troughenergy = decibel(props['troughenergy'], ref_energy)
1540 ax.plot(troughfreq, troughenergy, zorder=60, **peak_style)
1541 ax.text(troughfreq, troughenergy - 3,
1542 f'{troughenergy:.0f}dB @ {troughfreq:.0f}Hz',
1543 ha='center', va='top', zorder=100, fontsize=fontsize)
1544 ax.set_xlim(min_freq, max_freq)
1545 ax.set_xscale('log')
1546 ax.set_ylim(min_db, 2)
1547 ax.set_xlabel('Frequency [Hz]')
1548 ax.set_ylabel('Energy [dB]')
1549 return ref_energy
1552def save_pulse_fish(eod_props, unit, basename, **kwargs):
1553 """Save properties of pulse EODs to file.
1555 Parameters
1556 ----------
1557 eod_props: list of dict
1558 Properties of EODs as returned by `analyze_wave()` and
1559 `analyze_pulse()`. Only properties of pulse fish are saved.
1560 unit: string
1561 Unit of the waveform data.
1562 basename: string or stream
1563 If string, path and basename of file.
1564 If `basename` does not have an extension,
1565 '-pulsefish' and a file extension are appended.
1566 If stream, write pulse fish properties into this stream.
1567 kwargs:
1568 Arguments passed on to `TableData.write()`.
1570 Returns
1571 -------
1572 filename: Path or None
1573 Path and full name of the written file in case of `basename`
1574 being a string. Otherwise, the file name and extension that
1575 would have been appended to a basename.
1576 None if no pulse fish are contained in eod_props and
1577 consequently no file was written.
1579 See Also
1580 --------
1581 load_pulse_fish()
1582 """
1583 pulse_props = [p for p in eod_props if p['type'] == 'pulse']
1584 if len(pulse_props) == 0:
1585 return None
1586 td = TableData()
1587 if 'twin' in pulse_props[0] or 'samplerate' in pulse_props[0] or \
1588 'nfft' in pulse_props[0]:
1589 td.append_section('recording')
1590 if 'twin' in pulse_props[0]:
1591 td.append('twin', 's', '%7.2f', value=pulse_props)
1592 td.append('window', 's', '%7.2f', value=pulse_props)
1593 td.append('winclipped', '%', '%.2f',
1594 value=pulse_props, fac=100)
1595 if 'samplerate' in pulse_props[0]:
1596 td.append('samplerate', 'kHz', '%.3f', value=pulse_props,
1597 fac=0.001)
1598 if 'nfft' in pulse_props[0]:
1599 td.append('nfft', '', '%d', value=pulse_props)
1600 td.append('dfreq', 'Hz', '%.2f', value=pulse_props)
1601 td.append_section('waveform')
1602 td.append('index', '', '%d', value=pulse_props)
1603 td.append('n', '', '%d', value=pulse_props)
1604 td.append('EODf', 'Hz', '%7.2f', value=pulse_props)
1605 td.append('period', 'ms', '%7.2f', value=pulse_props, fac=1000)
1606 td.append('aoffs', unit, '%.5f', value=pulse_props)
1607 td.append('pos-ampl', unit, '%.5f', value=pulse_props)
1608 td.append('neg-ampl', unit, '%.5f', value=pulse_props)
1609 td.append('max-ampl', unit, '%.5f', value=pulse_props)
1610 td.append('p-p-amplitude', unit, '%.5f', value=pulse_props)
1611 td.append('p-p-dist', 'ms', '%.3f', value=pulse_props, fac=1000)
1612 if 'noise' in pulse_props[0]:
1613 td.append('noise', '%', '%.2f', value=pulse_props, fac=100)
1614 td.append('rmserror', '%', '%.2f', value=pulse_props, fac=100)
1615 if 'clipped' in pulse_props[0]:
1616 td.append('clipped', '%', '%.2f', value=pulse_props, fac=100)
1617 td.append('flipped', '', '%d', value=pulse_props)
1618 td.append('startendthresh', '%', '%.2f', value=pulse_props, fac=100)
1619 td.append('peakthresh', '%', '%.2f', value=pulse_props, fac=100)
1620 td.append('tstart', 'ms', '%.3f', value=pulse_props, fac=1000)
1621 td.append('tend', 'ms', '%.3f', value=pulse_props, fac=1000)
1622 td.append('width', 'ms', '%.3f', value=pulse_props, fac=1000)
1623 td.append('totalarea', f'{unit}*ms', '%.4f', value=pulse_props, fac=1000)
1624 td.append('pos-area', '%', '%.2f', value=pulse_props, fac=100)
1625 td.append('neg-area', '%', '%.2f', value=pulse_props, fac=100)
1626 td.append('polaritybalance', '%', '%.2f', value=pulse_props, fac=100)
1627 td.append('median', 'ms', '%.3f', value=pulse_props, fac=1000)
1628 td.append('quartile1', 'ms', '%.3f', value=pulse_props, fac=1000)
1629 td.append('quartile3', 'ms', '%.3f', value=pulse_props, fac=1000)
1630 td.append('iq-range', 'ms', '%.3f', value=pulse_props, fac=1000)
1631 td.append('tau', 'ms', '%.3f', value=pulse_props, fac=1000)
1632 td.append('taustart', 'ms', '%.3f', value=pulse_props, fac=1000)
1633 td.append('firstphase', '', '%d', value=pulse_props)
1634 td.append('lastphase', '', '%d', value=pulse_props)
1635 td.append_section('spectrum')
1636 td.append('peakfreq', 'Hz', '%.2f', value=pulse_props)
1637 td.append('peakenergy', f'{unit}^2s/Hz', '%.3g', value=pulse_props)
1638 td.append('troughfreq', 'Hz', '%.2f', value=pulse_props)
1639 td.append('troughenergy', f'{unit}^2s/Hz', '%.3g', value=pulse_props)
1640 td.append('energyatt5', 'dB', '%.2f', value=pulse_props)
1641 td.append('energyatt50', 'dB', '%.2f', value=pulse_props)
1642 td.append('lowcutoff', 'Hz', '%.2f', value=pulse_props)
1643 td.append('highcutoff', 'Hz', '%.2f', value=pulse_props)
1644 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
1645 fp = '-pulsefish' if not ext else ''
1646 return td.write_file_stream(basename, fp, **kwargs)
1649def load_pulse_fish(file_path):
1650 """Load properties of pulse EODs from file.
1652 All times are scaled to seconds, all frequencies to Hertz, and all
1653 percentages to fractions.
1655 Parameters
1656 ----------
1657 file_path: string
1658 Path of the file to be loaded.
1660 Returns
1661 -------
1662 eod_props: list of dict
1663 Properties of EODs.
1665 Raises
1666 ------
1667 FileNotFoundError:
1668 If `file_path` does not exist.
1670 See Also
1671 --------
1672 save_pulse_fish()
1674 """
1675 data = TableData(file_path)
1676 eod_props = data.dicts()
1677 for props in eod_props:
1678 if 'winclipped' in props:
1679 props['winclipped'] /= 100
1680 if 'samplerate' in props:
1681 props['samplerate'] *= 1000
1682 if 'nfft' in props:
1683 props['nfft'] = int(props['nfft'])
1684 if 'clipped' in props:
1685 props['clipped'] /= 100
1686 props['type'] = 'pulse'
1687 props['index'] = int(props['index'])
1688 props['n'] = int(props['n'])
1689 props['totalarea'] /= 1000
1690 props['pos-area'] /= 100
1691 props['neg-area'] /= 100
1692 props['polaritybalance'] /= 100
1693 props['median'] /= 1000
1694 props['quartile1'] /= 1000
1695 props['quartile3'] /= 1000
1696 props['iq-range'] /= 1000
1697 props['firstphase'] = int(props['firstphase'])
1698 props['lastphase'] = int(props['lastphase'])
1699 props['period'] /= 1000
1700 props['noise'] /= 100
1701 props['startendthresh'] /= 100
1702 props['peakthresh'] /= 100
1703 props['tstart'] /= 1000
1704 props['tend'] /= 1000
1705 props['p-p-dist'] /= 1000
1706 props['width'] /= 1000
1707 props['tau'] /= 1000
1708 props['taustart'] /= 1000
1709 props['rmserror'] /= 100
1710 return eod_props
1713def save_pulse_spectrum(spec_data, unit, idx, basename, **kwargs):
1714 """Save energy spectrum of pulse EOD to file.
1716 Parameters
1717 ----------
1718 spec_data: 2D array of floats
1719 Energy spectrum of single pulse as returned by `analyze_pulse()`.
1720 unit: string
1721 Unit of the waveform data.
1722 idx: int or None
1723 Index of fish.
1724 basename: string or stream
1725 If string, path and basename of file.
1726 If `basename` does not have an extension,
1727 '-pulsespectrum', the fish index, and a file extension are appended.
1728 If stream, write pulse spectrum into this stream.
1729 kwargs:
1730 Arguments passed on to `TableData.write()`.
1732 Returns
1733 -------
1734 filename: Path
1735 Path and full name of the written file in case of `basename`
1736 being a string. Otherwise, the file name and extension that
1737 would have been appended to a basename.
1739 See Also
1740 --------
1741 load_pulse_spectrum()
1742 """
1743 td = TableData(spec_data[:, :2], ['frequency', 'energy'],
1744 ['Hz', f'{unit}^2s/Hz'], ['%.2f', '%.4e'])
1745 if spec_data.shape[1] > 2:
1746 td.append('gaussian', f'{unit}^2s/Hz', '%.4e', value=spec_data[:, 2])
1747 fp = ''
1748 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
1749 if not ext:
1750 fp = '-pulsespectrum'
1751 if idx is not None:
1752 fp += f'-{idx}'
1753 return td.write_file_stream(basename, fp, **kwargs)
1756def load_pulse_spectrum(file_path):
1757 """Load energy spectrum of pulse EOD from file.
1759 Parameters
1760 ----------
1761 file_path: string
1762 Path of the file to be loaded.
1764 Returns
1765 -------
1766 spec: 2D array of floats
1767 Energy spectrum of single pulse: frequency, energy
1769 Raises
1770 ------
1771 FileNotFoundError:
1772 If `file_path` does not exist.
1774 See Also
1775 --------
1776 save_pulse_spectrum()
1777 """
1778 data = TableData(file_path)
1779 spec = data.array()
1780 return spec
1783def save_pulse_phases(phases, unit, idx, basename, **kwargs):
1784 """Save phase properties of pulse EOD to file.
1786 Parameters
1787 ----------
1788 phases: dict
1789 Dictionary with
1791 - "indices": indices of each phase
1792 (1 is P1, i.e. the largest positive peak)
1793 - "times": times of each phase relative to P1 in seconds
1794 - "amplitudes": amplitudes of each phase
1795 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
1796 - "widths": widths of each phase at `width_frac` height
1797 - "areas": areas of each phase
1798 - "relareas": areas of the phases relative to total area
1799 - "zeros": time of zero crossing towards next phase in seconds
1801 as returned by `analyze_pulse_phases()` and `analyze_pulse()`.
1802 unit: string
1803 Unit of the waveform data.
1804 idx: int or None
1805 Index of fish.
1806 basename: string or stream
1807 If string, path and basename of file.
1808 If `basename` does not have an extension,
1809 '-pulsephases', the fish index, and a file extension are appended.
1810 If stream, write pulse phases into this stream.
1811 kwargs:
1812 Arguments passed on to `TableData.write()`.
1814 Returns
1815 -------
1816 filename: Path
1817 Path and full name of the written file in case of `basename`
1818 being a string. Otherwise, the file name and extension that
1819 would have been appended to a basename.
1821 See Also
1822 --------
1823 load_pulse_phases()
1824 """
1825 if len(phases) == 0:
1826 return None
1827 td = TableData()
1828 td.append('type', '', '%s', value=['P']*len(phases['indices']))
1829 td.append('index', '', '%.0f', value=phases['indices'])
1830 td.append('time', 'ms', '%.4f', value=phases['times'], fac=1000)
1831 td.append('amplitude', unit, '%.5f', value=phases['amplitudes'])
1832 td.append('relampl', '%', '%.2f', value=phases['relamplitudes'], fac=100)
1833 td.append('width', 'ms', '%.4f', value=phases['widths'], fac=1000)
1834 td.append('area', f'{unit}*ms', '%.4f', value=phases['areas'], fac=1000)
1835 td.append('relarea', '%', '%.2f', value=phases['relareas'], fac=100)
1836 td.append('zeros', 'ms', '%.4f', value=phases['zeros'], fac=1000)
1837 fp = ''
1838 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
1839 if not ext:
1840 fp = '-pulsephases'
1841 if idx is not None:
1842 fp += f'-{idx}'
1843 return td.write_file_stream(basename, fp, **kwargs)
1846def load_pulse_phases(file_path):
1847 """Load phase properties of pulse EOD from file.
1849 Parameters
1850 ----------
1851 file_path: string
1852 Path of the file to be loaded.
1854 Returns
1855 -------
1856 phases: dict
1857 Dictionary with
1859 - "indices": indices of each phase
1860 (1 is P1, i.e. the largest positive peak)
1861 - "times": times of each phase relative to P1 in seconds
1862 - "amplitudes": amplitudes of each phase
1863 - "relamplitudes": amplitudes normalized to amplitude of P1 phase
1864 - "widths": widths of each phase at `width_frac` height
1865 - "areas": areas of each phase
1866 - "relareas": areas of the phases relative to total area
1867 - "zeros": time of zero crossing towards next phase in seconds
1869 unit: string
1870 Unit of phase amplitudes.
1872 Raises
1873 ------
1874 FileNotFoundError:
1875 If `file_path` does not exist.
1877 See Also
1878 --------
1879 save_pulse_phases()
1880 """
1881 data = TableData(file_path)
1882 phases = dict(indices=data['index'].astype(int),
1883 times=data['time']*0.001,
1884 amplitudes=data['amplitude'],
1885 relamplitudes=data['relampl']*0.01,
1886 widths=data['width']*0.001,
1887 areas=data['area']*0.001,
1888 relareas=data['relarea']*0.01,
1889 zeros=data['zeros']*0.001)
1890 return phases, data.unit('amplitude')
1893def save_pulse_gaussians(pulse, unit, idx, basename, **kwargs):
1894 """Save Gaussian phase properties of pulse EOD to file.
1896 Parameters
1897 ----------
1898 pulse: dict
1899 Dictionary with
1901 - "times": phase times in seconds,
1902 - "amplitudes": amplitudes, and
1903 - "stdevs": standard deviations in seconds
1905 of Gaussians fitted to the pulse waveform as returned by
1906 `decompose_pulse()` and `analyze_pulse()`.
1907 unit: string
1908 Unit of the waveform data.
1909 idx: int or None
1910 Index of fish.
1911 basename: string or stream
1912 If string, path and basename of file.
1913 If `basename` does not have an extension,
1914 '-pulsegaussians', the fish index, and a file extension are appended.
1915 If stream, write pulse phases into this stream.
1916 kwargs:
1917 Arguments passed on to `TableData.write()`.
1919 Returns
1920 -------
1921 filename: Path
1922 Path and full name of the written file in case of `basename`
1923 being a string. Otherwise, the file name and extension that
1924 would have been appended to a basename.
1926 See Also
1927 --------
1928 load_pulse_gaussians()
1930 """
1931 if len(pulse) == 0:
1932 return None
1933 td = TableData(pulse,
1934 units=['ms', unit, 'ms'],
1935 formats=['%.3f', '%.5f', '%.3f'])
1936 td['times'] *= 1000
1937 td['stdevs'] *= 1000
1938 fp = ''
1939 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
1940 if not ext:
1941 fp = '-pulsegaussians'
1942 if idx is not None:
1943 fp += f'-{idx}'
1944 return td.write_file_stream(basename, fp, **kwargs)
1947def load_pulse_gaussians(file_path):
1948 """Load Gaussian phase properties of pulse EOD from file.
1950 Parameters
1951 ----------
1952 file_path: string
1953 Path of the file to be loaded.
1955 Returns
1956 -------
1957 pulse: dict
1958 Dictionary with
1960 - "times": phase times in seconds,
1961 - "amplitudes": amplitudes, and
1962 - "stdevs": standard deviations in seconds
1964 of Gaussians fitted to the pulse waveform.
1965 Use the functions provided in thunderfish.fakefish to simulate
1966 pulse fish EODs from this data.
1967 unit: string
1968 Unit of Gaussian amplitudes.
1970 Raises
1971 ------
1972 FileNotFoundError:
1973 If `file_path` does not exist.
1975 See Also
1976 --------
1977 save_pulse_gaussians()
1979 """
1980 data = TableData(file_path)
1981 pulses = data.dict()
1982 pulses['times'] = 0.001*np.array(data['times'])
1983 pulses['amplitudes'] = np.array(data['amplitudes'])
1984 pulses['stdevs'] = 0.001*np.array(data['stdevs'])
1985 return pulses, data.unit('amplitudes')
1988def save_pulse_times(pulse_times, idx, basename, **kwargs):
1989 """Save times of pulse EOD to file.
1991 Parameters
1992 ----------
1993 pulse_times: dict or array of floats
1994 Times of EOD pulses. Either as array of times or
1995 `props['peaktimes']` or `props['times']` as returned by
1996 `analyze_pulse()`.
1997 idx: int or None
1998 Index of fish.
1999 basename: string or stream
2000 If string, path and basename of file.
2001 If `basename` does not have an extension,
2002 '-pulsetimes', the fish index, and a file extension are appended.
2003 If stream, write pulse times into this stream.
2004 kwargs:
2005 Arguments passed on to `TableData.write()`.
2007 Returns
2008 -------
2009 filename: Path
2010 Path and full name of the written file in case of `basename`
2011 being a string. Otherwise, the file name and extension that
2012 would have been appended to a basename.
2014 See Also
2015 --------
2016 load_pulse_times()
2017 """
2018 if isinstance(pulse_times, dict):
2019 props = pulse_times
2020 pulse_times = props.get('times', [])
2021 pulse_times = props.get('peaktimes', pulse_times)
2022 if len(pulse_times) == 0:
2023 return None
2024 td = TableData()
2025 td.append('time', 's', '%.4f', value=pulse_times)
2026 fp = ''
2027 ext = Path(basename).suffix if not hasattr(basename, 'write') else ''
2028 if not ext:
2029 fp = '-pulsetimes'
2030 if idx is not None:
2031 fp += f'-{idx}'
2032 return td.write_file_stream(basename, fp, **kwargs)
2035def load_pulse_times(file_path):
2036 """Load times of pulse EOD from file.
2038 Parameters
2039 ----------
2040 file_path: string
2041 Path of the file to be loaded.
2043 Returns
2044 -------
2045 pulse_times: array of floats
2046 Times of pulse EODs in seconds.
2048 Raises
2049 ------
2050 FileNotFoundError:
2051 If `file_path` does not exist.
2053 See Also
2054 --------
2055 save_pulse_times()
2056 """
2057 data = TableData(file_path)
2058 pulse_times = data.array()[:, 0]
2059 return pulse_times
2062def main():
2063 import matplotlib.pyplot as plt
2064 from thunderlab.eventdetection import snippets
2065 from .fakefish import pulsefish_eods, export_pulsefish
2066 from .eodanalysis import plot_eod_waveform
2068 print('Analysis of pulse-type EOD waveforms.')
2070 deltaf = 2 # frequency resolution
2071 eodf = 80 # EOD frequency
2073 # data:
2074 rate = 100_000
2075 data = pulsefish_eods('Triphasic', eodf, rate, 5.0,
2076 jitter_cv=0.01, noise_std=0.01)
2077 unit = 'mV'
2078 eod_idx, _ = detect_peaks(data, 1.0)
2079 eod_times = eod_idx/rate
2081 # mean EOD:
2082 ns = int(0.003*rate)
2083 snips = snippets(data, eod_idx, start=-ns, stop=ns)
2084 mean_eod = np.mean(snips, 0)
2086 # analyse EOD:
2087 mean_eod, props, peaks, pulses, energy = \
2088 analyze_pulse(mean_eod, rate, eod_times,
2089 freq_resolution=deltaf, fit_gaussians=False)
2091 # write out as python code:
2092 export_pulsefish(pulses, '', 'Triphasic spec')
2094 # full spectrum:
2095 freq, power = pulsetrain_spectrum(eod_times, mean_eod, None,
2096 freq_resolution=deltaf, fade_frac=0.05)
2097 dfreq, dpower = psd(data, rate, freq_resolution=deltaf)
2099 # plot:
2100 fig, axs = plt.subplots(1, 2, layout='constrained')
2101 plot_eod_waveform(axs[0], mean_eod, props, peaks, unit=unit)
2102 axs[0].set_title(f'pulse fish: EODf = {props["EODf"]:.1f} Hz')
2103 ref = plot_pulse_spectrum(axs[1], energy, props)
2104 fac = 6300/deltaf # conversion power to energy spectrum (independent of recording duration and sampling rate, dependent on EOD frequency)
2105 print(f'FFT window={1/freq[1]:.2f}s, number of pulses={len(eod_times)}')
2106 axs[1].plot(freq, decibel(power, ref*fac), 'k')
2107 axs[1].plot(dfreq, decibel(dpower, ref*fac), 'gray')
2108 plt.show()
2111if __name__ == '__main__':
2112 main()