Coverage for src/thunderfish/fakefish.py: 98%

1"""Simulate EOD waveforms. 

2 

3 

4## Species names 

5 

6- `species_name`: translate species ids to full species names. 

7- `abbrv_genus()`: abbreviate genus in a species name. 

8 

9 

10## Muscial intervals 

11 

12- `musical_intervals`: names and frequency ratios of musical intervals 

13 

14 

15## Wavefish 

16 

17- `wavefish_spectrum()`: amplitudes and phases of a wavefish EOD. 

18- `wavefish_eods()`: simulate EOD waveform of a wave-type fish. 

19- `normalize_wavefish()`: normalize amplitudes and phases of EOD wave-type waveform. 

20- `export_wavefish()`: serialize wavefish parameter to file. 

21- `chirps()`: simulate frequency trace with chirps. 

22- `rises()`: simulate frequency trace with rises. 

23 

24 

25## Pulsefish 

26 

27- `pulsefish_eods()`: simulate EOD waveform of a pulse-type fish. 

28- `normalize_pulsefish()`: normalize times and stdevs of pulse-type EOD waveform. 

29- `export_pulsefish()`: serialize pulsefish parameter to file. 

30 

31 

32## Interactive waveform generation 

33 

34- `generate_waveform()`: interactively generate audio file with simulated EOD waveforms. 

35""" 

36 

37import sys 

38import numpy as np 

39 

40 

41species_name = dict(Sine='Sinewave', 

42 Alepto='Apteronotus leptorhynchus', 

43 Arostratus='Apteronotus rostratus', 

44 Eigenmannia='Eigenmannia spec.', 

45 Sternarchella='Sternarchella terminalis', 

46 Sternopygus='Sternopygus dariensis') 

47"""Translate species ids used by wavefish_harmonics and pulsefish_eodpeaks to full species names. 

48""" 

49 

50 

51def abbrv_genus(name): 

52 """Abbreviate genus in a species name. 

53 

54 Parameters 

55 ---------- 

56 name: string 

57 Full species name of the form 'Genus species subspecies' 

58 

59 Returns 

60 ------- 

61 name: string 

62 The species name with abbreviated genus, i.e. 'G. species subspecies' 

63 """ 

64 ns = name.split() 

65 return ns[0][0] + '. ' + ' '.join(ns[1:]) 

66 

67 

68musical_intervals = { 

69 'unison': (1/1, 1, 1, 0), 

70 'minor second': (16/15, 16, 15, 1), 

71 'major second': (9/8, 9, 8, 2), 

72 'minor third': (6/5, 6, 5, 3), 

73 'major third': (5/4, 5, 4, 4), 

74 'forth': (4/3, 4, 3, 5), 

75 'tritone': (45/32, 45, 32, 6), # =1.406, half way between forth and fifth: 17/6/2=1.4167, sqrt(2)=1.4142 

76 'fifth': (3/2, 3, 2, 7), 

77 'minor sixth': (8/5, 8, 5, 8), 

78 'major sixth': (5/3, 5, 3, 9), 

79 'subminor seventh': (7/4, 7, 4, 9.5), 

80 'minor seventh': (9/5, 9, 5, 10), 

81 'major seventh': (15/8, 15, 8, 11), 

82 'octave': (2/1, 2, 1, 12), 

83} 

84"""Name, frequency ratio, nominator, denominator, and index of musical intervals 

85""" 

86 

87# Amplitudes and phases of various wavefish species: 

88 

89Sine_harmonics = dict(amplitudes=(1.0,), phases=(0.5*np.pi,)) 

90 

91Apteronotus_leptorhynchus_harmonics = \ 

92 dict(amplitudes=(0.90062, 0.15311, 0.072049, 0.012609, 0.011708), 

93 phases=(1.3623, 2.3246, 0.9869, 2.6492, -2.6885)) 

94 

95Apteronotus_rostratus_harmonics = \ 

96 dict(amplitudes=(0.64707, 0.43874, 0.063592, 0.07379, 0.040199, 0.023073, 

97 0.0097678), 

98 phases=(2.2988, 0.78876, -1.316, 2.2416, 2.0413, 1.1022, 

99 -2.0513)) 

100 

101Eigenmannia_harmonics = \ 

102 dict(amplitudes=(1.0087, 0.23201, 0.060524, 0.020175, 0.010087, 0.0080699), 

103 phases=(1.3414, 1.3228, 2.9242, 2.8157, 2.6871, -2.8415)) 

104 

105Sternarchella_terminalis_harmonics = \ 

