source: fact/tools/pyscripts/sandbox/dneise/spike_studies/drs_spikes.py@ 20115

Last change on this file since 20115 was 13632, checked in by neise, 13 years ago
initial commit of my tests and so on
  • Property svn:executable set to *
File size: 3.5 KB
Line 
1#!/usr/bin/python -tt
2#
3# Werner Lustermann
4# ETH Zurich
5#
6import numpy as np
7
8import fir_filter as fir
9
10class DRSSpikes(object):
11 """ remove spikes (single or double false readings) from DRS4 data
12 Strategy:
13 * filter the data, removing the signal, thus spike(s) are clearly visible
14 * search single and double spikes
15 * replace the spike by a value derived from the neighbors
16
17 """
18
19 def __init__(self, threshold=7.,
20 single_pattern=np.array( [-0.5, 1.0, -0.5]) ,
21 double_pattern=np.array([-1., 1., 1., -1.]),
22 user_action=lambda DRSSpikes_instance: None,
23 debug = False):
24 """ initialize spike filter
25 template_single: template of a single slice spike
26 template_double: template of a two slice spike
27
28 """
29
30 self.threshold = threshold
31 self.single_pattern = single_pattern * threshold
32 self.double_pattern = double_pattern * threshold
33
34 self.remove_signal = fir.RemoveSignal()
35
36 self.user_action = user_action
37 self.debug = debug
38
39 def __call__(self, data):
40
41 self.row, self.col = data.shape
42 indicator = self.remove_signal(data)
43 a = indicator.flatten()
44 singles = []
45 doubles = []
46
47 # a spike in the first or last channel is considered as a filter artefact
48 candidates = np.where(a[1:-2] > self.threshold)
49 # candidates = np.where(a[1:1022] > self.threshold)
50 cc = candidates[0]
51 #print 'cc: ', cc
52 #: find single spikes
53 p = self.single_pattern * np.sign( self.single_pattern )
54 for i, can in enumerate( zip(a[cc], a[cc+1], a[cc+2]) ):
55 #print 'can : p', can, p
56 can = can * np.sign(self.single_pattern)
57 if np.all(can > p):
58 singles.append(cc[i])
59
60 #: find double spikes
61 p = self.double_pattern * np.sign( self.double_pattern )
62 for i, can in enumerate( zip(a[cc], a[cc+1], a[cc+2], a[cc+3]) ):
63 #print 'data: ', [data[0,cc[i]+k] for k in range(3)]
64 #print 'can : p', can, p
65 can = can * np.sign(self.double_pattern)
66 if np.all(can > p):
67 doubles.append(cc[i])
68
69 if self.debug:
70 print 'singles: ', singles
71 print 'doubles: ', doubles
72
73 # make relevant information a member,
74 self.singles = singles
75 self.doubles = doubles
76 self.indicator = a
77 self.candidates = cc
78
79 self.user_action(self)
80
81 data = self.remove_single_spikes(singles, data)
82 data = self.remove_double_spikes(doubles, data)
83 return data
84
85 def remove_single_spikes(self, singles, data):
86 data = data.flatten()
87 for spike in singles:
88 data[spike] = (data[spike-1] + data[spike+1]) / 2.
89 return data.reshape(self.row, self.col)
90
91 def remove_double_spikes(self, doubles, data):
92 data = data.flatten()
93 for spike in doubles:
94 data[spike:spike+2] = (data[spike-1] + data[spike+2]) / 2.
95 return data.reshape(self.row, self.col)
96
97
98def _test():
99
100 a = np.ones((3,12)) * 3.
101 a[0,3] = 7.
102 a[1,7] = 14.
103 a[1,8] = 14.
104 a[2,4] = 50.
105 a[2,5] = 45.
106 a[2,8] = 20.
107
108 print a
109
110 SpikeRemover = DRSSpikes(3., debug=True)
111 print 'single spike pattern ', SpikeRemover.single_pattern
112 print 'double spike pattern ', SpikeRemover.double_pattern
113 afilt = SpikeRemover(a)
114 print afilt
115
116if __name__ == '__main__':
117 """ test the class """
118 _test()
119
120
Note: See TracBrowser for help on using the repository browser.