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