-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy paththreshold.py
160 lines (112 loc) · 4.34 KB
/
threshold.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
import math
import numpy as np
__author__="Andrew"
__date__ ="$Jan 31, 2011 11:04:09 PM$"
from numpy import *
def complement(a, b):
'''zero biased complement function'''
return b-a-1
def fall_threshold (trace, lo, hi):
'''
Wrapper for rise_threshold
'''
l = len(trace)
#send reversed trace
a,b,c,d,e,f= rise_threshold(trace[::-1],lo,hi)
return a, complement (b,l), c, complement (d,l), complement (e,l), complement (f,l)
def rise_threshold(trace, lo, hi):
''' find the threshold crossing points for a trace
Arguments:
-- trace: a 1-D numpy array containing a single peak
-- lo : the low threshold as a decimal between 0 and 1
-- hi : the high threshold as above
Returns: -- nearest low value,
-- index of nearest low value
-- nearest high value
-- index of nearest high value
-- rise times between thresholds
interpolates between bracketing points if no direct hit
'''
max = trace.max()
min = trace.min()
l = len(trace)
if math.fabs(max) < math.fabs(min) :
peak_negative = True
#make trace positive to simplify maths
trace = - trace
max = trace.max()
min = trace.min()
else:
peak_negative = False
max_idx = trace.argmax()
print "peak at index %i of %i elements" %(max_idx, l)
# find index of max point
#slice array and use only rising phase
if max_idx < l:
rise_trace = trace[:max_idx+1:] #include max point
#print rise_trace
else:
rise_trace = trace
rise_l = len(rise_trace)
range = max - min
lo_threshold = range * lo
hi_threshold = range * hi
#[::-1] reverses the array, counting down from peak. Ensures closest value to peak
near_lo, near_lo_idx = find_nearest (rise_trace[::-1],lo_threshold)
near_hi, near_hi_idx = find_nearest (rise_trace[::-1],hi_threshold)
near_hi_idx = rise_l - near_hi_idx - 1 #array zero biassed but len not
near_lo_idx = rise_l - near_lo_idx - 1 #array zero biassed but len not
#interpolate
if near_hi > hi_threshold:
left = rise_trace[near_hi_idx - 1]
right = near_hi
#(point-threshold) / spread = time to sub from point to get threshold cross
hi_crossing = near_hi_idx - float(right - hi_threshold)/(right-left)
elif near_hi < hi_threshold:
left = near_hi
right = rise_trace[near_hi_idx+1]
hi_crossing = near_hi_idx + float(hi_threshold-left)/(right-left)
elif near_hi == hi_threshold:
hi_crossing = near_hi_idx
if near_lo > lo_threshold:
left = rise_trace[near_lo_idx - 1]
right = near_lo
lo_crossing = near_lo_idx- float(right - lo_threshold)/(right-left)
elif near_lo < lo_threshold:
left = near_lo
right = rise_trace[near_lo_idx+1]
lo_crossing = near_lo_idx + float(lo_threshold-left)/(right-left)
elif near_lo == lo_threshold:
lo_crossing = near_lo_idx
if peak_negative:
#re-invert values that are returned
near_lo = - near_lo
near_hi = - near_hi
return near_lo, near_lo_idx, near_hi, near_hi_idx, lo_crossing,hi_crossing
def find_nearest(array, value):
"""the closest match between the value and the values in the array"""
idx = (np.abs(array - value)).argmin()
return array[idx], idx
if __name__ == "__main__":
low = 0.1
high = 0.9
#find max and min in data
data = np.array([0,0,0,0,0,0.01,0.075,.221,.9,1,0.9,1.,1,0,0,0,0,0])
data_set = []
fac = [1,20,-10,.5]
for n in fac:
data_set.append(data*n)
for da in data_set:
a,b,c,d,e,f = rise_threshold (da,low,high)
aa,bb,cc,dd,ee,ff = fall_threshold (da,low,high)
#print a,b,c,d
print da
print "%i-%i%% rise time = %f time units" %(int(low*100),int(high*100),e-f)
print "%i-%i%% fall time = %f time units" %(int(low*100),int(high*100),ee-ff)
print a,da[b],c,da[d]
print aa,da[bb],cc,da[dd]
a,b,c,d,e,f = rise_threshold (data,.5,.6) #.6 is not used
aa,bb,cc,dd,ee,ff = fall_threshold (data,.5,.6)
print a,b,c,d,e,f
print aa,bb,cc,dd,ee,ff
print "FWHM = %f" %(ee-e)