Coverage for src/thunderfish/thunderfish.py: 0%
921 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-03-23 22:57 +0000
1"""# thunderfish
3Automatically detect and analyze all EOD waveforms in short recordings
4and generated summary plots and data tables.
6Run it from the thunderfish development directory as:
7```sh
8> python3 -m thunderfish.thunderfish audiofile.wav
9```
10Or install thunderfish
11```sh
12> sudo pip3 install .
13```
14Then you can run it directly from every directory:
15```sh
16> thunderfish audiofile.wav
17```
18"""
20import sys
21import os
22import glob
23import io
24import zipfile
25import argparse
26import traceback
27import numpy as np
28import matplotlib.pyplot as plt
29import matplotlib.gridspec as gridspec
30import matplotlib.ticker as ticker
31import matplotlib.lines as ml
32from matplotlib.transforms import Bbox
33from matplotlib.backends.backend_pdf import PdfPages
34from multiprocessing import Pool, freeze_support, cpu_count
35from audioio import play, fade, load_audio
36from thunderlab.configfile import ConfigFile
37from thunderlab.dataloader import load_data
38from thunderlab.powerspectrum import decibel, plot_decibel_psd, multi_psd
39from thunderlab.powerspectrum import add_multi_psd_config, multi_psd_args
40from thunderlab.tabledata import TableData, add_write_table_config, write_table_args
41from .version import __version__, __year__
42from .bestwindow import add_clip_config, add_best_window_config
43from .bestwindow import clip_args, best_window_args
44from .bestwindow import analysis_window, plot_data_window
45from .checkpulse import check_pulse, add_check_pulse_config, check_pulse_args
46from .pulses import extract_pulsefish
47from .harmonics import add_psd_peak_detection_config, add_harmonic_groups_config
48from .harmonics import harmonic_groups, harmonic_groups_args, psd_peak_detection_args
49from .harmonics import colors_markers, plot_harmonic_groups
50from .consistentfishes import consistent_fishes
51from .eodanalysis import eod_waveform, waveeod_waveform, analyze_wave, analyze_pulse
52from .eodanalysis import clipped_fraction
53from .eodanalysis import plot_eod_recording, plot_pulse_eods
54from .eodanalysis import plot_eod_waveform, plot_eod_snippets
55from .eodanalysis import plot_pulse_spectrum, plot_wave_spectrum
56from .eodanalysis import add_eod_analysis_config, eod_waveform_args
57from .eodanalysis import analyze_wave_args, analyze_pulse_args
58from .eodanalysis import add_species_config
59from .eodanalysis import wave_quality, wave_quality_args, add_eod_quality_config
60from .eodanalysis import pulse_quality, pulse_quality_args
61from .eodanalysis import save_eod_waveform, save_wave_eodfs, save_wave_fish, save_pulse_fish
62from .eodanalysis import save_wave_spectrum, save_pulse_spectrum, save_pulse_peaks, save_pulse_times
63from .eodanalysis import load_eod_waveform, load_wave_eodfs, load_wave_fish, load_pulse_fish
64from .eodanalysis import load_wave_spectrum, load_pulse_spectrum, load_pulse_peaks
65from .eodanalysis import save_analysis, load_analysis, load_recording
66from .eodanalysis import parse_filename, file_types
69# Colors taken from https://github.com/bendalab/plottools/blob/master/src/plottools/colors.py (colors_vivid)
71trace_color = '#D71000'
72""" Default color for traces."""
74data_color = '#1040C0'
75""" Default color for the full data trace."""
77fit_color = '#E09080'
78""" Default color for fits."""
80spectrum_color = '#1040C0'
81""" Default color for power spectra."""
83data_styles = dict(dstyle=dict(color=data_color, lw=1.4),
84 wstyle=dict(color=trace_color, lw=1.4))
86trace_style = dict(color=trace_color, lw=2)
88spectrum_style = dict(color=spectrum_color, lw=2)
90eod_styles = dict(mstyle=dict(color=trace_color, lw=3),
91 pstyle=dict(facecolor='#00A050', alpha=0.2,
92 edgecolor='none'),
93 nstyle=dict(facecolor='#1040C0', alpha=0.2,
94 edgecolor='none'),
95 sstyle=dict(color='0.8'),
96 fstyle=dict(color=fit_color, lw=6),
97 zstyle=dict(color='0.3', lw=1))
99snippet_style = dict(scaley=False, lw=0.5, color='0.6')
101wave_spec_styles = dict(mstyle=dict(color=spectrum_color, markersize=10),
102 sstyle=dict(color=spectrum_color, alpha=0.5, lw=2))
104pulse_spec_styles = dict(sstyle=dict(color=spectrum_color, lw=3),
105 pstyle=dict(ls='', marker='o', markersize=12,
106 color=spectrum_color, mec='none', mew=0,
107 alpha=0.6),
108 cstyle=dict(ls='-', color='0.5', lw=1),
109 att5_color='0.8',
110 att50_color='0.9')
112def configuration():
113 """Assemble configuration parameter for thunderfish.
115 Returns
116 -------
117 cfg: ConfigFile
118 Configuration parameters.
119 """
120 cfg = ConfigFile()
121 add_multi_psd_config(cfg)
122 cfg.add('frequencyThreshold', 1.0, 'Hz',
123 'The fundamental frequency of each fish needs to be detected in each power spectrum within this threshold.')
124 # TODO: make this threshold dependent on frequency resolution!
125 cfg.add('minPSDAverages', 3, '', 'Minimum number of fft averages for estimating the power spectrum.') # needed by fishfinder
126 add_psd_peak_detection_config(cfg)
127 add_harmonic_groups_config(cfg)
128 add_clip_config(cfg)
129 cfg.add('unwrapData', False, '', 'Unwrap clipped voltage traces.')
130 add_best_window_config(cfg, win_size=8.0, w_cv_ampl=10.0)
131 add_check_pulse_config(cfg)
132 add_eod_analysis_config(cfg, min_pulse_win=0.004)
133 del cfg['eodSnippetFac']
134 del cfg['eodMinSnippet']
135 del cfg['eodMinSem']
136 add_eod_quality_config(cfg)
137 add_species_config(cfg)
138 add_write_table_config(cfg, table_format='csv', unit_style='row',
139 align_columns=True, shrink_width=False)
140 return cfg
143def save_configuration(cfg, config_file):
144 """Save configuration parameter for thunderfish to a file.
146 Parameters
147 ----------
148 cfg: ConfigFile
149 Configuration parameters and their values.
150 config_file: string
151 Name of the configuration file to be loaded.
152 """
153 ext = os.path.splitext(config_file)[1]
154 if ext != os.extsep + 'cfg':
155 print('configuration file name must have .cfg as extension!')
156 else:
157 print('write configuration to %s ...' % config_file)
158 del cfg['fileColumnNumbers']
159 del cfg['fileShrinkColumnWidth']
160 del cfg['fileMissing']
161 del cfg['fileLaTeXLabelCommand']
162 del cfg['fileLaTeXMergeStd']
163 cfg.dump(config_file)
166def detect_eods(data, rate, min_clip, max_clip, name, mode,
167 verbose, plot_level, cfg):
168 """Detect EODs of all fish present in the data.
170 Parameters
171 ----------
172 data: array of floats
173 The recording in which to detect EODs.
174 rate: float
175 Sampling rate of the dataset.
176 min_clip: float
177 Minimum amplitude that is not clipped.
178 max_clip: float
179 Maximum amplitude that is not clipped.
180 name: string
181 Name of the recording (e.g. its filename).
182 mode: string
183 Characters in the string indicate what and how to analyze:
184 - 'w': analyze wavefish
185 - 'p': analyze pulsefish
186 - 'P': analyze only the pulsefish with the largest amplitude (not implemented yet)
187 verbose: int
188 Print out information about EOD detection if greater than zero.
189 plot_level : int
190 Similar to verbosity levels, but with plots.
191 cfg: ConfigFile
192 Configuration parameters.
194 Returns
195 -------
196 psd_data: list of 2D arrays
197 List of power spectra (frequencies and power) of the analysed data
198 for different frequency resolutions.
199 wave_eodfs: list of 2D arrays
200 Frequency and power of fundamental frequency/harmonics of all wave fish.
201 wave_indices: array of int
202 Indices of wave fish mapping from wave_eodfs to eod_props.
203 If negative, then that EOD frequency has no waveform described in eod_props.
204 eod_props: list of dict
205 Lists of EOD properties as returned by analyze_pulse() and analyze_wave()
206 for each waveform in mean_eods.
207 mean_eods: list of 2-D arrays with time, mean, sem, and fit.
208 Averaged EOD waveforms of pulse and wave fish.
209 spec_data: list of 2_D arrays
210 For each pulsefish a power spectrum of the single pulse and for
211 each wavefish the relative amplitudes and phases of the harmonics.
212 peak_data: list of 2_D arrays
213 For each pulse fish a list of peak properties
214 (index, time, and amplitude), empty array for wave fish.
215 power_thresh: 2 D array or None
216 Frequency (first column) and power (second column) of threshold
217 derived from single pulse spectra to discard false wave fish.
218 None if no pulse fish was detected.
219 skip_reason: list of string
220 Reasons, why an EOD was discarded.
221 """
222 dfreq = np.nan
223 nfft = 0
224 psd_data = [[]]
225 wave_eodfs = []
226 wave_indices = []
227 if 'w' in mode:
228 # detect wave fish:
229 psd_data = multi_psd(data, rate, **multi_psd_args(cfg))
230 dfreq = np.mean(np.diff(psd_data[0][:,0]))
231 nfft = int(rate/dfreq)
232 h_kwargs = psd_peak_detection_args(cfg)
233 h_kwargs.update(harmonic_groups_args(cfg))
234 wave_eodfs_list = []
235 for i, psd in enumerate(psd_data):
236 wave_eodfs = harmonic_groups(psd[:,0], psd[:,1], verbose-1, **h_kwargs)[0]
237 if verbose > 0 and len(psd_data) > 1:
238 numpsdresolutions = cfg.value('numberPSDResolutions')
239 print('fundamental frequencies detected in power spectrum of window %d at resolution %d:'
240 % (i//numpsdresolutions, i%numpsdresolutions))
241 if len(wave_eodfs) > 0:
242 print(' ' + ' '.join(['%.1f' % freq[0, 0] for freq in wave_eodfs]))
243 else:
244 print(' none')
245 wave_eodfs_list.append(wave_eodfs)
246 wave_eodfs = consistent_fishes(wave_eodfs_list,
247 df_th=cfg.value('frequencyThreshold'))
248 if verbose > 0:
249 if len(wave_eodfs) > 0:
250 print('found %2d EOD frequencies consistent in all power spectra:' % len(wave_eodfs))
251 print(' ' + ' '.join(['%.1f' % freq[0, 0] for freq in wave_eodfs]))
252 else:
253 print('no fundamental frequencies are consistent in all power spectra')
255 # analysis results:
256 eod_props = []
257 mean_eods = []
258 spec_data = []
259 peak_data = []
260 power_thresh = None
261 skip_reason = []
262 max_pulse_amplitude = 0.0
263 zoom_window = []
265 if 'p' in mode:
266 # detect pulse fish:
267 _, eod_times, eod_peaktimes, zoom_window, _ = \
268 extract_pulsefish(data, rate, max_clip, verbose=verbose-1,
269 plot_level=plot_level,
270 save_path=os.path.splitext(os.path.basename(name))[0])
272 #eod_times = []
273 #eod_peaktimes = []
274 if verbose > 0:
275 if len(eod_times) > 0:
276 print('found %2d pulsefish EODs' % len(eod_times))
277 else:
278 print('no pulsefish EODs found')
280 # analyse eod waveform of pulse-fish:
281 min_freq_res = cfg.value('frequencyResolution')
282 for k, (eod_ts, eod_pts) in enumerate(zip(eod_times, eod_peaktimes)):
283 mean_eod, eod_times0 = \
284 eod_waveform(data, rate, eod_ts, win_fac=0.8,
285 min_win=cfg.value('eodMinPulseSnippet'),
286 min_sem=False, **eod_waveform_args(cfg))
287 mean_eod, props, peaks, power = analyze_pulse(mean_eod, eod_times0,
288 freq_resolution=min_freq_res,
289 **analyze_pulse_args(cfg))
290 if len(peaks) == 0:
291 print('error: no peaks in pulse EOD detected')
292 continue
293 clipped_frac = clipped_fraction(data, rate, eod_times0,
294 mean_eod, min_clip, max_clip)
295 props['peaktimes'] = eod_pts # XXX that should go into analyze pulse
296 props['index'] = len(eod_props)
297 props['clipped'] = clipped_frac
298 props['samplerate'] = rate
299 props['nfft'] = nfft
300 props['dfreq'] = dfreq
302 # add good waveforms only:
303 skips, msg, skipped_clipped = pulse_quality(props, **pulse_quality_args(cfg))
305 if len(skips) == 0:
306 eod_props.append(props)
307 mean_eods.append(mean_eod)
308 spec_data.append(power)
309 peak_data.append(peaks)
310 if verbose > 0:
311 print('take %6.1fHz pulse fish: %s' % (props['EODf'], msg))
312 else:
313 skip_reason += ['%.1fHz pulse fish %s' % (props['EODf'], skips)]
314 if verbose > 0:
315 print('skip %6.1fHz pulse fish: %s (%s)' %
316 (props['EODf'], skips, msg))
318 # threshold for wave fish peaks based on single pulse spectra:
319 if len(skips) == 0 or skipped_clipped:
320 if max_pulse_amplitude < props['p-p-amplitude']:
321 max_pulse_amplitude = props['p-p-amplitude']
322 i0 = np.argmin(np.abs(mean_eod[:,0]))
323 i1 = len(mean_eod) - i0
324 pulse_data = np.zeros(len(data))
325 for t in props['peaktimes']:
326 idx = int(t*rate)
327 ii0 = i0 if idx-i0 >= 0 else idx
328 ii1 = i1 if idx+i1 < len(pulse_data) else len(pulse_data)-1-idx
329 pulse_data[idx-ii0:idx+ii1] = mean_eod[i0-ii0:i0+ii1,1]
330 pulse_psd = multi_psd(pulse_data, rate, **multi_psd_args(cfg))
331 pulse_power = pulse_psd[0][:,1]
332 pulse_power *= len(data)/rate/props['period']/len(props['peaktimes'])
333 pulse_power *= 5.0
334 if power_thresh is None:
335 power_thresh = pulse_psd[0]
336 power_thresh[:,1] = pulse_power
337 else:
338 power_thresh[:,1] += pulse_power
340 # remove wavefish below pulse fish power:
341 if 'w' in mode and power_thresh is not None:
342 n = len(wave_eodfs)
343 maxh = 3 # XXX make parameter
344 df = power_thresh[1,0] - power_thresh[0,0]
345 for k, fish in enumerate(reversed(wave_eodfs)):
346 idx = np.array(fish[:maxh,0]//df, dtype=int)
347 for offs in range(-2, 3):
348 nbelow = np.sum(fish[:maxh,1] < power_thresh[idx+offs,1])
349 if nbelow > 0:
350 wave_eodfs.pop(n-1-k)
351 if verbose > 0:
352 print('skip %6.1fHz wave fish: %2d harmonics are below pulsefish threshold' % (fish[0,0], nbelow))
353 break
355 if 'w' in mode:
356 # analyse EOD waveform of all wavefish:
357 powers = np.array([np.sum(fish[:,1]) for fish in wave_eodfs])
358 power_indices = np.argsort(-powers)
359 wave_indices = np.zeros(len(wave_eodfs), dtype=int) - 3
360 for k, idx in enumerate(power_indices):
361 fish = wave_eodfs[idx]
362 """
363 eod_times = np.arange(0.0, len(data)/rate, 1.0/fish[0,0])
364 mean_eod, eod_times = \
365 eod_waveform(data, rate, eod_times, win_fac=3.0, min_win=0.0,
366 min_sem=(k==0), **eod_waveform_args(cfg))
367 """
368 mean_eod, eod_times = \
369 waveeod_waveform(data, rate, fish[0, 0], win_fac=3.0,
370 unfilter_cutoff=cfg.value('unfilterCutoff'))
371 if len(mean_eod) == 0:
372 continue
373 mean_eod, props, sdata, error_str = \
374 analyze_wave(mean_eod, fish, **analyze_wave_args(cfg))
375 if error_str:
376 print(name + ': ' + error_str)
377 clipped_frac = clipped_fraction(data, rate, eod_times,
378 mean_eod, min_clip, max_clip)
379 props['n'] = len(eod_times)
380 props['index'] = len(eod_props)
381 props['clipped'] = clipped_frac
382 props['samplerate'] = rate
383 props['nfft'] = nfft
384 props['dfreq'] = dfreq
385 # remove wave fish that are smaller than the largest pulse fish:
386 if props['p-p-amplitude'] < 0.01*max_pulse_amplitude:
387 rm_indices = power_indices[k:]
388 if verbose > 0:
389 print('skip %6.1fHz wave fish: power=%5.1fdB, p-p amplitude=%5.1fdB smaller than pulse fish=%5.1dB - 20dB' %
390 (props['EODf'], decibel(powers[idx]),
391 decibel(props['p-p-amplitude']), decibel(max_pulse_amplitude)))
392 for idx in rm_indices[1:]:
393 print('skip %6.1fHz wave fish: power=%5.1fdB even smaller' %
394 (wave_eodfs[idx][0,0], decibel(powers[idx])))
395 wave_eodfs = [eodfs for idx, eodfs in enumerate(wave_eodfs)
396 if idx not in rm_indices]
397 wave_indices = np.array([idcs for idx, idcs in enumerate(wave_indices)
398 if idx not in rm_indices], dtype=int)
399 break
400 # add good waveforms only:
401 remove, skips, msg = wave_quality(props, sdata[1:,3], **wave_quality_args(cfg))
402 if len(skips) == 0:
403 wave_indices[idx] = props['index']
404 eod_props.append(props)
405 mean_eods.append(mean_eod)
406 spec_data.append(sdata)
407 peak_data.append([])
408 if verbose > 0:
409 print('take %6.1fHz wave fish: %s' % (props['EODf'], msg))
410 else:
411 wave_indices[idx] = -2 if remove else -1
412 skip_reason += ['%.1fHz wave fish %s' % (props['EODf'], skips)]
413 if verbose > 0:
414 print('%-6s %6.1fHz wave fish: %s (%s)' %
415 ('remove' if remove else 'skip', props['EODf'], skips, msg))
416 wave_eodfs = [eodfs for idx, eodfs in zip(wave_indices, wave_eodfs) if idx > -2]
417 wave_indices = np.array([idx for idx in wave_indices if idx > -2], dtype=int)
418 return (psd_data, wave_eodfs, wave_indices, eod_props, mean_eods,
419 spec_data, peak_data, power_thresh, skip_reason, zoom_window)
422def remove_eod_files(output_basename, verbose, cfg):
423 """Remove all files from previous runs of thunderfish
424 """
425 ff = cfg.value('fileFormat')
426 if ff == 'py':
427 fext = 'py'
428 else:
429 fext = TableData.extensions[ff]
430 # remove all files from previous runs of thunderfish:
431 for fn in glob.glob('%s*.%s' % (output_basename, fext)):
432 os.remove(fn)
433 if verbose > 0:
434 print('removed file %s' % fn)
437def plot_style():
438 """Set style of plots.
439 """
440 plt.rcParams['figure.facecolor'] = 'white'
441 plt.rcParams['axes.facecolor'] = 'none'
442 plt.rcParams['xtick.direction'] = 'out'
443 plt.rcParams['ytick.direction'] = 'out'
446def axes_style(ax):
447 """Fix style of axes.
449 Parameters
450 ----------
451 ax: matplotlib axes
452 """
453 ax.spines['top'].set_visible(False)
454 ax.spines['right'].set_visible(False)
455 ax.get_xaxis().tick_bottom()
456 ax.get_yaxis().tick_left()
459def plot_eods(base_name, message_filename,
460 raw_data, rate, channel, idx0, idx1, clipped,
461 psd_data, wave_eodfs, wave_indices, mean_eods, eod_props,
462 peak_data, spec_data, indices, unit, zoom_window,
463 n_snippets=10, power_thresh=None, label_power=True,
464 all_eods=False, spec_plots='auto', skip_bad=True,
465 log_freq=False, min_freq=0.0, max_freq=3000.0,
466 interactive=True, verbose=0):
467 """Creates an output plot for the thunderfish program.
469 This output contains the raw trace where the analysis window is
470 marked, the power-spectrum of this analysis window where the
471 detected fish are marked, plots of averaged EOD plots, and
472 spectra of the EOD waveforms.
474 Parameters
475 ----------
476 base_name: string
477 Basename of audio_file.
478 message_filename: string or None
479 Path to meta-data message.
480 raw_data: array
481 Dataset.
482 rate: float
483 Sampling rate of the dataset.
484 channel: int or None
485 Channel of the recording to be put into the plot title.
486 If None, do not write the channel into the title.
487 idx0: float
488 Index of the beginning of the analysis window in the dataset.
489 idx1: float
490 Index of the end of the analysis window in the dataset.
491 clipped: float
492 Fraction of clipped amplitudes.
493 psd_data: 2D array
494 Power spectrum (frequencies and power) of the analysed data.
495 wave_eodfs: array
496 Frequency and power of fundamental frequency/harmonics of several fish.
497 wave_indices: array of int
498 Indices of wave fish mapping from wave_eodfs to eod_props.
499 If negative, then that EOD frequency has no waveform described in eod_props.
500 mean_eods: list of 2-D arrays with time, mean and std.
501 Mean trace for the mean EOD plot.
502 eod_props: list of dict
503 Properties for each waveform in mean_eods.
504 peak_data: list of 2_D arrays
505 For each pulsefish a list of peak properties
506 (index, time, and amplitude).
507 spec_data: list of 2_D arrays
508 For each pulsefish a power spectrum of the single pulse and for
509 each wavefish the relative amplitudes and phases of the harmonics.
510 indices: list of int or None
511 Indices of the fish in eod_props to be plotted.
512 If None try to plot all.
513 unit: string
514 Unit of the trace and the mean EOD.
515 n_snippets: int
516 Number of EOD waveform snippets to be plotted. If zero do not plot any.
517 power_thresh: 2 D array or None
518 Frequency (first column) and power (second column) of threshold
519 derived from single pulse spectra to discard false wave fish.
520 label_power: boolean
521 If `True` put the power in decibel in addition to the frequency
522 into the legend.
523 all_eods: bool
524 Plot all EOD waveforms.
525 spec_plots: bool or 'auto'
526 Plot amplitude spectra of EOD waveforms.
527 If 'auto', plot them if there is a singel waveform only.
528 skip_bad: bool
529 Skip harmonic groups without index (entry in indices is negative).
530 log_freq: boolean
531 Logarithmic (True) or linear (False) frequency axis of power spectrum of recording.
532 min_freq: float
533 Limits of frequency axis of power spectrum of recording
534 are set to `(min_freq, max_freq)` if `max_freq` is greater than zero
535 max_freq: float
536 Limits of frequency axis of power spectrum of recording
537 are set to `(min_freq, max_freq)` and limits of power axis are computed
538 from powers below max_freq if `max_freq` is greater than zero
539 interactive: bool
540 If True install some keyboard interaction.
541 verbose: int
542 Print out information about data to be plotted if greater than zero.
544 Returns
545 -------
546 fig: plt.figure
547 Figure with the plots.
548 """
550 def keypress(event):
551 if event.key in 'pP':
552 if idx1 > idx0:
553 playdata = 1.0 * raw_data[idx0:idx1]
554 else:
555 playdata = 1.0 * raw_data[:]
556 fade(playdata, rate, 0.1)
557 play(playdata, rate, blocking=False)
558 if event.key in 'mM' and message_filename:
559 # play voice message:
560 msg, msg_rate = load_audio(message_filename)
561 play(msg, msg_rate, blocking=False)
563 def recording_format_coord(x, y):
564 return 'full recording: x=%.3f s, y=%.3f' % (x, y)
566 def recordingzoom_format_coord(x, y):
567 return 'recording zoom-in: x=%.3f s, y=%.3f' % (x, y)
569 def psd_format_coord(x, y):
570 return 'power spectrum: x=%.1f Hz, y=%.1f dB' % (x, y)
572 def meaneod_format_coord(x, y):
573 return 'mean EOD waveform: x=%.2f ms, y=%.3f' % (x, y)
575 def ampl_format_coord(x, y):
576 return u'amplitude spectrum: x=%.0f, y=%.2f' % (x, y)
578 def phase_format_coord(x, y):
579 return u'phase spectrum: x=%.0f, y=%.2f \u03c0' % (x, y/np.pi)
581 def pulsepsd_format_coord(x, y):
582 return 'single pulse power spectrum: x=%.1f Hz, y=%.1f dB' % (x, y)
584 # count number of fish types to be plotted:
585 if indices is None:
586 indices = np.arange(len(eod_props))
587 else:
588 indices = np.array(indices, dtype=int)
589 nwave = 0
590 npulse = 0
591 for idx in indices:
592 if eod_props[idx]['type'] == 'pulse':
593 npulse += 1
594 elif eod_props[idx]['type'] == 'wave':
595 nwave += 1
596 neods = nwave + npulse
598 if verbose > 0:
599 print('plot: %2d waveforms: %2d wave fish, %2d pulse fish and %2d EOD frequencies.'
600 % (len(indices), nwave, npulse, len(wave_eodfs)))
602 # size and positions:
603 if spec_plots == 'auto':
604 spec_plots = len(indices) == 1
605 large_plots = spec_plots or len(indices) <= 2
606 width = 14.0
607 height = 10.0
608 if all_eods and len(indices) > 0:
609 nrows = len(indices) if spec_plots else (len(indices)+1)//2
610 if large_plots:
611 height = 6.0 + 4.0*nrows
612 else:
613 height = 6.4 + 1.9*nrows
614 leftx = 1.0/width
615 midx = 0.5 + leftx
616 fullwidth = 1.0-1.4/width
617 halfwidth = 0.5-1.4/width
618 pheight = 3.0/height
620 # figure:
621 plot_style()
622 fig = plt.figure(figsize=(width, height))
623 if interactive:
624 fig.canvas.mpl_connect('key_press_event', keypress)
626 # plot title:
627 title = base_name
628 if channel is not None:
629 title += ' c%d' % channel
630 ax = fig.add_axes([0.2/width, 1.0-0.6/height, 1.0-0.4/width, 0.55/height])
631 ax.text(0.0, 1.0, title, fontsize=22, va='top')
632 ax.text(1.0, 1.0, 'thunderfish by Benda-Lab', fontsize=16, ha='right', va='top')
633 ax.text(1.0, 0.0, 'version %s' % __version__, fontsize=16, ha='right', va='bottom')
634 ax.set_frame_on(False)
635 ax.set_axis_off()
636 ax.set_navigate(False)
638 # layout of recording and psd plots:
639 #force_both = True # set to True for debugging pulse and wave detection
640 force_both = False
641 posy = 1.0 - 4.0/height
642 axr = None
643 axp = None
644 legend_inside = True
645 legendwidth = 2.2/width if label_power else 1.7/width
646 if neods == 0:
647 axr = fig.add_axes([leftx, posy, fullwidth, pheight]) # top, wide
648 if len(psd_data) > 0:
649 axp = fig.add_axes([leftx, 2.0/height, fullwidth, pheight]) # bottom, wide
650 else:
651 if npulse == 0 and nwave > 2 and psd_data is not None and \
652 len(psd_data) > 0 and not force_both:
653 axp = fig.add_axes([leftx, posy, fullwidth-legendwidth, pheight]) # top, wide
654 legend_inside = False
655 elif (npulse > 0 or psd_data is None or len(psd_data) == 0) \
656 and len(wave_eodfs) == 0 and not force_both:
657 axr = fig.add_axes([leftx, posy, fullwidth, pheight]) # top, wide
658 else:
659 axr = fig.add_axes([leftx, posy, halfwidth, pheight]) # top left
660 label_power = False
661 legendwidth = 2.2/width
662 axp = fig.add_axes([midx, posy, halfwidth, pheight]) # top, right
664 # best window data:
665 data = raw_data[idx0:idx1] if idx1 > idx0 else raw_data
667 # plot recording
668 pulse_colors, pulse_markers = colors_markers()
669 pulse_colors = pulse_colors[3:]
670 pulse_markers = pulse_markers[3:]
671 if axr is not None:
672 axes_style(axr)
673 twidth = 0.1
674 if len(indices) > 0:
675 if eod_props[indices[0]]['type'] == 'wave':
676 twidth = 5.0/eod_props[indices[0]]['EODf']
677 else:
678 if len(wave_eodfs) > 0:
679 twidth = 3.0/eod_props[indices[0]]['EODf']
680 else:
681 twidth = 10.0/eod_props[indices[0]]['EODf']
682 twidth = (1+twidth//0.005)*0.005
683 if data is not None and len(data) > 0:
684 plot_eod_recording(axr, data, rate, unit, twidth,
685 idx0/rate, trace_style)
686 plot_pulse_eods(axr, data, rate,
687 zoom_window, twidth, eod_props,
688 idx0/rate, colors=pulse_colors,
689 markers=pulse_markers, frameon=True,
690 loc='upper right')
691 if axr.get_legend() is not None:
692 axr.get_legend().get_frame().set_color('white')
693 axr.set_title('Recording', fontsize=14, y=1.05)
694 axr.format_coord = recordingzoom_format_coord
696 # plot psd
697 wave_colors, wave_markers = colors_markers()
698 if axp is not None:
699 axes_style(axp)
700 if power_thresh is not None:
701 axp.plot(power_thresh[:,0], decibel(power_thresh[:,1]), '#CCCCCC', lw=1)
702 if len(wave_eodfs) > 0:
703 kwargs = {}
704 if len(wave_eodfs) > 1:
705 title = '%d EOD frequencies' % len(wave_eodfs)
706 kwargs = {'title': title if len(wave_eodfs) > 2 else None }
707 if legend_inside:
708 kwargs.update({'bbox_to_anchor': (1.05, 1.1),
709 'loc': 'upper right', 'legend_rows': 10,
710 'frameon': True})
711 else:
712 kwargs.update({'bbox_to_anchor': (1.02, 1.1),
713 'loc': 'upper left', 'legend_rows': 14,
714 'labelspacing': 0.6, 'frameon': False})
715 plot_harmonic_groups(axp, wave_eodfs, wave_indices, max_groups=0,
716 skip_bad=skip_bad,
717 sort_by_freq=True, label_power=label_power,
718 colors=wave_colors, markers=wave_markers,
719 **kwargs)
720 if legend_inside:
721 axp.get_legend().get_frame().set_color('white')
722 if psd_data is not None and len(psd_data) > 0:
723 plot_decibel_psd(axp, psd_data[:,0], psd_data[:,1],
724 log_freq=log_freq, min_freq=min_freq,
725 max_freq=max_freq, ymarg=5.0,
726 sstyle=spectrum_style)
727 axp.yaxis.set_major_locator(ticker.MaxNLocator(6))
728 if len(wave_eodfs) == 1:
729 axp.get_legend().set_visible(False)
730 label = '%6.1f Hz' % wave_eodfs[0][0, 0]
731 axp.set_title('Powerspectrum: %s' % label, y=1.05, fontsize=14)
732 else:
733 axp.set_title('Powerspectrum', y=1.05, fontsize=14)
734 axp.format_coord = psd_format_coord
736 # get fish labels from legends:
737 if axp is not None:
738 w, _ = axp.get_legend_handles_labels()
739 eodf_labels = [wi.get_label().split()[0] for wi in w]
740 legend_wave_eodfs = np.array([float(f) if f[0] != '(' else np.nan for f in eodf_labels])
741 if axr is not None:
742 p, _ = axr.get_legend_handles_labels()
743 eodf_labels = [pi.get_label().split()[0] for pi in p]
744 legend_pulse_eodfs = np.array([float(f) if f[0] != '(' else np.nan for f in eodf_labels])
746 # layout:
747 sheight = 1.4/height
748 sstep = 1.6/height
749 max_plots = len(indices)
750 if not all_eods:
751 if large_plots:
752 max_plots = 1 if spec_plots else 2
753 else:
754 max_plots = 4
755 if large_plots:
756 pstep = pheight + 1.0/height
757 ty = 1.08
758 my = 1.10
759 ny = 6
760 else:
761 posy -= 0.2/height
762 pheight = 1.3/height
763 pstep = 1.9/height
764 ty = 1.10
765 my = 1.16
766 ny = 4
767 posy -= pstep
769 # sort indices by p-p amplitude:
770 pp_ampls = [eod_props[idx]['p-p-amplitude'] for idx in indices]
771 pp_indices = np.argsort(pp_ampls)[::-1]
773 # plot EOD waveform and spectra:
774 for k, idx in enumerate(indices[pp_indices]):
775 if k >= max_plots:
776 break
777 # plot EOD waveform:
778 mean_eod = mean_eods[idx]
779 props = eod_props[idx]
780 peaks = peak_data[idx]
781 lx = leftx if spec_plots or k%2 == 0 else midx
782 ax = fig.add_axes([lx, posy, halfwidth, pheight])
783 axes_style(ax)
784 ax.yaxis.set_major_locator(ticker.MaxNLocator(ny))
785 if len(indices) > 1:
786 ax.text(0.3, ty, '{EODf:.1f} Hz {type} fish'.format(**props),
787 transform=ax.transAxes, fontsize=14, zorder=20)
788 mx = 0.25
789 else:
790 ax.text(-0.1, ty, '{EODf:.1f} Hz {type} fish'.format(**props),
791 transform=ax.transAxes, fontsize=14, zorder=20)
792 ax.text(0.5, ty, 'Averaged EOD', ha='center',
793 transform=ax.transAxes, fontsize=14, zorder=20)
794 mx = -0.14
795 eodf = props['EODf']
796 if props['type'] == 'wave':
797 if axp is not None:
798 wk = np.nanargmin(np.abs(legend_wave_eodfs - eodf))
799 ma = ml.Line2D([mx], [my], color=w[wk].get_color(), marker=w[wk].get_marker(),
800 markersize=w[wk].get_markersize(), mec='none', clip_on=False,
801 label=w[wk].get_label(), transform=ax.transAxes)
802 ax.add_line(ma)
803 else:
804 if axr is not None and len(legend_pulse_eodfs) > 0:
805 pk = np.argmin(np.abs(legend_pulse_eodfs - eodf))
806 ma = ml.Line2D([mx], [my], color=p[pk].get_color(), marker=p[pk].get_marker(),
807 markersize=p[pk].get_markersize(), mec='none', clip_on=False,
808 label=p[pk].get_label(), transform=ax.transAxes)
809 ax.add_line(ma)
810 plot_eod_waveform(ax, mean_eod, props, peaks, unit, **eod_styles)
811 if props['type'] == 'pulse' and 'times' in props:
812 plot_eod_snippets(ax, data, rate, mean_eod[0,0], mean_eod[-1,0],
813 props['times'], n_snippets, props['flipped'],
814 snippet_style)
815 if not large_plots and k < max_plots-2:
816 ax.set_xlabel('')
817 ax.format_coord = meaneod_format_coord
819 # plot spectra:
820 if spec_plots:
821 spec = spec_data[idx]
822 if props['type'] == 'pulse':
823 ax = fig.add_axes([midx, posy, halfwidth, pheight])
824 axes_style(ax)
825 plot_pulse_spectrum(ax, spec, props, **pulse_spec_styles)
826 ax.set_title('Single pulse spectrum', fontsize=14, y=1.05)
827 ax.format_coord = pulsepsd_format_coord
828 else:
829 axa = fig.add_axes([midx, posy+sstep, halfwidth, sheight])
830 axes_style(axa)
831 axp = fig.add_axes([midx, posy, halfwidth, sheight])
832 axes_style(axp)
833 plot_wave_spectrum(axa, axp, spec, props, unit,
834 **wave_spec_styles)
835 axa.set_title('Amplitude and phase spectrum', fontsize=14, y=1.05)
836 axa.set_xticklabels([])
837 axa.yaxis.set_major_locator(ticker.MaxNLocator(4))
838 axa.format_coord = ampl_format_coord
839 axp.format_coord = phase_format_coord
841 if spec_plots or k%2 == 1:
842 posy -= pstep
844 # whole trace:
845 ax = fig.add_axes([leftx, 0.6/height, fullwidth, 0.9/height])
846 axes_style(ax)
847 if raw_data is not None and len(raw_data) > 0:
848 plot_data_window(ax, raw_data, rate, unit, idx0, idx1, clipped,
849 **data_styles)
850 ax.format_coord = recording_format_coord
852 return fig
855def plot_eod_subplots(base_name, subplots, raw_data, rate, idx0, idx1,
856 clipped, psd_data, wave_eodfs, wave_indices, mean_eods,
857 eod_props, peak_data, spec_data, unit, zoom_window,
858 n_snippets=10, power_thresh=None, label_power=True,
859 skip_bad=True, log_freq=False,
860 min_freq=0.0, max_freq=3000.0, save=True):
861 """Plot time traces and spectra into separate windows or files.
863 Parameters
864 ----------
865 base_name: string
866 Basename of audio_file.
867 subplots: string
868 Specifies which subplots to plot:
869 r) recording with best window, t) data trace with detected pulse fish,
870 p) power spectrum with detected wave fish, w/W) mean EOD waveform,
871 s/S) EOD spectrum, e/E) EOD waveform and spectra. With capital letters
872 all fish are saved into a single pdf file, with small letters each fish
873 is saved into a separate file.
874 raw_data: array
875 Dataset.
876 rate: float
877 Sampling rate of the dataset.
878 idx0: float
879 Index of the beginning of the analysis window in the dataset.
880 idx1: float
881 Index of the end of the analysis window in the dataset.
882 clipped: float
883 Fraction of clipped amplitudes.
884 psd_data: 2D array
885 Power spectrum (frequencies and power) of the analysed data.
886 wave_eodfs: array
887 Frequency and power of fundamental frequency/harmonics of several fish.
888 wave_indices: array of int
889 Indices of wave fish mapping from wave_eodfs to eod_props.
890 If negative, then that EOD frequency has no waveform described in eod_props.
891 mean_eods: list of 2-D arrays with time, mean and std.
892 Mean trace for the mean EOD plot.
893 eod_props: list of dict
894 Properties for each waveform in mean_eods.
895 peak_data: list of 2_D arrays
896 For each pulsefish a list of peak properties
897 (index, time, and amplitude).
898 spec_data: list of 2_D arrays
899 For each pulsefish a power spectrum of the single pulse and for
900 each wavefish the relative amplitudes and phases of the harmonics.
901 unit: string
902 Unit of the trace and the mean EOD.
903 n_snippets: int
904 Number of EOD waveform snippets to be plotted. If zero do not plot any.
905 power_thresh: 2 D array or None
906 Frequency (first column) and power (second column) of threshold
907 derived from single pulse spectra to discard false wave fish.
908 label_power: boolean
909 If `True` put the power in decibel in addition to the frequency
910 into the legend.
911 skip_bad: bool
912 Skip harmonic groups without index (entry in indices is negative).
913 log_freq: boolean
914 Logarithmic (True) or linear (False) frequency axis of power spectrum of recording.
915 min_freq: float
916 Limits of frequency axis of power spectrum of recording
917 are set to `(min_freq, max_freq)` if `max_freq` is greater than zero
918 max_freq: float
919 Limits of frequency axis of power spectrum of recording
920 are set to `(min_freq, max_freq)` and limits of power axis are computed
921 from powers below max_freq if `max_freq` is greater than zero
922 save: bool
923 If True save plots to files instead of showing them.
924 """
925 plot_style()
926 if 'r' in subplots:
927 fig, ax = plt.subplots(figsize=(10, 2))
928 fig.subplots_adjust(left=0.07, right=0.99, bottom=0.22, top=0.95)
929 plot_data_window(ax, raw_data, rate, unit, idx0, idx1, clipped,
930 **data_styles)
931 ax.yaxis.set_major_locator(ticker.MaxNLocator(5))
932 axes_style(ax)
933 if save:
934 fig.savefig(base_name + '-recording.pdf')
935 if 't' in subplots:
936 fig, ax = plt.subplots(figsize=(10, 6))
937 twidth = 0.1
938 if len(eod_props) > 0:
939 if eod_props[0]['type'] == 'wave':
940 twidth = 5.0/eod_props[0]['EODf']
941 else:
942 if len(wave_eodfs) > 0:
943 twidth = 3.0/eod_props[0]['EODf']
944 else:
945 twidth = 10.0/eod_props[0]['EODf']
946 twidth = (1+twidth//0.005)*0.005
947 pulse_colors, pulse_markers = colors_markers()
948 pulse_colors = pulse_colors[3:]
949 pulse_markers = pulse_markers[3:]
950 plot_eod_recording(ax, raw_data[idx0:idx1], rate, unit,
951 twidth, idx0/rate, trace_style)
952 plot_pulse_eods(ax, raw_data[idx0:idx1], rate, zoom_window,
953 twidth, eod_props, idx0/rate,
954 colors=pulse_colors, markers=pulse_markers,
955 frameon=True, loc='upper right')
956 if ax.get_legend() is not None:
957 ax.get_legend().get_frame().set_color('white')
958 axes_style(ax)
959 if save:
960 fig.savefig(base_name + '-trace.pdf')
961 if 'p' in subplots:
962 fig, ax = plt.subplots(figsize=(10, 5))
963 fig.subplots_adjust(left=0.08, right=0.975, bottom=0.11, top=0.9)
964 axes_style(ax)
965 if power_thresh is not None:
966 ax.plot(power_thresh[:,0], decibel(power_thresh[:,1]), '#CCCCCC', lw=1)
967 if len(wave_eodfs) > 0:
968 kwargs = {}
969 if len(wave_eodfs) > 1:
970 title = '%d EOD frequencies' % len(wave_eodfs)
971 kwargs = {'title': title if len(wave_eodfs) > 2 else None }
972 if len(wave_eodfs) > 2:
973 fig.subplots_adjust(left=0.08, right=0.78, bottom=0.11, top=0.9)
974 kwargs.update({'bbox_to_anchor': (1.01, 1.1),
975 'loc': 'upper left', 'legend_rows': 14,
976 'labelspacing': 0.6})
977 else:
978 kwargs.update({'bbox_to_anchor': (1.05, 1.1),
979 'loc': 'upper right', 'legend_rows': 10})
980 wave_colors, wave_markers = colors_markers()
981 plot_harmonic_groups(ax, wave_eodfs, wave_indices, max_groups=0,
982 skip_bad=skip_bad,
983 sort_by_freq=True, label_power=label_power,
984 colors=wave_colors, markers=wave_markers,
985 frameon=False, **kwargs)
986 plot_decibel_psd(ax, psd_data[:,0], psd_data[:,1], log_freq=log_freq,
987 min_freq=min_freq, max_freq=max_freq, ymarg=5.0,
988 sstyle=spectrum_style)
989 ax.yaxis.set_major_locator(ticker.MaxNLocator(6))
990 if len(wave_eodfs) == 1:
991 ax.get_legend().set_visible(False)
992 label = '%6.1f Hz' % wave_eodfs[0][0, 0]
993 ax.set_title('Powerspectrum: %s' % label, y=1.05)
994 else:
995 ax.set_title('Powerspectrum', y=1.05)
996 if save:
997 fig.savefig(base_name + '-psd.pdf')
998 if 'w' in subplots or 'W' in subplots:
999 mpdf = None
1000 if 'W' in subplots:
1001 mpdf = PdfPages(base_name + '-waveforms.pdf')
1002 for meod, props, peaks in zip(mean_eods, eod_props, peak_data):
1003 if meod is None:
1004 continue
1005 fig, ax = plt.subplots(figsize=(5, 3))
1006 fig.subplots_adjust(left=0.18, right=0.98, bottom=0.15, top=0.9)
1007 if not props is None:
1008 ax.set_title('{index:d}: {EODf:.1f} Hz {type} fish'.format(**props))
1009 plot_eod_waveform(ax, meod, props, peaks, unit, **eod_styles)
1010 data = raw_data[idx0:idx1] if idx1 > idx0 else raw_data
1011 if not props is None and props['type'] == 'pulse' and \
1012 'times' in props:
1013 plot_eod_snippets(ax, data, rate, meod[0,0],
1014 meod[-1,0], props['times'],
1015 n_snippets, False, snippet_style)
1016 ax.yaxis.set_major_locator(ticker.MaxNLocator(6))
1017 axes_style(ax)
1018 if mpdf is None:
1019 if save:
1020 fig.savefig(base_name + '-waveform-%d.pdf' % props['index'])
1021 else:
1022 mpdf.savefig(fig)
1023 if mpdf is not None:
1024 mpdf.close()
1025 if 's' in subplots or 'S' in subplots:
1026 mpdf = None
1027 if 'S' in subplots:
1028 mpdf = PdfPages(base_name + '-spectrum.pdf')
1029 for props, peaks, spec in zip(eod_props, peak_data, spec_data):
1030 if spec is None:
1031 continue
1032 if props is not None and props['type'] == 'pulse':
1033 fig, ax = plt.subplots(figsize=(5, 3.5))
1034 fig.subplots_adjust(left=0.15, right=0.967, bottom=0.16, top=0.88)
1035 axes_style(ax)
1036 ax.set_title('{index:d}: {EODf:.1f} Hz {type} fish'.format(**props), y=1.07)
1037 plot_pulse_spectrum(ax, spec, props, **pulse_spec_styles)
1038 else:
1039 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(5, 3.5))
1040 fig.subplots_adjust(left=0.15, right=0.97, bottom=0.16, top=0.88, hspace=0.4)
1041 axes_style(ax1)
1042 axes_style(ax2)
1043 if not props is None:
1044 ax1.set_title('{index:d}: {EODf:.1f} Hz {type} fish'.format(**props), y=1.15)
1045 plot_wave_spectrum(ax1, ax2, spec, props, unit,
1046 **wave_spec_styles)
1047 ax1.set_xticklabels([])
1048 ax1.yaxis.set_major_locator(ticker.MaxNLocator(4))
1049 if mpdf is None:
1050 if save:
1051 fig.savefig(base_name + '-spectrum-%d.pdf' % props['index'])
1052 else:
1053 mpdf.savefig(fig)
1054 if mpdf is not None:
1055 mpdf.close()
1056 if 'e' in subplots or 'E' in subplots:
1057 mpdf = None
1058 if 'E' in subplots:
1059 mpdf = PdfPages(base_name + '-eods.pdf')
1060 for meod, props, peaks, spec in zip(mean_eods, eod_props, peak_data, spec_data):
1061 if meod is None or spec is None:
1062 continue
1063 fig = plt.figure(figsize=(10, 3.5))
1064 gs = gridspec.GridSpec(nrows=2, ncols=2, left=0.09, right=0.98,
1065 bottom=0.16, top=0.88, wspace=0.4, hspace=0.4)
1066 ax1 = fig.add_subplot(gs[:,0])
1067 if not props is None:
1068 ax1.set_title('{index:d}: {EODf:.1f} Hz {type} fish'.format(**props), y=1.07)
1069 plot_eod_waveform(ax1, meod, props, peaks, unit, **eod_styles)
1070 data = raw_data[idx0:idx1] if idx1 > idx0 else raw_data
1071 if not props is None and props['type'] == 'pulse' and 'times' in props:
1072 plot_eod_snippets(ax1, data, rate, meod[0,0],
1073 meod[-1,0], props['times'],
1074 n_snippets, False, snippet_style)
1075 ax1.yaxis.set_major_locator(ticker.MaxNLocator(6))
1076 axes_style(ax1)
1077 if not props is None and props['type'] == 'pulse':
1078 ax2 = fig.add_subplot(gs[:,1])
1079 axes_style(ax2)
1080 plot_pulse_spectrum(ax2, spec, props, **pulse_spec_styles)
1081 ax2.set_title('Single pulse spectrum', y=1.07)
1082 else:
1083 ax2 = fig.add_subplot(gs[0,1])
1084 ax3 = fig.add_subplot(gs[1,1])
1085 axes_style(ax2)
1086 axes_style(ax3)
1087 plot_wave_spectrum(ax2, ax3, spec, props, unit,
1088 **wave_spec_styles)
1089 ax2.set_title('Amplitude and phase spectrum', y=1.15)
1090 ax2.set_xticklabels([])
1091 ax2.yaxis.set_major_locator(ticker.MaxNLocator(4))
1092 if mpdf is None:
1093 if save:
1094 fig.savefig(base_name + '-eod-%d.pdf' % props['index'])
1095 else:
1096 mpdf.savefig(fig)
1097 if mpdf is not None:
1098 mpdf.close()
1099 if not save:
1100 plt.show()
1101 plt.close('all')
1104def thunderfish_plot(files, data_path=None, load_kwargs={},
1105 all_eods=False, spec_plots='auto', skip_bad=True,
1106 save_plot=False, multi_pdf=None,
1107 save_subplots='', log_freq=False, min_freq=0.0,
1108 max_freq=3000.0, output_folder='.',
1109 keep_path=False, verbose=0):
1110 """Generate plots from saved analysis results.
1112 Parameters
1113 ----------
1114 files: list of str
1115 Analysis files from a single recording.
1116 data_path: str
1117 Path where to find the raw data.
1118 load_kwargs: dict
1119 Key-word arguments for the `load_data()` function.
1120 all_eods: bool
1121 If True, plot all EOD waveforms.
1122 spec_plots: bool or 'auto'
1123 Plot amplitude spectra of EOD waveforms.
1124 If 'auto', plot them if there is a singel waveform only.
1125 skip_bad: bool
1126 Skip harmonic groups without index in the spectrum plot.
1127 save_plot: bool
1128 If True, save plots as pdf file.
1129 multi_pdf: matplotlib.PdfPages or None
1130 PdfPages instance in which to save plots.
1131 save_subplots: string
1132 If not empty, specifies subplots to be saved as separate pdf
1133 files: r) recording with best window, t) data trace with
1134 detected pulse fish, p) power spectrum with detected wave
1135 fish, w/W) mean EOD waveform, s/S) EOD spectrum, e/E) EOD
1136 waveform and spectra. Capital letters produce a single
1137 multipage pdf containing plots of all detected fish.
1138 log_freq: boolean
1139 Logarithmic (True) or linear (False) frequency axis of
1140 power spectrum of recording.
1141 min_freq: float
1142 Limits of frequency axis of power spectrum of recording
1143 are set to `(min_freq, max_freq)`, if `max_freq` is greater than zero
1144 max_freq: float
1145 Limits of frequency axis of power spectrum of recording
1146 are set to `(min_freq, max_freq)` and limits of power axis are computed
1147 from powers below max_freq, if `max_freq` is greater than zero
1148 output_folder: string
1149 Folder where to save results.
1150 keep_path: bool
1151 Add relative path of data files to output path.
1152 verbose: int
1153 Verbosity level (for debugging).
1154 """
1155 if len(save_subplots) == 0:
1156 save_subplots = 'rtpwsed' # plot everything
1157 # load analysis results:
1158 mean_eods, wave_eodfs, wave_indices, eod_props, spec_data, \
1159 peak_data, base_name, channel, unit = load_analysis(files)
1160 if len(mean_eods) == 0 or all(me is None for me in mean_eods):
1161 save_subplots = save_subplots.replace('w', '')
1162 save_subplots = save_subplots.replace('W', '')
1163 save_subplots = save_subplots.replace('e', '')
1164 save_subplots = save_subplots.replace('E', '')
1165 save_subplots = save_subplots.replace('d', '')
1166 if len(spec_data) == 0 or all(sd is None for sd in spec_data):
1167 save_subplots = save_subplots.replace('s', '')
1168 save_subplots = save_subplots.replace('S', '')
1169 save_subplots = save_subplots.replace('e', '')
1170 save_subplots = save_subplots.replace('E', '')
1171 save_subplots = save_subplots.replace('d', '')
1172 clipped = 0.0
1173 if len(eod_props) > 0 and not eod_props[0] is None and \
1174 'winclipped' in eod_props[0]:
1175 clipped = eod_props[0]['winclipped']
1176 zoom_window = [1.2, 1.3]
1177 # load recording:
1178 psd_data = None
1179 if base_name:
1180 name = os.path.basename(base_name) if data_path and data_path != '.' else base_name
1181 data_path = os.path.join(data_path, name)
1182 data, rate, idx0, idx1, data_path = \
1183 load_recording(data_path, channel, load_kwargs,
1184 eod_props, verbose-1)
1185 if data is None:
1186 save_subplots = save_subplots.replace('r', '')
1187 save_subplots = save_subplots.replace('t', '')
1188 save_subplots = save_subplots.replace('d', '')
1189 if verbose > 0:
1190 print('loaded', data_path)
1191 if len(eod_props) > 0 and not eod_props[0] is None and \
1192 'dfreq' in eod_props[0] and data is not None and len(data) > 0:
1193 psd_data = multi_psd(data[idx0:idx1], rate,
1194 1.1*eod_props[0]['dfreq'])[0]
1195 if psd_data is not None and len(psd_data) > 0:
1196 for idx, fish in zip(wave_indices, wave_eodfs):
1197 if idx < 0:
1198 for k in range(len(fish)):
1199 fish[k,1] = psd_data[np.argmin(np.abs(psd_data[:,0] - fish[k,0])),1]
1200 if psd_data is None:
1201 save_subplots = save_subplots.replace('p', '')
1202 save_subplots = save_subplots.replace('d', '')
1203 # file name for output files:
1204 fn = base_name if keep_path else os.path.basename(base_name)
1205 output_basename = os.path.join(output_folder, fn)
1206 if channel >= 0:
1207 output_basename += f'-c{channel}'
1208 # make directory if necessary:
1209 if keep_path:
1210 outpath = os.path.dirname(output_basename)
1211 if not os.path.exists(outpath):
1212 if verbose > 0:
1213 print('mkdir %s' % outpath)
1214 os.makedirs(outpath)
1215 # plot:
1216 if len(save_subplots) == 0 or 'd' in save_subplots:
1217 fig = plot_eods(os.path.basename(base_name), None, data, rate,
1218 channel, idx0, idx1, clipped, psd_data,
1219 wave_eodfs, wave_indices, mean_eods,
1220 eod_props, peak_data, spec_data, None, unit,
1221 zoom_window, 10, None, True, all_eods,
1222 spec_plots, skip_bad, log_freq, min_freq,
1223 max_freq, interactive=not save_plot,
1224 verbose=verbose-1)
1225 if save_plot:
1226 if multi_pdf is not None:
1227 multi_pdf.savefig(fig)
1228 else:
1229 fig.savefig(output_basename + '.pdf')
1230 else:
1231 fig.canvas.manager.set_window_title('thunderfish')
1232 plt.show()
1233 plt.close()
1234 save_subplots = save_subplots.replace('d', '')
1235 if len(save_subplots) > 0:
1236 plot_eod_subplots(output_basename, save_subplots, data, rate,
1237 idx0, idx1, clipped, psd_data, wave_eodfs,
1238 wave_indices, mean_eods, eod_props,
1239 peak_data, spec_data, unit, zoom_window, 10,
1240 None, True, skip_bad, log_freq, min_freq,
1241 max_freq, save_plot)
1242 return None
1245def thunderfish(filename, load_kwargs, cfg, channel=0,
1246 time=None, time_file=False,
1247 mode='wp', log_freq=False, min_freq=0.0, max_freq=3000,
1248 save_data=False, zip_file=False,
1249 all_eods=False, spec_plots='auto', skip_bad=True,
1250 save_plot=False, multi_pdf=None, save_subplots='',
1251 output_folder='.', keep_path=False,
1252 verbose=0, plot_level=0):
1253 """Automatically detect and analyze all EOD waveforms in a short recording.
1255 Parameters
1256 ----------
1257 filename: string
1258 Path of the data file to be analyzed.
1259 load_kwargs: dict
1260 Key-word arguments for the `load_data()` function.
1261 cfg: dict
1262 channel: int
1263 Channel to be analyzed.
1264 time: string, float, or None
1265 Start time of analysis window: "beginning", "center", "end",
1266 "best", or time in seconds (as float or string). If not None
1267 overwrites "windowPosition" in cofiguration file.
1268 time_file: bool
1269 If `True` add time of analysis window to output file names.
1270 mode: 'w', 'p', 'P', 'wp', or 'wP'
1271 Analyze wavefish ('w'), all pulse fish ('p'), or largest pulse
1272 fish only ('P').
1273 log_freq: boolean
1274 Logarithmic (True) or linear (False) frequency axis of
1275 power spectrum of recording.
1276 min_freq: float
1277 Limits of frequency axis of power spectrum of recording
1278 are set to `(min_freq, max_freq)`, if `max_freq` is greater than zero
1279 max_freq: float
1280 Limits of frequency axis of power spectrum of recording
1281 are set to `(min_freq, max_freq)` and limits of power axis are computed
1282 from powers below max_freq, if `max_freq` is greater than zero
1283 save_data: bool
1284 If True save analysis results in files. If False, just plot the data.
1285 zip_data: bool
1286 If True, store all analysis results in a single zip file.
1287 all_eods: bool
1288 If True, plot all EOD waveforms.
1289 spec_plots: bool or 'auto'
1290 Plot amplitude spectra of EOD waveforms.
1291 If 'auto', plot them if there is a singel waveform only.
1292 skip_bad: bool
1293 Skip harmonic groups without index in the spectrum plot.
1294 save_plot: bool
1295 If True, save plots as pdf file.
1296 multi_pdf: matplotlib.PdfPages or None
1297 PdfPages instance in which to save plots.
1298 save_subplots: string
1299 If not empty, specifies subplots to be saved as separate pdf
1300 files: r) recording with best window, t) data trace with
1301 detected pulse fish, p) power spectrum with detected wave
1302 fish, w/W) mean EOD waveform, s/S) EOD spectrum, e/E) EOD
1303 waveform and spectra. Capital letters produce a single
1304 multipage pdf containing plots of all detected fish.
1305 output_folder: string
1306 Folder where to save results.
1307 keep_path: bool
1308 Add relative path of data files to output path.
1309 verbose: int
1310 Verbosity level (for debugging).
1311 plot_level: int
1312 Plot intermediate results.
1314 Returns
1315 -------
1316 msg: string or None
1317 In case of errors, an error message.
1318 """
1319 # check data file:
1320 if len(filename) == 0:
1321 return 'you need to specify a file containing some data'
1323 # file names:
1324 fn = filename if keep_path else os.path.basename(filename)
1325 outfilename = os.path.splitext(fn)[0]
1326 messagefilename = os.path.splitext(fn)[0] + '-message.wav'
1327 if not os.path.isfile(messagefilename):
1328 messagefilename = None
1330 # load data:
1331 try:
1332 all_data, rate, unit, ampl_max = load_data(filename,
1333 verbose=verbose,
1334 **load_kwargs)
1335 except IOError as e:
1336 return '%s: failed to open file: %s' % (filename, str(e))
1337 # select channel:
1338 channels = all_data.shape[1]
1339 chan_list = [channel]
1340 if channel < 0:
1341 chan_list = range(channels)
1342 elif channel >= channels:
1343 return '%s: invalid channel %d (%d channels)' % (filename, channel, channels)
1344 # process all channels:
1345 for chan in chan_list:
1346 raw_data = all_data[:,chan]
1347 if len(raw_data) <= 1:
1348 return '%s: empty data file' % filename
1349 if verbose >= 0 and len(chan_list) > 1:
1350 print(' channel %d' % chan)
1352 # analysis window:
1353 win_pos = cfg.value('windowPosition')
1354 if time is not None:
1355 win_pos = time
1356 data, idx0, idx1, clipped, min_clip, max_clip = \
1357 analysis_window(raw_data, rate, ampl_max, win_pos, cfg,
1358 plot_level>0)
1359 found_bestwindow = idx1 > 0
1360 if not found_bestwindow:
1361 return '%s: not enough data for requested window length. You may want to adjust the windowSize parameter in the configuration file.' % filename
1363 # detect EODs in the data:
1364 psd_data, wave_eodfs, wave_indices, eod_props, \
1365 mean_eods, spec_data, peak_data, power_thresh, skip_reason, zoom_window = \
1366 detect_eods(data, rate, min_clip, max_clip, filename, mode,
1367 verbose, plot_level, cfg)
1368 if not found_bestwindow:
1369 wave_eodfs = []
1370 wave_indices = []
1371 eod_props = []
1372 mean_eods = []
1374 # add analysis window to EOD properties:
1375 for props in eod_props:
1376 props['twin'] = idx0/rate
1377 props['window'] = (idx1 - idx0)/rate
1378 props['winclipped'] = clipped
1380 # warning message in case no fish has been found:
1381 if found_bestwindow and not eod_props :
1382 msg = ', '.join(skip_reason)
1383 if msg:
1384 print(filename + ': no fish found: %s' % msg)
1385 else:
1386 print(filename + ': no fish found.')
1388 # file name for output files:
1389 output_basename = os.path.join(output_folder, outfilename)
1390 if channels > 1:
1391 if channels > 100:
1392 output_basename += '-c%03d' % chan
1393 elif channels > 10:
1394 output_basename += '-c%02d' % chan
1395 else:
1396 output_basename += '-c%d' % chan
1397 if time_file:
1398 output_basename += '-t%.0fs' % (idx0/rate)
1399 # make directory if necessary:
1400 if keep_path and found_bestwindow:
1401 outpath = os.path.dirname(output_basename)
1402 if not os.path.exists(outpath):
1403 if verbose > 0:
1404 print('mkdir %s' % outpath)
1405 os.makedirs(outpath)
1406 # save results to files:
1407 if save_data:
1408 remove_eod_files(output_basename, verbose, cfg)
1409 if found_bestwindow:
1410 save_analysis(output_basename, zip_file, eod_props,
1411 mean_eods, spec_data, peak_data,
1412 wave_eodfs, wave_indices, unit, verbose,
1413 **write_table_args(cfg))
1414 # summary plots:
1415 if save_plot or not save_data:
1416 n_snippets = 10
1417 if len(save_subplots) == 0 or 'd' in save_subplots:
1418 chl = chan if channels > 1 else None
1419 fig = plot_eods(outfilename, messagefilename,
1420 raw_data, rate, chl, idx0, idx1,
1421 clipped, psd_data[0], wave_eodfs,
1422 wave_indices, mean_eods, eod_props,
1423 peak_data, spec_data, None, unit,
1424 zoom_window, n_snippets, power_thresh,
1425 True, all_eods, spec_plots, skip_bad,
1426 log_freq, min_freq, max_freq,
1427 interactive=not save_plot,
1428 verbose=verbose)
1429 if save_plot:
1430 if multi_pdf is not None:
1431 multi_pdf.savefig(fig)
1432 else:
1433 fig.savefig(output_basename + '.pdf')
1434 else:
1435 fig.canvas.manager.set_window_title('thunderfish')
1436 plt.show()
1437 plt.close()
1438 save_subplots = save_subplots.replace('d', '')
1439 if len(save_subplots) > 0:
1440 plot_eod_subplots(output_basename, save_subplots,
1441 raw_data, rate, idx0, idx1, clipped,
1442 psd_data[0], wave_eodfs,
1443 wave_indices, mean_eods, eod_props,
1444 peak_data, spec_data, unit,
1445 zoom_window, n_snippets,
1446 power_thresh, True, skip_bad,
1447 log_freq, min_freq, max_freq,
1448 save_plot)
1449 return None
1452def run_thunderfish(file_args):
1453 """Helper function for mutlithreading Pool().map().
1454 """
1455 results = file_args[1][0]
1456 verbose = file_args[1][-1] if results else file_args[1][-2]+1
1457 if verbose > 1:
1458 print('='*70)
1459 try:
1460 if results:
1461 thunderfish_plot(file_args[0], *file_args[1][1:])
1462 else:
1463 if verbose > 0:
1464 print('analyze recording %s ...' % file_args[0])
1465 msg = thunderfish(file_args[0], *file_args[1][1:])
1466 if msg:
1467 print(msg)
1468 except (KeyboardInterrupt, SystemExit):
1469 print('\nthunderfish interrupted by user... exit now.')
1470 sys.exit(0)
1471 except:
1472 print(traceback.format_exc())
1475def main(cargs=None):
1476 # config file name:
1477 cfgfile = __package__ + '.cfg'
1479 # command line arguments:
1480 if cargs is None:
1481 cargs = sys.argv[1:]
1482 parser = argparse.ArgumentParser(add_help=False,
1483 description='Analyze EOD waveforms of weakly electric fish.',
1484 epilog='version %s by Benda-Lab (2015-%s)' % (__version__, __year__))
1485 parser.add_argument('-h', '--help', action='store_true',
1486 help='show this help message and exit')
1487 parser.add_argument('--version', action='version', version=__version__)
1488 parser.add_argument('-v', action='count', dest='verbose', default=0,
1489 help='verbosity level. Increase by specifying -v multiple times, or like -vvv')
1490 parser.add_argument('-V', action='count', dest='plot_level', default=0,
1491 help='level for debugging plots. Increase by specifying -V multiple times, or like -VVV')
1492 parser.add_argument('-c', dest='save_config', action='store_true',
1493 help='save configuration to file {0} after reading all configuration files'.format(cfgfile))
1494 parser.add_argument('--channel', default=0, type=int,
1495 help='channel to be analyzed (defaults to first channel, negative channel selects all channels)')
1496 parser.add_argument('-t', dest='time', default=None, type=str, metavar='TIME',
1497 help='start time of analysis window in recording: "beginning", "center", "end", "best", or time in seconds (overwrites "windowPosition" in cofiguration file)')
1498 parser.add_argument('-u', dest='unwrap', action='store_true',
1499 help='unwrap clipped files, toggles unwrap setting of config file.')
1500 parser.add_argument('-T', dest='time_file', action='store_true',
1501 help='add start time of analysis file to output file names')
1502 parser.add_argument('-m', dest='mode', default='wp', type=str,
1503 choices=['w', 'p', 'wp'],
1504 help='extract wave "w" and/or pulse "p" fish EODs')
1505 parser.add_argument('-a', dest='all_eods', action='store_true',
1506 help='show all EOD waveforms in the summary plot')
1507 parser.add_argument('-S', dest='spec_plots', action='store_true',
1508 help='plot spectra for all EOD waveforms in the summary plot')
1509 parser.add_argument('-b', dest='skip_bad', action='store_false',
1510 help='indicate bad EODs in legend of power spectrum')
1511 parser.add_argument('-l', dest='log_freq', type=float, metavar='MINFREQ',
1512 nargs='?', const=100.0, default=0.0,
1513 help='logarithmic frequency axis in power spectrum with optional minimum frequency (defaults to 100 Hz)')
1514 parser.add_argument('-p', dest='save_plot', action='store_true',
1515 help='save output plots as pdf files')
1516 parser.add_argument('-M', dest='multi_pdf', default='', type=str, metavar='PDFFILE',
1517 help='save all summary plots of all recordings in a multi page pdf file. Disables parallel jobs.')
1518 parser.add_argument('-P', dest='save_subplots', default='', type=str, metavar='rtpwsed',
1519 help='save subplots as separate pdf files: r) recording with analysis window, t) data trace with detected pulse fish, p) power spectrum with detected wave fish, w/W) mean EOD waveform, s/S) EOD spectrum, e/E) EOD waveform and spectra, d) the default summary plot. Capital letters produce a single multipage pdf containing plots of all detected fish')
1520 parser.add_argument('-d', dest='rawdata_path', default='.', type=str, metavar='PATH',
1521 help='path to raw EOD recordings needed for plotting based on analysis results')
1522 parser.add_argument('-j', dest='jobs', nargs='?', type=int, default=None, const=0,
1523 help='number of jobs run in parallel. Without argument use all CPU cores.')
1524 parser.add_argument('-s', dest='save_data', action='store_true',
1525 help='save analysis results to files')
1526 parser.add_argument('-z', dest='zip_file', action='store_true',
1527 help='save analysis results in a single zip file')
1528 parser.add_argument('-f', dest='format', default='auto', type=str,
1529 choices=TableData.formats + ['py'],
1530 help='file format used for saving analysis results, defaults to the format specified in the configuration file or "csv"')
1531 parser.add_argument('-o', dest='outpath', default='.', type=str,
1532 help='path where to store results and figures (defaults to current working directory)')
1533 parser.add_argument('-k', dest='keep_path', action='store_true',
1534 help='keep path of input file when saving analysis files, i.e. append path of input file to OUTPATH')
1535 parser.add_argument('-i', dest='load_kwargs', default=[],
1536 action='append', metavar='KWARGS',
1537 help='key-word arguments for the data loader function')
1538 parser.add_argument('file', nargs='*', default='', type=str,
1539 help='name of a file with time series data of an EOD recording, may include wildcards')
1540 args = parser.parse_args(cargs)
1542 # help:
1543 if args.help:
1544 parser.print_help()
1545 print('')
1546 print('examples:')
1547 print('- analyze the single file data.wav interactively:')
1548 print(' > thunderfish data.wav')
1549 print('- extract wavefish only:')
1550 print(' > thunderfish -m w data.wav')
1551 print('- automatically analyze all wav files in the current working directory and save analysis results and plot to files:')
1552 print(' > thunderfish -s -p *.wav')
1553 print('- analyze all wav files in the river1/ directory, use all CPUs, and write files directly to "results/":')
1554 print(' > thunderfish -j -s -p -o results/ river1/*.wav')
1555 print('- analyze all wav files in the river1/ directory and write files to "results/river1/":')
1556 print(' > thunderfish -s -p -o results/ -k river1/*.wav')
1557 print('- write configuration file:')
1558 print(' > thunderfish -c')
1559 parser.exit()
1561 # set verbosity level from command line:
1562 verbose = args.verbose
1563 plot_level = args.plot_level
1564 if verbose < plot_level:
1565 verbose = plot_level
1567 # interactive plot:
1568 plt.rcParams['keymap.quit'] = 'ctrl+w, alt+q, q'
1569 plt.ioff()
1571 # expand wildcard patterns:
1572 files = []
1573 if os.name == 'nt':
1574 for fn in args.file:
1575 files.extend(glob.glob(fn))
1576 else:
1577 files = [f for f in args.file if '-message' not in f]
1579 # save configuration:
1580 if args.save_config:
1581 file_name = files[0] if len(files) else ''
1582 cfg = configuration()
1583 cfg.load_files(cfgfile, file_name, 4, verbose)
1584 save_configuration(cfg, cfgfile)
1585 exit()
1586 elif len(files) == 0:
1587 parser.error('you need to specify at least one file for the analysis')
1589 # configure:
1590 cfg = configuration()
1591 cfg.load_files(cfgfile, files[0], 4, verbose)
1592 if args.format != 'auto':
1593 cfg.set('fileFormat', args.format)
1594 if args.unwrap:
1595 cfg.set('unwrapData', not cfg.value('unwrapData'))
1597 # plot parameter:
1598 spec_plots = 'auto'
1599 if args.spec_plots:
1600 spec_plots = True
1602 # multi-page pdfs:
1603 multi_pdf = None
1604 if len(args.multi_pdf) > 0:
1605 args.save_plot = True
1606 args.jobs = None # PdfPages does not work yet with mutliprocessing
1607 ext = os.path.splitext(args.multi_pdf)[1]
1608 if ext != os.extsep + 'pdf':
1609 args.multi_pdf += os.extsep + 'pdf'
1610 multi_pdf = PdfPages(args.multi_pdf)
1612 # create output folder:
1613 if args.save_data or args.save_plot:
1614 if not os.path.exists(args.outpath):
1615 if verbose > 1:
1616 print('mkdir %s' % args.outpath)
1617 os.makedirs(args.outpath)
1619 # kwargs for data loader:
1620 load_kwargs = {}
1621 for s in args.load_kwargs:
1622 for kw in s.split(','):
1623 kws = kw.split(':')
1624 if len(kws) == 2:
1625 load_kwargs[kws[0].strip()] = kws[1].strip()
1627 # frequency limits for power spectrum:
1628 min_freq = 0.0
1629 max_freq = 3000.0
1630 log_freq = args.log_freq
1631 if log_freq > 0.0:
1632 min_freq = log_freq
1633 max_freq = min_freq*20
1634 if max_freq < 2000:
1635 max_freq = 2000
1636 log_freq = True
1637 else:
1638 log_freq = False
1640 # check if all input files are results:
1641 exts = TableData.ext_formats.values()
1642 results = True
1643 # check and group by recording:
1644 result_files = []
1645 for f in sorted(files):
1646 _, base_name, _, _, ftype, _, ext = parse_filename(f)
1647 if ext == 'zip' or (ext in exts and ftype in file_types):
1648 if len(result_files) == 0 or \
1649 not result_files[-1][-1].startswith(base_name):
1650 result_files.append([f])
1651 else:
1652 result_files[-1].append(f)
1653 else:
1654 results = False
1655 break
1656 if results:
1657 files = result_files
1659 # adjust verbosity:
1660 v = verbose
1661 if len(files) > 1:
1662 v += 1
1664 # run on pool:
1665 pool_args = (results, load_kwargs, cfg, args.channel, args.time,
1666 args.time_file, args.mode, log_freq, min_freq,
1667 max_freq, args.save_data, args.zip_file,
1668 args.all_eods, spec_plots, args.skip_bad,
1669 args.save_plot, multi_pdf, args.save_subplots,
1670 args.outpath, args.keep_path, v-1, plot_level)
1671 if results:
1672 pool_args = (results, args.rawdata_path, load_kwargs,
1673 args.all_eods, spec_plots, args.skip_bad,
1674 args.save_plot, multi_pdf, args.save_subplots,
1675 log_freq, min_freq, max_freq, args.outpath,
1676 args.keep_path, v)
1677 if args.jobs is not None and (args.save_data or args.save_plot) and len(files) > 1:
1678 cpus = cpu_count() if args.jobs == 0 else args.jobs
1679 if verbose > 1:
1680 print('run on %d cpus' % cpus)
1681 p = Pool(cpus)
1682 p.map(run_thunderfish, zip(files, [pool_args]*len(files)))
1683 else:
1684 list(map(run_thunderfish, zip(files, [pool_args]*len(files))))
1685 if multi_pdf is not None:
1686 multi_pdf.close()
1689if __name__ == '__main__':
1690 freeze_support() # needed by multiprocessing for some weired windows stuff
1691 main()