106 dict(amplitudes=(0.11457, 0.4401, 0.41055, 0.20132, 0.061364, 0.011389, 

107 0.0057985), 

108 phases=(-2.7106, 2.4472, 1.6829, 0.79085, 0.119, -0.82355, 

109 -1.9956)) 

110 

111Sternopygus_dariensis_harmonics = \ 

112 dict(amplitudes=(0.98843, 0.41228, 0.047848, 0.11048, 0.022801, 0.030706, 

113 0.019018), 

114 phases=(1.4153, 1.3141, 3.1062, -2.3961, -1.9524, 0.54321, 

115 1.6844)) 

116 

117wavefish_harmonics = dict(Sine=Sine_harmonics, 

118 Alepto=Apteronotus_leptorhynchus_harmonics, 

119 Arostratus=Apteronotus_rostratus_harmonics, 

120 Eigenmannia=Eigenmannia_harmonics, 

121 Sternarchella=Sternarchella_terminalis_harmonics, 

122 Sternopygus=Sternopygus_dariensis_harmonics) 

123"""Amplitudes and phases of EOD waveforms of various species of wave-type electric fish.""" 

124 

125 

126def wavefish_spectrum(fish): 

127 """Amplitudes and phases of a wavefish EOD. 

128 

129 Parameters 

130 ---------- 

131 fish: string, dict or tuple of lists/arrays 

132 Specify relative amplitudes and phases of the fundamental and its harmonics. 

133 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary. 

134 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys. 

135 If tuple then the first element is the list of amplitudes and 

136 the second one the list of relative phases in radians. 

137 

138 Returns 

139 ------- 

140 amplitudes: array of floats 

141 Amplitudes of the fundamental and its harmonics. 

142 phases: array of floats 

143 Phases in radians of the fundamental and its harmonics. 

144 

145 Raises 

146 ------ 

147 KeyError: 

148 Unknown fish. 

149 IndexError: 

150 Amplitudes and phases differ in length. 

151 """ 

152 if isinstance(fish, (tuple, list)): 

153 amplitudes = fish[0] 

154 phases = fish[1] 

155 elif isinstance(fish, dict): 

156 amplitudes = fish['amplitudes'] 

157 phases = fish['phases'] 

158 else: 

159 if not fish in wavefish_harmonics: 

160 raise KeyError('unknown wavefish. Choose one of ' + 

161 ', '.join(wavefish_harmonics.keys())) 

162 amplitudes = wavefish_harmonics[fish]['amplitudes'] 

163 phases = wavefish_harmonics[fish]['phases'] 

164 if len(amplitudes) != len(phases): 

165 raise IndexError('need exactly as many phases as amplitudes') 

166 # remove NaNs: 

167 for k in reversed(range(len(amplitudes))): 

168 if np.isfinite(amplitudes[k]) or np.isfinite(phases[k]): 

169 amplitudes = amplitudes[:k+1] 

170 phases = phases[:k+1] 

171 break 

172 return amplitudes, phases 

173 

174 

175def wavefish_eods(fish='Eigenmannia', frequency=100.0, rate=44100.0, 

176 duration=1.0, phase0=0.0, noise_std=0.05): 

177 """Simulate EOD waveform of a wave-type fish. 

178  

179 The waveform is constructed by superimposing sinewaves of integral 

180 multiples of the fundamental frequency - the fundamental and its 

181 harmonics. The fundamental frequency of the EOD (EODf) is given by 

182 `frequency`. With `fish` relative amplitudes and phases of the 

183 fundamental and its harmonics are specified. 

184 

185 The generated waveform is `duration` seconds long and is sampled with 

186 `rate` Hertz. Gaussian white noise with a standard deviation of 

187 `noise_std` is added to the generated waveform. 

188 

189 Parameters 

190 ---------- 

191 fish: string, dict or tuple of lists/arrays 

192 Specify relative amplitudes and phases of the fundamental and its harmonics. 

193 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary. 

194 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys. 

195 If tuple then the first element is the list of amplitudes and 

196 the second one the list of relative phases in radians. 

197 frequency: float or array of floats 

198 EOD frequency of the fish in Hertz. Either fixed number or array for 

199 time-dependent frequencies. 

200 rate: float 

201 Sampling rate in Hertz. 

202 duration: float 

203 Duration of the generated data in seconds. Only used if frequency is scalar. 

204 phase0: float 

205 Phase offset of the EOD waveform in radians. 

206 noise_std: float 

207 Standard deviation of additive Gaussian white noise. 

208 

209 Returns 

210 ------- 

211 data: array of floats 

212 Generated data of a wave-type fish. 

213 

214 Raises 

215 ------ 

216 KeyError: 

217 Unknown fish. 

218 IndexError: 

219 Amplitudes and phases differ in length. 

220 """ 

