##// END OF EJS Templates
jroproc_spectra.py: math module is not used anymore
Ivan Valdez -
r722:6a6dc64e1b2f
parent child
Show More
@@ -1,878 +1,877
1 import numpy
1 import numpy
2 import math
3
2
4 from jroproc_base import ProcessingUnit, Operation
3 from jroproc_base import ProcessingUnit, Operation
5 from schainpy.model.data.jrodata import Spectra
4 from schainpy.model.data.jrodata import Spectra
6 from schainpy.model.data.jrodata import hildebrand_sekhon
5 from schainpy.model.data.jrodata import hildebrand_sekhon
7
6
8 class SpectraProc(ProcessingUnit):
7 class SpectraProc(ProcessingUnit):
9
8
10 def __init__(self):
9 def __init__(self):
11
10
12 ProcessingUnit.__init__(self)
11 ProcessingUnit.__init__(self)
13
12
14 self.buffer = None
13 self.buffer = None
15 self.firstdatatime = None
14 self.firstdatatime = None
16 self.profIndex = 0
15 self.profIndex = 0
17 self.dataOut = Spectra()
16 self.dataOut = Spectra()
18 self.id_min = None
17 self.id_min = None
19 self.id_max = None
18 self.id_max = None
20
19
21 def __updateSpecFromVoltage(self):
20 def __updateSpecFromVoltage(self):
22
21
23 self.dataOut.timeZone = self.dataIn.timeZone
22 self.dataOut.timeZone = self.dataIn.timeZone
24 self.dataOut.dstFlag = self.dataIn.dstFlag
23 self.dataOut.dstFlag = self.dataIn.dstFlag
25 self.dataOut.errorCount = self.dataIn.errorCount
24 self.dataOut.errorCount = self.dataIn.errorCount
26 self.dataOut.useLocalTime = self.dataIn.useLocalTime
25 self.dataOut.useLocalTime = self.dataIn.useLocalTime
27
26
28 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
27 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
29 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
28 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
30 self.dataOut.channelList = self.dataIn.channelList
29 self.dataOut.channelList = self.dataIn.channelList
31 self.dataOut.heightList = self.dataIn.heightList
30 self.dataOut.heightList = self.dataIn.heightList
32 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
31 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
33
32
34 self.dataOut.nBaud = self.dataIn.nBaud
33 self.dataOut.nBaud = self.dataIn.nBaud
35 self.dataOut.nCode = self.dataIn.nCode
34 self.dataOut.nCode = self.dataIn.nCode
36 self.dataOut.code = self.dataIn.code
35 self.dataOut.code = self.dataIn.code
37 self.dataOut.nProfiles = self.dataOut.nFFTPoints
36 self.dataOut.nProfiles = self.dataOut.nFFTPoints
38
37
39 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
38 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
40 self.dataOut.utctime = self.firstdatatime
39 self.dataOut.utctime = self.firstdatatime
41 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
40 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
42 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
41 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
43 self.dataOut.flagShiftFFT = False
42 self.dataOut.flagShiftFFT = False
44
43
45 self.dataOut.nCohInt = self.dataIn.nCohInt
44 self.dataOut.nCohInt = self.dataIn.nCohInt
46 self.dataOut.nIncohInt = 1
45 self.dataOut.nIncohInt = 1
47
46
48 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
47 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
49
48
50 self.dataOut.frequency = self.dataIn.frequency
49 self.dataOut.frequency = self.dataIn.frequency
51 self.dataOut.realtime = self.dataIn.realtime
50 self.dataOut.realtime = self.dataIn.realtime
52
51
53 self.dataOut.azimuth = self.dataIn.azimuth
52 self.dataOut.azimuth = self.dataIn.azimuth
54 self.dataOut.zenith = self.dataIn.zenith
53 self.dataOut.zenith = self.dataIn.zenith
55
54
56 self.dataOut.beam.codeList = self.dataIn.beam.codeList
55 self.dataOut.beam.codeList = self.dataIn.beam.codeList
57 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
56 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
58 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
57 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
59
58
60 def __getFft(self):
59 def __getFft(self):
61 """
60 """
62 Convierte valores de Voltaje a Spectra
61 Convierte valores de Voltaje a Spectra
63
62
64 Affected:
63 Affected:
65 self.dataOut.data_spc
64 self.dataOut.data_spc
66 self.dataOut.data_cspc
65 self.dataOut.data_cspc
67 self.dataOut.data_dc
66 self.dataOut.data_dc
68 self.dataOut.heightList
67 self.dataOut.heightList
69 self.profIndex
68 self.profIndex
70 self.buffer
69 self.buffer
71 self.dataOut.flagNoData
70 self.dataOut.flagNoData
72 """
71 """
73 fft_volt = numpy.fft.fft(self.buffer,n=self.dataOut.nFFTPoints,axis=1)
72 fft_volt = numpy.fft.fft(self.buffer,n=self.dataOut.nFFTPoints,axis=1)
74 fft_volt = fft_volt.astype(numpy.dtype('complex'))
73 fft_volt = fft_volt.astype(numpy.dtype('complex'))
75 dc = fft_volt[:,0,:]
74 dc = fft_volt[:,0,:]
76
75
77 #calculo de self-spectra
76 #calculo de self-spectra
78 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
77 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
79 spc = fft_volt * numpy.conjugate(fft_volt)
78 spc = fft_volt * numpy.conjugate(fft_volt)
80 spc = spc.real
79 spc = spc.real
81
80
82 blocksize = 0
81 blocksize = 0
83 blocksize += dc.size
82 blocksize += dc.size
84 blocksize += spc.size
83 blocksize += spc.size
85
84
86 cspc = None
85 cspc = None
87 pairIndex = 0
86 pairIndex = 0
88 if self.dataOut.pairsList != None:
87 if self.dataOut.pairsList != None:
89 #calculo de cross-spectra
88 #calculo de cross-spectra
90 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
89 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
91 for pair in self.dataOut.pairsList:
90 for pair in self.dataOut.pairsList:
92 if pair[0] not in self.dataOut.channelList:
91 if pair[0] not in self.dataOut.channelList:
93 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
92 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
94 if pair[1] not in self.dataOut.channelList:
93 if pair[1] not in self.dataOut.channelList:
95 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
94 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
96
95
97 cspc[pairIndex,:,:] = fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:])
96 cspc[pairIndex,:,:] = fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:])
98 pairIndex += 1
97 pairIndex += 1
99 blocksize += cspc.size
98 blocksize += cspc.size
100
99
101 self.dataOut.data_spc = spc
100 self.dataOut.data_spc = spc
102 self.dataOut.data_cspc = cspc
101 self.dataOut.data_cspc = cspc
103 self.dataOut.data_dc = dc
102 self.dataOut.data_dc = dc
104 self.dataOut.blockSize = blocksize
103 self.dataOut.blockSize = blocksize
105 self.dataOut.flagShiftFFT = True
104 self.dataOut.flagShiftFFT = True
106
105
107 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None):
106 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None):
108
107
109 self.dataOut.flagNoData = True
108 self.dataOut.flagNoData = True
110
109
111 if self.dataIn.type == "Spectra":
110 if self.dataIn.type == "Spectra":
112 self.dataOut.copy(self.dataIn)
111 self.dataOut.copy(self.dataIn)
113 return True
112 return True
114
113
115 if self.dataIn.type == "Voltage":
114 if self.dataIn.type == "Voltage":
116
115
117 if nFFTPoints == None:
116 if nFFTPoints == None:
118 raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
117 raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
119
118
120 if nProfiles == None:
119 if nProfiles == None:
121 nProfiles = nFFTPoints
120 nProfiles = nFFTPoints
122
121
123 if ippFactor == None:
122 if ippFactor == None:
124 ippFactor = 1
123 ippFactor = 1
125
124
126 self.dataOut.ippFactor = ippFactor
125 self.dataOut.ippFactor = ippFactor
127
126
128 self.dataOut.nFFTPoints = nFFTPoints
127 self.dataOut.nFFTPoints = nFFTPoints
129 self.dataOut.pairsList = pairsList
128 self.dataOut.pairsList = pairsList
130
129
131 if self.buffer is None:
130 if self.buffer is None:
132 self.buffer = numpy.zeros( (self.dataIn.nChannels,
131 self.buffer = numpy.zeros( (self.dataIn.nChannels,
133 nProfiles,
132 nProfiles,
134 self.dataIn.nHeights),
133 self.dataIn.nHeights),
135 dtype='complex')
134 dtype='complex')
136
135
137 if self.dataIn.flagDataAsBlock:
136 if self.dataIn.flagDataAsBlock:
138 #data dimension: [nChannels, nProfiles, nSamples]
137 #data dimension: [nChannels, nProfiles, nSamples]
139 nVoltProfiles = self.dataIn.data.shape[1]
138 nVoltProfiles = self.dataIn.data.shape[1]
140 nVoltProfiles = self.dataIn.nProfiles
139 nVoltProfiles = self.dataIn.nProfiles
141
140
142 if nVoltProfiles == nProfiles:
141 if nVoltProfiles == nProfiles:
143 self.buffer = self.dataIn.data.copy()
142 self.buffer = self.dataIn.data.copy()
144 self.profIndex = nVoltProfiles
143 self.profIndex = nVoltProfiles
145
144
146 elif nVoltProfiles < nProfiles:
145 elif nVoltProfiles < nProfiles:
147
146
148 if self.profIndex == 0:
147 if self.profIndex == 0:
149 self.id_min = 0
148 self.id_min = 0
150 self.id_max = nVoltProfiles
149 self.id_max = nVoltProfiles
151
150
152 self.buffer[:,self.id_min:self.id_max,:] = self.dataIn.data
151 self.buffer[:,self.id_min:self.id_max,:] = self.dataIn.data
153 self.profIndex += nVoltProfiles
152 self.profIndex += nVoltProfiles
154 self.id_min += nVoltProfiles
153 self.id_min += nVoltProfiles
155 self.id_max += nVoltProfiles
154 self.id_max += nVoltProfiles
156 else:
155 else:
157 raise ValueError, "The type object %s has %d profiles, it should be equal to %d profiles"%(self.dataIn.type,self.dataIn.data.shape[1],nProfiles)
156 raise ValueError, "The type object %s has %d profiles, it should be equal to %d profiles"%(self.dataIn.type,self.dataIn.data.shape[1],nProfiles)
158 self.dataOut.flagNoData = True
157 self.dataOut.flagNoData = True
159 return 0
158 return 0
160 else:
159 else:
161 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
160 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
162 self.profIndex += 1
161 self.profIndex += 1
163
162
164 if self.firstdatatime == None:
163 if self.firstdatatime == None:
165 self.firstdatatime = self.dataIn.utctime
164 self.firstdatatime = self.dataIn.utctime
166
165
167 if self.profIndex == nProfiles:
166 if self.profIndex == nProfiles:
168 self.__updateSpecFromVoltage()
167 self.__updateSpecFromVoltage()
169 self.__getFft()
168 self.__getFft()
170
169
171 self.dataOut.flagNoData = False
170 self.dataOut.flagNoData = False
172 self.firstdatatime = None
171 self.firstdatatime = None
173 self.profIndex = 0
172 self.profIndex = 0
174
173
175 return True
174 return True
176
175
177 raise ValueError, "The type of input object '%s' is not valid"%(self.dataIn.type)
176 raise ValueError, "The type of input object '%s' is not valid"%(self.dataIn.type)
178
177
179 def __selectPairs(self, channelList=None):
178 def __selectPairs(self, channelList=None):
180
179
181 if channelList == None:
180 if channelList == None:
182 return
181 return
183
182
184 pairsIndexListSelected = []
183 pairsIndexListSelected = []
185 for pairIndex in self.dataOut.pairsIndexList:
184 for pairIndex in self.dataOut.pairsIndexList:
186 #First pair
185 #First pair
187 if self.dataOut.pairsList[pairIndex][0] not in channelList:
186 if self.dataOut.pairsList[pairIndex][0] not in channelList:
188 continue
187 continue
189 #Second pair
188 #Second pair
190 if self.dataOut.pairsList[pairIndex][1] not in channelList:
189 if self.dataOut.pairsList[pairIndex][1] not in channelList:
191 continue
190 continue
192
191
193 pairsIndexListSelected.append(pairIndex)
192 pairsIndexListSelected.append(pairIndex)
194
193
195 if not pairsIndexListSelected:
194 if not pairsIndexListSelected:
196 self.dataOut.data_cspc = None
195 self.dataOut.data_cspc = None
197 self.dataOut.pairsList = []
196 self.dataOut.pairsList = []
198 return
197 return
199
198
200 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
199 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
201 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
200 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
202
201
203 return
202 return
204
203
205 def selectChannels(self, channelList):
204 def selectChannels(self, channelList):
206
205
207 channelIndexList = []
206 channelIndexList = []
208
207
209 for channel in channelList:
208 for channel in channelList:
210 if channel not in self.dataOut.channelList:
209 if channel not in self.dataOut.channelList:
211 raise ValueError, "Error selecting channels: The value %d in channelList is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList))
210 raise ValueError, "Error selecting channels: The value %d in channelList is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList))
212
211
213 index = self.dataOut.channelList.index(channel)
212 index = self.dataOut.channelList.index(channel)
214 channelIndexList.append(index)
213 channelIndexList.append(index)
215
214
216 self.selectChannelsByIndex(channelIndexList)
215 self.selectChannelsByIndex(channelIndexList)
217
216
218 def selectChannelsByIndex(self, channelIndexList):
217 def selectChannelsByIndex(self, channelIndexList):
219 """
218 """
220 Selecciona un bloque de datos en base a canales segun el channelIndexList
219 Selecciona un bloque de datos en base a canales segun el channelIndexList
221
220
222 Input:
221 Input:
223 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
222 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
224
223
225 Affected:
224 Affected:
226 self.dataOut.data_spc
225 self.dataOut.data_spc
227 self.dataOut.channelIndexList
226 self.dataOut.channelIndexList
228 self.dataOut.nChannels
227 self.dataOut.nChannels
229
228
230 Return:
229 Return:
231 None
230 None
232 """
231 """
233
232
234 for channelIndex in channelIndexList:
233 for channelIndex in channelIndexList:
235 if channelIndex not in self.dataOut.channelIndexList:
234 if channelIndex not in self.dataOut.channelIndexList:
236 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList)
235 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList)
237
236
238 # nChannels = len(channelIndexList)
237 # nChannels = len(channelIndexList)
239
238
240 data_spc = self.dataOut.data_spc[channelIndexList,:]
239 data_spc = self.dataOut.data_spc[channelIndexList,:]
241 data_dc = self.dataOut.data_dc[channelIndexList,:]
240 data_dc = self.dataOut.data_dc[channelIndexList,:]
242
241
243 self.dataOut.data_spc = data_spc
242 self.dataOut.data_spc = data_spc
244 self.dataOut.data_dc = data_dc
243 self.dataOut.data_dc = data_dc
245
244
246 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
245 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
247 # self.dataOut.nChannels = nChannels
246 # self.dataOut.nChannels = nChannels
248
247
249 self.__selectPairs(self.dataOut.channelList)
248 self.__selectPairs(self.dataOut.channelList)
250
249
251 return 1
250 return 1
252
251
253 def selectHeights(self, minHei, maxHei):
252 def selectHeights(self, minHei, maxHei):
254 """
253 """
255 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
254 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
256 minHei <= height <= maxHei
255 minHei <= height <= maxHei
257
256
258 Input:
257 Input:
259 minHei : valor minimo de altura a considerar
258 minHei : valor minimo de altura a considerar
260 maxHei : valor maximo de altura a considerar
259 maxHei : valor maximo de altura a considerar
261
260
262 Affected:
261 Affected:
263 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
262 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
264
263
265 Return:
264 Return:
266 1 si el metodo se ejecuto con exito caso contrario devuelve 0
265 1 si el metodo se ejecuto con exito caso contrario devuelve 0
267 """
266 """
268
267
269 if (minHei > maxHei):
268 if (minHei > maxHei):
270 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
269 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
271
270
272 if (minHei < self.dataOut.heightList[0]):
271 if (minHei < self.dataOut.heightList[0]):
273 minHei = self.dataOut.heightList[0]
272 minHei = self.dataOut.heightList[0]
274
273
275 if (maxHei > self.dataOut.heightList[-1]):
274 if (maxHei > self.dataOut.heightList[-1]):
276 maxHei = self.dataOut.heightList[-1]
275 maxHei = self.dataOut.heightList[-1]
277
276
278 minIndex = 0
277 minIndex = 0
279 maxIndex = 0
278 maxIndex = 0
280 heights = self.dataOut.heightList
279 heights = self.dataOut.heightList
281
280
282 inda = numpy.where(heights >= minHei)
281 inda = numpy.where(heights >= minHei)
283 indb = numpy.where(heights <= maxHei)
282 indb = numpy.where(heights <= maxHei)
284
283
285 try:
284 try:
286 minIndex = inda[0][0]
285 minIndex = inda[0][0]
287 except:
286 except:
288 minIndex = 0
287 minIndex = 0
289
288
290 try:
289 try:
291 maxIndex = indb[0][-1]
290 maxIndex = indb[0][-1]
292 except:
291 except:
293 maxIndex = len(heights)
292 maxIndex = len(heights)
294
293
295 self.selectHeightsByIndex(minIndex, maxIndex)
294 self.selectHeightsByIndex(minIndex, maxIndex)
296
295
297 return 1
296 return 1
298
297
299 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
298 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
300 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
299 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
301
300
302 if hei_ref != None:
301 if hei_ref != None:
303 newheis = numpy.where(self.dataOut.heightList>hei_ref)
302 newheis = numpy.where(self.dataOut.heightList>hei_ref)
304
303
305 minIndex = min(newheis[0])
304 minIndex = min(newheis[0])
306 maxIndex = max(newheis[0])
305 maxIndex = max(newheis[0])
307 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
306 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
308 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
307 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
309
308
310 # determina indices
309 # determina indices
311 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
310 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
312 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
311 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
313 beacon_dB = numpy.sort(avg_dB)[-nheis:]
312 beacon_dB = numpy.sort(avg_dB)[-nheis:]
314 beacon_heiIndexList = []
313 beacon_heiIndexList = []
315 for val in avg_dB.tolist():
314 for val in avg_dB.tolist():
316 if val >= beacon_dB[0]:
315 if val >= beacon_dB[0]:
317 beacon_heiIndexList.append(avg_dB.tolist().index(val))
316 beacon_heiIndexList.append(avg_dB.tolist().index(val))
318
317
319 #data_spc = data_spc[:,:,beacon_heiIndexList]
318 #data_spc = data_spc[:,:,beacon_heiIndexList]
320 data_cspc = None
319 data_cspc = None
321 if self.dataOut.data_cspc is not None:
320 if self.dataOut.data_cspc is not None:
322 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
321 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
323 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
322 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
324
323
325 data_dc = None
324 data_dc = None
326 if self.dataOut.data_dc is not None:
325 if self.dataOut.data_dc is not None:
327 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
326 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
328 #data_dc = data_dc[:,beacon_heiIndexList]
327 #data_dc = data_dc[:,beacon_heiIndexList]
329
328
330 self.dataOut.data_spc = data_spc
329 self.dataOut.data_spc = data_spc
331 self.dataOut.data_cspc = data_cspc
330 self.dataOut.data_cspc = data_cspc
332 self.dataOut.data_dc = data_dc
331 self.dataOut.data_dc = data_dc
333 self.dataOut.heightList = heightList
332 self.dataOut.heightList = heightList
334 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
333 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
335
334
336 return 1
335 return 1
337
336
338
337
339 def selectHeightsByIndex(self, minIndex, maxIndex):
338 def selectHeightsByIndex(self, minIndex, maxIndex):
340 """
339 """
341 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
340 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
342 minIndex <= index <= maxIndex
341 minIndex <= index <= maxIndex
343
342
344 Input:
343 Input:
345 minIndex : valor de indice minimo de altura a considerar
344 minIndex : valor de indice minimo de altura a considerar
346 maxIndex : valor de indice maximo de altura a considerar
345 maxIndex : valor de indice maximo de altura a considerar
347
346
348 Affected:
347 Affected:
349 self.dataOut.data_spc
348 self.dataOut.data_spc
350 self.dataOut.data_cspc
349 self.dataOut.data_cspc
351 self.dataOut.data_dc
350 self.dataOut.data_dc
352 self.dataOut.heightList
351 self.dataOut.heightList
353
352
354 Return:
353 Return:
355 1 si el metodo se ejecuto con exito caso contrario devuelve 0
354 1 si el metodo se ejecuto con exito caso contrario devuelve 0
356 """
355 """
357
356
358 if (minIndex < 0) or (minIndex > maxIndex):
357 if (minIndex < 0) or (minIndex > maxIndex):
359 raise ValueError, "Error selecting heights by index: Index range in (%d,%d) is not valid" % (minIndex, maxIndex)
358 raise ValueError, "Error selecting heights by index: Index range in (%d,%d) is not valid" % (minIndex, maxIndex)
360
359
361 if (maxIndex >= self.dataOut.nHeights):
360 if (maxIndex >= self.dataOut.nHeights):
362 maxIndex = self.dataOut.nHeights-1
361 maxIndex = self.dataOut.nHeights-1
363
362
364 #Spectra
363 #Spectra
365 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
364 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
366
365
367 data_cspc = None
366 data_cspc = None
368 if self.dataOut.data_cspc is not None:
367 if self.dataOut.data_cspc is not None:
369 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
368 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
370
369
371 data_dc = None
370 data_dc = None
372 if self.dataOut.data_dc is not None:
371 if self.dataOut.data_dc is not None:
373 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
372 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
374
373
375 self.dataOut.data_spc = data_spc
374 self.dataOut.data_spc = data_spc
376 self.dataOut.data_cspc = data_cspc
375 self.dataOut.data_cspc = data_cspc
377 self.dataOut.data_dc = data_dc
376 self.dataOut.data_dc = data_dc
378
377
379 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
378 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
380
379
381 return 1
380 return 1
382
381
383 def removeDC(self, mode = 2):
382 def removeDC(self, mode = 2):
384 jspectra = self.dataOut.data_spc
383 jspectra = self.dataOut.data_spc
385 jcspectra = self.dataOut.data_cspc
384 jcspectra = self.dataOut.data_cspc
386
385
387
386
388 num_chan = jspectra.shape[0]
387 num_chan = jspectra.shape[0]
389 num_hei = jspectra.shape[2]
388 num_hei = jspectra.shape[2]
390
389
391 if jcspectra is not None:
390 if jcspectra is not None:
392 jcspectraExist = True
391 jcspectraExist = True
393 num_pairs = jcspectra.shape[0]
392 num_pairs = jcspectra.shape[0]
394 else: jcspectraExist = False
393 else: jcspectraExist = False
395
394
396 freq_dc = jspectra.shape[1]/2
395 freq_dc = jspectra.shape[1]/2
397 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
396 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
398
397
399 if ind_vel[0]<0:
398 if ind_vel[0]<0:
400 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
399 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
401
400
402 if mode == 1:
401 if mode == 1:
403 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
402 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
404
403
405 if jcspectraExist:
404 if jcspectraExist:
406 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
405 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
407
406
408 if mode == 2:
407 if mode == 2:
409
408
410 vel = numpy.array([-2,-1,1,2])
409 vel = numpy.array([-2,-1,1,2])
411 xx = numpy.zeros([4,4])
410 xx = numpy.zeros([4,4])
412
411
413 for fil in range(4):
412 for fil in range(4):
414 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
413 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
415
414
416 xx_inv = numpy.linalg.inv(xx)
415 xx_inv = numpy.linalg.inv(xx)
417 xx_aux = xx_inv[0,:]
416 xx_aux = xx_inv[0,:]
418
417
419 for ich in range(num_chan):
418 for ich in range(num_chan):
420 yy = jspectra[ich,ind_vel,:]
419 yy = jspectra[ich,ind_vel,:]
421 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
420 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
422
421
423 junkid = jspectra[ich,freq_dc,:]<=0
422 junkid = jspectra[ich,freq_dc,:]<=0
424 cjunkid = sum(junkid)
423 cjunkid = sum(junkid)
425
424
426 if cjunkid.any():
425 if cjunkid.any():
427 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
426 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
428
427
429 if jcspectraExist:
428 if jcspectraExist:
430 for ip in range(num_pairs):
429 for ip in range(num_pairs):
431 yy = jcspectra[ip,ind_vel,:]
430 yy = jcspectra[ip,ind_vel,:]
432 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
431 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
433
432
434
433
435 self.dataOut.data_spc = jspectra
434 self.dataOut.data_spc = jspectra
436 self.dataOut.data_cspc = jcspectra
435 self.dataOut.data_cspc = jcspectra
437
436
438 return 1
437 return 1
439
438
440 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
439 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
441
440
442 jspectra = self.dataOut.data_spc
441 jspectra = self.dataOut.data_spc
443 jcspectra = self.dataOut.data_cspc
442 jcspectra = self.dataOut.data_cspc
444 jnoise = self.dataOut.getNoise()
443 jnoise = self.dataOut.getNoise()
445 num_incoh = self.dataOut.nIncohInt
444 num_incoh = self.dataOut.nIncohInt
446
445
447 num_channel = jspectra.shape[0]
446 num_channel = jspectra.shape[0]
448 num_prof = jspectra.shape[1]
447 num_prof = jspectra.shape[1]
449 num_hei = jspectra.shape[2]
448 num_hei = jspectra.shape[2]
450
449
451 #hei_interf
450 #hei_interf
452 if hei_interf is None:
451 if hei_interf is None:
453 count_hei = num_hei/2 #Como es entero no importa
452 count_hei = num_hei/2 #Como es entero no importa
454 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
453 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
455 hei_interf = numpy.asarray(hei_interf)[0]
454 hei_interf = numpy.asarray(hei_interf)[0]
456 #nhei_interf
455 #nhei_interf
457 if (nhei_interf == None):
456 if (nhei_interf == None):
458 nhei_interf = 5
457 nhei_interf = 5
459 if (nhei_interf < 1):
458 if (nhei_interf < 1):
460 nhei_interf = 1
459 nhei_interf = 1
461 if (nhei_interf > count_hei):
460 if (nhei_interf > count_hei):
462 nhei_interf = count_hei
461 nhei_interf = count_hei
463 if (offhei_interf == None):
462 if (offhei_interf == None):
464 offhei_interf = 0
463 offhei_interf = 0
465
464
466 ind_hei = range(num_hei)
465 ind_hei = range(num_hei)
467 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
466 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
468 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
467 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
469 mask_prof = numpy.asarray(range(num_prof))
468 mask_prof = numpy.asarray(range(num_prof))
470 num_mask_prof = mask_prof.size
469 num_mask_prof = mask_prof.size
471 comp_mask_prof = [0, num_prof/2]
470 comp_mask_prof = [0, num_prof/2]
472
471
473
472
474 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
473 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
475 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
474 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
476 jnoise = numpy.nan
475 jnoise = numpy.nan
477 noise_exist = jnoise[0] < numpy.Inf
476 noise_exist = jnoise[0] < numpy.Inf
478
477
479 #Subrutina de Remocion de la Interferencia
478 #Subrutina de Remocion de la Interferencia
480 for ich in range(num_channel):
479 for ich in range(num_channel):
481 #Se ordena los espectros segun su potencia (menor a mayor)
480 #Se ordena los espectros segun su potencia (menor a mayor)
482 power = jspectra[ich,mask_prof,:]
481 power = jspectra[ich,mask_prof,:]
483 power = power[:,hei_interf]
482 power = power[:,hei_interf]
484 power = power.sum(axis = 0)
483 power = power.sum(axis = 0)
485 psort = power.ravel().argsort()
484 psort = power.ravel().argsort()
486
485
487 #Se estima la interferencia promedio en los Espectros de Potencia empleando
486 #Se estima la interferencia promedio en los Espectros de Potencia empleando
488 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
487 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
489
488
490 if noise_exist:
489 if noise_exist:
491 # tmp_noise = jnoise[ich] / num_prof
490 # tmp_noise = jnoise[ich] / num_prof
492 tmp_noise = jnoise[ich]
491 tmp_noise = jnoise[ich]
493 junkspc_interf = junkspc_interf - tmp_noise
492 junkspc_interf = junkspc_interf - tmp_noise
494 #junkspc_interf[:,comp_mask_prof] = 0
493 #junkspc_interf[:,comp_mask_prof] = 0
495
494
496 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
495 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
497 jspc_interf = jspc_interf.transpose()
496 jspc_interf = jspc_interf.transpose()
498 #Calculando el espectro de interferencia promedio
497 #Calculando el espectro de interferencia promedio
499 noiseid = numpy.where(jspc_interf <= tmp_noise/ math.sqrt(num_incoh))
498 noiseid = numpy.where(jspc_interf <= tmp_noise/ numpy.sqrt(num_incoh))
500 noiseid = noiseid[0]
499 noiseid = noiseid[0]
501 cnoiseid = noiseid.size
500 cnoiseid = noiseid.size
502 interfid = numpy.where(jspc_interf > tmp_noise/ math.sqrt(num_incoh))
501 interfid = numpy.where(jspc_interf > tmp_noise/ numpy.sqrt(num_incoh))
503 interfid = interfid[0]
502 interfid = interfid[0]
504 cinterfid = interfid.size
503 cinterfid = interfid.size
505
504
506 if (cnoiseid > 0): jspc_interf[noiseid] = 0
505 if (cnoiseid > 0): jspc_interf[noiseid] = 0
507
506
508 #Expandiendo los perfiles a limpiar
507 #Expandiendo los perfiles a limpiar
509 if (cinterfid > 0):
508 if (cinterfid > 0):
510 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
509 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
511 new_interfid = numpy.asarray(new_interfid)
510 new_interfid = numpy.asarray(new_interfid)
512 new_interfid = {x for x in new_interfid}
511 new_interfid = {x for x in new_interfid}
513 new_interfid = numpy.array(list(new_interfid))
512 new_interfid = numpy.array(list(new_interfid))
514 new_cinterfid = new_interfid.size
513 new_cinterfid = new_interfid.size
515 else: new_cinterfid = 0
514 else: new_cinterfid = 0
516
515
517 for ip in range(new_cinterfid):
516 for ip in range(new_cinterfid):
518 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
517 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
519 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
518 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
520
519
521
520
522 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
521 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
523
522
524 #Removiendo la interferencia del punto de mayor interferencia
523 #Removiendo la interferencia del punto de mayor interferencia
525 ListAux = jspc_interf[mask_prof].tolist()
524 ListAux = jspc_interf[mask_prof].tolist()
526 maxid = ListAux.index(max(ListAux))
525 maxid = ListAux.index(max(ListAux))
527
526
528
527
529 if cinterfid > 0:
528 if cinterfid > 0:
530 for ip in range(cinterfid*(interf == 2) - 1):
529 for ip in range(cinterfid*(interf == 2) - 1):
531 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/math.sqrt(num_incoh))).nonzero()
530 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/numpy.sqrt(num_incoh))).nonzero()
532 cind = len(ind)
531 cind = len(ind)
533
532
534 if (cind > 0):
533 if (cind > 0):
535 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/math.sqrt(num_incoh))
534 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/numpy.sqrt(num_incoh))
536
535
537 ind = numpy.array([-2,-1,1,2])
536 ind = numpy.array([-2,-1,1,2])
538 xx = numpy.zeros([4,4])
537 xx = numpy.zeros([4,4])
539
538
540 for id1 in range(4):
539 for id1 in range(4):
541 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
540 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
542
541
543 xx_inv = numpy.linalg.inv(xx)
542 xx_inv = numpy.linalg.inv(xx)
544 xx = xx_inv[:,0]
543 xx = xx_inv[:,0]
545 ind = (ind + maxid + num_mask_prof)%num_mask_prof
544 ind = (ind + maxid + num_mask_prof)%num_mask_prof
546 yy = jspectra[ich,mask_prof[ind],:]
545 yy = jspectra[ich,mask_prof[ind],:]
547 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
546 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
548
547
549
548
550 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/math.sqrt(num_incoh))).nonzero()
549 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/numpy.sqrt(num_incoh))).nonzero()
551 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/math.sqrt(num_incoh))
550 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/numpy.sqrt(num_incoh))
552
551
553 #Remocion de Interferencia en el Cross Spectra
552 #Remocion de Interferencia en el Cross Spectra
554 if jcspectra is None: return jspectra, jcspectra
553 if jcspectra is None: return jspectra, jcspectra
555 num_pairs = jcspectra.size/(num_prof*num_hei)
554 num_pairs = jcspectra.size/(num_prof*num_hei)
556 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
555 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
557
556
558 for ip in range(num_pairs):
557 for ip in range(num_pairs):
559
558
560 #-------------------------------------------
559 #-------------------------------------------
561
560
562 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
561 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
563 cspower = cspower[:,hei_interf]
562 cspower = cspower[:,hei_interf]
564 cspower = cspower.sum(axis = 0)
563 cspower = cspower.sum(axis = 0)
565
564
566 cspsort = cspower.ravel().argsort()
565 cspsort = cspower.ravel().argsort()
567 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
566 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
568 junkcspc_interf = junkcspc_interf.transpose()
567 junkcspc_interf = junkcspc_interf.transpose()
569 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
568 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
570
569
571 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
570 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
572
571
573 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
572 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
574 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
573 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
575 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
574 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
576
575
577 for iprof in range(num_prof):
576 for iprof in range(num_prof):
578 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
577 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
579 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
578 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
580
579
581 #Removiendo la Interferencia
580 #Removiendo la Interferencia
582 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
581 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
583
582
584 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
583 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
585 maxid = ListAux.index(max(ListAux))
584 maxid = ListAux.index(max(ListAux))
586
585
587 ind = numpy.array([-2,-1,1,2])
586 ind = numpy.array([-2,-1,1,2])
588 xx = numpy.zeros([4,4])
587 xx = numpy.zeros([4,4])
589
588
590 for id1 in range(4):
589 for id1 in range(4):
591 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
590 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
592
591
593 xx_inv = numpy.linalg.inv(xx)
592 xx_inv = numpy.linalg.inv(xx)
594 xx = xx_inv[:,0]
593 xx = xx_inv[:,0]
595
594
596 ind = (ind + maxid + num_mask_prof)%num_mask_prof
595 ind = (ind + maxid + num_mask_prof)%num_mask_prof
597 yy = jcspectra[ip,mask_prof[ind],:]
596 yy = jcspectra[ip,mask_prof[ind],:]
598 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
597 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
599
598
600 #Guardar Resultados
599 #Guardar Resultados
601 self.dataOut.data_spc = jspectra
600 self.dataOut.data_spc = jspectra
602 self.dataOut.data_cspc = jcspectra
601 self.dataOut.data_cspc = jcspectra
603
602
604 return 1
603 return 1
605
604
606 def setRadarFrequency(self, frequency=None):
605 def setRadarFrequency(self, frequency=None):
607
606
608 if frequency != None:
607 if frequency != None:
609 self.dataOut.frequency = frequency
608 self.dataOut.frequency = frequency
610
609
611 return 1
610 return 1
612
611
613 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
612 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
614 #validacion de rango
613 #validacion de rango
615 if minHei == None:
614 if minHei == None:
616 minHei = self.dataOut.heightList[0]
615 minHei = self.dataOut.heightList[0]
617
616
618 if maxHei == None:
617 if maxHei == None:
619 maxHei = self.dataOut.heightList[-1]
618 maxHei = self.dataOut.heightList[-1]
620
619
621 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
620 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
622 print 'minHei: %.2f is out of the heights range'%(minHei)
621 print 'minHei: %.2f is out of the heights range'%(minHei)
623 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
622 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
624 minHei = self.dataOut.heightList[0]
623 minHei = self.dataOut.heightList[0]
625
624
626 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
625 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
627 print 'maxHei: %.2f is out of the heights range'%(maxHei)
626 print 'maxHei: %.2f is out of the heights range'%(maxHei)
628 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
627 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
629 maxHei = self.dataOut.heightList[-1]
628 maxHei = self.dataOut.heightList[-1]
630
629
631 # validacion de velocidades
630 # validacion de velocidades
632 velrange = self.dataOut.getVelRange(1)
631 velrange = self.dataOut.getVelRange(1)
633
632
634 if minVel == None:
633 if minVel == None:
635 minVel = velrange[0]
634 minVel = velrange[0]
636
635
637 if maxVel == None:
636 if maxVel == None:
638 maxVel = velrange[-1]
637 maxVel = velrange[-1]
639
638
640 if (minVel < velrange[0]) or (minVel > maxVel):
639 if (minVel < velrange[0]) or (minVel > maxVel):
641 print 'minVel: %.2f is out of the velocity range'%(minVel)
640 print 'minVel: %.2f is out of the velocity range'%(minVel)
642 print 'minVel is setting to %.2f'%(velrange[0])
641 print 'minVel is setting to %.2f'%(velrange[0])
643 minVel = velrange[0]
642 minVel = velrange[0]
644
643
645 if (maxVel > velrange[-1]) or (maxVel < minVel):
644 if (maxVel > velrange[-1]) or (maxVel < minVel):
646 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
645 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
647 print 'maxVel is setting to %.2f'%(velrange[-1])
646 print 'maxVel is setting to %.2f'%(velrange[-1])
648 maxVel = velrange[-1]
647 maxVel = velrange[-1]
649
648
650 # seleccion de indices para rango
649 # seleccion de indices para rango
651 minIndex = 0
650 minIndex = 0
652 maxIndex = 0
651 maxIndex = 0
653 heights = self.dataOut.heightList
652 heights = self.dataOut.heightList
654
653
655 inda = numpy.where(heights >= minHei)
654 inda = numpy.where(heights >= minHei)
656 indb = numpy.where(heights <= maxHei)
655 indb = numpy.where(heights <= maxHei)
657
656
658 try:
657 try:
659 minIndex = inda[0][0]
658 minIndex = inda[0][0]
660 except:
659 except:
661 minIndex = 0
660 minIndex = 0
662
661
663 try:
662 try:
664 maxIndex = indb[0][-1]
663 maxIndex = indb[0][-1]
665 except:
664 except:
666 maxIndex = len(heights)
665 maxIndex = len(heights)
667
666
668 if (minIndex < 0) or (minIndex > maxIndex):
667 if (minIndex < 0) or (minIndex > maxIndex):
669 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
668 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
670
669
671 if (maxIndex >= self.dataOut.nHeights):
670 if (maxIndex >= self.dataOut.nHeights):
672 maxIndex = self.dataOut.nHeights-1
671 maxIndex = self.dataOut.nHeights-1
673
672
674 # seleccion de indices para velocidades
673 # seleccion de indices para velocidades
675 indminvel = numpy.where(velrange >= minVel)
674 indminvel = numpy.where(velrange >= minVel)
676 indmaxvel = numpy.where(velrange <= maxVel)
675 indmaxvel = numpy.where(velrange <= maxVel)
677 try:
676 try:
678 minIndexVel = indminvel[0][0]
677 minIndexVel = indminvel[0][0]
679 except:
678 except:
680 minIndexVel = 0
679 minIndexVel = 0
681
680
682 try:
681 try:
683 maxIndexVel = indmaxvel[0][-1]
682 maxIndexVel = indmaxvel[0][-1]
684 except:
683 except:
685 maxIndexVel = len(velrange)
684 maxIndexVel = len(velrange)
686
685
687 #seleccion del espectro
686 #seleccion del espectro
688 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
687 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
689 #estimacion de ruido
688 #estimacion de ruido
690 noise = numpy.zeros(self.dataOut.nChannels)
689 noise = numpy.zeros(self.dataOut.nChannels)
691
690
692 for channel in range(self.dataOut.nChannels):
691 for channel in range(self.dataOut.nChannels):
693 daux = data_spc[channel,:,:]
692 daux = data_spc[channel,:,:]
694 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
693 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
695
694
696 self.dataOut.noise_estimation = noise.copy()
695 self.dataOut.noise_estimation = noise.copy()
697
696
698 return 1
697 return 1
699
698
700 class IncohInt(Operation):
699 class IncohInt(Operation):
701
700
702
701
703 __profIndex = 0
702 __profIndex = 0
704 __withOverapping = False
703 __withOverapping = False
705
704
706 __byTime = False
705 __byTime = False
707 __initime = None
706 __initime = None
708 __lastdatatime = None
707 __lastdatatime = None
709 __integrationtime = None
708 __integrationtime = None
710
709
711 __buffer_spc = None
710 __buffer_spc = None
712 __buffer_cspc = None
711 __buffer_cspc = None
713 __buffer_dc = None
712 __buffer_dc = None
714
713
715 __dataReady = False
714 __dataReady = False
716
715
717 __timeInterval = None
716 __timeInterval = None
718
717
719 n = None
718 n = None
720
719
721
720
722
721
723 def __init__(self):
722 def __init__(self):
724
723
725 Operation.__init__(self)
724 Operation.__init__(self)
726 # self.isConfig = False
725 # self.isConfig = False
727
726
728 def setup(self, n=None, timeInterval=None, overlapping=False):
727 def setup(self, n=None, timeInterval=None, overlapping=False):
729 """
728 """
730 Set the parameters of the integration class.
729 Set the parameters of the integration class.
731
730
732 Inputs:
731 Inputs:
733
732
734 n : Number of coherent integrations
733 n : Number of coherent integrations
735 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
734 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
736 overlapping :
735 overlapping :
737
736
738 """
737 """
739
738
740 self.__initime = None
739 self.__initime = None
741 self.__lastdatatime = 0
740 self.__lastdatatime = 0
742
741
743 self.__buffer_spc = 0
742 self.__buffer_spc = 0
744 self.__buffer_cspc = 0
743 self.__buffer_cspc = 0
745 self.__buffer_dc = 0
744 self.__buffer_dc = 0
746
745
747 self.__profIndex = 0
746 self.__profIndex = 0
748 self.__dataReady = False
747 self.__dataReady = False
749 self.__byTime = False
748 self.__byTime = False
750
749
751 if n is None and timeInterval is None:
750 if n is None and timeInterval is None:
752 raise ValueError, "n or timeInterval should be specified ..."
751 raise ValueError, "n or timeInterval should be specified ..."
753
752
754 if n is not None:
753 if n is not None:
755 self.n = int(n)
754 self.n = int(n)
756 else:
755 else:
757 self.__integrationtime = int(timeInterval) #if (type(timeInterval)!=integer) -> change this line
756 self.__integrationtime = int(timeInterval) #if (type(timeInterval)!=integer) -> change this line
758 self.n = None
757 self.n = None
759 self.__byTime = True
758 self.__byTime = True
760
759
761 def putData(self, data_spc, data_cspc, data_dc):
760 def putData(self, data_spc, data_cspc, data_dc):
762
761
763 """
762 """
764 Add a profile to the __buffer_spc and increase in one the __profileIndex
763 Add a profile to the __buffer_spc and increase in one the __profileIndex
765
764
766 """
765 """
767
766
768 self.__buffer_spc += data_spc
767 self.__buffer_spc += data_spc
769
768
770 if data_cspc is None:
769 if data_cspc is None:
771 self.__buffer_cspc = None
770 self.__buffer_cspc = None
772 else:
771 else:
773 self.__buffer_cspc += data_cspc
772 self.__buffer_cspc += data_cspc
774
773
775 if data_dc is None:
774 if data_dc is None:
776 self.__buffer_dc = None
775 self.__buffer_dc = None
777 else:
776 else:
778 self.__buffer_dc += data_dc
777 self.__buffer_dc += data_dc
779
778
780 self.__profIndex += 1
779 self.__profIndex += 1
781
780
782 return
781 return
783
782
784 def pushData(self):
783 def pushData(self):
785 """
784 """
786 Return the sum of the last profiles and the profiles used in the sum.
785 Return the sum of the last profiles and the profiles used in the sum.
787
786
788 Affected:
787 Affected:
789
788
790 self.__profileIndex
789 self.__profileIndex
791
790
792 """
791 """
793
792
794 data_spc = self.__buffer_spc
793 data_spc = self.__buffer_spc
795 data_cspc = self.__buffer_cspc
794 data_cspc = self.__buffer_cspc
796 data_dc = self.__buffer_dc
795 data_dc = self.__buffer_dc
797 n = self.__profIndex
796 n = self.__profIndex
798
797
799 self.__buffer_spc = 0
798 self.__buffer_spc = 0
800 self.__buffer_cspc = 0
799 self.__buffer_cspc = 0
801 self.__buffer_dc = 0
800 self.__buffer_dc = 0
802 self.__profIndex = 0
801 self.__profIndex = 0
803
802
804 return data_spc, data_cspc, data_dc, n
803 return data_spc, data_cspc, data_dc, n
805
804
806 def byProfiles(self, *args):
805 def byProfiles(self, *args):
807
806
808 self.__dataReady = False
807 self.__dataReady = False
809 avgdata_spc = None
808 avgdata_spc = None
810 avgdata_cspc = None
809 avgdata_cspc = None
811 avgdata_dc = None
810 avgdata_dc = None
812
811
813 self.putData(*args)
812 self.putData(*args)
814
813
815 if self.__profIndex == self.n:
814 if self.__profIndex == self.n:
816
815
817 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
816 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
818 self.n = n
817 self.n = n
819 self.__dataReady = True
818 self.__dataReady = True
820
819
821 return avgdata_spc, avgdata_cspc, avgdata_dc
820 return avgdata_spc, avgdata_cspc, avgdata_dc
822
821
823 def byTime(self, datatime, *args):
822 def byTime(self, datatime, *args):
824
823
825 self.__dataReady = False
824 self.__dataReady = False
826 avgdata_spc = None
825 avgdata_spc = None
827 avgdata_cspc = None
826 avgdata_cspc = None
828 avgdata_dc = None
827 avgdata_dc = None
829
828
830 self.putData(*args)
829 self.putData(*args)
831
830
832 if (datatime - self.__initime) >= self.__integrationtime:
831 if (datatime - self.__initime) >= self.__integrationtime:
833 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
832 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
834 self.n = n
833 self.n = n
835 self.__dataReady = True
834 self.__dataReady = True
836
835
837 return avgdata_spc, avgdata_cspc, avgdata_dc
836 return avgdata_spc, avgdata_cspc, avgdata_dc
838
837
839 def integrate(self, datatime, *args):
838 def integrate(self, datatime, *args):
840
839
841 if self.__profIndex == 0:
840 if self.__profIndex == 0:
842 self.__initime = datatime
841 self.__initime = datatime
843
842
844 if self.__byTime:
843 if self.__byTime:
845 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
844 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
846 else:
845 else:
847 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
846 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
848
847
849 if not self.__dataReady:
848 if not self.__dataReady:
850 return None, None, None, None
849 return None, None, None, None
851
850
852 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
851 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
853
852
854 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
853 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
855
854
856 if n==1:
855 if n==1:
857 return
856 return
858
857
859 dataOut.flagNoData = True
858 dataOut.flagNoData = True
860
859
861 if not self.isConfig:
860 if not self.isConfig:
862 self.setup(n, timeInterval, overlapping)
861 self.setup(n, timeInterval, overlapping)
863 self.isConfig = True
862 self.isConfig = True
864
863
865 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
864 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
866 dataOut.data_spc,
865 dataOut.data_spc,
867 dataOut.data_cspc,
866 dataOut.data_cspc,
868 dataOut.data_dc)
867 dataOut.data_dc)
869
868
870 if self.__dataReady:
869 if self.__dataReady:
871
870
872 dataOut.data_spc = avgdata_spc
871 dataOut.data_spc = avgdata_spc
873 dataOut.data_cspc = avgdata_cspc
872 dataOut.data_cspc = avgdata_cspc
874 dataOut.data_dc = avgdata_dc
873 dataOut.data_dc = avgdata_dc
875
874
876 dataOut.nIncohInt *= self.n
875 dataOut.nIncohInt *= self.n
877 dataOut.utctime = avgdatatime
876 dataOut.utctime = avgdatatime
878 dataOut.flagNoData = False
877 dataOut.flagNoData = False
General Comments 0
You need to be logged in to leave comments. Login now