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