221 # get relative amplitude and phases: 

222 amplitudes, phases = wavefish_spectrum(fish) 

223 # compute phase: 

224 if np.isscalar(frequency): 

225 phase = np.arange(0, duration, 1.0/rate) 

226 phase *= frequency 

227 else: 

228 phase = np.cumsum(frequency)/rate 

229 # generate EOD: 

230 data = np.zeros(len(phase)) 

231 for har, (ampl, phi) in enumerate(zip(amplitudes, phases)): 

232 if np.isfinite(ampl) and np.isfinite(phi): 

233 data += ampl * np.sin(2*np.pi*(har+1)*phase + phi - (har+1)*phase0) 

234 # add noise: 

235 data += noise_std * np.random.randn(len(data)) 

236 return data 

237 

238 

239def normalize_wavefish(fish, mode='peak'): 

240 """Normalize amplitudes and phases of wave-type EOD waveform. 

241 

242 The amplitudes and phases of the Fourier components are adjusted 

243 such that the resulting EOD waveform has a peak-to-peak amplitude 

244 of two and the peak of the waveform is at time zero (mode is set 

245 to 'peak') or that the fundamental has an amplitude of one and a 

246 phase of 0 (mode is set to 'zero'). 

247 

248 Parameters 

249 ---------- 

250 fish: string, dict or tuple of lists/arrays 

251 Specify relative amplitudes and phases of the fundamental and its harmonics. 

252 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary. 

253 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys. 

254 If tuple then the first element is the list of amplitudes and 

255 the second one the list of relative phases in radians. 

256 mode: 'peak' or 'zero' 

257 How to normalize amplitude and phases: 

258 - 'peak': normalize waveform to a peak-to-peak amplitude of two 

259 and shift it such that its peak is at time zero. 

260 - 'zero': scale amplitude of fundamental to one and its phase to zero. 

261 

262 Returns 

263 ------- 

264 amplitudes: array of floats 

265 Adjusted amplitudes of the fundamental and its harmonics. 

266 phases: array of floats 

267 Adjusted phases in radians of the fundamental and its harmonics. 

268 

269 """ 

270 # get relative amplitude and phases: 

271 amplitudes, phases = wavefish_spectrum(fish) 

272 if mode == 'zero': 

273 newamplitudes = np.array(amplitudes)/amplitudes[0] 

274 newphases = np.array([p+(k+1)*(-phases[0]) for k, p in enumerate(phases)]) 

275 newphases %= 2.0*np.pi 

276 newphases[newphases>np.pi] -= 2.0*np.pi 

277 return newamplitudes, newphases 

278 else: 

279 # generate waveform: 

280 eodf = 100.0 

281 rate = 100000.0 

282 data = wavefish_eods(fish, eodf, rate, 2.0/eodf, noise_std=0.0) 

283 # normalize amplitudes: 

284 ampl = 0.5*(np.max(data) - np.min(data)) 

285 newamplitudes = np.array(amplitudes)/ampl 

286 # shift phases: 

287 deltat = np.argmax(data[:int(rate/eodf)])/rate 

288 deltap = 2.0*np.pi*deltat*eodf 

289 newphases = np.array([p+(k+1)*deltap for k, p in enumerate(phases)]) 

290 newphases %= 2.0*np.pi 

291 newphases[newphases>np.pi] -= 2.0*np.pi 

292 # return: 

293 return newamplitudes, newphases 

294 

295 

296def export_wavefish(fish, name='Unknown_harmonics', file=None): 

297 """Serialize wavefish parameter to python code. 

298 

299 Add output to the wavefish_harmonics dictionary! 

300 

301 Parameters 

302 ---------- 

303 fish: string, dict or tuple of lists/arrays 

304 Specify relative amplitudes and phases of the fundamental and its harmonics. 

305 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary. 

306 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys. 

307 If tuple then the first element is the list of amplitudes and 

308 the second one the list of relative phases in radians. 

309 name: string 

310 Name of the dictionary to be written. 

311 file: string or file or None 

312 File name or open file object where to write wavefish dictionary. 

313 

314 Returns 

315 ------- 

316 fish: dict 

317 Dictionary with amplitudes and phases. 

318 """ 

319 # get relative amplitude and phases: 

320 amplitudes, phases = wavefish_spectrum(fish) 

321 # write out dictionary: 

322 if file is None: 

323 file = sys.stdout 

324 try: 

325 file.write('') 

326 closeit = False 

327 except AttributeError: 

