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

1""" 

2Check for pulse-type weakly electric fish 

3 

4Check whether a pulse-type or a wave-type weakly electric fish is present in a recording. 

5 

6## Main functions 

7 

8- `check_pulse()`: checks for pulse-type fish. 

9 

10## Configuration parameter 

11 

12- `add_check_pulse_config()`: add parameters for `check_pulse()` to configuration. 

13- `check_pulse_args()`: retrieve parameters for `check_pulse()` from configuration. 

14 

15""" 

16 

17import numpy as np 

18from thunderlab.eventdetection import percentile_threshold, detect_peaks, trim 

19 

20 

21def check_pulse(data, sem, samplerate, thresh_fac=0.8, percentile=0.0, 

22 sem_fac=0.05, pulse_thresh=0.15, verbose=0): 

23 """Detects if a fish is pulse- or wave-type based on the proportion of the time distance 

24 between a peak and its following trough, relative to the time between consecutive peaks. 

25 

26 

27 Parameters 

28 ---------- 

29 data: 1-D array 

30 The data to be analyzed. 

31 sem: 1-D array or None 

32 Standard error of the mean corresponding to data. 

33 samplerate: float 

34 Sampling rate of the data in Hertz. 

35 percentile: float 

36 The interpercentile range is computed at percentile and 100.0-percentile. 

37 If zero take maxmimum minus minimum. 

38 thresh_fac: float 

39 The threshold for peak detection is the inter-percentile-range 

40 multiplied by this factor. 

41 sem_fac: float 

42 Base analysis on peaks and troughs with the corresponding standard error of the mean 

43 of the data below `sem_fac` times the interpercentile range. 

44 pulse_thresh: float 

45 A positive number setting the minimum distance between peaks and troughs. 

46 verbose: int 

47 If > 1, print information in the command line. 

48 

49 Returns 

50 ------- 

51 pulse_fish: bool 

52 True if algorithm suggests a pulse-type fish. 

53 peak_ratio: float 

54 Returns a float between 0. and 1. which gives the proportion of peak-2-trough, 

55 from peak-2-peak time distance, i.e. pulse width relative to pulse interval. 

56 """ 

57 

58 def ratio(peak_idx, trough_idx): 

59 if len(peak_idx) < 2: 

60 return 1.0 if len(peak_idx)+len(trough_idx) < 1 else 0.0 

61 # ratio of peak-to-trough to peak-to-peak time distances: 

62 ratios = np.abs((trough_idx - peak_idx))[:-1].astype(float) / np.diff(peak_idx) 

63 # fix for cases where trough of eod comes before peak: 

64 ratios[ratios > 0.5] = 1.0 - ratios[ratios > 0.5] 

65 return np.median(ratios) 

66 

67 # threshold for peak detection: 

68 pp_ampl = percentile_threshold(data, thresh_fac=1.0, percentile=percentile) 

69 threshold = thresh_fac*pp_ampl 

70 if verbose > 1: 

71 print(' check_pulse(): amplitude threshold for peak detection is %g' % threshold) 

72 # detect large peaks and troughs: 

73 peak_idx, trough_idx = detect_peaks(data, threshold) 

74 # XXX check here for less than 2 peaks + troughs XXX 

75 if verbose > 1: 

76 print(' check_pulse(): detected %d peaks %d troughs' % 

77 (len(peak_idx), len(trough_idx))) 

78 peak_idx, trough_idx = trim(peak_idx, trough_idx) 

79 if verbose > 1: 

80 print(' check_pulse(): trim %d peaks %d troughs' % 

81 (len(peak_idx), len(trough_idx))) 

82 # take only peaks and troughs of reliable sem: 

83 if sem is not None: 

84 athresh = sem_fac*pp_ampl 

85 sel = (sem[peak_idx]<athresh) & (sem[trough_idx]<athresh) 

86 peak_idx = peak_idx[sel] 

87 trough_idx = trough_idx[sel] 

88 if verbose > 1: 

89 print(' check_pulse(): sem left %d peaks %d troughs' % 

90 (len(peak_idx), len(trough_idx))) 

