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