328 file = open(file, 'w') 

329 closeit = True 

330 n = 6 

331 file.write(name + ' = \\\n') 

332 file.write(' dict(amplitudes=(') 

333 file.write(', '.join([f'{a:.5g}' for a in amplitudes[:n]])) 

334 for k in range(n, len(amplitudes), n): 

335 file.write(',\n') 

336 file.write(' ' * (9+12)) 

337 file.write(', '.join([f'{a:.5g}' for a in amplitudes[k:k+n]])) 

338 file.write('),\n') 

339 file.write(' ' * 9 + 'phases=(') 

340 file.write(', '.join(['{p:.5g}' for p in phases[:n]])) 

341 for k in range(n, len(phases), n): 

342 file.write(',\n') 

343 file.write(' ' * (9+8)) 

344 file.write(', '.join([f'{p:.5g}' for p in phases[k:k+n]])) 

345 file.write('))\n') 

346 if closeit: 

347 file.close() 

348 # return dictionary: 

349 harmonics = dict(amplitudes=amplitudes, 

350 phases=phases) 

351 return harmonics 

352 

353 

354def chirps(eodf=100.0, rate=44100.0, duration=1.0, chirp_freq=5.0, 

355 chirp_size=100.0, chirp_width=0.01, chirp_kurtosis=1.0, chirp_contrast=0.05): 

356 """Simulate frequency trace with chirps. 

357 

358 A chirp is modeled as a Gaussian frequency modulation. 

359 The first chirp is placed at 0.5/chirp_freq. 

360 

361 Parameters 

362 ---------- 

363 eodf: float 

364 EOD frequency of the fish in Hertz. 

365 rate: float 

366 Sampling rate in Hertz. 

367 duration: float 

368 Duration of the generated data in seconds. 

369 chirp_freq: float 

370 Frequency of occurance of chirps in Hertz. 

371 chirp_size: float 

372 Size of the chirp (maximum frequency increase above eodf) in Hertz. 

373 chirp_width: float 

374 Width of the chirp at 10% height in seconds. 

375 chirp_kurtosis: float: 

376 Shape of the chirp. =1: Gaussian, >1: more rectangular, <1: more peaked. 

377 chirp_contrast: float 

378 Maximum amplitude reduction of EOD during chirp. 

379 

380 Returns 

381 ------- 

382 frequency: array of floats 

383 Generated frequency trace that can be passed on to wavefish_eods(). 

384 amplitude: array of floats 

385 Generated amplitude modulation that can be used to multiply the trace generated by 

386 wavefish_eods(). 

387 """ 

388 # baseline eod frequency and amplitude modulation: 

389 n = len(np.arange(0, duration, 1.0/rate)) 

390 frequency = eodf * np.ones(n) 

391 am = np.ones(n) 

392 # time points for chirps: 

393 chirp_period = 1.0/chirp_freq 

394 chirp_times = np.arange(0.5*chirp_period, duration, chirp_period) 

395 # chirp frequency waveform: 

396 chirp_t = np.arange(-2.0*chirp_width, 2.0*chirp_width, 1./rate) 

397 chirp_sig = 0.5*chirp_width / (2.0*np.log(10.0))**(0.5/chirp_kurtosis) 

398 gauss = np.exp(-0.5*((chirp_t/chirp_sig)**2.0)**chirp_kurtosis) 

399 # add chirps on baseline eodf: 

400 for ct in chirp_times: 

401 index = int(ct*rate) 

402 i0 = index - len(gauss)//2 

403 i1 = i0 + len(gauss) 

404 gi0 = 0 

405 gi1 = len(gauss) 

406 if i0 < 0: 

407 gi0 -= i0 

408 i0 = 0 

409 if i1 >= len(frequency): 

410 gi1 -= i1 - len(frequency) 

411 i1 = len(frequency) 

412 frequency[i0:i1] += chirp_size * gauss[gi0:gi1] 

413 am[i0:i1] -= chirp_contrast * gauss[gi0:gi1] 

414 return frequency, am 

415 

416 

417def rises(eodf=100.0, rate=44100.0, duration=1.0, rise_freq=0.1, 

418 rise_size=10.0, rise_tau=1.0, decay_tau=10.0): 

