Coverage for src/thunderfish/eodanalysis.py: 11%

1""" 

2Analyse EOD waveforms. 

3 

4## EOD detection 

5 

6- `detect_eods()`: detect EODs of all fish present in the data. 

7 

8## EOD waveform analysis 

9 

10- `eod_waveform()`: compute an averaged EOD waveform. 

11- `adjust_eodf()`: adjust EOD frequencies to a standard temperature. 

12 

13## Undo high-pass filter 

14 

15- `unfilter()`: apply inverse high-pass filter on time series. 

16- `unfilter_coeff()`: apply inverse high-pass filter on Fourier coefficients. 

17 

18## Similarity of EOD waveforms 

19 

20- `wave_similarity()`: root-mean squared difference between two wave fish EODs. 

21- `pulse_similarity()`: root-mean squared difference between two pulse fish EODs. 

22- `load_species_waveforms()`: load template EOD waveforms for species matching. 

23 

24## Quality assessment 

25 

26- `clipped_fraction()`: compute fraction of clipped EOD waveform snippets. 

27 

28## Visualization 

29 

30- `plot_eod_recording()`: plot a zoomed in range of the recorded trace. 

31- `zoom_eod_recording()`: set optimal time limits for EOD recording. 

32- `plot_eod_snippets()`: plot a few EOD waveform snippets. 

33 

34## Storage 

35 

36- `save_eod_waveform()`: save mean EOD waveform to file. 

37- `load_eod_waveform()`: load EOD waveform from file. 

38- `save_analysis(): save EOD analysis results to files. 

39- `load_analysis()`: load EOD analysis files. 

40- `load_recording()`: load recording. 

41 

42## Configuration 

43 

44- `add_eod_analysis_config()`: add parameters for EOD analysis functions to configuration. 

45- `eod_waveform_args()`: retrieve parameters for `eod_waveform()` from configuration. 

46- `add_species_config()`: add parameters needed for assigning EOD waveforms to species. 

47- `add_eod_quality_config()`: add parameters for `wave_quality()` and `pulse_quality()` to configuration. 

48""" 

49 

50import os 

51import io 

52import numpy as np 

53import matplotlib.pyplot as plt 

54 

55from pathlib import Path 

56from zipfile import ZipFile 

57from audioio import get_str 

58from thunderlab.eventdetection import detect_peaks, snippets 

59from thunderlab.fourier import normalize_fourier_coeffs 

60from thunderlab.powerspectrum import decibel, spectrogram, spectrum_args 

61from thunderlab.tabledata import TableData 

62from thunderlab.dataloader import DataLoader 

63 

64from .files import parse_filename, open_from_zip, close_zip 

65from .fakefish import pulsefish_spectrum 

66from .fakefish import normalize_pulsefish, export_pulsefish 

67from .fakefish import normalize_wavefish, export_wavefish 

68from .harmonics import harmonic_groups_args, psd_peak_detection_args 

69from .harmonics import harmonic_groups, closest, consistent 

70from .harmonics import plot_selected_groups 

71from .pulses import extract_pulsefish, extract_pulsefish_args 

72from .pulseanalysis import analyze_pulse, analyze_pulse_args 

73from .pulseanalysis import pulsetrain_spectrum 

74from .pulseanalysis import pulse_quality, pulse_quality_args 

75from .pulseanalysis import save_pulse_fish, load_pulse_fish 

76from .pulseanalysis import save_pulse_spectrum, load_pulse_spectrum 

77from .pulseanalysis import save_pulse_phases, load_pulse_phases 

78from .pulseanalysis import save_pulse_gaussians, load_pulse_gaussians 

79from .pulseanalysis import save_pulse_times, load_pulse_times 

80from .waveanalysis import extract_wave, analyze_wave, wave_quality 

81from .waveanalysis import extract_wave_args, analyze_wave_args 

82from .waveanalysis import wave_quality_args 

83from .waveanalysis import save_wave_eodfs, load_wave_eodfs 

84from .waveanalysis import save_wave_fish, load_wave_fish 

85from .waveanalysis import save_wave_phases, load_wave_phases 

86from .waveanalysis import save_wave_spectrum, load_wave_spectrum 

87 

88 

89def detect_eods(data, rate, power_freqs, power_times, powers, 

90 min_clip, max_clip, name, mode, verbose, plot_level, cfg): 

