Coverage for src/thunderfish/eodanalysis.py: 32%
1310 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 16:21 +0000
1"""
2Analysis of EOD waveforms.
4## EOD waveform analysis
6- `eod_waveform()`: compute an averaged EOD waveform.
7- `analyze_wave()`: analyze the EOD waveform of a wave fish.
8- `analyze_pulse()`: analyze the EOD waveform of a pulse fish.
9- `adjust_eodf()`: adjust EOD frequencies to a standard temperature.
11## Similarity of EOD waveforms
13- `wave_similarity()`: root-mean squared difference between two wave fish EODs.
14- `pulse_similarity()`: root-mean squared difference between two pulse fish EODs.
15- `load_species_waveforms()`: load template EOD waveforms for species matching.
17## Quality assessment
19- `clipped_fraction()`: compute fraction of clipped EOD waveform snippets.
20- `wave_quality()`: asses quality of EOD waveform of a wave fish.
21- `pulse_quality()`: asses quality of EOD waveform of a pulse fish.
23## Visualization
25- `plot_eod_recording()`: plot a zoomed in range of the recorded trace.
26- `plot_pulse_eods()`: mark pulse EODs in a plot of an EOD recording.
27- `plot_eod_snippets()`: plot a few EOD waveform snippets.
28- `plot_eod_waveform()`: plot and annotate the averaged EOD-waveform with standard error.
29- `plot_wave_spectrum()`: plot and annotate spectrum of wave EODs.
30- `plot_pulse_spectrum()`: plot and annotate spectrum of single pulse EOD.
32## Storage
34- `save_eod_waveform()`: save mean EOD waveform to file.
35- `load_eod_waveform()`: load EOD waveform from file.
36- `save_wave_eodfs()`: save frequencies of wave EODs to file.
37- `load_wave_eodfs()`: load frequencies of wave EODs from file.
38- `save_wave_fish()`: save properties of wave EODs to file.
39- `load_wave_fish()`: load properties of wave EODs from file.
40- `save_pulse_fish()`: save properties of pulse EODs to file.
41- `load_pulse_fish()`: load properties of pulse EODs from file.
42- `save_wave_spectrum()`: save amplitude and phase spectrum of wave EOD to file.
43- `load_wave_spectrum()`: load amplitude and phase spectrum of wave EOD from file.
44- `save_pulse_spectrum()`: save power spectrum of pulse EOD to file.
45- `load_pulse_spectrum()`: load power spectrum of pulse EOD from file.
46- `save_pulse_peaks()`: save peak properties of pulse EOD to file.
47- `load_pulse_peaks()`: load peak properties of pulse EOD from file.
48- `save_pulse_times()`: save times of pulse EOD to file.
49- `load_pulse_times()`: load times of pulse EOD from file.
50- `parse_filename()`: parse components of an EOD analysis file name.
51- `save_analysis(): save EOD analysis results to files.
52- `load_analysis()`: load EOD analysis files.
53- `load_recording()`: load recording.
55## Fit functions
57- `fourier_series()`: Fourier series of sine waves with amplitudes and phases.
58- `exp_decay()`: exponential decay.
60## Filter functions
62- `unfilter()`: apply inverse low-pass filter on data.
64## Configuration
66- `add_eod_analysis_config()`: add parameters for EOD analysis functions to configuration.
67- `eod_waveform_args()`: retrieve parameters for `eod_waveform()` from configuration.
68- `analyze_wave_args()`: retrieve parameters for `analyze_wave()` from configuration.
69- `analyze_pulse_args()`: retrieve parameters for `analyze_pulse()` from configuration.
70- `add_species_config()`: add parameters needed for assigning EOD waveforms to species.
71- `add_eod_quality_config()`: add parameters needed for assesing the quality of an EOD waveform.
72- `wave_quality_args()`: retrieve parameters for `wave_quality()` from configuration.
73- `pulse_quality_args()`: retrieve parameters for `pulse_quality()` from configuration.
74"""
76import os
77import io
78import glob
79import zipfile
80import numpy as np
81from scipy.optimize import curve_fit
82import matplotlib.patches as mpatches
83import matplotlib.pyplot as plt
84from thunderlab.eventdetection import percentile_threshold, detect_peaks, snippets, peak_width
85from thunderlab.eventdetection import threshold_crossings, threshold_crossing_times, merge_events
86from thunderlab.powerspectrum import next_power_of_two, nfft, decibel
87from thunderlab.tabledata import TableData
88from thunderlab.dataloader import load_data
89from .harmonics import fundamental_freqs_and_power
92def eod_waveform(data, samplerate, eod_times, win_fac=2.0, min_win=0.01,
93 min_sem=False, max_eods=None, unfilter_cutoff=0.0):
94 """Detect EODs in the given data, extract data snippets around each EOD,
95 and compute a mean waveform with standard error.
97 Retrieving the EOD waveform of a wave fish works under the following
98 conditions: (i) at a signal-to-noise ratio \\(SNR = P_s/P_n\\),
99 i.e. the power \\(P_s\\) of the EOD of interest relative to the
100 largest other EOD \\(P_n\\), we need to average over at least \\(n >
101 (SNR \\cdot c_s^2)^{-1}\\) snippets to bring the standard error of the
102 averaged EOD waveform down to \\(c_s\\) relative to its
103 amplitude. For a s.e.m. less than 5% ( \\(c_s=0.05\\) ) and an SNR of
104 -10dB (the signal is 10 times smaller than the noise, \\(SNR=0.1\\) ) we
105 get \\(n > 0.00025^{-1} = 4000\\) data snippets - a recording a
106 couple of seconds long. (ii) Very important for wave fish is that
107 they keep their frequency constant. Slight changes in the EOD
108 frequency will corrupt the average waveform. If the period of the
109 waveform changes by \\(c_f=\\Delta T/T\\), then after \\(n =
110 1/c_f\\) periods moved the modified waveform through a whole period.
111 This is in the range of hundreds or thousands waveforms.
113 NOTE: we need to take into account a possible error in the estimate
114 of the EOD period. This will limit the maximum number of snippets to
115 be averaged.
117 If `min_sem` is set, the algorithm checks for a global minimum of
118 the s.e.m. as a function of snippet number. If there is one then
119 the average is computed for this number of snippets, otherwise all
120 snippets are taken from the provided data segment. Note that this
121 check only works for the strongest EOD in a recording. For weaker
122 EOD the s.e.m. always decays with snippet number (empirical
123 observation).
125 TODO: use power spectra to check for changes in EOD frequency!
127 Parameters
128 ----------
129 data: 1-D array of float
130 The data to be analysed.
131 samplerate: float
132 Sampling rate of the data in Hertz.
133 eod_times: 1-D array of float
134 Array of EOD times in seconds over which the waveform should be
135 averaged.
136 WARNING: The first data point must be at time zero!
137 win_fac: float
138 The snippet size is the EOD period times `win_fac`. The EOD period
139 is determined as the minimum interval between EOD times.
140 min_win: float
141 The minimum size of the snippets in seconds.
142 min_sem: bool
143 If set, check for minimum in s.e.m. to set the maximum numbers
144 of EODs to be used for computing the average waveform.
145 max_eods: int or None
146 Maximum number of EODs to be used for averaging.
147 unfilter_cutoff: float
148 If not zero, the cutoff frequency for an inverse high-pass filter
149 applied to the mean EOD waveform.
151 Returns
152 -------
153 mean_eod: 2-D array
154 Average of the EOD snippets. First column is time in seconds,
155 second column the mean eod, third column the standard error.
156 eod_times: 1-D array
157 Times of EOD peaks in seconds that have been actually used to calculate the
158 averaged EOD waveform.
159 """
160 # indices of EOD times:
161 eod_idx = np.round(eod_times * samplerate).astype(int)
163 # window size:
164 period = np.min(np.diff(eod_times))
165 win = 0.5*win_fac*period
166 if 2*win < min_win:
167 win = 0.5*min_win
168 win_inx = int(win * samplerate)
170 # extract snippets:
171 eod_times = eod_times[(eod_idx >= win_inx) & (eod_idx < len(data)-win_inx)]
172 eod_idx = eod_idx[(eod_idx >= win_inx) & (eod_idx < len(data)-win_inx)]
173 if max_eods and max_eods > 0 and len(eod_idx) > max_eods:
174 dn = (len(eod_idx) - max_eods)//2
175 eod_times = eod_times[dn:dn+max_eods]
176 eod_idx = eod_idx[dn:dn+max_eods]
177 eod_snippets = snippets(data, eod_idx, -win_inx, win_inx)
178 if len(eod_snippets) == 0:
179 return np.zeros((0, 3)), eod_times
181 # optimal number of snippets:
182 step = 10
183 if min_sem and len(eod_snippets) > step:
184 sems = [np.mean(np.std(eod_snippets[:k], axis=0, ddof=1)/np.sqrt(k))
185 for k in range(step, len(eod_snippets), step)]
186 idx = np.argmin(sems)
187 # there is a local minimum:
188 if idx > 0 and idx < len(sems)-1:
189 maxn = step*(idx+1)
190 eod_snippets = eod_snippets[:maxn]
191 eod_times = eod_times[:maxn]
193 # mean and std of snippets:
194 mean_eod = np.zeros((len(eod_snippets[0]), 3))
195 mean_eod[:,1] = np.mean(eod_snippets, axis=0)
196 if len(eod_snippets) > 1:
197 mean_eod[:,2] = np.std(eod_snippets, axis=0, ddof=1)/np.sqrt(len(eod_snippets))
199 # apply inverse filter:
200 if unfilter_cutoff and unfilter_cutoff > 0.0:
201 unfilter(mean_eod[:,1], samplerate, unfilter_cutoff)
203 # time axis:
204 mean_eod[:,0] = (np.arange(len(mean_eod)) - win_inx) / samplerate
206 return mean_eod, eod_times
209def unfilter(data, samplerate, cutoff):
210 """Apply inverse high-pass filter on data.
212 Assumes high-pass filter
213 \\[ \\tau \\dot y = -y + \\tau \\dot x \\]
214 has been applied on the original data \\(x\\), where
215 \\(\\tau=(2\\pi f_{cutoff})^{-1}\\) is the time constant of the
216 filter. To recover \\(x\\) the ODE
217 \\[ \\tau \\dot x = y + \\tau \\dot y \\]
218 is applied on the filtered data \\(y\\).
220 Parameters
221 ----------
222 data: ndarray
223 High-pass filtered original data.
224 samplerate: float
225 Sampling rate of `data` in Hertz.
226 cutoff: float
227 Cutoff frequency \\(f_{cutoff}\\) of the high-pass filter in Hertz.
229 Returns
230 -------
231 data: ndarray
232 Recovered original data.
233 """
234 tau = 0.5/np.pi/cutoff
235 fac = tau*samplerate
236 data -= np.mean(data)
237 d0 = data[0]
238 x = d0
239 for k in range(len(data)):
240 d1 = data[k]
241 x += (d1 - d0) + d0/fac
242 data[k] = x
243 d0 = d1
244 return data
247def fourier_series(t, freq, *ap):
248 """Fourier series of sine waves with amplitudes and phases.
250 x(t) = sum_{i=0}^n ap[2*i]*sin(2 pi (i+1) freq t + ap[2*i+1])
252 Parameters
253 ----------
254 t: float or array
255 Time.
256 freq: float
257 Fundamental frequency.
258 *ap: list of floats
259 The amplitudes and phases (in rad) of the fundamental and harmonics.
261 Returns
262 -------
263 x: float or array
264 The Fourier series evaluated at times `t`.
265 """
266 omega = 2.0*np.pi*freq
267 x = 0.0
268 for i, (a, p) in enumerate(zip(ap[0:-1:2], ap[1::2])):
269 x += a*np.sin((i+1)*omega*t+p)
270 return x
273def analyze_wave(eod, freq, n_harm=10, power_n_harmonics=0,
274 n_harmonics=3, flip_wave='none'):
275 """Analyze the EOD waveform of a wave fish.
277 Parameters
278 ----------
279 eod: 2-D array
280 The eod waveform. First column is time in seconds, second
281 column the EOD waveform, third column, if present, is the
282 standard error of the EOD waveform, Further columns are
283 optional but not used.
284 freq: float or 2-D array
285 The frequency of the EOD or the list of harmonics (rows) with
286 frequency and peak height (columns) as returned from
287 `harmonics.harmonic_groups()`.
288 n_harm: int
289 Maximum number of harmonics used for the Fourier decomposition.
290 power_n_harmonics: int
291 Sum over the first `power_n_harmonics` harmonics for computing
292 the total power. If 0 sum over all harmonics.
293 n_harmonics: int
294 The maximum power of higher harmonics is computed from
295 harmonics higher than the maximum harmonics within the first
296 three harmonics plus `n_harmonics`.
297 flip_wave: 'auto', 'none', 'flip'
298 - 'auto' flip waveform such that the larger extremum is positive.
299 - 'flip' flip waveform.
300 - 'none' do not flip waveform.
302 Returns
303 -------
304 meod: 2-D array of floats
305 The eod waveform. First column is time in seconds, second
306 column the eod waveform. Further columns are kept from the
307 input `eod`. And a column is added with the fit of the fourier
308 series to the waveform.
309 props: dict
310 A dictionary with properties of the analyzed EOD waveform.
312 - type: set to 'wave'.
313 - EODf: is set to the EOD fundamental frequency.
314 - p-p-amplitude: peak-to-peak amplitude of the Fourier fit.
315 - flipped: True if the waveform was flipped.
316 - amplitude: amplitude factor of the Fourier fit.
317 - noise: root-mean squared standard error mean of the averaged
318 EOD waveform relative to the p-p amplitude.
319 - rmserror: root-mean-square error between Fourier-fit and
320 EOD waveform relative to the p-p amplitude. If larger than
321 about 0.05 the data are bad.
322 - ncrossings: number of zero crossings per period
323 - peakwidth: width of the peak at the averaged amplitude relative
324 to EOD period.
325 - troughwidth: width of the trough at the averaged amplitude
326 relative to EOD period.
327 - leftpeak: time from positive zero crossing to peak relative
328 to EOD period.
329 - rightpeak: time from peak to negative zero crossing relative to
330 EOD period.
331 - lefttrough: time from negative zero crossing to trough relative to
332 EOD period.
333 - righttrough: time from trough to positive zero crossing relative to
334 EOD period.
335 - p-p-distance: time between peak and trough relative to EOD period.
336 - relpeakampl: amplitude of peak or trough, whichever is larger, relative to p-p amplitude.
337 - power: summed power of all harmonics of the extracted EOD waveform
338 in decibel relative to one.
339 - datapower: summed power of all harmonics of the original data in
340 decibel relative to one. Only if `freq` is a list of harmonics.
341 - thd: total harmonic distortion, i.e. square root of sum of
342 amplitudes squared of harmonics relative to amplitude
343 of fundamental.
344 - dbdiff: smoothness of power spectrum as standard deviation of
345 differences in decibel power.
346 - maxdb: maximum power of higher harmonics relative to peak power
347 in decibel.
349 spec_data: 2-D array of floats
350 First size columns are from the spectrum of the extracted
351 waveform. First column is the index of the harmonics, second
352 column its frequency, third column its amplitude, fourth
353 column its amplitude relative to the fundamental, fifth column
354 is power of harmonics relative to fundamental in decibel, and
355 sixth column the phase shift relative to the fundamental.
356 If `freq` is a list of harmonics, a seventh column is added to
357 `spec_data` that contains the powers of the harmonics from the
358 original power spectrum of the raw data. Rows are the
359 harmonics, first row is the fundamental frequency with index
360 0, relative amplitude of one, relative power of 0dB, and phase
361 shift of zero.
362 error_str: string
363 If fitting of the fourier series failed,
364 this is reported in this string.
366 Raises
367 ------
368 IndexError:
369 EOD data is less than one period long.
370 """
371 error_str = ''
373 freq0 = freq
374 if hasattr(freq, 'shape'):
375 freq0 = freq[0][0]
377 # storage:
378 meod = np.zeros((eod.shape[0], eod.shape[1]+1))
379 meod[:,:-1] = eod
381 # subtract mean and flip:
382 period = 1.0/freq0
383 pinx = int(np.ceil(period/(meod[1,0]-meod[0,0])))
384 maxn = (len(meod)//pinx)*pinx
385 if maxn < pinx: maxn = len(meod)
386 offs = (len(meod) - maxn)//2
387 meod[:,1] -= np.mean(meod[offs:offs+pinx,1])
388 flipped = False
389 if 'flip' in flip_wave or ('auto' in flip_wave and -np.min(meod[:,1]) > np.max(meod[:,1])):
390 meod[:,1] = -meod[:,1]
391 flipped = True
393 # move peak of waveform to zero:
394 offs = len(meod)//4
395 maxinx = offs+np.argmax(meod[offs:3*offs,1])
396 meod[:,0] -= meod[maxinx,0]
398 # indices of exactly one or two periods around peak:
399 if len(meod) < pinx:
400 raise IndexError('data need to contain at least one EOD period')
401 if len(meod) >= 2*pinx:
402 i0 = maxinx - pinx if maxinx >= pinx else 0
403 i1 = i0 + 2*pinx
404 if i1 > len(meod):
405 i1 = len(meod)
406 i0 = i1 - 2*pinx
407 else:
408 i0 = maxinx - pinx//2 if maxinx >= pinx//2 else 0
409 i1 = i0 + pinx
411 # subtract mean:
412 meod[:,1] -= np.mean(meod[i0:i1,1])
414 # zero crossings:
415 ui, di = threshold_crossings(meod[:,1], 0.0)
416 ut, dt = threshold_crossing_times(meod[:,0], meod[:,1], 0.0, ui, di)
417 ut, dt = merge_events(ut, dt, 0.02/freq0)
418 ncrossings = int(np.round((len(ut) + len(dt))/(meod[-1,0]-meod[0,0])/freq0))
419 if np.any(ut<0.0):
420 up_time = ut[ut<0.0][-1]
421 else:
422 up_time = 0.0
423 error_str += '%.1f Hz wave fish: no upward zero crossing. ' % freq0
424 if np.any(dt>0.0):
425 down_time = dt[dt>0.0][0]
426 else:
427 down_time = 0.0
428 error_str += '%.1f Hz wave fish: no downward zero crossing. ' % freq0
429 peak_width = down_time - up_time
430 trough_width = period - peak_width
431 peak_time = 0.0
432 trough_time = meod[maxinx+np.argmin(meod[maxinx:maxinx+pinx,1]),0]
433 phase1 = peak_time - up_time
434 phase2 = down_time - peak_time
435 phase3 = trough_time - down_time
436 phase4 = up_time + period - trough_time
437 distance = trough_time - peak_time
439 # fit fourier series:
440 ampl = 0.5*(np.max(meod[:,1])-np.min(meod[:,1]))
441 while n_harm > 1:
442 params = [freq0]
443 for i in range(1, n_harm+1):
444 params.extend([ampl/i, 0.0])
445 try:
446 popt, pcov = curve_fit(fourier_series, meod[i0:i1,0],
447 meod[i0:i1,1], params, maxfev=2000)
448 break
449 except (RuntimeError, TypeError):
450 error_str += '%.1f Hz wave fish: fit of fourier series failed for %d harmonics. ' % (freq0, n_harm)
451 n_harm //= 2
452 # store fourier fit:
453 meod[:,-1] = fourier_series(meod[:,0], *popt)
454 # make all amplitudes positive:
455 for i in range(n_harm):
456 if popt[i*2+1] < 0.0:
457 popt[i*2+1] *= -1.0
458 popt[i*2+2] += np.pi
459 # phases relative to fundamental:
460 # phi0 = 2*pi*f0*dt <=> dt = phi0/(2*pi*f0)
461 # phik = 2*pi*i*f0*dt = i*phi0
462 phi0 = popt[2]
463 for i in range(n_harm):
464 popt[i*2+2] -= (i + 1)*phi0
465 # all phases in the range -pi to pi:
466 popt[i*2+2] %= 2.0*np.pi
467 if popt[i*2+2] > np.pi:
468 popt[i*2+2] -= 2.0*np.pi
469 # shift time axis:
470 # meod[:,0] += phi0/2/np.pi/freq0
471 # store fourier spectrum:
472 if hasattr(freq, 'shape'):
473 n = n_harm
474 n += np.sum(freq[:,0] > (n_harm+0.5)*freq[0,0])
475 spec_data = np.zeros((n, 7))
476 spec_data[:,:] = np.nan
477 k = 0
478 for i in range(n_harm):
479 while k < len(freq) and freq[k,0] < (i+0.5)*freq0:
480 k += 1
481 if k >= len(freq):
482 break
483 if freq[k,0] < (i+1.5)*freq0:
484 spec_data[i,6] = freq[k,1]
485 k += 1
486 for i in range(n_harm, n):
487 if k >= len(freq):
488 break
489 spec_data[i,:2] = [np.round(freq[k,0]/freq0)-1, freq[k,0]]
490 spec_data[i,6] = freq[k,1]
491 k += 1
492 else:
493 spec_data = np.zeros((n_harm, 6))
494 for i in range(n_harm):
495 spec_data[i,:6] = [i, (i+1)*freq0, popt[i*2+1], popt[i*2+1]/popt[1],
496 decibel((popt[i*2+1]/popt[1])**2.0), popt[i*2+2]]
497 # smoothness of power spectrum:
498 db_powers = decibel(spec_data[:n_harm,2]**2)
499 db_diff = np.std(np.diff(db_powers))
500 # maximum relative power of higher harmonics:
501 p_max = np.argmax(db_powers[:3])
502 db_powers -= db_powers[p_max]
503 if len(db_powers[p_max+n_harmonics:]) == 0:
504 max_harmonics_power = -100.0
505 else:
506 max_harmonics_power = np.max(db_powers[p_max+n_harmonics:])
507 # total harmonic distortion:
508 thd = np.sqrt(np.nansum(spec_data[1:,3]))
510 # peak-to-peak and trough amplitudes:
511 ppampl = np.max(meod[i0:i1,1]) - np.min(meod[i0:i1,1])
512 relpeakampl = min(np.max(meod[i0:i1,1]), np.abs(np.min(meod[i0:i1,1])))/ppampl
514 # variance and fit error:
515 rmssem = np.sqrt(np.mean(meod[i0:i1,2]**2.0))/ppampl if eod.shape[1] > 2 else None
516 rmserror = np.sqrt(np.mean((meod[i0:i1,1] - meod[i0:i1,-1])**2.0))/ppampl
518 # store results:
519 props = {}
520 props['type'] = 'wave'
521 props['EODf'] = freq0
522 props['p-p-amplitude'] = ppampl
523 props['flipped'] = flipped
524 props['amplitude'] = 0.5*ppampl # remove it
525 props['rmserror'] = rmserror
526 if rmssem:
527 props['noise'] = rmssem
528 props['ncrossings'] = ncrossings
529 props['peakwidth'] = peak_width/period
530 props['troughwidth'] = trough_width/period
531 props['leftpeak'] = phase1/period
532 props['rightpeak'] = phase2/period
533 props['lefttrough'] = phase3/period
534 props['righttrough'] = phase4/period
535 props['p-p-distance'] = distance/period
536 props['relpeakampl'] = relpeakampl
537 pnh = power_n_harmonics if power_n_harmonics > 0 else n_harm
538 pnh = min(n_harm, pnh)
539 props['power'] = decibel(np.sum(spec_data[:pnh,2]**2.0))
540 if hasattr(freq, 'shape'):
541 props['datapower'] = decibel(np.sum(freq[:pnh,1]))
542 props['thd'] = thd
543 props['dbdiff'] = db_diff
544 props['maxdb'] = max_harmonics_power
546 return meod, props, spec_data, error_str
549def exp_decay(t, tau, ampl, offs):
550 """Exponential decay function.
552 x(t) = ampl*exp(-t/tau) + offs
554 Parameters
555 ----------
556 t: float or array
557 Time.
558 tau: float
559 Time constant of exponential decay.
560 ampl: float
561 Amplitude of exponential decay, i.e. initial value minus
562 steady-state value.
563 offs: float
564 Steady-state value.
566 Returns
567 -------
568 x: float or array
569 The exponential decay evaluated at times `t`.
571 """
572 return offs + ampl*np.exp(-t/tau)
575def analyze_pulse(eod, eod_times=None, min_pulse_win=0.001,
576 peak_thresh_fac=0.01, min_dist=50.0e-6,
577 width_frac=0.5, fit_frac = 0.5, freq_resolution=1.0,
578 flip_pulse='none', ipi_cv_thresh=0.5,
579 ipi_percentile=30.0):
580 """Analyze the EOD waveform of a pulse fish.
582 Parameters
583 ----------
584 eod: 2-D array
585 The eod waveform. First column is time in seconds, second
586 column the EOD waveform, third column, if present, is the
587 standard error of the EOD waveform, Further columns are
588 optional but not used.
589 eod_times: 1-D array or None
590 List of times of detected EOD peaks.
591 min_pulse_win: float
592 The minimum size of cut-out EOD waveform.
593 peak_thresh_fac: float
594 Set the threshold for peak detection to the maximum pulse
595 amplitude times this factor.
596 min_dist: float
597 Minimum distance between peak and troughs of the pulse.
598 width_frac: float
599 The width of a peak is measured at this fraction of a peak's
600 height (0-1).
601 fit_frac: float or None
602 An exponential is fitted to the tail of the last peak/trough
603 starting where the waveform falls below this fraction of the
604 peak's height (0-1).
605 freq_resolution: float
606 The frequency resolution of the power spectrum of the single pulse.
607 flip_pulse: 'auto', 'none', 'flip'
608 - 'auto' flip waveform such that the first large extremum is positive.
609 - 'flip' flip waveform.
610 - 'none' do not flip waveform.
611 ipi_cv_thresh: float
612 If the coefficient of variation of the interpulse intervals
613 are smaller than this threshold, then the EOD frequency is
614 computed as the inverse of the mean of all interpulse
615 intervals. Otherwise only intervals smaller than a certain
616 quantile are used.
617 ipi_percentile: float
618 When computing the EOD frequency, period, mean and standard
619 deviation of interpulse intervals from a subset of the
620 interpulse intervals, only intervals smaller than this
621 percentile (between 0 and 100) are used.
623 Returns
624 -------
625 meod: 2-D array of floats
626 The eod waveform. First column is time in seconds,
627 second column the eod waveform.
628 Further columns are kept from the input `eod`.
629 As a last column the fit to the tail of the last peak is appended.
630 props: dict
631 A dictionary with properties of the analyzed EOD waveform.
633 - type: set to 'pulse'.
634 - EODf: the inverse of the median interval between `eod_times`,
635 if provided.
636 - period: the median interval between `eod_times`, if provided.
637 - IPI-mean: the mean interval between `eod_times`, if provided.
638 - IPI-std: the standard deviation of the intervals between
639 `eod_times`, if provided.
640 - max-ampl: the amplitude of the largest positive peak (P1).
641 - min-ampl: the amplitude of the largest negative peak (P2).
642 - p-p-amplitude: peak-to-peak amplitude of the EOD waveform.
643 - noise: root-mean squared standard error mean of the averaged
644 EOD waveform relative to the p-p amplitude.
645 - tstart: time in seconds where the pulse starts,
646 i.e. crosses the threshold for the first time.
647 - tend: time in seconds where the pulse ends,
648 i.e. crosses the threshold for the last time.
649 - width: total width of the pulse in seconds (tend-tstart).
650 - P2-P1-dist: distance between P2 and P1 in seconds.
651 - tau: time constant of exponential decay of pulse tail in seconds.
652 - firstpeak: index of the first peak in the pulse (i.e. -1 for P-1)
653 - lastpeak: index of the last peak in the pulse (i.e. 3 for P3)
654 - peakfreq: frequency at peak power of the single pulse spectrum
655 in Hertz.
656 - peakpower: peak power of the single pulse spectrum in decibel.
657 - poweratt5: how much the average power below 5 Hz is attenuated
658 relative to the peak power in decibel.
659 - poweratt50: how much the average power below 5 Hz is attenuated
660 relative to the peak power in decibel.
661 - lowcutoff: frequency at which the power reached half of the
662 peak power relative to the initial power in Hertz.
663 - flipped: True if the waveform was flipped.
664 - n: number of pulses analyzed (i.e. `len(eod_times)`), if provided.
665 - times: the times of the detected EOD pulses (i.e. `eod_times`),
666 if provided.
668 Empty if waveform is not a pulse EOD.
669 peaks: 2-D array
670 For each peak and trough (rows) of the EOD waveform
671 5 columns: the peak index (1 is P1, i.e. the largest positive peak),
672 time relative to largest positive peak, amplitude,
673 amplitude normalized to largest postive peak,
674 and width of peak/trough at half height.
675 Empty if waveform is not a pulse EOD.
676 power: 2-D array
677 The power spectrum of a single pulse. First column are the
678 frequencies, second column the power in x^2/Hz such that the
679 integral equals the variance. Empty if waveform is not a
680 pulse EOD.
682 """
683 # storage:
684 meod = np.zeros((eod.shape[0], eod.shape[1]+1))
685 meod[:,:eod.shape[1]] = eod
686 meod[:,-1] = np.nan
687 toffs = 0
689 # cut out stable estimate if standard deviation is available:
690 if eod.shape[1] > 2 and np.max(meod[:,2]) > 3*np.min(meod[:,2]):
691 n = len(meod)
692 idx0 = np.argmax(np.abs(meod[n//10:9*n//10,1])) + n//10
693 toffs += meod[idx0,0]
694 meod[:,0] -= meod[idx0,0]
695 # minimum in standard deviation:
696 lstd_idx = np.argmin(meod[:idx0-5,2])
697 rstd_idx = np.argmin(meod[idx0+5:,2]) + idx0
698 # central, left, and right maximum of standard deviation:
699 max_std = np.max(meod[lstd_idx:rstd_idx,2])
700 l_std = np.max(meod[:len(meod)//4,2])
701 r_std = np.max(meod[len(meod)*3//4:,2])
702 lidx = 0
703 ridx = len(meod)
704 if l_std > max_std and lstd_idx > lidx:
705 lidx = lstd_idx - np.argmax(meod[lstd_idx:0:-1,2] >= 0.25*l_std + 0.75*meod[lstd_idx,2])
706 if r_std > max_std and rstd_idx < ridx:
707 ridx = rstd_idx + np.argmax(meod[rstd_idx:,2] >= 0.25*r_std + 0.75*meod[rstd_idx,2])
708 #plt.plot(meod[:,0], meod[:,1])
709 #plt.plot(meod[:,0], meod[:,2], '-r')
710 #plt.plot([meod[lidx,0], meod[lidx,0]], [-0.1, 0.1], '-k')
711 #plt.plot([meod[ridx-1,0], meod[ridx-1,0]], [-0.1, 0.1], '-b')
712 #plt.show()
713 meod = meod[lidx:ridx,:]
715 # subtract mean computed from the ends of the snippet:
716 n = len(meod)//20 if len(meod) >= 20 else 1
717 meod[:,1] -= 0.5*(np.mean(meod[:n,1]) + np.mean(meod[-n:,1]))
719 # largest positive and negative peak:
720 flipped = False
721 max_idx = np.argmax(meod[:,1])
722 max_ampl = np.abs(meod[max_idx,1])
723 min_idx = np.argmin(meod[:,1])
724 min_ampl = np.abs(meod[min_idx,1])
725 amplitude = np.max((max_ampl, min_ampl))
726 if max_ampl > 0.2*amplitude and min_ampl > 0.2*amplitude:
727 # two major peaks:
728 if 'flip' in flip_pulse or ('auto' in flip_pulse and min_idx < max_idx):
729 # flip:
730 meod[:,1] = -meod[:,1]
731 peak_idx = min_idx
732 min_idx = max_idx
733 max_idx = peak_idx
734 flipped = True
735 elif 'flip' in flip_pulse or ('auto' in flip_pulse and min_ampl > 0.2*amplitude):
736 # flip:
737 meod[:,1] = -meod[:,1]
738 peak_idx = min_idx
739 min_idx = max_idx
740 max_idx = peak_idx
741 flipped = True
742 max_ampl = np.abs(meod[max_idx,1])
743 min_ampl = np.abs(meod[min_idx,1])
745 # move peak of waveform to zero:
746 toffs += meod[max_idx,0]
747 meod[:,0] -= meod[max_idx,0]
749 # minimum threshold for peak detection:
750 n = len(meod[:,1])//10 if len(meod) >= 20 else 2
751 thl_max = np.max(meod[:n,1])
752 thl_min = np.min(meod[:n,1])
753 thr_max = np.max(meod[-n:,1])
754 thr_min = np.min(meod[-n:,1])
755 min_thresh = 2*np.max([thl_max, thr_max]) - np.min([thl_min, thr_min])
756 if min_thresh > 0.5*(max_ampl + min_ampl):
757 min_thresh = 0.5*(max_ampl + min_ampl)
758 fit_frac = None
759 # threshold for peak detection:
760 threshold = max_ampl*peak_thresh_fac
761 if threshold < min_thresh:
762 threshold = min_thresh
763 if threshold <= 0.0:
764 return meod, {}, [], []
766 # cut out relevant signal:
767 lidx = np.argmax(np.abs(meod[:,1]) > threshold)
768 ridx = len(meod) - 1 - np.argmax(np.abs(meod[::-1,1]) > threshold)
769 t0 = meod[lidx,0]
770 t1 = meod[ridx,0]
771 width = t1 - t0
772 if width < min_pulse_win:
773 width = min_pulse_win
774 dt = meod[1,0] - meod[0,0]
775 width_idx = int(np.round(width/dt))
776 # expand width:
777 leidx = lidx - width_idx//2
778 if leidx < 0:
779 leidx = 0
780 reidx = ridx + width_idx//2
781 if reidx >= len(meod):
782 reidx = len(meod)
783 meod = meod[leidx:reidx,:]
784 lidx -= leidx
785 ridx -= leidx
786 max_idx -= leidx
787 min_idx -= leidx
788 tau = None
789 dist = 0.0
790 peaks = []
792 # amplitude and variance:
793 ppampl = max_ampl + min_ampl
794 rmssem = np.sqrt(np.mean(meod[:,2]**2.0))/ppampl if eod.shape[1] > 2 else None
796 # find smaller peaks:
797 peak_idx, trough_idx = detect_peaks(meod[:,1], threshold)
799 if len(peak_idx) > 0:
800 # and their width:
801 peak_widths = peak_width(meod[:,0], meod[:,1], peak_idx, trough_idx,
802 peak_frac=width_frac, base='max')
803 trough_widths = peak_width(meod[:,0], -meod[:,1], trough_idx, peak_idx,
804 peak_frac=width_frac, base='max')
805 # combine peaks and troughs:
806 pt_idx = np.concatenate((peak_idx, trough_idx))
807 pt_widths = np.concatenate((peak_widths, trough_widths))
808 pts_idx = np.argsort(pt_idx)
809 peak_list = pt_idx[pts_idx]
810 width_list = pt_widths[pts_idx]
811 # remove multiple peaks that are too close: XXX replace by Dexters function that keeps the maximum peak
812 rmidx = [(k, k+1) for k in np.where(np.diff(meod[peak_list,0]) < min_dist)[0]]
813 # flatten and keep maximum peak:
814 rmidx = np.unique([k for kk in rmidx for k in kk if peak_list[k] != max_idx])
815 # delete:
816 if len(rmidx) > 0:
817 peak_list = np.delete(peak_list, rmidx)
818 width_list = np.delete(width_list, rmidx)
819 if len(peak_list) == 0:
820 return meod, {}, [], []
821 # find P1:
822 p1i = np.argmax(peak_list == max_idx)
823 # truncate peaks to the left: XXX REALLY? WHY?
824 offs = 0 if p1i <= 2 else p1i - 2
825 peak_list = peak_list[offs:]
826 width_list = width_list[offs:]
827 p1i -= offs
828 # store peaks:
829 peaks = np.zeros((len(peak_list), 5))
830 for i, pi in enumerate(peak_list):
831 peaks[i,:] = [i+1-p1i, meod[pi,0], meod[pi,1], meod[pi,1]/max_ampl, width_list[i]]
832 # P2 - P1 distance:
833 dist = peaks[p1i+1,1] - peaks[p1i,1] if p1i+1 < len(peaks) else 0.0
834 # fit exponential to last peak/trough:
835 pi = peak_list[-1]
836 # positive or negative decay:
837 sign = 1.0
838 if np.sum(meod[pi:,1] < -0.5*threshold) > np.sum(meod[pi:,1] > 0.5*threshold):
839 sign = -1.0
840 if sign*meod[pi,1] < 0.0:
841 pi += np.argmax(sign*meod[pi:,1])
842 pi_ampl = np.abs(meod[pi,1])
843 n = len(meod[pi:])
844 # no sufficiently large initial value:
845 if fit_frac and pi_ampl*fit_frac <= 0.5*threshold:
846 fit_frac = False
847 # no sufficiently long decay:
848 if n < 10:
849 fit_frac = False
850 # not decaying towards zero:
851 max_line = pi_ampl - (pi_ampl-threshold)*np.arange(n)/n + 1e-8
852 if np.any(np.abs(meod[pi+2:,1]) > max_line[2:]):
853 fit_frac = False
854 if fit_frac:
855 thresh = meod[pi,1]*fit_frac
856 inx = pi + np.argmax(sign*meod[pi:,1] < sign*thresh)
857 thresh = meod[inx,1]*np.exp(-1.0)
858 tau_inx = np.argmax(sign*meod[inx:,1] < sign*thresh)
859 if tau_inx < 2:
860 tau_inx = 2
861 rridx = inx + 6*tau_inx
862 if rridx > len(meod)-1:
863 tau = None
864 else:
865 tau = meod[inx+tau_inx,0]-meod[inx,0]
866 params = [tau, meod[inx,1]-meod[rridx,1], meod[rridx,1]]
867 try:
868 popt, pcov = curve_fit(exp_decay, meod[inx:rridx,0]-meod[inx,0],
869 meod[inx:rridx,1], params,
870 bounds=([0.0, -np.inf, -np.inf], np.inf))
871 except TypeError:
872 popt, pcov = curve_fit(exp_decay, meod[inx:rridx,0]-meod[inx,0],
873 meod[inx:rridx,1], params)
874 if popt[0] > 1.2*tau:
875 tau_inx = int(np.round(popt[0]/dt))
876 rridx = inx + 6*tau_inx
877 if rridx > len(meod)-1:
878 rridx = len(meod)-1
879 try:
880 popt, pcov = curve_fit(exp_decay, meod[inx:rridx,0]-meod[inx,0],
881 meod[inx:rridx,1], popt,
882 bounds=([0.0, -np.inf, -np.inf], np.inf))
883 except TypeError:
884 popt, pcov = curve_fit(exp_decay, meod[inx:rridx,0]-meod[inx,0],
885 meod[inx:rridx,1], popt)
886 tau = popt[0]
887 meod[inx:rridx,-1] = exp_decay(meod[inx:rridx,0]-meod[inx,0], *popt)
889 # power spectrum of single pulse:
890 samplerate = 1.0/(meod[1,0]-meod[0,0])
891 n_fft = nfft(samplerate, freq_resolution)
893 n0 = max_idx
894 n1 = len(meod)-max_idx
895 n = 2*max(n0, n1)
896 if n_fft < n:
897 n_fft = next_power_of_two(n)
898 data = np.zeros(n_fft)
899 data[n_fft//2-n0:n_fft//2+n1] = meod[:,1]
900 nr = n//4
901 data[n_fft//2-n0:n_fft//2-n0+nr] *= np.arange(nr)/nr
902 data[n_fft//2+n1-nr:n_fft//2+n1] *= np.arange(nr)[::-1]/nr
903 freqs = np.fft.rfftfreq(n_fft, 1.0/samplerate)
904 fourier = np.fft.rfft(data)/n_fft/freqs[1]
905 power = np.abs(fourier)**2.0
906 ppower = np.zeros((len(power), 2))
908 ppower[:,0] = freqs
909 ppower[:,1] = power
910 maxpower = np.max(power)
911 att5 = decibel(np.mean(power[freqs<5.0]), maxpower)
912 att50 = decibel(np.mean(power[freqs<50.0]), maxpower)
913 lowcutoff = freqs[decibel(power, maxpower) > 0.5*att5][0]
915 # analyze pulse timing:
916 if eod_times is not None:
917 inter_pulse_intervals = np.diff(eod_times)
918 ipi_cv = np.std(inter_pulse_intervals)/np.mean(inter_pulse_intervals)
919 if ipi_cv < ipi_cv_thresh:
920 period = np.median(inter_pulse_intervals)
921 ipi_mean = np.mean(inter_pulse_intervals)
922 ipi_std = np.std(inter_pulse_intervals)
923 else:
924 intervals = inter_pulse_intervals[inter_pulse_intervals <
925 np.percentile(inter_pulse_intervals, ipi_percentile)]
926 period = np.median(intervals)
927 ipi_mean = np.mean(intervals)
928 ipi_std = np.std(intervals)
930 # store properties:
931 props = {}
932 props['type'] = 'pulse'
933 if eod_times is not None:
934 props['EODf'] = 1.0/period
935 props['period'] = period
936 props['IPI-mean'] = ipi_mean
937 props['IPI-std'] = ipi_std
938 props['max-ampl'] = max_ampl
939 props['min-ampl'] = min_ampl
940 props['p-p-amplitude'] = ppampl
941 if rmssem:
942 props['noise'] = rmssem
943 props['tstart'] = t0
944 props['tend'] = t1
945 props['width'] = t1-t0
946 props['P2-P1-dist'] = dist
947 if tau:
948 props['tau'] = tau
949 props['firstpeak'] = peaks[0, 0] if len(peaks) > 0 else 1
950 props['lastpeak'] = peaks[-1, 0] if len(peaks) > 0 else 1
951 props['peakfreq'] = freqs[np.argmax(power)]
952 props['peakpower'] = decibel(maxpower)
953 props['poweratt5'] = att5
954 props['poweratt50'] = att50
955 props['lowcutoff'] = lowcutoff
956 props['flipped'] = flipped
957 if eod_times is not None:
958 props['n'] = len(eod_times)
959 props['times'] = eod_times + toffs
961 return meod, props, peaks, ppower
964def adjust_eodf(eodf, temp, temp_adjust=25.0, q10=1.62):
965 """Adjust EOD frequencies to a standard temperature using Q10.
967 Parameters
968 ----------
969 eodf: float or ndarray
970 EOD frequencies.
971 temp: float
972 Temperature in degree celsisus at which EOD frequencies in
973 `eodf` were measured.
974 temp_adjust: float
975 Standard temperature in degree celsisus to which EOD
976 frequencies are adjusted.
977 q10: float
978 Q10 value describing temperature dependence of EOD
979 frequencies. The default of 1.62 is from Dunlap, Smith, Yetka
980 (2000) Brain Behav Evol, measured for Apteronotus
981 lepthorhynchus in the lab.
983 Returns
984 -------
985 eodf_corrected: float or array
986 EOD frequencies adjusted to `temp_adjust` using `q10`.
987 """
988 return eodf*q10**((temp_adjust - temp) / 10.0)
991def load_species_waveforms(species_file='none'):
992 """Load template EOD waveforms for species matching.
994 Parameters
995 ----------
996 species_file: string
997 Name of file containing species definitions. The content of
998 this file is as follows:
1000 - Empty lines and line starting with a hash ('#') are skipped.
1001 - A line with the key-word 'wavefish' marks the beginning of the
1002 table for wave fish.
1003 - A line with the key-word 'pulsefish' marks the beginning of the
1004 table for pulse fish.
1005 - Each line in a species table has three fields,
1006 separated by colons (':'):
1008 1. First field is an abbreviation of the species name.
1009 2. Second field is the filename of the recording containing the
1010 EOD waveform.
1011 3. The optional third field is the EOD frequency of the EOD waveform.
1013 The EOD frequency is used to normalize the time axis of a
1014 wave fish EOD to one EOD period. If it is not specified in
1015 the third field, it is taken from the corresponding
1016 *-wavespectrum-* file, if present. Otherwise the species is
1017 excluded and a warning is issued.
1019 Example file content:
1020 ``` plain
1021 Wavefish
1022 Aptero : F_91009L20-eodwaveform-0.csv : 823Hz
1023 Eigen : C_91008L01-eodwaveform-0.csv
1025 Pulsefish
1026 Gymnotus : pulsefish/gymnotus.csv
1027 Brachy : H_91009L11-eodwaveform-0.csv
1028 ```
1030 Returns
1031 -------
1032 wave_names: list of strings
1033 List of species names of wave-type fish.
1034 wave_eods: list of 2-D arrays
1035 List of EOD waveforms of wave-type fish corresponding to
1036 `wave_names`. First column of a waveform is time in seconds,
1037 second column is the EOD waveform. The waveforms contain
1038 exactly one EOD period.
1039 pulse_names: list of strings
1040 List of species names of pulse-type fish.
1041 pulse_eods: list of 2-D arrays
1042 List of EOD waveforms of pulse-type fish corresponding to
1043 `pulse_names`. First column of a waveform is time in seconds,
1044 second column is the EOD waveform.
1045 """
1046 if len(species_file) == 0 or species_file == 'none' or \
1047 not os.path.isfile(species_file):
1048 return [], [], [], []
1049 wave_names = []
1050 wave_eods = []
1051 pulse_names = []
1052 pulse_eods = []
1053 fish_type = 'wave'
1054 with open(species_file, 'r') as sf:
1055 for line in sf:
1056 line = line.strip()
1057 if len(line) == 0 or line[0] == '#':
1058 continue
1059 if line.lower() == 'wavefish':
1060 fish_type = 'wave'
1061 elif line.lower() == 'pulsefish':
1062 fish_type = 'pulse'
1063 else:
1064 ls = [s.strip() for s in line.split(':')]
1065 if len(ls) < 2:
1066 continue
1067 name = ls[0]
1068 waveform_file = ls[1]
1069 eod = TableData(waveform_file).array()
1070 eod[:,0] *= 0.001
1071 if fish_type == 'wave':
1072 eodf = None
1073 if len(ls) > 2:
1074 eodf = float(ls[2].replace('Hz', '').strip())
1075 else:
1076 spectrum_file = waveform_file.replace('eodwaveform', 'wavespectrum')
1077 try:
1078 wave_spec = TableData(spectrum_file)
1079 eodf = wave_spec[0, 1]
1080 except FileNotFoundError:
1081 pass
1082 if eodf is None:
1083 print('warning: unknown EOD frequency of %s. Skip.' % name)
1084 continue
1085 eod[:,0] *= eodf
1086 wave_names.append(name)
1087 wave_eods.append(eod[:,:2])
1088 elif fish_type == 'pulse':
1089 pulse_names.append(name)
1090 pulse_eods.append(eod[:,:2])
1091 return wave_names, wave_eods, pulse_names, pulse_eods
1094def wave_similarity(eod1, eod2, eod1f=1.0, eod2f=1.0):
1095 """Root-mean squared difference between two wave fish EODs.
1097 Compute the root-mean squared difference between two wave fish
1098 EODs over one period. The better sampled signal is down-sampled to
1099 the worse sampled one. Amplitude are normalized to peak-to-peak
1100 amplitude before computing rms difference. Also compute the rms
1101 difference between the one EOD and the other one inverted in
1102 amplitude. The smaller of the two rms values is returned.
1104 Parameters
1105 ----------
1106 eod1: 2-D array
1107 Time and amplitude of reference EOD.
1108 eod2: 2-D array
1109 Time and amplitude of EOD that is to be compared to `eod1`.
1110 eod1f: float
1111 EOD frequency of `eod1` used to transform the time axis of `eod1`
1112 to multiples of the EOD period. If already normalized to EOD period,
1113 as for example by the `load_species_waveforms()` function, then
1114 set the EOD frequency to one (default).
1115 eod2f: float
1116 EOD frequency of `eod2` used to transform the time axis of `eod2`
1117 to multiples of the EOD period. If already normalized to EOD period,
1118 as for example by the `load_species_waveforms()` function, then
1119 set the EOD frequency to one (default).
1121 Returns
1122 -------
1123 rmse: float
1124 Root-mean-squared difference between the two EOD waveforms relative to
1125 their standard deviation over one period.
1126 """
1127 # copy:
1128 eod1 = np.array(eod1[:,:2])
1129 eod2 = np.array(eod2[:,:2])
1130 # scale to multiples of EOD period:
1131 eod1[:,0] *= eod1f
1132 eod2[:,0] *= eod2f
1133 # make eod1 the waveform with less samples per period:
1134 n1 = int(1.0/(eod1[1,0]-eod1[0,0]))
1135 n2 = int(1.0/(eod2[1,0]-eod2[0,0]))
1136 if n1 > n2:
1137 eod1, eod2 = eod2, eod1
1138 n1, n2 = n2, n1
1139 # one period around time zero:
1140 i0 = np.argmin(np.abs(eod1[:,0]))
1141 k0 = i0-n1//2
1142 if k0 < 0:
1143 k0 = 0
1144 k1 = k0 + n1 + 1
1145 if k1 >= len(eod1):
1146 k1 = len(eod1)
1147 # cut out one period from the reference EOD around maximum:
1148 i = k0 + np.argmax(eod1[k0:k1,1])
1149 k0 = i-n1//2
1150 if k0 < 0:
1151 k0 = 0
1152 k1 = k0 + n1 + 1
1153 if k1 >= len(eod1):
1154 k1 = len(eod1)
1155 eod1 = eod1[k0:k1,:]
1156 # normalize amplitudes of first EOD:
1157 eod1[:,1] -= np.min(eod1[:,1])
1158 eod1[:,1] /= np.max(eod1[:,1])
1159 sigma = np.std(eod1[:,1])
1160 # set time zero to maximum of second EOD:
1161 i0 = np.argmin(np.abs(eod2[:,0]))
1162 k0 = i0-n2//2
1163 if k0 < 0:
1164 k0 = 0
1165 k1 = k0 + n2 + 1
1166 if k1 >= len(eod2):
1167 k1 = len(eod2)
1168 i = k0 + np.argmax(eod2[k0:k1,1])
1169 eod2[:,0] -= eod2[i,0]
1170 # interpolate eod2 to the time base of eod1:
1171 eod2w = np.interp(eod1[:,0], eod2[:,0], eod2[:,1])
1172 # normalize amplitudes of second EOD:
1173 eod2w -= np.min(eod2w)
1174 eod2w /= np.max(eod2w)
1175 # root-mean-square difference:
1176 rmse1 = np.sqrt(np.mean((eod1[:,1] - eod2w)**2))/sigma
1177 # root-mean-square difference of the flipped signal:
1178 i = k0 + np.argmin(eod2[k0:k1,1])
1179 eod2[:,0] -= eod2[i,0]
1180 eod2w = np.interp(eod1[:,0], eod2[:,0], -eod2[:,1])
1181 eod2w -= np.min(eod2w)
1182 eod2w /= np.max(eod2w)
1183 rmse2 = np.sqrt(np.mean((eod1[:,1] - eod2w)**2))/sigma
1184 # take the smaller value:
1185 rmse = min(rmse1, rmse2)
1186 return rmse
1189def pulse_similarity(eod1, eod2, wfac=10.0):
1190 """Root-mean squared difference between two pulse fish EODs.
1192 Compute the root-mean squared difference between two pulse fish
1193 EODs over `wfac` times the distance between minimum and maximum of
1194 the wider EOD. The waveforms are normalized to their maxima prior
1195 to computing the rms difference. Also compute the rms difference
1196 between the one EOD and the other one inverted in amplitude. The
1197 smaller of the two rms values is returned.
1199 Parameters
1200 ----------
1201 eod1: 2-D array
1202 Time and amplitude of reference EOD.
1203 eod2: 2-D array
1204 Time and amplitude of EOD that is to be compared to `eod1`.
1205 wfac: float
1206 Multiply the distance between minimum and maximum by this factor
1207 to get the window size over which to compute the rms difference.
1209 Returns
1210 -------
1211 rmse: float
1212 Root-mean-squared difference between the two EOD waveforms
1213 relative to their standard deviation over the analysis window.
1214 """
1215 # copy:
1216 eod1 = np.array(eod1[:,:2])
1217 eod2 = np.array(eod2[:,:2])
1218 # width of the pulses:
1219 imin1 = np.argmin(eod1[:,1])
1220 imax1 = np.argmax(eod1[:,1])
1221 w1 = np.abs(eod1[imax1,0]-eod1[imin1,0])
1222 imin2 = np.argmin(eod2[:,1])
1223 imax2 = np.argmax(eod2[:,1])
1224 w2 = np.abs(eod2[imax2,0]-eod2[imin2,0])
1225 w = wfac*max(w1, w2)
1226 # cut out signal from the reference EOD:
1227 n = int(w/(eod1[1,0]-eod1[0,0]))
1228 i0 = imax1-n//2
1229 if i0 < 0:
1230 i0 = 0
1231 i1 = imax1+n//2+1
1232 if i1 >= len(eod1):
1233 i1 = len(eod1)
1234 eod1[:,0] -= eod1[imax1,0]
1235 eod1 = eod1[i0:i1,:]
1236 # normalize amplitude of first EOD:
1237 eod1[:,1] /= np.max(eod1[:,1])
1238 sigma = np.std(eod1[:,1])
1239 # interpolate eod2 to the time base of eod1:
1240 eod2[:,0] -= eod2[imax2,0]
1241 eod2w = np.interp(eod1[:,0], eod2[:,0], eod2[:,1])
1242 # normalize amplitude of second EOD:
1243 eod2w /= np.max(eod2w)
1244 # root-mean-square difference:
1245 rmse1 = np.sqrt(np.mean((eod1[:,1] - eod2w)**2))/sigma
1246 # root-mean-square difference of the flipped signal:
1247 eod2[:,0] -= eod2[imin2,0]
1248 eod2w = np.interp(eod1[:,0], eod2[:,0], -eod2[:,1])
1249 eod2w /= np.max(eod2w)
1250 rmse2 = np.sqrt(np.mean((eod1[:,1] - eod2w)**2))/sigma
1251 # take the smaller value:
1252 rmse = min(rmse1, rmse2)
1253 return rmse
1256def clipped_fraction(data, samplerate, eod_times, mean_eod,
1257 min_clip=-np.inf, max_clip=np.inf):
1258 """Compute fraction of clipped EOD waveform snippets.
1260 Cut out snippets at each `eod_times` based on time axis of
1261 `mean_eod`. Check which fraction of snippets exceeds clipping
1262 amplitude `min_clip` and `max_clip`.
1264 Parameters
1265 ----------
1266 data: 1-D array of float
1267 The data to be analysed.
1268 samplerate: float
1269 Sampling rate of the data in Hertz.
1270 eod_times: 1-D array of float
1271 Array of EOD times in seconds.
1272 mean_eod: 2-D array with time, mean, sem, and fit.
1273 Averaged EOD waveform of wave fish. Only the time axis is used
1274 to set width of snippets.
1275 min_clip: float
1276 Minimum amplitude that is not clipped.
1277 max_clip: float
1278 Maximum amplitude that is not clipped.
1280 Returns
1281 -------
1282 clipped_frac: float
1283 Fraction of snippets that are clipped.
1284 """
1285 # snippets:
1286 idx0 = np.argmin(np.abs(mean_eod[:,0])) # index of time zero
1287 w0 = -idx0
1288 w1 = len(mean_eod[:,0]) - idx0
1289 eod_idx = np.round(eod_times * samplerate).astype(int)
1290 eod_snippets = snippets(data, eod_idx, w0, w1)
1291 # fraction of clipped snippets:
1292 clipped_frac = np.sum(np.any((eod_snippets > max_clip) |
1293 (eod_snippets < min_clip), axis=1))\
1294 / len(eod_snippets)
1295 return clipped_frac
1298def wave_quality(props, harm_relampl=None, min_freq=0.0,
1299 max_freq=2000.0, max_clipped_frac=0.1,
1300 max_crossings=4, max_rms_sem=0.0, max_rms_error=0.05,
1301 min_power=-100.0, max_thd=0.0, max_db_diff=20.0,
1302 max_harmonics_db=-5.0, max_relampl_harm1=0.0,
1303 max_relampl_harm2=0.0, max_relampl_harm3=0.0):
1304 """Assess the quality of an EOD waveform of a wave fish.
1306 Parameters
1307 ----------
1308 props: dict
1309 A dictionary with properties of the analyzed EOD waveform
1310 as returned by `analyze_wave()`.
1311 harm_relampl: 1-D array of floats or None
1312 Relative amplitude of at least the first 3 harmonics without
1313 the fundamental.
1314 min_freq: float
1315 Minimum EOD frequency (`props['EODf']`).
1316 max_freq: float
1317 Maximum EOD frequency (`props['EODf']`).
1318 max_clipped_frac: float
1319 If larger than zero, maximum allowed fraction of clipped data
1320 (`props['clipped']`).
1321 max_crossings: int
1322 If larger than zero, maximum number of zero crossings per EOD period
1323 (`props['ncrossings']`).
1324 max_rms_sem: float
1325 If larger than zero, maximum allowed standard error of the
1326 data relative to p-p amplitude (`props['noise']`).
1327 max_rms_error: float
1328 If larger than zero, maximum allowed root-mean-square error
1329 between EOD waveform and Fourier fit relative to p-p amplitude
1330 (`props['rmserror']`).
1331 min_power: float
1332 Minimum power of the EOD in dB (`props['power']`).
1333 max_thd: float
1334 If larger than zero, then maximum total harmonic distortion
1335 (`props['thd']`).
1336 max_db_diff: float
1337 If larger than zero, maximum standard deviation of differences between
1338 logarithmic powers of harmonics in decibel (`props['dbdiff']`).
1339 Low values enforce smoother power spectra.
1340 max_harmonics_db:
1341 Maximum power of higher harmonics relative to peak power in
1342 decibel (`props['maxdb']`).
1343 max_relampl_harm1: float
1344 If larger than zero, maximum allowed amplitude of first harmonic
1345 relative to fundamental.
1346 max_relampl_harm2: float
1347 If larger than zero, maximum allowed amplitude of second harmonic
1348 relative to fundamental.
1349 max_relampl_harm3: float
1350 If larger than zero, maximum allowed amplitude of third harmonic
1351 relative to fundamental.
1353 Returns
1354 -------
1355 remove: bool
1356 If True then this is most likely not an electric fish. Remove
1357 it from both the waveform properties and the list of EOD
1358 frequencies. If False, keep it in the list of EOD
1359 frequencies, but remove it from the waveform properties if
1360 `skip_reason` is not empty.
1361 skip_reason: string
1362 An empty string if the waveform is good, otherwise a string
1363 indicating the failure.
1364 msg: string
1365 A textual representation of the values tested.
1366 """
1367 remove = False
1368 msg = []
1369 skip_reason = []
1370 # EOD frequency:
1371 if 'EODf' in props:
1372 eodf = props['EODf']
1373 msg += ['EODf=%5.1fHz' % eodf]
1374 if eodf < min_freq or eodf > max_freq:
1375 remove = True
1376 skip_reason += ['invalid EODf=%5.1fHz (minimumFrequency=%5.1fHz, maximumFrequency=%5.1f))' %
1377 (eodf, min_freq, max_freq)]
1378 # clipped fraction:
1379 if 'clipped' in props:
1380 clipped_frac = props['clipped']
1381 msg += ['clipped=%3.0f%%' % (100.0*clipped_frac)]
1382 if max_clipped_frac > 0 and clipped_frac >= max_clipped_frac:
1383 skip_reason += ['clipped=%3.0f%% (maximumClippedFraction=%3.0f%%)' %
1384 (100.0*clipped_frac, 100.0*max_clipped_frac)]
1385 # too many zero crossings:
1386 if 'ncrossings' in props:
1387 ncrossings = props['ncrossings']
1388 msg += ['zero crossings=%d' % ncrossings]
1389 if max_crossings > 0 and ncrossings > max_crossings:
1390 skip_reason += ['too many zero crossings=%d (maximumCrossings=%d)' %
1391 (ncrossings, max_crossings)]
1392 # noise:
1393 rms_sem = None
1394 if 'rmssem' in props:
1395 rms_sem = props['rmssem']
1396 if 'noise' in props:
1397 rms_sem = props['noise']
1398 if rms_sem is not None:
1399 msg += ['rms sem waveform=%6.2f%%' % (100.0*rms_sem)]
1400 if max_rms_sem > 0.0 and rms_sem >= max_rms_sem:
1401 skip_reason += ['noisy waveform s.e.m.=%6.2f%% (max %6.2f%%)' %
1402 (100.0*rms_sem, 100.0*max_rms_sem)]
1403 # fit error:
1404 if 'rmserror' in props:
1405 rms_error = props['rmserror']
1406 msg += ['rmserror=%6.2f%%' % (100.0*rms_error)]
1407 if max_rms_error > 0.0 and rms_error >= max_rms_error:
1408 skip_reason += ['noisy rmserror=%6.2f%% (maximumVariance=%6.2f%%)' %
1409 (100.0*rms_error, 100.0*max_rms_error)]
1410 # wave power:
1411 if 'power' in props:
1412 power = props['power']
1413 msg += ['power=%6.1fdB' % power]
1414 if power < min_power:
1415 skip_reason += ['small power=%6.1fdB (minimumPower=%6.1fdB)' %
1416 (power, min_power)]
1417 # total harmonic distortion:
1418 if 'thd' in props:
1419 thd = props['thd']
1420 msg += ['thd=%5.1f%%' % (100.0*thd)]
1421 if max_thd > 0.0 and thd > max_thd:
1422 skip_reason += ['large THD=%5.1f%% (maxximumTotalHarmonicDistortion=%5.1f%%)' %
1423 (100.0*thd, 100.0*max_thd)]
1424 # smoothness of spectrum:
1425 if 'dbdiff' in props:
1426 db_diff = props['dbdiff']
1427 msg += ['dBdiff=%5.1fdB' % db_diff]
1428 if max_db_diff > 0.0 and db_diff > max_db_diff:
1429 remove = True
1430 skip_reason += ['not smooth s.d. diff=%5.1fdB (maxximumPowerDifference=%5.1fdB)' %
1431 (db_diff, max_db_diff)]
1432 # maximum power of higher harmonics:
1433 if 'maxdb' in props:
1434 max_harmonics = props['maxdb']
1435 msg += ['max harmonics=%5.1fdB' % max_harmonics]
1436 if max_harmonics > max_harmonics_db:
1437 remove = True
1438 skip_reason += ['maximum harmonics=%5.1fdB too strong (maximumHarmonicsPower=%5.1fdB)' %
1439 (max_harmonics, max_harmonics_db)]
1440 # relative amplitude of harmonics:
1441 if harm_relampl is not None:
1442 for k, max_relampl in enumerate([max_relampl_harm1, max_relampl_harm2, max_relampl_harm3]):
1443 if k >= len(harm_relampl):
1444 break
1445 msg += ['ampl%d=%5.1f%%' % (k+1, 100.0*harm_relampl[k])]
1446 if max_relampl > 0.0 and k < len(harm_relampl) and harm_relampl[k] >= max_relampl:
1447 num_str = ['First', 'Second', 'Third']
1448 skip_reason += ['distorted ampl%d=%5.1f%% (maximum%sHarmonicAmplitude=%5.1f%%)' %
1449 (k+1, 100.0*harm_relampl[k], num_str[k], 100.0*max_relampl)]
1450 return remove, ', '.join(skip_reason), ', '.join(msg)
1453def pulse_quality(props, max_clipped_frac=0.1, max_rms_sem=0.0):
1454 """Assess the quality of an EOD waveform of a pulse fish.
1456 Parameters
1457 ----------
1458 props: dict
1459 A dictionary with properties of the analyzed EOD waveform
1460 as returned by `analyze_pulse()`.
1461 max_clipped_frac: float
1462 Maximum allowed fraction of clipped data.
1463 max_rms_sem: float
1464 If not zero, maximum allowed standard error of the data
1465 relative to p-p amplitude.
1467 Returns
1468 -------
1469 skip_reason: string
1470 An empty string if the waveform is good, otherwise a string
1471 indicating the failure.
1472 msg: string
1473 A textual representation of the values tested.
1474 skipped_clipped: bool
1475 True if waveform was skipped because of clipping.
1476 """
1477 msg = []
1478 skip_reason = []
1479 skipped_clipped = False
1480 # clipped fraction:
1481 if 'clipped' in props:
1482 clipped_frac = props['clipped']
1483 msg += ['clipped=%3.0f%%' % (100.0*clipped_frac)]
1484 if clipped_frac >= max_clipped_frac:
1485 skip_reason += ['clipped=%3.0f%% (maximumClippedFraction=%3.0f%%)' %
1486 (100.0*clipped_frac, 100.0*max_clipped_frac)]
1487 skipped_clipped = True
1488 # noise:
1489 rms_sem = None
1490 if 'rmssem' in props:
1491 rms_sem = props['rmssem']
1492 if 'noise' in props:
1493 rms_sem = props['noise']
1494 if rms_sem is not None:
1495 msg += ['rms sem waveform=%6.2f%%' % (100.0*rms_sem)]
1496 if max_rms_sem > 0.0 and rms_sem >= max_rms_sem:
1497 skip_reason += ['noisy waveform s.e.m.=%6.2f%% (maximumRMSNoise=%6.2f%%)' %
1498 (100.0*rms_sem, 100.0*max_rms_sem)]
1499 return ', '.join(skip_reason), ', '.join(msg), skipped_clipped
1502def plot_eod_recording(ax, data, samplerate, unit=None, width=0.1,
1503 toffs=0.0, kwargs={'lw': 2, 'color': 'red'}):
1504 """Plot a zoomed in range of the recorded trace.
1506 Parameters
1507 ----------
1508 ax: matplotlib axes
1509 Axes used for plotting.
1510 data: 1D ndarray
1511 Recorded data to be plotted.
1512 samplerate: float
1513 Sampling rate of the data in Hertz.
1514 unit: string
1515 Optional unit of the data used for y-label.
1516 width: float
1517 Width of data segment to be plotted in seconds.
1518 toffs: float
1519 Time of first data value in seconds.
1520 kwargs: dict
1521 Arguments passed on to the plot command for the recorded trace.
1522 """
1523 widx2 = int(width*samplerate)//2
1524 i0 = len(data)//2 - widx2
1525 i0 = (i0//widx2)*widx2
1526 i1 = i0 + 2*widx2
1527 if i0 < 0:
1528 i0 = 0
1529 if i1 >= len(data):
1530 i1 = len(data)
1531 time = np.arange(len(data))/samplerate + toffs
1532 tunit = 'sec'
1533 if np.abs(time[i0]) < 1.0 and np.abs(time[i1]) < 1.0:
1534 time *= 1000.0
1535 tunit = 'ms'
1536 ax.plot(time, data, **kwargs)
1537 ax.set_xlim(time[i0], time[i1])
1539 ax.set_xlabel('Time [%s]' % tunit)
1540 ymin = np.min(data[i0:i1])
1541 ymax = np.max(data[i0:i1])
1542 dy = ymax - ymin
1543 ax.set_ylim(ymin-0.05*dy, ymax+0.05*dy)
1544 if len(unit) == 0 or unit == 'a.u.':
1545 ax.set_ylabel('Amplitude')
1546 else:
1547 ax.set_ylabel('Amplitude [%s]' % unit)
1550def plot_pulse_eods(ax, data, samplerate, zoom_window, width, eod_props,
1551 toffs=0.0, colors=None, markers=None, marker_size=10,
1552 legend_rows=8, **kwargs):
1553 """Mark pulse EODs in a plot of an EOD recording.
1555 Parameters
1556 ----------
1557 ax: matplotlib axes
1558 Axes used for plotting.
1559 data: 1D ndarray
1560 Recorded data (these are not plotted!).
1561 samplerate: float
1562 Sampling rate of the data in Hertz.
1563 zoom_window: tuple of floats
1564 Start and end time of the data to be plotted in seconds.
1565 width: float
1566 Minimum width of the data to be plotted in seconds.
1567 eod_props: list of dictionaries
1568 Lists of EOD properties as returned by `analyze_pulse()`
1569 and `analyze_wave()`. From the entries with 'type' ==
1570 'pulse' the properties 'EODf' and 'times' are used. 'EODf'
1571 is the averaged EOD frequency, and 'times' is a list of
1572 detected EOD pulse times.
1573 toffs: float
1574 Time of first data value in seconds that will be added
1575 to the pulse times in `eod_props`.
1576 colors: list of colors or None
1577 If not None list of colors for plotting each cluster
1578 markers: list of markers or None
1579 If not None list of markers for plotting each cluster
1580 marker_size: float
1581 Size of markers used to mark the pulses.
1582 legend_rows: int
1583 Maximum number of rows to be used for the legend.
1584 kwargs:
1585 Key word arguments for the legend of the plot.
1586 """
1587 k = 0
1588 for eod in eod_props:
1589 if eod['type'] != 'pulse':
1590 continue
1591 if 'times' not in eod:
1592 continue
1594 width = np.min([width, np.diff(zoom_window)[0]])
1595 while len(eod['peaktimes'][(eod['peaktimes']>(zoom_window[1]-width)) & (eod['peaktimes']<(zoom_window[1]))]) == 0:
1596 width *= 2
1597 if zoom_window[1] - width < 0:
1598 width = width/2
1599 break
1601 x = eod['peaktimes'] + toffs
1602 pidx = np.round(eod['peaktimes']*samplerate).astype(int)
1603 y = data[pidx[(pidx>0)&(pidx<len(data))]]
1604 x = x[(pidx>0)&(pidx<len(data))]
1605 color_kwargs = {}
1606 #if colors is not None:
1607 # color_kwargs['color'] = colors[k%len(colors)]
1608 if marker_size is not None:
1609 color_kwargs['ms'] = marker_size
1610 label = '%6.1f Hz' % eod['EODf']
1611 if legend_rows > 5 and k >= legend_rows:
1612 label = None
1613 if markers is None:
1614 ax.plot(x, y, 'o', label=label, zorder=-1, **color_kwargs)
1615 else:
1616 ax.plot(x, y, linestyle='none', label=label,
1617 marker=markers[k%len(markers)], mec=None, mew=0.0,
1618 zorder=-1, **color_kwargs)
1619 k += 1
1621 # legend:
1622 if k > 1:
1623 if legend_rows > 0:
1624 if legend_rows > 5:
1625 ncol = 1
1626 else:
1627 ncol = (len(idx)-1) // legend_rows + 1
1628 ax.legend(numpoints=1, ncol=ncol, **kwargs)
1629 else:
1630 ax.legend(numpoints=1, **kwargs)
1632 # reset window so at least one EOD of each cluster is visible
1633 if len(zoom_window)>0:
1634 ax.set_xlim(max(toffs,toffs+zoom_window[1]-width), toffs+zoom_window[1])
1636 i0 = max(0,int((zoom_window[1]-width)*samplerate))
1637 i1 = int(zoom_window[1]*samplerate)
1639 ymin = np.min(data[i0:i1])
1640 ymax = np.max(data[i0:i1])
1641 dy = ymax - ymin
1642 ax.set_ylim(ymin-0.05*dy, ymax+0.05*dy)
1645def plot_eod_snippets(ax, data, samplerate, tmin, tmax, eod_times,
1646 n_snippets=10, flip=False,
1647 kwargs={'zorder': -5, 'scaley': False,
1648 'lw': 0.5, 'color': '#CCCCCC'}):
1649 """Plot a few EOD waveform snippets.
1651 Parameters
1652 ----------
1653 ax: matplotlib axes
1654 Axes used for plotting.
1655 data: 1D ndarray
1656 Recorded data from which the snippets are taken.
1657 samplerate: float
1658 Sampling rate of the data in Hertz.
1659 tmin: float
1660 Start time of each snippet.
1661 tmax: float
1662 End time of each snippet.
1663 eod_times: 1-D array
1664 EOD peak times from which a few are selected to be plotted.
1665 n_snippets: int
1666 Number of snippets to be plotted. If zero do not plot anything.
1667 flip: bool
1668 If True flip the snippets upside down.
1669 kwargs: dict
1670 Arguments passed on to the plot command for plotting the snippets.
1671 """
1672 if n_snippets <= 0:
1673 return
1674 i0 = int(tmin*samplerate)
1675 i1 = int(tmax*samplerate)
1676 time = 1000.0*np.arange(i0, i1)/samplerate
1677 step = len(eod_times)//n_snippets
1678 if step < 1:
1679 step = 1
1680 for t in eod_times[n_snippets//2::step]:
1681 idx = int(np.round(t*samplerate))
1682 if idx+i0 < 0 or idx+i1 >= len(data):
1683 continue
1684 snippet = data[idx+i0:idx+i1]
1685 if flip:
1686 snippet *= -1
1687 ax.plot(time, snippet - np.mean(snippet[:len(snippet)//4]), **kwargs)
1690def plot_eod_waveform(ax, eod_waveform, props, peaks=None, unit=None,
1691 mkwargs={'zorder': 10, 'lw': 2, 'color': 'red'},
1692 skwargs={'zorder': -10, 'color': '#CCCCCC'},
1693 fkwargs={'zorder': 5, 'lw': 6, 'color': 'steelblue'},
1694 zkwargs={'zorder': -5, 'lw': 1, 'color': '#AAAAAA'}):
1695 """Plot mean EOD, its standard error, and an optional fit to the EOD.
1697 Parameters
1698 ----------
1699 ax: matplotlib axes
1700 Axes used for plotting.
1701 eod_waveform: 2-D array
1702 EOD waveform. First column is time in seconds, second column
1703 the (mean) eod waveform. The optional third column is the
1704 standard error, and the optional fourth column is a fit on the
1705 waveform.
1706 props: dict
1707 A dictionary with properties of the analyzed EOD waveform as
1708 returned by `analyze_wave()` and `analyze_pulse()`.
1709 peaks: 2_D arrays or None
1710 List of peak properties (index, time, and amplitude) of a EOD pulse
1711 as returned by `analyze_pulse()`.
1712 unit: string
1713 Optional unit of the data used for y-label.
1714 mkwargs: dict
1715 Arguments passed on to the plot command for the mean EOD.
1716 skwargs: dict
1717 Arguments passed on to the fill_between command for the
1718 standard error of the EOD.
1719 fkwargs: dict
1720 Arguments passed on to the plot command for the fitted EOD.
1721 zkwargs: dict
1722 Arguments passed on to the plot command for the zero line.
1723 """
1724 ax.autoscale(True)
1725 time = 1000.0 * eod_waveform[:,0]
1726 # plot zero line:
1727 ax.axhline(0.0, **zkwargs)
1728 # plot fit:
1729 if eod_waveform.shape[1] > 3:
1730 ax.plot(time, eod_waveform[:,3], **fkwargs)
1731 # plot waveform:
1732 mean_eod = eod_waveform[:,1]
1733 ax.plot(time, mean_eod, **mkwargs)
1734 # plot standard error:
1735 if eod_waveform.shape[1] > 2:
1736 std_eod = eod_waveform[:,2]
1737 if np.mean(std_eod)/(np.max(mean_eod) - np.min(mean_eod)) > 0.1:
1738 ax.autoscale_view(False)
1739 ax.autoscale(False)
1740 ax.fill_between(time, mean_eod + std_eod, mean_eod - std_eod,
1741 **skwargs)
1742 # ax height dimensions:
1743 pixely = np.abs(np.diff(ax.get_window_extent().get_points()[:,1]))[0]
1744 ymin, ymax = ax.get_ylim()
1745 unity = ymax - ymin
1746 dyu = np.abs(unity)/pixely
1747 font_size = plt.rcParams['font.size']*dyu
1748 # annotate fit:
1749 tau = None if props is None else props.get('tau', None)
1750 ty = 0.0
1751 if tau is not None and eod_waveform.shape[1] > 3:
1752 if tau < 0.001:
1753 label = u'\u03c4=%.0f\u00b5s' % (1.e6*tau)
1754 else:
1755 label = u'\u03c4=%.2fms' % (1.e3*tau)
1756 inx = np.argmin(np.isnan(eod_waveform[:,3]))
1757 x = eod_waveform[inx,0] + 1.5*tau
1758 ty = 0.7*eod_waveform[inx,3]
1759 if np.abs(ty) < 0.5*font_size:
1760 ty = 0.5*font_size*np.sign(ty)
1761 va = 'bottom' if ty > 0.0 else 'top'
1762 ax.text(1000.0*x, ty, label, ha='left', va=va, zorder=20)
1763 # annotate peaks:
1764 if peaks is not None and len(peaks)>0:
1765 for p in peaks:
1766 ax.scatter(1000.0*p[1], p[2], s=80, clip_on=False, zorder=0, alpha=0.4,
1767 c=mkwargs['color'], edgecolors=mkwargs['color'])
1768 label = u'P%d' % p[0]
1769 if p[0] != 1:
1770 if np.abs(p[1]) < 0.001:
1771 ts = u'%.0f\u00b5s' % (1.0e6*p[1])
1772 elif np.abs(p[1]) < 0.01:
1773 ts = u'%.2fms' % (1.0e3*p[1])
1774 else:
1775 ts = u'%.3gms' % (1.0e3*p[1])
1776 if np.abs(p[3]) < 0.05:
1777 ps = u'%.1f%%' % (100*p[3])
1778 else:
1779 ps = u'%.0f%%' % (100*p[3])
1780 label += u'(%s @ %s)' % (ps, ts)
1781 va = 'bottom'
1782 dy = 0.4*font_size
1783 if p[0] % 2 == 0:
1784 va = 'top'
1785 dy = -dy
1786 if p[0] == 1:
1787 dy = 0.0
1788 """
1789 if p[2] <= np.min(peaks[:,2]):
1790 dy = -0.8*font_size
1791 va = 'bottom'
1792 """
1793 if p[2] + dy < ymin + 1.3*font_size:
1794 dy = ymin + 1.3*font_size - p[2]
1795 sign = np.sign(p[2])
1796 if p[0] == np.max(peaks[:,0]) and ty*p[2] > 0.0 and \
1797 sign*p[2]+dy < sign*ty+1.2*font_size:
1798 dy = ty + sign*1.2*font_size - p[2]
1799 dx = 0.05*time[-1]
1800 if p[1] >= 0.0:
1801 ax.text(1000.0*p[1]+dx, p[2]+dy, label,
1802 ha='left', va=va, zorder=20)
1803 else:
1804 ax.text(1000.0*p[1]-dx, p[2]+dy, label,
1805 ha='right', va=va, zorder=20)
1806 # annotate plot:
1807 if unit is None or len(unit) == 0 or unit == 'a.u.':
1808 unit = ''
1809 if props is not None:
1810 props['unit'] = unit
1811 label = 'p-p amplitude = {p-p-amplitude:.3g} {unit}\n'.format(**props)
1812 if 'n' in props:
1813 props['eods'] = 'EODs' if props['n'] > 1 else 'EOD'
1814 label += 'n = {n} {eods}\n'.format(**props)
1815 if props['flipped']:
1816 label += 'flipped\n'
1817 if -eod_waveform[0,0] < 0.6*eod_waveform[-1,0]:
1818 ax.text(0.97, 0.97, label, transform=ax.transAxes,
1819 va='top', ha='right', zorder=20)
1820 else:
1821 ax.text(0.03, 0.97, label, transform=ax.transAxes,
1822 va='top', zorder=20)
1823 # axis:
1824 if props is not None and props['type'] == 'wave':
1825 lim = 750.0/props['EODf']
1826 ax.set_xlim([-lim, +lim])
1827 else:
1828 ax.set_xlim(time[0], time[-1])
1829 ax.set_xlabel('Time [msec]')
1830 if unit:
1831 ax.set_ylabel('Amplitude [%s]' % unit)
1832 else:
1833 ax.set_ylabel('Amplitude')
1836def plot_wave_spectrum(axa, axp, spec, props, unit=None,
1837 color='b', lw=2, markersize=10):
1838 """Plot and annotate spectrum of wave EOD.
1840 Parameters
1841 ----------
1842 axa: matplotlib axes
1843 Axes for amplitude plot.
1844 axp: matplotlib axes
1845 Axes for phase plot.
1846 spec: 2-D array
1847 The amplitude spectrum of a single pulse as returned by
1848 `analyze_wave()`. First column is the index of the harmonics,
1849 second column its frequency, third column its amplitude,
1850 fourth column its amplitude relative to the fundamental, fifth
1851 column is power of harmonics relative to fundamental in
1852 decibel, and sixth column the phase shift relative to the
1853 fundamental.
1854 props: dict
1855 A dictionary with properties of the analyzed EOD waveform as
1856 returned by `analyze_wave()`.
1857 unit: string
1858 Optional unit of the data used for y-label.
1859 color:
1860 Color for line and points of spectrum.
1861 lw: float
1862 Linewidth for spectrum.
1863 markersize: float
1864 Size of points on spectrum.
1865 """
1866 n = min(9, np.sum(np.isfinite(spec[:,2])))
1867 # amplitudes:
1868 markers, stemlines, _ = axa.stem(spec[:n,0]+1, spec[:n,2], basefmt='none')
1869 plt.setp(markers, color=color, markersize=markersize, clip_on=False)
1870 plt.setp(stemlines, color=color, lw=lw)
1871 axa.set_xlim(0.5, n+0.5)
1872 axa.set_ylim(bottom=0)
1873 axa.xaxis.set_major_locator(plt.MultipleLocator(1))
1874 axa.tick_params('x', direction='out')
1875 if unit:
1876 axa.set_ylabel('Amplitude [%s]' % unit)
1877 else:
1878 axa.set_ylabel('Amplitude')
1879 # phases:
1880 phases = spec[:n,5]
1881 phases[phases<0.0] = phases[phases<0.0] + 2.0*np.pi
1882 markers, stemlines, _ = axp.stem(spec[:n,0]+1, phases[:n], basefmt='none')
1883 plt.setp(markers, color=color, markersize=markersize, clip_on=False)
1884 plt.setp(stemlines, color=color, lw=lw)
1885 axp.set_xlim(0.5, n+0.5)
1886 axp.xaxis.set_major_locator(plt.MultipleLocator(1))
1887 axp.tick_params('x', direction='out')
1888 axp.set_ylim(0, 2.0*np.pi)
1889 axp.set_yticks([0, np.pi, 2.0*np.pi])
1890 axp.set_yticklabels([u'0', u'\u03c0', u'2\u03c0'])
1891 axp.set_xlabel('Harmonics')
1892 axp.set_ylabel('Phase')
1895def plot_pulse_spectrum(ax, power, props, min_freq=1.0, max_freq=10000.0,
1896 color='b', lw=3, markersize=80):
1897 """Plot and annotate spectrum of single pulse EOD.
1899 Parameters
1900 ----------
1901 ax: matplotlib axes
1902 Axes used for plotting.
1903 power: 2-D array
1904 The power spectrum of a single pulse as returned by `analyze_pulse()`.
1905 First column are the frequencies, second column the power.
1906 props: dict
1907 A dictionary with properties of the analyzed EOD waveform as
1908 returned by `analyze_pulse()`.
1909 min_freq: float
1910 Minimun frequency of the spectrum to be plotted (logscale!).
1911 max_freq: float
1912 Maximun frequency of the spectrum to be plotted (logscale!).
1913 color:
1914 Color for line and points of spectrum.
1915 lw: float
1916 Linewidth for spectrum.
1917 markersize: float
1918 Size of points on spectrum.
1919 """
1920 box = mpatches.Rectangle((1,-60), 49, 60, linewidth=0, facecolor='#DDDDDD',
1921 zorder=1)
1922 ax.add_patch(box)
1923 att = props['poweratt50']
1924 if att < -5.0:
1925 ax.text(10.0, att+1.0, '%.0f dB' % att, ha='left', va='bottom', zorder=10)
1926 else:
1927 ax.text(10.0, att-1.0, '%.0f dB' % att, ha='left', va='top', zorder=10)
1928 box = mpatches.Rectangle((1,-60), 4, 60, linewidth=0, facecolor='#CCCCCC',
1929 zorder=2)
1930 ax.add_patch(box)
1931 att = props['poweratt5']
1932 if att < -5.0:
1933 ax.text(4.0, att+1.0, '%.0f dB' % att, ha='right', va='bottom', zorder=10)
1934 else:
1935 ax.text(4.0, att-1.0, '%.0f dB' % att, ha='right', va='top', zorder=10)
1936 lowcutoff = props['lowcutoff']
1937 if lowcutoff >= min_freq:
1938 ax.plot([lowcutoff, lowcutoff, 1.0], [-60.0, 0.5*att, 0.5*att], '#BBBBBB',
1939 zorder=3)
1940 ax.text(1.2*lowcutoff, 0.5*att-1.0, '%.0f Hz' % lowcutoff, ha='left', va='top', zorder=10)
1941 db = decibel(power[:,1])
1942 smax = np.nanmax(db)
1943 ax.plot(power[:,0], db - smax, color, lw=lw, zorder=4)
1944 peakfreq = props['peakfreq']
1945 if peakfreq >= min_freq:
1946 ax.scatter([peakfreq], [0.0], c=color, edgecolors=color, s=markersize, alpha=0.4, zorder=5)
1947 ax.text(peakfreq*1.2, 1.0, '%.0f Hz' % peakfreq, va='bottom', zorder=10)
1948 ax.set_xlim(min_freq, max_freq)
1949 ax.set_xscale('log')
1950 ax.set_ylim(-60.0, 2.0)
1951 ax.set_xlabel('Frequency [Hz]')
1952 ax.set_ylabel('Power [dB]')
1955def save_eod_waveform(mean_eod, unit, idx, basename, **kwargs):
1956 """Save mean EOD waveform to file.
1958 Parameters
1959 ----------
1960 mean_eod: 2D array of floats
1961 Averaged EOD waveform as returned by `eod_waveform()`,
1962 `analyze_wave()`, and `analyze_pulse()`.
1963 unit: string
1964 Unit of the waveform data.
1965 idx: int or None
1966 Index of fish.
1967 basename: string or stream
1968 If string, path and basename of file.
1969 '-eodwaveform', the fish index, and a file extension are appended.
1970 If stream, write EOD waveform data into this stream.
1971 kwargs:
1972 Arguments passed on to `TableData.write()`.
1974 Returns
1975 -------
1976 filename: string
1977 Path and full name of the written file in case of `basename`
1978 being a string. Otherwise, the file name and extension that
1979 would have been appended to a basename.
1981 See Also
1982 --------
1983 load_eod_waveform()
1984 """
1985 td = TableData(mean_eod[:,:3]*[1000.0, 1.0, 1.0], ['time', 'mean', 'sem'],
1986 ['ms', unit, unit], ['%.3f', '%.5f', '%.5f'])
1987 if mean_eod.shape[1] > 3:
1988 td.append('fit', unit, '%.5f', mean_eod[:,3])
1989 fp = '-eodwaveform'
1990 if idx is not None:
1991 fp += '-%d' % idx
1992 return td.write_file_stream(basename, fp, **kwargs)
1995def load_eod_waveform(file_path):
1996 """Load EOD waveform from file.
1998 Parameters
1999 ----------
2000 file_path: string
2001 Path of the file to be loaded.
2003 Returns
2004 -------
2005 mean_eod: 2D array of floats
2006 Averaged EOD waveform: time in seconds, mean, standard deviation, fit.
2007 unit: string
2008 Unit of EOD waveform.
2010 Raises
2011 ------
2012 FileNotFoundError:
2013 If `file_path` does not exist.
2015 See Also
2016 --------
2017 save_eod_waveform()
2018 """
2019 data = TableData(file_path)
2020 mean_eod = data.array()
2021 mean_eod[:,0] *= 0.001
2022 return mean_eod, data.unit('mean')
2025def save_wave_eodfs(wave_eodfs, wave_indices, basename, **kwargs):
2026 """Save frequencies of wave EODs to file.
2028 Parameters
2029 ----------
2030 wave_eodfs: list of 2D arrays
2031 Each item is a matrix with the frequencies and powers
2032 (columns) of the fundamental and harmonics (rows) as returned
2033 by `harmonics.harmonic_groups()`.
2034 wave_indices: array
2035 Indices identifying each fish or NaN.
2036 If None no index column is inserted.
2037 basename: string or stream
2038 If string, path and basename of file.
2039 '-waveeodfs' and a file extension are appended.
2040 If stream, write EOD frequencies data into this stream.
2041 kwargs:
2042 Arguments passed on to `TableData.write()`.
2044 Returns
2045 -------
2046 filename: string
2047 Path and full name of the written file in case of `basename`
2048 being a string. Otherwise, the file name and extension that
2049 would have been appended to a basename.
2051 See Also
2052 --------
2053 load_wave_eodfs()
2055 """
2056 eodfs = fundamental_freqs_and_power(wave_eodfs)
2057 td = TableData()
2058 if wave_indices is not None:
2059 td.append('index', '', '%d', [wi if wi >= 0 else np.nan for wi in wave_indices])
2060 td.append('EODf', 'Hz', '%7.2f', eodfs[:,0])
2061 td.append('datapower', 'dB', '%7.2f', eodfs[:,1])
2062 fp = '-waveeodfs'
2063 return td.write_file_stream(basename, fp, **kwargs)
2066def load_wave_eodfs(file_path):
2067 """Load frequencies of wave EODs from file.
2069 Parameters
2070 ----------
2071 file_path: string
2072 Path of the file to be loaded.
2074 Returns
2075 -------
2076 eodfs: 2D array of floats
2077 EODfs and power of wave type fish.
2078 Indices can contain NaNs.
2079 indices: array of ints
2080 Corresponding indices of fish, can contain negative numbers to
2081 indicate frequencies without fish.
2083 Raises
2084 ------
2085 FileNotFoundError:
2086 If `file_path` does not exist.
2088 See Also
2089 --------
2090 save_wave_eodfs()
2091 """
2092 data = TableData(file_path)
2093 eodfs = data.array()
2094 if 'index' in data:
2095 indices = data[:,'index']
2096 indices[~np.isfinite(indices)] = -1
2097 indices = np.array(indices, dtype=int)
2098 eodfs = eodfs[:,1:]
2099 else:
2100 indices = np.zeros(data.rows(), dtype=int) - 1
2101 return eodfs, indices
2104def save_wave_fish(eod_props, unit, basename, **kwargs):
2105 """Save properties of wave EODs to file.
2107 Parameters
2108 ----------
2109 eod_props: list of dict
2110 Properties of EODs as returned by `analyze_wave()` and
2111 `analyze_pulse()`. Only properties of wave fish are saved.
2112 unit: string
2113 Unit of the waveform data.
2114 basename: string or stream
2115 If string, path and basename of file.
2116 '-wavefish' and a file extension are appended.
2117 If stream, write wave fish properties into this stream.
2118 kwargs:
2119 Arguments passed on to `TableData.write()`.
2121 Returns
2122 -------
2123 filename: string or None
2124 Path and full name of the written file in case of `basename`
2125 being a string. Otherwise, the file name and extension that
2126 would have been appended to a basename.
2127 None if no wave fish are contained in eod_props and
2128 consequently no file was written.
2130 See Also
2131 --------
2132 load_wave_fish()
2133 """
2134 wave_props = [p for p in eod_props if p['type'] == 'wave']
2135 if len(wave_props) == 0:
2136 return None
2137 td = TableData()
2138 if 'twin' in wave_props[0] or 'samplerate' in wave_props[0] or \
2139 'nfft' in wave_props[0]:
2140 td.append_section('recording')
2141 if 'twin' in wave_props[0]:
2142 td.append('twin', 's', '%7.2f', wave_props)
2143 td.append('window', 's', '%7.2f', wave_props)
2144 td.append('winclipped', '%', '%.2f', wave_props, 100.0)
2145 if 'samplerate' in wave_props[0]:
2146 td.append('samplerate', 'kHz', '%.3f', wave_props, 0.001)
2147 if 'nfft' in wave_props[0]:
2148 td.append('nfft', '', '%d', wave_props)
2149 td.append('dfreq', 'Hz', '%.2f', wave_props)
2150 td.append_section('waveform')
2151 td.append('index', '', '%d', wave_props)
2152 td.append('EODf', 'Hz', '%7.2f', wave_props)
2153 td.append('p-p-amplitude', unit, '%.5f', wave_props)
2154 td.append('power', 'dB', '%7.2f', wave_props)
2155 if 'datapower' in wave_props[0]:
2156 td.append('datapower', 'dB', '%7.2f', wave_props)
2157 td.append('thd', '%', '%.2f', wave_props, 100.0)
2158 td.append('dbdiff', 'dB', '%7.2f', wave_props)
2159 td.append('maxdb', 'dB', '%7.2f', wave_props)
2160 if 'noise' in wave_props[0]:
2161 td.append('noise', '%', '%.1f', wave_props, 100.0)
2162 td.append('rmserror', '%', '%.2f', wave_props, 100.0)
2163 if 'clipped' in wave_props[0]:
2164 td.append('clipped', '%', '%.1f', wave_props, 100.0)
2165 td.append('flipped', '', '%d', wave_props)
2166 td.append('n', '', '%5d', wave_props)
2167 td.append_section('timing')
2168 td.append('ncrossings', '', '%d', wave_props)
2169 td.append('peakwidth', '%', '%.2f', wave_props, 100.0)
2170 td.append('troughwidth', '%', '%.2f', wave_props, 100.0)
2171 td.append('leftpeak', '%', '%.2f', wave_props, 100.0)
2172 td.append('rightpeak', '%', '%.2f', wave_props, 100.0)
2173 td.append('lefttrough', '%', '%.2f', wave_props, 100.0)
2174 td.append('righttrough', '%', '%.2f', wave_props, 100.0)
2175 td.append('p-p-distance', '%', '%.2f', wave_props, 100.0)
2176 td.append('relpeakampl', '%', '%.2f', wave_props, 100.0)
2177 fp = '-wavefish'
2178 return td.write_file_stream(basename, fp, **kwargs)
2181def load_wave_fish(file_path):
2182 """Load properties of wave EODs from file.
2184 All times are scaled to seconds, all frequencies to Hertz and all
2185 percentages to fractions.
2187 Parameters
2188 ----------
2189 file_path: string
2190 Path of the file to be loaded.
2192 Returns
2193 -------
2194 eod_props: list of dict
2195 Properties of EODs.
2197 Raises
2198 ------
2199 FileNotFoundError:
2200 If `file_path` does not exist.
2202 See Also
2203 --------
2204 save_wave_fish()
2206 """
2207 data = TableData(file_path)
2208 eod_props = data.dicts()
2209 for props in eod_props:
2210 if 'winclipped' in props:
2211 props['winclipped'] /= 100
2212 if 'samplerate' in props:
2213 props['samplerate'] *= 1000
2214 if 'nfft' in props:
2215 props['nfft'] = int(props['nfft'])
2216 props['index'] = int(props['index'])
2217 props['n'] = int(props['n'])
2218 props['type'] = 'wave'
2219 props['thd'] /= 100
2220 props['noise'] /= 100
2221 props['rmserror'] /= 100
2222 if 'clipped' in props:
2223 props['clipped'] /= 100
2224 props['ncrossings'] = int(props['ncrossings'])
2225 props['peakwidth'] /= 100
2226 props['troughwidth'] /= 100
2227 props['leftpeak'] /= 100
2228 props['rightpeak'] /= 100
2229 props['lefttrough'] /= 100
2230 props['righttrough'] /= 100
2231 props['p-p-distance'] /= 100
2232 props['relpeakampl'] /= 100
2233 return eod_props
2236def save_pulse_fish(eod_props, unit, basename, **kwargs):
2237 """Save properties of pulse EODs to file.
2239 Parameters
2240 ----------
2241 eod_props: list of dict
2242 Properties of EODs as returned by `analyze_wave()` and
2243 `analyze_pulse()`. Only properties of pulse fish are saved.
2244 unit: string
2245 Unit of the waveform data.
2246 basename: string or stream
2247 If string, path and basename of file.
2248 '-pulsefish' and a file extension are appended.
2249 If stream, write pulse fish properties into this stream.
2250 kwargs:
2251 Arguments passed on to `TableData.write()`.
2253 Returns
2254 -------
2255 filename: string or None
2256 Path and full name of the written file in case of `basename`
2257 being a string. Otherwise, the file name and extension that
2258 would have been appended to a basename.
2259 None if no pulse fish are contained in eod_props and
2260 consequently no file was written.
2262 See Also
2263 --------
2264 load_pulse_fish()
2265 """
2266 pulse_props = [p for p in eod_props if p['type'] == 'pulse']
2267 if len(pulse_props) == 0:
2268 return None
2269 td = TableData()
2270 if 'twin' in pulse_props[0] or 'samplerate' in pulse_props[0] or \
2271 'nfft' in pulse_props[0]:
2272 td.append_section('recording')
2273 if 'twin' in pulse_props[0]:
2274 td.append('twin', 's', '%7.2f', pulse_props)
2275 td.append('window', 's', '%7.2f', pulse_props)
2276 td.append('winclipped', '%', '%.2f', pulse_props, 100.0)
2277 if 'samplerate' in pulse_props[0]:
2278 td.append('samplerate', 'kHz', '%.3f', pulse_props, 0.001)
2279 if 'nfft' in pulse_props[0]:
2280 td.append('nfft', '', '%d', pulse_props)
2281 td.append('dfreq', 'Hz', '%.2f', pulse_props)
2282 td.append_section('waveform')
2283 td.append('index', '', '%d', pulse_props)
2284 td.append('EODf', 'Hz', '%7.2f', pulse_props)
2285 td.append('period', 'ms', '%7.2f', pulse_props, 1000.0)
2286 td.append('max-ampl', unit, '%.5f', pulse_props)
2287 td.append('min-ampl', unit, '%.5f', pulse_props)
2288 td.append('p-p-amplitude', unit, '%.5f', pulse_props)
2289 if 'noise' in pulse_props[0]:
2290 td.append('noise', '%', '%.1f', pulse_props, 100.0)
2291 if 'clipped' in pulse_props[0]:
2292 td.append('clipped', '%', '%.1f', pulse_props, 100.0)
2293 td.append('flipped', '', '%d', pulse_props)
2294 td.append('tstart', 'ms', '%.3f', pulse_props, 1000.0)
2295 td.append('tend', 'ms', '%.3f', pulse_props, 1000.0)
2296 td.append('width', 'ms', '%.3f', pulse_props, 1000.0)
2297 td.append('P2-P1-dist', 'ms', '%.3f', pulse_props, 1000.0)
2298 td.append('tau', 'ms', '%.3f', pulse_props, 1000.0)
2299 td.append('firstpeak', '', '%d', pulse_props)
2300 td.append('lastpeak', '', '%d', pulse_props)
2301 td.append('n', '', '%d', pulse_props)
2302 td.append_section('power spectrum')
2303 td.append('peakfreq', 'Hz', '%.2f', pulse_props)
2304 td.append('peakpower', 'dB', '%.2f', pulse_props)
2305 td.append('poweratt5', 'dB', '%.2f', pulse_props)
2306 td.append('poweratt50', 'dB', '%.2f', pulse_props)
2307 td.append('lowcutoff', 'Hz', '%.2f', pulse_props)
2308 fp = '-pulsefish'
2309 return td.write_file_stream(basename, fp, **kwargs)
2312def load_pulse_fish(file_path):
2313 """Load properties of pulse EODs from file.
2315 All times are scaled to seconds, all frequencies to Hertz and all
2316 percentages to fractions.
2318 Parameters
2319 ----------
2320 file_path: string
2321 Path of the file to be loaded.
2323 Returns
2324 -------
2325 eod_props: list of dict
2326 Properties of EODs.
2328 Raises
2329 ------
2330 FileNotFoundError:
2331 If `file_path` does not exist.
2333 See Also
2334 --------
2335 save_pulse_fish()
2337 """
2338 data = TableData(file_path)
2339 eod_props = data.dicts()
2340 for props in eod_props:
2341 if 'winclipped' in props:
2342 props['winclipped'] /= 100
2343 if 'samplerate' in props:
2344 props['samplerate'] *= 1000
2345 if 'nfft' in props:
2346 props['nfft'] = int(props['nfft'])
2347 props['index'] = int(props['index'])
2348 props['n'] = int(props['n'])
2349 props['firstpeak'] = int(props['firstpeak'])
2350 props['lastpeak'] = int(props['lastpeak'])
2351 props['type'] = 'pulse'
2352 if 'clipped' in props:
2353 props['clipped'] /= 100
2354 props['period'] /= 1000
2355 props['noise'] /= 100
2356 props['tstart'] /= 1000
2357 props['tend'] /= 1000
2358 props['width'] /= 1000
2359 props['P2-P1-dist'] /= 1000
2360 props['tau'] /= 1000
2361 return eod_props
2364def save_wave_spectrum(spec_data, unit, idx, basename, **kwargs):
2365 """Save amplitude and phase spectrum of wave EOD to file.
2367 Parameters
2368 ----------
2369 spec_data: 2D array of floats
2370 Amplitude and phase spectrum of wave EOD as returned by
2371 `analyze_wave()`.
2372 unit: string
2373 Unit of the waveform data.
2374 idx: int or None
2375 Index of fish.
2376 basename: string or stream
2377 If string, path and basename of file.
2378 '-wavespectrum', the fish index, and a file extension are appended.
2379 If stream, write wave spectrum into this stream.
2380 kwargs:
2381 Arguments passed on to `TableData.write()`.
2383 Returns
2384 -------
2385 filename: string
2386 Path and full name of the written file in case of `basename`
2387 being a string. Otherwise, the file name and extension that
2388 would have been appended to a basename.
2390 See Also
2391 --------
2392 load_wave_spectrum()
2394 """
2395 td = TableData(spec_data[:,:6]*[1.0, 1.0, 1.0, 100.0, 1.0, 1.0],
2396 ['harmonics', 'frequency', 'amplitude', 'relampl', 'relpower', 'phase'],
2397 ['', 'Hz', unit, '%', 'dB', 'rad'],
2398 ['%.0f', '%.2f', '%.6f', '%10.2f', '%6.2f', '%8.4f'])
2399 if spec_data.shape[1] > 6:
2400 td.append('datapower', '%s^2/Hz' % unit, '%11.4e', spec_data[:,6])
2401 fp = '-wavespectrum'
2402 if idx is not None:
2403 fp += '-%d' % idx
2404 return td.write_file_stream(basename, fp, **kwargs)
2407def load_wave_spectrum(file_path):
2408 """Load amplitude and phase spectrum of wave EOD from file.
2410 Parameters
2411 ----------
2412 file_path: string
2413 Path of the file to be loaded.
2415 Returns
2416 -------
2417 spec: 2D array of floats
2418 Amplitude and phase spectrum of wave EOD:
2419 harmonics, frequency, amplitude, relative amplitude in dB,
2420 relative power in dB, phase, data power in unit squared.
2421 Can contain NaNs.
2422 unit: string
2423 Unit of amplitudes.
2425 Raises
2426 ------
2427 FileNotFoundError:
2428 If `file_path` does not exist.
2430 See Also
2431 --------
2432 save_wave_spectrum()
2433 """
2434 data = TableData(file_path)
2435 spec = data.array()
2436 spec[:,3] *= 0.01
2437 return spec, data.unit('amplitude')
2440def save_pulse_spectrum(spec_data, unit, idx, basename, **kwargs):
2441 """Save power spectrum of pulse EOD to file.
2443 Parameters
2444 ----------
2445 spec_data: 2D array of floats
2446 Power spectrum of single pulse as returned by `analyze_pulse()`.
2447 unit: string
2448 Unit of the waveform data.
2449 idx: int or None
2450 Index of fish.
2451 basename: string or stream
2452 If string, path and basename of file.
2453 '-pulsespectrum', the fish index, and a file extension are appended.
2454 If stream, write pulse spectrum into this stream.
2455 kwargs:
2456 Arguments passed on to `TableData.write()`.
2458 Returns
2459 -------
2460 filename: string
2461 Path and full name of the written file in case of `basename`
2462 being a string. Otherwise, the file name and extension that
2463 would have been appended to a basename.
2465 See Also
2466 --------
2467 load_pulse_spectrum()
2468 """
2469 td = TableData(spec_data[:,:2], ['frequency', 'power'],
2470 ['Hz', '%s^2/Hz' % unit], ['%.2f', '%.4e'])
2471 fp = '-pulsespectrum'
2472 if idx is not None:
2473 fp += '-%d' % idx
2474 return td.write_file_stream(basename, fp, **kwargs)
2477def load_pulse_spectrum(file_path):
2478 """Load power spectrum of pulse EOD from file.
2480 Parameters
2481 ----------
2482 file_path: string
2483 Path of the file to be loaded.
2485 Returns
2486 -------
2487 spec: 2D array of floats
2488 Power spectrum of single pulse: frequency, power
2490 Raises
2491 ------
2492 FileNotFoundError:
2493 If `file_path` does not exist.
2495 See Also
2496 --------
2497 save_pulse_spectrum()
2498 """
2499 data = TableData(file_path)
2500 spec = data.array()
2501 return spec
2504def save_pulse_peaks(peak_data, unit, idx, basename, **kwargs):
2505 """Save peak properties of pulse EOD to file.
2507 Parameters
2508 ----------
2509 peak_data: 2D array of floats
2510 Properties of peaks and troughs of pulse EOD as returned by
2511 `analyze_pulse()`.
2512 unit: string
2513 Unit of the waveform data.
2514 idx: int or None
2515 Index of fish.
2516 basename: string or stream
2517 If string, path and basename of file.
2518 '-pulsepeaks', the fish index, and a file extension are appended.
2519 If stream, write pulse peaks into this stream.
2520 kwargs:
2521 Arguments passed on to `TableData.write()`.
2523 Returns
2524 -------
2525 filename: string
2526 Path and full name of the written file in case of `basename`
2527 being a string. Otherwise, the file name and extension that
2528 would have been appended to a basename.
2530 See Also
2531 --------
2532 load_pulse_peaks()
2533 """
2534 if len(peak_data) == 0:
2535 return None
2536 td = TableData(peak_data[:,:5]*[1.0, 1000.0, 1.0, 100.0, 1000.0],
2537 ['P', 'time', 'amplitude', 'relampl', 'width'],
2538 ['', 'ms', unit, '%', 'ms'],
2539 ['%.0f', '%.3f', '%.5f', '%.2f', '%.3f'])
2540 fp = '-pulsepeaks'
2541 if idx is not None:
2542 fp += '-%d' % idx
2543 return td.write_file_stream(basename, fp, **kwargs)
2546def load_pulse_peaks(file_path):
2547 """Load peak properties of pulse EOD from file.
2549 Parameters
2550 ----------
2551 file_path: string
2552 Path of the file to be loaded.
2554 Returns
2555 -------
2556 peak_data: 2D array of floats
2557 Properties of peaks and troughs of pulse EOD:
2558 P, time, amplitude, relampl, width
2559 unit: string
2560 Unit of peak amplitudes.
2562 Raises
2563 ------
2564 FileNotFoundError:
2565 If `file_path` does not exist.
2567 See Also
2568 --------
2569 save_pulse_peaks()
2570 """
2571 data = TableData(file_path)
2572 peaks = data.array()
2573 peaks[:,1] *= 0.001
2574 peaks[:,3] *= 0.01
2575 peaks[:,4] *= 0.001
2576 return peaks, data.unit('amplitude')
2579def save_pulse_times(pulse_times, idx, basename, **kwargs):
2580 """Save times of pulse EOD to file.
2582 Parameters
2583 ----------
2584 pulse_times: dict or array of floats
2585 Times of EOD pulses. Either as array of times or
2586 `props['peaktimes']` or `props['times']` as returned by
2587 `analyze_pulse()`.
2588 idx: int or None
2589 Index of fish.
2590 basename: string or stream
2591 If string, path and basename of file.
2592 '-pulsetimes', the fish index, and a file extension are appended.
2593 If stream, write pulse times into this stream.
2594 kwargs:
2595 Arguments passed on to `TableData.write()`.
2597 Returns
2598 -------
2599 filename: string
2600 Path and full name of the written file in case of `basename`
2601 being a string. Otherwise, the file name and extension that
2602 would have been appended to a basename.
2604 See Also
2605 --------
2606 load_pulse_times()
2607 """
2608 if isinstance(pulse_times, dict):
2609 props = pulse_times
2610 pulse_times = props.get('times', [])
2611 pulse_times = props.get('peaktimes', pulse_times)
2612 if len(pulse_times) == 0:
2613 return None
2614 td = TableData()
2615 td.append('time', 's', '%.4f', pulse_times)
2616 fp = '-pulsetimes'
2617 if idx is not None:
2618 fp += '-%d' % idx
2619 return td.write_file_stream(basename, fp, **kwargs)
2622def load_pulse_times(file_path):
2623 """Load times of pulse EOD from file.
2625 Parameters
2626 ----------
2627 file_path: string
2628 Path of the file to be loaded.
2630 Returns
2631 -------
2632 pulse_times: array of floats
2633 Times of pulse EODs in seconds.
2635 Raises
2636 ------
2637 FileNotFoundError:
2638 If `file_path` does not exist.
2640 See Also
2641 --------
2642 save_pulse_times()
2643 """
2644 data = TableData(file_path)
2645 pulse_times = data.array()[:,0]
2646 return pulse_times
2649file_types = ['waveeodfs', 'wavefish', 'pulsefish', 'eodwaveform',
2650 'wavespectrum', 'pulsepeaks', 'pulsespectrum', 'pulsetimes']
2651"""List of all file types generated and supported by the `save_*` and `load_*` functions."""
2654def parse_filename(file_path):
2655 """Parse components of an EOD analysis file name.
2657 Analysis files generated by the `eodanalysis` module are named
2658 according to
2659 ```plain
2660 PATH/RECORDING-CHANNEL-TIME-FTYPE-N.EXT
2661 ```
2663 Parameters
2664 ----------
2665 file_path: string
2666 Path of the file to be parsed.
2668 Returns
2669 -------
2670 recording: string
2671 Path and basename of the recording, i.e. 'PATH/RECORDING'.
2672 A leading './' is removed.
2673 base_path: string
2674 Path and basename of the analysis results,
2675 i.e. 'PATH/RECORDING-CHANNEL-TIME'. A leading './' is removed.
2676 channel: int
2677 Channel of the recording
2678 ('CHANNEL' component of the file name if present).
2679 -1 if not present in `file_path`.
2680 time: float
2681 Start time of analysis window in seconds
2682 ('TIME' component of the file name if present).
2683 `None` if not present in `file_path`.
2684 ftype: string
2685 Type of analysis file (e.g. 'wavespectrum', 'pulsepeaks', etc.),
2686 ('FTYPE' component of the file name if present).
2687 See `file_types` for a list of all supported file types.
2688 Empty string if not present in `file_path`.
2689 index: int
2690 Index of the EOD.
2691 ('N' component of the file name if present).
2692 -1 if not present in `file_path`.
2693 ext: string
2694 File extension *without* leading period
2695 ('EXT' component of the file name).
2697 """
2698 name, ext = os.path.splitext(file_path)
2699 ext = ext[1:]
2700 parts = name.split('-')
2701 index = -1
2702 if len(parts) > 0 and parts[-1].isdigit():
2703 index = int(parts[-1])
2704 parts = parts[:-1]
2705 ftype = ''
2706 if len(parts) > 0:
2707 ftype = parts[-1]
2708 parts = parts[:-1]
2709 base_path = '-'.join(parts)
2710 if base_path.startswith('./'):
2711 base_path = base_path[2:]
2712 time = None
2713 if len(parts) > 0 and len(parts[-1]) > 0 and \
2714 parts[-1][0] == 't' and parts[-1][-1] == 's' and \
2715 parts[-1][1:-1].isdigit():
2716 time = float(parts[-1][1:-1])
2717 parts = parts[:-1]
2718 channel = -1
2719 if len(parts) > 0 and len(parts[-1]) > 0 and \
2720 parts[-1][0] == 'c' and parts[-1][1:].isdigit():
2721 channel = int(parts[-1][1:])
2722 parts = parts[:-1]
2723 recording = '-'.join(parts)
2724 if recording.startswith('./'):
2725 recording = recording[2:]
2726 return recording, base_path, channel, time, ftype, index, ext
2729def save_analysis(output_basename, zip_file, eod_props, mean_eods,
2730 spec_data, peak_data, wave_eodfs, wave_indices, unit,
2731 verbose, **kwargs):
2732 """Save EOD analysis results to files.
2734 Parameters
2735 ----------
2736 output_basename: string
2737 Path and basename of files to be saved.
2738 zip_file: bool
2739 If `True`, write all analysis results into a zip archive.
2740 eod_props: list of dict
2741 Properties of EODs as returned by `analyze_wave()` and
2742 `analyze_pulse()`.
2743 mean_eods: list of 2D array of floats
2744 Averaged EOD waveforms as returned by `eod_waveform()`,
2745 `analyze_wave()`, and `analyze_pulse()`.
2746 spec_data: list of 2D array of floats
2747 Power spectra of single pulses as returned by
2748 `analyze_pulse()`.
2749 peak_data: list of 2D array of floats
2750 Properties of peaks and troughs of pulse EODs as returned by
2751 `analyze_pulse()`.
2752 wave_eodfs: list of 2D array of float
2753 Each item is a matrix with the frequencies and powers
2754 (columns) of the fundamental and harmonics (rows) as returned
2755 by `harmonics.harmonic_groups()`.
2756 wave_indices: array of int
2757 Indices identifying each fish in `wave_eodfs` or NaN. unit:
2758 string Unit of the waveform data.
2759 verbose: int
2760 Verbosity level.
2761 kwargs:
2762 Arguments passed on to `TableData.write()`.
2763 """
2764 def write_file_zip(zf, save_func, output, *args, **kwargs):
2765 if zf is None:
2766 fp = save_func(*args, basename=output, **kwargs)
2767 if verbose > 0 and fp is not None:
2768 print('wrote file %s' % fp)
2769 else:
2770 with io.StringIO() as df:
2771 fp = save_func(*args, basename=df, **kwargs)
2772 if fp is not None:
2773 fp = output_basename + fp
2774 zf.writestr(os.path.basename(fp), df.getvalue())
2775 if verbose > 0:
2776 print('zipped file %s' % fp)
2779 if 'table_format' in kwargs and kwargs['table_format'] == 'py':
2780 with open(output_basename+'.py', 'w') as f:
2781 name = os.path.basename(output_basename)
2782 for k, sdata in enumerate(spec_data):
2783 # save wave fish only:
2784 if len(sdata)>0 and sdata.shape[1] > 2:
2785 fish = dict(amplitudes=sdata[:,3], phases=sdata[:,5])
2786 fish = normalize_wavefish(fish)
2787 export_wavefish(fish, name+'-%d_harmonics' % k, f)
2788 else:
2789 zf = None
2790 if zip_file:
2791 zf = zipfile.ZipFile(output_basename + '.zip', 'w')
2792 # all wave fish in wave_eodfs:
2793 if len(wave_eodfs) > 0:
2794 write_file_zip(zf, save_wave_eodfs, output_basename,
2795 wave_eodfs, wave_indices, **kwargs)
2796 # all wave and pulse fish:
2797 for i, (mean_eod, sdata, pdata, props) in enumerate(zip(mean_eods, spec_data, peak_data, eod_props)):
2798 write_file_zip(zf, save_eod_waveform, output_basename,
2799 mean_eod, unit, i, **kwargs)
2800 # power spectrum:
2801 if len(sdata)>0:
2802 if sdata.shape[1] == 2:
2803 write_file_zip(zf, save_pulse_spectrum, output_basename,
2804 sdata, unit, i, **kwargs)
2805 else:
2806 write_file_zip(zf, save_wave_spectrum, output_basename,
2807 sdata, unit, i, **kwargs)
2808 # peaks:
2809 write_file_zip(zf, save_pulse_peaks, output_basename,
2810 pdata, unit, i, **kwargs)
2811 # times:
2812 write_file_zip(zf, save_pulse_times, output_basename,
2813 props, i, **kwargs)
2814 # wave fish properties:
2815 write_file_zip(zf, save_wave_fish, output_basename,
2816 eod_props, unit, **kwargs)
2817 # pulse fish properties:
2818 write_file_zip(zf, save_pulse_fish, output_basename,
2819 eod_props, unit, **kwargs)
2822def load_analysis(file_pathes):
2823 """Load all EOD analysis files.
2825 Parameters
2826 ----------
2827 file_pathes: list of string
2828 Pathes of the analysis files of a single recording to be loaded.
2830 Returns
2831 -------
2832 mean_eods: list of 2D array of floats
2833 Averaged EOD waveforms: time in seconds, mean, standard deviation, fit.
2834 wave_eodfs: 2D array of floats
2835 EODfs and power of wave type fish.
2836 wave_indices: array of ints
2837 Corresponding indices of fish, can contain negative numbers to
2838 indicate frequencies without fish.
2839 eod_props: list of dict
2840 Properties of EODs. The 'index' property is an index into the
2841 reurned lists.
2842 spec_data: list of 2D array of floats
2843 Amplitude and phase spectrum of wave-type EODs with columns
2844 harmonics, frequency, amplitude, relative amplitude in dB,
2845 relative power in dB, phase, data power in unit squared.
2846 Power spectrum of single pulse-type EODs with columns frequency, power
2847 peak_data: list of 2D array of floats
2848 Properties of peaks and troughs of pulse-type EODs with columns
2849 P, time, amplitude, relampl, width
2850 recording: string
2851 Path and base name of the recording file.
2852 channel: int
2853 Analysed channel of the recording.
2854 unit: string
2855 Unit of EOD waveform.
2856 """
2857 recording = None
2858 channel = -1
2859 eod_props = []
2860 zf = None
2861 if len(file_pathes) == 1 and os.path.splitext(file_pathes[0])[1][1:] == 'zip':
2862 zf = zipfile.ZipFile(file_pathes[0])
2863 file_pathes = sorted(zf.namelist())
2864 # first, read wave- and pulse-fish summaries:
2865 pulse_fish = False
2866 wave_fish = False
2867 for f in file_pathes:
2868 recording, _, channel, _, ftype, _, _ = parse_filename(f)
2869 if zf is not None:
2870 f = io.TextIOWrapper(zf.open(f, 'r'))
2871 if ftype == 'wavefish':
2872 eod_props.extend(load_wave_fish(f))
2873 wave_fish = True
2874 elif ftype == 'pulsefish':
2875 eod_props.extend(load_pulse_fish(f))
2876 pulse_fish = True
2877 idx_offs = 0
2878 if wave_fish and not pulse_fish:
2879 idx_offs = sorted([ep['index'] for ep in eod_props])[0]
2880 # then load all other files:
2881 neods = len(eod_props)
2882 if neods < 1:
2883 neods = 1
2884 eod_props = [None]
2885 wave_eodfs = np.array([])
2886 wave_indices = np.array([])
2887 mean_eods = [None]*neods
2888 spec_data = [None]*neods
2889 peak_data = [None]*neods
2890 unit = None
2891 for f in file_pathes:
2892 recording, _, channel, _, ftype, idx, _ = parse_filename(f)
2893 if neods == 1 and idx > 0:
2894 idx = 0
2895 idx -= idx_offs
2896 if zf is not None:
2897 f = io.TextIOWrapper(zf.open(f, 'r'))
2898 if ftype == 'waveeodfs':
2899 wave_eodfs, wave_indices = load_wave_eodfs(f)
2900 elif ftype == 'eodwaveform':
2901 mean_eods[idx], unit = load_eod_waveform(f)
2902 elif ftype == 'wavespectrum':
2903 spec_data[idx], unit = load_wave_spectrum(f)
2904 elif ftype == 'pulsepeaks':
2905 peak_data[idx], unit = load_pulse_peaks(f)
2906 elif ftype == 'pulsetimes':
2907 pulse_times = load_pulse_times(f)
2908 eod_props[idx]['times'] = pulse_times
2909 eod_props[idx]['peaktimes'] = pulse_times
2910 elif ftype == 'pulsespectrum':
2911 spec_data[idx] = load_pulse_spectrum(f)
2912 # fix wave spectra:
2913 wave_eodfs = [fish.reshape(1, 2) if len(fish)>0 else fish
2914 for fish in wave_eodfs]
2915 if len(wave_eodfs) > 0 and len(spec_data) > 0:
2916 eodfs = []
2917 for idx, fish in zip(wave_indices, wave_eodfs):
2918 if idx >= 0:
2919 spec = spec_data[idx]
2920 specd = np.zeros((np.sum(np.isfinite(spec[:,-1])),
2921 2))
2922 specd[:,0] = spec[np.isfinite(spec[:,-1]),1]
2923 specd[:,1] = spec[np.isfinite(spec[:,-1]),-1]
2924 eodfs.append(specd)
2925 else:
2926 specd = np.zeros((10, 2))
2927 specd[:,0] = np.arange(len(specd))*fish[0,0]
2928 specd[:,1] = np.nan
2929 eodfs.append(specd)
2930 wave_eodfs = eodfs
2931 return mean_eods, wave_eodfs, wave_indices, eod_props, spec_data, \
2932 peak_data, recording, channel, unit
2935def load_recording(file_path, channel=0, load_kwargs={},
2936 eod_props=None, verbose=0):
2937 """Load recording.
2939 Parameters
2940 ----------
2941 file_path: string
2942 Full path of the file with the recorded data.
2943 Extension is optional. If absent, look for the first file
2944 with a reasonable extension.
2945 channel: int
2946 Channel of the recording to be returned.
2947 load_kwargs: dict
2948 Keyword arguments that are passed on to the
2949 format specific loading functions.
2950 eod_props: list of dict or None
2951 List of EOD properties from which start and end times of
2952 analysis window are extracted.
2953 verbose: int
2954 Verbosity level passed on to load function.
2956 Returns
2957 -------
2958 data: array of float
2959 Data of the requested `channel`.
2960 samplerate: float
2961 Sampling rate in Hertz.
2962 idx0: int
2963 Start index of the analysis window.
2964 idx1: int
2965 End index of the analysis window.
2966 data_file: str
2967 Full path and name of the loaded file inclusively extension.
2969 """
2970 data = None
2971 samplerate = 0.0
2972 idx0 = 0
2973 idx1 = 0
2974 data_file = ''
2975 if len(os.path.splitext(file_path)[1]) > 1:
2976 data_file = file_path
2977 else:
2978 data_files = glob.glob(file_path + os.extsep + '*')
2979 for dfile in data_files:
2980 if not os.path.splitext(dfile)[1][1:] in ['zip'] + list(TableData.ext_formats.values()):
2981 data_file = dfile
2982 break
2983 if os.path.exists(data_file):
2984 data, samplerate, unit, amax = load_data(data_file,
2985 verbose=verbose,
2986 **load_kwargs)
2987 idx0 = 0
2988 idx1 = len(data)
2989 if eod_props is not None and len(eod_props) > 0 and 'twin' in eod_props[0]:
2990 idx0 = int(eod_props[0]['twin']*samplerate)
2991 if len(eod_props) > 0 and 'window' in eod_props[0]:
2992 idx1 = idx0 + int(eod_props[0]['window']*samplerate)
2993 return data[:,channel], samplerate, idx0, idx1, data_file
2996def add_eod_analysis_config(cfg, thresh_fac=0.8, percentile=0.1,
2997 win_fac=2.0, min_win=0.01, max_eods=None,
2998 min_sem=False, unfilter_cutoff=0.0,
2999 flip_wave='none', flip_pulse='none',
3000 n_harm=10, min_pulse_win=0.001,
3001 peak_thresh_fac=0.01, min_dist=50.0e-6,
3002 width_frac = 0.5, fit_frac = 0.5,
3003 ipi_cv_thresh=0.5, ipi_percentile=30.0):
3004 """Add all parameters needed for the eod analysis functions as a new
3005 section to a configuration.
3007 Parameters
3008 ----------
3009 cfg: ConfigFile
3010 The configuration.
3012 See `eod_waveform()`, `analyze_wave()`, and `analyze_pulse()` for
3013 details on the remaining arguments.
3014 """
3015 cfg.add_section('EOD analysis:')
3016 cfg.add('eodSnippetFac', win_fac, '', 'The duration of EOD snippets is the EOD period times this factor.')
3017 cfg.add('eodMinSnippet', min_win, 's', 'Minimum duration of cut out EOD snippets.')
3018 cfg.add('eodMaxEODs', max_eods or 0, '', 'The maximum number of EODs used to compute the average EOD. If 0 use all EODs.')
3019 cfg.add('eodMinSem', min_sem, '', 'Use minimum of s.e.m. to set maximum number of EODs used to compute the average EOD.')
3020 cfg.add('unfilterCutoff', unfilter_cutoff, 'Hz', 'If non-zero remove effect of high-pass filter with this cut-off frequency.')
3021 cfg.add('flipWaveEOD', flip_wave, '', 'Flip EOD of wave fish to make largest extremum positive (flip, none, or auto).')
3022 cfg.add('flipPulseEOD', flip_pulse, '', 'Flip EOD of pulse fish to make the first large peak positive (flip, none, or auto).')
3023 cfg.add('eodHarmonics', n_harm, '', 'Number of harmonics fitted to the EOD waveform.')
3024 cfg.add('eodMinPulseSnippet', min_pulse_win, 's', 'Minimum duration of cut out EOD snippets for a pulse fish.')
3025 cfg.add('eodPeakThresholdFactor', peak_thresh_fac, '', 'Threshold for detection of peaks in pulse EODs as a fraction of the pulse amplitude.')
3026 cfg.add('eodMinimumDistance', min_dist, 's', 'Minimum distance between peaks and troughs in a EOD pulse.')
3027 cfg.add('eodPulseWidthFraction', width_frac, '', 'The width of a pulse is measured at this fraction of the pulse height.')
3028 cfg.add('eodExponentialFitFraction', fit_frac, '', 'An exponential function is fitted on the tail of a pulse starting at this fraction of the height of the last peak.')
3029 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.')
3030 cfg.add('ipiPercentile', ipi_percentile, '%', 'Use only interpulse intervals shorter than this percentile to compute EOD frequency.')
3033def eod_waveform_args(cfg):
3034 """Translates a configuration to the respective parameter names of
3035 the function `eod_waveform()`.
3037 The return value can then be passed as key-word arguments to this
3038 function.
3040 Parameters
3041 ----------
3042 cfg: ConfigFile
3043 The configuration.
3045 Returns
3046 -------
3047 a: dict
3048 Dictionary with names of arguments of the `eod_waveform()` function
3049 and their values as supplied by `cfg`.
3050 """
3051 a = cfg.map({'win_fac': 'eodSnippetFac',
3052 'min_win': 'eodMinSnippet',
3053 'max_eods': 'eodMaxEODs',
3054 'min_sem': 'eodMinSem',
3055 'unfilter_cutoff': 'unfilterCutoff'})
3056 return a
3059def analyze_wave_args(cfg):
3060 """Translates a configuration to the respective parameter names of
3061 the function `analyze_wave()`.
3063 The return value can then be passed as key-word arguments to this
3064 function.
3066 Parameters
3067 ----------
3068 cfg: ConfigFile
3069 The configuration.
3071 Returns
3072 -------
3073 a: dict
3074 Dictionary with names of arguments of the `analyze_wave()` function
3075 and their values as supplied by `cfg`.
3076 """
3077 a = cfg.map({'n_harm': 'eodHarmonics',
3078 'power_n_harmonics': 'powerNHarmonics',
3079 'flip_wave': 'flipWaveEOD'})
3080 return a
3083def analyze_pulse_args(cfg):
3084 """Translates a configuration to the respective parameter names of
3085 the function `analyze_pulse()`.
3087 The return value can then be passed as key-word arguments to this
3088 function.
3090 Parameters
3091 ----------
3092 cfg: ConfigFile
3093 The configuration.
3095 Returns
3096 -------
3097 a: dict
3098 Dictionary with names of arguments of the `analyze_pulse()` function
3099 and their values as supplied by `cfg`.
3100 """
3101 a = cfg.map({'min_pulse_win': 'eodMinPulseSnippet',
3102 'peak_thresh_fac': 'eodPeakThresholdFactor',
3103 'min_dist': 'eodMinimumDistance',
3104 'width_frac': 'eodPulseWidthFraction',
3105 'fit_frac': 'eodExponentialFitFraction',
3106 'flip_pulse': 'flipPulseEOD',
3107 'ipi_cv_thresh': 'ipiCVThresh',
3108 'ipi_percentile': 'ipiPercentile'})
3109 return a
3112def add_species_config(cfg, species_file='none', wave_max_rms=0.2,
3113 pulse_max_rms=0.2):
3114 """Add parameters needed for assigning EOD waveforms to species.
3116 Parameters
3117 ----------
3118 cfg: ConfigFile
3119 The configuration.
3120 species_file: string
3121 File path to a file containing species names and corresponding
3122 file names of EOD waveform templates. If 'none', no species
3123 assignemnt is performed.
3124 wave_max_rms: float
3125 Maximum allowed rms difference (relative to standard deviation
3126 of EOD waveform) to an EOD waveform template for assignment to
3127 a wave fish species.
3128 pulse_max_rms: float
3129 Maximum allowed rms difference (relative to standard deviation
3130 of EOD waveform) to an EOD waveform template for assignment to
3131 a pulse fish species.
3132 """
3133 cfg.add_section('Species assignment:')
3134 cfg.add('speciesFile', species_file, '', 'File path to a file containing species names and corresponding file names of EOD waveform templates.')
3135 cfg.add('maximumWaveSpeciesRMS', wave_max_rms, '', 'Maximum allowed rms difference (relative to standard deviation of EOD waveform) to an EOD waveform template for assignment to a wave fish species.')
3136 cfg.add('maximumPulseSpeciesRMS', pulse_max_rms, '', 'Maximum allowed rms difference (relative to standard deviation of EOD waveform) to an EOD waveform template for assignment to a pulse fish species.')
3139def add_eod_quality_config(cfg, max_clipped_frac=0.1, max_variance=0.0,
3140 max_rms_error=0.05, min_power=-100.0, max_thd=0.0,
3141 max_crossings=4, max_relampl_harm1=0.0,
3142 max_relampl_harm2=0.0, max_relampl_harm3=0.0):
3143 """Add parameters needed for assesing the quality of an EOD waveform.
3145 Parameters
3146 ----------
3147 cfg: ConfigFile
3148 The configuration.
3150 See `wave_quality()` and `pulse_quality()` for details on
3151 the remaining arguments.
3152 """
3153 cfg.add_section('Waveform selection:')
3154 cfg.add('maximumClippedFraction', max_clipped_frac, '', 'Take waveform of the fish with the highest power only if the fraction of clipped signals is below this value.')
3155 cfg.add('maximumVariance', max_variance, '', 'Skip waveform of fish if the standard error of the EOD waveform relative to the peak-to-peak amplitude is larger than this number. A value of zero allows any variance.')
3156 cfg.add('maximumRMSError', max_rms_error, '', 'Skip waveform of wave fish if the root-mean-squared error of the fit relative to the peak-to-peak amplitude is larger than this number.')
3157 cfg.add('minimumPower', min_power, 'dB', 'Skip waveform of wave fish if its power is smaller than this value.')
3158 cfg.add('maximumTotalHarmonicDistortion', max_thd, '', 'Skip waveform of wave fish if its total harmonic distortion is larger than this value. If set to zero do not check.')
3159 cfg.add('maximumCrossings', max_crossings, '', 'Maximum number of zero crossings per EOD period.')
3160 cfg.add('maximumFirstHarmonicAmplitude', max_relampl_harm1, '', 'Skip waveform of wave fish if the amplitude of the first harmonic is higher than this factor times the amplitude of the fundamental. If set to zero do not check.')
3161 cfg.add('maximumSecondHarmonicAmplitude', max_relampl_harm2, '', 'Skip waveform of wave fish if the ampltude of the second harmonic is higher than this factor times the amplitude of the fundamental. That is, the waveform appears to have twice the frequency than the fundamental. If set to zero do not check.')
3162 cfg.add('maximumThirdHarmonicAmplitude', max_relampl_harm3, '', 'Skip waveform of wave fish if the ampltude of the third harmonic is higher than this factor times the amplitude of the fundamental. If set to zero do not check.')
3165def wave_quality_args(cfg):
3166 """Translates a configuration to the respective parameter names of
3167 the function `wave_quality()`.
3169 The return value can then be passed as key-word arguments to this
3170 function.
3172 Parameters
3173 ----------
3174 cfg: ConfigFile
3175 The configuration.
3177 Returns
3178 -------
3179 a: dict
3180 Dictionary with names of arguments of the `wave_quality()` function
3181 and their values as supplied by `cfg`.
3182 """
3183 a = cfg.map({'max_clipped_frac': 'maximumClippedFraction',
3184 'max_rms_sem': 'maximumVariance',
3185 'max_rms_error': 'maximumRMSError',
3186 'min_power': 'minimumPower',
3187 'max_crossings': 'maximumCrossings',
3188 'min_freq': 'minimumFrequency',
3189 'max_freq': 'maximumFrequency',
3190 'max_thd': 'maximumTotalHarmonicDistortion',
3191 'max_db_diff': 'maximumPowerDifference',
3192 'max_harmonics_db': 'maximumHarmonicsPower',
3193 'max_relampl_harm1': 'maximumFirstHarmonicAmplitude',
3194 'max_relampl_harm2': 'maximumSecondHarmonicAmplitude',
3195 'max_relampl_harm3': 'maximumThirdHarmonicAmplitude'})
3196 return a
3199def pulse_quality_args(cfg):
3200 """Translates a configuration to the respective parameter names of
3201 the function `pulse_quality()`.
3203 The return value can then be passed as key-word arguments to this
3204 function.
3206 Parameters
3207 ----------
3208 cfg: ConfigFile
3209 The configuration.
3211 Returns
3212 -------
3213 a: dict
3214 Dictionary with names of arguments of the `pulse_quality()` function
3215 and their values as supplied by `cfg`.
3216 """
3217 a = cfg.map({'max_clipped_frac': 'maximumClippedFraction',
3218 'max_rms_sem': 'maximumRMSNoise'})
3219 return a
3222def main():
3223 import matplotlib.pyplot as plt
3224 from .fakefish import pulsefish_eods
3226 print('Analysis of EOD waveforms.')
3228 # data:
3229 samplerate = 44100.0
3230 data = pulsefish_eods('Triphasic', 83.0, samplerate, 5.0, noise_std=0.02)
3231 unit = 'mV'
3232 eod_idx, _ = detect_peaks(data, 1.0)
3233 eod_times = eod_idx/samplerate
3235 # analyse EOD:
3236 mean_eod, eod_times = eod_waveform(data, samplerate, eod_times)
3237 mean_eod, props, peaks, power = analyze_pulse(mean_eod, eod_times)
3239 # plot:
3240 fig, axs = plt.subplots(1, 2)
3241 plot_eod_waveform(axs[0], mean_eod, props, peaks, unit=unit)
3242 axs[0].set_title('{type} fish: EODf = {EODf:.1f} Hz'.format(**props))
3243 plot_pulse_spectrum(axs[1], power, props)
3244 plt.show()
3247if __name__ == '__main__':
3248 main()