419 """Simulate frequency trace with rises. 

420 

421 A rise is modeled as a double exponential frequency modulation. 

422 

423 Parameters 

424 ---------- 

425 eodf: float 

426 EOD frequency of the fish in Hertz. 

427 rate: float 

428 Sampling rate in Hertz. 

429 duration: float 

430 Duration of the generated data in seconds. 

431 rise_freq: float 

432 Frequency of occurance of rises in Hertz. 

433 rise_size: float 

434 Size of the rise (frequency increase above eodf) in Hertz. 

435 rise_tau: float 

436 Time constant of the frequency increase of the rise in seconds. 

437 decay_tau: float 

438 Time constant of the frequency decay of the rise in seconds. 

439 

440 Returns 

441 ------- 

442 data: array of floats 

443 Generated frequency trace that can be passed on to wavefish_eods(). 

444 """ 

445 n = len(np.arange(0, duration, 1.0/rate)) 

446 # baseline eod frequency: 

447 frequency = eodf * np.ones(n) 

448 # time points for rises: 

449 rise_period = 1.0/rise_freq 

450 rise_times = np.arange(0.5*rise_period, duration, rise_period) 

451 # rise frequency waveform: 

452 rise_t = np.arange(0.0, 5.0*decay_tau, 1./rate) 

453 rise = rise_size * (1.0-np.exp(-rise_t/rise_tau)) * np.exp(-rise_t/decay_tau) 

454 # add rises on baseline eodf: 

455 for r in rise_times: 

456 index = int(r*rate) 

457 if index+len(rise) > len(frequency): 

458 rise_index = len(frequency)-index 

459 frequency[index:index+rise_index] += rise[:rise_index] 

460 break 

461 else: 

462 frequency[index:index+len(rise)] += rise 

463 return frequency 

464 

465 

466# Positions, amplitudes and standard deviations of peaks of various pulsefish species: 

467 

468Monophasic_peaks = \ 

469 dict(times=(0,), 

470 amplitudes=(1,), 

471 stdevs=(0.0003,)) 

472 

473Biphasic_peaks = \ 

474 dict(times=(9e-05, 0.00049), 

475 amplitudes=(1.1922, -0.95374), 

476 stdevs=(0.0003, 0.00025)) 

477 

478Triphasic_peaks = \ 

479 dict(times=(3e-05, 0.00018, 0.00043), 

480 amplitudes=(1.2382, -0.9906, 0.12382), 

481 stdevs=(0.0001, 0.0001, 0.0002)) 

482 

483pulsefish_eodpeaks = dict(Monophasic=Monophasic_peaks, 

484 Biphasic=Biphasic_peaks, 

485 Triphasic=Triphasic_peaks) 

486"""Standard deviations, amplitudes and positions of Gaussians that 

487 make up EOD waveforms of pulse-type electric fish. 

488""" 

489 

490 

491def pulsefish_peaks(fish): 

492 """Position, amplitudes and standard deviations of peaks in pulsefish EOD waveforms. 

493 

494 Parameters 

495 ---------- 

496 fish: string, dict or tuple of floats/lists/arrays 

497 Specify positions, amplitudes and standard deviations Gaussians peaks that are 

498 superimposed to simulate EOD waveforms of pulse-type electric fishes.  

499 If string then take positions, amplitudes and standard deviations  

500 from the `pulsefish_eodpeaks` dictionary. 

501 If dictionary then take pulse properties from the 'times', 'amlitudes' 

502 and 'stdevs' keys. 

503 If tuple then the first element is the list of peak positions, 

504 the second is the list of corresponding amplitudes, and 

505 the third one the list of corresponding standard deviations. 

506 

507 Returns 

508 ------- 

509 times : array of floats 

510 Positions of the peaks. 

511 amplitudes : array of floats 

512 Amplitudes of the peaks. 

513 stdevs : array of floats 

514 Standard deviations of the peaks. 

515 

516 Raises 

517 ------ 

518 KeyError: 

519 Unknown fish. 

520 IndexError: 

521 Peak positions, amplitudes, or standard deviations differ in length. 

522 """ 

523 if isinstance(fish, (tuple, list)): 

524 peak_times = fish[0] 

525 peak_amplitudes = fish[1] 

526 peak_stdevs = fish[2] 

527 elif isinstance(fish, dict): 

528 peak_times = fish['times'] 

529 peak_amplitudes = fish['amplitudes'] 

530 peak_stdevs = fish['stdevs'] 

531 else: 

532 if not fish in pulsefish_eodpeaks: 

533 raise KeyError('unknown pulse-type fish. Choose one of ' + 

534 ', '.join(pulsefish_eodpeaks.keys())) 

535 peak_times = pulsefish_eodpeaks[fish]['times'] 

536 peak_amplitudes = pulsefish_eodpeaks[fish]['amplitudes'] 

537 peak_stdevs = pulsefish_eodpeaks[fish]['stdevs'] 

538 if len(peak_stdevs) != len(peak_amplitudes) or len(peak_stdevs) != len(peak_times): 

