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