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 

18 

19from thunderlab.eventdetection import percentile_threshold, detect_peaks, trim 

20 

21 

22def check_pulse(data, sem, rate, thresh_fac=0.8, percentile=0.0, 

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

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

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

26 

27 

28 Parameters 

29 ---------- 

30 data: 1-D array 

31 The data to be analyzed. 

32 sem: 1-D array or None 

33 Standard error of the mean corresponding to data. 

34 rate: float 

35 Sampling rate of the data in Hertz. 

36 percentile: float 

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

38 If zero take maxmimum minus minimum. 

39 thresh_fac: float 

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

41 multiplied by this factor. 

42 sem_fac: float 

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

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

45 pulse_thresh: float 

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

47 verbose: int 

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

49 

50 Returns 

51 ------- 

52 pulse_fish: bool 

53 True if algorithm suggests a pulse-type fish. 

54 peak_ratio: float 

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

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

57 """ 

58 

59 def ratio(peak_idx, trough_idx): 

60 if len(peak_idx) < 2: 

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

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

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

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

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

66 return np.median(ratios) 

67 

68 # threshold for peak detection: 

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

70 threshold = thresh_fac*pp_ampl 

71 if verbose > 1: 

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

73 # detect large peaks and troughs: 

74 peak_idx, trough_idx = detect_peaks(data, threshold) 

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

76 if verbose > 1: 

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

78 (len(peak_idx), len(trough_idx))) 

79 peak_idx, trough_idx = trim(peak_idx, trough_idx) 

80 if verbose > 1: 

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

82 (len(peak_idx), len(trough_idx))) 

83 # take only peaks and troughs of reliable sem: 

84 if sem is not None: 

85 athresh = sem_fac*pp_ampl 

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

87 peak_idx = peak_idx[sel] 

88 trough_idx = trough_idx[sel] 

89 if verbose > 1: 

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

91 (len(peak_idx), len(trough_idx))) 

92 # compute ratios: 

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

94 pulse_fish = peak_ratio < pulse_thresh 

95 if verbose > 0: 

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

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

98 return pulse_fish, peak_ratio 

99 

100 

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

102 sem_fac=0.05, pulse_thresh=0.1): 

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

104 

105 Parameters 

106 ---------- 

107 cfg: ConfigFile 

108 The configuration. 

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

110 """ 

111 

112 cfg.add_section('Classify waveforms:') 

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

114 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.') 

115 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.') 

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

117 

118 

119def check_pulse_args(cfg): 

120 """ Translates a configuration to the 

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

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

123 

124 Parameters 

125 ---------- 

126 cfg: ConfigFile 

127 The configuration. 

128 

129 Returns 

130 ------- 

131 a: dict 

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

133 and their values as supplied by `cfg`. 

134 """ 

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

136 'percentile': 'pulseWidthPercentile', 

137 'sem_fac': 'pulseWidthSEMFactor', 

138 'pulse_thresh': 'pulseWidthThresholdRatio'}) 

139 return a 

140 

141 

142if __name__ == "__main__": 

143 import sys 

144 import matplotlib.pyplot as plt 

145 from .bestwindow import best_window 

146 from .fakefish import wavefish_eods, pulsefish_eods 

147 

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

149 

150 # generate data: 

151 rate = 44100.0 

152 if len(sys.argv) < 2: 

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

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

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

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

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

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

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

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

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

162 else: # load data given by the user 

163 from thunderlab.dataloader import load_data 

164 

165 file_path = sys.argv[1] 

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

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

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

169 

170 # run pulse-width-based detector: 

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