539 raise IndexError('need exactly as many standard deviations as amplitudes and times') 

540 return peak_times, peak_amplitudes, peak_stdevs 

541 

542 

543def pulsefish_eods(fish='Biphasic', frequency=100.0, rate=44100.0, 

544 duration=1.0, noise_std=0.01, jitter_cv=0.1, 

545 first_pulse=None): 

546 """Simulate EOD waveform of a pulse-type fish. 

547 

548 Pulses are spaced by 1/frequency, jittered as determined by jitter_cv. Each pulse is 

549 a combination of Gaussian peaks, whose positions, amplitudes and widths are 

550 given by 'fish'. 

551 

552 The generated waveform is duration seconds long and is sampled with rate Hertz. 

553 Gaussian white noise with a standard deviation of noise_std is added to the generated 

554 pulse train. 

555 

556 Parameters 

557 ---------- 

558 fish: string, dict or tuple of floats/lists/arrays 

559 Specify positions, amplitudes and standard deviations Gaussians peaks that are 

560 superimposed to simulate EOD waveforms of pulse-type electric fishes.  

561 If string then take positions, amplitudes and standard deviations  

562 from the `pulsefish_eodpeaks` dictionary. 

563 If dictionary then take pulse properties from the 'times', 'amlitudes' 

564 and 'stdevs' keys. 

565 If tuple then the first element is the list of peak positions, 

566 the second is the list of corresponding amplitudes, and 

567 the third one the list of corresponding standard deviations. 

568 frequency: float 

569 EOD frequency of the fish in Hz. 

570 rate: float 

571 Sampling Rate in Hz. 

572 duration: float 

573 Duration of the generated data in seconds. 

574 noise_std: float 

575 Standard deviation of additive Gaussian white noise. 

576 jitter_cv: float 

577 Gaussian distributed jitter of pulse times as coefficient of variation 

578 of inter-pulse intervals. 

579 first_pulse: float or None 

580 The position of the first pulse. If None it is choosen automatically 

581 depending on pulse width, jitter, and frequency. 

582 

583 Returns 

584 ------- 

585 data: array of floats 

586 Generated data of a pulse-type fish. 

587 

588 Raises 

589 ------ 

590 KeyError: 

591 Unknown fish. 

592 IndexError: 

593 Peak positions, amplitudes, or standard deviations differ in length. 

594 """ 

595 # get peak properties: 

596 peak_times, peak_amplitudes, peak_stdevs = pulsefish_peaks(fish) 

597 # time axis for single pulse: 

598 min_time_inx = np.argmin(peak_times) 

599 max_time_inx = np.argmax(peak_times) 

600 tmax = max(np.abs(peak_times[min_time_inx]-4.0*peak_stdevs[min_time_inx]), 

601 np.abs(peak_times[max_time_inx]+4.0*peak_stdevs[max_time_inx])) 

602 x = np.arange(-tmax, tmax, 1.0/rate) 

603 pulse_duration = x[-1] - x[0] 

604 

605 # generate a single pulse: 

606 pulse = np.zeros(len(x)) 

607 for time, ampl, std in zip(peak_times, peak_amplitudes, peak_stdevs): 

608 pulse += ampl * np.exp(-0.5*((x-time)/std)**2) 

609 poffs = len(pulse)//2 

610 

611 # paste the pulse into the noise floor: 

612 time = np.arange(0, duration, 1.0/rate) 

613 data = np.random.randn(len(time)) * noise_std 

614 period = 1.0/frequency 

615 jitter_std = period * jitter_cv 

616 if first_pulse is None: 

617 first_pulse = np.max([pulse_duration, 3.0*jitter_std]) 

618 pulse_times = np.arange(first_pulse, duration, period ) 

619 pulse_times += jitter_std*np.random.randn(len(pulse_times)) 

620 pulse_indices = np.round(pulse_times * rate).astype(int) 

621 for inx in pulse_indices[(pulse_indices>=poffs)&(pulse_indices-poffs+len(pulse)<len(data))]: 

622 data[inx-poffs:inx-poffs+len(pulse)] += pulse 

623 return data 

624 

625 

626def normalize_pulsefish(fish): 

