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