91 """Detect EODs of all fish present in the data. 

92 

93 This is meant to detect EODs in short (many seconds) recordings, 

94 like for example with the fishfinder or dedicated head-tail 

95 measurements. The assumption is that the EOD of every fish is 

96 stable in amplitude and discharge rate within the data. 

97 

98 Parameters 

99 ---------- 

100 data: array of floats 

101 The recording in which to detect EODs. 

102 rate: float 

103 Sampling rate of the dataset. 

104 power_freqs: 1D array of float or None 

105 Frequencies for `powers`. Can be empty. 

106 power_times: 1D array of float 

107 Times for `powers`. Can be empty. 

108 powers: 2D array of float or None 

109 Spectrogram of the data, if available. First column are 

110 frequencies, second column times. 

111 min_clip: float 

112 Minimum amplitude that is not clipped. 

113 max_clip: float 

114 Maximum amplitude that is not clipped. 

115 name: string 

116 Name of the recording (e.g. its filename). 

117 mode: string 

118 Characters in the string indicate what and how to analyze: 

119 - 'w': analyze wavefish 

120 - 'p': analyze pulsefish 

121 - 'P': analyze only the pulsefish with the largest amplitude (not implemented yet)  

122 verbose: int 

123 Print out information about EOD detection if greater than zero. 

124 plot_level : int 

125 Similar to verbosity levels, but with plots.  

126 cfg: ConfigFile 

127 Configuration parameters. 

128 

129 Returns 

130 ------- 

131 power_freqs: 1D array of float 

132 Frequencies for `powers`. Can be empty. 

133 powers: 1D array of float 

134 Power spectrum of the data. 

135 wave_eodfs: list of 2D arrays 

136 Frequency and power of fundamental frequency/harmonics of all wave fish. 

137 wave_indices: array of int 

138 Indices of wave fish mapping from wave_eodfs to eod_props. 

139 If negative, then that EOD frequency has no waveform described in eod_props. 

140 eod_props: list of dict 

141 Lists of EOD properties as returned by analyze_pulse() and analyze_wave() 

142 for each waveform in mean_eods. 

143 mean_eods: list of 2-D arrays with time, mean, sem, and fit. 

144 Averaged EOD waveforms of pulse and wave fish. 

145 spec_data: list of 2_D arrays 

146 For each pulsefish a power spectrum of the single pulse and for 

147 each wavefish the relative amplitudes and phases of the harmonics. 

148 phase_data: list of dict 

149 For each pulse fish a dictionary with phase properties 

150 (indices, times, amplitudes, relamplitudes, widths, areas, relareas, zeros), 

151 empty dict for wave fish. 

152 pulse_data: list of dict 

153 For each pulse fish a dictionary with phase times, amplitudes and standard 

154 deviations of Gaussians fitted to the pulse waveform. Use the 

155 functions provided in thunderfish.fakefish to simulate pulse 

156 fish EODs from this data. 

157 power_thresh: 2 D array or None 

158 Frequency (first column) and power (second column) of threshold 

159 derived from single pulse spectra to discard false wave fish. 

160 None if no pulse fish was detected. 

161 skip_reason: list of string 

162 Reasons, why an EOD was discarded. 

163 

164 """ 

165 dfreq = np.nan 

166 nfft = 0 

167 wave_eodfs = [] 

168 wave_indices = [] 

169 if 'w' in mode: 

170 # detect wave fish: 

171 if power_freqs is None or power_times is None or powers is None or \ 

172 power_freqs[1] > cfg.value('frequencyResolution'): 

173 power_freqs, power_times, powers = \ 

174 spectrogram(data, rate, **spectrum_args(cfg)) 

175 dfreq = np.mean(np.diff(power_freqs)) 

176 nfft = int(rate/dfreq) 

177 h_kwargs = psd_peak_detection_args(cfg) 

178 h_kwargs.update(harmonic_groups_args(cfg)) 

179 wave_eodfs_list = [] 

180 for i, psd in enumerate(powers.T): 

181 wave_eodfs = harmonic_groups(power_freqs, psd, verbose - 1, 

182 **h_kwargs)[0] 

183 if verbose > 0 and powers.shape[1] > 1: 

184 print(f'{len(wave_eodfs)} fundamental frequencies detected in spectrum of window {i}:') 

185 if len(wave_eodfs) > 0: 

186 print(' ' + ' '.join([f'{freq[0, 0]:.1f}' for freq in wave_eodfs])) 

187 else: 

188 print(' none') 

189 wave_eodfs_list.append(wave_eodfs) 

