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