Coverage for src / thunderfish / waveanalysis.py: 34%

1""" 

2Analysis of wave-type EOD waveforms. 

3 

4## Analysis of wave-type EODs 

5 

6- `waveeod_waveform()`: retrieve average EOD waveform via Fourier transform. 

7- `analyze_wave()`: analyze the EOD waveform of a wave fish. 

8 

9## Visualization 

10 

11- `plot_wave_spectrum()`: plot and annotate spectrum of wave EODs. 

12 

13## Storage 

14 

15- `save_wave_eodfs()`: save frequencies of wave EODs to file. 

16- `load_wave_eodfs()`: load frequencies of wave EODs from file. 

17- `save_wave_fish()`: save properties of wave EODs to file. 

18- `load_wave_fish()`: load properties of wave EODs from file. 

19- `save_wave_spectrum()`: save amplitude and phase spectrum of wave EOD to file. 

20- `load_wave_spectrum()`: load amplitude and phase spectrum of wave EOD from file. 

21 

22## Fit functions 

23 

24- `fourier_series()`: Fourier series of sine waves with amplitudes and phases. 

25- `exp_decay()`: exponential decay. 

26""" 

27 

28import numpy as np 

29 

30try: 

31 from matplotlib.ticker import MultipleLocator 

32 from matplotlib.artist import setp 

33except ImportError: 

34 pass 

35 

36from pathlib import Path 

37from scipy.optimize import curve_fit 

38from numba import jit 

39from thunderlab.eventdetection import threshold_crossings, threshold_crossing_times, merge_events 

40from thunderlab.powerspectrum import decibel 

41from thunderlab.tabledata import TableData 

42 

43from .harmonics import fundamental_freqs_and_power 

44 

45 

46def waveeod_waveform(data, rate, freq, win_fac=2.0): 

47 """Retrieve average EOD waveform via Fourier transform. 

48 

49 TODO: use power spectra to check minimum data segment needed and 

50 check for changes in frequency over several segments! 

51 

52 TODO: return waveform with higher samplign rate? (important for 

53 2kHz waves on 24kHz sampling). But this seems to render some EODs 

54 inacceptable in the further thunderfish processing pipeline. 

55 

56 Parameters 

57 ---------- 

58 data: 1-D array of float 

59 The data to be analysed. 

60 rate: float 

61 Sampling rate of the data in Hertz. 

62 freq: float 

63 EOD frequency. 

64 win_fac: float 

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

66 is determined as the minimum interval between EOD times. 

67  

68 Returns 

69 ------- 

70 mean_eod: 2-D array 

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

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

73 eod_times: 1-D array 

74 Times of EOD peaks in seconds that have been actually used to 

75 calculate the averaged EOD waveform. 

76 

77 """ 

78 

79 @jit(nopython=True) 

80 def fourier_wave(data, rate, freq): 

81 """ 

82 extracting wave via fourier coefficients 

83 """ 

84 twave = np.arange(0, (1+win_fac)/freq, 1/rate) 

85 wave = np.zeros(len(twave)) 

86 t = np.arange(len(data))/rate 

87 for k in range(0, 31): 

88 Xk = np.trapz(data*np.exp(-1j*2*np.pi*k*freq*t), t)*2/t[-1] 

89 wave += np.real(Xk*np.exp(1j*2*np.pi*k*freq*twave)) 

90 return wave 

91 

92 @jit(nopython=True) 

93 def fourier_range(data, rate, f0, f1, df): 

94 wave = np.zeros(1) 

95 freq = f0 

96 for f in np.arange(f0, f1, df): 

97 w = fourier_wave(data, rate, f) 

98 if np.max(w) - np.min(w) > np.max(wave) - np.min(wave): 

99 wave = w 

100 freq = f 

101 return wave, freq 

102 

103 # TODO: parameterize! 

104 tsnippet = 2 

105 min_corr = 0.98 

106 min_ampl_frac = 0.5 

107 frange = 0.1 

108 fstep = 0.1 

109 waves = [] 

110 freqs = [] 

111 times = [] 

112 step = int(tsnippet*rate) 

