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