Coverage for src / thunderfish / fakefish.py: 96%
530 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-02-09 14:25 +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# EODs of various wavefish species:
93Sine_harmonics = dict(amplitudes=(1.0,), phases=(0.5*np.pi,),
94 species='Sinewave')
96Apteronotus_leptorhynchus_harmonics = \
97 dict(amplitudes=(0.90062, 0.15311, 0.072049, 0.012609, 0.011708),
98 phases=(1.3623, 2.3246, 0.9869, 2.6492, -2.6885),
99 species='Apteronotus leptorhynchus')
101Apteronotus_rostratus_harmonics = \
102 dict(amplitudes=(0.64707, 0.43874, 0.063592, 0.07379, 0.040199, 0.023073,
103 0.0097678),
104 phases=(2.2988, 0.78876, -1.316, 2.2416, 2.0413, 1.1022,
105 -2.0513),
106 species='Apteronotus rostratus')
108Eigenmannia_harmonics = \
109 dict(amplitudes=(1.0087, 0.23201, 0.060524, 0.020175, 0.010087, 0.0080699),
110 phases=(1.3414, 1.3228, 2.9242, 2.8157, 2.6871, -2.8415),
111 species='Eigenmannia')
113Sternarchella_terminalis_harmonics = \
114 dict(amplitudes=(0.11457, 0.4401, 0.41055, 0.20132, 0.061364, 0.011389,
115 0.0057985),
116 phases=(-2.7106, 2.4472, 1.6829, 0.79085, 0.119, -0.82355,
117 -1.9956),
118 species='Sternarchella terminalis')
120Sternopygus_dariensis_harmonics = \
121 dict(amplitudes=(0.98843, 0.41228, 0.047848, 0.11048, 0.022801, 0.030706,
122 0.019018),
123 phases=(1.4153, 1.3141, 3.1062, -2.3961, -1.9524, 0.54321,
124 1.6844),
125 species='Sternopygus dariensis')
127wavefish_harmonics = dict(Sine=Sine_harmonics,
128 Alepto=Apteronotus_leptorhynchus_harmonics,
129 Arostratus=Apteronotus_rostratus_harmonics,
130 Eigenmannia=Eigenmannia_harmonics,
131 Sternarchella=Sternarchella_terminalis_harmonics,
132 Sternopygus=Sternopygus_dariensis_harmonics)
133"""Amplitudes and phases of EOD waveforms of various species of wave-type electric fish."""
136def wavefish_spectrum(fish):
137 """Amplitudes and phases of a wavefish EOD.
139 Parameters
140 ----------
141 fish: string, dict or tuple of lists/arrays or 2-D array
142 Specify relative amplitudes and phases of the fundamental and its harmonics.
143 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary.
144 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys.
145 If tuple then the first element is the list of amplitudes and
146 the second one the list of relative phases in radians.
147 If 2-D array, as returned from waveanalysis.analyse_wave(),
148 then take relative amplitudes from third column and phases
149 from fifth colum.
151 Returns
152 -------
153 amplitudes: array of floats
154 Amplitudes of the fundamental and its harmonics.
155 phases: array of floats
156 Phases in radians of the fundamental and its harmonics.
158 Raises
159 ------
160 KeyError:
161 Unknown fish.
162 IndexError:
163 Amplitudes and phases differ in length.
164 """
165 if isinstance(fish, np.ndarray) and fish.ndim == 2:
166 amplitudes = fish[:, 3]
167 phases = fish[:, 5]
168 elif isinstance(fish, (tuple, list)):
169 amplitudes = fish[0]
170 phases = fish[1]
171 elif isinstance(fish, dict):
172 amplitudes = fish['amplitudes']
173 phases = fish['phases']
174 else:
175 if not fish in wavefish_harmonics:
176 raise KeyError('unknown wavefish. Choose one of ' +
177 ', '.join(wavefish_harmonics.keys()))
178 amplitudes = wavefish_harmonics[fish]['amplitudes']
179 phases = wavefish_harmonics[fish]['phases']
180 if len(amplitudes) != len(phases):
181 raise IndexError('need exactly as many phases as amplitudes')
182 # remove NaNs:
183 for k in reversed(range(len(amplitudes))):
184 if np.isfinite(amplitudes[k]) or np.isfinite(phases[k]):
185 amplitudes = amplitudes[:k+1]
186 phases = phases[:k+1]
187 break
188 return amplitudes, phases
191def wavefish_eods(fish='Eigenmannia', frequency=100.0, rate=44100.0,
192 duration=1.0, phase0=0.0, noise_std=0.05):
193 """Simulate EOD waveform of a wave-type fish.
195 The waveform is constructed by superimposing sinewaves of integral
196 multiples of the fundamental frequency - the fundamental and its
197 harmonics. The fundamental frequency of the EOD (EODf) is given by
198 `frequency`. With `fish` relative amplitudes and phases of the
199 fundamental and its harmonics are specified.
201 The generated waveform is `duration` seconds long and is sampled with
202 `rate` Hertz. Gaussian white noise with a standard deviation of
203 `noise_std` is added to the generated waveform.
205 Parameters
206 ----------
207 fish: string, dict or tuple of lists/arrays or 2-D array
208 Specify relative amplitudes and phases of the fundamental and its harmonics.
209 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary.
210 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys.
211 If tuple then the first element is the list of amplitudes and
212 the second one the list of relative phases in radians.
213 If 2-D array, as returned from waveanalysis.analyse_wave(),
214 then take relative amplitudes from third column and phases
215 from fifth colum.
216 frequency: float or array of floats
217 EOD frequency of the fish in Hertz. Either fixed number or array for
218 time-dependent frequencies.
219 rate: float
220 Sampling rate in Hertz.
221 duration: float
222 Duration of the generated data in seconds. Only used if frequency is scalar.
223 phase0: float
224 Phase offset of the EOD waveform in radians.
225 noise_std: float
226 Standard deviation of additive Gaussian white noise.
228 Returns
229 -------
230 data: array of floats
231 Generated data of a wave-type fish.
233 Raises
234 ------
235 KeyError:
236 Unknown fish.
237 IndexError:
238 Amplitudes and phases differ in length.
239 """
240 # get relative amplitude and phases:
241 amplitudes, phases = wavefish_spectrum(fish)
242 # compute phase:
243 if np.isscalar(frequency):
244 phase = np.arange(0, duration, 1.0/rate)
245 phase *= frequency
246 else:
247 phase = np.cumsum(frequency)/rate
248 # generate EOD:
249 data = np.zeros(len(phase))
250 for har, (ampl, phi) in enumerate(zip(amplitudes, phases)):
251 if np.isfinite(ampl) and np.isfinite(phi):
252 data += ampl * np.sin(2*np.pi*(har+1)*phase + phi - (har+1)*phase0)
253 # add noise:
254 data += noise_std * np.random.randn(len(data))
255 return data
258def normalize_wavefish(fish, mode='peak'):
259 """Normalize amplitudes and phases of wave-type EOD waveform.
261 The amplitudes and phases of the Fourier components are adjusted
262 such that the resulting EOD waveform has a peak-to-peak amplitude
263 of two and the peak of the waveform is at time zero (mode is set
264 to 'peak') or that the fundamental has an amplitude of one and a
265 phase of 0 (mode is set to 'zero').
267 Parameters
268 ----------
269 fish: string, dict or tuple of lists/arrays or 2-D array
270 Specify relative amplitudes and phases of the fundamental and its harmonics.
271 If string then take amplitudes and phases from the `wavefish_harmonics` dictionary.
272 If dictionary then take amplitudes and phases from the 'amlitudes' and 'phases' keys.
273 If tuple then the first element is the list of amplitudes and
274 the second one the list of relative phases in radians.
275 If 2-D array, as returned from waveanalysis.analyse_wave(),
276 then take relative amplitudes from third column and phases
277 from fifth colum.
278 mode: 'peak' or 'zero'
279 How to normalize amplitude and phases:
280 - 'peak': normalize waveform to a peak-to-peak amplitude of two
281 and shift it such that its peak is at time zero.
282 - 'zero': scale amplitude of fundamental to one and its phase to zero.
284 Returns
285 -------
286 amplitudes: array of floats
287 Adjusted amplitudes of the fundamental and its harmonics.
288 phases: array of floats
289 Adjusted phases in radians of the fundamental and its harmonics.
291 """
292 # get relative amplitude and phases:
293 amplitudes, phases = wavefish_spectrum(fish)
294 if mode == 'zero':
295 newamplitudes = np.array(amplitudes)/amplitudes[0]
296 newphases = np.array([p+(k+1)*(-phases[0]) for k, p in enumerate(phases)])
297 newphases %= 2.0*np.pi
298 newphases[newphases>np.pi] -= 2.0*np.pi
299 return newamplitudes, newphases
300 else:
301 # generate waveform:
302 eodf = 100.0
303 rate = 100000.0
304 data = wavefish_eods(fish, eodf, rate, 2.0/eodf, noise_std=0.0)
305 # normalize amplitudes:
306 ampl = 0.5*(np.max(data) - np.min(data))
307 newamplitudes = np.array(amplitudes)/ampl
308 # shift phases:
309 deltat = np.argmax(data[:int(rate/eodf)])/rate
310 deltap = 2.0*np.pi*deltat*eodf
311 newphases = np.array([p+(k+1)*deltap for k, p in enumerate(phases)])
312 newphases %= 2.0*np.pi
313 newphases[newphases>np.pi] -= 2.0*np.pi
314 # return:
315 return newamplitudes, newphases
318def export_wavefish(fish, name='Unknown_harmonics', species='', file=None):
319 """Serialize wavefish parameter to python code.
321 Add output to the wavefish_harmonics dictionary!
323 Parameters
324 ----------
325 fish: string, dict or tuple of lists/arrays or 2-D array
326 Specify relative amplitudes and phases of the fundamental and its harmonics and optional species name.
327 If string then take amplitudes, phases, and species from the `wavefish_harmonics` dictionary.
328 If dictionary then take amplitudes, phases, and species from the 'amlitudes', 'phases', and 'species' keys.
329 If tuple then the first element is the list of amplitudes and
330 the second one the list of relative phases in radians.
331 If 2-D array, as returned from waveanalysis.analyse_wave(),
332 then take relative amplitudes from third column and phases
333 from fifth colum.
334 name: str
335 Name of the dictionary to be written. If empty take species name.
336 species: str
337 Name of the fish species.
338 file: str or Path or file or None
339 File path or open file object where to write wavefish dictionary.
341 Returns
342 -------
343 fish: dict
344 Dictionary with amplitudes, phases, and species.
345 """
346 if fish is None or (isinstance(fish, dict) and len(fish) == 0):
347 return {}
348 # get relative amplitude and phases:
349 amplitudes, phases = wavefish_spectrum(fish)
350 harmonics = dict(amplitudes=amplitudes,
351 phases=phases)
352 # species:
353 if not species:
354 if isinstance(fish, dict):
355 species = fish.get('species', '')
356 elif isinstance(fish, str) and fish in wavefish_harmonics:
357 species = wavefish_harmonics[fish]['species']
358 # name:
359 if not name:
360 name = species
361 name = name.replace(' ', '_')
362 # write out dictionary:
363 if file is None:
364 file = sys.stdout
365 try:
366 file.write('')
367 closeit = False
368 except AttributeError:
369 file = open(file, 'w')
370 closeit = True
371 n = 6
372 file.write(name + ' = \\\n')
373 file.write(' dict(amplitudes=(')
374 file.write(', '.join([f'{a:.5g}' for a in amplitudes[:n]]))
375 for k in range(n, len(amplitudes), n):
376 file.write(',\n')
377 file.write(' ' * (9 + 12))
378 file.write(', '.join([f'{a:.5g}' for a in amplitudes[k:k + n]]))
379 file.write('),\n')
380 file.write(' ' * 9 + 'phases=(')
381 file.write(', '.join([f'{p:.5g}' for p in phases[:n]]))
382 for k in range(n, len(phases), n):
383 file.write(',\n')
384 file.write(' ' * (9 + 8))
385 file.write(', '.join([f'{p:.5g}' for p in phases[k:k + n]]))
386 if species:
387 file.write('),\n')
388 file.write(' ' * 9 + f'species=\'{species}\')\n')
389 harmonics['species'] = species
390 else:
391 file.write('))\n')
392 if closeit:
393 file.close()
394 return harmonics
397def chirps(eodf=100.0, rate=44100.0, duration=1.0, chirp_freq=5.0,
398 chirp_size=100.0, chirp_width=0.01, chirp_kurtosis=1.0, chirp_contrast=0.05):
399 """Simulate frequency trace with chirps.
401 A chirp is modeled as a Gaussian frequency modulation.
402 The first chirp is placed at 0.5/chirp_freq.
404 Parameters
405 ----------
406 eodf: float
407 EOD frequency of the fish in Hertz.
408 rate: float
409 Sampling rate in Hertz.
410 duration: float
411 Duration of the generated data in seconds.
412 chirp_freq: float
413 Frequency of occurance of chirps in Hertz.
414 chirp_size: float
415 Size of the chirp (maximum frequency increase above eodf) in Hertz.
416 chirp_width: float
417 Width of the chirp at 10% height in seconds.
418 chirp_kurtosis: float:
419 Shape of the chirp. =1: Gaussian, >1: more rectangular, <1: more peaked.
420 chirp_contrast: float
421 Maximum amplitude reduction of EOD during chirp.
423 Returns
424 -------
425 frequency: array of floats
426 Generated frequency trace that can be passed on to wavefish_eods().
427 amplitude: array of floats
428 Generated amplitude modulation that can be used to multiply the trace generated by
429 wavefish_eods().
430 """
431 # baseline eod frequency and amplitude modulation:
432 n = len(np.arange(0, duration, 1.0/rate))
433 frequency = eodf * np.ones(n)
434 am = np.ones(n)
435 # time points for chirps:
436 chirp_period = 1.0/chirp_freq
437 chirp_times = np.arange(0.5*chirp_period, duration, chirp_period)
438 # chirp frequency waveform:
439 chirp_t = np.arange(-2.0*chirp_width, 2.0*chirp_width, 1./rate)
440 chirp_sig = 0.5*chirp_width / (2.0*np.log(10.0))**(0.5/chirp_kurtosis)
441 gauss = np.exp(-0.5*((chirp_t/chirp_sig)**2.0)**chirp_kurtosis)
442 # add chirps on baseline eodf:
443 for ct in chirp_times:
444 index = int(ct*rate)
445 i0 = index - len(gauss)//2
446 i1 = i0 + len(gauss)
447 gi0 = 0
448 gi1 = len(gauss)
449 if i0 < 0:
450 gi0 -= i0
451 i0 = 0
452 if i1 >= len(frequency):
453 gi1 -= i1 - len(frequency)
454 i1 = len(frequency)
455 frequency[i0:i1] += chirp_size * gauss[gi0:gi1]
456 am[i0:i1] -= chirp_contrast * gauss[gi0:gi1]
457 return frequency, am
460def rises(eodf=100.0, rate=44100.0, duration=1.0, rise_freq=0.1,
461 rise_size=10.0, rise_tau=1.0, decay_tau=10.0):
462 """Simulate frequency trace with rises.
464 A rise is modeled as a double exponential frequency modulation.
466 Parameters
467 ----------
468 eodf: float
469 EOD frequency of the fish in Hertz.
470 rate: float
471 Sampling rate in Hertz.
472 duration: float
473 Duration of the generated data in seconds.
474 rise_freq: float
475 Frequency of occurance of rises in Hertz.
476 rise_size: float
477 Size of the rise (frequency increase above eodf) in Hertz.
478 rise_tau: float
479 Time constant of the frequency increase of the rise in seconds.
480 decay_tau: float
481 Time constant of the frequency decay of the rise in seconds.
483 Returns
484 -------
485 data: array of floats
486 Generated frequency trace that can be passed on to wavefish_eods().
487 """
488 n = len(np.arange(0, duration, 1.0/rate))
489 # baseline eod frequency:
490 frequency = eodf * np.ones(n)
491 # time points for rises:
492 rise_period = 1.0/rise_freq
493 rise_times = np.arange(0.5*rise_period, duration, rise_period)
494 # rise frequency waveform:
495 rise_t = np.arange(0.0, 5.0*decay_tau, 1./rate)
496 rise = rise_size * (1.0-np.exp(-rise_t/rise_tau)) * np.exp(-rise_t/decay_tau)
497 # add rises on baseline eodf:
498 for r in rise_times:
499 index = int(r*rate)
500 if index+len(rise) > len(frequency):
501 rise_index = len(frequency)-index
502 frequency[index:index+rise_index] += rise[:rise_index]
503 break
504 else:
505 frequency[index:index+len(rise)] += rise
506 return frequency
509# EODs of various pulsefish species:
511Monophasic_phases = \
512 dict(times=(0,),
513 amplitudes=(1,),
514 stdevs=(0.0003,),
515 species='Biphasic')
517Biphasic_phases = \
518 dict(times=(9e-05, 0.00049),
519 amplitudes=(1.1922, -0.95374),
520 stdevs=(0.0003, 0.00025),
521 species='Biphasic')
523Triphasic_phases = \
524 dict(times=(3e-05, 0.00018, 0.00043),
525 amplitudes=(1.2382, -0.9906, 0.12382),
526 stdevs=(0.0001, 0.0001, 0.0002),
527 species='Triphasic')
529pulsefish_eodphases = dict(Monophasic=Monophasic_phases,
530 Biphasic=Biphasic_phases,
531 Triphasic=Triphasic_phases)
532"""Positions, amplitudes, and standard deviations of Gaussians that
533 make up the EOD phases of pulse-type electric fish.
534"""
537def pulsefish_phases(fish):
538 """Position, amplitudes and standard deviations of EOD phases in pulsefish waveforms.
540 Parameters
541 ----------
542 fish: string, dict or tuple of floats/lists/arrays
543 Specify positions, amplitudes, and standard deviations
544 of Gaussians phases that are superimposed to simulate EOD waveforms
545 of pulse-type electric fishes.
546 If string then take positions, amplitudes and standard deviations
547 from the `pulsefish_eodphases` dictionary.
548 If dictionary then take pulse properties from the 'times', 'amlitudes'
549 and 'stdevs' keys.
550 If tuple then the first element is the list of phase positions,
551 the second is the list of corresponding amplitudes, and
552 the third one the list of corresponding standard deviations.
554 Returns
555 -------
556 times : array of floats
557 Positions of the phases.
558 amplitudes : array of floats
559 Amplitudes of the phases.
560 stdevs : array of floats
561 Standard deviations of the phases.
563 Raises
564 ------
565 KeyError:
566 Unknown fish.
567 IndexError:
568 Phase positions, amplitudes, or standard deviations differ in length.
569 """
570 if isinstance(fish, (tuple, list)):
571 times = fish[0]
572 amplitudes = fish[1]
573 stdevs = fish[2]
574 elif isinstance(fish, dict):
575 times = fish['times']
576 amplitudes = fish['amplitudes']
577 stdevs = fish['stdevs']
578 else:
579 if not fish in pulsefish_eodphases:
580 raise KeyError('unknown pulse-type fish. Choose one of ' +
581 ', '.join(pulsefish_eodphases.keys()))
582 phases = pulsefish_eodphases[fish]
583 times = phases['times']
584 amplitudes = phases['amplitudes']
585 stdevs = phases['stdevs']
586 if len(stdevs) != len(amplitudes) or len(stdevs) != len(times):
587 raise IndexError('need exactly as many standard deviations as amplitudes and times')
588 return times, amplitudes, stdevs
591def pulsefish_spectrum(fish, freqs):
592 """Analytically computed single-pulse spectrum.
594 The spectrum is the sum of the spectra of all of the phases. Each
595 phase is the spectrum of a Gaussian (which again is a Gaussian
596 centered at f=0) multiplied with a shift term to account for its
597 position relative to the first phase.
599 You get the amplitude spectral density by taking the absolute value.
600 The energy spectral density is the amplitude spectral density squared:
601 ```
602 spec = pulsefish_spectrum('Biphasic', freqs)
603 ampl = np.abs(spec)
604 energy = np.abs(spec)**2
605 ```
607 Parameters
608 ----------
609 fish: string, dict or tuple of floats/lists/arrays
610 Specify positions, amplitudes, and standard deviations
611 of Gaussians phases that are superimposed to simulate EOD waveforms
612 of pulse-type electric fishes.
613 If string then take positions, amplitudes and standard deviations
614 from the `pulsefish_eodphases` dictionary.
615 If dictionary then take pulse properties from the 'times', 'amlitudes'
616 and 'stdevs' keys.
617 If tuple then the first element is the list of phase positions,
618 the second is the list of corresponding amplitudes, and
619 the third one the list of corresponding standard deviations.
620 freqs: 1-D array of float
621 Frequencies at which the spectrum is computed.
623 Returns
624 -------
625 spectrum: 1-D array of complex
626 The one-sided complex-valued spectrum of the single pulse EOD.
627 The squared magnitude (the energy spectrum) has
628 the unit (x s)^2 = x^2 s/Hz.
630 """
631 times, ampls, stdevs = pulsefish_phases(fish)
632 spec = np.zeros(len(freqs), dtype=complex)
633 for dt, a, s in zip(times, ampls, stdevs):
634 gauss = a*np.sqrt(2*np.pi)*s*np.exp(-0.5*(2*np.pi*s*freqs)**2)
635 shift = np.exp(-2j*np.pi*freqs*dt)
636 spec += gauss*shift
637 spec *= np.sqrt(2) # because of one-sided spectrum
638 return spec
641def pulsefish_waveform(fish, t):
642 """Compute a single pulsefish EOD.
644 You may use pulsefish_parameter() to pass the paramters returned
645 by pulsefish_phases() to this function.
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 t: array of float
661 The time array over which the pulse waveform is evaluated.
663 Returns
664 -------
665 pulse: array of float
666 The pulse waveform for the times given in `t`.
668 """
669 times, ampls, stdevs = pulsefish_phases(fish)
670 pulse = np.zeros(len(t))
671 for dt, a, s in zip(times, ampls, stdevs):
672 pulse += a*np.exp(-0.5*((t - dt)/s)**2)
673 return pulse
676def pulsefish_eods(fish='Biphasic', frequency=100.0, rate=44100.0,
677 duration=1.0, noise_std=0.01, jitter_cv=0.1,
678 first_pulse=None):
679 """Simulate EODs of a pulse-type fish.
681 Pulses are spaced by 1/frequency, jittered as determined by
682 jitter_cv. Each pulse is a combination of Gaussian phases, whose
683 positions, amplitudes and widths are given by `fish`.
685 The generated waveform is duration seconds long and is sampled
686 with rate Hertz. Gaussian white noise with a standard deviation
687 of noise_std is added to the generated pulse train.
689 Parameters
690 ----------
691 fish: string, dict or tuple of floats/lists/arrays
692 Specify positions, amplitudes, and standard deviations
693 of Gaussians phases that are superimposed to simulate EOD waveforms
694 of pulse-type electric fishes.
695 If string then take positions, amplitudes and standard deviations
696 from the `pulsefish_eodphases` dictionary.
697 If dictionary then take pulse properties from the 'times', 'amlitudes'
698 and 'stdevs' keys.
699 If tuple then the first element is the list of phase positions,
700 the second is the list of corresponding amplitudes, and
701 the third one the list of corresponding standard deviations.
702 frequency: float
703 EOD frequency of the fish in Hz.
704 rate: float
705 Sampling Rate in Hz.
706 duration: float
707 Duration of the generated data in seconds.
708 noise_std: float
709 Standard deviation of additive Gaussian white noise.
710 jitter_cv: float
711 Gaussian distributed jitter of pulse times as coefficient of variation
712 of inter-pulse intervals.
713 first_pulse: float or None
714 The position of the first pulse. If None it is choosen automatically
715 depending on pulse width, jitter, and frequency.
717 Returns
718 -------
719 data: array of floats
720 Generated data of a pulse-type fish.
722 Raises
723 ------
724 KeyError:
725 Unknown fish.
726 IndexError:
727 Phase positions, amplitudes, or standard deviations differ in length.
729 """
730 # get phase properties:
731 phase_times, phase_amplitudes, phase_stdevs = pulsefish_phases(fish)
732 # time axis for single pulse:
733 min_time_inx = np.argmin(phase_times)
734 max_time_inx = np.argmax(phase_times)
735 tmax = max(np.abs(phase_times[min_time_inx]-4.0*phase_stdevs[min_time_inx]),
736 np.abs(phase_times[max_time_inx]+4.0*phase_stdevs[max_time_inx]))
737 x = np.arange(-tmax, tmax, 1.0/rate)
738 pulse_duration = x[-1] - x[0]
740 # generate a single pulse:
741 pulse = pulsefish_waveform((phase_times, phase_amplitudes, phase_stdevs), x)
742 #pulse = np.zeros(len(x))
743 #for time, ampl, std in zip(phase_times, phase_amplitudes, phase_stdevs):
744 # pulse += ampl * np.exp(-0.5*((x-time)/std)**2)
745 poffs = len(pulse)//2
747 # paste the pulse into the noise floor:
748 time = np.arange(0, duration, 1.0/rate)
749 data = np.random.randn(len(time)) * noise_std
750 period = 1.0/frequency
751 jitter_std = period * jitter_cv
752 if first_pulse is None:
753 first_pulse = np.max([pulse_duration, 3.0*jitter_std])
754 pulse_times = np.arange(first_pulse, duration, period )
755 pulse_times += jitter_std*np.random.randn(len(pulse_times))
756 pulse_indices = np.round(pulse_times * rate).astype(int)
757 for inx in pulse_indices[(pulse_indices>=poffs)&(pulse_indices-poffs+len(pulse)<len(data))]:
758 data[inx-poffs:inx-poffs+len(pulse)] += pulse
759 return data
762def normalize_pulsefish(fish):
763 """Normalize times and stdevs of pulse-type EOD waveform.
765 The positions and amplitudes of Gaussian phases are adjusted such
766 that the resulting EOD waveform has a maximum peak amplitude of one
767 and has the largest phase at time zero.
769 Parameters
770 ----------
771 fish: string, dict or tuple of floats/lists/arrays
772 Specify positions, amplitudes, and standard deviations
773 of Gaussians phases that are superimposed to simulate EOD waveforms
774 of pulse-type electric fishes.
775 If string then take positions, amplitudes and standard deviations
776 from the `pulsefish_eodphases` dictionary.
777 If dictionary then take pulse properties from the 'times', 'amlitudes'
778 and 'stdevs' keys.
779 If tuple then the first element is the list of phase positions,
780 the second is the list of corresponding amplitudes, and
781 the third one the list of corresponding standard deviations.
783 Returns
784 -------
785 fish: dict
786 Dictionary with adjusted times and standard deviations.
787 """
788 # get phase properties:
789 phase_times, phase_amplitudes, phase_stdevs = pulsefish_phases(fish)
790 # generate waveform:
791 eodf = 10.0
792 rate = 100000.0
793 first_pulse = 0.5/eodf
794 data = pulsefish_eods(fish, eodf, rate, 1.0/eodf, noise_std=0.0,
795 jitter_cv=0.0, first_pulse=first_pulse)
796 # maximum peak:
797 idx = np.argmax(np.abs(data))
798 # normalize amplitudes:
799 ampl = data[idx]
800 newamplitudes = np.array(phase_amplitudes)/ampl
801 # shift times:
802 deltat = idx/rate - first_pulse
803 newtimes = np.array(phase_times) - deltat
804 # store and return:
805 phases = dict(times=newtimes,
806 amplitudes=newamplitudes,
807 stdevs=phase_stdevs)
808 return phases
811def export_pulsefish(fish, name='Unknown_phases', species='', file=None):
812 """Serialize pulsefish parameter to python code.
814 Add output to the pulsefish_eodphases dictionary!
816 Parameters
817 ----------
818 fish: string, dict or tuple of floats/lists/arrays
819 Specify positions, amplitudes, and standard deviations
820 of Gaussians phases that are superimposed to simulate EOD waveforms
821 of pulse-type electric fishes.
822 If string then take positions, amplitudes, standard deviations, and
823 species from the `pulsefish_eodphases` dictionary.
824 If dictionary then take pulse properties and species from the 'times',
825 'amlitudes', 'stdevs' and 'species' keys.
826 If tuple then the first element is the list of phase positions,
827 the second is the list of corresponding amplitudes, and
828 the third one the list of corresponding standard deviations.
829 name: string
830 Name of the dictionary to be written. If empty take species name.
831 species: str
832 Name of the fish species.
833 file: str or Path or file or None
834 File name or open file object where to write pulsefish dictionary.
836 Returns
837 -------
838 fish: dict
839 Dictionary with phase times, amplitudes, standard deviations and
840 species.
841 """
842 if fish is None or (isinstance(fish, dict) and len(fish) == 0):
843 return {}
844 # get phase properties:
845 phase_times, phase_amplitudes, phase_stdevs = pulsefish_phases(fish)
846 phases = dict(times=phase_times,
847 amplitudes=phase_amplitudes,
848 stdevs=phase_stdevs)
849 # species:
850 if not species:
851 if isinstance(fish, dict):
852 species = fish.get('species', '')
853 elif isinstance(fish, str) and fish in pulsefish_eodphases:
854 species = pulsefish_eodphases[fish]['species']
855 # name:
856 if not name:
857 name = species
858 name = name.replace(' ', '_')
859 # write out dictionary:
860 if file is None:
861 file = sys.stdout
862 try:
863 file.write('')
864 closeit = False
865 except AttributeError:
866 file = open(file, 'w')
867 closeit = True
868 n = 6
869 file.write(name + ' = \\\n')
870 file.write(' dict(times=(')
871 file.write(', '.join([f'{a:.5g}' for a in phase_times[:n]]))
872 for k in range(n, len(phase_times), n):
873 file.write(',\n')
874 file.write(' ' * (9 + 12))
875 file.write(', '.join([f'{a:.5g}' for a in phase_times[k:k + n]]))
876 if len(phase_times) == 1:
877 file.write(',')
878 file.write('),\n')
879 file.write(' ' * 9 + 'amplitudes=(')
880 file.write(', '.join([f'{p:.5g}' for p in phase_amplitudes[:n]]))
881 for k in range(n, len(phase_amplitudes), n):
882 file.write(',\n')
883 file.write(' ' * (9 + 8))
884 file.write(', '.join([f'{p:.5g}' for p in phase_amplitudes[k:k + n]]))
885 if len(phase_amplitudes) == 1:
886 file.write(',')
887 file.write('),\n')
888 file.write(' ' * 9 + 'stdevs=(')
889 file.write(', '.join([f'{s:.5g}' for s in phase_stdevs[:n]]))
890 for k in range(n, len(phase_stdevs), n):
891 file.write(',\n')
892 file.write(' ' * (9 + 8))
893 file.write(', '.join([f'{s:.5g}' for s in phase_stdevs[k:k + n]]))
894 if len(phase_stdevs) == 1:
895 file.write(',')
896 if species:
897 file.write('),\n')
898 file.write(' ' * 9 + f'species=\'{species}\')\n')
899 phases['species'] = species
900 else:
901 file.write('))\n')
902 if closeit:
903 file.close()
904 return phases
907def generate_waveform(filename):
908 """Interactively generate audio file with simulated EOD waveforms.
910 Parameters needed to generate EOD waveforms are take from console input.
912 Parameters
913 ----------
914 filename: string
915 Name of file inclusively extension (e.g. '.wav')
916 used to store the simulated EOD waveforms.
917 """
918 import os
919 from audioio import write_audio
920 from thunderlab.consoleinput import read, select, save_inputs
921 # generate file:
922 rate = read('Sampling rate in Hz', '44100', float, 1.0)
923 duration = read('Duration in seconds', '10', float, 0.001)
924 nfish = read('Number of fish', '1', int, 1)
925 ndata = read('Number of electrodes', '1', int, 1)
926 fish_spread = 1
927 if ndata > 1:
928 fish_spread = read('Number of electrodes fish are spread over', '2', int, 1)
929 data = np.random.randn(int(duration*rate), ndata)*0.01
930 fish_indices = np.random.randint(ndata, size=nfish)
931 eodt = 'a'
932 eodf = 800.0
933 eoda = 1.0
934 eodsig = 'n'
935 pulse_jitter = 0.1
936 chirp_freq = 5.0
937 chirp_size = 100.0
938 chirp_width = 0.01
939 chirp_kurtosis = 1.0
940 rise_freq = 0.1
941 rise_size = 10.0
942 rise_tau = 1.0
943 rise_decay_tau = 10.0
944 for k in range(nfish):
945 print('')
946 fish = 'Fish %d: ' % (k+1)
947 eodt = select(fish + 'EOD type', eodt, ['a', 'e', '1', '2', '3'],
948 ['Apteronotus', 'Eigenmannia',
949 'monophasic pulse', 'biphasic pulse', 'triphasic pulse'])
950 eodf = read(fish + 'EOD frequency in Hz', '%g'%eodf, float, 1.0, 3000.0)
951 eoda = read(fish + 'EOD amplitude', '%g'%eoda, float, 0.0, 10.0)
952 if eodt in 'ae':
953 eodsig = select(fish + 'Add communication signals', eodsig, ['n', 'c', 'r'],
954 ['fixed EOD', 'chirps', 'rises'])
955 eodfreq = eodf
956 if eodsig == 'c':
957 chirp_freq = read('Number of chirps per second', '%g'%chirp_freq, float, 0.001)
958 chirp_size = read('Size of chirp in Hz', '%g'%chirp_size, float, 1.0)
959 chirp_width = 0.001*read('Width of chirp in ms', '%g'%(1000.0*chirp_width), float, 1.0)
960 eodfreq, _ = chirps(eodf, rate, duration,
961 chirp_freq, chirp_size, chirp_width, chirp_kurtosis)
962 elif eodsig == 'r':
963 rise_freq = read('Number of rises per second', '%g'%rise_freq, float, 0.00001)
964 rise_size = read('Size of rise in Hz', '%g'%rise_size, float, 0.01)
965 rise_tau = read('Time-constant of rise onset in seconds', '%g'%rise_tau, float, 0.01)
966 rise_decay_tau = read('Time-constant of rise decay in seconds', '%g'%rise_decay_tau, float, 0.01)
967 eodfreq = rises(eodf, rate, duration,
968 rise_freq, rise_size, rise_tau, rise_decay_tau)
969 if eodt == 'a':
970 fishdata = eoda*wavefish_eods('Alepto', eodfreq, rate, duration,
971 phase0=0.0, noise_std=0.0)
972 elif eodt == 'e':
973 fishdata = eoda*wavefish_eods('Eigenmannia', eodfreq, rate,
974 duration, phase0=0.0, noise_std=0.0)
975 else:
976 pulse_jitter = read(fish + 'CV of pulse jitter', '%g'%pulse_jitter, float, 0.0, 2.0)
977 if eodt == '1':
978 fishdata = eoda*pulsefish_eods('Monophasic', eodf, rate, duration,
979 jitter_cv=pulse_jitter, noise_std=0.0)
980 elif eodt == '2':
981 fishdata = eoda*pulsefish_eods('Biphasic', eodf, rate, duration,
982 jitter_cv=pulse_jitter, noise_std=0.0)
983 elif eodt == '3':
984 fishdata = eoda*pulsefish_eods('Triphasic', eodf, rate, duration,
985 jitter_cv=pulse_jitter, noise_std=0.0)
986 i = fish_indices[k]
987 for j in range(fish_spread):
988 data[:, (i+j)%ndata] += fishdata*(0.2**j)
990 maxdata = np.max(np.abs(data))
991 write_audio(filename, 0.9*data/maxdata, rate)
992 input_file = os.path.splitext(filename)[0] + '.inp'
993 save_inputs(input_file)
994 print(f'\nWrote fakefish data to file "{filename}".')
997def generate_testfiles(log_spec=True):
998 """Generate recordings of various pulse EODs and their spectrum.
1000 The spectrum is analytically computed and thus can be used to test
1001 analyis tools.
1003 Three files are generated for each pulse waveform:
1005 1. A wav file with the simulated recording.
1006 2. A csv file with the spectra (see below for details).
1007 3. A pdf file with a plot showing the averaged EOD pulse and the spectrum.
1009 The csv file contains six columns separated by semicolons of the
1010 single pulse spectra:
1012 1. `f`: the frequency components in Hertz
1013 2- `real`: real part of the Fourier spectrum
1014 3. `imag`: imaginary part of the Fourier spectrum
1015 4. `ampl`: amplitude spectrum, i.e. magnitude of Fourier spectrum
1016 5. `energy`: energy spectrum, i.e. squared amplitude
1017 6. `level`: energy sepctrum in decibel relative to maximum energy
1020 Parameters
1021 ----------
1022 log_spec: bool
1023 If True, plot decibel agains logarithmic frequency axis.
1024 If False, plot energy agains linear frequency axis.
1026 Returns
1027 -------
1028 files: list of str
1029 List of generated files.
1030 """
1031 import matplotlib.pyplot as plt
1032 from scipy.signal import find_peaks
1033 from audioio import write_audio
1034 from thunderlab.eventdetection import snippets
1035 from .pulseanalysis import pulse_spectrum
1037 np.seterr(all="ignore")
1038 rng = np.random.default_rng()
1039 rate = 50000.0
1040 duration = 20.0
1041 sigmat = 0.0002
1042 monophasic = dict(name='monophasic',
1043 times=(0,),
1044 amplitudes=(1.0,),
1045 stdevs=(sigmat,))
1046 biphasic30 = dict(name='biphasic30',
1047 times=(0, 2*sigmat),
1048 amplitudes=(1.0, -0.3),
1049 stdevs=(sigmat, sigmat))
1050 biphasic60 = dict(name='biphasic60',
1051 times=(0, 2*sigmat),
1052 amplitudes=(1.0, -0.6),
1053 stdevs=(sigmat, sigmat))
1054 biphasic100 = dict(name='biphasic100',
1055 times=(0, 2*sigmat),
1056 amplitudes=(1.0, -1.0),
1057 stdevs=(sigmat, sigmat))
1058 triphasic = dict(name='triphasic', **Triphasic_phases)
1059 files = []
1060 for pulse in [monophasic, biphasic30, biphasic60, biphasic100, triphasic]:
1061 print(pulse['name'], '...')
1062 # fake recording:
1063 eodf = rng.uniform(40.0, 120.0)
1064 eoda = rng.uniform(0.8, 8.0)
1065 data = eoda*pulsefish_eods(pulse, eodf, rate, duration,
1066 jitter_cv=0.01, noise_std=0.002)
1067 maxdata = np.max(np.abs(data))
1068 fac = 0.9/maxdata
1069 metadata = dict(gain=f'{1/fac:.3f}mV')
1070 write_audio(pulse['name'] + '.wav', fac*data, rate, metadata)
1071 print(f' wrote {pulse['name']}.wav')
1072 files.append(pulse['name'] + '.wav')
1073 # average EOD pulse:
1074 tmax = 0.002
1075 idxs, _ = find_peaks(data, prominence=0.9*eoda)
1076 iw = int(tmax*rate)
1077 snips = snippets(data, idxs, start=-iw, stop=iw)
1078 eod_data = np.mean(snips, 0)
1079 eod_time = (np.arange(len(eod_data)) - iw)/rate
1080 # analytic spectra:
1081 freqs = np.arange(0, rate/2, 1.0)
1082 spec = pulsefish_spectrum(pulse, freqs)
1083 spec *= eoda
1084 ampl = np.abs(spec)
1085 energy = ampl**2
1086 level = 10*np.log10(energy/np.max(energy))
1087 spec_data = np.zeros((len(freqs), 6))
1088 spec_data[:, 0] = freqs
1089 spec_data[:, 1] = np.real(spec)
1090 spec_data[:, 2] = np.imag(spec)
1091 spec_data[:, 3] = ampl
1092 spec_data[:, 4] = energy
1093 spec_data[:, 5] = level
1094 np.savetxt(pulse['name'] + '.csv', spec_data, fmt='%g', delimiter=';',
1095 header='f/Hz;real;imag;ampl;energy;level/dB')
1096 print(f' wrote {pulse['name']}.csv')
1097 files.append(pulse['name'] + '.csv')
1098 # numerical spectrum:
1099 nfreqs, nenergy = pulse_spectrum(eod_data, rate, 1.0, 0.05)
1100 nlevel = 10*np.log10(nenergy/np.max(energy))
1101 # check normalization:
1102 print(f' integral over analytic energy spectrum: {np.sum(energy)*freqs[1]:9.3e}')
1103 print(f' integral over numeric energy spectrum : {np.sum(nenergy)*nfreqs[1]:9.3e}')
1104 print(f' integral over squared signal : {np.sum(eod_data**2)/rate:9.3e}')
1105 # plot waveform:
1106 pa = np.sum(eod_data[eod_data >= 0])/rate
1107 na = np.sum(-eod_data[eod_data <= 0])/rate
1108 balance = (pa - na)/(pa + na)
1109 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3),
1110 layout='constrained')
1111 fig.suptitle(pulse['name'])
1112 ax1.axhline(0, color='gray')
1113 ax1.plot(1000*eod_time, eod_data, color='C0', label='data')
1114 ax1.text(0.1, 0.05, f'polarity balance = {100*balance:.0f}%',
1115 transform=ax1.transAxes)
1116 ax1.set_xlim(-1000*tmax, 1000*tmax)
1117 ax1.set_ylim(-1.1*eoda, 1.1*eoda)
1118 ax1.set_xlabel('time [ms]')
1119 ax1.set_ylabel('averaged EOD')
1120 # plot spectrum:
1121 ip = np.argmax(level)
1122 fmax = freqs[ip]
1123 fmaxpos = fmax if fmax > 1 else 1
1124 if log_spec:
1125 pmax = level[ip]
1126 ax2.plot(nfreqs, nlevel + 0.5, color='C1', label='numeric')
1127 ax2.plot(freqs, level, color='C3', label='analytic')
1128 ax2.plot(fmaxpos, pmax, 'o', color='C3')
1129 ax2.set_xlim(1, 1e4)
1130 ax2.set_xscale('log')
1131 ax2.set_ylim(-60, 10)
1132 ax2.set_ylabel('energy [dB]')
1133 dc = level[0]
1134 ax2.text(2, dc - 4, f'{dc:.0f}dB')
1135 ax2.text(fmaxpos*1.05, pmax + 1, f'{fmax:.0f}Hz')
1136 ax2.legend(loc='lower left', frameon=False)
1137 else:
1138 pmax = energy[ip]
1139 ax2.plot(nfreqs, nenergy, color='C1', label='numeric')
1140 ax2.plot(freqs, energy, color='C3', label='analytic')
1141 ax2.plot(fmaxpos, pmax, 'o', color='C3')
1142 ax2.set_xlim(0, 3000)
1143 ax2.set_ylim(bottom=0)
1144 ax2.set_ylabel('energy')
1145 ax2.text(fmaxpos + 200, pmax, f'{fmax:.0f}Hz', va='center')
1146 ax2.legend(loc='upper right', frameon=False)
1147 ax2.set_xlabel('frequency [Hz]')
1148 fig.savefig(pulse['name'] + '.pdf')
1149 plt.show()
1150 plt.close(fig)
1151 print(f' wrote {pulse['name']}.pdf')
1152 files.append(pulse['name'] + '.pdf')
1153 return files
1156def demo():
1157 import matplotlib.pyplot as plt
1158 rate = 40000.0 # in Hz
1159 duration = 10.0 # in sec
1161 inset_len = 0.01 # in sec
1162 inset_indices = int(inset_len*rate)
1163 ws_fac = 0.1 # whitespace factor or ylim (between 0. and 1.)
1165 # generate data:
1166 eodf = 400.0
1167 wavefish = wavefish_eods('Alepto', eodf, rate, duration, noise_std=0.02)
1168 eodf = 650.0
1169 wavefish += 0.5*wavefish_eods('Eigenmannia', eodf, rate, duration)
1171 pulsefish = pulsefish_eods('Biphasic', 80.0, rate, duration,
1172 noise_std=0.02, jitter_cv=0.1, first_pulse=inset_len/2)
1173 time = np.arange(len(wavefish))/rate
1175 fig, ax = plt.subplots(2, 2, figsize=(19, 10), layout='constrained')
1177 # get proper wavefish ylim
1178 ymin = np.min(wavefish)
1179 ymax = np.max(wavefish)
1180 dy = ws_fac*(ymax - ymin)
1181 ymin -= dy
1182 ymax += dy
1184 # complete wavefish:
1185 ax[0][0].set_title('Wavefish')
1186 ax[0][0].set_ylim(ymin, ymax)
1187 ax[0][0].plot(time, wavefish)
1189 # wavefish zoom in:
1190 ax[0][1].set_title('Wavefish ZOOM IN')
1191 ax[0][1].set_ylim(ymin, ymax)
1192 ax[0][1].plot(time[:inset_indices], wavefish[:inset_indices], '-o')
1194 # get proper pulsefish ylim
1195 ymin = np.min(pulsefish)
1196 ymax = np.max(pulsefish)
1197 dy = ws_fac*(ymax - ymin)
1198 ymin -= dy
1199 ymax += dy
1201 # complete pulsefish:
1202 ax[1][0].set_title('Pulsefish')
1203 ax[1][0].set_ylim(ymin, ymax)
1204 ax[1][0].plot(time, pulsefish)
1206 # pulsefish zoom in:
1207 ax[1][1].set_title('Pulsefish ZOOM IN')
1208 ax[1][1].set_ylim(ymin, ymax)
1209 ax[1][1].plot(time[:inset_indices], pulsefish[:inset_indices], '-o')
1211 for row in ax:
1212 for c_ax in row:
1213 c_ax.set_xlabel('Time [sec]')
1214 c_ax.set_ylabel('Amplitude')
1216 # chirps:
1217 chirps_freq = chirps(600.0, rate, duration)
1218 chirps_data = wavefish_eods('Alepto', chirps_freq, rate)
1220 # rises:
1221 rises_freq = rises(600.0, rate, duration, rise_size=20.0)
1222 rises_data = wavefish_eods('Alepto', rises_freq, rate)
1224 nfft = 256
1225 fig, ax = plt.subplots(2, 1, figsize=(19, 10), layout='constrained')
1226 ax[0].set_title('Chirps')
1227 ax[0].specgram(chirps_data, Fs=rate, NFFT=nfft, noverlap=nfft//16)
1228 time = np.arange(len(chirps_freq))/rate
1229 ax[0].plot(time[:-nfft//2], chirps_freq[nfft//2:], '-k', lw=2)
1230 ax[0].set_ylim(0.0, 3000.0)
1231 ax[0].set_ylabel('Frequency [Hz]')
1233 nfft = 4096
1234 ax[1].set_title('Rises')
1235 ax[1].specgram(rises_data, Fs=rate, NFFT=nfft, noverlap=nfft//2)
1236 time = np.arange(len(rises_freq))/rate
1237 ax[1].plot(time[:-nfft//4], rises_freq[nfft//4:], '-k', lw=2)
1238 ax[1].set_ylim(500.0, 700.0)
1239 ax[1].set_ylabel('Frequency [Hz]')
1240 ax[1].set_xlabel('Time [s]')
1242 plt.show()
1245def main(args=[]):
1246 from .version import __year__
1248 if len(args) > 0:
1249 if (len(args) == 1 and args[0].lower() != '-t') and args[0] != '-s':
1250 print('usage: fakefish [-h|--help] [-s audiofile] [-t] [-T]')
1251 print('')
1252 print('Without arguments, run a demo for illustrating fakefish functionality.')
1253 print('')
1254 print('-s audiofile: writes audiofile with user defined simulated electric fishes.')
1255 print('-t: write audiofiles for a number of pulse waveforms and corresponding analytic spectra in csv files for testing. Plot frequency axis logarithmically.')
1256 print('-T: write audiofiles for a number of pulse waveforms and corresponding analytic spectra in csv files for testing. Plot frequency axis linearly.')
1257 print('')
1258 print(f'by bendalab ({__year__})')
1259 elif args[0] == '-s':
1260 generate_waveform(args[1])
1261 elif args[0] == '-t':
1262 generate_testfiles(True)
1263 elif args[0] == '-T':
1264 generate_testfiles(False)
1265 else:
1266 demo()
1269if __name__ == '__main__':
1270 import sys
1271 main(sys.argv[1:])