627 """Normalize times and stdevs of pulse-type EOD waveform. 

628 

629 The positions and amplitudes of Gaussian peaks are adjusted such 

630 that the resulting EOD waveform has a maximum peak amplitude of one 

631 and has the largest peak at time zero. 

632 

633 Parameters 

634 ---------- 

635 fish: string, dict or tuple of floats/lists/arrays 

636 Specify positions, amplitudes and standard deviations Gaussians peaks that are 

637 superimposed to simulate EOD waveforms of pulse-type electric fishes.  

638 If string then take positions, amplitudes and standard deviations  

639 from the `pulsefish_eodpeaks` dictionary. 

640 If dictionary then take pulse properties from the 'times', 'amlitudes' 

641 and 'stdevs' keys. 

642 If tuple then the first element is the list of peak positions, 

643 the second is the list of corresponding amplitudes, and 

644 the third one the list of corresponding standard deviations. 

645 

646 Returns 

647 ------- 

648 fish: dict 

649 Dictionary with adjusted times and standard deviations. 

650 """ 

651 # get peak properties: 

652 peak_times, peak_amplitudes, peak_stdevs = pulsefish_peaks(fish) 

653 # generate waveform: 

654 eodf = 10.0 

655 rate = 100000.0 

656 first_pulse = 0.5/eodf 

657 data = pulsefish_eods(fish, eodf, rate, 1.0/eodf, noise_std=0.0, 

658 jitter_cv=0.0, first_pulse=first_pulse) 

659 # maximum peak: 

660 idx = np.argmax(np.abs(data)) 

661 # normalize amplitudes: 

662 ampl = data[idx] 

663 newamplitudes = np.array(peak_amplitudes)/ampl 

664 # shift times: 

665 deltat = idx/rate - first_pulse 

666 newtimes = np.array(peak_times) - deltat 

667 # store and return: 

668 peaks = dict(times=newtimes, 

669 amplitudes=newamplitudes, 

670 stdevs=peak_stdevs) 

671 return peaks 

672 

673 

674def export_pulsefish(fish, name='Unknown_peaks', file=None): 

675 """Serialize pulsefish parameter to python code. 

676 

677 Add output to the pulsefish_eodpeaks dictionary! 

678 

679 Parameters 

680 ---------- 

681 fish: string, dict or tuple of floats/lists/arrays 

682 Specify positions, amplitudes and standard deviations Gaussians peaks that are 

683 superimposed to simulate EOD waveforms of pulse-type electric fishes.  

684 If string then take positions, amplitudes and standard deviations  

685 from the `pulsefish_eodpeaks` dictionary. 

686 If dictionary then take pulse properties from the 'times', 'amlitudes' 

687 and 'stdevs' keys. 

688 If tuple then the first element is the list of peak positions, 

689 the second is the list of corresponding amplitudes, and 

690 the third one the list of corresponding standard deviations. 

691 name: string 

692 Name of the dictionary to be written. 

693 file: string or file or None 

694 File name or open file object where to write pulsefish dictionary. 

695 

696 Returns 

697 ------- 

698 fish: dict 

699 Dictionary with peak times, amplitudes and standard deviations. 

700 """ 

701 # get peak properties: 

702 peak_times, peak_amplitudes, peak_stdevs = pulsefish_peaks(fish) 

703 # write out dictionary: 

704 if file is None: 

705 file = sys.stdout 

706 try: 

707 file.write('') 

708 closeit = False 

709 except AttributeError: 

710 file = open(file, 'w') 

711 closeit = True 

712 n = 6 

713 file.write(name + ' = \\\n') 

714 file.write(' dict(times=(') 

715 file.write(', '.join([f'{a:.5g}' for a in peak_times[:n]])) 

716 for k in range(n, len(peak_times), n): 

717 file.write(',\n') 

718 file.write(' ' * (9+12)) 

719 file.write(', '.join([f'{a:.5g}' for a in peak_times[k:k+n]])) 

720 if len(peak_times) == 1: 

721 file.write(',') 

722 file.write('),\n') 

723 file.write(' ' * 9 + 'amplitudes=(') 

724 file.write(', '.join([f'{p:.5g}' for p in peak_amplitudes[:n]])) 

725 for k in range(n, len(peak_amplitudes), n): 

726 file.write(',\n') 

727 file.write(' ' * (9+8)) 

728 file.write(', '.join([f'{p:.5g}' for p in peak_amplitudes[k:k+n]])) 

729 if len(peak_amplitudes) == 1: 

730 file.write(',') 

731 file.write('),\n') 

732 file.write(' ' * 9 + 'stdevs=(') 

733 file.write(', '.join([f'{p:.5g}' for p in peak_stdevs[:n]])) 

734 for k in range(n, len(peak_stdevs), n): 

735 file.write(',\n') 

736 file.write(' ' * (9+8)) 

737 file.write(', '.join([f'{p:.5g}' for p in peak_stdevs[k:k+n]])) 

738 if len(peak_stdevs) == 1: 

