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