Coverage for src/thunderfish/chirp.py: 0%

1""" 

2Detection of chirps in weakly electric fish recordings. 

3 

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""" 

8 

9import numpy as np 

10import matplotlib.pyplot as plt 

11from thunderlab.powerspectrum import spectrogram 

12from thunderlab.eventdetection import std_threshold, detect_peaks, trim_to_peak 

13from .harmonics import harmonic_groups 

14 

15 

16def true_chirp_power_drop(chirp_time_idx, power, power_window=100): 

17 """ 

18 Chirp is only accepted as such if the power of the frequency drops down as expected. 

19 

20 :param chirp_time_idx: (array) indices of chirps. 

21 :param power: (array) power array containing for each timestamp the max value in power of a certain frequency range. 

22 :param power_window (int) datapoints arroung a detected chirp used to verify that there is a chirp. 

23 :return: chirp_time_idx: (array) indices of chirps that have been confirmed to be chirps. 

24 """ 

25 

26 true_chirp_time_idx = [] 

27 

28 for i in range(len(chirp_time_idx)): 

29 idx0 = int(chirp_time_idx[i] - power_window/2) 

30 if idx0 < 0: 

31 idx0 = 0 

32 idx1 = int(chirp_time_idx[i] + power_window/2) 

33 if idx1 > len(power): 

34 idx1 = len(power) 

35 

36 tmp_median = np.median(power[idx0:idx1]) 

37 tmp_std = np.std(power[idx0:idx1], ddof=1) 

38 

39 if np.min(power[idx0:idx1]) < tmp_median - 3*tmp_std: 

40 true_chirp_time_idx.append(chirp_time_idx[i]) 

41 

42 return np.array(true_chirp_time_idx) 

43 

44 

45def true_chirp_power_rise_above(chirp_time_idx, power_above): 

46 

47 median_power_above = np.median(power_above) 

48 std_power_above = np.std(power_above, ddof=1) 

49 

50 if median_power_above > 0.001: 

51 print('another fish disturbs the chirp approval! Have to rely on other algorithms.') 

52 return chirp_time_idx 

53 else: 

54 true_chirp_time_idx = [] 

55 

56 for i in range(len(chirp_time_idx)): 

57 if power_above[int(chirp_time_idx[i])] > median_power_above + 3*std_power_above: 

58 true_chirp_time_idx.append(chirp_time_idx[i]) 

59 

60 return true_chirp_time_idx 

61 

62 

63def chirp_detection(spectrum, freqs, time, fishlist=None, fundamentals=None, min_power= 0.005, freq_tolerance=1., chirp_th=1., 

64 plot_data_func=None): 

65 """ 

66 Detects chirps on the basis of a spectrogram. 

67 

68 :param spectrum: (2d-array) spectrum calulated with the `spectrogram()` function. 

69 :param freqs: (array) frequencies of the spectrum. 

70 :param time: (array) time of the nffts used in the spectrum. 

71 :param fishlist: (array) power und frequncy for each fundamental/harmonic of a detected fish. 

72 fishlist[fish][harmonic][frequency, power] 

73 :param min_power: (float) minimum power of the fundamental frequency for each fish to participate in chirp detection. 

74 :param freq_tolerance: (float) frequency tolerance in the spectrum to detect the power of a certain frequency. 

75 :param chirp_th: (float) minimum chirp duration to be accepted as a chirp. 

76 :param plot_data: (bool) If True: plots the process of chirp detection. 

77 :return:chirp_time: (array) array of times (in sec) where chirps have been detected. 

78 """ 

79 if not hasattr(fundamentals, '__len__'): 

80 if not hasattr(fishlist, '__len__'): 

81 print('fishlist or fundamentals missing as argument !!!') 

82 quit() 

83 else: 

84 fundamentals = [] 

85 for fish in fishlist: 

86 if fish[0][1] > min_power: 

87 fundamentals.append(fish[0][0]) 

88 

89 chirp_time = np.array([]) 

90 chirp_freq = np.array([]) 

91 

92 for enu, fundamental in enumerate(fundamentals): 

93 # extract power of only the part of the spectrum that has to be analysied for each fundamental and get the peak 

94 # power of every piont in time. 

95 power = np.max(spectrum[(freqs >= fundamental - freq_tolerance) & (freqs <= fundamental + freq_tolerance)], axis=0) 

96 power_above = np.max(spectrum[(freqs >= fundamental + 50.0 -freq_tolerance) & (freqs <= fundamental + 50.0 + freq_tolerance)], axis=0) 

97 #power = np.mean(spectrum[(freqs >= fundamental - freq_tolerance) & (freqs <= fundamental + freq_tolerance)], axis=0) 

98 # calculate the slope by calculating the difference in the power 

99 power_diff = np.diff(power) 

100 

101 # peak detection in the power_diff to detect drops in power indicating chrips 

