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