Coverage for src / thunderfish / chirp.py: 0%
87 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-02-09 14:25 +0000
1"""
2Detection of chirps in weakly electric fish recordings.
4- `chirp_analysis()`: calculates spectrogram, detects fishes and extracts chirp times(combined for all fishes).
5 !!! recommended for short recordings (up to 5 min) where only the chirp times shall be extracted !!!
6- `chirp_detection()`: extracts chirp times with help of given spectrogram and fishlist.
7"""
9import numpy as np
10import matplotlib.pyplot as plt
12from thunderlab.powerspectrum import spectrogram
13from thunderlab.eventdetection import std_threshold, detect_peaks, trim_to_peak
15from .harmonics import harmonic_groups
18def true_chirp_power_drop(chirp_time_idx, power, power_window=100):
19 """
20 Chirp is only accepted as such if the power of the frequency drops down as expected.
22 :param chirp_time_idx: (array) indices of chirps.
23 :param power: (array) power array containing for each timestamp the max value in power of a certain frequency range.
24 :param power_window (int) datapoints arroung a detected chirp used to verify that there is a chirp.
25 :return: chirp_time_idx: (array) indices of chirps that have been confirmed to be chirps.
26 """
28 true_chirp_time_idx = []
30 for i in range(len(chirp_time_idx)):
31 idx0 = int(chirp_time_idx[i] - power_window/2)
32 if idx0 < 0:
33 idx0 = 0
34 idx1 = int(chirp_time_idx[i] + power_window/2)
35 if idx1 > len(power):
36 idx1 = len(power)
38 tmp_median = np.median(power[idx0:idx1])
39 tmp_std = np.std(power[idx0:idx1], ddof=1)
41 if np.min(power[idx0:idx1]) < tmp_median - 3*tmp_std:
42 true_chirp_time_idx.append(chirp_time_idx[i])
44 return np.array(true_chirp_time_idx)
47def true_chirp_power_rise_above(chirp_time_idx, power_above):
49 median_power_above = np.median(power_above)
50 std_power_above = np.std(power_above, ddof=1)
52 if median_power_above > 0.001:
53 print('another fish disturbs the chirp approval! Have to rely on other algorithms.')
54 return chirp_time_idx
55 else:
56 true_chirp_time_idx = []
58 for i in range(len(chirp_time_idx)):
59 if power_above[int(chirp_time_idx[i])] > median_power_above + 3*std_power_above:
60 true_chirp_time_idx.append(chirp_time_idx[i])
62 return true_chirp_time_idx
65def chirp_detection(spectrum, freqs, time, fishlist=None, fundamentals=None, min_power= 0.005, freq_tolerance=1., chirp_th=1.,
66 plot_data_func=None):
67 """
68 Detects chirps on the basis of a spectrogram.
70 :param spectrum: (2d-array) spectrum calulated with the `spectrogram()` function.
71 :param freqs: (array) frequencies of the spectrum.
72 :param time: (array) time of the nffts used in the spectrum.
73 :param fishlist: (array) power und frequncy for each fundamental/harmonic of a detected fish.
74 fishlist[fish][harmonic][frequency, power]
75 :param min_power: (float) minimum power of the fundamental frequency for each fish to participate in chirp detection.
76 :param freq_tolerance: (float) frequency tolerance in the spectrum to detect the power of a certain frequency.
77 :param chirp_th: (float) minimum chirp duration to be accepted as a chirp.
78 :param plot_data: (bool) If True: plots the process of chirp detection.
79 :return:chirp_time: (array) array of times (in sec) where chirps have been detected.
80 """
81 if not hasattr(fundamentals, '__len__'):
82 if not hasattr(fishlist, '__len__'):
83 print('fishlist or fundamentals missing as argument !!!')
84 quit()
85 else:
86 fundamentals = []
87 for fish in fishlist:
88 if fish[0][1] > min_power:
89 fundamentals.append(fish[0][0])
91 chirp_time = np.array([])
92 chirp_freq = np.array([])
94 for enu, fundamental in enumerate(fundamentals):
95 # extract power of only the part of the spectrum that has to be analysied for each fundamental and get the peak
96 # power of every piont in time.
97 power = np.max(spectrum[(freqs >= fundamental - freq_tolerance) & (freqs <= fundamental + freq_tolerance)], axis=0)
98 power_above = np.max(spectrum[(freqs >= fundamental + 50.0 -freq_tolerance) & (freqs <= fundamental + 50.0 + freq_tolerance)], axis=0)
99 #power = np.mean(spectrum[(freqs >= fundamental - freq_tolerance) & (freqs <= fundamental + freq_tolerance)], axis=0)
100 # calculate the slope by calculating the difference in the power
101 power_diff = np.diff(power)
103 # peak detection in the power_diff to detect drops in power indicating chrips
104 threshold = std_threshold(power_diff)
105 peaks, troughs = detect_peaks(power_diff, threshold)
106 troughs, peaks = trim_to_peak(troughs, peaks) # reversed troughs and peaks in output and input to get trim_to_troughs
108 # exclude peaks and troughs with to much time diff to be a chirp
109 # ToDO: not nice !!!
110 peaks = peaks[(troughs - peaks) < chirp_th]
111 troughs = troughs[(troughs - peaks) < chirp_th]
113 if len(troughs) > 0:
114 # chirps times defined as the mean time between the troughs and peaks
115 chirp_time_idx = np.mean([troughs, peaks], axis=0)
117 # exclude detected chirps if the powervalue doesn't drop far enought
118 chirp_time_idx = true_chirp_power_drop(chirp_time_idx, power)
120 # chirp_time_idx = true_chirp_power_rise_above(chirp_time_idx, power_above)
121 # add times of detected chirps to the list.
122 chirp_time = np.concatenate((chirp_time, np.array([time[int(i)] for i in chirp_time_idx])))
123 chirp_freq = np.concatenate((chirp_freq, np.array(fundamental* np.ones(len(chirp_time_idx)))))
125 else:
126 chirp_time = np.array([])
127 chirp_freq = np.array([])
129 if plot_data_func:
130 plot_data_func(enu, chirp_time, time, power, power_above, power_diff, fundamental)
132 return chirp_time, chirp_freq
135def chirp_detection_plot(enu, chirp_time, time, power, power2, power_diff, fundamental):
136 """
137 plots the process of chirp detection.
139 :param enu: (int) indication which fish in list is processed.
140 :param chirp_time: (array) timestamps when chirps have been detected.
141 :param time: (array) time array.
142 :param power: (array) power of a certain frequency band.
143 :param power_diff: (array) slope of the power array.
144 :param fundamental: (float) fundamental frequency around which the algorithm looked for chirps.
145 """
146 try:
147 ax
148 except NameError:
149 fig, ax = plt.subplots()
150 colors = ['r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue']
151 # if enu == 0:
152 # fig, ax = plt.subplots()
153 # colors = ['r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue']
154 ax.plot(chirp_time, np.zeros(len(chirp_time)), 'o', markersize=10, color=colors[enu], alpha=0.8, label='chirps')
155 ax.set_xlabel('time in sec')
156 ax.set_ylabel('power')
158 ax.plot(time, power, colors[enu], marker='.', label='%.1f Hz' % fundamental)
159 ax.plot(time, power2, colors[enu+1], label='%.1f Hz' % (fundamental+50.0))
160 ax.plot(time[:len(power_diff)], power_diff, colors[enu], label='slope')
161 ax.legend(loc='upper right', bbox_to_anchor=(1, 1), frameon=False)
164def chirp_analysis(data, rate):
165 """
166 Performs all steps to detect chirps in a given dataset. This includes spectrogram calculation, fish detection and
167 analysing of specific frequency bands.
168 For further documentation see functions chirp_spectrogram() and chirp_detection().
169 !!! recommended for short recordings (up to 5 min) where only the chirp times shall be extracted !!!
171 :param data: (array) data.
172 :param rate: (float) smaplerate of the data.
173 :param min_power: (float) minimal power of the fish fundamental to include this fish in chirp detection.
174 """
175 freqs, time, spectrum = spectrogram(data, rate, freq_resolution=2., overlap_frac=0.95)
177 power = np.mean(spectrum, axis=1) # spectrum[:, t0:t1] to only let spectrum of certain time....
179 fishlist = harmonic_groups(freqs, power)[0]
181 chirp_time, chirp_freq = chirp_detection(spectrum, freqs, time, fishlist, plot_data_func=chirp_detection_plot)
183 plt.show()
185 return chirp_time, chirp_freq
188if __name__ == '__main__':
189 ###
190 # If you want to test the code I propose to use the file '60427L05.WAV' of the transect
191 # '2016_04_27__downstream_stonewall_at_pool' made in colombia, 2016.
192 ###
193 import sys
194 import matplotlib.pyplot as plt
195 from thunderlab.dataloader import load_data
197 data_file = sys.argv[1]
198 raw_data, rate, unit, amax = load_data(data_file)
200 chirp_time, chirp_freq = chirp_analysis(raw_data[:,0], rate)
202 # power = np.mean(spectrum[:, t:t + nffts_per_psd], axis=1)