102 threshold = std_threshold(power_diff) 

103 peaks, troughs = detect_peaks(power_diff, threshold) 

104 troughs, peaks = trim_to_peak(troughs, peaks) # reversed troughs and peaks in output and input to get trim_to_troughs 

105 

106 # exclude peaks and troughs with to much time diff to be a chirp 

107 # ToDO: not nice !!! 

108 peaks = peaks[(troughs - peaks) < chirp_th] 

109 troughs = troughs[(troughs - peaks) < chirp_th] 

110 

111 if len(troughs) > 0: 

112 # chirps times defined as the mean time between the troughs and peaks 

113 chirp_time_idx = np.mean([troughs, peaks], axis=0) 

114 

115 # exclude detected chirps if the powervalue doesn't drop far enought 

116 chirp_time_idx = true_chirp_power_drop(chirp_time_idx, power) 

117 

118 # chirp_time_idx = true_chirp_power_rise_above(chirp_time_idx, power_above) 

119 # add times of detected chirps to the list. 

120 chirp_time = np.concatenate((chirp_time, np.array([time[int(i)] for i in chirp_time_idx]))) 

121 chirp_freq = np.concatenate((chirp_freq, np.array(fundamental* np.ones(len(chirp_time_idx))))) 

122 

123 else: 

124 chirp_time = np.array([]) 

125 chirp_freq = np.array([]) 

126 

127 if plot_data_func: 

128 plot_data_func(enu, chirp_time, time, power, power_above, power_diff, fundamental) 

129 

130 return chirp_time, chirp_freq 

131 

132 

133def chirp_detection_plot(enu, chirp_time, time, power, power2, power_diff, fundamental): 

134 """ 

135 plots the process of chirp detection. 

136 

137 :param enu: (int) indication which fish in list is processed. 

138 :param chirp_time: (array) timestamps when chirps have been detected. 

139 :param time: (array) time array. 

140 :param power: (array) power of a certain frequency band. 

141 :param power_diff: (array) slope of the power array. 

142 :param fundamental: (float) fundamental frequency around which the algorithm looked for chirps. 

143 """ 

144 try: 

145 ax 

146 except NameError: 

147 fig, ax = plt.subplots() 

148 colors = ['r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue'] 

149 # if enu == 0: 

150 # fig, ax = plt.subplots() 

151 # colors = ['r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue', 'r', 'g', 'k', 'blue'] 

152 ax.plot(chirp_time, np.zeros(len(chirp_time)), 'o', markersize=10, color=colors[enu], alpha=0.8, label='chirps') 

153 ax.set_xlabel('time in sec') 

154 ax.set_ylabel('power') 

155 

156 ax.plot(time, power, colors[enu], marker='.', label='%.1f Hz' % fundamental) 

157 ax.plot(time, power2, colors[enu+1], label='%.1f Hz' % (fundamental+50.0)) 

158 ax.plot(time[:len(power_diff)], power_diff, colors[enu], label='slope') 

159 ax.legend(loc='upper right', bbox_to_anchor=(1, 1), frameon=False) 

160 

161 

162def chirp_analysis(data, samplerate): 

163 """ 

164 Performs all steps to detect chirps in a given dataset. This includes spectrogram calculation, fish detection and 

165 analysing of specific frequency bands. 

166 For further documentation see functions chirp_spectrogram() and chirp_detection(). 

167 !!! recommended for short recordings (up to 5 min) where only the chirp times shall be extracted !!! 

168 

169 :param data: (array) data. 

170 :param samplerate: (float) smaplerate of the data. 

171 :param min_power: (float) minimal power of the fish fundamental to include this fish in chirp detection. 

172 """ 

173 freqs, time, spectrum = spectrogram(data, samplerate, freq_resolution=2., overlap_frac=0.95) 

174 

175 power = np.mean(spectrum, axis=1) # spectrum[:, t0:t1] to only let spectrum of certain time.... 

176 

177 fishlist = harmonic_groups(freqs, power)[0] 

178 

179 chirp_time, chirp_freq = chirp_detection(spectrum, freqs, time, fishlist, plot_data_func=chirp_detection_plot) 

180 

181 plt.show() 

182 

183 return chirp_time, chirp_freq 

184 

185 

186if __name__ == '__main__': 

187 ### 

188 # If you want to test the code I propose to use the file '60427L05.WAV' of the transect 

189 # '2016_04_27__downstream_stonewall_at_pool' made in colombia, 2016. 

190 ### 

191 import sys 

192 import matplotlib.pyplot as plt 

193 from thunderlab.dataloader import load_data 

194 

195 data_file = sys.argv[1] 

196 raw_data, samplerate, unit, amax = load_data(data_file) 

197 

198 chirp_time, chirp_freq = chirp_analysis(raw_data[:,0], samplerate) 

199 

200 # power = np.mean(spectrum[:, t:t + nffts_per_psd], axis=1)