Coverage for src/thunderfish/thunderfish.py: 0%
905 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 16:21 +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, 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, samplerate, 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 samplerate: 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, samplerate, **multi_psd_args(cfg))
188 dfreq = np.mean(np.diff(psd_data[0][:,0]))
189 nfft = int(samplerate/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, samplerate, 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, samplerate, 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, samplerate, 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'] = samplerate
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*samplerate)
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, samplerate, **multi_psd_args(cfg))
286 pulse_power = pulse_psd[0][:,1]
287 pulse_power *= len(data)/samplerate/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 eod_times = np.arange(0.0, len(data)/samplerate, 1.0/fish[0,0])
318 mean_eod, eod_times = \
319 eod_waveform(data, samplerate, eod_times, win_fac=3.0, min_win=0.0,
320 min_sem=(k==0), **eod_waveform_args(cfg))
321 mean_eod, props, sdata, error_str = \
322 analyze_wave(mean_eod, fish, **analyze_wave_args(cfg))
323 if error_str:
324 print(name + ': ' + error_str)
325 clipped_frac = clipped_fraction(data, samplerate, eod_times,
326 mean_eod, min_clip, max_clip)
327 props['n'] = len(eod_times)
328 props['index'] = len(eod_props)
329 props['clipped'] = clipped_frac
330 props['samplerate'] = samplerate
331 props['nfft'] = nfft
332 props['dfreq'] = dfreq
333 # remove wave fish that are smaller than the largest pulse fish:
334 if props['p-p-amplitude'] < 0.01*max_pulse_amplitude:
335 rm_indices = power_indices[k:]
336 if verbose > 0:
337 print('skip %6.1fHz wave fish: power=%5.1fdB, p-p amplitude=%5.1fdB smaller than pulse fish=%5.1dB - 20dB' %
338 (props['EODf'], decibel(powers[idx]),
339 decibel(props['p-p-amplitude']), decibel(max_pulse_amplitude)))
340 for idx in rm_indices[1:]:
341 print('skip %6.1fHz wave fish: power=%5.1fdB even smaller' %
342 (wave_eodfs[idx][0,0], decibel(powers[idx])))
343 wave_eodfs = [eodfs for idx, eodfs in enumerate(wave_eodfs)
344 if idx not in rm_indices]
345 wave_indices = np.array([idcs for idx, idcs in enumerate(wave_indices)
346 if idx not in rm_indices], dtype=int)
347 break
348 # add good waveforms only:
349 remove, skips, msg = wave_quality(props, sdata[1:,3], **wave_quality_args(cfg))
350 if len(skips) == 0:
351 wave_indices[idx] = props['index']
352 eod_props.append(props)
353 mean_eods.append(mean_eod)
354 spec_data.append(sdata)
355 peak_data.append([])
356 if verbose > 0:
357 print('take %6.1fHz wave fish: %s' % (props['EODf'], msg))
358 else:
359 wave_indices[idx] = -2 if remove else -1
360 skip_reason += ['%.1fHz wave fish %s' % (props['EODf'], skips)]
361 if verbose > 0:
362 print('%-6s %6.1fHz wave fish: %s (%s)' %
363 ('remove' if remove else 'skip', props['EODf'], skips, msg))
364 wave_eodfs = [eodfs for idx, eodfs in zip(wave_indices, wave_eodfs) if idx > -2]
365 wave_indices = np.array([idx for idx in wave_indices if idx > -2], dtype=int)
366 return (psd_data, wave_eodfs, wave_indices, eod_props, mean_eods,
367 spec_data, peak_data, power_thresh, skip_reason, zoom_window)
370def remove_eod_files(output_basename, verbose, cfg):
371 """Remove all files from previous runs of thunderfish
372 """
373 ff = cfg.value('fileFormat')
374 if ff == 'py':
375 fext = 'py'
376 else:
377 fext = TableData.extensions[ff]
378 # remove all files from previous runs of thunderfish:
379 for fn in glob.glob('%s*.%s' % (output_basename, fext)):
380 os.remove(fn)
381 if verbose > 0:
382 print('removed file %s' % fn)
385def plot_style():
386 """Set style of plots.
387 """
388 plt.rcParams['figure.facecolor'] = 'white'
389 plt.rcParams['axes.facecolor'] = 'none'
390 plt.rcParams['xtick.direction'] = 'out'
391 plt.rcParams['ytick.direction'] = 'out'
394def axes_style(ax):
395 """Fix style of axes.
397 Parameters
398 ----------
399 ax: matplotlib axes
400 """
401 ax.spines['top'].set_visible(False)
402 ax.spines['right'].set_visible(False)
403 ax.get_xaxis().tick_bottom()
404 ax.get_yaxis().tick_left()
407def plot_eods(base_name, message_filename,
408 raw_data, samplerate, channel, idx0, idx1, clipped,
409 psd_data, wave_eodfs, wave_indices, mean_eods, eod_props,
410 peak_data, spec_data, indices, unit, zoom_window,
411 n_snippets=10, power_thresh=None, label_power=True,
412 all_eods=False, spec_plots='auto', skip_bad=True,
413 log_freq=False, min_freq=0.0, max_freq=3000.0,
414 interactive=True, verbose=0):
415 """Creates an output plot for the thunderfish program.
417 This output contains the raw trace where the analysis window is
418 marked, the power-spectrum of this analysis window where the
419 detected fish are marked, plots of averaged EOD plots, and
420 spectra of the EOD waveforms.
422 Parameters
423 ----------
424 base_name: string
425 Basename of audio_file.
426 message_filename: string or None
427 Path to meta-data message.
428 raw_data: array
429 Dataset.
430 samplerate: float
431 Sampling rate of the dataset.
432 channel: int or None
433 Channel of the recording to be put into the plot title.
434 If None, do not write the channel into the title.
435 idx0: float
436 Index of the beginning of the analysis window in the dataset.
437 idx1: float
438 Index of the end of the analysis window in the dataset.
439 clipped: float
440 Fraction of clipped amplitudes.
441 psd_data: 2D array
442 Power spectrum (frequencies and power) of the analysed data.
443 wave_eodfs: array
444 Frequency and power of fundamental frequency/harmonics of several fish.
445 wave_indices: array of int
446 Indices of wave fish mapping from wave_eodfs to eod_props.
447 If negative, then that EOD frequency has no waveform described in eod_props.
448 mean_eods: list of 2-D arrays with time, mean and std.
449 Mean trace for the mean EOD plot.
450 eod_props: list of dict
451 Properties for each waveform in mean_eods.
452 peak_data: list of 2_D arrays
453 For each pulsefish a list of peak properties
454 (index, time, and amplitude).
455 spec_data: list of 2_D arrays
456 For each pulsefish a power spectrum of the single pulse and for
457 each wavefish the relative amplitudes and phases of the harmonics.
458 indices: list of int or None
459 Indices of the fish in eod_props to be plotted.
460 If None try to plot all.
461 unit: string
462 Unit of the trace and the mean EOD.
463 n_snippets: int
464 Number of EOD waveform snippets to be plotted. If zero do not plot any.
465 power_thresh: 2 D array or None
466 Frequency (first column) and power (second column) of threshold
467 derived from single pulse spectra to discard false wave fish.
468 label_power: boolean
469 If `True` put the power in decibel in addition to the frequency
470 into the legend.
471 all_eods: bool
472 Plot all EOD waveforms.
473 spec_plots: bool or 'auto'
474 Plot amplitude spectra of EOD waveforms.
475 If 'auto', plot them if there is a singel waveform only.
476 skip_bad: bool
477 Skip harmonic groups without index (entry in indices is negative).
478 log_freq: boolean
479 Logarithmic (True) or linear (False) frequency axis of power spectrum of recording.
480 min_freq: float
481 Limits of frequency axis of power spectrum of recording
482 are set to `(min_freq, max_freq)` if `max_freq` is greater than zero
483 max_freq: float
484 Limits of frequency axis of power spectrum of recording
485 are set to `(min_freq, max_freq)` and limits of power axis are computed
486 from powers below max_freq if `max_freq` is greater than zero
487 interactive: bool
488 If True install some keyboard interaction.
489 verbose: int
490 Print out information about data to be plotted if greater than zero.
492 Returns
493 -------
494 fig: plt.figure
495 Figure with the plots.
496 """
498 def keypress(event):
499 if event.key in 'pP':
500 if idx1 > idx0:
501 playdata = 1.0 * raw_data[idx0:idx1]
502 else:
503 playdata = 1.0 * raw_data[:]
504 fade(playdata, samplerate, 0.1)
505 play(playdata, samplerate, blocking=False)
506 if event.key in 'mM' and message_filename:
507 # play voice message:
508 msg, msg_rate = load_audio(message_filename)
509 play(msg, msg_rate, blocking=False)
511 def recording_format_coord(x, y):
512 return 'full recording: x=%.3f s, y=%.3f' % (x, y)
514 def recordingzoom_format_coord(x, y):
515 return 'recording zoom-in: x=%.3f s, y=%.3f' % (x, y)
517 def psd_format_coord(x, y):
518 return 'power spectrum: x=%.1f Hz, y=%.1f dB' % (x, y)
520 def meaneod_format_coord(x, y):
521 return 'mean EOD waveform: x=%.2f ms, y=%.3f' % (x, y)
523 def ampl_format_coord(x, y):
524 return u'amplitude spectrum: x=%.0f, y=%.2f' % (x, y)
526 def phase_format_coord(x, y):
527 return u'phase spectrum: x=%.0f, y=%.2f \u03c0' % (x, y/np.pi)
529 def pulsepsd_format_coord(x, y):
530 return 'single pulse power spectrum: x=%.1f Hz, y=%.1f dB' % (x, y)
532 # count number of fish types to be plotted:
533 if indices is None:
534 indices = np.arange(len(eod_props))
535 else:
536 indices = np.array(indices, dtype=int)
537 nwave = 0
538 npulse = 0
539 for idx in indices:
540 if eod_props[idx]['type'] == 'pulse':
541 npulse += 1
542 elif eod_props[idx]['type'] == 'wave':
543 nwave += 1
544 neods = nwave + npulse
546 if verbose > 0:
547 print('plot: %2d waveforms: %2d wave fish, %2d pulse fish and %2d EOD frequencies.'
548 % (len(indices), nwave, npulse, len(wave_eodfs)))
550 # size and positions:
551 if spec_plots == 'auto':
552 spec_plots = len(indices) == 1
553 large_plots = spec_plots or len(indices) <= 2
554 width = 14.0
555 height = 10.0
556 if all_eods and len(indices) > 0:
557 nrows = len(indices) if spec_plots else (len(indices)+1)//2
558 if large_plots:
559 height = 6.0 + 4.0*nrows
560 else:
561 height = 6.4 + 1.9*nrows
562 leftx = 1.0/width
563 midx = 0.5 + leftx
564 fullwidth = 1.0-1.4/width
565 halfwidth = 0.5-1.4/width
566 pheight = 3.0/height
568 # figure:
569 plot_style()
570 fig = plt.figure(figsize=(width, height))
571 if interactive:
572 fig.canvas.mpl_connect('key_press_event', keypress)
574 # plot title:
575 title = base_name
576 if channel is not None:
577 title += ' c%d' % channel
578 ax = fig.add_axes([0.2/width, 1.0-0.6/height, 1.0-0.4/width, 0.55/height])
579 ax.text(0.0, 1.0, title, fontsize=22, va='top')
580 ax.text(1.0, 1.0, 'thunderfish by Benda-Lab', fontsize=16, ha='right', va='top')
581 ax.text(1.0, 0.0, 'version %s' % __version__, fontsize=16, ha='right', va='bottom')
582 ax.set_frame_on(False)
583 ax.set_axis_off()
584 ax.set_navigate(False)
586 # layout of recording and psd plots:
587 #force_both = True # set to True for debugging pulse and wave detection
588 force_both = False
589 posy = 1.0 - 4.0/height
590 axr = None
591 axp = None
592 legend_inside = True
593 legendwidth = 2.2/width if label_power else 1.7/width
594 if neods == 0:
595 axr = fig.add_axes([leftx, posy, fullwidth, pheight]) # top, wide
596 if len(psd_data) > 0:
597 axp = fig.add_axes([leftx, 2.0/height, fullwidth, pheight]) # bottom, wide
598 else:
599 if npulse == 0 and nwave > 2 and psd_data is not None and \
600 len(psd_data) > 0 and not force_both:
601 axp = fig.add_axes([leftx, posy, fullwidth-legendwidth, pheight]) # top, wide
602 legend_inside = False
603 elif (npulse > 0 or psd_data is None or len(psd_data) == 0) \
604 and len(wave_eodfs) == 0 and not force_both:
605 axr = fig.add_axes([leftx, posy, fullwidth, pheight]) # top, wide
606 else:
607 axr = fig.add_axes([leftx, posy, halfwidth, pheight]) # top left
608 label_power = False
609 legendwidth = 2.2/width
610 axp = fig.add_axes([midx, posy, halfwidth, pheight]) # top, right
612 # best window data:
613 data = raw_data[idx0:idx1] if idx1 > idx0 else raw_data
615 # plot recording
616 pulse_colors, pulse_markers = colors_markers()
617 pulse_colors = pulse_colors[3:]
618 pulse_markers = pulse_markers[3:]
619 if axr is not None:
620 axes_style(axr)
621 twidth = 0.1
622 if len(indices) > 0:
623 if eod_props[indices[0]]['type'] == 'wave':
624 twidth = 5.0/eod_props[indices[0]]['EODf']
625 else:
626 if len(wave_eodfs) > 0:
627 twidth = 3.0/eod_props[indices[0]]['EODf']
628 else:
629 twidth = 10.0/eod_props[indices[0]]['EODf']
630 twidth = (1+twidth//0.005)*0.005
631 if data is not None and len(data) > 0:
632 plot_eod_recording(axr, data, samplerate, unit, twidth,
633 idx0/samplerate)
634 plot_pulse_eods(axr, data, samplerate,
635 zoom_window, twidth, eod_props,
636 idx0/samplerate, colors=pulse_colors,
637 markers=pulse_markers, frameon=True,
638 loc='upper right')
639 if axr.get_legend() is not None:
640 axr.get_legend().get_frame().set_color('white')
641 axr.set_title('Recording', fontsize=14, y=1.05)
642 axr.format_coord = recordingzoom_format_coord
644 # plot psd
645 wave_colors, wave_markers = colors_markers()
646 if axp is not None:
647 axes_style(axp)
648 if power_thresh is not None:
649 axp.plot(power_thresh[:,0], decibel(power_thresh[:,1]), '#CCCCCC', lw=1)
650 if len(wave_eodfs) > 0:
651 kwargs = {}
652 if len(wave_eodfs) > 1:
653 title = '%d EOD frequencies' % len(wave_eodfs)
654 kwargs = {'title': title if len(wave_eodfs) > 2 else None }
655 if legend_inside:
656 kwargs.update({'bbox_to_anchor': (1.05, 1.1),
657 'loc': 'upper right', 'legend_rows': 10,
658 'frameon': True})
659 else:
660 kwargs.update({'bbox_to_anchor': (1.02, 1.1),
661 'loc': 'upper left', 'legend_rows': 14,
662 'labelspacing': 0.6, 'frameon': False})
663 plot_harmonic_groups(axp, wave_eodfs, wave_indices, max_groups=0,
664 skip_bad=skip_bad,
665 sort_by_freq=True, label_power=label_power,
666 colors=wave_colors, markers=wave_markers,
667 **kwargs)
668 if legend_inside:
669 axp.get_legend().get_frame().set_color('white')
670 if psd_data is not None and len(psd_data) > 0:
671 plot_decibel_psd(axp, psd_data[:,0], psd_data[:,1], log_freq=log_freq,
672 min_freq=min_freq, max_freq=max_freq, ymarg=5.0, color='blue')
673 axp.yaxis.set_major_locator(ticker.MaxNLocator(6))
674 if len(wave_eodfs) == 1:
675 axp.get_legend().set_visible(False)
676 label = '%6.1f Hz' % wave_eodfs[0][0, 0]
677 axp.set_title('Powerspectrum: %s' % label, y=1.05, fontsize=14)
678 else:
679 axp.set_title('Powerspectrum', y=1.05, fontsize=14)
680 axp.format_coord = psd_format_coord
682 # get fish labels from legends:
683 if axp is not None:
684 w, _ = axp.get_legend_handles_labels()
685 eodf_labels = [wi.get_label().split()[0] for wi in w]
686 legend_wave_eodfs = np.array([float(f) if f[0] != '(' else np.nan for f in eodf_labels])
687 if axr is not None:
688 p, _ = axr.get_legend_handles_labels()
689 eodf_labels = [pi.get_label().split()[0] for pi in p]
690 legend_pulse_eodfs = np.array([float(f) if f[0] != '(' else np.nan for f in eodf_labels])
692 # layout:
693 sheight = 1.4/height
694 sstep = 1.6/height
695 max_plots = len(indices)
696 if not all_eods:
697 if large_plots:
698 max_plots = 1 if spec_plots else 2
699 else:
700 max_plots = 4
701 if large_plots:
702 pstep = pheight + 1.0/height
703 ty = 1.08
704 my = 1.10
705 ny = 6
706 else:
707 posy -= 0.2/height
708 pheight = 1.3/height
709 pstep = 1.9/height
710 ty = 1.10
711 my = 1.16
712 ny = 4
713 posy -= pstep
715 # sort indices by p-p amplitude:
716 pp_ampls = [eod_props[idx]['p-p-amplitude'] for idx in indices]
717 pp_indices = np.argsort(pp_ampls)[::-1]
719 # plot EOD waveform and spectra:
720 for k, idx in enumerate(indices[pp_indices]):
721 if k >= max_plots:
722 break
723 # plot EOD waveform:
724 mean_eod = mean_eods[idx]
725 props = eod_props[idx]
726 peaks = peak_data[idx]
727 lx = leftx if spec_plots or k%2 == 0 else midx
728 ax = fig.add_axes([lx, posy, halfwidth, pheight])
729 axes_style(ax)
730 ax.yaxis.set_major_locator(ticker.MaxNLocator(ny))
731 if len(indices) > 1:
732 ax.text(0.3, ty, '{EODf:.1f} Hz {type} fish'.format(**props),
733 transform=ax.transAxes, fontsize=14, zorder=20)
734 mx = 0.25
735 else:
736 ax.text(-0.1, ty, '{EODf:.1f} Hz {type} fish'.format(**props),
737 transform=ax.transAxes, fontsize=14, zorder=20)
738 ax.text(0.5, ty, 'Averaged EOD', ha='center',
739 transform=ax.transAxes, fontsize=14, zorder=20)
740 mx = -0.14
741 eodf = props['EODf']
742 if props['type'] == 'wave':
743 if axp is not None:
744 wk = np.nanargmin(np.abs(legend_wave_eodfs - eodf))
745 ma = ml.Line2D([mx], [my], color=w[wk].get_color(), marker=w[wk].get_marker(),
746 markersize=w[wk].get_markersize(), mec='none', clip_on=False,
747 label=w[wk].get_label(), transform=ax.transAxes)
748 ax.add_line(ma)
749 else:
750 if axr is not None and len(legend_pulse_eodfs) > 0:
751 pk = np.argmin(np.abs(legend_pulse_eodfs - eodf))
752 ma = ml.Line2D([mx], [my], color=p[pk].get_color(), marker=p[pk].get_marker(),
753 markersize=p[pk].get_markersize(), mec='none', clip_on=False,
754 label=p[pk].get_label(), transform=ax.transAxes)
755 ax.add_line(ma)
756 plot_eod_waveform(ax, mean_eod, props, peaks, unit)
757 if props['type'] == 'pulse' and 'times' in props:
758 plot_eod_snippets(ax, data, samplerate, mean_eod[0,0], mean_eod[-1,0],
759 props['times'], n_snippets, props['flipped'])
760 if not large_plots and k < max_plots-2:
761 ax.set_xlabel('')
762 ax.format_coord = meaneod_format_coord
764 # plot spectra:
765 if spec_plots:
766 spec = spec_data[idx]
767 if props['type'] == 'pulse':
768 ax = fig.add_axes([midx, posy, halfwidth, pheight])
769 axes_style(ax)
770 plot_pulse_spectrum(ax, spec, props)
771 ax.set_title('Single pulse spectrum', fontsize=14, y=1.05)
772 ax.format_coord = pulsepsd_format_coord
773 else:
774 axa = fig.add_axes([midx, posy+sstep, halfwidth, sheight])
775 axes_style(axa)
776 axp = fig.add_axes([midx, posy, halfwidth, sheight])
777 axes_style(axp)
778 plot_wave_spectrum(axa, axp, spec, props, unit)
779 axa.set_title('Amplitude and phase spectrum', fontsize=14, y=1.05)
780 axa.set_xticklabels([])
781 axa.yaxis.set_major_locator(ticker.MaxNLocator(4))
782 axa.format_coord = ampl_format_coord
783 axp.format_coord = phase_format_coord
785 if spec_plots or k%2 == 1:
786 posy -= pstep
788 # whole trace:
789 ax = fig.add_axes([leftx, 0.6/height, fullwidth, 0.9/height])
790 axes_style(ax)
791 if raw_data is not None and len(raw_data) > 0:
792 plot_data_window(ax, raw_data, samplerate, unit, idx0, idx1, clipped)
793 ax.format_coord = recording_format_coord
795 return fig
798def plot_eod_subplots(base_name, subplots, raw_data, samplerate, idx0, idx1,
799 clipped, psd_data, wave_eodfs, wave_indices, mean_eods,
800 eod_props, peak_data, spec_data, unit, zoom_window,
801 n_snippets=10, power_thresh=None, label_power=True,
802 skip_bad=True, log_freq=False,
803 min_freq=0.0, max_freq=3000.0, save=True):
804 """Plot time traces and spectra into separate windows or files.
806 Parameters
807 ----------
808 base_name: string
809 Basename of audio_file.
810 subplots: string
811 Specifies which subplots to plot:
812 r) recording with best window, t) data trace with detected pulse fish,
813 p) power spectrum with detected wave fish, w/W) mean EOD waveform,
814 s/S) EOD spectrum, e/E) EOD waveform and spectra. With capital letters
815 all fish are saved into a single pdf file, with small letters each fish
816 is saved into a separate file.
817 raw_data: array
818 Dataset.
819 samplerate: float
820 Sampling rate of the dataset.
821 idx0: float
822 Index of the beginning of the analysis window in the dataset.
823 idx1: float
824 Index of the end of the analysis window in the dataset.
825 clipped: float
826 Fraction of clipped amplitudes.
827 psd_data: 2D array
828 Power spectrum (frequencies and power) of the analysed data.
829 wave_eodfs: array
830 Frequency and power of fundamental frequency/harmonics of several fish.
831 wave_indices: array of int
832 Indices of wave fish mapping from wave_eodfs to eod_props.
833 If negative, then that EOD frequency has no waveform described in eod_props.
834 mean_eods: list of 2-D arrays with time, mean and std.
835 Mean trace for the mean EOD plot.
836 eod_props: list of dict
837 Properties for each waveform in mean_eods.
838 peak_data: list of 2_D arrays
839 For each pulsefish a list of peak properties
840 (index, time, and amplitude).
841 spec_data: list of 2_D arrays
842 For each pulsefish a power spectrum of the single pulse and for
843 each wavefish the relative amplitudes and phases of the harmonics.
844 unit: string
845 Unit of the trace and the mean EOD.
846 n_snippets: int
847 Number of EOD waveform snippets to be plotted. If zero do not plot any.
848 power_thresh: 2 D array or None
849 Frequency (first column) and power (second column) of threshold
850 derived from single pulse spectra to discard false wave fish.
851 label_power: boolean
852 If `True` put the power in decibel in addition to the frequency
853 into the legend.
854 skip_bad: bool
855 Skip harmonic groups without index (entry in indices is negative).
856 log_freq: boolean
857 Logarithmic (True) or linear (False) frequency axis of power spectrum of recording.
858 min_freq: float
859 Limits of frequency axis of power spectrum of recording
860 are set to `(min_freq, max_freq)` if `max_freq` is greater than zero
861 max_freq: float
862 Limits of frequency axis of power spectrum of recording
863 are set to `(min_freq, max_freq)` and limits of power axis are computed
864 from powers below max_freq if `max_freq` is greater than zero
865 save: bool
866 If True save plots to files instead of showing them.
867 """
868 plot_style()
869 if 'r' in subplots:
870 fig, ax = plt.subplots(figsize=(10, 2))
871 fig.subplots_adjust(left=0.07, right=0.99, bottom=0.22, top=0.95)
872 plot_data_window(ax, raw_data, samplerate, unit, idx0, idx1, clipped)
873 ax.yaxis.set_major_locator(ticker.MaxNLocator(5))
874 axes_style(ax)
875 if save:
876 fig.savefig(base_name + '-recording.pdf')
877 if 't' in subplots:
878 fig, ax = plt.subplots(figsize=(10, 6))
879 twidth = 0.1
880 if len(eod_props) > 0:
881 if eod_props[0]['type'] == 'wave':
882 twidth = 5.0/eod_props[0]['EODf']
883 else:
884 if len(wave_eodfs) > 0:
885 twidth = 3.0/eod_props[0]['EODf']
886 else:
887 twidth = 10.0/eod_props[0]['EODf']
888 twidth = (1+twidth//0.005)*0.005
889 pulse_colors, pulse_markers = colors_markers()
890 pulse_colors = pulse_colors[3:]
891 pulse_markers = pulse_markers[3:]
892 plot_eod_recording(ax, raw_data[idx0:idx1], samplerate, unit,
893 twidth, idx0/samplerate)
894 plot_pulse_eods(ax, raw_data[idx0:idx1], samplerate,
895 zoom_window, twidth, eod_props,
896 idx0/samplerate, colors=pulse_colors,
897 markers=pulse_markers, frameon=True,
898 loc='upper right')
899 if ax.get_legend() is not None:
900 ax.get_legend().get_frame().set_color('white')
901 axes_style(ax)
902 if save:
903 fig.savefig(base_name + '-trace.pdf')
904 if 'p' in subplots:
905 fig, ax = plt.subplots(figsize=(10, 5))
906 fig.subplots_adjust(left=0.08, right=0.975, bottom=0.11, top=0.9)
907 axes_style(ax)
908 if power_thresh is not None:
909 ax.plot(power_thresh[:,0], decibel(power_thresh[:,1]), '#CCCCCC', lw=1)
910 if len(wave_eodfs) > 0:
911 kwargs = {}
912 if len(wave_eodfs) > 1:
913 title = '%d EOD frequencies' % len(wave_eodfs)
914 kwargs = {'title': title if len(wave_eodfs) > 2 else None }
915 if len(wave_eodfs) > 2:
916 fig.subplots_adjust(left=0.08, right=0.78, bottom=0.11, top=0.9)
917 kwargs.update({'bbox_to_anchor': (1.01, 1.1),
918 'loc': 'upper left', 'legend_rows': 14,
919 'labelspacing': 0.6})
920 else:
921 kwargs.update({'bbox_to_anchor': (1.05, 1.1),
922 'loc': 'upper right', 'legend_rows': 10})
923 wave_colors, wave_markers = colors_markers()
924 plot_harmonic_groups(ax, wave_eodfs, wave_indices, max_groups=0,
925 skip_bad=skip_bad,
926 sort_by_freq=True, label_power=label_power,
927 colors=wave_colors, markers=wave_markers,
928 frameon=False, **kwargs)
929 plot_decibel_psd(ax, psd_data[:,0], psd_data[:,1], log_freq=log_freq,
930 min_freq=min_freq, max_freq=max_freq, ymarg=5.0, color='blue')
931 ax.yaxis.set_major_locator(ticker.MaxNLocator(6))
932 if len(wave_eodfs) == 1:
933 ax.get_legend().set_visible(False)
934 label = '%6.1f Hz' % wave_eodfs[0][0, 0]
935 ax.set_title('Powerspectrum: %s' % label, y=1.05)
936 else:
937 ax.set_title('Powerspectrum', y=1.05)
938 if save:
939 fig.savefig(base_name + '-psd.pdf')
940 if 'w' in subplots or 'W' in subplots:
941 mpdf = None
942 if 'W' in subplots:
943 mpdf = PdfPages(base_name + '-waveforms.pdf')
944 for meod, props, peaks in zip(mean_eods, eod_props, peak_data):
945 if meod is None:
946 continue
947 fig, ax = plt.subplots(figsize=(5, 3))
948 fig.subplots_adjust(left=0.18, right=0.98, bottom=0.15, top=0.9)
949 if not props is None:
950 ax.set_title('{index:d}: {EODf:.1f} Hz {type} fish'.format(**props))
951 plot_eod_waveform(ax, meod, props, peaks, unit)
952 data = raw_data[idx0:idx1] if idx1 > idx0 else raw_data
953 if not props is None and props['type'] == 'pulse' and \
954 'times' in props:
955 plot_eod_snippets(ax, data, samplerate, meod[0,0], meod[-1,0],
956 props['times'], n_snippets)
957 ax.yaxis.set_major_locator(ticker.MaxNLocator(6))
958 axes_style(ax)
959 if mpdf is None:
960 if save:
961 fig.savefig(base_name + '-waveform-%d.pdf' % props['index'])
962 else:
963 mpdf.savefig(fig)
964 if mpdf is not None:
965 mpdf.close()
966 if 's' in subplots or 'S' in subplots:
967 mpdf = None
968 if 'S' in subplots:
969 mpdf = PdfPages(base_name + '-spectrum.pdf')
970 for props, peaks, spec in zip(eod_props, peak_data, spec_data):
971 if spec is None:
972 continue
973 if props is not None and props['type'] == 'pulse':
974 fig, ax = plt.subplots(figsize=(5, 3.5))
975 fig.subplots_adjust(left=0.15, right=0.967, bottom=0.16, top=0.88)
976 axes_style(ax)
977 ax.set_title('{index:d}: {EODf:.1f} Hz {type} fish'.format(**props), y=1.07)
978 plot_pulse_spectrum(ax, spec, props)
979 else:
980 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(5, 3.5))
981 fig.subplots_adjust(left=0.15, right=0.97, bottom=0.16, top=0.88, hspace=0.4)
982 axes_style(ax1)
983 axes_style(ax2)
984 if not props is None:
985 ax1.set_title('{index:d}: {EODf:.1f} Hz {type} fish'.format(**props), y=1.15)
986 plot_wave_spectrum(ax1, ax2, spec, props, unit)
987 ax1.set_xticklabels([])
988 ax1.yaxis.set_major_locator(ticker.MaxNLocator(4))
989 if mpdf is None:
990 if save:
991 fig.savefig(base_name + '-spectrum-%d.pdf' % props['index'])
992 else:
993 mpdf.savefig(fig)
994 if mpdf is not None:
995 mpdf.close()
996 if 'e' in subplots or 'E' in subplots:
997 mpdf = None
998 if 'E' in subplots:
999 mpdf = PdfPages(base_name + '-eods.pdf')
1000 for meod, props, peaks, spec in zip(mean_eods, eod_props, peak_data, spec_data):
1001 if meod is None or spec is None:
1002 continue
1003 fig = plt.figure(figsize=(10, 3.5))
1004 gs = gridspec.GridSpec(nrows=2, ncols=2, left=0.09, right=0.98,
1005 bottom=0.16, top=0.88, wspace=0.4, hspace=0.4)
1006 ax1 = fig.add_subplot(gs[:,0])
1007 if not props is None:
1008 ax1.set_title('{index:d}: {EODf:.1f} Hz {type} fish'.format(**props), y=1.07)
1009 plot_eod_waveform(ax1, meod, props, peaks, unit)
1010 data = raw_data[idx0:idx1] if idx1 > idx0 else raw_data
1011 if not props is None and props['type'] == 'pulse' and 'times' in props:
1012 plot_eod_snippets(ax1, data, samplerate, meod[0,0], meod[-1,0],
1013 props['times'], n_snippets)
1014 ax1.yaxis.set_major_locator(ticker.MaxNLocator(6))
1015 axes_style(ax1)
1016 if not props is None and props['type'] == 'pulse':
1017 ax2 = fig.add_subplot(gs[:,1])
1018 axes_style(ax2)
1019 plot_pulse_spectrum(ax2, spec, props)
1020 ax2.set_title('Single pulse spectrum', y=1.07)
1021 else:
1022 ax2 = fig.add_subplot(gs[0,1])
1023 ax3 = fig.add_subplot(gs[1,1])
1024 axes_style(ax2)
1025 axes_style(ax3)
1026 plot_wave_spectrum(ax2, ax3, spec, props, unit)
1027 ax2.set_title('Amplitude and phase spectrum', y=1.15)
1028 ax2.set_xticklabels([])
1029 ax2.yaxis.set_major_locator(ticker.MaxNLocator(4))
1030 if mpdf is None:
1031 if save:
1032 fig.savefig(base_name + '-eod-%d.pdf' % props['index'])
1033 else:
1034 mpdf.savefig(fig)
1035 if mpdf is not None:
1036 mpdf.close()
1037 if not save:
1038 plt.show()
1039 plt.close('all')
1042def thunderfish_plot(files, data_path=None, load_kwargs={},
1043 all_eods=False, spec_plots='auto', skip_bad=True,
1044 save_plot=False, multi_pdf=None,
1045 save_subplots='', log_freq=False, min_freq=0.0,
1046 max_freq=3000.0, output_folder='.',
1047 keep_path=False, verbose=0):
1048 """Generate plots from saved analysis results.
1050 Parameters
1051 ----------
1052 files: list of str
1053 Analysis files from a single recording.
1054 data_path: str
1055 Path where to find the raw data.
1056 load_kwargs: dict
1057 Key-word arguments for the `load_data()` function.
1058 all_eods: bool
1059 If True, plot all EOD waveforms.
1060 spec_plots: bool or 'auto'
1061 Plot amplitude spectra of EOD waveforms.
1062 If 'auto', plot them if there is a singel waveform only.
1063 skip_bad: bool
1064 Skip harmonic groups without index in the spectrum plot.
1065 save_plot: bool
1066 If True, save plots as pdf file.
1067 multi_pdf: matplotlib.PdfPages or None
1068 PdfPages instance in which to save plots.
1069 save_subplots: string
1070 If not empty, specifies subplots to be saved as separate pdf
1071 files: r) recording with best window, t) data trace with
1072 detected pulse fish, p) power spectrum with detected wave
1073 fish, w/W) mean EOD waveform, s/S) EOD spectrum, e/E) EOD
1074 waveform and spectra. Capital letters produce a single
1075 multipage pdf containing plots of all detected fish.
1076 log_freq: boolean
1077 Logarithmic (True) or linear (False) frequency axis of
1078 power spectrum of recording.
1079 min_freq: float
1080 Limits of frequency axis of power spectrum of recording
1081 are set to `(min_freq, max_freq)`, if `max_freq` is greater than zero
1082 max_freq: float
1083 Limits of frequency axis of power spectrum of recording
1084 are set to `(min_freq, max_freq)` and limits of power axis are computed
1085 from powers below max_freq, if `max_freq` is greater than zero
1086 output_folder: string
1087 Folder where to save results.
1088 keep_path: bool
1089 Add relative path of data files to output path.
1090 verbose: int
1091 Verbosity level (for debugging).
1092 """
1093 if len(save_subplots) == 0:
1094 save_subplots = 'rtpwsed' # plot everything
1095 # load analysis results:
1096 mean_eods, wave_eodfs, wave_indices, eod_props, spec_data, \
1097 peak_data, base_name, channel, unit = load_analysis(files)
1098 if len(mean_eods) == 0 or all(me is None for me in mean_eods):
1099 save_subplots = save_subplots.replace('w', '')
1100 save_subplots = save_subplots.replace('W', '')
1101 save_subplots = save_subplots.replace('e', '')
1102 save_subplots = save_subplots.replace('E', '')
1103 save_subplots = save_subplots.replace('d', '')
1104 if len(spec_data) == 0 or all(sd is None for sd in spec_data):
1105 save_subplots = save_subplots.replace('s', '')
1106 save_subplots = save_subplots.replace('S', '')
1107 save_subplots = save_subplots.replace('e', '')
1108 save_subplots = save_subplots.replace('E', '')
1109 save_subplots = save_subplots.replace('d', '')
1110 clipped = 0.0
1111 if len(eod_props) > 0 and not eod_props[0] is None and \
1112 'winclipped' in eod_props[0]:
1113 clipped = eod_props[0]['winclipped']
1114 zoom_window = [1.2, 1.3]
1115 # load recording:
1116 psd_data = None
1117 if base_name:
1118 name = os.path.basename(base_name) if data_path and data_path != '.' else base_name
1119 data_path = os.path.join(data_path, name)
1120 data, samplerate, idx0, idx1, data_path = \
1121 load_recording(data_path, channel, load_kwargs,
1122 eod_props, verbose-1)
1123 if data is None:
1124 save_subplots = save_subplots.replace('r', '')
1125 save_subplots = save_subplots.replace('t', '')
1126 save_subplots = save_subplots.replace('d', '')
1127 if verbose > 0:
1128 print('loaded', data_path)
1129 if len(eod_props) > 0 and not eod_props[0] is None and \
1130 'dfreq' in eod_props[0] and data is not None and len(data) > 0:
1131 psd_data = multi_psd(data[idx0:idx1],
1132 samplerate,
1133 1.1*eod_props[0]['dfreq'])[0]
1134 if psd_data is not None and len(psd_data) > 0:
1135 for idx, fish in zip(wave_indices, wave_eodfs):
1136 if idx < 0:
1137 for k in range(len(fish)):
1138 fish[k,1] = psd_data[np.argmin(np.abs(psd_data[:,0] - fish[k,0])),1]
1139 if psd_data is None:
1140 save_subplots = save_subplots.replace('p', '')
1141 save_subplots = save_subplots.replace('d', '')
1142 # file name for output files:
1143 fn = base_name if keep_path else os.path.basename(base_name)
1144 output_basename = os.path.join(output_folder, fn)
1145 if channel >= 0:
1146 output_basename += f'-c{channel}'
1147 # make directory if necessary:
1148 if keep_path:
1149 outpath = os.path.dirname(output_basename)
1150 if not os.path.exists(outpath):
1151 if verbose > 0:
1152 print('mkdir %s' % outpath)
1153 os.makedirs(outpath)
1154 # plot:
1155 if len(save_subplots) == 0 or 'd' in save_subplots:
1156 fig = plot_eods(os.path.basename(base_name), None,
1157 data, samplerate,
1158 channel, idx0, idx1, clipped, psd_data,
1159 wave_eodfs, wave_indices, mean_eods,
1160 eod_props, peak_data, spec_data, None, unit,
1161 zoom_window, 10, None, True, all_eods,
1162 spec_plots, skip_bad, log_freq, min_freq,
1163 max_freq, interactive=not save_plot,
1164 verbose=verbose-1)
1165 if save_plot:
1166 if multi_pdf is not None:
1167 multi_pdf.savefig(fig)
1168 else:
1169 fig.savefig(output_basename + '.pdf')
1170 else:
1171 fig.canvas.manager.set_window_title('thunderfish')
1172 plt.show()
1173 plt.close()
1174 save_subplots = save_subplots.replace('d', '')
1175 if len(save_subplots) > 0:
1176 plot_eod_subplots(output_basename, save_subplots, data,
1177 samplerate, idx0, idx1, clipped,
1178 psd_data, wave_eodfs, wave_indices,
1179 mean_eods, eod_props, peak_data,
1180 spec_data, unit, zoom_window, 10, None,
1181 True, skip_bad, log_freq, min_freq,
1182 max_freq, save_plot)
1183 return None
1186def thunderfish(filename, load_kwargs, cfg, channel=0,
1187 time=None, time_file=False,
1188 mode='wp', log_freq=False, min_freq=0.0, max_freq=3000,
1189 save_data=False, zip_file=False,
1190 all_eods=False, spec_plots='auto', skip_bad=True,
1191 save_plot=False, multi_pdf=None, save_subplots='',
1192 output_folder='.', keep_path=False,
1193 verbose=0, plot_level=0):
1194 """Automatically detect and analyze all EOD waveforms in a short recording.
1196 Parameters
1197 ----------
1198 filename: string
1199 Path of the data file to be analyzed.
1200 load_kwargs: dict
1201 Key-word arguments for the `load_data()` function.
1202 cfg: dict
1203 channel: int
1204 Channel to be analyzed.
1205 time: string, float, or None
1206 Start time of analysis window: "beginning", "center", "end",
1207 "best", or time in seconds (as float or string). If not None
1208 overwrites "windowPosition" in cofiguration file.
1209 time_file: bool
1210 If `True` add time of analysis window to output file names.
1211 mode: 'w', 'p', 'P', 'wp', or 'wP'
1212 Analyze wavefish ('w'), all pulse fish ('p'), or largest pulse
1213 fish only ('P').
1214 log_freq: boolean
1215 Logarithmic (True) or linear (False) frequency axis of
1216 power spectrum of recording.
1217 min_freq: float
1218 Limits of frequency axis of power spectrum of recording
1219 are set to `(min_freq, max_freq)`, if `max_freq` is greater than zero
1220 max_freq: float
1221 Limits of frequency axis of power spectrum of recording
1222 are set to `(min_freq, max_freq)` and limits of power axis are computed
1223 from powers below max_freq, if `max_freq` is greater than zero
1224 save_data: bool
1225 If True save analysis results in files. If False, just plot the data.
1226 zip_data: bool
1227 If True, store all analysis results in a single zip file.
1228 all_eods: bool
1229 If True, plot all EOD waveforms.
1230 spec_plots: bool or 'auto'
1231 Plot amplitude spectra of EOD waveforms.
1232 If 'auto', plot them if there is a singel waveform only.
1233 skip_bad: bool
1234 Skip harmonic groups without index in the spectrum plot.
1235 save_plot: bool
1236 If True, save plots as pdf file.
1237 multi_pdf: matplotlib.PdfPages or None
1238 PdfPages instance in which to save plots.
1239 save_subplots: string
1240 If not empty, specifies subplots to be saved as separate pdf
1241 files: r) recording with best window, t) data trace with
1242 detected pulse fish, p) power spectrum with detected wave
1243 fish, w/W) mean EOD waveform, s/S) EOD spectrum, e/E) EOD
1244 waveform and spectra. Capital letters produce a single
1245 multipage pdf containing plots of all detected fish.
1246 output_folder: string
1247 Folder where to save results.
1248 keep_path: bool
1249 Add relative path of data files to output path.
1250 verbose: int
1251 Verbosity level (for debugging).
1252 plot_level: int
1253 Plot intermediate results.
1255 Returns
1256 -------
1257 msg: string or None
1258 In case of errors, an error message.
1259 """
1260 # check data file:
1261 if len(filename) == 0:
1262 return 'you need to specify a file containing some data'
1264 # file names:
1265 fn = filename if keep_path else os.path.basename(filename)
1266 outfilename = os.path.splitext(fn)[0]
1267 messagefilename = os.path.splitext(fn)[0] + '-message.wav'
1268 if not os.path.isfile(messagefilename):
1269 messagefilename = None
1271 # load data:
1272 try:
1273 all_data, samplerate, unit, ampl_max = load_data(filename,
1274 verbose=verbose,
1275 **load_kwargs)
1276 except IOError as e:
1277 return '%s: failed to open file: %s' % (filename, str(e))
1278 # select channel:
1279 channels = all_data.shape[1]
1280 chan_list = [channel]
1281 if channel < 0:
1282 chan_list = range(channels)
1283 elif channel >= channels:
1284 return '%s: invalid channel %d (%d channels)' % (filename, channel, channels)
1285 # process all channels:
1286 for chan in chan_list:
1287 raw_data = all_data[:,chan]
1288 if len(raw_data) <= 1:
1289 return '%s: empty data file' % filename
1290 if verbose >= 0 and len(chan_list) > 1:
1291 print(' channel %d' % chan)
1293 # analysis window:
1294 win_pos = cfg.value('windowPosition')
1295 if time is not None:
1296 win_pos = time
1297 data, idx0, idx1, clipped, min_clip, max_clip = \
1298 analysis_window(raw_data, samplerate, ampl_max, win_pos,
1299 cfg, plot_level>0)
1300 found_bestwindow = idx1 > 0
1301 if not found_bestwindow:
1302 return '%s: not enough data for requested window length. You may want to adjust the windowSize parameter in the configuration file.' % filename
1304 # detect EODs in the data:
1305 psd_data, wave_eodfs, wave_indices, eod_props, \
1306 mean_eods, spec_data, peak_data, power_thresh, skip_reason, zoom_window = \
1307 detect_eods(data, samplerate, min_clip, max_clip, filename,
1308 mode, verbose, plot_level, cfg)
1309 if not found_bestwindow:
1310 wave_eodfs = []
1311 wave_indices = []
1312 eod_props = []
1313 mean_eods = []
1315 # add analysis window to EOD properties:
1316 for props in eod_props:
1317 props['twin'] = idx0/samplerate
1318 props['window'] = (idx1 - idx0)/samplerate
1319 props['winclipped'] = clipped
1321 # warning message in case no fish has been found:
1322 if found_bestwindow and not eod_props :
1323 msg = ', '.join(skip_reason)
1324 if msg:
1325 print(filename + ': no fish found: %s' % msg)
1326 else:
1327 print(filename + ': no fish found.')
1329 # file name for output files:
1330 output_basename = os.path.join(output_folder, outfilename)
1331 if channels > 1:
1332 if channels > 100:
1333 output_basename += '-c%03d' % chan
1334 elif channels > 10:
1335 output_basename += '-c%02d' % chan
1336 else:
1337 output_basename += '-c%d' % chan
1338 if time_file:
1339 output_basename += '-t%.0fs' % (idx0/samplerate)
1340 # make directory if necessary:
1341 if keep_path and found_bestwindow:
1342 outpath = os.path.dirname(output_basename)
1343 if not os.path.exists(outpath):
1344 if verbose > 0:
1345 print('mkdir %s' % outpath)
1346 os.makedirs(outpath)
1347 # save results to files:
1348 if save_data:
1349 remove_eod_files(output_basename, verbose, cfg)
1350 if found_bestwindow:
1351 save_analysis(output_basename, zip_file, eod_props,
1352 mean_eods, spec_data, peak_data,
1353 wave_eodfs, wave_indices, unit, verbose,
1354 **write_table_args(cfg))
1355 # summary plots:
1356 if save_plot or not save_data:
1357 n_snippets = 10
1358 if len(save_subplots) == 0 or 'd' in save_subplots:
1359 chl = chan if channels > 1 else None
1360 fig = plot_eods(outfilename, messagefilename,
1361 raw_data, samplerate,
1362 chl, idx0, idx1, clipped, psd_data[0],
1363 wave_eodfs, wave_indices, mean_eods,
1364 eod_props, peak_data, spec_data, None,
1365 unit, zoom_window, n_snippets,
1366 power_thresh, True, all_eods,
1367 spec_plots, skip_bad, log_freq,
1368 min_freq, max_freq, interactive=not
1369 save_plot, verbose=verbose)
1370 if save_plot:
1371 if multi_pdf is not None:
1372 multi_pdf.savefig(fig)
1373 else:
1374 fig.savefig(output_basename + '.pdf')
1375 else:
1376 fig.canvas.manager.set_window_title('thunderfish')
1377 plt.show()
1378 plt.close()
1379 save_subplots = save_subplots.replace('d', '')
1380 if len(save_subplots) > 0:
1381 plot_eod_subplots(output_basename, save_subplots,
1382 raw_data, samplerate, idx0, idx1,
1383 clipped, psd_data[0], wave_eodfs,
1384 wave_indices, mean_eods, eod_props,
1385 peak_data, spec_data, unit,
1386 zoom_window, n_snippets,
1387 power_thresh, True, skip_bad,
1388 log_freq, min_freq, max_freq,
1389 save_plot)
1390 return None
1393def run_thunderfish(file_args):
1394 """Helper function for mutlithreading Pool().map().
1395 """
1396 results = file_args[1][0]
1397 verbose = file_args[1][-1] if results else file_args[1][-2]+1
1398 if verbose > 1:
1399 print('='*70)
1400 try:
1401 if results:
1402 thunderfish_plot(file_args[0], *file_args[1][1:])
1403 else:
1404 if verbose > 0:
1405 print('analyze recording %s ...' % file_args[0])
1406 msg = thunderfish(file_args[0], *file_args[1][1:])
1407 if msg:
1408 print(msg)
1409 except (KeyboardInterrupt, SystemExit):
1410 print('\nthunderfish interrupted by user... exit now.')
1411 sys.exit(0)
1412 except:
1413 print(traceback.format_exc())
1416def main(cargs=None):
1417 # config file name:
1418 cfgfile = __package__ + '.cfg'
1420 # command line arguments:
1421 if cargs is None:
1422 cargs = sys.argv[1:]
1423 parser = argparse.ArgumentParser(add_help=False,
1424 description='Analyze EOD waveforms of weakly electric fish.',
1425 epilog='version %s by Benda-Lab (2015-%s)' % (__version__, __year__))
1426 parser.add_argument('-h', '--help', action='store_true',
1427 help='show this help message and exit')
1428 parser.add_argument('--version', action='version', version=__version__)
1429 parser.add_argument('-v', action='count', dest='verbose', default=0,
1430 help='verbosity level. Increase by specifying -v multiple times, or like -vvv')
1431 parser.add_argument('-V', action='count', dest='plot_level', default=0,
1432 help='level for debugging plots. Increase by specifying -V multiple times, or like -VVV')
1433 parser.add_argument('-c', dest='save_config', action='store_true',
1434 help='save configuration to file {0} after reading all configuration files'.format(cfgfile))
1435 parser.add_argument('--channel', default=0, type=int,
1436 help='channel to be analyzed (defaults to first channel, negative channel selects all channels)')
1437 parser.add_argument('-t', dest='time', default=None, type=str, metavar='TIME',
1438 help='start time of analysis window in recording: "beginning", "center", "end", "best", or time in seconds (overwrites "windowPosition" in cofiguration file)')
1439 parser.add_argument('-u', dest='unwrap', action='store_true',
1440 help='unwrap clipped files, toggles unwrap setting of config file.')
1441 parser.add_argument('-T', dest='time_file', action='store_true',
1442 help='add start time of analysis file to output file names')
1443 parser.add_argument('-m', dest='mode', default='wp', type=str,
1444 choices=['w', 'p', 'wp'],
1445 help='extract wave "w" and/or pulse "p" fish EODs')
1446 parser.add_argument('-a', dest='all_eods', action='store_true',
1447 help='show all EOD waveforms in the summary plot')
1448 parser.add_argument('-S', dest='spec_plots', action='store_true',
1449 help='plot spectra for all EOD waveforms in the summary plot')
1450 parser.add_argument('-b', dest='skip_bad', action='store_false',
1451 help='indicate bad EODs in legend of power spectrum')
1452 parser.add_argument('-l', dest='log_freq', type=float, metavar='MINFREQ',
1453 nargs='?', const=100.0, default=0.0,
1454 help='logarithmic frequency axis in power spectrum with optional minimum frequency (defaults to 100 Hz)')
1455 parser.add_argument('-p', dest='save_plot', action='store_true',
1456 help='save output plots as pdf files')
1457 parser.add_argument('-M', dest='multi_pdf', default='', type=str, metavar='PDFFILE',
1458 help='save all summary plots of all recordings in a multi page pdf file. Disables parallel jobs.')
1459 parser.add_argument('-P', dest='save_subplots', default='', type=str, metavar='rtpwsed',
1460 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')
1461 parser.add_argument('-d', dest='rawdata_path', default='.', type=str, metavar='PATH',
1462 help='path to raw EOD recordings needed for plotting based on analysis results')
1463 parser.add_argument('-j', dest='jobs', nargs='?', type=int, default=None, const=0,
1464 help='number of jobs run in parallel. Without argument use all CPU cores.')
1465 parser.add_argument('-s', dest='save_data', action='store_true',
1466 help='save analysis results to files')
1467 parser.add_argument('-z', dest='zip_file', action='store_true',
1468 help='save analysis results in a single zip file')
1469 parser.add_argument('-f', dest='format', default='auto', type=str,
1470 choices=TableData.formats + ['py'],
1471 help='file format used for saving analysis results, defaults to the format specified in the configuration file or "csv"')
1472 parser.add_argument('-o', dest='outpath', default='.', type=str,
1473 help='path where to store results and figures (defaults to current working directory)')
1474 parser.add_argument('-k', dest='keep_path', action='store_true',
1475 help='keep path of input file when saving analysis files, i.e. append path of input file to OUTPATH')
1476 parser.add_argument('-i', dest='load_kwargs', default=[],
1477 action='append', metavar='KWARGS',
1478 help='key-word arguments for the data loader function')
1479 parser.add_argument('file', nargs='*', default='', type=str,
1480 help='name of a file with time series data of an EOD recording, may include wildcards')
1481 args = parser.parse_args(cargs)
1483 # help:
1484 if args.help:
1485 parser.print_help()
1486 print('')
1487 print('examples:')
1488 print('- analyze the single file data.wav interactively:')
1489 print(' > thunderfish data.wav')
1490 print('- extract wavefish only:')
1491 print(' > thunderfish -m w data.wav')
1492 print('- automatically analyze all wav files in the current working directory and save analysis results and plot to files:')
1493 print(' > thunderfish -s -p *.wav')
1494 print('- analyze all wav files in the river1/ directory, use all CPUs, and write files directly to "results/":')
1495 print(' > thunderfish -j -s -p -o results/ river1/*.wav')
1496 print('- analyze all wav files in the river1/ directory and write files to "results/river1/":')
1497 print(' > thunderfish -s -p -o results/ -k river1/*.wav')
1498 print('- write configuration file:')
1499 print(' > thunderfish -c')
1500 parser.exit()
1502 # set verbosity level from command line:
1503 verbose = args.verbose
1504 plot_level = args.plot_level
1505 if verbose < plot_level:
1506 verbose = plot_level
1508 # interactive plot:
1509 plt.rcParams['keymap.quit'] = 'ctrl+w, alt+q, q'
1510 plt.ioff()
1512 # expand wildcard patterns:
1513 files = []
1514 if os.name == 'nt':
1515 for fn in args.file:
1516 files.extend(glob.glob(fn))
1517 else:
1518 files = [f for f in args.file if '-message' not in f]
1520 # save configuration:
1521 if args.save_config:
1522 file_name = files[0] if len(files) else ''
1523 cfg = configuration()
1524 cfg.load_files(cfgfile, file_name, 4, verbose)
1525 save_configuration(cfg, cfgfile)
1526 exit()
1527 elif len(files) == 0:
1528 parser.error('you need to specify at least one file for the analysis')
1530 # configure:
1531 cfg = configuration()
1532 cfg.load_files(cfgfile, files[0], 4, verbose)
1533 if args.format != 'auto':
1534 cfg.set('fileFormat', args.format)
1535 if args.unwrap:
1536 cfg.set('unwrapData', not cfg.value('unwrapData'))
1538 # plot parameter:
1539 spec_plots = 'auto'
1540 if args.spec_plots:
1541 spec_plots = True
1543 # multi-page pdfs:
1544 multi_pdf = None
1545 if len(args.multi_pdf) > 0:
1546 args.save_plot = True
1547 args.jobs = None # PdfPages does not work yet with mutliprocessing
1548 ext = os.path.splitext(args.multi_pdf)[1]
1549 if ext != os.extsep + 'pdf':
1550 args.multi_pdf += os.extsep + 'pdf'
1551 multi_pdf = PdfPages(args.multi_pdf)
1553 # create output folder:
1554 if args.save_data or args.save_plot:
1555 if not os.path.exists(args.outpath):
1556 if verbose > 1:
1557 print('mkdir %s' % args.outpath)
1558 os.makedirs(args.outpath)
1560 # kwargs for data loader:
1561 load_kwargs = {}
1562 for s in args.load_kwargs:
1563 for kw in s.split(','):
1564 kws = kw.split(':')
1565 if len(kws) == 2:
1566 load_kwargs[kws[0].strip()] = kws[1].strip()
1568 # frequency limits for power spectrum:
1569 min_freq = 0.0
1570 max_freq = 3000.0
1571 log_freq = args.log_freq
1572 if log_freq > 0.0:
1573 min_freq = log_freq
1574 max_freq = min_freq*20
1575 if max_freq < 2000:
1576 max_freq = 2000
1577 log_freq = True
1578 else:
1579 log_freq = False
1581 # check if all input files are results:
1582 exts = TableData.ext_formats.values()
1583 results = True
1584 # check and group by recording:
1585 result_files = []
1586 for f in sorted(files):
1587 _, base_name, _, _, ftype, _, ext = parse_filename(f)
1588 if ext == 'zip' or (ext in exts and ftype in file_types):
1589 if len(result_files) == 0 or \
1590 not result_files[-1][-1].startswith(base_name):
1591 result_files.append([f])
1592 else:
1593 result_files[-1].append(f)
1594 else:
1595 results = False
1596 break
1597 if results:
1598 files = result_files
1600 # adjust verbosity:
1601 v = verbose
1602 if len(files) > 1:
1603 v += 1
1605 # run on pool:
1606 pool_args = (results, load_kwargs, cfg, args.channel, args.time,
1607 args.time_file, args.mode, log_freq, min_freq,
1608 max_freq, args.save_data, args.zip_file,
1609 args.all_eods, spec_plots, args.skip_bad,
1610 args.save_plot, multi_pdf, args.save_subplots,
1611 args.outpath, args.keep_path, v-1, plot_level)
1612 if results:
1613 pool_args = (results, args.rawdata_path, load_kwargs,
1614 args.all_eods, spec_plots, args.skip_bad,
1615 args.save_plot, multi_pdf, args.save_subplots,
1616 log_freq, min_freq, max_freq, args.outpath,
1617 args.keep_path, v)
1618 if args.jobs is not None and (args.save_data or args.save_plot) and len(files) > 1:
1619 cpus = cpu_count() if args.jobs == 0 else args.jobs
1620 if verbose > 1:
1621 print('run on %d cpus' % cpus)
1622 p = Pool(cpus)
1623 p.map(run_thunderfish, zip(files, [pool_args]*len(files)))
1624 else:
1625 list(map(run_thunderfish, zip(files, [pool_args]*len(files))))
1626 if multi_pdf is not None:
1627 multi_pdf.close()
1630if __name__ == '__main__':
1631 freeze_support() # needed by multiprocessing for some weired windows stuff
1632 main()