739 file.write(',') 

740 file.write('))\n') 

741 if closeit: 

742 file.close() 

743 # return dictionary: 

744 peaks = dict(times=peak_times, 

745 amplitudes=peak_amplitudes, 

746 stdevs=peak_stdevs) 

747 return peaks 

748 

749 

750def generate_waveform(filename): 

751 """Interactively generate audio file with simulated EOD waveforms. 

752 

753 Parameters needed to generate EOD waveforms are take from console input. 

754 

755 Parameters 

756 ---------- 

757 filename: string 

758 Name of file inclusively extension (e.g. '.wav') 

759 used to store the simulated EOD waveforms. 

760 """ 

761 import os 

762 from audioio import write_audio 

763 from thunderlab.consoleinput import read, select, save_inputs 

764 # generate file: 

765 rate = read('Sampling rate in Hz', '44100', float, 1.0) 

766 duration = read('Duration in seconds', '10', float, 0.001) 

767 nfish = read('Number of fish', '1', int, 1) 

768 ndata = read('Number of electrodes', '1', int, 1) 

769 fish_spread = 1 

770 if ndata > 1: 

771 fish_spread = read('Number of electrodes fish are spread over', '2', int, 1) 

772 data = np.random.randn(int(duration*rate), ndata)*0.01 

773 fish_indices = np.random.randint(ndata, size=nfish) 

774 eodt = 'a' 

775 eodf = 800.0 

776 eoda = 1.0 

777 eodsig = 'n' 

778 pulse_jitter = 0.1 

779 chirp_freq = 5.0 

780 chirp_size = 100.0 

781 chirp_width = 0.01 

782 chirp_kurtosis = 1.0 

783 rise_freq = 0.1 

784 rise_size = 10.0 

785 rise_tau = 1.0 

786 rise_decay_tau = 10.0 

787 for k in range(nfish): 

788 print('') 

789 fish = 'Fish %d: ' % (k+1) 

790 eodt = select(fish + 'EOD type', eodt, ['a', 'e', '1', '2', '3'], 

791 ['Apteronotus', 'Eigenmannia', 

792 'monophasic pulse', 'biphasic pulse', 'triphasic pulse']) 

793 eodf = read(fish + 'EOD frequency in Hz', '%g'%eodf, float, 1.0, 3000.0) 

794 eoda = read(fish + 'EOD amplitude', '%g'%eoda, float, 0.0, 10.0) 

795 if eodt in 'ae': 

796 eodsig = select(fish + 'Add communication signals', eodsig, ['n', 'c', 'r'], 

797 ['fixed EOD', 'chirps', 'rises']) 

798 eodfreq = eodf 

799 if eodsig == 'c': 

800 chirp_freq = read('Number of chirps per second', '%g'%chirp_freq, float, 0.001) 

801 chirp_size = read('Size of chirp in Hz', '%g'%chirp_size, float, 1.0) 

802 chirp_width = 0.001*read('Width of chirp in ms', '%g'%(1000.0*chirp_width), float, 1.0) 

803 eodfreq, _ = chirps(eodf, rate, duration, 

804 chirp_freq, chirp_size, chirp_width, chirp_kurtosis) 

805 elif eodsig == 'r': 

806 rise_freq = read('Number of rises per second', '%g'%rise_freq, float, 0.00001) 

807 rise_size = read('Size of rise in Hz', '%g'%rise_size, float, 0.01) 

808 rise_tau = read('Time-constant of rise onset in seconds', '%g'%rise_tau, float, 0.01) 

809 rise_decay_tau = read('Time-constant of rise decay in seconds', '%g'%rise_decay_tau, float, 0.01) 

810 eodfreq = rises(eodf, rate, duration, 

811 rise_freq, rise_size, rise_tau, rise_decay_tau) 

812 if eodt == 'a': 

813 fishdata = eoda*wavefish_eods('Alepto', eodfreq, rate, duration, 

814 phase0=0.0, noise_std=0.0) 

815 elif eodt == 'e': 

816 fishdata = eoda*wavefish_eods('Eigenmannia', eodfreq, rate, 

817 duration, phase0=0.0, noise_std=0.0) 

818 else: 

819 pulse_jitter = read(fish + 'CV of pulse jitter', '%g'%pulse_jitter, float, 0.0, 2.0) 

820 if eodt == '1': 

821 fishdata = eoda*pulsefish_eods('Monophasic', eodf, rate, duration, 

822 jitter_cv=pulse_jitter, noise_std=0.0) 

823 elif eodt == '2': 

824 fishdata = eoda*pulsefish_eods('Biphasic', eodf, rate, duration,