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