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