Coverage for src/thunderfish/fakefish.py: 96%
540 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-19 20:05 +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_phases()`: position, amplitudes and standard deviations of EOD phases in pulsefish waveforms.
23- `pulsefish_spectrum()`: analytically computed single-pulse spectrum.
24- `pulsefish_waveform()`: compute a single pulsefish EOD.
25- `pulsefish_eods()`: simulate EODs of a pulse-type fish.
26- `normalize_pulsefish()`: normalize times and stdevs of pulse-type EOD waveform.
27- `export_pulsefish()`: serialize pulsefish parameter to python code.
30## Waveform generation
32- `generate_testfiles()`: generate recordings of various pulse EODs and their spectrum.
33- `generate_waveform()`: interactively generate audio file with simulated EOD waveforms.
34"""
36import sys
37import numpy as np
39from thunderlab.fourier import normalize_fourier_coeffs, fourier_synthesis
42species_name = dict(Sine='Sinewave',
43 Alepto='Apteronotus leptorhynchus',
44 Arostratus='Apteronotus rostratus',
45 Eigenmannia='Eigenmannia spec.',
46 Sternarchella='Sternarchella terminalis',
47 Sternopygus='Sternopygus dariensis')
48"""Translate species ids used by wavefish_harmonics and pulsefish_eodphases to full species names.
49"""
52def abbrv_genus(name):
53 """Abbreviate genus in a species name.
55 Parameters
56 ----------
57 name: string
58 Full species name of the form 'Genus species subspecies'
60 Returns
61 -------
62 name: string
63 The species name with abbreviated genus, i.e. 'G. species subspecies'
64 """
65 ns = name.split()
66 return ns[0][0] + '. ' + ' '.join(ns[1:])
69# EODs of various wavefish species:
71Sine_harmonics = dict(amplitudes=(0, 1.0,), phases=(0, -0.5*np.pi,),
72 species='Sinewave')
74Apteronotus_leptorhynchus_harmonics = \
75 dict(amplitudes=(0, 0.85271, 0.27593, 0.078503, 0.01434,
76 0.0067499, 0.0028604, 0.00032096, 0.00040814,
77 0.00015664, 2.2418e-05, 8.5942e-06, 1.5948e-05,
78 1.4839e-05, 3.6825e-05, 2.4283e-05, 2.251e-05,
79 2.3141e-05, 5.6629e-05, 4.7806e-05, 2.1825e-05),
80 phases=(0, -0.1131, 0.75971, -1.4472, -0.68193,
81 -1.1429, 1.9871, -1.0602, 1.6677, 0.50438,
82 1.3706, 0.069786, -2.911, 1.2664, -1.8063,
83 -1.9486, -3.0924, -1.2103, -1.0031, -0.040183,
84 2.8477),
85 species='Apteronotus leptorhynchus')
87Apteronotus_rostratus_harmonics = \
88 dict(amplitudes=(0, 0.64707, 0.43874, 0.063592, 0.07379, 0.040199,
89 0.023073, 0.0097678),
90 phases=(0, 0.72800,
91 -0.78204, -2.88680, 0.67080, 0.47050, -0.46860,
92 -3.62210),
93 species='Apteronotus rostratus')
95Eigenmannia_harmonics = \
96 dict(amplitudes=(0, 1.0087, 0.23201, 0.060524, 0.020175,
97 0.010087, 0.0080699),
98 phases=(0, -0.22940, -0.24800, 1.35340, 1.24490, 1.11630, -4.41230),
99 species='Eigenmannia')
101Sternarchella_terminalis_harmonics = \
102 dict(amplitudes=(0, 0.11457, 0.4401, 0.41055, 0.20132, 0.061364, 0.011389,
103 0.0057985),
104 phases=(0, -4.28140, 0.87640, 0.11210, -0.77995, -1.45180, -2.39435,
105 -3.56640),
106 species='Sternarchella terminalis')
108Sternopygus_dariensis_harmonics = \
109 dict(amplitudes=(0, 0.98843, 0.41228, 0.047848, 0.11048, 0.022801, 0.030706,
110 0.019018),
111 phases=(0, -0.15550, -0.25670, 1.53540, -3.96690, -3.52320, -1.02759,
112 0.11360),
113 species='Sternopygus dariensis')
115wavefish_harmonics = dict(Sine=Sine_harmonics,
116 Alepto=Apteronotus_leptorhynchus_harmonics,
117 Arostratus=Apteronotus_rostratus_harmonics,
118 Eigenmannia=Eigenmannia_harmonics,
119 Sternarchella=Sternarchella_terminalis_harmonics,
120 Sternopygus=Sternopygus_dariensis_harmonics)
121"""Amplitudes and phases of EOD waveforms of various species of wave-type electric fish."""
124def wavefish_spectrum(fish):
125 """Fourier coefficients (amplitude and phases) of a wavefish EOD.
127 Parameters
128 ----------
129 fish: string, dict or tuple of lists/arrays or 2-D array of float or 1-D array of complex
130 Specify relative amplitudes and phases of the fundamental and
131 its harmonics.
132 If string then take amplitudes and phases from the
133 `wavefish_harmonics` dictionary.
134 If dictionary then take amplitudes and phases from the
135 'amlitudes' and 'phases' keys.
136 If tuple then the first element is the list of amplitudes and
137 the second one the list of relative phases in radians.
138 If 2-D array with less than five columns, then first column is
139 amplitudes and second one phases is radians.
140 If 2-D array with five or more columns, as returned from
141 waveanalysis.analyse_wave(), then take relative amplitudes
142 from third column and phases from fifth colum.
143 If 1-D array of complex or string, take the coefficients as is.
144 In all other formats, the first element is the amplitude of
145 the fundamental frequency.
147 Returns
148 -------
149 coeffs: 1D array of complex
150 For each harmonics the complex valued Fourier coefficient.
151 The first one is the offset and the second one is the coefficient
152 of the fundamental.
154 Raises
155 ------
156 KeyError:
157 Unknown fish.
158 IndexError:
159 Amplitudes and phases differ in length.
161 """
162 if isinstance(fish, np.ndarray) and fish.ndim == 1 and \
163 fish.dtype == complex:
164 return fish
165 elif isinstance(fish, np.ndarray) and fish.ndim == 2 and \
166 fish.dtype == float:
167 if fish.shape[1] < 5:
168 amplitudes = fish[:, 0]
169 phases = fish[:, 1]
170 else:
171 amplitudes = fish[:, 3]
172 phases = fish[:, 5]
173 elif isinstance(fish, (tuple, list)):
174 amplitudes = fish[0]
175 phases = fish[1]
176 elif isinstance(fish, dict):
177 amplitudes = fish['amplitudes']
178 phases = fish['phases']
179 else:
180 if not fish in wavefish_harmonics:
181 raise KeyError('unknown wavefish. Choose one of ' +
182 ', '.join(wavefish_harmonics.keys()))
183 amplitudes = wavefish_harmonics[fish]['amplitudes']
184 phases = wavefish_harmonics[fish]['phases']
185 coeffs = np.zeros(len(amplitudes), dtype=complex)
186 for k in range(len(amplitudes)):
187 coeffs[k] = amplitudes[k]*np.exp(1j*phases[k])
188 return coeffs
189 if len(amplitudes) != len(phases):
190 raise IndexError('need exactly as many phases as amplitudes')
191 # remove NaNs:
192 for k in reversed(range(len(amplitudes))):
193 if np.isfinite(amplitudes[k]) and np.isfinite(phases[k]):
194 amplitudes = amplitudes[:k+1]
195 phases = phases[:k+1]
196 break
197 # make coefficients:
198 coeffs = np.zeros(len(amplitudes) + 1, dtype=complex)
199 for k in range(len(amplitudes)):
200 coeffs[k + 1] = amplitudes[k]*np.exp(1j*phases[k])
201 return coeffs
204def wavefish_eods(fish='Eigenmannia', frequency=100.0, rate=44100.0,
205 duration=1.0, phase0=0.0, noise_std=0.0):
206 """Simulate EOD waveform of a wave-type fish.
208 The waveform is constructed by superimposing sinewaves of integral
209 multiples of the fundamental frequency - the fundamental and its
210 harmonics. The fundamental frequency of the EOD (EODf) is given by
211 `frequency`. With `fish` relative amplitudes and phases of the
212 fundamental and its harmonics are specified.
214 The generated waveform is `duration` seconds long and is sampled with
215 `rate` Hertz. Gaussian white noise with a standard deviation of
216 `noise_std` is added to the generated waveform.
218 Parameters
219 ----------
220 fish: string, dict or tuple of lists/arrays or 2-D array float or 1-D array of complex
221 See wavefish_spectrum() for details.
222 frequency: float or array of floats
223 EOD frequency of the fish in Hertz. Either fixed number or array for
224 time-dependent frequencies.
225 rate: float
226 Sampling rate in Hertz.
227 duration: float
228 Duration of the generated data in seconds. Only used if frequency is scalar.
229 phase0: float
230 Phase offset of the EOD waveform in radians.
231 noise_std: float
232 Standard deviation of additive Gaussian white noise.
234 Returns
235 -------
236 data: array of floats
237 Generated data of a wave-type fish.
239 Raises
240 ------
241 KeyError:
242 Unknown fish.
243 IndexError:
244 Amplitudes and phases differ in length.
245 """
246 # get relative amplitude and phases:
247 coeffs = wavefish_spectrum(fish)
248 # compute phase:
249 if np.isscalar(frequency):
250 phase = np.arange(0, duration, 1.0/rate)
251 phase *= frequency
252 else:
253 phase = np.cumsum(frequency)/rate
254 # generate EOD:
255 iomega = 1j*2*np.pi*phase
256 data = np.zeros(len(phase))
257 for h, c in enumerate(coeffs):
258 data += np.real(c*np.exp(iomega*h + 1j*h*phase0))
259 # add noise:
260 if noise_std > 0:
261 data += noise_std*np.random.randn(len(data))
262 return data
265def normalize_wavefish(fish, mode='peak'):
266 """Normalize amplitudes and phases of wave-type EOD waveform.
268 The amplitudes and phases of the Fourier components are adjusted
269 such that the resulting EOD waveform has a peak-to-peak amplitude
270 of two and the peak of the waveform is at time zero (mode is set
271 to 'peak') or that the fundamental has an amplitude of one and a
272 phase of 0 (mode is set to 'zero').
274 Parameters
275 ----------
276 fish: string, dict or tuple of lists/arrays or 2-D array float or 1-D array of complex
277 See wavefish_spectrum() for details.
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 coeffs = wavefish_spectrum(fish)
293 if mode == 'zero':
294 coeffs = normalize_fourier_coeffs(coeffs)/np.abs(coeffs[1])
295 return coeffs
296 else:
297 # generate waveform:
298 eodf = 100.0
299 rate = 100000.0
300 n = int(2/eodf*rate) + 1
301 data = fourier_synthesis(eodf, coeffs, rate, n)
302 # normalize amplitudes:
303 ampl = 0.5*(np.max(data) - np.min(data))
304 coeffs /= ampl
305 # shift phases:
306 deltat = np.argmax(data[:int(rate/eodf)])/rate
307 deltap = 2j*np.pi*deltat*eodf
308 coeffs *= np.exp(deltap*np.arange(len(coeffs)))
309 return coeffs
312def export_wavefish(fish, name='Unknown_harmonics', species='', file=None):
313 """Serialize wavefish parameter to python code.
315 Add output to the wavefish_harmonics dictionary!
317 Parameters
318 ----------
319 fish: string, dict or tuple of lists/arrays or 2-D array float or 1-D array of complex
320 See wavefish_spectrum() for details.
321 name: str
322 Name of the dictionary to be written. If empty take species name.
323 species: str
324 Name of the fish species.
325 file: str or Path or file or None
326 File path or open file object where to write wavefish dictionary.
328 Returns
329 -------
330 fish: dict
331 Dictionary with amplitudes, phases, and species.
332 """
333 if fish is None or (isinstance(fish, dict) and len(fish) == 0):
334 return {}
335 # get relative amplitude and phases:
336 coeffs = wavefish_spectrum(fish)
337 amplitudes = np.abs(coeffs)
338 phases = np.angle(coeffs)
339 harmonics = dict(amplitudes=amplitudes,
340 phases=phases)
341 # species:
342 if not species:
343 if isinstance(fish, dict):
344 species = fish.get('species', '')
345 elif isinstance(fish, str) and fish in wavefish_harmonics:
346 species = wavefish_harmonics[fish]['species']
347 # name:
348 if not name:
349 name = species
350 name = name.replace(' ', '_')
351 # write out dictionary:
352 if file is None:
353 file = sys.stdout
354 try:
355 file.write('')
356 closeit = False
357 except AttributeError:
358 file = open(file, 'w')
359 closeit = True
360 n = 5
361 file.write(name + ' = \\\n')
362 file.write(' dict(amplitudes=(')
363 file.write(', '.join([f'{a:.5g}' for a in amplitudes[:n]]))
364 for k in range(n, len(amplitudes), n):
365 file.write(',\n')
366 file.write(' ' * (9 + 12))
367 file.write(', '.join([f'{a:.5g}' for a in amplitudes[k:k + n]]))
368 file.write('),\n')
369 file.write(' ' * 9 + 'phases=(')
370 file.write(', '.join([f'{p:.5g}' for p in phases[:n]]))
371 for k in range(n, len(phases), n):
372 file.write(',\n')
373 file.write(' ' * (9 + 8))
374 file.write(', '.join([f'{p:.5g}' for p in phases[k:k + n]]))
375 if species:
376 file.write('),\n')
377 file.write(' ' * 9 + f'species=\'{species}\')\n')
378 harmonics['species'] = species
379 else:
380 file.write('))\n')
381 if closeit:
382 file.close()
383 return harmonics
386def chirps(eodf=100.0, rate=44100.0, duration=1.0, chirp_freq=5.0,
387 chirp_size=100.0, chirp_width=0.01, chirp_kurtosis=1.0, chirp_contrast=0.05):
388 """Simulate frequency trace with chirps.
390 A chirp is modeled as a Gaussian frequency modulation.
391 The first chirp is placed at 0.5/chirp_freq.
393 Parameters
394 ----------
395 eodf: float
396 EOD frequency of the fish in Hertz.
397 rate: float
398 Sampling rate in Hertz.
399 duration: float
400 Duration of the generated data in seconds.
401 chirp_freq: float
402 Frequency of occurance of chirps in Hertz.
403 chirp_size: float
404 Size of the chirp (maximum frequency increase above eodf) in Hertz.
405 chirp_width: float
406 Width of the chirp at 10% height in seconds.
407 chirp_kurtosis: float:
408 Shape of the chirp. =1: Gaussian, >1: more rectangular, <1: more peaked.
409 chirp_contrast: float
410 Maximum amplitude reduction of EOD during chirp.
412 Returns
413 -------
414 frequency: array of floats
415 Generated frequency trace that can be passed on to wavefish_eods().
416 amplitude: array of floats
417 Generated amplitude modulation that can be used to multiply the trace generated by
418 wavefish_eods().
419 """
420 # baseline eod frequency and amplitude modulation:
421 n = len(np.arange(0, duration, 1.0/rate))
422 frequency = eodf * np.ones(n)
423 am = np.ones(n)
424 # time points for chirps:
425 chirp_period = 1.0/chirp_freq
426 chirp_times = np.arange(0.5*chirp_period, duration, chirp_period)
427 # chirp frequency waveform:
428 chirp_t = np.arange(-2.0*chirp_width, 2.0*chirp_width, 1./rate)
429 chirp_sig = 0.5*chirp_width / (2.0*np.log(10.0))**(0.5/chirp_kurtosis)
430 gauss = np.exp(-0.5*((chirp_t/chirp_sig)**2.0)**chirp_kurtosis)
431 # add chirps on baseline eodf:
432 for ct in chirp_times:
433 index = int(ct*rate)
434 i0 = index - len(gauss)//2
435 i1 = i0 + len(gauss)
436 gi0 = 0
437 gi1 = len(gauss)
438 if i0 < 0:
439 gi0 -= i0
440 i0 = 0
441 if i1 >= len(frequency):
442 gi1 -= i1 - len(frequency)
443 i1 = len(frequency)
444 frequency[i0:i1] += chirp_size * gauss[gi0:gi1]
445 am[i0:i1] -= chirp_contrast * gauss[gi0:gi1]
446 return frequency, am
449def rises(eodf=100.0, rate=44100.0, duration=1.0, rise_freq=0.1,
450 rise_size=10.0, rise_tau=1.0, decay_tau=10.0):
451 """Simulate frequency trace with rises.
453 A rise is modeled as a double exponential frequency modulation.
455 Parameters
456 ----------
457 eodf: float
458 EOD frequency of the fish in Hertz.
459 rate: float
460 Sampling rate in Hertz.
461 duration: float
462 Duration of the generated data in seconds.
463 rise_freq: float
464 Frequency of occurance of rises in Hertz.
465 rise_size: float
466 Size of the rise (frequency increase above eodf) in Hertz.
467 rise_tau: float
468 Time constant of the frequency increase of the rise in seconds.
469 decay_tau: float
470 Time constant of the frequency decay of the rise in seconds.
472 Returns
473 -------
474 data: array of floats
475 Generated frequency trace that can be passed on to wavefish_eods().
476 """
477 n = len(np.arange(0, duration, 1.0/rate))
478 # baseline eod frequency:
479 frequency = eodf * np.ones(n)
480 # time points for rises:
481 rise_period = 1.0/rise_freq
482 rise_times = np.arange(0.5*rise_period, duration, rise_period)
483 # rise frequency waveform:
484 rise_t = np.arange(0.0, 5.0*decay_tau, 1./rate)
485 rise = rise_size * (1.0-np.exp(-rise_t/rise_tau)) * np.exp(-rise_t/decay_tau)
486 # add rises on baseline eodf:
487 for r in rise_times:
488 index = int(r*rate)
489 if index+len(rise) > len(frequency):
490 rise_index = len(frequency)-index
491 frequency[index:index+rise_index] += rise[:rise_index]
492 break
493 else:
494 frequency[index:index+len(rise)] += rise
495 return frequency
498# EODs of various pulsefish species:
500Monophasic_phases = \
501 dict(times=(0,),
502 amplitudes=(1,),
503 stdevs=(0.0003,),
504 species='Biphasic')
506Biphasic_phases = \
507 dict(times=(9e-05, 0.00049),
508 amplitudes=(1.1922, -0.95374),
509 stdevs=(0.0003, 0.00025),
510 species='Biphasic')
512Triphasic_phases = \
513 dict(times=(3e-05, 0.00018, 0.00043),
514 amplitudes=(1.2382, -0.9906, 0.12382),
515 stdevs=(0.0001, 0.0001, 0.0002),
516 species='Triphasic')
518pulsefish_eodphases = dict(Monophasic=Monophasic_phases,
519 Biphasic=Biphasic_phases,
520 Triphasic=Triphasic_phases)
521"""Positions, amplitudes, and standard deviations of Gaussians that
522 make up the EOD phases of pulse-type electric fish.
523"""
526def pulsefish_phases(fish):
527 """Position, amplitudes and standard deviations of EOD phases in pulsefish waveforms.
529 Parameters
530 ----------
531 fish: string, dict or tuple of floats/lists/arrays
532 Specify positions, amplitudes, and standard deviations
533 of Gaussians phases that are superimposed to simulate EOD waveforms
534 of pulse-type electric fishes.
535 If string then take positions, amplitudes and standard deviations
536 from the `pulsefish_eodphases` dictionary.
537 If dictionary then take pulse properties from the 'times', 'amlitudes'
538 and 'stdevs' keys.
539 If tuple then the first element is the list of phase positions,
540 the second is the list of corresponding amplitudes, and
541 the third one the list of corresponding standard deviations.
543 Returns
544 -------
545 times : array of floats
546 Positions of the phases.
547 amplitudes : array of floats
548 Amplitudes of the phases.
549 stdevs : array of floats
550 Standard deviations of the phases.
552 Raises
553 ------
554 KeyError:
555 Unknown fish.
556 IndexError:
557 Phase positions, amplitudes, or standard deviations differ in length.
558 """
559 if isinstance(fish, (tuple, list)):
560 times = fish[0]
561 amplitudes = fish[1]
562 stdevs = fish[2]
563 elif isinstance(fish, dict):
564 times = fish['times']
565 amplitudes = fish['amplitudes']
566 stdevs = fish['stdevs']
567 else:
568 if not fish in pulsefish_eodphases:
569 raise KeyError('unknown pulse-type fish. Choose one of ' +
570 ', '.join(pulsefish_eodphases.keys()))
571 phases = pulsefish_eodphases[fish]
572 times = phases['times']
573 amplitudes = phases['amplitudes']
574 stdevs = phases['stdevs']
575 if len(stdevs) != len(amplitudes) or len(stdevs) != len(times):
576 raise IndexError('need exactly as many standard deviations as amplitudes and times')
577 return times, amplitudes, stdevs
580def pulsefish_spectrum(fish, freqs):
581 """Analytically computed single-pulse spectrum.
583 The spectrum is the sum of the spectra of all of the phases. Each
584 phase is the spectrum of a Gaussian (which again is a Gaussian
585 centered at f=0) multiplied with a shift term to account for its
586 position relative to the first phase.
588 You get the amplitude spectral density by taking the absolute value.
589 The energy spectral density is the amplitude spectral density squared:
590 ```
591 spec = pulsefish_spectrum('Biphasic', freqs)
592 ampl = np.abs(spec)
593 energy = np.abs(spec)**2
594 ```
596 Parameters
597 ----------
598 fish: string, dict or tuple of floats/lists/arrays
599 Specify positions, amplitudes, and standard deviations
600 of Gaussians phases that are superimposed to simulate EOD waveforms
601 of pulse-type electric fishes.
602 If string then take positions, amplitudes and standard deviations
603 from the `pulsefish_eodphases` dictionary.
604 If dictionary then take pulse properties from the 'times', 'amlitudes'
605 and 'stdevs' keys.
606 If tuple then the first element is the list of phase positions,
607 the second is the list of corresponding amplitudes, and
608 the third one the list of corresponding standard deviations.
609 freqs: 1-D array of float
610 Frequencies at which the spectrum is computed.
612 Returns
613 -------
614 spectrum: 1-D array of complex
615 The one-sided complex-valued spectrum of the single pulse EOD.
616 The squared magnitude (the energy spectrum) has
617 the unit (x s)^2 = x^2 s/Hz.
619 """
620 times, ampls, stdevs = pulsefish_phases(fish)
621 spec = np.zeros(len(freqs), dtype=complex)
622 for dt, a, s in zip(times, ampls, stdevs):
623 gauss = a*np.sqrt(2*np.pi)*s*np.exp(-0.5*(2*np.pi*s*freqs)**2)
624 shift = np.exp(-2j*np.pi*freqs*dt)
625 spec += gauss*shift
626 spec *= np.sqrt(2) # because of one-sided spectrum
627 return spec
630def pulsefish_waveform(fish, t):
631 """Compute a single pulsefish EOD.
633 You may use pulsefish_parameter() to pass the paramters returned
634 by pulsefish_phases() to this function.
636 Parameters
637 ----------
638 fish: string, dict or tuple of floats/lists/arrays
639 Specify positions, amplitudes, and standard deviations
640 of Gaussians phases that are superimposed to simulate EOD waveforms
641 of pulse-type electric fishes.
642 If string then take positions, amplitudes and standard deviations
643 from the `pulsefish_eodphases` dictionary.
644 If dictionary then take pulse properties from the 'times', 'amlitudes'
645 and 'stdevs' keys.
646 If tuple then the first element is the list of phase positions,
647 the second is the list of corresponding amplitudes, and
648 the third one the list of corresponding standard deviations.
649 t: array of float
650 The time array over which the pulse waveform is evaluated.
652 Returns
653 -------
654 pulse: array of float
655 The pulse waveform for the times given in `t`.
657 """
658 times, ampls, stdevs = pulsefish_phases(fish)
659 pulse = np.zeros(len(t))
660 for dt, a, s in zip(times, ampls, stdevs):
661 pulse += a*np.exp(-0.5*((t - dt)/s)**2)
662 return pulse
665def pulsefish_eods(fish='Biphasic', frequency=100.0, rate=44100.0,
666 duration=1.0, noise_std=0.01, jitter_cv=0.1,
667 first_pulse=None):
668 """Simulate EODs of a pulse-type fish.
670 Pulses are spaced by 1/frequency, jittered as determined by
671 jitter_cv. Each pulse is a combination of Gaussian phases, whose
672 positions, amplitudes and widths are given by `fish`.
674 The generated waveform is duration seconds long and is sampled
675 with rate Hertz. Gaussian white noise with a standard deviation
676 of noise_std is added to the generated pulse train.
678 Parameters
679 ----------
680 fish: string, dict or tuple of floats/lists/arrays
681 Specify positions, amplitudes, and standard deviations
682 of Gaussians phases that are superimposed to simulate EOD waveforms
683 of pulse-type electric fishes.
684 If string then take positions, amplitudes and standard deviations
685 from the `pulsefish_eodphases` dictionary.
686 If dictionary then take pulse properties from the 'times', 'amlitudes'
687 and 'stdevs' keys.
688 If tuple then the first element is the list of phase positions,
689 the second is the list of corresponding amplitudes, and
690 the third one the list of corresponding standard deviations.
691 frequency: float
692 EOD frequency of the fish in Hz.
693 rate: float
694 Sampling Rate in Hz.
695 duration: float
696 Duration of the generated data in seconds.
697 noise_std: float
698 Standard deviation of additive Gaussian white noise.
699 jitter_cv: float
700 Gaussian distributed jitter of pulse times as coefficient of variation
701 of inter-pulse intervals.
702 first_pulse: float or None
703 The position of the first pulse. If None it is choosen automatically
704 depending on pulse width, jitter, and frequency.
706 Returns
707 -------
708 data: array of floats
709 Generated data of a pulse-type fish.
711 Raises
712 ------
713 KeyError:
714 Unknown fish.
715 IndexError:
716 Phase positions, amplitudes, or standard deviations differ in length.
718 """
719 # get phase properties:
720 phase_times, phase_amplitudes, phase_stdevs = pulsefish_phases(fish)
721 # time axis for single pulse:
722 min_time_inx = np.argmin(phase_times)
723 max_time_inx = np.argmax(phase_times)
724 tmax = max(np.abs(phase_times[min_time_inx]-4.0*phase_stdevs[min_time_inx]),
725 np.abs(phase_times[max_time_inx]+4.0*phase_stdevs[max_time_inx]))
726 x = np.arange(-tmax, tmax, 1.0/rate)
727 pulse_duration = x[-1] - x[0]
729 # generate a single pulse:
730 pulse = pulsefish_waveform((phase_times, phase_amplitudes, phase_stdevs), x)
731 #pulse = np.zeros(len(x))
732 #for time, ampl, std in zip(phase_times, phase_amplitudes, phase_stdevs):
733 # pulse += ampl * np.exp(-0.5*((x-time)/std)**2)
734 poffs = len(pulse)//2
736 # paste the pulse into the noise floor:
737 time = np.arange(0, duration, 1.0/rate)
738 data = np.random.randn(len(time)) * noise_std
739 period = 1.0/frequency
740 jitter_std = period * jitter_cv
741 if first_pulse is None:
742 first_pulse = np.max([pulse_duration, 3.0*jitter_std])
743 pulse_times = np.arange(first_pulse, duration, period )
744 pulse_times += jitter_std*np.random.randn(len(pulse_times))
745 pulse_indices = np.round(pulse_times * rate).astype(int)
746 for inx in pulse_indices[(pulse_indices>=poffs)&(pulse_indices-poffs+len(pulse)<len(data))]:
747 data[inx-poffs:inx-poffs+len(pulse)] += pulse
748 return data
751def normalize_pulsefish(fish):
752 """Normalize times and stdevs of pulse-type EOD waveform.
754 The positions and amplitudes of Gaussian phases are adjusted such
755 that the resulting EOD waveform has a maximum peak amplitude of one
756 and has the largest phase at time zero.
758 Parameters
759 ----------
760 fish: string, dict or tuple of floats/lists/arrays
761 Specify positions, amplitudes, and standard deviations
762 of Gaussians phases that are superimposed to simulate EOD waveforms
763 of pulse-type electric fishes.
764 If string then take positions, amplitudes and standard deviations
765 from the `pulsefish_eodphases` dictionary.
766 If dictionary then take pulse properties from the 'times', 'amlitudes'
767 and 'stdevs' keys.
768 If tuple then the first element is the list of phase positions,
769 the second is the list of corresponding amplitudes, and
770 the third one the list of corresponding standard deviations.
772 Returns
773 -------
774 fish: dict
775 Dictionary with adjusted times and standard deviations.
776 """
777 # get phase properties:
778 phase_times, phase_amplitudes, phase_stdevs = pulsefish_phases(fish)
779 # generate waveform:
780 eodf = 10.0
781 rate = 100000.0
782 first_pulse = 0.5/eodf
783 data = pulsefish_eods(fish, eodf, rate, 1.0/eodf, noise_std=0.0,
784 jitter_cv=0.0, first_pulse=first_pulse)
785 # maximum peak:
786 idx = np.argmax(np.abs(data))
787 # normalize amplitudes:
788 ampl = data[idx]
789 newamplitudes = np.array(phase_amplitudes)/ampl
790 # shift times:
791 deltat = idx/rate - first_pulse
792 newtimes = np.array(phase_times) - deltat
793 # store and return:
794 phases = dict(times=newtimes,
795 amplitudes=newamplitudes,
796 stdevs=phase_stdevs)
797 return phases
800def export_pulsefish(fish, name='Unknown_phases', species='', file=None):
801 """Serialize pulsefish parameter to python code.
803 Add output to the pulsefish_eodphases dictionary!
805 Parameters
806 ----------
807 fish: string, dict or tuple of floats/lists/arrays
808 Specify positions, amplitudes, and standard deviations
809 of Gaussians phases that are superimposed to simulate EOD waveforms
810 of pulse-type electric fishes.
811 If string then take positions, amplitudes, standard deviations, and
812 species from the `pulsefish_eodphases` dictionary.
813 If dictionary then take pulse properties and species from the 'times',
814 'amlitudes', 'stdevs' and 'species' keys.
815 If tuple then the first element is the list of phase positions,
816 the second is the list of corresponding amplitudes, and
817 the third one the list of corresponding standard deviations.
818 name: string
819 Name of the dictionary to be written. If empty take species name.
820 species: str
821 Name of the fish species.
822 file: str or Path or file or None
823 File name or open file object where to write pulsefish dictionary.
825 Returns
826 -------
827 fish: dict
828 Dictionary with phase times, amplitudes, standard deviations and
829 species.
830 """
831 if fish is None or (isinstance(fish, dict) and len(fish) == 0):
832 return {}
833 # get phase properties:
834 phase_times, phase_amplitudes, phase_stdevs = pulsefish_phases(fish)
835 phases = dict(times=phase_times,
836 amplitudes=phase_amplitudes,
837 stdevs=phase_stdevs)
838 # species:
839 if not species:
840 if isinstance(fish, dict):
841 species = fish.get('species', '')
842 elif isinstance(fish, str) and fish in pulsefish_eodphases:
843 species = pulsefish_eodphases[fish]['species']
844 # name:
845 if not name:
846 name = species
847 name = name.replace(' ', '_')
848 # write out dictionary:
849 if file is None:
850 file = sys.stdout
851 try:
852 file.write('')
853 closeit = False
854 except AttributeError:
855 file = open(file, 'w')
856 closeit = True
857 n = 6
858 file.write(name + ' = \\\n')
859 file.write(' dict(times=(')
860 file.write(', '.join([f'{a:.5g}' for a in phase_times[:n]]))
861 for k in range(n, len(phase_times), n):
862 file.write(',\n')
863 file.write(' ' * (9 + 12))
864 file.write(', '.join([f'{a:.5g}' for a in phase_times[k:k + n]]))
865 if len(phase_times) == 1:
866 file.write(',')
867 file.write('),\n')
868 file.write(' ' * 9 + 'amplitudes=(')
869 file.write(', '.join([f'{p:.5g}' for p in phase_amplitudes[:n]]))
870 for k in range(n, len(phase_amplitudes), n):
871 file.write(',\n')
872 file.write(' ' * (9 + 8))
873 file.write(', '.join([f'{p:.5g}' for p in phase_amplitudes[k:k + n]]))
874 if len(phase_amplitudes) == 1:
875 file.write(',')
876 file.write('),\n')
877 file.write(' ' * 9 + 'stdevs=(')
878 file.write(', '.join([f'{s:.5g}' for s in phase_stdevs[:n]]))
879 for k in range(n, len(phase_stdevs), n):
880 file.write(',\n')
881 file.write(' ' * (9 + 8))
882 file.write(', '.join([f'{s:.5g}' for s in phase_stdevs[k:k + n]]))
883 if len(phase_stdevs) == 1:
884 file.write(',')
885 if species:
886 file.write('),\n')
887 file.write(' ' * 9 + f'species=\'{species}\')\n')
888 phases['species'] = species
889 else:
890 file.write('))\n')
891 if closeit:
892 file.close()
893 return phases
896def generate_waveform(filename):
897 """Interactively generate audio file with simulated EOD waveforms.
899 Parameters needed to generate EOD waveforms are take from console input.
901 Parameters
902 ----------
903 filename: string
904 Name of file inclusively extension (e.g. '.wav')
905 used to store the simulated EOD waveforms.
906 """
907 import os
908 from audioio import write_audio
909 from thunderlab.consoleinput import read, select, save_inputs
910 # generate file:
911 rate = read('Sampling rate in Hz', '44100', float, 1.0)
912 duration = read('Duration in seconds', '10', float, 0.001)
913 nfish = read('Number of fish', '1', int, 1)
914 ndata = read('Number of electrodes', '1', int, 1)
915 fish_spread = 1
916 if ndata > 1:
917 fish_spread = read('Number of electrodes fish are spread over', '2', int, 1)
918 data = np.random.randn(int(duration*rate), ndata)*0.01
919 fish_indices = np.random.randint(ndata, size=nfish)
920 eodt = 'a'
921 eodf = 800.0
922 eoda = 1.0
923 eodsig = 'n'
924 pulse_jitter = 0.1
925 chirp_freq = 5.0
926 chirp_size = 100.0
927 chirp_width = 0.01
928 chirp_kurtosis = 1.0
929 rise_freq = 0.1
930 rise_size = 10.0
931 rise_tau = 1.0
932 rise_decay_tau = 10.0
933 for k in range(nfish):
934 print('')
935 fish = 'Fish %d: ' % (k+1)
936 eodt = select(fish + 'EOD type', eodt, ['a', 'e', '1', '2', '3'],
937 ['Apteronotus', 'Eigenmannia',
938 'monophasic pulse', 'biphasic pulse', 'triphasic pulse'])
939 eodf = read(fish + 'EOD frequency in Hz', '%g'%eodf, float, 1.0, 3000.0)
940 eoda = read(fish + 'EOD amplitude', '%g'%eoda, float, 0.0, 10.0)
941 if eodt in 'ae':
942 eodsig = select(fish + 'Add communication signals', eodsig, ['n', 'c', 'r'],
943 ['fixed EOD', 'chirps', 'rises'])
944 eodfreq = eodf
945 if eodsig == 'c':
946 chirp_freq = read('Number of chirps per second', '%g'%chirp_freq, float, 0.001)
947 chirp_size = read('Size of chirp in Hz', '%g'%chirp_size, float, 1.0)
948 chirp_width = 0.001*read('Width of chirp in ms', '%g'%(1000.0*chirp_width), float, 1.0)
949 eodfreq, _ = chirps(eodf, rate, duration,
950 chirp_freq, chirp_size, chirp_width, chirp_kurtosis)
951 elif eodsig == 'r':
952 rise_freq = read('Number of rises per second', '%g'%rise_freq, float, 0.00001)
953 rise_size = read('Size of rise in Hz', '%g'%rise_size, float, 0.01)
954 rise_tau = read('Time-constant of rise onset in seconds', '%g'%rise_tau, float, 0.01)
955 rise_decay_tau = read('Time-constant of rise decay in seconds', '%g'%rise_decay_tau, float, 0.01)
956 eodfreq = rises(eodf, rate, duration,
957 rise_freq, rise_size, rise_tau, rise_decay_tau)
958 if eodt == 'a':
959 fishdata = eoda*wavefish_eods('Alepto', eodfreq, rate, duration,
960 phase0=0.0, noise_std=0.0)
961 elif eodt == 'e':
962 fishdata = eoda*wavefish_eods('Eigenmannia', eodfreq, rate,
963 duration, phase0=0.0, noise_std=0.0)
964 else:
965 pulse_jitter = read(fish + 'CV of pulse jitter', '%g'%pulse_jitter, float, 0.0, 2.0)
966 if eodt == '1':
967 fishdata = eoda*pulsefish_eods('Monophasic', eodf, rate, duration,
968 jitter_cv=pulse_jitter, noise_std=0.0)
969 elif eodt == '2':
970 fishdata = eoda*pulsefish_eods('Biphasic', eodf, rate, duration,
971 jitter_cv=pulse_jitter, noise_std=0.0)
972 elif eodt == '3':
973 fishdata = eoda*pulsefish_eods('Triphasic', eodf, rate, duration,
974 jitter_cv=pulse_jitter, noise_std=0.0)
975 i = fish_indices[k]
976 for j in range(fish_spread):
977 data[:, (i+j)%ndata] += fishdata*(0.2**j)
979 maxdata = np.max(np.abs(data))
980 write_audio(filename, 0.9*data/maxdata, rate)
981 input_file = os.path.splitext(filename)[0] + '.inp'
982 save_inputs(input_file)
983 print(f'\nWrote fakefish data to file "{filename}".')
986def generate_testfiles(log_spec=True):
987 """Generate recordings of various pulse EODs and their spectrum.
989 The spectrum is analytically computed and thus can be used to test
990 analyis tools.
992 Three files are generated for each pulse waveform:
994 1. A wav file with the simulated recording.
995 2. A csv file with the spectra (see below for details).
996 3. A pdf file with a plot showing the averaged EOD pulse and the spectrum.
998 The csv file contains six columns separated by semicolons of the
999 single pulse spectra:
1001 1. `f`: the frequency components in Hertz
1002 2- `real`: real part of the Fourier spectrum
1003 3. `imag`: imaginary part of the Fourier spectrum
1004 4. `ampl`: amplitude spectrum, i.e. magnitude of Fourier spectrum
1005 5. `energy`: energy spectrum, i.e. squared amplitude
1006 6. `level`: energy sepctrum in decibel relative to maximum energy
1009 Parameters
1010 ----------
1011 log_spec: bool
1012 If True, plot decibel agains logarithmic frequency axis.
1013 If False, plot energy agains linear frequency axis.
1015 Returns
1016 -------
1017 files: list of str
1018 List of generated files.
1019 """
1020 import matplotlib.pyplot as plt
1021 from scipy.signal import find_peaks
1022 from audioio import write_audio
1023 from thunderlab.eventdetection import snippets
1024 from .pulseanalysis import pulse_spectrum
1026 np.seterr(all="ignore")
1027 rng = np.random.default_rng()
1028 rate = 50000.0
1029 duration = 20.0
1030 sigmat = 0.0002
1031 monophasic = dict(name='monophasic',
1032 times=(0,),
1033 amplitudes=(1.0,),
1034 stdevs=(sigmat,))
1035 biphasic30 = dict(name='biphasic30',
1036 times=(0, 2*sigmat),
1037 amplitudes=(1.0, -0.3),
1038 stdevs=(sigmat, sigmat))
1039 biphasic60 = dict(name='biphasic60',
1040 times=(0, 2*sigmat),
1041 amplitudes=(1.0, -0.6),
1042 stdevs=(sigmat, sigmat))
1043 biphasic100 = dict(name='biphasic100',
1044 times=(0, 2*sigmat),
1045 amplitudes=(1.0, -1.0),
1046 stdevs=(sigmat, sigmat))
1047 triphasic = dict(name='triphasic', **Triphasic_phases)
1048 files = []
1049 for pulse in [monophasic, biphasic30, biphasic60, biphasic100, triphasic]:
1050 print(pulse['name'], '...')
1051 # fake recording:
1052 eodf = rng.uniform(40.0, 120.0)
1053 eoda = rng.uniform(0.8, 8.0)
1054 data = eoda*pulsefish_eods(pulse, eodf, rate, duration,
1055 jitter_cv=0.01, noise_std=0.002)
1056 maxdata = np.max(np.abs(data))
1057 fac = 0.9/maxdata
1058 metadata = dict(gain=f'{1/fac:.3f}mV')
1059 write_audio(pulse['name'] + '.wav', fac*data, rate, metadata)
1060 print(f' wrote {pulse['name']}.wav')
1061 files.append(pulse['name'] + '.wav')
1062 # average EOD pulse:
1063 tmax = 0.002
1064 idxs, _ = find_peaks(data, prominence=0.9*eoda)
1065 iw = int(tmax*rate)
1066 snips = snippets(data, idxs, start=-iw, stop=iw)
1067 eod_data = np.mean(snips, 0)
1068 eod_time = (np.arange(len(eod_data)) - iw)/rate
1069 # analytic spectra:
1070 freqs = np.arange(0, rate/2, 1.0)
1071 spec = pulsefish_spectrum(pulse, freqs)
1072 spec *= eoda
1073 ampl = np.abs(spec)
1074 energy = ampl**2
1075 level = 10*np.log10(energy/np.max(energy))
1076 spec_data = np.zeros((len(freqs), 6))
1077 spec_data[:, 0] = freqs
1078 spec_data[:, 1] = np.real(spec)
1079 spec_data[:, 2] = np.imag(spec)
1080 spec_data[:, 3] = ampl
1081 spec_data[:, 4] = energy
1082 spec_data[:, 5] = level
1083 np.savetxt(pulse['name'] + '.csv', spec_data, fmt='%g', delimiter=';',
1084 header='f/Hz;real;imag;ampl;energy;level/dB')
1085 print(f' wrote {pulse['name']}.csv')
1086 files.append(pulse['name'] + '.csv')
1087 # numerical spectrum:
1088 nfreqs, nenergy = pulse_spectrum(eod_data, rate, 1.0, 0.05)
1089 nlevel = 10*np.log10(nenergy/np.max(energy))
1090 # check normalization:
1091 print(f' integral over analytic energy spectrum: {np.sum(energy)*freqs[1]:9.3e}')
1092 print(f' integral over numeric energy spectrum : {np.sum(nenergy)*nfreqs[1]:9.3e}')
1093 print(f' integral over squared signal : {np.sum(eod_data**2)/rate:9.3e}')
1094 # plot waveform:
1095 pa = np.sum(eod_data[eod_data >= 0])/rate
1096 na = np.sum(-eod_data[eod_data <= 0])/rate
1097 balance = (pa - na)/(pa + na)
1098 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3),
1099 layout='constrained')
1100 fig.suptitle(pulse['name'])
1101 ax1.axhline(0, color='gray')
1102 ax1.plot(1000*eod_time, eod_data, color='C0', label='data')
1103 ax1.text(0.1, 0.05, f'polarity balance = {100*balance:.0f}%',
1104 transform=ax1.transAxes)
1105 ax1.set_xlim(-1000*tmax, 1000*tmax)
1106 ax1.set_ylim(-1.1*eoda, 1.1*eoda)
1107 ax1.set_xlabel('time [ms]')
1108 ax1.set_ylabel('averaged EOD')
1109 # plot spectrum:
1110 ip = np.argmax(level)
1111 fmax = freqs[ip]
1112 fmaxpos = fmax if fmax > 1 else 1
1113 if log_spec:
1114 pmax = level[ip]
1115 ax2.plot(nfreqs, nlevel + 0.5, color='C1', label='numeric')
1116 ax2.plot(freqs, level, color='C3', label='analytic')
1117 ax2.plot(fmaxpos, pmax, 'o', color='C3')
1118 ax2.set_xlim(1, 1e4)
1119 ax2.set_xscale('log')
1120 ax2.set_ylim(-60, 10)
1121 ax2.set_ylabel('energy [dB]')
1122 dc = level[0]
1123 ax2.text(2, dc - 4, f'{dc:.0f}dB')
1124 ax2.text(fmaxpos*1.05, pmax + 1, f'{fmax:.0f}Hz')
1125 ax2.legend(loc='lower left', frameon=False)
1126 else:
1127 pmax = energy[ip]
1128 ax2.plot(nfreqs, nenergy, color='C1', label='numeric')
1129 ax2.plot(freqs, energy, color='C3', label='analytic')
1130 ax2.plot(fmaxpos, pmax, 'o', color='C3')
1131 ax2.set_xlim(0, 3000)
1132 ax2.set_ylim(bottom=0)
1133 ax2.set_ylabel('energy')
1134 ax2.text(fmaxpos + 200, pmax, f'{fmax:.0f}Hz', va='center')
1135 ax2.legend(loc='upper right', frameon=False)
1136 ax2.set_xlabel('frequency [Hz]')
1137 fig.savefig(pulse['name'] + '.pdf')
1138 plt.show()
1139 plt.close(fig)
1140 print(f' wrote {pulse['name']}.pdf')
1141 files.append(pulse['name'] + '.pdf')
1142 return files
1145def demo():
1146 import matplotlib.pyplot as plt
1147 rate = 40000.0 # in Hz
1148 duration = 10.0 # in sec
1150 inset_len = 0.01 # in sec
1151 inset_indices = int(inset_len*rate)
1152 ws_fac = 0.1 # whitespace factor or ylim (between 0. and 1.)
1154 # generate data:
1155 eodf = 400.0
1156 wavefish = wavefish_eods('Alepto', eodf, rate, duration, noise_std=0.02)
1157 eodf = 650.0
1158 wavefish += 0.5*wavefish_eods('Eigenmannia', eodf, rate, duration)
1160 pulsefish = pulsefish_eods('Biphasic', 80.0, rate, duration,
1161 noise_std=0.02, jitter_cv=0.1, first_pulse=inset_len/2)
1162 time = np.arange(len(wavefish))/rate
1164 fig, ax = plt.subplots(2, 2, figsize=(19, 10), layout='constrained')
1166 # get proper wavefish ylim
1167 ymin = np.min(wavefish)
1168 ymax = np.max(wavefish)
1169 dy = ws_fac*(ymax - ymin)
1170 ymin -= dy
1171 ymax += dy
1173 # complete wavefish:
1174 ax[0][0].set_title('Wavefish')
1175 ax[0][0].set_ylim(ymin, ymax)
1176 ax[0][0].plot(time, wavefish)
1178 # wavefish zoom in:
1179 ax[0][1].set_title('Wavefish ZOOM IN')
1180 ax[0][1].set_ylim(ymin, ymax)
1181 ax[0][1].plot(time[:inset_indices], wavefish[:inset_indices], '-o')
1183 # get proper pulsefish ylim
1184 ymin = np.min(pulsefish)
1185 ymax = np.max(pulsefish)
1186 dy = ws_fac*(ymax - ymin)
1187 ymin -= dy
1188 ymax += dy
1190 # complete pulsefish:
1191 ax[1][0].set_title('Pulsefish')
1192 ax[1][0].set_ylim(ymin, ymax)
1193 ax[1][0].plot(time, pulsefish)
1195 # pulsefish zoom in:
1196 ax[1][1].set_title('Pulsefish ZOOM IN')
1197 ax[1][1].set_ylim(ymin, ymax)
1198 ax[1][1].plot(time[:inset_indices], pulsefish[:inset_indices], '-o')
1200 for row in ax:
1201 for c_ax in row:
1202 c_ax.set_xlabel('Time [sec]')
1203 c_ax.set_ylabel('Amplitude')
1205 # chirps:
1206 chirps_freq = chirps(600.0, rate, duration)
1207 chirps_data = wavefish_eods('Alepto', chirps_freq, rate)
1209 # rises:
1210 rises_freq = rises(600.0, rate, duration, rise_size=20.0)
1211 rises_data = wavefish_eods('Alepto', rises_freq, rate)
1213 nfft = 256
1214 fig, ax = plt.subplots(2, 1, figsize=(19, 10), layout='constrained')
1215 ax[0].set_title('Chirps')
1216 ax[0].specgram(chirps_data, Fs=rate, NFFT=nfft, noverlap=nfft//16)
1217 time = np.arange(len(chirps_freq))/rate
1218 ax[0].plot(time[:-nfft//2], chirps_freq[nfft//2:], '-k', lw=2)
1219 ax[0].set_ylim(0.0, 3000.0)
1220 ax[0].set_ylabel('Frequency [Hz]')
1222 nfft = 4096
1223 ax[1].set_title('Rises')
1224 ax[1].specgram(rises_data, Fs=rate, NFFT=nfft, noverlap=nfft//2)
1225 time = np.arange(len(rises_freq))/rate
1226 ax[1].plot(time[:-nfft//4], rises_freq[nfft//4:], '-k', lw=2)
1227 ax[1].set_ylim(500.0, 700.0)
1228 ax[1].set_ylabel('Frequency [Hz]')
1229 ax[1].set_xlabel('Time [s]')
1231 plt.show()
1234def main(args=[]):
1235 from .version import __year__
1237 if len(args) > 0:
1238 if (len(args) == 1 and args[0].lower() != '-t') and args[0] != '-s':
1239 print('usage: fakefish [-h|--help] [-s audiofile] [-t] [-T]')
1240 print('')
1241 print('Without arguments, run a demo for illustrating fakefish functionality.')
1242 print('')
1243 print('-s audiofile: writes audiofile with user defined simulated electric fishes.')
1244 print('-t: write audiofiles for a number of pulse waveforms and corresponding analytic spectra in csv files for testing. Plot frequency axis logarithmically.')
1245 print('-T: write audiofiles for a number of pulse waveforms and corresponding analytic spectra in csv files for testing. Plot frequency axis linearly.')
1246 print('')
1247 print(f'by bendalab ({__year__})')
1248 elif args[0] == '-s':
1249 generate_waveform(args[1])
1250 elif args[0] == '-t':
1251 generate_testfiles(True)
1252 elif args[0] == '-T':
1253 generate_testfiles(False)
1254 else:
1255 demo()
1258if __name__ == '__main__':
1259 import sys
1260 main(sys.argv[1:])