113 for i in range(0, len(data) - step//2, step//2): 

114 w, f = fourier_range(data[i:i + step], rate, freq - frange, 

115 freq + frange + fstep/2, fstep) 

116 waves.append(w) 

117 freqs.append(f) 

118 """ 

119 waves.append(np.zeros(1)) 

120 freqs.append(freq) 

121 for f in np.arange(freq - frange, freq + frange + fstep/2, fstep): 

122 w = fourier_wave(data[i:i + step], rate, f) 

123 if np.max(w) - np.min(w) > np.max(waves[-1]) - np.min(waves[-1]): 

124 waves[-1] = w 

125 freqs[-1] = f 

126 """ 

127 times.append(np.arange(i/rate, (i + step)/rate, 1/freqs[-1])) 

128 eod_freq = np.mean(freqs) 

129 mean_eod = np.zeros((0, 3)) 

130 eod_times = np.zeros((0)) 

131 if len(waves) == 0: 

132 return mean_eod, eod_times 

133 for k in range(len(waves)): 

134 period = int(np.ceil(rate/freqs[k])) 

135 i = np.argmax(waves[k][:period]) 

136 waves[k] = waves[k][i:] 

137 n = np.min([len(w) for w in waves]) 

138 waves = np.array([w[:n] for w in waves]) 

139 # only snippets that are similar: 

140 if len(waves) > 1: 

141 corr = np.corrcoef(waves) 

142 nmax = np.argmax(np.sum(corr > min_corr, axis=1)) 

143 if nmax <= 1: 

144 nmax = 2 

145 select = np.sum(corr > min_corr, axis=1) >= nmax 

146 waves = waves[select] 

147 times = [times[k] for k in range(len(times)) if select[k]] 

148 if len(waves) == 0: 

149 return mean_eod, eod_times 

150 # only the largest snippets: 

151 ampls = np.std(waves, axis=1) 

152 select = ampls >= min_ampl_frac*np.max(ampls) 

153 waves = waves[select] 

154 times = [times[k] for k in range(len(times)) if select[k]] 

155 if len(waves) == 0: 

156 return mean_eod, eod_times 

157 """ 

158 #plt.plot(freqs) 

159 plt.plot(waves.T) 

160 plt.show() 

161 """ 

162 mean_eod = np.zeros((n, 3)) 

163 mean_eod[:, 0] = np.arange(len(mean_eod))/rate 

164 mean_eod[:, 1] = np.mean(waves, axis=0) 

165 mean_eod[:, 2] = np.std(waves, axis=0) 

166 eod_times = np.concatenate(times) 

167 return mean_eod, eod_times 

168 

169 

170def fourier_series(t, freq, *ap): 

171 """Fourier series of sine waves with amplitudes and phases. 

172 

173 x(t) = sum_{i=0}^n ap[2*i]*sin(2 pi (i+1) freq t + ap[2*i+1]) 

174  

175 Parameters 

176 ---------- 

177 t: float or array 

178 Time. 

179 freq: float 

180 Fundamental frequency. 

181 *ap: list of floats 

182 The amplitudes and phases (in rad) of the fundamental and harmonics. 

183  

184 Returns 

185 ------- 

186 x: float or array 

187 The Fourier series evaluated at times `t`. 

188 """ 

189 omega = 2.0*np.pi*freq 

190 x = 0.0 

191 for i, (a, p) in enumerate(zip(ap[0:-1:2], ap[1::2])): 

192 x += a*np.sin((i+1)*omega*t+p) 

193 return x 

194 

195 

196def analyze_wave(eod, ratetime, freq, n_harm=10, power_n_harmonics=0, 

197 n_harmonics=3, flip_wave='none'): 

198 """Analyze the EOD waveform of a wave fish. 

199  

200 Parameters 

201 ---------- 

202 eod: 1-D or 2-D array 

203 The eod waveform to be analyzed. 

204 If an 1-D array, then this is the waveform and you 

205 need to also pass a time array or sampling rate in `ratetime`. 

206 If a 2-D array, then first column is time in seconds, second 

207 column the EOD waveform, third column, if present, is the 

208 standard error of the EOD waveform. Further columns are 

209 optional and are not used. 

210 ratetime: None or float or array of float 

211 If a 1-D array is passed on to `eod` then either the sampling 

212 rate in Hertz or the time array corresponding to `eod`. 

213 freq: float or 2-D array 

214 The frequency of the EOD or the list of harmonics (rows) with 

215 frequency and peak height (columns) as returned from 

216 `harmonics.harmonic_groups()`. 

217 n_harm: int 

218 Maximum number of harmonics used for the Fourier decomposition. 

219 power_n_harmonics: int 

220 Sum over the first `power_n_harmonics` harmonics for computing 

221 the total power. If 0 sum over all harmonics. 

222 n_harmonics: int 

223 The maximum power of higher harmonics is computed from 

224 harmonics higher than the maximum harmonics within the first 

225 three harmonics plus `n_harmonics`. 

226 flip_wave: 'auto', 'none', 'flip' 

227 - 'auto' flip waveform such that the larger extremum is positive. 

228 - 'flip' flip waveform. 

229 - 'none' do not flip waveform. 

230  

231 Returns 

232 ------- 

233 meod: 2-D array of floats 

234 The eod waveform. First column is time in seconds, second 

235 column the eod waveform. Further columns are kept from the 

236 input `eod`. And a column is added with the fit of the fourier 

237 series to the waveform. 

238 props: dict 

239 A dictionary with properties of the analyzed EOD waveform. 

240 

241 - type: set to 'wave'. 

242 - EODf: is set to the EOD fundamental frequency. 

243 - p-p-amplitude: peak-to-peak amplitude of the Fourier fit. 

244 - flipped: True if the waveform was flipped. 

245 - amplitude: amplitude factor of the Fourier fit. 

246 - noise: root-mean squared standard error mean of the averaged 

247 EOD waveform relative to the p-p amplitude. 

248 - rmserror: root-mean-square error between Fourier-fit and 

249 EOD waveform relative to the p-p amplitude. If larger than 

250 about 0.05 the data are bad. 

251 - ncrossings: number of zero crossings per period 

252 - peakwidth: width of the peak at the averaged amplitude relative 

253 to EOD period. 

254 - troughwidth: width of the trough at the averaged amplitude 

255 relative to EOD period. 

256 - minwidth: peakwidth or troughwidth, whichever is smaller. 

257 - leftpeak: time from positive zero crossing to peak relative 

258 to EOD period. 

259 - rightpeak: time from peak to negative zero crossing relative to 

260 EOD period. 

261 - lefttrough: time from negative zero crossing to trough relative to 

262 EOD period. 

263 - righttrough: time from trough to positive zero crossing relative to 

264 EOD period. 

265 - p-p-distance: time between peak and trough relative to EOD period. 

266 - min-p-p-distance: p-p-distance or EOD period minus p-p-distance, 

267 whichever is smaller, relative to EOD period. 

268 - relpeakampl: amplitude of peak or trough, whichever is larger, 

269 relative to p-p amplitude. 

270 - power: summed power of all harmonics of the extracted EOD waveform 

271 in decibel relative to one. 

272 - datapower: summed power of all harmonics of the original data in 

273 decibel relative to one. Only if `freq` is a list of harmonics. 

274 - thd: total harmonic distortion, i.e. square root of sum of 

275 amplitudes squared of harmonics relative to amplitude 

276 of fundamental.  

277 - dbdiff: smoothness of power spectrum as standard deviation of 

278 differences in decibel power. 

279 - maxdb: maximum power of higher harmonics relative to peak power 

280 in decibel. 

281 

282 spec_data: 2-D array of floats 

283 First six columns are from the spectrum of the extracted 

284 waveform. First column is the index of the harmonics, second 

285 column its frequency, third column its amplitude, fourth 

286 column its amplitude relative to the fundamental, fifth column 

287 is power of harmonics relative to fundamental in decibel, and 

288 sixth column the phase shift relative to the fundamental. 

289 If `freq` is a list of harmonics, a seventh column is added to 

290 `spec_data` that contains the powers of the harmonics from the 

291 original power spectrum of the raw data. Rows are the 

292 harmonics, first row is the fundamental frequency with index 

293 0, relative amplitude of one, relative power of 0dB, and phase 

294 shift of zero. 

295 error_str: string 

296 If fitting of the fourier series failed, 

297 this is reported in this string. 

298 

299 Raises 

300 ------ 

301 IndexError: 

302 EOD data is less than one period long. 

303 """ 

304 if eod.ndim == 2: 

305 if eod.shape[1] > 2: 

306 eeod = eod 

307 else: 

308 eeod = np.column_stack((eod, np.zeros(len(eod)))) 

309 else: 

310 if isinstance(ratetime, (list, tuple, np.ndarray)): 

311 time = ratetime 

312 else: 

313 time = np.arange(len(eod))/ratetime 

314 eeod = np.zeros((len(eod), 3)) 

315 eeod[:, 0] = time 

316 eeod[:, 1] = eod 

317 # storage: 

318 meod = np.zeros((eeod.shape[0], eeod.shape[1] + 1)) 

319 meod[:, :eeod.shape[1]] = eeod 

320 meod[:, -1] = np.nan 

321 

322 freq0 = freq 

323 if hasattr(freq, 'shape'): 

324 freq0 = freq[0][0] 

325 

326 error_str = '' 

327 

328 # subtract mean and flip: 

329 period = 1.0/freq0 

330 pinx = int(np.ceil(period/(meod[1,0]-meod[0,0]))) 

331 maxn = (len(meod)//pinx)*pinx 

332 if maxn < pinx: maxn = len(meod) 

333 offs = (len(meod) - maxn)//2 

334 meod[:, 1] -= np.mean(meod[offs:offs+pinx,1]) 

335 flipped = False 

336 if 'flip' in flip_wave or ('auto' in flip_wave and -np.min(meod[:, 1]) > np.max(meod[:, 1])): 

337 meod[:, 1] = -meod[:, 1] 

338 flipped = True 

339 

340 # move peak of waveform to zero: 

341 offs = len(meod)//4 

342 maxinx = offs+np.argmax(meod[offs:3*offs,1]) 

343 meod[:, 0] -= meod[maxinx,0] 

344 

345 # indices of exactly one or two periods around peak: 

346 if len(meod) < pinx: 

347 raise IndexError('data need to contain at least one EOD period') 

348 if len(meod) >= 2*pinx: 

349 i0 = maxinx - pinx if maxinx >= pinx else 0 

350 i1 = i0 + 2*pinx 

351 if i1 > len(meod): 

352 i1 = len(meod) 

353 i0 = i1 - 2*pinx 

354 else: 

355 i0 = maxinx - pinx//2 if maxinx >= pinx//2 else 0 

356 i1 = i0 + pinx 

357 

358 # subtract mean: 

359 meod[:, 1] -= np.mean(meod[i0:i1,1]) 

360 

361 # zero crossings: 

362 ui, di = threshold_crossings(meod[:, 1], 0.0) 

363 ut, dt = threshold_crossing_times(meod[:, 0], meod[:, 1], 0.0, ui, di) 

364 ut, dt = merge_events(ut, dt, 0.02/freq0) 

365 ncrossings = int(np.round((len(ut) + len(dt))/(meod[-1,0]-meod[0,0])/freq0)) 

366 if np.any(ut<0.0): 

367 up_time = ut[ut<0.0][-1] 

368 else: 

369 up_time = 0.0 

370 error_str += '%.1f Hz wave fish: no upward zero crossing. ' % freq0 

371 if np.any(dt>0.0): 

372 down_time = dt[dt>0.0][0] 

373 else: 

374 down_time = 0.0 

375 error_str += '%.1f Hz wave fish: no downward zero crossing. ' % freq0 

376 peak_width = down_time - up_time 

377 trough_width = period - peak_width 

378 peak_time = 0.0 

379 trough_time = meod[maxinx+np.argmin(meod[maxinx:maxinx+pinx,1]),0] 

380 phase1 = peak_time - up_time 

381 phase2 = down_time - peak_time 

382 phase3 = trough_time - down_time 

383 phase4 = up_time + period - trough_time 

384 distance = trough_time - peak_time 

385 min_distance = distance 

386 if distance > period/2: 

387 min_distance = period - distance 

388 

389 # fit fourier series: 

390 ampl = 0.5*(np.max(meod[:, 1])-np.min(meod[:, 1])) 

391 while n_harm > 1: 

392 params = [freq0] 

393 for i in range(1, n_harm+1): 

394 params.extend([ampl/i, 0.0]) 

395 try: 

396 popt, pcov = curve_fit(fourier_series, meod[i0:i1,0], 

397 meod[i0:i1,1], params, maxfev=2000) 

398 break 

399 except (RuntimeError, TypeError): 

400 error_str += '%.1f Hz wave fish: fit of fourier series failed for %d harmonics. ' % (freq0, n_harm) 

401 n_harm //= 2 

402 # store fourier fit: 

403 meod[:, -1] = fourier_series(meod[:, 0], *popt) 

404 # make all amplitudes positive: 

405 for i in range(n_harm): 

406 if popt[i*2+1] < 0.0: 

407 popt[i*2+1] *= -1.0 

408 popt[i*2+2] += np.pi 

409 # phases relative to fundamental: 

410 # phi0 = 2*pi*f0*dt <=> dt = phi0/(2*pi*f0) 

411 # phik = 2*pi*i*f0*dt = i*phi0 

412 phi0 = popt[2] 

413 for i in range(n_harm): 

414 popt[i*2+2] -= (i + 1)*phi0 

415 # all phases in the range -pi to pi: 

416 popt[i*2+2] %= 2*np.pi 

417 if popt[i*2+2] > np.pi: 

418 popt[i*2+2] -= 2*np.pi 

419 # store fourier spectrum: 

420 if hasattr(freq, 'shape'): 

421 n = n_harm 

422 n += np.sum(freq[:, 0] > (n_harm+0.5)*freq[0,0]) 

423 spec_data = np.zeros((n, 7)) 

424 spec_data[:, :] = np.nan 

425 k = 0 

426 for i in range(n_harm): 

427 while k < len(freq) and freq[k,0] < (i+0.5)*freq0: 

428 k += 1 

429 if k >= len(freq): 

430 break 

431 if freq[k,0] < (i+1.5)*freq0: 

432 spec_data[i,6] = freq[k,1] 

433 k += 1 

434 for i in range(n_harm, n): 

435 if k >= len(freq): 

436 break 

437 spec_data[i,:2] = [np.round(freq[k,0]/freq0)-1, freq[k,0]] 

438 spec_data[i,6] = freq[k,1] 

439 k += 1 

440 else: 

441 spec_data = np.zeros((n_harm, 6)) 

442 for i in range(n_harm): 

443 spec_data[i,:6] = [i, (i+1)*freq0, popt[i*2+1], popt[i*2+1]/popt[1], 

444 decibel((popt[i*2+1]/popt[1])**2.0), popt[i*2+2]] 

445 # smoothness of power spectrum: 

446 db_powers = decibel(spec_data[:n_harm,2]**2) 

447 db_diff = np.std(np.diff(db_powers)) 

448 # maximum relative power of higher harmonics: 

449 p_max = np.argmax(db_powers[:3]) 

450 db_powers -= db_powers[p_max] 

451 if len(db_powers[p_max+n_harmonics:]) == 0: 

452 max_harmonics_power = -100.0 

453 else: 

454 max_harmonics_power = np.max(db_powers[p_max+n_harmonics:]) 

455 # total harmonic distortion: 

456 thd = np.sqrt(np.nansum(spec_data[1:, 3])) 

457 

458 # peak-to-peak and trough amplitudes: 

459 ppampl = np.max(meod[i0:i1,1]) - np.min(meod[i0:i1,1]) 

460 relpeakampl = max(np.max(meod[i0:i1,1]), np.abs(np.min(meod[i0:i1,1])))/ppampl 

461 

462 # variance and fit error: 

463 rmssem = np.sqrt(np.mean(meod[i0:i1,2]**2.0))/ppampl if meod.shape[1] > 2 else None 

464 rmserror = np.sqrt(np.mean((meod[i0:i1,1] - meod[i0:i1,-1])**2.0))/ppampl 

465 

466 # store results: 

467 props = {} 

468 props['type'] = 'wave' 

469 props['EODf'] = freq0 

470 props['p-p-amplitude'] = ppampl 

471 props['flipped'] = flipped 

472 props['amplitude'] = 0.5*ppampl # remove it 

473 props['rmserror'] = rmserror 

474 if rmssem: 

475 props['noise'] = rmssem 

476 props['ncrossings'] = ncrossings 

477 props['peakwidth'] = peak_width/period 

478 props['troughwidth'] = trough_width/period 

479 props['minwidth'] = min(peak_width, trough_width)/period 

480 props['leftpeak'] = phase1/period 

481 props['rightpeak'] = phase2/period 

482 props['lefttrough'] = phase3/period 

483 props['righttrough'] = phase4/period 

484 props['p-p-distance'] = distance/period 

485 props['min-p-p-distance'] = min_distance/period 

486 props['relpeakampl'] = relpeakampl 

487 pnh = power_n_harmonics if power_n_harmonics > 0 else n_harm 

488 pnh = min(n_harm, pnh) 

489 props['power'] = decibel(np.sum(spec_data[:pnh,2]**2.0)) 

490 if hasattr(freq, 'shape'): 

491 props['datapower'] = decibel(np.sum(freq[:pnh,1])) 

492 props['thd'] = thd 

493 props['dbdiff'] = db_diff 

494 props['maxdb'] = max_harmonics_power 

495 

496 return meod, props, spec_data, error_str 

497 

498 

499def plot_wave_spectrum(axa, axp, spec, props, unit=None, 

500 ampl_style=dict(ls='', marker='o', color='tab:blue', markersize=6), 

501 ampl_stem_style=dict(color='tab:blue', alpha=0.5, lw=2), 

502 phase_style=dict(ls='', marker='p', color='tab:blue', markersize=6), 

503 phase_stem_style=dict(color='tab:blue', alpha=0.5, lw=2)): 

504 """Plot and annotate spectrum of wave EOD. 

505 

506 Parameters 

507 ---------- 

508 axa: matplotlib axes 

509 Axes for amplitude plot. 

510 axp: matplotlib axes 

511 Axes for phase plot. 

512 spec: 2-D array 

513 The amplitude spectrum of a single pulse as returned by 

514 `analyze_wave()`. First column is the index of the harmonics, 

515 second column its frequency, third column its amplitude, 

516 fourth column its amplitude relative to the fundamental, fifth 

517 column is power of harmonics relative to fundamental in 

518 decibel, and sixth column the phase shift relative to the 

519 fundamental. 

520 props: dict 

521 A dictionary with properties of the analyzed EOD waveform as 

522 returned by `analyze_wave()`. 

523 unit: string 

524 Optional unit of the data used for y-label. 

525 ampl_style: dict 

526 Properties of the markers of the amplitude plot. 

527 ampl_stem_style: dict 

528 Properties of the stems of the amplitude plot. 

529 phase_style: dict 

530 Properties of the markers of the phase plot. 

531 phase_stem_style: dict 

532 Properties of the stems of the phase plot. 

533 """ 

534 n = min(9, np.sum(np.isfinite(spec[:, 2]))) 

535 # amplitudes: 

536 markers, stemlines, _ = axa.stem(spec[:n, 0] + 1, spec[:n, 2], 

537 basefmt='none') 

538 setp(markers, clip_on=False, **ampl_style) 

539 setp(stemlines, **ampl_stem_style) 

540 axa.set_xlim(0.5, n + 0.5) 

541 axa.set_ylim(bottom=0) 

542 axa.xaxis.set_major_locator(MultipleLocator(1)) 

543 axa.tick_params('x', direction='out') 

544 if unit: 

545 axa.set_ylabel(f'Amplitude [{unit}]') 

546 else: 

547 axa.set_ylabel('Amplitude') 

548 # phases: 

549 phases = spec[:n, 5] 

550 phases[phases<0.0] = phases[phases<0.0] + 2*np.pi 

551 markers, stemlines, _ = axp.stem(spec[:n, 0] + 1, phases[:n], 

552 basefmt='none') 

553 setp(markers, clip_on=False, **phase_style) 

554 setp(stemlines, **phase_stem_style) 

555 axp.set_xlim(0.5, n + 0.5) 

556 axp.xaxis.set_major_locator(MultipleLocator(1)) 

557 axp.tick_params('x', direction='out') 

558 axp.set_ylim(0, 2*np.pi) 

559 axp.set_yticks([0, np.pi, 2*np.pi]) 

560 axp.set_yticklabels(['0', '\u03c0', '2\u03c0']) 

561 axp.set_xlabel('Harmonics') 

562 axp.set_ylabel('Phase') 

563 

564 

565def save_wave_eodfs(wave_eodfs, wave_indices, basename, **kwargs): 

566 """Save frequencies of wave EODs to file. 

567 

568 Parameters 

569 ---------- 

570 wave_eodfs: list of 2D arrays 

571 Each item is a matrix with the frequencies and powers 

572 (columns) of the fundamental and harmonics (rows) as returned 

573 by `harmonics.harmonic_groups()`. 

574 wave_indices: array 

575 Indices identifying each fish or NaN. 

576 If None no index column is inserted. 

577 basename: string or stream 

578 If string, path and basename of file. 

579 If `basename` does not have an extension, 

580 '-waveeodfs' and a file extension according to `kwargs` are appended. 

581 If stream, write EOD frequencies data into this stream. 

582 kwargs: 

583 Arguments passed on to `TableData.write()`. 

584 

585 Returns 

586 ------- 

587 filename: Path 

588 Path and full name of the written file in case of `basename` 

589 being a string. Otherwise, the file name and extension that 

590 would have been appended to a basename. 

591 

592 See Also 

593 -------- 

594 load_wave_eodfs() 

595 

596 """ 

597 eodfs = fundamental_freqs_and_power(wave_eodfs) 

598 td = TableData() 

599 if wave_indices is not None: 

600 td.append('index', '', '%d', 

601 value=[wi if wi >= 0 else np.nan 

602 for wi in wave_indices]) 

603 td.append('EODf', 'Hz', '%7.2f', value=eodfs[:, 0]) 

604 td.append('datapower', 'dB', '%7.2f', value=eodfs[:, 1]) 

605 ext = Path(basename).suffix if not hasattr(basename, 'write') else '' 

606 fp = '-waveeodfs' if not ext else '' 

607 return td.write_file_stream(basename, fp, **kwargs) 

608 

609 

610def load_wave_eodfs(file_path): 

611 """Load frequencies of wave EODs from file. 

612 

613 Parameters 

614 ---------- 

615 file_path: string 

616 Path of the file to be loaded. 

617 

618 Returns 

619 ------- 

620 eodfs: 2D array of floats 

621 EODfs and power of wave type fish. 

622 indices: array of ints 

623 Corresponding indices of fish, can contain negative numbers to 

624 indicate frequencies without fish. 

625 

626 Raises 

627 ------ 

628 FileNotFoundError: 

629 If `file_path` does not exist. 

630 

631 See Also 

632 -------- 

633 save_wave_eodfs() 

634 """ 

635 data = TableData(file_path) 

636 eodfs = data.array() 

637 if 'index' in data: 

638 indices = data[:, 'index'] 

639 indices[~np.isfinite(indices)] = -1 

640 indices = np.array(indices, dtype=int) 

641 eodfs = eodfs[:, 1:] 

642 else: 

643 indices = np.zeros(data.rows(), dtype=int) - 1 

644 return eodfs, indices 

645 

646 

647def save_wave_fish(eod_props, unit, basename, **kwargs): 

648 """Save properties of wave EODs to file. 

649 

650 Parameters 

651 ---------- 

652 eod_props: list of dict 

653 Properties of EODs as returned by `analyze_wave()` and 

654 `analyze_pulse()`. Only properties of wave fish are saved. 

655 unit: string 

656 Unit of the waveform data. 

657 basename: string or stream 

658 If string, path and basename of file. 

659 If `basename` does not have an extension, 

660 '-wavefish' and a file extension are appended. 

661 If stream, write wave fish properties into this stream. 

662 kwargs: 

663 Arguments passed on to `TableData.write()`. 

664 

665 Returns 

666 ------- 

667 filename: Path or None 

668 Path and full name of the written file in case of `basename` 

669 being a string. Otherwise, the file name and extension that 

670 would have been appended to a basename. 

671 None if no wave fish are contained in eod_props and 

672 consequently no file was written. 

673 

674 See Also 

675 -------- 

676 load_wave_fish() 

677 """ 

678 wave_props = [p for p in eod_props if p['type'] == 'wave'] 

679 if len(wave_props) == 0: 

680 return None 

681 td = TableData() 

682 if 'twin' in wave_props[0] or 'rate' in wave_props[0] or \ 

683 'nfft' in wave_props[0]: 

684 td.append_section('recording') 

685 if 'twin' in wave_props[0]: 

686 td.append('twin', 's', '%7.2f', value=wave_props) 

687 td.append('window', 's', '%7.2f', value=wave_props) 

688 td.append('winclipped', '%', '%.2f', 

689 value=wave_props, fac=100.0) 

690 if 'samplerate' in wave_props[0]: 

691 td.append('samplerate', 'kHz', '%.3f', 

692 value=wave_props, fac=0.001) 

693 if 'nfft' in wave_props[0]: 

694 td.append('nfft', '', '%d', value=wave_props) 

695 td.append('dfreq', 'Hz', '%.2f', value=wave_props) 

696 td.append_section('waveform') 

697 td.append('index', '', '%d', value=wave_props) 

698 td.append('EODf', 'Hz', '%7.2f', value=wave_props) 

699 td.append('p-p-amplitude', unit, '%.5f', value=wave_props) 

700 td.append('power', 'dB', '%7.2f', value=wave_props) 

701 if 'datapower' in wave_props[0]: 

702 td.append('datapower', 'dB', '%7.2f', value=wave_props) 

703 td.append('thd', '%', '%.2f', value=wave_props, fac=100.0) 

704 td.append('dbdiff', 'dB', '%7.2f', value=wave_props) 

705 td.append('maxdb', 'dB', '%7.2f', value=wave_props) 

706 if 'noise' in wave_props[0]: 

707 td.append('noise', '%', '%.1f', value=wave_props, fac=100.0) 

708 td.append('rmserror', '%', '%.2f', value=wave_props, fac=100.0) 

709 if 'clipped' in wave_props[0]: 

710 td.append('clipped', '%', '%.1f', value=wave_props, fac=100.0) 

711 td.append('flipped', '', '%d', value=wave_props) 

712 td.append('n', '', '%5d', value=wave_props) 

713 td.append_section('timing') 

714 td.append('ncrossings', '', '%d', value=wave_props) 

715 td.append('peakwidth', '%', '%.2f', value=wave_props, fac=100.0) 

716 td.append('troughwidth', '%', '%.2f', value=wave_props, fac=100.0) 

717 td.append('minwidth', '%', '%.2f', value=wave_props, fac=100.0) 

718 td.append('leftpeak', '%', '%.2f', value=wave_props, fac=100.0) 

719 td.append('rightpeak', '%', '%.2f', value=wave_props, fac=100.0) 

720 td.append('lefttrough', '%', '%.2f', value=wave_props, fac=100.0) 

721 td.append('righttrough', '%', '%.2f', value=wave_props, fac=100.0) 

722 td.append('p-p-distance', '%', '%.2f', value=wave_props, fac=100.0) 

723 td.append('min-p-p-distance', '%', '%.2f', 

724 value=wave_props, fac=100.0) 

725 td.append('relpeakampl', '%', '%.2f', value=wave_props, fac=100.0) 

726 ext = Path(basename).suffix if not hasattr(basename, 'write') else '' 

727 fp = '-wavefish' if not ext else '' 

728 return td.write_file_stream(basename, fp, **kwargs) 

729 

730 

731def load_wave_fish(file_path): 

732 """Load properties of wave EODs from file. 

733 

734 All times are scaled to seconds, all frequencies to Hertz and all 

735 percentages to fractions. 

736 

737 Parameters 

738 ---------- 

739 file_path: string 

740 Path of the file to be loaded. 

741 

742 Returns 

743 ------- 

744 eod_props: list of dict 

745 Properties of EODs. 

746 

747 Raises 

748 ------ 

749 FileNotFoundError: 

750 If `file_path` does not exist. 

751 

752 See Also