Coverage for src / thunderfish / fakefish.py: 79%
483 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-15 17:50 +0000
1"""Simulate EOD waveforms.
4## Species names
6- `species_name`: translate species ids to full species names.
7- `abbrv_genus()`: abbreviate genus in a species name.
10## Muscial intervals
12- `musical_intervals`: names and frequency ratios of musical intervals
15## Wavefish
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.
25## Pulsefish
27- `pulsefish_phases()`: position, amplitudes and standard deviations of EOD phases in pulsefish waveforms.
28- `pulsefish_spectrum()`: analytically computed single-pulse spectrum.
29- `pulsefish_waveform()`: compute a single pulsefish EOD.
30- `pulsefish_eods()`: simulate EODs of a pulse-type fish.
31- `normalize_pulsefish()`: normalize times and stdevs of pulse-type EOD waveform.
32- `export_pulsefish()`: serialize pulsefish parameter to python code.
35## Waveform generation
37- `generate_testfiles()`: generate recordings of various pulse EODs and their spectrum.
38- `generate_waveform()`: interactively generate audio file with simulated EOD waveforms.
39"""
41import sys
42import numpy as np
45species_name = dict(Sine='Sinewave',
46 Alepto='Apteronotus leptorhynchus',
47 Arostratus='Apteronotus rostratus',
48 Eigenmannia='Eigenmannia spec.',
49 Sternarchella='Sternarchella terminalis',
50 Sternopygus='Sternopygus dariensis')
51"""Translate species ids used by wavefish_harmonics and pulsefish_eodphases to full species names.
52"""
55def abbrv_genus(name):
56 """Abbreviate genus in a species name.
58 Parameters
59 ----------
60 name: string
61 Full species name of the form 'Genus species subspecies'
63 Returns
64 -------
65 name: string
66 The species name with abbreviated genus, i.e. 'G. species subspecies'
67 """
68 ns = name.split()
69 return ns[0][0] + '. ' + ' '.join(ns[1:])
72musical_intervals = {
73 'unison': (1/1, 1, 1, 0),
74 'minor second': (16/15, 16, 15, 1),
75 'major second': (9/8, 9, 8, 2),
76 'minor third': (6/5, 6, 5, 3),
77 'major third': (5/4, 5, 4, 4),
78 'forth': (4/3, 4, 3, 5),
79 'tritone': (45/32, 45, 32, 6), # =1.406, half way between forth and fifth: 17/6/2=1.4167, sqrt(2)=1.4142
80 'fifth': (3/2, 3, 2, 7),
81 'minor sixth': (8/5, 8, 5, 8),
82 'major sixth': (5/3, 5, 3, 9),
83 'subminor seventh': (7/4, 7, 4, 9.5),
84 'minor seventh': (9/5, 9, 5, 10),
85 'major seventh': (15/8, 15, 8, 11),
86 'octave': (2/1, 2, 1, 12),
87}
88"""Name, frequency ratio, nominator, denominator, and index of musical intervals
89"""
91# Amplitudes and phases of various wavefish species:
93Sine_harmonics = dict(amplitudes=(1.0,), phases=(0.5*np.pi,))
95Apteronotus_leptorhynchus_harmonics = \
96 dict(amplitudes=(0.90062, 0.15311, 0.072049, 0.012609, 0.011708),
97 phases=(1.3623, 2.3246, 0.9869, 2.6492, -2.6885))
99Apteronotus_rostratus_harmonics = \
100 dict(amplitudes=(0.64707, 0.43874, 0.063592, 0.07379, 0.040199, 0.023073,
101 0.0097678),
102 phases=(2.2988, 0.78876, -1.316, 2.2416, 2.0413, 1.1022,
103 -2.0513))
105Eigenmannia_harmonics = \
106 dict(amplitudes=(1.0087, 0.23201, 0.060524, 0.020175, 0.010087, 0.0080699),
107 phases=(1.3414, 1.3228, 2.9242, 2.8157, 2.6871, -2.8415))
109Sternarchella_terminalis_harmonics = \
110 dict(amplitudes=(0.11457, 0.4401, 0.41055, 0.20132, 0.061364, 0.011389,
111 0.0057985),
112 phases=(-2.7106, 2.4472, 1.6829, 0.79085, 0.119, -0.82355,
113 -1.9956))
115Sternopygus_dariensis_harmonics = \
116 dict(amplitudes=(0.98843, 0.41228, 0.047848, 0.11048, 0.022801, 0.030706,
117 0.019018),
118 phases=(1.4153, 1.3141, 3.1062, -2.3961, -1.9524, 0.54321,
119 1.6844))
121wavefish_harmonics = dict(Sine=Sine_harmonics,
122 Alepto=Apteronotus_leptorhynchus_harmonics,
123 Arostratus=Apteronotus_rostratus_harmonics,
124 Eigenmannia=Eigenmannia_harmonics,
125 Sternarchella=Sternarchella_terminalis_harmonics,
126 Sternopygus=Sternopygus_dariensis_harmonics)
127"""Amplitudes and phases of EOD waveforms of various species of wave-type electric fish."""
130def wavefish_spectrum(fish):
131 """Amplitudes and phases of a wavefish EOD.
133 Parameters
134 ----------
135 fish: string, dict or tuple of lists/arrays
136 Specify relative amplitudes and phases of the fundamental and its harmonics.
137 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary.
138 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys.
139 If tuple then the first element is the list of amplitudes and
140 the second one the list of relative phases in radians.
142 Returns
143 -------
144 amplitudes: array of floats
145 Amplitudes of the fundamental and its harmonics.
146 phases: array of floats
147 Phases in radians of the fundamental and its harmonics.
149 Raises
150 ------
151 KeyError:
152 Unknown fish.
153 IndexError:
154 Amplitudes and phases differ in length.
155 """
156 if isinstance(fish, (tuple, list)):
157 amplitudes = fish[0]
158 phases = fish[1]
159 elif isinstance(fish, dict):
160 amplitudes = fish['amplitudes']
161 phases = fish['phases']
162 else:
163 if not fish in wavefish_harmonics:
164 raise KeyError('unknown wavefish. Choose one of ' +
165 ', '.join(wavefish_harmonics.keys()))
166 amplitudes = wavefish_harmonics[fish]['amplitudes']
167 phases = wavefish_harmonics[fish]['phases']
168 if len(amplitudes) != len(phases):
169 raise IndexError('need exactly as many phases as amplitudes')
170 # remove NaNs:
171 for k in reversed(range(len(amplitudes))):
172 if np.isfinite(amplitudes[k]) or np.isfinite(phases[k]):
173 amplitudes = amplitudes[:k+1]
174 phases = phases[:k+1]
175 break
176 return amplitudes, phases
179def wavefish_eods(fish='Eigenmannia', frequency=100.0, rate=44100.0,
180 duration=1.0, phase0=0.0, noise_std=0.05):
181 """Simulate EOD waveform of a wave-type fish.
183 The waveform is constructed by superimposing sinewaves of integral
184 multiples of the fundamental frequency - the fundamental and its
185 harmonics. The fundamental frequency of the EOD (EODf) is given by
186 `frequency`. With `fish` relative amplitudes and phases of the
187 fundamental and its harmonics are specified.
189 The generated waveform is `duration` seconds long and is sampled with
190 `rate` Hertz. Gaussian white noise with a standard deviation of
191 `noise_std` is added to the generated waveform.
193 Parameters
194 ----------
195 fish: string, dict or tuple of lists/arrays
196 Specify relative amplitudes and phases of the fundamental and its harmonics.
197 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary.
198 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys.
199 If tuple then the first element is the list of amplitudes and
200 the second one the list of relative phases in radians.
201 frequency: float or array of floats
202 EOD frequency of the fish in Hertz. Either fixed number or array for
203 time-dependent frequencies.
204 rate: float
205 Sampling rate in Hertz.
206 duration: float
207 Duration of the generated data in seconds. Only used if frequency is scalar.
208 phase0: float
209 Phase offset of the EOD waveform in radians.
210 noise_std: float
211 Standard deviation of additive Gaussian white noise.
213 Returns
214 -------
215 data: array of floats
216 Generated data of a wave-type fish.
218 Raises
219 ------
220 KeyError:
221 Unknown fish.
222 IndexError:
223 Amplitudes and phases differ in length.
224 """
225 # get relative amplitude and phases:
226 amplitudes, phases = wavefish_spectrum(fish)
227 # compute phase:
228 if np.isscalar(frequency):
229 phase = np.arange(0, duration, 1.0/rate)
230 phase *= frequency
231 else:
232 phase = np.cumsum(frequency)/rate
233 # generate EOD:
234 data = np.zeros(len(phase))
235 for har, (ampl, phi) in enumerate(zip(amplitudes, phases)):
236 if np.isfinite(ampl) and np.isfinite(phi):
237 data += ampl * np.sin(2*np.pi*(har+1)*phase + phi - (har+1)*phase0)
238 # add noise:
239 data += noise_std * np.random.randn(len(data))
240 return data
243def normalize_wavefish(fish, mode='peak'):
244 """Normalize amplitudes and phases of wave-type EOD waveform.
246 The amplitudes and phases of the Fourier components are adjusted
247 such that the resulting EOD waveform has a peak-to-peak amplitude
248 of two and the peak of the waveform is at time zero (mode is set
249 to 'peak') or that the fundamental has an amplitude of one and a
250 phase of 0 (mode is set to 'zero').
252 Parameters
253 ----------
254 fish: string, dict or tuple of lists/arrays
255 Specify relative amplitudes and phases of the fundamental and its harmonics.
256 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary.
257 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys.
258 If tuple then the first element is the list of amplitudes and
259 the second one the list of relative phases in radians.
260 mode: 'peak' or 'zero'
261 How to normalize amplitude and phases:
262 - 'peak': normalize waveform to a peak-to-peak amplitude of two
263 and shift it such that its peak is at time zero.
264 - 'zero': scale amplitude of fundamental to one and its phase to zero.
266 Returns
267 -------
268 amplitudes: array of floats
269 Adjusted amplitudes of the fundamental and its harmonics.
270 phases: array of floats
271 Adjusted phases in radians of the fundamental and its harmonics.
273 """
274 # get relative amplitude and phases:
275 amplitudes, phases = wavefish_spectrum(fish)
276 if mode == 'zero':
277 newamplitudes = np.array(amplitudes)/amplitudes[0]
278 newphases = np.array([p+(k+1)*(-phases[0]) for k, p in enumerate(phases)])
279 newphases %= 2.0*np.pi
280 newphases[newphases>np.pi] -= 2.0*np.pi
281 return newamplitudes, newphases
282 else:
283 # generate waveform:
284 eodf = 100.0
285 rate = 100000.0
286 data = wavefish_eods(fish, eodf, rate, 2.0/eodf, noise_std=0.0)
287 # normalize amplitudes:
288 ampl = 0.5*(np.max(data) - np.min(data))
289 newamplitudes = np.array(amplitudes)/ampl
290 # shift phases:
291 deltat = np.argmax(data[:int(rate/eodf)])/rate
292 deltap = 2.0*np.pi*deltat*eodf
293 newphases = np.array([p+(k+1)*deltap for k, p in enumerate(phases)])
294 newphases %= 2.0*np.pi
295 newphases[newphases>np.pi] -= 2.0*np.pi
296 # return:
297 return newamplitudes, newphases
300def export_wavefish(fish, name='Unknown_harmonics', file=None):
301 """Serialize wavefish parameter to python code.
303 Add output to the wavefish_harmonics dictionary!
305 Parameters
306 ----------
307 fish: string, dict or tuple of lists/arrays
308 Specify relative amplitudes and phases of the fundamental and its harmonics.
309 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary.
310 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys.
311 If tuple then the first element is the list of amplitudes and
312 the second one the list of relative phases in radians.
313 name: string
314 Name of the dictionary to be written.
315 file: string or file or None
316 File name or open file object where to write wavefish dictionary.
318 Returns
319 -------
320 fish: dict
321 Dictionary with amplitudes and phases.
322 """
323 # get relative amplitude and phases:
324 amplitudes, phases = wavefish_spectrum(fish)
325 # write out dictionary:
326 if file is None:
327 file = sys.stdout
328 try:
329 file.write('')
330 closeit = False
331 except AttributeError:
332 file = open(file, 'w')
333 closeit = True
334 n = 6
335 file.write(name + ' = \\\n')
336 file.write(' dict(amplitudes=(')
337 file.write(', '.join([f'{a:.5g}' for a in amplitudes[:n]]))
338 for k in range(n, len(amplitudes), n):
339 file.write(',\n')
340 file.write(' ' * (9+12))
341 file.write(', '.join([f'{a:.5g}' for a in amplitudes[k:k+n]]))
342 file.write('),\n')
343 file.write(' ' * 9 + 'phases=(')
344 file.write(', '.join(['{p:.5g}' for p in phases[:n]]))
345 for k in range(n, len(phases), n):
346 file.write(',\n')
347 file.write(' ' * (9+8))
348 file.write(', '.join([f'{p:.5g}' for p in phases[k:k+n]]))
349 file.write('))\n')
350 if closeit:
351 file.close()
352 # return dictionary:
353 harmonics = dict(amplitudes=amplitudes,
354 phases=phases)
355 return harmonics
358def chirps(eodf=100.0, rate=44100.0, duration=1.0, chirp_freq=5.0,
359 chirp_size=100.0, chirp_width=0.01, chirp_kurtosis=1.0, chirp_contrast=0.05):
360 """Simulate frequency trace with chirps.
362 A chirp is modeled as a Gaussian frequency modulation.
363 The first chirp is placed at 0.5/chirp_freq.
365 Parameters
366 ----------
367 eodf: float
368 EOD frequency of the fish in Hertz.
369 rate: float
370 Sampling rate in Hertz.
371 duration: float
372 Duration of the generated data in seconds.
373 chirp_freq: float
374 Frequency of occurance of chirps in Hertz.
375 chirp_size: float
376 Size of the chirp (maximum frequency increase above eodf) in Hertz.
377 chirp_width: float
378 Width of the chirp at 10% height in seconds.
379 chirp_kurtosis: float:
380 Shape of the chirp. =1: Gaussian, >1: more rectangular, <1: more peaked.
381 chirp_contrast: float
382 Maximum amplitude reduction of EOD during chirp.
384 Returns
385 -------
386 frequency: array of floats
387 Generated frequency trace that can be passed on to wavefish_eods().
388 amplitude: array of floats
389 Generated amplitude modulation that can be used to multiply the trace generated by
390 wavefish_eods().
391 """
392 # baseline eod frequency and amplitude modulation:
393 n = len(np.arange(0, duration, 1.0/rate))
394 frequency = eodf * np.ones(n)
395 am = np.ones(n)
396 # time points for chirps:
397 chirp_period = 1.0/chirp_freq
398 chirp_times = np.arange(0.5*chirp_period, duration, chirp_period)
399 # chirp frequency waveform:
400 chirp_t = np.arange(-2.0*chirp_width, 2.0*chirp_width, 1./rate)
401 chirp_sig = 0.5*chirp_width / (2.0*np.log(10.0))**(0.5/chirp_kurtosis)
402 gauss = np.exp(-0.5*((chirp_t/chirp_sig)**2.0)**chirp_kurtosis)
403 # add chirps on baseline eodf:
404 for ct in chirp_times:
405 index = int(ct*rate)
406 i0 = index - len(gauss)//2
407 i1 = i0 + len(gauss)
408 gi0 = 0
409 gi1 = len(gauss)
410 if i0 < 0:
411 gi0 -= i0
412 i0 = 0
413 if i1 >= len(frequency):
414 gi1 -= i1 - len(frequency)
415 i1 = len(frequency)
416 frequency[i0:i1] += chirp_size * gauss[gi0:gi1]
417 am[i0:i1] -= chirp_contrast * gauss[gi0:gi1]
418 return frequency, am
421def rises(eodf=100.0, rate=44100.0, duration=1.0, rise_freq=0.1,
422 rise_size=10.0, rise_tau=1.0, decay_tau=10.0):
423 """Simulate frequency trace with rises.
425 A rise is modeled as a double exponential frequency modulation.
427 Parameters
428 ----------
429 eodf: float
430 EOD frequency of the fish in Hertz.
431 rate: float
432 Sampling rate in Hertz.
433 duration: float
434 Duration of the generated data in seconds.
435 rise_freq: float
436 Frequency of occurance of rises in Hertz.
437 rise_size: float
438 Size of the rise (frequency increase above eodf) in Hertz.
439 rise_tau: float
440 Time constant of the frequency increase of the rise in seconds.
441 decay_tau: float
442 Time constant of the frequency decay of the rise in seconds.
444 Returns
445 -------
446 data: array of floats
447 Generated frequency trace that can be passed on to wavefish_eods().
448 """
449 n = len(np.arange(0, duration, 1.0/rate))
450 # baseline eod frequency:
451 frequency = eodf * np.ones(n)
452 # time points for rises:
453 rise_period = 1.0/rise_freq
454 rise_times = np.arange(0.5*rise_period, duration, rise_period)
455 # rise frequency waveform:
456 rise_t = np.arange(0.0, 5.0*decay_tau, 1./rate)
457 rise = rise_size * (1.0-np.exp(-rise_t/rise_tau)) * np.exp(-rise_t/decay_tau)
458 # add rises on baseline eodf:
459 for r in rise_times:
460 index = int(r*rate)
461 if index+len(rise) > len(frequency):
462 rise_index = len(frequency)-index
463 frequency[index:index+rise_index] += rise[:rise_index]
464 break
465 else:
466 frequency[index:index+len(rise)] += rise
467 return frequency
470# Positions, amplitudes and standard deviations of EOD phases of various pulsefish species:
472Monophasic_phases = \
473 dict(times=(0,),
474 amplitudes=(1,),
475 stdevs=(0.0003,))
477Biphasic_phases = \
478 dict(times=(9e-05, 0.00049),
479 amplitudes=(1.1922, -0.95374),
480 stdevs=(0.0003, 0.00025))
482Triphasic_phases = \
483 dict(times=(3e-05, 0.00018, 0.00043),
484 amplitudes=(1.2382, -0.9906, 0.12382),
485 stdevs=(0.0001, 0.0001, 0.0002))
487pulsefish_eodphases = dict(Monophasic=Monophasic_phases,
488 Biphasic=Biphasic_phases,
489 Triphasic=Triphasic_phases)
490"""Positions, amplitudes, and standard deviations of Gaussians that
491 make up the EOD phases of pulse-type electric fish.
492"""
495def pulsefish_phases(fish):
496 """Position, amplitudes and standard deviations of EOD phases in pulsefish waveforms.
498 Parameters
499 ----------
500 fish: string, dict or tuple of floats/lists/arrays
501 Specify positions, amplitudes, and standard deviations
502 of Gaussians phases that are superimposed to simulate EOD waveforms
503 of pulse-type electric fishes.
504 If string then take positions, amplitudes and standard deviations
505 from the `pulsefish_eodphases` dictionary.
506 If dictionary then take pulse properties from the 'times', 'amlitudes'
507 and 'stdevs' keys.
508 If tuple then the first element is the list of phase positions,
509 the second is the list of corresponding amplitudes, and
510 the third one the list of corresponding standard deviations.
512 Returns
513 -------
514 times : array of floats
515 Positions of the phases.
516 amplitudes : array of floats
517 Amplitudes of the phases.
518 stdevs : array of floats
519 Standard deviations of the phases.
521 Raises
522 ------
523 KeyError:
524 Unknown fish.
525 IndexError:
526 Phase positions, amplitudes, or standard deviations differ in length.
527 """
528 if isinstance(fish, (tuple, list)):
529 times = fish[0]
530 amplitudes = fish[1]
531 stdevs = fish[2]
532 elif isinstance(fish, dict):
533 times = fish['times']
534 amplitudes = fish['amplitudes']
535 stdevs = fish['stdevs']
536 else:
537 if not fish in pulsefish_eodphases:
538 raise KeyError('unknown pulse-type fish. Choose one of ' +
539 ', '.join(pulsefish_eodphases.keys()))
540 phases = pulsefish_eodphases[fish]
541 times = phases['times']
542 amplitudes = phases['amplitudes']
543 stdevs = phases['stdevs']
544 if len(stdevs) != len(amplitudes) or len(stdevs) != len(times):
545 raise IndexError('need exactly as many standard deviations as amplitudes and times')
546 return times, amplitudes, stdevs
549def pulsefish_spectrum(fish, freqs):
550 """Analytically computed single-pulse spectrum.
552 The spectrum is the sum of the spectra of all of the phases. Each
553 phase is the spectrum of a Gaussian (which again is a Gaussian
554 centered at f=0) multiplied with a shift term to account for its
555 position relative to the first phase.
557 You get the amplitude spectral density by taking the absolute value.
558 The energy spectral density is the amplitude spectral density squared:
559 ```
560 spec = pulsefish_spectrum('Biphasic', freqs)
561 ampl = np.abs(spec)
562 energy = np.abs(spec)**2
563 ```
565 Parameters
566 ----------
567 fish: string, dict or tuple of floats/lists/arrays
568 Specify positions, amplitudes, and standard deviations
569 of Gaussians phases that are superimposed to simulate EOD waveforms
570 of pulse-type electric fishes.
571 If string then take positions, amplitudes and standard deviations
572 from the `pulsefish_eodphases` dictionary.
573 If dictionary then take pulse properties from the 'times', 'amlitudes'
574 and 'stdevs' keys.
575 If tuple then the first element is the list of phase positions,
576 the second is the list of corresponding amplitudes, and
577 the third one the list of corresponding standard deviations.
578 freqs: 1-D array of float
579 Frequencies at which the spectrum is computed.
581 Returns
582 -------
583 spectrum: 1-D array of complex
584 The one-sided complex-valued spectrum of the single pulse EOD.
585 The squared magnitude (the energy spectrum) has
586 the unit (x s)^2 = x^2 s/Hz.
588 """
589 times, ampls, stdevs = pulsefish_phases(fish)
590 spec = np.zeros(len(freqs), dtype=complex)
591 for dt, a, s in zip(times, ampls, stdevs):
592 gauss = a*np.sqrt(2*np.pi)*s*np.exp(-0.5*(2*np.pi*s*freqs)**2)
593 shift = np.exp(-2j*np.pi*freqs*dt)
594 spec += gauss*shift
595 spec *= np.sqrt(2) # because of one-sided spectrum
596 return spec
599def pulsefish_waveform(fish, t):
600 """Compute a single pulsefish EOD.
602 You may use pulsefish_parameter() to pass the paramters returned
603 by pulsefish_phases() to this function.
605 Parameters
606 ----------
607 fish: string, dict or tuple of floats/lists/arrays
608 Specify positions, amplitudes, and standard deviations
609 of Gaussians phases that are superimposed to simulate EOD waveforms
610 of pulse-type electric fishes.
611 If string then take positions, amplitudes and standard deviations
612 from the `pulsefish_eodphases` dictionary.
613 If dictionary then take pulse properties from the 'times', 'amlitudes'
614 and 'stdevs' keys.
615 If tuple then the first element is the list of phase positions,
616 the second is the list of corresponding amplitudes, and
617 the third one the list of corresponding standard deviations.
618 t: array of float
619 The time array over which the pulse waveform is evaluated.
621 Returns
622 -------
623 pulse: array of float
624 The pulse waveform for the times given in `t`.
626 """
627 times, ampls, stdevs = pulsefish_phases(fish)
628 pulse = np.zeros(len(t))
629 for dt, a, s in zip(times, ampls, stdevs):
630 pulse += a*np.exp(-0.5*((t - dt)/s)**2)
631 return pulse
634def pulsefish_eods(fish='Biphasic', frequency=100.0, rate=44100.0,
635 duration=1.0, noise_std=0.01, jitter_cv=0.1,
636 first_pulse=None):
637 """Simulate EODs of a pulse-type fish.
639 Pulses are spaced by 1/frequency, jittered as determined by
640 jitter_cv. Each pulse is a combination of Gaussian phases, whose
641 positions, amplitudes and widths are given by `fish`.
643 The generated waveform is duration seconds long and is sampled
644 with rate Hertz. Gaussian white noise with a standard deviation
645 of noise_std is added to the generated pulse train.
647 Parameters
648 ----------
649 fish: string, dict or tuple of floats/lists/arrays
650 Specify positions, amplitudes, and standard deviations
651 of Gaussians phases that are superimposed to simulate EOD waveforms
652 of pulse-type electric fishes.
653 If string then take positions, amplitudes and standard deviations
654 from the `pulsefish_eodphases` dictionary.
655 If dictionary then take pulse properties from the 'times', 'amlitudes'
656 and 'stdevs' keys.
657 If tuple then the first element is the list of phase positions,
658 the second is the list of corresponding amplitudes, and
659 the third one the list of corresponding standard deviations.
660 frequency: float
661 EOD frequency of the fish in Hz.
662 rate: float
663 Sampling Rate in Hz.
664 duration: float
665 Duration of the generated data in seconds.
666 noise_std: float
667 Standard deviation of additive Gaussian white noise.
668 jitter_cv: float
669 Gaussian distributed jitter of pulse times as coefficient of variation
670 of inter-pulse intervals.
671 first_pulse: float or None
672 The position of the first pulse. If None it is choosen automatically
673 depending on pulse width, jitter, and frequency.
675 Returns
676 -------
677 data: array of floats
678 Generated data of a pulse-type fish.
680 Raises
681 ------
682 KeyError:
683 Unknown fish.
684 IndexError:
685 Phase positions, amplitudes, or standard deviations differ in length.
687 """
688 # get phase properties:
689 phase_times, phase_amplitudes, phase_stdevs = pulsefish_phases(fish)
690 # time axis for single pulse:
691 min_time_inx = np.argmin(phase_times)
692 max_time_inx = np.argmax(phase_times)
693 tmax = max(np.abs(phase_times[min_time_inx]-4.0*phase_stdevs[min_time_inx]),
694 np.abs(phase_times[max_time_inx]+4.0*phase_stdevs[max_time_inx]))
695 x = np.arange(-tmax, tmax, 1.0/rate)
696 pulse_duration = x[-1] - x[0]
698 # generate a single pulse:
699 pulse = pulsefish_waveform((phase_times, phase_amplitudes, phase_stdevs), x)
700 #pulse = np.zeros(len(x))
701 #for time, ampl, std in zip(phase_times, phase_amplitudes, phase_stdevs):
702 # pulse += ampl * np.exp(-0.5*((x-time)/std)**2)
703 poffs = len(pulse)//2
705 # paste the pulse into the noise floor:
706 time = np.arange(0, duration, 1.0/rate)
707 data = np.random.randn(len(time)) * noise_std
708 period = 1.0/frequency
709 jitter_std = period * jitter_cv
710 if first_pulse is None:
711 first_pulse = np.max([pulse_duration, 3.0*jitter_std])
712 pulse_times = np.arange(first_pulse, duration, period )
713 pulse_times += jitter_std*np.random.randn(len(pulse_times))
714 pulse_indices = np.round(pulse_times * rate).astype(int)
715 for inx in pulse_indices[(pulse_indices>=poffs)&(pulse_indices-poffs+len(pulse)<len(data))]:
716 data[inx-poffs:inx-poffs+len(pulse)] += pulse
717 return data
720def normalize_pulsefish(fish):
721 """Normalize times and stdevs of pulse-type EOD waveform.
723 The positions and amplitudes of Gaussian phases are adjusted such
724 that the resulting EOD waveform has a maximum peak amplitude of one
725 and has the largest phase at time zero.
727 Parameters
728 ----------
729 fish: string, dict or tuple of floats/lists/arrays
730 Specify positions, amplitudes, and standard deviations
731 of Gaussians phases that are superimposed to simulate EOD waveforms
732 of pulse-type electric fishes.
733 If string then take positions, amplitudes and standard deviations
734 from the `pulsefish_eodphases` dictionary.
735 If dictionary then take pulse properties from the 'times', 'amlitudes'
736 and 'stdevs' keys.
737 If tuple then the first element is the list of phase positions,
738 the second is the list of corresponding amplitudes, and
739 the third one the list of corresponding standard deviations.
741 Returns
742 -------
743 fish: dict
744 Dictionary with adjusted times and standard deviations.
745 """
746 # get phase properties:
747 phase_times, phase_amplitudes, phase_stdevs = pulsefish_phases(fish)
748 # generate waveform:
749 eodf = 10.0
750 rate = 100000.0
751 first_pulse = 0.5/eodf
752 data = pulsefish_eods(fish, eodf, rate, 1.0/eodf, noise_std=0.0,
753 jitter_cv=0.0, first_pulse=first_pulse)
754 # maximum peak:
755 idx = np.argmax(np.abs(data))
756 # normalize amplitudes:
757 ampl = data[idx]
758 newamplitudes = np.array(phase_amplitudes)/ampl
759 # shift times:
760 deltat = idx/rate - first_pulse
761 newtimes = np.array(phase_times) - deltat
762 # store and return:
763 phases = dict(times=newtimes,
764 amplitudes=newamplitudes,
765 stdevs=phase_stdevs)
766 return phases
769def export_pulsefish(fish, name='Unknown_phases', file=None):
770 """Serialize pulsefish parameter to python code.
772 Add output to the pulsefish_eodphases dictionary!
774 Parameters
775 ----------
776 fish: string, dict or tuple of floats/lists/arrays
777 Specify positions, amplitudes, and standard deviations
778 of Gaussians phases that are superimposed to simulate EOD waveforms
779 of pulse-type electric fishes.
780 If string then take positions, amplitudes and standard deviations
781 from the `pulsefish_eodphases` dictionary.
782 If dictionary then take pulse properties from the 'times', 'amlitudes'
783 and 'stdevs' keys.
784 If tuple then the first element is the list of phase positions,
785 the second is the list of corresponding amplitudes, and
786 the third one the list of corresponding standard deviations.
787 name: string
788 Name of the dictionary to be written.
789 file: string or file or None
790 File name or open file object where to write pulsefish dictionary.
792 Returns
793 -------
794 fish: dict
795 Dictionary with phase times, amplitudes and standard deviations.
796 """
797 # get phase properties:
798 phase_times, phase_amplitudes, phase_stdevs = pulsefish_phases(fish)
799 # write out dictionary:
800 if file is None:
801 file = sys.stdout
802 try:
803 file.write('')
804 closeit = False
805 except AttributeError:
806 file = open(file, 'w')
807 closeit = True
808 n = 6
809 file.write(name + ' = \\\n')
810 file.write(' dict(times=(')
811 file.write(', '.join([f'{a:.5g}' for a in phase_times[:n]]))
812 for k in range(n, len(phase_times), n):
813 file.write(',\n')
814 file.write(' ' * (9+12))
815 file.write(', '.join([f'{a:.5g}' for a in phase_times[k:k+n]]))
816 if len(phase_times) == 1:
817 file.write(',')
818 file.write('),\n')
819 file.write(' ' * 9 + 'amplitudes=(')
820 file.write(', '.join([f'{p:.5g}' for p in phase_amplitudes[:n]]))
821 for k in range(n, len(phase_amplitudes), n):
822 file.write(',\n')
823 file.write(' ' * (9+8))
824 file.write(', '.join([f'{p:.5g}' for p in phase_amplitudes[k:k+n]]))
825 if len(phase_amplitudes) == 1:
826 file.write(',')
827 file.write('),\n')
828 file.write(' ' * 9 + 'stdevs=(')
829 file.write(', '.join([f'{p:.5g}' for p in phase_stdevs[:n]]))
830 for k in range(n, len(phase_stdevs), n):
831 file.write(',\n')
832 file.write(' ' * (9+8))
833 file.write(', '.join([f'{p:.5g}' for p in phase_stdevs[k:k+n]]))
834 if len(phase_stdevs) == 1:
835 file.write(',')
836 file.write('))\n')
837 if closeit:
838 file.close()
839 # return dictionary:
840 phases = dict(times=phase_times,
841 amplitudes=phase_amplitudes,
842 stdevs=phase_stdevs)
843 return phases
846def generate_waveform(filename):
847 """Interactively generate audio file with simulated EOD waveforms.
849 Parameters needed to generate EOD waveforms are take from console input.
851 Parameters
852 ----------
853 filename: string
854 Name of file inclusively extension (e.g. '.wav')
855 used to store the simulated EOD waveforms.
856 """
857 import os
858 from audioio import write_audio
859 from thunderlab.consoleinput import read, select, save_inputs
860 # generate file:
861 rate = read('Sampling rate in Hz', '44100', float, 1.0)
862 duration = read('Duration in seconds', '10', float, 0.001)
863 nfish = read('Number of fish', '1', int, 1)
864 ndata = read('Number of electrodes', '1', int, 1)
865 fish_spread = 1
866 if ndata > 1:
867 fish_spread = read('Number of electrodes fish are spread over', '2', int, 1)
868 data = np.random.randn(int(duration*rate), ndata)*0.01
869 fish_indices = np.random.randint(ndata, size=nfish)
870 eodt = 'a'
871 eodf = 800.0
872 eoda = 1.0
873 eodsig = 'n'
874 pulse_jitter = 0.1
875 chirp_freq = 5.0
876 chirp_size = 100.0
877 chirp_width = 0.01
878 chirp_kurtosis = 1.0
879 rise_freq = 0.1
880 rise_size = 10.0
881 rise_tau = 1.0
882 rise_decay_tau = 10.0
883 for k in range(nfish):
884 print('')
885 fish = 'Fish %d: ' % (k+1)
886 eodt = select(fish + 'EOD type', eodt, ['a', 'e', '1', '2', '3'],
887 ['Apteronotus', 'Eigenmannia',
888 'monophasic pulse', 'biphasic pulse', 'triphasic pulse'])
889 eodf = read(fish + 'EOD frequency in Hz', '%g'%eodf, float, 1.0, 3000.0)
890 eoda = read(fish + 'EOD amplitude', '%g'%eoda, float, 0.0, 10.0)
891 if eodt in 'ae':
892 eodsig = select(fish + 'Add communication signals', eodsig, ['n', 'c', 'r'],
893 ['fixed EOD', 'chirps', 'rises'])
894 eodfreq = eodf
895 if eodsig == 'c':
896 chirp_freq = read('Number of chirps per second', '%g'%chirp_freq, float, 0.001)
897 chirp_size = read('Size of chirp in Hz', '%g'%chirp_size, float, 1.0)
898 chirp_width = 0.001*read('Width of chirp in ms', '%g'%(1000.0*chirp_width), float, 1.0)
899 eodfreq, _ = chirps(eodf, rate, duration,
900 chirp_freq, chirp_size, chirp_width, chirp_kurtosis)
901 elif eodsig == 'r':
902 rise_freq = read('Number of rises per second', '%g'%rise_freq, float, 0.00001)
903 rise_size = read('Size of rise in Hz', '%g'%rise_size, float, 0.01)
904 rise_tau = read('Time-constant of rise onset in seconds', '%g'%rise_tau, float, 0.01)
905 rise_decay_tau = read('Time-constant of rise decay in seconds', '%g'%rise_decay_tau, float, 0.01)
906 eodfreq = rises(eodf, rate, duration,
907 rise_freq, rise_size, rise_tau, rise_decay_tau)
908 if eodt == 'a':
909 fishdata = eoda*wavefish_eods('Alepto', eodfreq, rate, duration,
910 phase0=0.0, noise_std=0.0)
911 elif eodt == 'e':
912 fishdata = eoda*wavefish_eods('Eigenmannia', eodfreq, rate,
913 duration, phase0=0.0, noise_std=0.0)
914 else:
915 pulse_jitter = read(fish + 'CV of pulse jitter', '%g'%pulse_jitter, float, 0.0, 2.0)
916 if eodt == '1':
917 fishdata = eoda*pulsefish_eods('Monophasic', eodf, rate, duration,
918 jitter_cv=pulse_jitter, noise_std=0.0)
919 elif eodt == '2':
920 fishdata = eoda*pulsefish_eods('Biphasic', eodf, rate, duration,
921 jitter_cv=pulse_jitter, noise_std=0.0)
922 elif eodt == '3':
923 fishdata = eoda*pulsefish_eods('Triphasic', eodf, rate, duration,
924 jitter_cv=pulse_jitter, noise_std=0.0)
925 i = fish_indices[k]
926 for j in range(fish_spread):
927 data[:, (i+j)%ndata] += fishdata*(0.2**j)
929 maxdata = np.max(np.abs(data))
930 write_audio(filename, 0.9*data/maxdata, rate)
931 input_file = os.path.splitext(filename)[0] + '.inp'
932 save_inputs(input_file)
933 print(f'\nWrote fakefish data to file "{filename}".')
936def generate_testfiles():
937 """Generate recordings of various pulse EODs and their spectrum.
939 The spectrum is analytically computed and thus can be used to test
940 analyis tools.
942 Three files are generated for each pulse waveform:
944 1. A wav file with the simulated recording.
945 2. A csv file with the spectra (see below for details).
946 3. A pdf file with a plot showing the averaged EOD pulse and the spectrum.
948 The csv file contains six columns separated by semicolons of the
949 single pulse spectra:
951 1. `f`: the frequency components in Hertz
952 2- `real`: real part of the Fourier spectrum
953 3. `imag`: imaginary part of the Fourier spectrum
954 4. `ampl`: amplitude spectrum, i.e. magnitude of Fourier spectrum
955 5. `energy`: energy spectrum, i.e. squared amplitude
956 6. `level`: energy sepctrum in decibel relative to maximum energy
958 """
959 import matplotlib.pyplot as plt
960 from scipy.signal import find_peaks
961 from audioio import write_audio
962 from thunderlab.eventdetection import snippets
963 from .eodanalysis import pulse_spectrum
965 np.seterr(all="ignore")
966 rng = np.random.default_rng()
967 rate = 50000.0
968 duration = 20.0
969 sigmat = 0.0002
970 monophasic = dict(name='monophasic',
971 times=(0,),
972 amplitudes=(1.0,),
973 stdevs=(sigmat,))
974 biphasic30 = dict(name='biphasic30',
975 times=(0, 2*sigmat),
976 amplitudes=(1.0, -0.3),
977 stdevs=(sigmat, sigmat))
978 biphasic60 = dict(name='biphasic60',
979 times=(0, 2*sigmat),
980 amplitudes=(1.0, -0.6),
981 stdevs=(sigmat, sigmat))
982 biphasic100 = dict(name='biphasic100',
983 times=(0, 2*sigmat),
984 amplitudes=(1.0, -1.0),
985 stdevs=(sigmat, sigmat))
986 triphasic = dict(name='triphasic', **Triphasic_phases)
987 for pulse in [monophasic, biphasic30, biphasic60, biphasic100, triphasic]:
988 print(pulse['name'], '...')
989 # fake recording:
990 eodf = rng.uniform(40.0, 120.0)
991 eoda = rng.uniform(0.8, 8.0)
992 data = eoda*pulsefish_eods(pulse, eodf, rate, duration,
993 jitter_cv=0.01, noise_std=0.002)
994 maxdata = np.max(np.abs(data))
995 fac = 0.9/maxdata
996 metadata = dict(gain=f'{1/fac:.3f}mV')
997 write_audio(pulse['name'] + '.wav', fac*data, rate, metadata)
998 print(f' wrote {pulse['name']}.wav')
999 # average EOD pulse:
1000 tmax = 0.002
1001 idxs, _ = find_peaks(data, prominence=0.9*eoda)
1002 iw = int(tmax*rate)
1003 snips = snippets(data, idxs, start=-iw, stop=iw)
1004 eod_data = np.mean(snips, 0)
1005 eod_time = (np.arange(len(eod_data)) - iw)/rate
1006 # analytic spectra:
1007 freqs = np.arange(0, rate/2, 1.0)
1008 spec = pulsefish_spectrum(pulse, freqs)
1009 spec *= eoda
1010 ampl = np.abs(spec)
1011 energy = ampl**2
1012 level = 10*np.log10(energy/np.max(energy))
1013 spec_data = np.zeros((len(freqs), 6))
1014 spec_data[:, 0] = freqs
1015 spec_data[:, 1] = np.real(spec)
1016 spec_data[:, 2] = np.imag(spec)
1017 spec_data[:, 3] = ampl
1018 spec_data[:, 4] = energy
1019 spec_data[:, 5] = level
1020 np.savetxt(pulse['name'] + '.csv', spec_data, fmt='%g', delimiter=';',
1021 header='f/Hz;real;imag;ampl;energy;level/dB')
1022 print(f' wrote {pulse['name']}.csv')
1023 # numerical spectrum:
1024 nfreqs, nenergy = pulse_spectrum(eod_data, rate, 1.0, 0.05)
1025 nlevel = 10*np.log10(nenergy/np.max(energy))
1026 # check normalization:
1027 print(f' integral over analytic energy spectrum: {np.sum(energy)*freqs[1]:9.3e}')
1028 print(f' integral over numeric energy spectrum : {np.sum(nenergy)*nfreqs[1]:9.3e}')
1029 print(f' integral over squared signal : {np.sum(eod_data**2)/rate:9.3e}')
1030 # plot waveform:
1031 pa = np.sum(eod_data[eod_data >= 0])/rate
1032 na = np.sum(-eod_data[eod_data <= 0])/rate
1033 balance = (pa - na)/(pa + na)
1034 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3),
1035 layout='constrained')
1036 fig.suptitle(pulse['name'])
1037 ax1.axhline(0, color='gray')
1038 ax1.plot(1000*eod_time, eod_data, color='C0', label='data')
1039 ax1.text(0.1, 0.05, f'polarity balance = {100*balance:.0f}%',
1040 transform=ax1.transAxes)
1041 ax1.set_xlim(-1000*tmax, 1000*tmax)
1042 ax1.set_ylim(-1.1*eoda, 1.1*eoda)
1043 ax1.set_xlabel('time [ms]')
1044 ax1.set_ylabel('averaged EOD')
1045 # plot spectrum:
1046 ip = np.argmax(level)
1047 fmax = freqs[ip]
1048 pmax = level[ip]
1049 fmaxpos = fmax if fmax > 1 else 1
1050 ax2.plot(nfreqs, nlevel + 0.5, color='C1', label='numeric')
1051 ax2.plot(freqs, level, color='C3', label='analytic')
1052 ax2.plot(fmaxpos, pmax, 'o', color='C3')
1053 ax2.set_xlim(1, 1e4)
1054 ax2.set_ylim(-60, 10)
1055 ax2.set_xscale('log')
1056 ax2.set_xlabel('frequency [Hz]')
1057 ax2.set_ylabel('energy [dB]')
1058 ax2.legend(loc='lower left', frameon=False)
1059 dc = level[0]
1060 ax2.text(2, dc - 4, f'{dc:.0f}dB')
1061 ax2.text(fmaxpos*1.05, pmax + 1, f'{fmax:.0f}Hz')
1062 fig.savefig(pulse['name'] + '.pdf')
1063 plt.show()
1064 plt.close(fig)
1065 print(f' wrote {pulse['name']}.pdf')
1068def demo():
1069 import matplotlib.pyplot as plt
1070 rate = 40000.0 # in Hz
1071 duration = 10.0 # in sec
1073 inset_len = 0.01 # in sec
1074 inset_indices = int(inset_len*rate)
1075 ws_fac = 0.1 # whitespace factor or ylim (between 0. and 1.)
1077 # generate data:
1078 eodf = 400.0
1079 wavefish = wavefish_eods('Alepto', eodf, rate, duration, noise_std=0.02)
1080 eodf = 650.0
1081 wavefish += 0.5*wavefish_eods('Eigenmannia', eodf, rate, duration)
1083 pulsefish = pulsefish_eods('Biphasic', 80.0, rate, duration,
1084 noise_std=0.02, jitter_cv=0.1, first_pulse=inset_len/2)
1085 time = np.arange(len(wavefish))/rate
1087 fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(19, 10))
1089 # get proper wavefish ylim
1090 ymin = np.min(wavefish)
1091 ymax = np.max(wavefish)
1092 dy = ws_fac*(ymax - ymin)
1093 ymin -= dy
1094 ymax += dy
1096 # complete wavefish:
1097 ax[0][0].set_title('Wavefish')
1098 ax[0][0].set_ylim(ymin, ymax)
1099 ax[0][0].plot(time, wavefish)
1101 # wavefish zoom in:
1102 ax[0][1].set_title('Wavefish ZOOM IN')
1103 ax[0][1].set_ylim(ymin, ymax)
1104 ax[0][1].plot(time[:inset_indices], wavefish[:inset_indices], '-o')
1106 # get proper pulsefish ylim
1107 ymin = np.min(pulsefish)
1108 ymax = np.max(pulsefish)
1109 dy = ws_fac*(ymax - ymin)
1110 ymin -= dy
1111 ymax += dy
1113 # complete pulsefish:
1114 ax[1][0].set_title('Pulsefish')
1115 ax[1][0].set_ylim(ymin, ymax)
1116 ax[1][0].plot(time, pulsefish)
1118 # pulsefish zoom in:
1119 ax[1][1].set_title('Pulsefish ZOOM IN')
1120 ax[1][1].set_ylim(ymin, ymax)
1121 ax[1][1].plot(time[:inset_indices], pulsefish[:inset_indices], '-o')
1123 for row in ax:
1124 for c_ax in row:
1125 c_ax.set_xlabel('Time [sec]')
1126 c_ax.set_ylabel('Amplitude')
1128 plt.tight_layout()
1130 # chirps:
1131 chirps_freq = chirps(600.0, rate, duration)
1132 chirps_data = wavefish_eods('Alepto', chirps_freq, rate)
1134 # rises:
1135 rises_freq = rises(600.0, rate, duration, rise_size=20.0)
1136 rises_data = wavefish_eods('Alepto', rises_freq, rate)
1138 nfft = 256
1139 fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(19, 10))
1140 ax[0].set_title('Chirps')
1141 ax[0].specgram(chirps_data, Fs=rate, NFFT=nfft, noverlap=nfft//16)
1142 time = np.arange(len(chirps_freq))/rate
1143 ax[0].plot(time[:-nfft//2], chirps_freq[nfft//2:], '-k', lw=2)
1144 ax[0].set_ylim(0.0, 3000.0)
1145 ax[0].set_ylabel('Frequency [Hz]')
1147 nfft = 4096
1148 ax[1].set_title('Rises')
1149 ax[1].specgram(rises_data, Fs=rate, NFFT=nfft, noverlap=nfft//2)
1150 time = np.arange(len(rises_freq))/rate
1151 ax[1].plot(time[:-nfft//4], rises_freq[nfft//4:], '-k', lw=2)
1152 ax[1].set_ylim(500.0, 700.0)
1153 ax[1].set_ylabel('Frequency [Hz]')
1154 ax[1].set_xlabel('Time [s]')
1155 plt.tight_layout()
1157 plt.show()
1160def main(args=[]):
1161 from .version import __year__
1163 if len(args) > 0:
1164 if (len(args) == 1 and args[0] != '-t') and args[0] != '-s':
1165 print('usage: fakefish [-h|--help] [-s audiofile] [-t]')
1166 print('')
1167 print('Without arguments, run a demo for illustrating fakefish functionality.')
1168 print('')
1169 print('-s audiofile: writes audiofile with user defined simulated electric fishes.')
1170 print('-t: write audiofiles for a number of pulse waveforms and corresponding analytic spectra in csv files for testing.')
1171 print('')
1172 print(f'by bendalab ({__year__})')
1173 elif args[0] == '-s':
1174 generate_waveform(args[1])
1175 elif args[0] == '-t':
1176 generate_testfiles()
1177 else:
1178 demo()
1181if __name__ == '__main__':
1182 import sys
1183 main(sys.argv[1:])