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