Coverage for src/thunderfish/fakefish.py: 98%
379 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 16:21 +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## Wavefish
12- `wavefish_spectrum()`: amplitudes and phases of a wavefish EOD.
13- `wavefish_eods()`: simulate EOD waveform of a wave-type fish.
14- `normalize_wavefish()`: normalize amplitudes and phases of EOD wave-type waveform.
15- `export_wavefish()`: serialize wavefish parameter to file.
16- `chirps()`: simulate frequency trace with chirps.
17- `rises()`: simulate frequency trace with rises.
20## Pulsefish
22- `pulsefish_eods()`: simulate EOD waveform of a pulse-type fish.
23- `normalize_pulsefish()`: normalize times and stdevs of pulse-type EOD waveform.
24- `export_pulsefish()`: serialize pulsefish parameter to file.
27## Interactive waveform generation
29- `generate_waveform()`: interactively generate audio file with simulated EOD waveforms.
30"""
32import sys
33import numpy as np
36species_name = dict(Sine='Sinewave',
37 Alepto='Apteronotus leptorhynchus',
38 Arostratus='Apteronotus rostratus',
39 Eigenmannia='Eigenmannia spec.',
40 Sternarchella='Sternarchella terminalis',
41 Sternopygus='Sternopygus dariensis')
42"""Translate species ids used by wavefish_harmonics and pulsefish_eodpeaks to full species names.
43"""
46def abbrv_genus(name):
47 """Abbreviate genus in a species name.
49 Parameters
50 ----------
51 name: string
52 Full species name of the form 'Genus species subspecies'
54 Returns
55 -------
56 name: string
57 The species name with abbreviated genus, i.e. 'G. species subspecies'
58 """
59 ns = name.split()
60 return ns[0][0] + '. ' + ' '.join(ns[1:])
63# Amplitudes and phases of various wavefish species:
65Sine_harmonics = dict(amplitudes=(1.0,), phases=(0.5*np.pi,))
67Apteronotus_leptorhynchus_harmonics = \
68 dict(amplitudes=(0.90062, 0.15311, 0.072049, 0.012609, 0.011708),
69 phases=(1.3623, 2.3246, 0.9869, 2.6492, -2.6885))
71Apteronotus_rostratus_harmonics = \
72 dict(amplitudes=(0.64707, 0.43874, 0.063592, 0.07379, 0.040199, 0.023073,
73 0.0097678),
74 phases=(2.2988, 0.78876, -1.316, 2.2416, 2.0413, 1.1022,
75 -2.0513))
77Eigenmannia_harmonics = \
78 dict(amplitudes=(1.0087, 0.23201, 0.060524, 0.020175, 0.010087, 0.0080699),
79 phases=(1.3414, 1.3228, 2.9242, 2.8157, 2.6871, -2.8415))
81Sternarchella_terminalis_harmonics = \
82 dict(amplitudes=(0.11457, 0.4401, 0.41055, 0.20132, 0.061364, 0.011389,
83 0.0057985),
84 phases=(-2.7106, 2.4472, 1.6829, 0.79085, 0.119, -0.82355,
85 -1.9956))
87Sternopygus_dariensis_harmonics = \
88 dict(amplitudes=(0.98843, 0.41228, 0.047848, 0.11048, 0.022801, 0.030706,
89 0.019018),
90 phases=(1.4153, 1.3141, 3.1062, -2.3961, -1.9524, 0.54321,
91 1.6844))
93wavefish_harmonics = dict(Sine=Sine_harmonics,
94 Alepto=Apteronotus_leptorhynchus_harmonics,
95 Arostratus=Apteronotus_rostratus_harmonics,
96 Eigenmannia=Eigenmannia_harmonics,
97 Sternarchella=Sternarchella_terminalis_harmonics,
98 Sternopygus=Sternopygus_dariensis_harmonics)
99"""Amplitudes and phases of EOD waveforms of various species of wave-type electric fish."""
102def wavefish_spectrum(fish):
103 """Amplitudes and phases of a wavefish EOD.
105 Parameters
106 ----------
107 fish: string, dict or tuple of lists/arrays
108 Specify relative amplitudes and phases of the fundamental and its harmonics.
109 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary.
110 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys.
111 If tuple then the first element is the list of amplitudes and
112 the second one the list of relative phases in radians.
114 Returns
115 -------
116 amplitudes: array of floats
117 Amplitudes of the fundamental and its harmonics.
118 phases: array of floats
119 Phases in radians of the fundamental and its harmonics.
121 Raises
122 ------
123 KeyError:
124 Unknown fish.
125 IndexError:
126 Amplitudes and phases differ in length.
127 """
128 if isinstance(fish, (tuple, list)):
129 amplitudes = fish[0]
130 phases = fish[1]
131 elif isinstance(fish, dict):
132 amplitudes = fish['amplitudes']
133 phases = fish['phases']
134 else:
135 if not fish in wavefish_harmonics:
136 raise KeyError('unknown wavefish. Choose one of ' +
137 ', '.join(wavefish_harmonics.keys()))
138 amplitudes = wavefish_harmonics[fish]['amplitudes']
139 phases = wavefish_harmonics[fish]['phases']
140 if len(amplitudes) != len(phases):
141 raise IndexError('need exactly as many phases as amplitudes')
142 # remove NaNs:
143 for k in reversed(range(len(amplitudes))):
144 if np.isfinite(amplitudes[k]) or np.isfinite(phases[k]):
145 amplitudes = amplitudes[:k+1]
146 phases = phases[:k+1]
147 break
148 return amplitudes, phases
151def wavefish_eods(fish='Eigenmannia', frequency=100.0, samplerate=44100.0,
152 duration=1.0, phase0=0.0, noise_std=0.05):
153 """Simulate EOD waveform of a wave-type fish.
155 The waveform is constructed by superimposing sinewaves of integral
156 multiples of the fundamental frequency - the fundamental and its
157 harmonics. The fundamental frequency of the EOD (EODf) is given by
158 `frequency`. With `fish` relative amplitudes and phases of the
159 fundamental and its harmonics are specified.
161 The generated waveform is `duration` seconds long and is sampled with
162 `samplerate` Hertz. Gaussian white noise with a standard deviation of
163 `noise_std` is added to the generated waveform.
165 Parameters
166 ----------
167 fish: string, dict or tuple of lists/arrays
168 Specify relative amplitudes and phases of the fundamental and its harmonics.
169 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary.
170 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys.
171 If tuple then the first element is the list of amplitudes and
172 the second one the list of relative phases in radians.
173 frequency: float or array of floats
174 EOD frequency of the fish in Hertz. Either fixed number or array for
175 time-dependent frequencies.
176 samplerate: float
177 Sampling rate in Hertz.
178 duration: float
179 Duration of the generated data in seconds. Only used if frequency is scalar.
180 phase0: float
181 Phase offset of the EOD waveform in radians.
182 noise_std: float
183 Standard deviation of additive Gaussian white noise.
185 Returns
186 -------
187 data: array of floats
188 Generated data of a wave-type fish.
190 Raises
191 ------
192 KeyError:
193 Unknown fish.
194 IndexError:
195 Amplitudes and phases differ in length.
196 """
197 # get relative amplitude and phases:
198 amplitudes, phases = wavefish_spectrum(fish)
199 # compute phase:
200 if np.isscalar(frequency):
201 phase = np.arange(0, duration, 1.0/samplerate)
202 phase *= frequency
203 else:
204 phase = np.cumsum(frequency)/samplerate
205 # generate EOD:
206 data = np.zeros(len(phase))
207 for har, (ampl, phi) in enumerate(zip(amplitudes, phases)):
208 if np.isfinite(ampl) and np.isfinite(phi):
209 data += ampl * np.sin(2*np.pi*(har+1)*phase + phi - (har+1)*phase0)
210 # add noise:
211 data += noise_std * np.random.randn(len(data))
212 return data
215def normalize_wavefish(fish, mode='peak'):
216 """Normalize amplitudes and phases of wave-type EOD waveform.
218 The amplitudes and phases of the Fourier components are adjusted
219 such that the resulting EOD waveform has a peak-to-peak amplitude
220 of two and the peak of the waveform is at time zero (mode is set
221 to 'peak') or that the fundamental has an amplitude of one and a
222 phase of 0 (mode is set to 'zero').
224 Parameters
225 ----------
226 fish: string, dict or tuple of lists/arrays
227 Specify relative amplitudes and phases of the fundamental and its harmonics.
228 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary.
229 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys.
230 If tuple then the first element is the list of amplitudes and
231 the second one the list of relative phases in radians.
232 mode: 'peak' or 'zero'
233 How to normalize amplitude and phases:
234 - 'peak': normalize waveform to a peak-to-peak amplitude of two
235 and shift it such that its peak is at time zero.
236 - 'zero': scale amplitude of fundamental to one and its phase to zero.
238 Returns
239 -------
240 amplitudes: array of floats
241 Adjusted amplitudes of the fundamental and its harmonics.
242 phases: array of floats
243 Adjusted phases in radians of the fundamental and its harmonics.
245 """
246 # get relative amplitude and phases:
247 amplitudes, phases = wavefish_spectrum(fish)
248 if mode == 'zero':
249 newamplitudes = np.array(amplitudes)/amplitudes[0]
250 newphases = np.array([p+(k+1)*(-phases[0]) for k, p in enumerate(phases)])
251 newphases %= 2.0*np.pi
252 newphases[newphases>np.pi] -= 2.0*np.pi
253 return newamplitudes, newphases
254 else:
255 # generate waveform:
256 eodf = 100.0
257 rate = 100000.0
258 data = wavefish_eods(fish, eodf, rate, 2.0/eodf, noise_std=0.0)
259 # normalize amplitudes:
260 ampl = 0.5*(np.max(data) - np.min(data))
261 newamplitudes = np.array(amplitudes)/ampl
262 # shift phases:
263 deltat = np.argmax(data[:int(rate/eodf)])/rate
264 deltap = 2.0*np.pi*deltat*eodf
265 newphases = np.array([p+(k+1)*deltap for k, p in enumerate(phases)])
266 newphases %= 2.0*np.pi
267 newphases[newphases>np.pi] -= 2.0*np.pi
268 # return:
269 return newamplitudes, newphases
272def export_wavefish(fish, name='Unknown_harmonics', file=None):
273 """Serialize wavefish parameter to python code.
275 Add output to the wavefish_harmonics dictionary!
277 Parameters
278 ----------
279 fish: string, dict or tuple of lists/arrays
280 Specify relative amplitudes and phases of the fundamental and its harmonics.
281 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary.
282 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys.
283 If tuple then the first element is the list of amplitudes and
284 the second one the list of relative phases in radians.
285 name: string
286 Name of the dictionary to be written.
287 file: string or file or None
288 File name or open file object where to write wavefish dictionary.
290 Returns
291 -------
292 fish: dict
293 Dictionary with amplitudes and phases.
294 """
295 # get relative amplitude and phases:
296 amplitudes, phases = wavefish_spectrum(fish)
297 # write out dictionary:
298 if file is None:
299 file = sys.stdout
300 try:
301 file.write('')
302 closeit = False
303 except AttributeError:
304 file = open(file, 'w')
305 closeit = True
306 n = 6
307 file.write(name + ' = \\\n')
308 file.write(' dict(amplitudes=(')
309 file.write(', '.join([f'{a:.5g}' for a in amplitudes[:n]]))
310 for k in range(n, len(amplitudes), n):
311 file.write(',\n')
312 file.write(' ' * (9+12))
313 file.write(', '.join([f'{a:.5g}' for a in amplitudes[k:k+n]]))
314 file.write('),\n')
315 file.write(' ' * 9 + 'phases=(')
316 file.write(', '.join(['{p:.5g}' for p in phases[:n]]))
317 for k in range(n, len(phases), n):
318 file.write(',\n')
319 file.write(' ' * (9+8))
320 file.write(', '.join([f'{p:.5g}' for p in phases[k:k+n]]))
321 file.write('))\n')
322 if closeit:
323 file.close()
324 # return dictionary:
325 harmonics = dict(amplitudes=amplitudes,
326 phases=phases)
327 return harmonics
330def chirps(eodf=100.0, samplerate=44100.0, duration=1.0, chirp_freq=5.0,
331 chirp_size=100.0, chirp_width=0.01, chirp_kurtosis=1.0, chirp_contrast=0.05):
332 """Simulate frequency trace with chirps.
334 A chirp is modeled as a Gaussian frequency modulation.
335 The first chirp is placed at 0.5/chirp_freq.
337 Parameters
338 ----------
339 eodf: float
340 EOD frequency of the fish in Hertz.
341 samplerate: float
342 Sampling rate in Hertz.
343 duration: float
344 Duration of the generated data in seconds.
345 chirp_freq: float
346 Frequency of occurance of chirps in Hertz.
347 chirp_size: float
348 Size of the chirp (maximum frequency increase above eodf) in Hertz.
349 chirp_width: float
350 Width of the chirp at 10% height in seconds.
351 chirp_kurtosis: float:
352 Shape of the chirp. =1: Gaussian, >1: more rectangular, <1: more peaked.
353 chirp_contrast: float
354 Maximum amplitude reduction of EOD during chirp.
356 Returns
357 -------
358 frequency: array of floats
359 Generated frequency trace that can be passed on to wavefish_eods().
360 amplitude: array of floats
361 Generated amplitude modulation that can be used to multiply the trace generated by
362 wavefish_eods().
363 """
364 # baseline eod frequency and amplitude modulation:
365 n = len(np.arange(0, duration, 1.0/samplerate))
366 frequency = eodf * np.ones(n)
367 am = np.ones(n)
368 # time points for chirps:
369 chirp_period = 1.0/chirp_freq
370 chirp_times = np.arange(0.5*chirp_period, duration, chirp_period)
371 # chirp frequency waveform:
372 chirp_t = np.arange(-2.0*chirp_width, 2.0*chirp_width, 1./samplerate)
373 chirp_sig = 0.5*chirp_width / (2.0*np.log(10.0))**(0.5/chirp_kurtosis)
374 gauss = np.exp(-0.5*((chirp_t/chirp_sig)**2.0)**chirp_kurtosis)
375 # add chirps on baseline eodf:
376 for ct in chirp_times:
377 index = int(ct*samplerate)
378 i0 = index - len(gauss)//2
379 i1 = i0 + len(gauss)
380 gi0 = 0
381 gi1 = len(gauss)
382 if i0 < 0:
383 gi0 -= i0
384 i0 = 0
385 if i1 >= len(frequency):
386 gi1 -= i1 - len(frequency)
387 i1 = len(frequency)
388 frequency[i0:i1] += chirp_size * gauss[gi0:gi1]
389 am[i0:i1] -= chirp_contrast * gauss[gi0:gi1]
390 return frequency, am
393def rises(eodf=100.0, samplerate=44100.0, duration=1.0, rise_freq=0.1,
394 rise_size=10.0, rise_tau=1.0, decay_tau=10.0):
395 """Simulate frequency trace with rises.
397 A rise is modeled as a double exponential frequency modulation.
399 Parameters
400 ----------
401 eodf: float
402 EOD frequency of the fish in Hertz.
403 samplerate: float
404 Sampling rate in Hertz.
405 duration: float
406 Duration of the generated data in seconds.
407 rise_freq: float
408 Frequency of occurance of rises in Hertz.
409 rise_size: float
410 Size of the rise (frequency increase above eodf) in Hertz.
411 rise_tau: float
412 Time constant of the frequency increase of the rise in seconds.
413 decay_tau: float
414 Time constant of the frequency decay of the rise in seconds.
416 Returns
417 -------
418 data: array of floats
419 Generated frequency trace that can be passed on to wavefish_eods().
420 """
421 n = len(np.arange(0, duration, 1.0/samplerate))
422 # baseline eod frequency:
423 frequency = eodf * np.ones(n)
424 # time points for rises:
425 rise_period = 1.0/rise_freq
426 rise_times = np.arange(0.5*rise_period, duration, rise_period)
427 # rise frequency waveform:
428 rise_t = np.arange(0.0, 5.0*decay_tau, 1./samplerate)
429 rise = rise_size * (1.0-np.exp(-rise_t/rise_tau)) * np.exp(-rise_t/decay_tau)
430 # add rises on baseline eodf:
431 for r in rise_times:
432 index = int(r*samplerate)
433 if index+len(rise) > len(frequency):
434 rise_index = len(frequency)-index
435 frequency[index:index+rise_index] += rise[:rise_index]
436 break
437 else:
438 frequency[index:index+len(rise)] += rise
439 return frequency
442# Positions, amplitudes and standard deviations of peaks of various pulsefish species:
444Monophasic_peaks = \
445 dict(times=(0,),
446 amplitudes=(1,),
447 stdevs=(0.0003,))
449Biphasic_peaks = \
450 dict(times=(9e-05, 0.00049),
451 amplitudes=(1.1922, -0.95374),
452 stdevs=(0.0003, 0.00025))
454Triphasic_peaks = \
455 dict(times=(3e-05, 0.00018, 0.00043),
456 amplitudes=(1.2382, -0.9906, 0.12382),
457 stdevs=(0.0001, 0.0001, 0.0002))
459pulsefish_eodpeaks = dict(Monophasic=Monophasic_peaks,
460 Biphasic=Biphasic_peaks,
461 Triphasic=Triphasic_peaks)
462"""Standard deviations, amplitudes and positions of Gaussians that
463 make up EOD waveforms of pulse-type electric fish.
464"""
467def pulsefish_peaks(fish):
468 """Position, amplitudes and standard deviations of peaks in pulsefish EOD waveforms.
470 Parameters
471 ----------
472 fish: string, dict or tuple of floats/lists/arrays
473 Specify positions, amplitudes and standard deviations Gaussians peaks that are
474 superimposed to simulate EOD waveforms of pulse-type electric fishes.
475 If string then take positions, amplitudes and standard deviations
476 from the `pulsefish_eodpeaks` dictionary.
477 If dictionary then take pulse properties from the 'times', 'amlitudes'
478 and 'stdevs' keys.
479 If tuple then the first element is the list of peak positions,
480 the second is the list of corresponding amplitudes, and
481 the third one the list of corresponding standard deviations.
483 Returns
484 -------
485 times : array of floats
486 Positions of the peaks.
487 amplitudes : array of floats
488 Amplitudes of the peaks.
489 stdevs : array of floats
490 Standard deviations of the peaks.
492 Raises
493 ------
494 KeyError:
495 Unknown fish.
496 IndexError:
497 Peak positions, amplitudes, or standard deviations differ in length.
498 """
499 if isinstance(fish, (tuple, list)):
500 peak_times = fish[0]
501 peak_amplitudes = fish[1]
502 peak_stdevs = fish[2]
503 elif isinstance(fish, dict):
504 peak_times = fish['times']
505 peak_amplitudes = fish['amplitudes']
506 peak_stdevs = fish['stdevs']
507 else:
508 if not fish in pulsefish_eodpeaks:
509 raise KeyError('unknown pulse-type fish. Choose one of ' +
510 ', '.join(pulsefish_eodpeaks.keys()))
511 peak_times = pulsefish_eodpeaks[fish]['times']
512 peak_amplitudes = pulsefish_eodpeaks[fish]['amplitudes']
513 peak_stdevs = pulsefish_eodpeaks[fish]['stdevs']
514 if len(peak_stdevs) != len(peak_amplitudes) or len(peak_stdevs) != len(peak_times):
515 raise IndexError('need exactly as many standard deviations as amplitudes and times')
516 return peak_times, peak_amplitudes, peak_stdevs
519def pulsefish_eods(fish='Biphasic', frequency=100.0, samplerate=44100.0,
520 duration=1.0, noise_std=0.01, jitter_cv=0.1,
521 first_pulse=None):
522 """Simulate EOD waveform of a pulse-type fish.
524 Pulses are spaced by 1/frequency, jittered as determined by jitter_cv. Each pulse is
525 a combination of Gaussian peaks, whose positions, amplitudes and widths are
526 given by 'fish'.
528 The generated waveform is duration seconds long and is sampled with samplerate Hertz.
529 Gaussian white noise with a standard deviation of noise_std is added to the generated
530 pulse train.
532 Parameters
533 ----------
534 fish: string, dict or tuple of floats/lists/arrays
535 Specify positions, amplitudes and standard deviations Gaussians peaks that are
536 superimposed to simulate EOD waveforms of pulse-type electric fishes.
537 If string then take positions, amplitudes and standard deviations
538 from the `pulsefish_eodpeaks` dictionary.
539 If dictionary then take pulse properties from the 'times', 'amlitudes'
540 and 'stdevs' keys.
541 If tuple then the first element is the list of peak positions,
542 the second is the list of corresponding amplitudes, and
543 the third one the list of corresponding standard deviations.
544 frequency: float
545 EOD frequency of the fish in Hz.
546 samplerate: float
547 Sampling Rate in Hz.
548 duration: float
549 Duration of the generated data in seconds.
550 noise_std: float
551 Standard deviation of additive Gaussian white noise.
552 jitter_cv: float
553 Gaussian distributed jitter of pulse times as coefficient of variation
554 of inter-pulse intervals.
555 first_pulse: float or None
556 The position of the first pulse. If None it is choosen automatically
557 depending on pulse width, jitter, and frequency.
559 Returns
560 -------
561 data: array of floats
562 Generated data of a pulse-type fish.
564 Raises
565 ------
566 KeyError:
567 Unknown fish.
568 IndexError:
569 Peak positions, amplitudes, or standard deviations differ in length.
570 """
571 # get peak properties:
572 peak_times, peak_amplitudes, peak_stdevs = pulsefish_peaks(fish)
573 # time axis for single pulse:
574 min_time_inx = np.argmin(peak_times)
575 max_time_inx = np.argmax(peak_times)
576 tmax = max(np.abs(peak_times[min_time_inx]-4.0*peak_stdevs[min_time_inx]),
577 np.abs(peak_times[max_time_inx]+4.0*peak_stdevs[max_time_inx]))
578 x = np.arange(-tmax, tmax, 1.0/samplerate)
579 pulse_duration = x[-1] - x[0]
581 # generate a single pulse:
582 pulse = np.zeros(len(x))
583 for time, ampl, std in zip(peak_times, peak_amplitudes, peak_stdevs):
584 pulse += ampl * np.exp(-0.5*((x-time)/std)**2)
585 poffs = len(pulse)//2
587 # paste the pulse into the noise floor:
588 time = np.arange(0, duration, 1.0/samplerate)
589 data = np.random.randn(len(time)) * noise_std
590 period = 1.0/frequency
591 jitter_std = period * jitter_cv
592 if first_pulse is None:
593 first_pulse = np.max([pulse_duration, 3.0*jitter_std])
594 pulse_times = np.arange(first_pulse, duration, period )
595 pulse_times += jitter_std*np.random.randn(len(pulse_times))
596 pulse_indices = np.round(pulse_times * samplerate).astype(int)
597 for inx in pulse_indices[(pulse_indices>=poffs)&(pulse_indices-poffs+len(pulse)<len(data))]:
598 data[inx-poffs:inx-poffs+len(pulse)] += pulse
599 return data
602def normalize_pulsefish(fish):
603 """Normalize times and stdevs of pulse-type EOD waveform.
605 The positions and amplitudes of Gaussian peaks are adjusted such
606 that the resulting EOD waveform has a maximum peak amplitude of one
607 and has the largest peak at time zero.
609 Parameters
610 ----------
611 fish: string, dict or tuple of floats/lists/arrays
612 Specify positions, amplitudes and standard deviations Gaussians peaks that are
613 superimposed to simulate EOD waveforms of pulse-type electric fishes.
614 If string then take positions, amplitudes and standard deviations
615 from the `pulsefish_eodpeaks` dictionary.
616 If dictionary then take pulse properties from the 'times', 'amlitudes'
617 and 'stdevs' keys.
618 If tuple then the first element is the list of peak positions,
619 the second is the list of corresponding amplitudes, and
620 the third one the list of corresponding standard deviations.
622 Returns
623 -------
624 fish: dict
625 Dictionary with adjusted times and standard deviations.
626 """
627 # get peak properties:
628 peak_times, peak_amplitudes, peak_stdevs = pulsefish_peaks(fish)
629 # generate waveform:
630 eodf = 10.0
631 rate = 100000.0
632 first_pulse = 0.5/eodf
633 data = pulsefish_eods(fish, eodf, rate, 1.0/eodf, noise_std=0.0,
634 jitter_cv=0.0, first_pulse=first_pulse)
635 # maximum peak:
636 idx = np.argmax(np.abs(data))
637 # normalize amplitudes:
638 ampl = data[idx]
639 newamplitudes = np.array(peak_amplitudes)/ampl
640 # shift times:
641 deltat = idx/rate - first_pulse
642 newtimes = np.array(peak_times) - deltat
643 # store and return:
644 peaks = dict(times=newtimes,
645 amplitudes=newamplitudes,
646 stdevs=peak_stdevs)
647 return peaks
650def export_pulsefish(fish, name='Unknown_peaks', file=None):
651 """Serialize pulsefish parameter to python code.
653 Add output to the pulsefish_eodpeaks dictionary!
655 Parameters
656 ----------
657 fish: string, dict or tuple of floats/lists/arrays
658 Specify positions, amplitudes and standard deviations Gaussians peaks that are
659 superimposed to simulate EOD waveforms of pulse-type electric fishes.
660 If string then take positions, amplitudes and standard deviations
661 from the `pulsefish_eodpeaks` dictionary.
662 If dictionary then take pulse properties from the 'times', 'amlitudes'
663 and 'stdevs' keys.
664 If tuple then the first element is the list of peak positions,
665 the second is the list of corresponding amplitudes, and
666 the third one the list of corresponding standard deviations.
667 name: string
668 Name of the dictionary to be written.
669 file: string or file or None
670 File name or open file object where to write pulsefish dictionary.
672 Returns
673 -------
674 fish: dict
675 Dictionary with peak times, amplitudes and standard deviations.
676 """
677 # get peak properties:
678 peak_times, peak_amplitudes, peak_stdevs = pulsefish_peaks(fish)
679 # write out dictionary:
680 if file is None:
681 file = sys.stdout
682 try:
683 file.write('')
684 closeit = False
685 except AttributeError:
686 file = open(file, 'w')
687 closeit = True
688 n = 6
689 file.write(name + ' = \\\n')
690 file.write(' dict(times=(')
691 file.write(', '.join([f'{a:.5g}' for a in peak_times[:n]]))
692 for k in range(n, len(peak_times), n):
693 file.write(',\n')
694 file.write(' ' * (9+12))
695 file.write(', '.join([f'{a:.5g}' for a in peak_times[k:k+n]]))
696 if len(peak_times) == 1:
697 file.write(',')
698 file.write('),\n')
699 file.write(' ' * 9 + 'amplitudes=(')
700 file.write(', '.join([f'{p:.5g}' for p in peak_amplitudes[:n]]))
701 for k in range(n, len(peak_amplitudes), n):
702 file.write(',\n')
703 file.write(' ' * (9+8))
704 file.write(', '.join([f'{p:.5g}' for p in peak_amplitudes[k:k+n]]))
705 if len(peak_amplitudes) == 1:
706 file.write(',')
707 file.write('),\n')
708 file.write(' ' * 9 + 'stdevs=(')
709 file.write(', '.join([f'{p:.5g}' for p in peak_stdevs[:n]]))
710 for k in range(n, len(peak_stdevs), n):
711 file.write(',\n')
712 file.write(' ' * (9+8))
713 file.write(', '.join([f'{p:.5g}' for p in peak_stdevs[k:k+n]]))
714 if len(peak_stdevs) == 1:
715 file.write(',')
716 file.write('))\n')
717 if closeit:
718 file.close()
719 # return dictionary:
720 peaks = dict(times=peak_times,
721 amplitudes=peak_amplitudes,
722 stdevs=peak_stdevs)
723 return peaks
726def generate_waveform(filename):
727 """Interactively generate audio file with simulated EOD waveforms.
729 Parameters needed to generate EOD waveforms are take from console input.
731 Parameters
732 ----------
733 filename: string
734 Name of file inclusively extension (e.g. '.wav')
735 used to store the simulated EOD waveforms.
736 """
737 import os
738 from audioio import write_audio
739 from thunderlab.consoleinput import read, select, save_inputs
740 # generate file:
741 samplerate = read('Sampling rate in Hz', '44100', float, 1.0)
742 duration = read('Duration in seconds', '10', float, 0.001)
743 nfish = read('Number of fish', '1', int, 1)
744 ndata = read('Number of electrodes', '1', int, 1)
745 fish_spread = 1
746 if ndata > 1:
747 fish_spread = read('Number of electrodes fish are spread over', '2', int, 1)
748 data = np.random.randn(int(duration*samplerate), ndata)*0.01
749 fish_indices = np.random.randint(ndata, size=nfish)
750 eodt = 'a'
751 eodf = 800.0
752 eoda = 1.0
753 eodsig = 'n'
754 pulse_jitter = 0.1
755 chirp_freq = 5.0
756 chirp_size = 100.0
757 chirp_width = 0.01
758 chirp_kurtosis = 1.0
759 rise_freq = 0.1
760 rise_size = 10.0
761 rise_tau = 1.0
762 rise_decay_tau = 10.0
763 for k in range(nfish):
764 print('')
765 fish = 'Fish %d: ' % (k+1)
766 eodt = select(fish + 'EOD type', eodt, ['a', 'e', '1', '2', '3'],
767 ['Apteronotus', 'Eigenmannia',
768 'monophasic pulse', 'biphasic pulse', 'triphasic pulse'])
769 eodf = read(fish + 'EOD frequency in Hz', '%g'%eodf, float, 1.0, 3000.0)
770 eoda = read(fish + 'EOD amplitude', '%g'%eoda, float, 0.0, 10.0)
771 if eodt in 'ae':
772 eodsig = select(fish + 'Add communication signals', eodsig, ['n', 'c', 'r'],
773 ['fixed EOD', 'chirps', 'rises'])
774 eodfreq = eodf
775 if eodsig == 'c':
776 chirp_freq = read('Number of chirps per second', '%g'%chirp_freq, float, 0.001)
777 chirp_size = read('Size of chirp in Hz', '%g'%chirp_size, float, 1.0)
778 chirp_width = 0.001*read('Width of chirp in ms', '%g'%(1000.0*chirp_width), float, 1.0)
779 eodfreq, _ = chirps(eodf, samplerate, duration,
780 chirp_freq, chirp_size, chirp_width, chirp_kurtosis)
781 elif eodsig == 'r':
782 rise_freq = read('Number of rises per second', '%g'%rise_freq, float, 0.00001)
783 rise_size = read('Size of rise in Hz', '%g'%rise_size, float, 0.01)
784 rise_tau = read('Time-constant of rise onset in seconds', '%g'%rise_tau, float, 0.01)
785 rise_decay_tau = read('Time-constant of rise decay in seconds', '%g'%rise_decay_tau, float, 0.01)
786 eodfreq = rises(eodf, samplerate, duration,
787 rise_freq, rise_size, rise_tau, rise_decay_tau)
788 if eodt == 'a':
789 fishdata = eoda*wavefish_eods('Alepto', eodfreq, samplerate, duration,
790 phase0=0.0, noise_std=0.0)
791 elif eodt == 'e':
792 fishdata = eoda*wavefish_eods('Eigenmannia', eodfreq, samplerate,
793 duration, phase0=0.0, noise_std=0.0)
794 else:
795 pulse_jitter = read(fish + 'CV of pulse jitter', '%g'%pulse_jitter, float, 0.0, 2.0)
796 if eodt == '1':
797 fishdata = eoda*pulsefish_eods('Monophasic', eodf, samplerate, duration,
798 jitter_cv=pulse_jitter, noise_std=0.0)
799 elif eodt == '2':
800 fishdata = eoda*pulsefish_eods('Biphasic', eodf, samplerate, duration,
801 jitter_cv=pulse_jitter, noise_std=0.0)
802 elif eodt == '3':
803 fishdata = eoda*pulsefish_eods('Triphasic', eodf, samplerate, duration,
804 jitter_cv=pulse_jitter, noise_std=0.0)
805 i = fish_indices[k]
806 for j in range(fish_spread):
807 data[:, (i+j)%ndata] += fishdata*(0.2**j)
809 maxdata = np.max(np.abs(data))
810 write_audio(filename, 0.9*data/maxdata, samplerate)
811 input_file = os.path.splitext(filename)[0] + '.inp'
812 save_inputs(input_file)
813 print(f'\nWrote fakefish data to file "{filename}".')
816def demo():
817 import matplotlib.pyplot as plt
818 samplerate = 40000.0 # in Hz
819 duration = 10.0 # in sec
821 inset_len = 0.01 # in sec
822 inset_indices = int(inset_len*samplerate)
823 ws_fac = 0.1 # whitespace factor or ylim (between 0. and 1.)
825 # generate data:
826 eodf = 400.0
827 wavefish = wavefish_eods('Alepto', eodf, samplerate, duration, noise_std=0.02)
828 eodf = 650.0
829 wavefish += 0.5*wavefish_eods('Eigenmannia', eodf, samplerate, duration)
831 pulsefish = pulsefish_eods('Biphasic', 80.0, samplerate, duration,
832 noise_std=0.02, jitter_cv=0.1, first_pulse=inset_len/2)
833 time = np.arange(len(wavefish))/samplerate
835 fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(19, 10))
837 # get proper wavefish ylim
838 ymin = np.min(wavefish)
839 ymax = np.max(wavefish)
840 dy = ws_fac*(ymax - ymin)
841 ymin -= dy
842 ymax += dy
844 # complete wavefish:
845 ax[0][0].set_title('Wavefish')
846 ax[0][0].set_ylim(ymin, ymax)
847 ax[0][0].plot(time, wavefish)
849 # wavefish zoom in:
850 ax[0][1].set_title('Wavefish ZOOM IN')
851 ax[0][1].set_ylim(ymin, ymax)
852 ax[0][1].plot(time[:inset_indices], wavefish[:inset_indices], '-o')
854 # get proper pulsefish ylim
855 ymin = np.min(pulsefish)
856 ymax = np.max(pulsefish)
857 dy = ws_fac*(ymax - ymin)
858 ymin -= dy
859 ymax += dy
861 # complete pulsefish:
862 ax[1][0].set_title('Pulsefish')
863 ax[1][0].set_ylim(ymin, ymax)
864 ax[1][0].plot(time, pulsefish)
866 # pulsefish zoom in:
867 ax[1][1].set_title('Pulsefish ZOOM IN')
868 ax[1][1].set_ylim(ymin, ymax)
869 ax[1][1].plot(time[:inset_indices], pulsefish[:inset_indices], '-o')
871 for row in ax:
872 for c_ax in row:
873 c_ax.set_xlabel('Time [sec]')
874 c_ax.set_ylabel('Amplitude')
876 plt.tight_layout()
878 # chirps:
879 chirps_freq = chirps(600.0, samplerate, duration)
880 chirps_data = wavefish_eods('Alepto', chirps_freq, samplerate)
882 # rises:
883 rises_freq = rises(600.0, samplerate, duration, rise_size=20.0)
884 rises_data = wavefish_eods('Alepto', rises_freq, samplerate)
886 nfft = 256
887 fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(19, 10))
888 ax[0].set_title('Chirps')
889 ax[0].specgram(chirps_data, Fs=samplerate, NFFT=nfft, noverlap=nfft//16)
890 time = np.arange(len(chirps_freq))/samplerate
891 ax[0].plot(time[:-nfft//2], chirps_freq[nfft//2:], '-k', lw=2)
892 ax[0].set_ylim(0.0, 3000.0)
893 ax[0].set_ylabel('Frequency [Hz]')
895 nfft = 4096
896 ax[1].set_title('Rises')
897 ax[1].specgram(rises_data, Fs=samplerate, NFFT=nfft, noverlap=nfft//2)
898 time = np.arange(len(rises_freq))/samplerate
899 ax[1].plot(time[:-nfft//4], rises_freq[nfft//4:], '-k', lw=2)
900 ax[1].set_ylim(500.0, 700.0)
901 ax[1].set_ylabel('Frequency [Hz]')
902 ax[1].set_xlabel('Time [s]')
903 plt.tight_layout()
905 plt.show()
908def main(args=[]):
909 from .version import __year__
911 if len(args) > 0:
912 if len(args) == 1 or args[0] != '-s':
913 print('usage: fakefish [-h|--help] [-s audiofile]')
914 print('')
915 print('Without arguments, run a demo for illustrating fakefish functionality.')
916 print('')
917 print('-s audiofile: writes audiofile with user defined simulated electric fishes.')
918 print('')
919 print(f'by bendalab ({__year__})')
920 else:
921 generate_waveform(args[1])
922 else:
923 demo()
926if __name__ == '__main__':
927 import sys
928 main(sys.argv[1:])