190 min_closest = max(1, (len(wave_eodfs_list) + 1) // 3) 

191 freq_thresh = cfg.value('frequencyThreshold') 

192 wave_eodfs, wave_windows = \ 

193 closest(wave_eodfs_list, 

194 df_thresh=freq_thresh, close_thresh=1*freq_thresh, 

195 min_closest=min_closest) 

196 if len(wave_windows) > 0: 

197 p0 = np.min(wave_windows[:, 0]) 

198 p1 = np.max(wave_windows[:, 1]) 

199 powers = np.mean(powers[:, p0:p1 + 1], 1) 

200 else: 

201 powers = np.mean(powers, 1) 

202 if verbose > 0: 

203 if verbose > 1: 

204 print() 

205 if len(wave_eodfs) > 0: 

206 fstr = 'ies' if len(wave_eodfs) > 1 else 'y' 

207 print(f'found {len(wave_eodfs):2d} EOD frequenc{fstr} sufficiently close in all spectra:') 

208 if verbose > 1: 

209 for freq, win in zip(wave_eodfs, wave_windows): 

210 print(f' {freq[0, 0]:6.1f}Hz in spectra {win[0]} - {win[1]}') 

211 else: 

212 print(' ' + ' '.join([f'{freq[0, 0]:.1f}' for freq in wave_eodfs])) 

213 else: 

214 print('no fundamental frequencies are sufficiently close in all spectra') 

215 if plot_level > 0: 

216 fig, ax = plt.subplots(layout='constrained') 

217 ax.set_title(f'{len(wave_eodfs)} closest EOD frequencies') 

218 plot_selected_groups(ax, wave_eodfs_list, wave_eodfs, wave_windows) 

219 ax.set_xlabel('index of spectrum segment') 

220 plt.show() 

221 

222 # analysis results: 

223 eod_props = [] 

224 mean_eods = [] 

225 spec_data = [] 

226 phase_data = [] 

227 pulseeod_data = [] 

228 power_thresh = None 

229 skip_reason = [] 

230 max_pulse_amplitude = 0.0 

231 

232 if 'p' in mode: 

233 # detect pulse fish: 

234 eods, eod_times, eod_peaktimes, _ = \ 

235 extract_pulsefish(data, rate, 

236 verbose=verbose - 1, 

237 plot_level=plot_level, 

238 **extract_pulsefish_args(cfg)) 

239 if verbose > 0: 

240 if len(eod_times) > 0: 

241 print(f'found {len(eod_times):2d} pulsefish EODs') 

242 else: 

243 print('no pulsefish EODs found') 

244 if plot_level > 0 and len(eod_times) > 0: 

245 fig, ax = plt.subplots() 

246 for eod, times in zip(eods, eod_times): 

247 ax.plot(eod[:, 0], eod[:, 1], label=f'{1/np.median(np.diff(times)):.1f}Hz') 

248 ax.legend() 

249 plt.show() 

250 

251 # analyse eod waveform of pulse-fish: 

252 for mean_eod, eod_ts, eod_pts in zip(eods, eod_times, eod_peaktimes): 

253 """ 

254 mean_eod, eod_times0 = \ 

255 eod_waveform(data, rate, eod_ts, win_fac=0.8, 

256 min_win=cfg.value('eodMinPulseSnippet'), 

257 min_sem=False, **eod_waveform_args(cfg)) 

258 """ 

259 unfilter_cutoff = cfg.value('unfilterCutoff') 

260 if unfilter_cutoff and unfilter_cutoff > 0: 

261 unfilter(mean_eod[:, 1], frate, unfilter_cutoff) 

262 mean_eod, props, phases, pulse, power = \ 

263 analyze_pulse(mean_eod, None, eod_ts, verbose=verbose-1, 

264 **analyze_pulse_args(cfg)) 

265 if len(phases) == 0: 

266 if verbose > 0: 

267 print('no phases in pulse EOD detected') 

268 continue 

269 clipped_frac = clipped_fraction(data, rate, eod_ts, 

270 mean_eod, min_clip, max_clip) 

271 props['peaktimes'] = eod_ts # XXX that should go into analyze pulse 

272 props['index'] = len(eod_props) 

273 props['clipped'] = clipped_frac 

274 props['samplerate'] = rate 

275 props['nfft'] = nfft 

276 props['dfreq'] = dfreq 

277 

278 # add good waveforms only: 

279 skips, msg, skipped_clipped = \ 

280 pulse_quality(props, **pulse_quality_args(cfg)) 

281 

282 if len(skips) == 0: 

283 eod_props.append(props) 

284 mean_eods.append(mean_eod) 

285 spec_data.append(power) 

286 phase_data.append(phases) 

287 pulseeod_data.append(pulse) 

288 if verbose > 0: 

289 print(f'take {props['EODf']:7.2f}Hz pulse fish: {msg}') 

290 else: 

291 skip_reason += [f'{props['EODf']:.2f}Hz pulse fish {skips}'] 

292 if verbose > 0: 

293 print(f'skip {props['EODf']:7.2f}Hz pulse fish: {skips} ({msg})') 

294 

295 # threshold for wave fish peaks based on single pulse spectra: 

296 if len(skips) == 0 or skipped_clipped: 

297 if max_pulse_amplitude < props['ppampl']: 

298 max_pulse_amplitude = props['ppampl'] 

299 

300 pulse_freqs, pulse_power = \ 

301 pulsetrain_spectrum(eod_pts, mean_eod, None, len(data)/rate, rate, 

302 fade_frac=0.05, **spectrum_args(cfg)) 

303 pulse_power *= len(data)/rate/props['period']/len(props['peaktimes']) 

304 pulse_power *= 5 

305 if power_thresh is None: 

306 power_thresh = np.zeros((len(pulse_freqs), 2)) 

307 power_thresh[:, 0] = pulse_freqs 

308 power_thresh[:, 1] = pulse_power 

309 else: 

310 power_thresh[:, 1] += pulse_power 

311 

312 # remove wavefish below pulse fish power: 

313 if 'w' in mode and power_thresh is not None: 

314 n = len(wave_eodfs) 

315 maxh = 3 # XXX TODO make parameter 

316 df = power_thresh[1,0] - power_thresh[0,0] 

317 for k, fish in enumerate(reversed(wave_eodfs)): 

318 idx = np.array(fish[:maxh,0]//df, dtype=int) 

319 for offs in range(-2, 3): 

320 nbelow = np.sum(fish[:maxh,1] < power_thresh[idx+offs,1]) 

321 if nbelow > 0: 

322 wave_eodfs.pop(n-1-k) 

323 if verbose > 0: 

324 print(f'skip {fish[0,0]:7.2f}Hz wave fish: {nbelow:2d} harmonics are below pulsefish threshold') 

325 break 

326 

327 if 'w' in mode: 

328 # analyse EOD waveform of all wavefish: 

329 if verbose > 1: 

330 print() 

331 fish_powers = np.array([np.sum(fish[:,1]) for fish in wave_eodfs]) 

332 power_indices = np.argsort(-fish_powers) 

333 wave_indices = np.zeros(len(wave_eodfs), dtype=int) - 3 

334 for k, idx in enumerate(power_indices): 

335 fish = wave_eodfs[idx] 

336 window = wave_windows[idx] 

337 iw = int(rate/power_freqs[1])//2 

338 i0 = int(power_times[window[0]]*rate) - iw 

339 i1 = int(power_times[window[1]]*rate) + iw 

340 coeffs, mean_eod, eod_freq, times, n_eods, eskips = \ 

341 extract_wave(data[i0:i1], rate, fish[0, 0], 

342 power_freqs[1], verbose=verbose - 1, 

343 plot_level=plot_level, 

344 **extract_wave_args(cfg)) 

345 if len(mean_eod) == 0 or len(eskips) > 0: 

346 # waveform extraction failed 

347 continue 

348 for h in range(len(fish)): 

349 fish[h, 0] = (h + 1)*eod_freq 

350 unfilter_cutoff = cfg.value('unfilterCutoff') 

351 if unfilter_cutoff: 

352 coeffs = unfilter_coeff(eod_freq, coeffs, unfilter_cutoff) 

353 if plot_level > 0: 

354 w = fourier_synthesis(eod_freq, coeffs, 

355 1/mean_eod[1, 0], len(mean_eod)) 

356 fig, ax = plt.subplots(layout='constrained') 

357 ax.set_title(f'EODf={eod_freq:.1f}Hz, unfilter high-pass filter $f_{{cutoff}}={unfilter_cutoff:.0f}$Hz') 

358 ax.plot(1000*mean_eod[:, 0], mean_eod[:, 1], 

359 label='original') 

360 ax.plot(1000*mean_eod[:, 0], w, label='unfiltered') 

361 ax.set_xlabel('time [ms]') 

362 ax.legend() 

363 plt.show() 

364 mean_eod, props, phases, sdata = \ 

365 analyze_wave(mean_eod, None, fish, coeffs, 

366 **analyze_wave_args(cfg)) 

367 eod_times = np.arange(i0/rate, i1/rate, 1/fish[0, 0]) 

368 clipped_frac = clipped_fraction(data[i0:i1], rate, eod_times, 

369 mean_eod, min_clip, max_clip) 

370 props['n'] = n_eods 

371 props['nsegments'] = len(times) 

372 props['index'] = len(eod_props) 

373 props['clipped'] = clipped_frac 

374 props['samplerate'] = rate 

375 props['nfft'] = nfft 

376 props['dfreq'] = dfreq 

377 # remove wave fish that are smaller than the largest pulse fish: 

378 if props['ppampl'] < 0.01*max_pulse_amplitude: 

379 rm_indices = power_indices[k:] 

380 if verbose > 0: 

381 print(f'skip {props['EODf']:7.2f}Hz wave fish: power={decibel(fish_powers[idx]):5.1f}dB, p-p amplitude={decibel(props['ppampl']):5.1f}dB smaller than pulse fish={decibel(max_pulse_amplitude):5.1f}dB - 20dB') 

382 for idx in rm_indices[1:]: 

383 print(f'skip {wave_eodfs[idx][0,0]:7.2f}Hz wave fish: power={decibel(fish_powers[idx]):5.1f}dB even smaller') 

384 if verbose > 1: 

385 print() 

386 wave_eodfs = [eodfs for idx, eodfs in enumerate(wave_eodfs) 

387 if idx not in rm_indices] 

388 wave_indices = np.array([idcs for idx, idcs in enumerate(wave_indices) 

389 if idx not in rm_indices], dtype=int) 

390 break 

391 # add good waveforms only: 

392 remove, skips, msg = wave_quality(props, sdata[1:, 3], 

393 **wave_quality_args(cfg)) 

394 if len(eskips) + len(skips) == 0: 

395 wave_indices[idx] = props['index'] 

396 eod_props.append(props) 

397 mean_eods.append(mean_eod) 

398 spec_data.append(sdata) 

399 phase_data.append(phases) 

400 pulseeod_data.append(dict()) 

401 if verbose > 0: 

402 print(f'take {props['EODf']:7.2f}Hz wave fish:', msg) 

403 else: 

404 if len(skips) == 0: 

405 skips = eskips 

406 elif len(eskips) > 0: 

407 skips = eskips + ', ' + skips 

408 wave_indices[idx] = -2 if remove else -1 

409 skip_reason += [f'{props["EODf"]:.1f}Hz wave fish {skips}'] 

410 if verbose > 0: 

411 rstr = 'remove' if remove else 'skip' 

412 print(f'{rstr:<7s} {props["EODf"]:7.2f}Hz wave fish: {skips} ({msg})') 

413 if verbose > 1: 

414 print() 

415 wave_eodfs = [eodfs for idx, eodfs in zip(wave_indices, wave_eodfs) if idx > -2] 

416 wave_indices = np.array([idx for idx in wave_indices if idx > -2], dtype=int) 

417 return (power_freqs, powers, wave_eodfs, wave_indices, eod_props, 

418 mean_eods, spec_data, phase_data, pulseeod_data, 

419 power_thresh, skip_reason) 

420 

421 

422def eod_waveform(data, rate, eod_times, win_fac=2.0, min_win=0.01, 

423 min_sem=False, max_eods=None): 

424 """Extract data snippets around each EOD, and compute a mean waveform with standard error. 

425 

426 Retrieving the EOD waveform of a wave fish works under the following 

427 conditions: (i) at a signal-to-noise ratio \\(SNR = P_s/P_n\\), 

428 i.e. the power \\(P_s\\) of the EOD of interest relative to the 

429 largest other EOD \\(P_n\\), we need to average over at least \\(n > 

430 (SNR \\cdot c_s^2)^{-1}\\) snippets to bring the standard error of the 

431 averaged EOD waveform down to \\(c_s\\) relative to its 

432 amplitude. For a s.e.m. less than 5% ( \\(c_s=0.05\\) ) and an SNR of 

433 -10dB (the signal is 10 times smaller than the noise, \\(SNR=0.1\\) ) we 

434 get \\(n > 0.00025^{-1} = 4000\\) data snippets - a recording a 

435 couple of seconds long. (ii) Very important for wave fish is that 

436 they keep their frequency constant. Slight changes in the EOD 

437 frequency will corrupt the average waveform. If the period of the 

438 waveform changes by \\(c_f=\\Delta T/T\\), then after \\(n = 

439 1/c_f\\) periods moved the modified waveform through a whole period. 

440 This is in the range of hundreds or thousands waveforms. 

441 

442 NOTE: we need to take into account a possible error in the estimate 

443 of the EOD period. This will limit the maximum number of snippets to 

444 be averaged. 

445 

446 If `min_sem` is set, the algorithm checks for a global minimum of 

447 the s.e.m. as a function of snippet number. If there is one then 

448 the average is computed for this number of snippets, otherwise all 

449 snippets are taken from the provided data segment. Note that this 

450 check only works for the strongest EOD in a recording. For weaker 

451 EOD the s.e.m. always decays with snippet number (empirical 

452 observation). 

453 

454 TODO: use power spectra to check for changes in EOD frequency! 

455 

456 Parameters 

457 ---------- 

458 data: 1-D array of float 

459 The data to be analysed. 

460 rate: float 

461 Sampling rate of the data in Hertz. 

462 eod_times: 1-D array of float 

463 Array of EOD times in seconds over which the waveform should be 

464 averaged. 

465 WARNING: The first data point must be at time zero! 

466 win_fac: float 

467 The snippet size is the EOD period times `win_fac`. The EOD period 

468 is determined as the minimum interval between EOD times. 

469 min_win: float 

470 The minimum size of the snippets in seconds. 

471 min_sem: bool 

472 If set, check for minimum in s.e.m. to set the maximum numbers 

473 of EODs to be used for computing the average waveform. 

474 max_eods: int or None 

475 Maximum number of EODs to be used for averaging. 

476 unfilter_cutoff: float 

477 If not zero, the cutoff frequency for an inverse high-pass filter 

478 applied to the mean EOD waveform. 

479  

480 Returns 

481 ------- 

482 mean_eod: 2-D array 

483 Average of the EOD snippets. First column is time in seconds, 

484 second column the mean eod, third column the standard error. 

485 eod_times: 1-D array 

486 Times of EOD peaks in seconds that have been actually used to calculate the 

487 averaged EOD waveform. 

488 """ 

489 # indices of EOD times: 

490 eod_idx = np.round(eod_times*rate).astype(int) 

491 

492 # window size: 

493 period = np.min(np.diff(eod_times)) 

494 win = 0.5*win_fac*period 

495 if 2*win < min_win: 

496 win = 0.5*min_win 

497 win_inx = int(win*rate) 

498 

499 # extract snippets: 

500 eod_times = eod_times[(eod_idx >= win_inx) & (eod_idx < len(data)-win_inx)] 

501 eod_idx = eod_idx[(eod_idx >= win_inx) & (eod_idx < len(data)-win_inx)] 

502 if max_eods and max_eods > 0 and len(eod_idx) > max_eods: 

503 dn = (len(eod_idx) - max_eods)//2 

504 eod_times = eod_times[dn:dn+max_eods] 

505 eod_idx = eod_idx[dn:dn+max_eods] 

506 eod_snippets = snippets(data, eod_idx, -win_inx, win_inx) 

507 if len(eod_snippets) == 0: 

508 return np.zeros((0, 3)), eod_times 

509 

510 # optimal number of snippets: 

511 step = 10 

512 if min_sem and len(eod_snippets) > step: 

513 sems = [np.mean(np.std(eod_snippets[:k], axis=0, ddof=1)/np.sqrt(k)) 

514 for k in range(step, len(eod_snippets), step)] 

515 idx = np.argmin(sems) 

516 # there is a local minimum: 

517 if idx > 0 and idx < len(sems)-1: 

518 maxn = step*(idx+1) 

519 eod_snippets = eod_snippets[:maxn] 

520 eod_times = eod_times[:maxn] 

521 

522 # mean and std of snippets: 

523 mean_eod = np.zeros((len(eod_snippets[0]), 3)) 

524 mean_eod[:, 1] = np.mean(eod_snippets, axis=0) 

525 if len(eod_snippets) > 1: 

526 mean_eod[:, 2] = np.std(eod_snippets, axis=0, ddof=1)/np.sqrt(len(eod_snippets)) 

527 

528 # time axis: 

529 mean_eod[:, 0] = (np.arange(len(mean_eod)) - win_inx) / rate 

530 

531 return mean_eod, eod_times 

532 

533 

534def adjust_eodf(eodf, temp, temp_adjust=25.0, q10=1.62): 

535 """Adjust EOD frequencies to a standard temperature using Q10. 

536 

537 Parameters 

538 ---------- 

539 eodf: float or ndarray 

540 EOD frequencies. 

541 temp: float 

542 Temperature in degree celsisus at which EOD frequencies in 

543 `eodf` were measured. 

544 temp_adjust: float 

545 Standard temperature in degree celsisus to which EOD 

546 frequencies are adjusted. 

547 q10: float 

548 Q10 value describing temperature dependence of EOD 

549 frequencies. The default of 1.62 is from Dunlap, Smith, Yetka 

550 (2000) Brain Behav Evol, measured for Apteronotus 

551 lepthorhynchus in the lab. 

552 

553 Returns 

554 ------- 

555 eodf_corrected: float or array 

556 EOD frequencies adjusted to `temp_adjust` using `q10`. 

557 """ 

558 return eodf*q10**((temp_adjust - temp) / 10.0) 

559 

560 

561def unfilter(data, rate, fcutoff): 

562 """Apply inverse high-pass filter on time series. 

563 

564 Assumes high-pass filter 

565 \\[ \\tau \\dot y = -y + \\tau \\dot x \\] 

566 has been applied on the original data \\(x\\), where 

567 \\(\\tau=(2\\pi f_{cutoff})^{-1}\\) is the time constant of the 

568 filter. To recover \\(x\\) the ODE 

569 \\[ \\tau \\dot x = y + \\tau \\dot y \\] 

570 is applied on the filtered data \\(y\\). 

571 

572 Parameters 

573 ---------- 

574 data: 1D ndarray of float 

575 High-pass filtered original data. 

576 rate: float 

577 Sampling rate of `data` in Hertz. 

578 fcutoff: float 

579 Cutoff frequency \\(f_{cutoff}\\) of the high-pass filter in Hertz. 

580 

581 Returns 

582 ------- 

583 data: ndarray 

584 Recovered original data. 

585 """ 

586 tau = 0.5/np.pi/fcutoff 

587 fac = tau*rate 

588 data -= np.mean(data) 

589 d0 = data[0] 

590 x = d0 

591 for k in range(len(data)): 

592 d1 = data[k] 

593 x += (d1 - d0) + d0/fac 

594 data[k] = x 

595 d0 = d1 

596 return data 

597 

598 

599def unfilter_coeff(freq, coeffs, fcutoff): 

600 """Apply inverse high-pass filter on Fourier coefficients. 

601 

602 Assumes first-order high-pass filter 

603 \\[ \\tau \\dot y = -y + \\tau \\dot x \\] 

604 has been applied on the original data \\(x\\), where 

605 \\(\\tau=(2\\pi f_{cutoff})^{-1}\\) is the time constant of the 

606 filter. The transfer function of this high-pass filter is 

607 \\[ H(\\omega) = \\frac{i\\omega\\tau}{1 + i\\omega\\tau} \\] 

608 

609 To undo the effect of the high-pass filter, the complex-valued 

610 Fourier coefficients \\(c_k\\) of the waveform for frequencies \\(k f_0\\) 

611 are multiplied with the inverse transfer function 

612 \\[ IH(\\omega) = \\frac{1}{H(\\omega)} = 1 - i\\frac{1}{\\omega\\tau} \\] 

613 which simplifies for the harmonics to 

614 \\[ IH(kf_0) = 1 - i\\frac{f_{cutoff}}{k f_0} \\] 

615  

616 

617 Parameters 

618 ---------- 

619 freq: float 

620 Fundamental frequency \\(f_0\\). 

621 coeffs: 1D array of complex 

622 For each harmonics the complex valued Fourier coefficient 

623 as, for example, returned by `fourier_coeffs()`. 

624 The first one is the offset. 

625 fcutoff: float 

626 Cutoff frequency \\(f_{cutoff}\\) of the high-pass filter in Hertz. 

627 

628 Returns 

629 ------- 

630 coeffs: 1D array of complex 

631 Unfiltered Fourier coefficients. Normalized to zero phase shift 

632 of the fundamental and no offset. 

633 """ 

634 ihp_coeffs = np.zeros(len(coeffs), dtype=complex) 

635 ihp_coeffs[1:] = 1 - 1j/(np.arange(1, len(coeffs))*freq/fcutoff) 

636 coeffs *= ihp_coeffs 

637 coeffs = normalize_fourier_coeffs(coeffs) 

638 return coeffs 

639 

640 

641def load_species_waveforms(species_file='none'): 

642 """Load template EOD waveforms for species matching. 

643  

644 Parameters 

645 ---------- 

646 species_file: str 

647 Name of file containing species definitions. The content of 

648 this file is as follows: 

649  

650 - Empty lines and line starting with a hash ('#') are skipped. 

651 - A line with the key-word 'wavefish' marks the beginning of the 

652 table for wave fish. 

653 - A line with the key-word 'pulsefish' marks the beginning of the 

654 table for pulse fish. 

655 - Each line in a species table has three fields, 

656 separated by colons (':'): 

657  

658 1. First field is an abbreviation of the species name. 

659 2. Second field is the filename of the recording containing the 

660 EOD waveform. 

661 3. The optional third field is the EOD frequency of the EOD waveform. 

662 

663 The EOD frequency is used to normalize the time axis of a 

664 wave fish EOD to one EOD period. If it is not specified in 

665 the third field, it is taken from the corresponding 

666 *-wavespectrum-* file, if present. Otherwise the species is 

667 excluded and a warning is issued. 

668 

669 Example file content: 

670 ``` plain 

671 Wavefish 

672 Aptero : F_91009L20-eodwaveform-0.csv : 823Hz 

673 Eigen : C_91008L01-eodwaveform-0.csv 

674 

675 Pulsefish 

676 Gymnotus : pulsefish/gymnotus.csv 

677 Brachy : H_91009L11-eodwaveform-0.csv 

678 ``` 

679  

680 Returns 

681 ------- 

682 wave_names: list of str 

683 List of species names of wave-type fish. 

684 wave_eods: list of 2-D arrays 

685 List of EOD waveforms of wave-type fish corresponding to 

686 `wave_names`. First column of a waveform is time in seconds, 

687 second column is the EOD waveform. The waveforms contain 

688 exactly one EOD period. 

689 pulse_names: list of str 

690 List of species names of pulse-type fish. 

691 pulse_eods: list of 2-D arrays 

692 List of EOD waveforms of pulse-type fish corresponding to 

693 `pulse_names`. First column of a waveform is time in seconds, 

694 second column is the EOD waveform. 

695 """ 

696 if not Path(species_file).is_file(): 

697 return [], [], [], [] 

698 wave_names = [] 

699 wave_eods = [] 

700 pulse_names = [] 

701 pulse_eods = [] 

702 fish_type = 'wave' 

703 with open(species_file, 'r') as sf: 

704 for line in sf: 

705 line = line.strip() 

706 if len(line) == 0 or line[0] == '#': 

707 continue 

708 if line.lower() == 'wavefish': 

709 fish_type = 'wave' 

710 elif line.lower() == 'pulsefish': 

711 fish_type = 'pulse' 

712 else: 

713 ls = [s.strip() for s in line.split(':')] 

714 if len(ls) < 2: 

715 continue 

716 name = ls[0] 

717 waveform_file = ls[1] 

718 eod = TableData(waveform_file).array() 

719 eod[:, 0] *= 0.001 

720 if fish_type == 'wave': 

721 eodf = None 

722 if len(ls) > 2: 

723 eodf = float(ls[2].replace('Hz', '').strip()) 

724 else: 

725 spectrum_file = waveform_file.replace('eodwaveform', 'wavespectrum') 

726 try: 

727 wave_spec = TableData(spectrum_file) 

728 eodf = wave_spec[0, 1] 

729 except FileNotFoundError: 

730 pass 

731 if eodf is None: 

732 print('warning: unknown EOD frequency of %s. Skip.' % name) 

733 continue 

734 eod[:, 0] *= eodf 

735 wave_names.append(name) 

736 wave_eods.append(eod[:, :2]) 

737 elif fish_type == 'pulse': 

738 pulse_names.append(name) 

739 pulse_eods.append(eod[:, :2]) 

740 return wave_names, wave_eods, pulse_names, pulse_eods 

741 

742 

743def wave_similarity(eod1, eod2, eod1f=1.0, eod2f=1.0): 

744 """Root-mean squared difference between two wave fish EODs. 

745 

746 Compute the root-mean squared difference between two wave fish 

747 EODs over one period. The better sampled signal is down-sampled to 

748 the worse sampled one. Amplitude are normalized to peak-to-peak 

749 amplitude before computing rms difference. Also compute the rms 

750 difference between the one EOD and the other one inverted in 

751 amplitude. The smaller of the two rms values is returned. 

752 

753 Parameters 

754 ---------- 

755 eod1: 2-D array 

756 Time and amplitude of reference EOD. 

757 eod2: 2-D array 

758 Time and amplitude of EOD that is to be compared to `eod1`. 

759 eod1f: float 

760 EOD frequency of `eod1` used to transform the time axis of `eod1` 

761 to multiples of the EOD period. If already normalized to EOD period, 

762 as for example by the `load_species_waveforms()` function, then 

763 set the EOD frequency to one (default). 

764 eod2f: float 

765 EOD frequency of `eod2` used to transform the time axis of `eod2` 

766 to multiples of the EOD period. If already normalized to EOD period, 

767 as for example by the `load_species_waveforms()` function, then 

768 set the EOD frequency to one (default). 

769 

770 Returns 

771 ------- 

772 rmse: float 

773 Root-mean-squared difference between the two EOD waveforms relative to 

774 their standard deviation over one period. 

775 """ 

776 # copy: 

777 eod1 = np.array(eod1[:, :2]) 

778 eod2 = np.array(eod2[:, :2]) 

779 # scale to multiples of EOD period: 

780 eod1[:, 0] *= eod1f 

781 eod2[:, 0] *= eod2f 

782 # make eod1 the waveform with less samples per period: 

783 n1 = int(1.0/(eod1[1,0]-eod1[0,0])) 

784 n2 = int(1.0/(eod2[1,0]-eod2[0,0])) 

785 if n1 > n2: 

786 eod1, eod2 = eod2, eod1 

787 n1, n2 = n2, n1 

788 # one period around time zero: 

789 i0 = np.argmin(np.abs(eod1[:, 0])) 

790 k0 = i0-n1//2 

791 if k0 < 0: 

792 k0 = 0 

793 k1 = k0 + n1 + 1 

794 if k1 >= len(eod1): 

795 k1 = len(eod1) 

796 # cut out one period from the reference EOD around maximum: 

797 i = k0 + np.argmax(eod1[k0:k1,1]) 

798 k0 = i-n1//2 

799 if k0 < 0: 

800 k0 = 0 

801 k1 = k0 + n1 + 1 

802 if k1 >= len(eod1): 

803 k1 = len(eod1) 

804 eod1 = eod1[k0:k1,:] 

805 # normalize amplitudes of first EOD: 

806 eod1[:, 1] -= np.min(eod1[:, 1]) 

807 eod1[:, 1] /= np.max(eod1[:, 1]) 

808 sigma = np.std(eod1[:, 1]) 

809 # set time zero to maximum of second EOD: 

810 i0 = np.argmin(np.abs(eod2[:, 0]))