91 # compute ratios: 

92 peak_ratio = np.mean([ratio(peak_idx, trough_idx), ratio(trough_idx, peak_idx)]) 

93 pulse_fish = peak_ratio < pulse_thresh 

94 if verbose > 0: 

95 print(' check_pulse(): classified as %s fish: pulse-width-ratio is %.3f' % 

96 ('pulse' if pulse_fish else 'wave', peak_ratio)) 

97 return pulse_fish, peak_ratio 

98 

99 

100def add_check_pulse_config(cfg, thresh_fac=0.4, percentile=0.0, 

101 sem_fac=0.05, pulse_thresh=0.1): 

102 """ Add parameter needed for `check_pulse()` as a new section to a configuration. 

103 

104 Parameters 

105 ---------- 

106 cfg: ConfigFile 

107 The configuration. 

108 See check_pulse() for details on the remaining arguments. 

109 """ 

110 

111 cfg.add_section('Classify waveforms:') 

112 cfg.add('pulseWidthPercentile', percentile, '%', 'Threshold for detecing peaks is based on interpercentile range. If zero use maximum minus minimum amplitude.') 

113 cfg.add('pulseWidthThresholdFactor', thresh_fac, '', 'The threshold for peak detection is this factor multiplied with the interpercentile range. Should be smaller than 0.5 if only a single pulse is contained in the data.') 

114 cfg.add('pulseWidthSEMFactor', sem_fac, '', 'Peaks and troughs where the standarad error of the mean of the data is larger than the interpercentile range times this factor are discarded from the analysis.') 

115 cfg.add('pulseWidthThresholdRatio', pulse_thresh, '', 'A pulsefish is detected if the width of the pulses relative to the intervals is smaller than this threshold.') 

116 

117 

118def check_pulse_args(cfg): 

119 """ Translates a configuration to the 

120 respective parameter names of the function `check_pulse()`. 

121 The return value can then be passed as key-word arguments to this function. 

122 

123 Parameters 

124 ---------- 

125 cfg: ConfigFile 

126 The configuration. 

127 

128 Returns 

129 ------- 

130 a: dict 

131 Dictionary with names of arguments of the check_pulse() function 

132 and their values as supplied by `cfg`. 

133 """ 

134 a = cfg.map({'thresh_fac': 'pulseWidthThresholdFactor', 

135 'percentile': 'pulseWidthPercentile', 

136 'sem_fac': 'pulseWidthSEMFactor', 

137 'pulse_thresh': 'pulseWidthThresholdRatio'}) 

138 return a 

139 

140 

141if __name__ == "__main__": 

142 import sys 

143 import matplotlib.pyplot as plt 

144 from .bestwindow import best_window 

145 from .fakefish import wavefish_eods, pulsefish_eods 

146 

147 print("\nChecking checkpulse module ...\n") 

148 

149 # generate data: 

150 rate = 44100.0 

151 if len(sys.argv) < 2: 

152 data = wavefish_eods('Eigenmannia', 80.0, rate, 8.0) 

153 elif sys.argv[1] == '-w': 

154 data = wavefish_eods('Alepto', 600.0, rate, 8.0) 

155 elif sys.argv[1] == '-m': 

156 data = pulsefish_eods('monophasic', 80.0, rate, 8.0) 

157 elif sys.argv[1] == '-b': 

158 data = pulsefish_eods('biphasic', 80.0, rate, 8.0) 

159 elif sys.argv[1] == '-t': 

160 data = pulsefish_eods('triphasic', 80.0, rate, 8.0) 

161 else: # load data given by the user 

162 from thunderlab.dataloader import load_data 

163 

164 file_path = sys.argv[1] 

165 print("loading %s ...\n" % file_path) 

166 rawdata, rate, unit, amax = load_data(sys.argv[1]) 

167 data, _ = best_window(rawdata[:,0], rate) 

168 

169 # run pulse-width-based detector: 

170 pulse_fish, ratio = check_pulse(data, None, rate)