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