##// END OF EJS Templates
schain-jc
RSS
Include ippFactor in getTimeInterval (Spectra)
Juan C. Espinoza -
r1160:c32e52f19bdc
parent child
Show More
Side by Side Unified
Context file:
r1160:c32e52f19bdc Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1251 +1,1251
1 '''
1 '''
2
2
3 $Author: murco $
3 $Author: murco $
4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 '''
5 '''
6
6
7 import copy
7 import copy
8 import numpy
8 import numpy
9 import datetime
9 import datetime
10
10
11 from jroheaderIO import SystemHeader, RadarControllerHeader
11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 from schainpy import cSchain
12 from schainpy import cSchain
13
13
14
14
15 def getNumpyDtype(dataTypeCode):
15 def getNumpyDtype(dataTypeCode):
16
16
17 if dataTypeCode == 0:
17 if dataTypeCode == 0:
18 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
18 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
19 elif dataTypeCode == 1:
19 elif dataTypeCode == 1:
20 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
20 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
21 elif dataTypeCode == 2:
21 elif dataTypeCode == 2:
22 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
22 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
23 elif dataTypeCode == 3:
23 elif dataTypeCode == 3:
24 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
24 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
25 elif dataTypeCode == 4:
25 elif dataTypeCode == 4:
26 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
26 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
27 elif dataTypeCode == 5:
27 elif dataTypeCode == 5:
28 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
28 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
29 else:
29 else:
30 raise ValueError, 'dataTypeCode was not defined'
30 raise ValueError, 'dataTypeCode was not defined'
31
31
32 return numpyDtype
32 return numpyDtype
33
33
34
34
35 def getDataTypeCode(numpyDtype):
35 def getDataTypeCode(numpyDtype):
36
36
37 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
37 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
38 datatype = 0
38 datatype = 0
39 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
39 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
40 datatype = 1
40 datatype = 1
41 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
41 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
42 datatype = 2
42 datatype = 2
43 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
43 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
44 datatype = 3
44 datatype = 3
45 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
45 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
46 datatype = 4
46 datatype = 4
47 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
47 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
48 datatype = 5
48 datatype = 5
49 else:
49 else:
50 datatype = None
50 datatype = None
51
51
52 return datatype
52 return datatype
53
53
54
54
55 def hildebrand_sekhon(data, navg):
55 def hildebrand_sekhon(data, navg):
56 """
56 """
57 This method is for the objective determination of the noise level in Doppler spectra. This
57 This method is for the objective determination of the noise level in Doppler spectra. This
58 implementation technique is based on the fact that the standard deviation of the spectral
58 implementation technique is based on the fact that the standard deviation of the spectral
59 densities is equal to the mean spectral density for white Gaussian noise
59 densities is equal to the mean spectral density for white Gaussian noise
60
60
61 Inputs:
61 Inputs:
62 Data : heights
62 Data : heights
63 navg : numbers of averages
63 navg : numbers of averages
64
64
65 Return:
65 Return:
66 -1 : any error
66 -1 : any error
67 anoise : noise's level
67 anoise : noise's level
68 """
68 """
69
69
70 sortdata = numpy.sort(data, axis=None)
70 sortdata = numpy.sort(data, axis=None)
71 # lenOfData = len(sortdata)
71 # lenOfData = len(sortdata)
72 # nums_min = lenOfData*0.2
72 # nums_min = lenOfData*0.2
73 #
73 #
74 # if nums_min <= 5:
74 # if nums_min <= 5:
75 # nums_min = 5
75 # nums_min = 5
76 #
76 #
77 # sump = 0.
77 # sump = 0.
78 #
78 #
79 # sumq = 0.
79 # sumq = 0.
80 #
80 #
81 # j = 0
81 # j = 0
82 #
82 #
83 # cont = 1
83 # cont = 1
84 #
84 #
85 # while((cont==1)and(j<lenOfData)):
85 # while((cont==1)and(j<lenOfData)):
86 #
86 #
87 # sump += sortdata[j]
87 # sump += sortdata[j]
88 #
88 #
89 # sumq += sortdata[j]**2
89 # sumq += sortdata[j]**2
90 #
90 #
91 # if j > nums_min:
91 # if j > nums_min:
92 # rtest = float(j)/(j-1) + 1.0/navg
92 # rtest = float(j)/(j-1) + 1.0/navg
93 # if ((sumq*j) > (rtest*sump**2)):
93 # if ((sumq*j) > (rtest*sump**2)):
94 # j = j - 1
94 # j = j - 1
95 # sump = sump - sortdata[j]
95 # sump = sump - sortdata[j]
96 # sumq = sumq - sortdata[j]**2
96 # sumq = sumq - sortdata[j]**2
97 # cont = 0
97 # cont = 0
98 #
98 #
99 # j += 1
99 # j += 1
100 #
100 #
101 # lnoise = sump /j
101 # lnoise = sump /j
102 #
102 #
103 # return lnoise
103 # return lnoise
104
104
105 return cSchain.hildebrand_sekhon(sortdata, navg)
105 return cSchain.hildebrand_sekhon(sortdata, navg)
106
106
107
107
108 class Beam:
108 class Beam:
109
109
110 def __init__(self):
110 def __init__(self):
111 self.codeList = []
111 self.codeList = []
112 self.azimuthList = []
112 self.azimuthList = []
113 self.zenithList = []
113 self.zenithList = []
114
114
115
115
116 class GenericData(object):
116 class GenericData(object):
117
117
118 flagNoData = True
118 flagNoData = True
119
119
120 def copy(self, inputObj=None):
120 def copy(self, inputObj=None):
121
121
122 if inputObj == None:
122 if inputObj == None:
123 return copy.deepcopy(self)
123 return copy.deepcopy(self)
124
124
125 for key in inputObj.__dict__.keys():
125 for key in inputObj.__dict__.keys():
126
126
127 attribute = inputObj.__dict__[key]
127 attribute = inputObj.__dict__[key]
128
128
129 # If this attribute is a tuple or list
129 # If this attribute is a tuple or list
130 if type(inputObj.__dict__[key]) in (tuple, list):
130 if type(inputObj.__dict__[key]) in (tuple, list):
131 self.__dict__[key] = attribute[:]
131 self.__dict__[key] = attribute[:]
132 continue
132 continue
133
133
134 # If this attribute is another object or instance
134 # If this attribute is another object or instance
135 if hasattr(attribute, '__dict__'):
135 if hasattr(attribute, '__dict__'):
136 self.__dict__[key] = attribute.copy()
136 self.__dict__[key] = attribute.copy()
137 continue
137 continue
138
138
139 self.__dict__[key] = inputObj.__dict__[key]
139 self.__dict__[key] = inputObj.__dict__[key]
140
140
141 def deepcopy(self):
141 def deepcopy(self):
142
142
143 return copy.deepcopy(self)
143 return copy.deepcopy(self)
144
144
145 def isEmpty(self):
145 def isEmpty(self):
146
146
147 return self.flagNoData
147 return self.flagNoData
148
148
149
149
150 class JROData(GenericData):
150 class JROData(GenericData):
151
151
152 # m_BasicHeader = BasicHeader()
152 # m_BasicHeader = BasicHeader()
153 # m_ProcessingHeader = ProcessingHeader()
153 # m_ProcessingHeader = ProcessingHeader()
154
154
155 systemHeaderObj = SystemHeader()
155 systemHeaderObj = SystemHeader()
156
156
157 radarControllerHeaderObj = RadarControllerHeader()
157 radarControllerHeaderObj = RadarControllerHeader()
158
158
159 # data = None
159 # data = None
160
160
161 type = None
161 type = None
162
162
163 datatype = None # dtype but in string
163 datatype = None # dtype but in string
164
164
165 # dtype = None
165 # dtype = None
166
166
167 # nChannels = None
167 # nChannels = None
168
168
169 # nHeights = None
169 # nHeights = None
170
170
171 nProfiles = None
171 nProfiles = None
172
172
173 heightList = None
173 heightList = None
174
174
175 channelList = None
175 channelList = None
176
176
177 flagDiscontinuousBlock = False
177 flagDiscontinuousBlock = False
178
178
179 useLocalTime = False
179 useLocalTime = False
180
180
181 utctime = None
181 utctime = None
182
182
183 timeZone = None
183 timeZone = None
184
184
185 dstFlag = None
185 dstFlag = None
186
186
187 errorCount = None
187 errorCount = None
188
188
189 blocksize = None
189 blocksize = None
190
190
191 # nCode = None
191 # nCode = None
192 #
192 #
193 # nBaud = None
193 # nBaud = None
194 #
194 #
195 # code = None
195 # code = None
196
196
197 flagDecodeData = False # asumo q la data no esta decodificada
197 flagDecodeData = False # asumo q la data no esta decodificada
198
198
199 flagDeflipData = False # asumo q la data no esta sin flip
199 flagDeflipData = False # asumo q la data no esta sin flip
200
200
201 flagShiftFFT = False
201 flagShiftFFT = False
202
202
203 # ippSeconds = None
203 # ippSeconds = None
204
204
205 # timeInterval = None
205 # timeInterval = None
206
206
207 nCohInt = None
207 nCohInt = None
208
208
209 # noise = None
209 # noise = None
210
210
211 windowOfFilter = 1
211 windowOfFilter = 1
212
212
213 # Speed of ligth
213 # Speed of ligth
214 C = 3e8
214 C = 3e8
215
215
216 frequency = 49.92e6
216 frequency = 49.92e6
217
217
218 realtime = False
218 realtime = False
219
219
220 beacon_heiIndexList = None
220 beacon_heiIndexList = None
221
221
222 last_block = None
222 last_block = None
223
223
224 blocknow = None
224 blocknow = None
225
225
226 azimuth = None
226 azimuth = None
227
227
228 zenith = None
228 zenith = None
229
229
230 beam = Beam()
230 beam = Beam()
231
231
232 profileIndex = None
232 profileIndex = None
233
233
234 def getNoise(self):
234 def getNoise(self):
235
235
236 raise NotImplementedError
236 raise NotImplementedError
237
237
238 def getNChannels(self):
238 def getNChannels(self):
239
239
240 return len(self.channelList)
240 return len(self.channelList)
241
241
242 def getChannelIndexList(self):
242 def getChannelIndexList(self):
243
243
244 return range(self.nChannels)
244 return range(self.nChannels)
245
245
246 def getNHeights(self):
246 def getNHeights(self):
247
247
248 return len(self.heightList)
248 return len(self.heightList)
249
249
250 def getHeiRange(self, extrapoints=0):
250 def getHeiRange(self, extrapoints=0):
251
251
252 heis = self.heightList
252 heis = self.heightList
253 # deltah = self.heightList[1] - self.heightList[0]
253 # deltah = self.heightList[1] - self.heightList[0]
254 #
254 #
255 # heis.append(self.heightList[-1])
255 # heis.append(self.heightList[-1])
256
256
257 return heis
257 return heis
258
258
259 def getDeltaH(self):
259 def getDeltaH(self):
260
260
261 delta = self.heightList[1] - self.heightList[0]
261 delta = self.heightList[1] - self.heightList[0]
262
262
263 return delta
263 return delta
264
264
265 def getltctime(self):
265 def getltctime(self):
266
266
267 if self.useLocalTime:
267 if self.useLocalTime:
268 return self.utctime - self.timeZone * 60
268 return self.utctime - self.timeZone * 60
269
269
270 return self.utctime
270 return self.utctime
271
271
272 def getDatatime(self):
272 def getDatatime(self):
273
273
274 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
274 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
275 return datatimeValue
275 return datatimeValue
276
276
277 def getTimeRange(self):
277 def getTimeRange(self):
278
278
279 datatime = []
279 datatime = []
280
280
281 datatime.append(self.ltctime)
281 datatime.append(self.ltctime)
282 datatime.append(self.ltctime + self.timeInterval + 1)
282 datatime.append(self.ltctime + self.timeInterval + 1)
283
283
284 datatime = numpy.array(datatime)
284 datatime = numpy.array(datatime)
285
285
286 return datatime
286 return datatime
287
287
288 def getFmaxTimeResponse(self):
288 def getFmaxTimeResponse(self):
289
289
290 period = (10**-6) * self.getDeltaH() / (0.15)
290 period = (10**-6) * self.getDeltaH() / (0.15)
291
291
292 PRF = 1. / (period * self.nCohInt)
292 PRF = 1. / (period * self.nCohInt)
293
293
294 fmax = PRF
294 fmax = PRF
295
295
296 return fmax
296 return fmax
297
297
298 def getFmax(self):
298 def getFmax(self):
299 PRF = 1. / (self.ippSeconds * self.nCohInt)
299 PRF = 1. / (self.ippSeconds * self.nCohInt)
300
300
301 fmax = PRF
301 fmax = PRF
302 return fmax
302 return fmax
303
303
304 def getVmax(self):
304 def getVmax(self):
305
305
306 _lambda = self.C / self.frequency
306 _lambda = self.C / self.frequency
307
307
308 vmax = self.getFmax() * _lambda / 2
308 vmax = self.getFmax() * _lambda / 2
309
309
310 return vmax
310 return vmax
311
311
312 def get_ippSeconds(self):
312 def get_ippSeconds(self):
313 '''
313 '''
314 '''
314 '''
315 return self.radarControllerHeaderObj.ippSeconds
315 return self.radarControllerHeaderObj.ippSeconds
316
316
317 def set_ippSeconds(self, ippSeconds):
317 def set_ippSeconds(self, ippSeconds):
318 '''
318 '''
319 '''
319 '''
320
320
321 self.radarControllerHeaderObj.ippSeconds = ippSeconds
321 self.radarControllerHeaderObj.ippSeconds = ippSeconds
322
322
323 return
323 return
324
324
325 def get_dtype(self):
325 def get_dtype(self):
326 '''
326 '''
327 '''
327 '''
328 return getNumpyDtype(self.datatype)
328 return getNumpyDtype(self.datatype)
329
329
330 def set_dtype(self, numpyDtype):
330 def set_dtype(self, numpyDtype):
331 '''
331 '''
332 '''
332 '''
333
333
334 self.datatype = getDataTypeCode(numpyDtype)
334 self.datatype = getDataTypeCode(numpyDtype)
335
335
336 def get_code(self):
336 def get_code(self):
337 '''
337 '''
338 '''
338 '''
339 return self.radarControllerHeaderObj.code
339 return self.radarControllerHeaderObj.code
340
340
341 def set_code(self, code):
341 def set_code(self, code):
342 '''
342 '''
343 '''
343 '''
344 self.radarControllerHeaderObj.code = code
344 self.radarControllerHeaderObj.code = code
345
345
346 return
346 return
347
347
348 def get_ncode(self):
348 def get_ncode(self):
349 '''
349 '''
350 '''
350 '''
351 return self.radarControllerHeaderObj.nCode
351 return self.radarControllerHeaderObj.nCode
352
352
353 def set_ncode(self, nCode):
353 def set_ncode(self, nCode):
354 '''
354 '''
355 '''
355 '''
356 self.radarControllerHeaderObj.nCode = nCode
356 self.radarControllerHeaderObj.nCode = nCode
357
357
358 return
358 return
359
359
360 def get_nbaud(self):
360 def get_nbaud(self):
361 '''
361 '''
362 '''
362 '''
363 return self.radarControllerHeaderObj.nBaud
363 return self.radarControllerHeaderObj.nBaud
364
364
365 def set_nbaud(self, nBaud):
365 def set_nbaud(self, nBaud):
366 '''
366 '''
367 '''
367 '''
368 self.radarControllerHeaderObj.nBaud = nBaud
368 self.radarControllerHeaderObj.nBaud = nBaud
369
369
370 return
370 return
371
371
372 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
372 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
373 channelIndexList = property(
373 channelIndexList = property(
374 getChannelIndexList, "I'm the 'channelIndexList' property.")
374 getChannelIndexList, "I'm the 'channelIndexList' property.")
375 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
375 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
376 #noise = property(getNoise, "I'm the 'nHeights' property.")
376 #noise = property(getNoise, "I'm the 'nHeights' property.")
377 datatime = property(getDatatime, "I'm the 'datatime' property")
377 datatime = property(getDatatime, "I'm the 'datatime' property")
378 ltctime = property(getltctime, "I'm the 'ltctime' property")
378 ltctime = property(getltctime, "I'm the 'ltctime' property")
379 ippSeconds = property(get_ippSeconds, set_ippSeconds)
379 ippSeconds = property(get_ippSeconds, set_ippSeconds)
380 dtype = property(get_dtype, set_dtype)
380 dtype = property(get_dtype, set_dtype)
381 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
381 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
382 code = property(get_code, set_code)
382 code = property(get_code, set_code)
383 nCode = property(get_ncode, set_ncode)
383 nCode = property(get_ncode, set_ncode)
384 nBaud = property(get_nbaud, set_nbaud)
384 nBaud = property(get_nbaud, set_nbaud)
385
385
386
386
387 class Voltage(JROData):
387 class Voltage(JROData):
388
388
389 # data es un numpy array de 2 dmensiones (canales, alturas)
389 # data es un numpy array de 2 dmensiones (canales, alturas)
390 data = None
390 data = None
391
391
392 def __init__(self):
392 def __init__(self):
393 '''
393 '''
394 Constructor
394 Constructor
395 '''
395 '''
396
396
397 self.useLocalTime = True
397 self.useLocalTime = True
398
398
399 self.radarControllerHeaderObj = RadarControllerHeader()
399 self.radarControllerHeaderObj = RadarControllerHeader()
400
400
401 self.systemHeaderObj = SystemHeader()
401 self.systemHeaderObj = SystemHeader()
402
402
403 self.type = "Voltage"
403 self.type = "Voltage"
404
404
405 self.data = None
405 self.data = None
406
406
407 # self.dtype = None
407 # self.dtype = None
408
408
409 # self.nChannels = 0
409 # self.nChannels = 0
410
410
411 # self.nHeights = 0
411 # self.nHeights = 0
412
412
413 self.nProfiles = None
413 self.nProfiles = None
414
414
415 self.heightList = None
415 self.heightList = None
416
416
417 self.channelList = None
417 self.channelList = None
418
418
419 # self.channelIndexList = None
419 # self.channelIndexList = None
420
420
421 self.flagNoData = True
421 self.flagNoData = True
422
422
423 self.flagDiscontinuousBlock = False
423 self.flagDiscontinuousBlock = False
424
424
425 self.utctime = None
425 self.utctime = None
426
426
427 self.timeZone = None
427 self.timeZone = None
428
428
429 self.dstFlag = None
429 self.dstFlag = None
430
430
431 self.errorCount = None
431 self.errorCount = None
432
432
433 self.nCohInt = None
433 self.nCohInt = None
434
434
435 self.blocksize = None
435 self.blocksize = None
436
436
437 self.flagDecodeData = False # asumo q la data no esta decodificada
437 self.flagDecodeData = False # asumo q la data no esta decodificada
438
438
439 self.flagDeflipData = False # asumo q la data no esta sin flip
439 self.flagDeflipData = False # asumo q la data no esta sin flip
440
440
441 self.flagShiftFFT = False
441 self.flagShiftFFT = False
442
442
443 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
443 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
444
444
445 self.profileIndex = 0
445 self.profileIndex = 0
446
446
447 def getNoisebyHildebrand(self, channel=None):
447 def getNoisebyHildebrand(self, channel=None):
448 """
448 """
449 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
449 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
450
450
451 Return:
451 Return:
452 noiselevel
452 noiselevel
453 """
453 """
454
454
455 if channel != None:
455 if channel != None:
456 data = self.data[channel]
456 data = self.data[channel]
457 nChannels = 1
457 nChannels = 1
458 else:
458 else:
459 data = self.data
459 data = self.data
460 nChannels = self.nChannels
460 nChannels = self.nChannels
461
461
462 noise = numpy.zeros(nChannels)
462 noise = numpy.zeros(nChannels)
463 power = data * numpy.conjugate(data)
463 power = data * numpy.conjugate(data)
464
464
465 for thisChannel in range(nChannels):
465 for thisChannel in range(nChannels):
466 if nChannels == 1:
466 if nChannels == 1:
467 daux = power[:].real
467 daux = power[:].real
468 else:
468 else:
469 daux = power[thisChannel, :].real
469 daux = power[thisChannel, :].real
470 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
470 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
471
471
472 return noise
472 return noise
473
473
474 def getNoise(self, type=1, channel=None):
474 def getNoise(self, type=1, channel=None):
475
475
476 if type == 1:
476 if type == 1:
477 noise = self.getNoisebyHildebrand(channel)
477 noise = self.getNoisebyHildebrand(channel)
478
478
479 return noise
479 return noise
480
480
481 def getPower(self, channel=None):
481 def getPower(self, channel=None):
482
482
483 if channel != None:
483 if channel != None:
484 data = self.data[channel]
484 data = self.data[channel]
485 else:
485 else:
486 data = self.data
486 data = self.data
487
487
488 power = data * numpy.conjugate(data)
488 power = data * numpy.conjugate(data)
489 powerdB = 10 * numpy.log10(power.real)
489 powerdB = 10 * numpy.log10(power.real)
490 powerdB = numpy.squeeze(powerdB)
490 powerdB = numpy.squeeze(powerdB)
491
491
492 return powerdB
492 return powerdB
493
493
494 def getTimeInterval(self):
494 def getTimeInterval(self):
495
495
496 timeInterval = self.ippSeconds * self.nCohInt
496 timeInterval = self.ippSeconds * self.nCohInt
497
497
498 return timeInterval
498 return timeInterval
499
499
500 noise = property(getNoise, "I'm the 'nHeights' property.")
500 noise = property(getNoise, "I'm the 'nHeights' property.")
501 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
501 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
502
502
503
503
504 class Spectra(JROData):
504 class Spectra(JROData):
505
505
506 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
506 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
507 data_spc = None
507 data_spc = None
508
508
509 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
509 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
510 data_cspc = None
510 data_cspc = None
511
511
512 # data dc es un numpy array de 2 dmensiones (canales, alturas)
512 # data dc es un numpy array de 2 dmensiones (canales, alturas)
513 data_dc = None
513 data_dc = None
514
514
515 # data power
515 # data power
516 data_pwr = None
516 data_pwr = None
517
517
518 nFFTPoints = None
518 nFFTPoints = None
519
519
520 # nPairs = None
520 # nPairs = None
521
521
522 pairsList = None
522 pairsList = None
523
523
524 nIncohInt = None
524 nIncohInt = None
525
525
526 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
526 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
527
527
528 nCohInt = None # se requiere para determinar el valor de timeInterval
528 nCohInt = None # se requiere para determinar el valor de timeInterval
529
529
530 ippFactor = None
530 ippFactor = None
531
531
532 profileIndex = 0
532 profileIndex = 0
533
533
534 plotting = "spectra"
534 plotting = "spectra"
535
535
536 def __init__(self):
536 def __init__(self):
537 '''
537 '''
538 Constructor
538 Constructor
539 '''
539 '''
540
540
541 self.useLocalTime = True
541 self.useLocalTime = True
542
542
543 self.radarControllerHeaderObj = RadarControllerHeader()
543 self.radarControllerHeaderObj = RadarControllerHeader()
544
544
545 self.systemHeaderObj = SystemHeader()
545 self.systemHeaderObj = SystemHeader()
546
546
547 self.type = "Spectra"
547 self.type = "Spectra"
548
548
549 # self.data = None
549 # self.data = None
550
550
551 # self.dtype = None
551 # self.dtype = None
552
552
553 # self.nChannels = 0
553 # self.nChannels = 0
554
554
555 # self.nHeights = 0
555 # self.nHeights = 0
556
556
557 self.nProfiles = None
557 self.nProfiles = None
558
558
559 self.heightList = None
559 self.heightList = None
560
560
561 self.channelList = None
561 self.channelList = None
562
562
563 # self.channelIndexList = None
563 # self.channelIndexList = None
564
564
565 self.pairsList = None
565 self.pairsList = None
566
566
567 self.flagNoData = True
567 self.flagNoData = True
568
568
569 self.flagDiscontinuousBlock = False
569 self.flagDiscontinuousBlock = False
570
570
571 self.utctime = None
571 self.utctime = None
572
572
573 self.nCohInt = None
573 self.nCohInt = None
574
574
575 self.nIncohInt = None
575 self.nIncohInt = None
576
576
577 self.blocksize = None
577 self.blocksize = None
578
578
579 self.nFFTPoints = None
579 self.nFFTPoints = None
580
580
581 self.wavelength = None
581 self.wavelength = None
582
582
583 self.flagDecodeData = False # asumo q la data no esta decodificada
583 self.flagDecodeData = False # asumo q la data no esta decodificada
584
584
585 self.flagDeflipData = False # asumo q la data no esta sin flip
585 self.flagDeflipData = False # asumo q la data no esta sin flip
586
586
587 self.flagShiftFFT = False
587 self.flagShiftFFT = False
588
588
589 self.ippFactor = 1
589 self.ippFactor = 1
590
590
591 #self.noise = None
591 #self.noise = None
592
592
593 self.beacon_heiIndexList = []
593 self.beacon_heiIndexList = []
594
594
595 self.noise_estimation = None
595 self.noise_estimation = None
596
596
597 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
597 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
598 """
598 """
599 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
599 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
600
600
601 Return:
601 Return:
602 noiselevel
602 noiselevel
603 """
603 """
604
604
605 noise = numpy.zeros(self.nChannels)
605 noise = numpy.zeros(self.nChannels)
606
606
607 for channel in range(self.nChannels):
607 for channel in range(self.nChannels):
608 daux = self.data_spc[channel,
608 daux = self.data_spc[channel,
609 xmin_index:xmax_index, ymin_index:ymax_index]
609 xmin_index:xmax_index, ymin_index:ymax_index]
610 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
610 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
611
611
612 return noise
612 return noise
613
613
614 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
614 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
615
615
616 if self.noise_estimation is not None:
616 if self.noise_estimation is not None:
617 # this was estimated by getNoise Operation defined in jroproc_spectra.py
617 # this was estimated by getNoise Operation defined in jroproc_spectra.py
618 return self.noise_estimation
618 return self.noise_estimation
619 else:
619 else:
620 noise = self.getNoisebyHildebrand(
620 noise = self.getNoisebyHildebrand(
621 xmin_index, xmax_index, ymin_index, ymax_index)
621 xmin_index, xmax_index, ymin_index, ymax_index)
622 return noise
622 return noise
623
623
624 def getFreqRangeTimeResponse(self, extrapoints=0):
624 def getFreqRangeTimeResponse(self, extrapoints=0):
625
625
626 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
626 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
627 freqrange = deltafreq * \
627 freqrange = deltafreq * \
628 (numpy.arange(self.nFFTPoints + extrapoints) -
628 (numpy.arange(self.nFFTPoints + extrapoints) -
629 self.nFFTPoints / 2.) - deltafreq / 2
629 self.nFFTPoints / 2.) - deltafreq / 2
630
630
631 return freqrange
631 return freqrange
632
632
633 def getAcfRange(self, extrapoints=0):
633 def getAcfRange(self, extrapoints=0):
634
634
635 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
635 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
636 freqrange = deltafreq * \
636 freqrange = deltafreq * \
637 (numpy.arange(self.nFFTPoints + extrapoints) -
637 (numpy.arange(self.nFFTPoints + extrapoints) -
638 self.nFFTPoints / 2.) - deltafreq / 2
638 self.nFFTPoints / 2.) - deltafreq / 2
639
639
640 return freqrange
640 return freqrange
641
641
642 def getFreqRange(self, extrapoints=0):
642 def getFreqRange(self, extrapoints=0):
643
643
644 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
644 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
645 freqrange = deltafreq * \
645 freqrange = deltafreq * \
646 (numpy.arange(self.nFFTPoints + extrapoints) -
646 (numpy.arange(self.nFFTPoints + extrapoints) -
647 self.nFFTPoints / 2.) - deltafreq / 2
647 self.nFFTPoints / 2.) - deltafreq / 2
648
648
649 return freqrange
649 return freqrange
650
650
651 def getVelRange(self, extrapoints=0):
651 def getVelRange(self, extrapoints=0):
652
652
653 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
653 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
654 velrange = deltav * (numpy.arange(self.nFFTPoints +
654 velrange = deltav * (numpy.arange(self.nFFTPoints +
655 extrapoints) - self.nFFTPoints / 2.) # - deltav/2
655 extrapoints) - self.nFFTPoints / 2.) # - deltav/2
656
656
657 return velrange
657 return velrange
658
658
659 def getNPairs(self):
659 def getNPairs(self):
660
660
661 return len(self.pairsList)
661 return len(self.pairsList)
662
662
663 def getPairsIndexList(self):
663 def getPairsIndexList(self):
664
664
665 return range(self.nPairs)
665 return range(self.nPairs)
666
666
667 def getNormFactor(self):
667 def getNormFactor(self):
668
668
669 pwcode = 1
669 pwcode = 1
670
670
671 if self.flagDecodeData:
671 if self.flagDecodeData:
672 pwcode = numpy.sum(self.code[0]**2)
672 pwcode = numpy.sum(self.code[0]**2)
673 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
673 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
674 normFactor = self.nProfiles * self.nIncohInt * \
674 normFactor = self.nProfiles * self.nIncohInt * \
675 self.nCohInt * pwcode * self.windowOfFilter
675 self.nCohInt * pwcode * self.windowOfFilter
676
676
677 return normFactor
677 return normFactor
678
678
679 def getFlagCspc(self):
679 def getFlagCspc(self):
680
680
681 if self.data_cspc is None:
681 if self.data_cspc is None:
682 return True
682 return True
683
683
684 return False
684 return False
685
685
686 def getFlagDc(self):
686 def getFlagDc(self):
687
687
688 if self.data_dc is None:
688 if self.data_dc is None:
689 return True
689 return True
690
690
691 return False
691 return False
692
692
693 def getTimeInterval(self):
693 def getTimeInterval(self):
694
694
695 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
695 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
696
696
697 return timeInterval
697 return timeInterval
698
698
699 def getPower(self):
699 def getPower(self):
700
700
701 factor = self.normFactor
701 factor = self.normFactor
702 z = self.data_spc / factor
702 z = self.data_spc / factor
703 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
703 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
704 avg = numpy.average(z, axis=1)
704 avg = numpy.average(z, axis=1)
705
705
706 return 10 * numpy.log10(avg)
706 return 10 * numpy.log10(avg)
707
707
708 def getCoherence(self, pairsList=None, phase=False):
708 def getCoherence(self, pairsList=None, phase=False):
709
709
710 z = []
710 z = []
711 if pairsList is None:
711 if pairsList is None:
712 pairsIndexList = self.pairsIndexList
712 pairsIndexList = self.pairsIndexList
713 else:
713 else:
714 pairsIndexList = []
714 pairsIndexList = []
715 for pair in pairsList:
715 for pair in pairsList:
716 if pair not in self.pairsList:
716 if pair not in self.pairsList:
717 raise ValueError, "Pair %s is not in dataOut.pairsList" % (
717 raise ValueError, "Pair %s is not in dataOut.pairsList" % (
718 pair)
718 pair)
719 pairsIndexList.append(self.pairsList.index(pair))
719 pairsIndexList.append(self.pairsList.index(pair))
720 for i in range(len(pairsIndexList)):
720 for i in range(len(pairsIndexList)):
721 pair = self.pairsList[pairsIndexList[i]]
721 pair = self.pairsList[pairsIndexList[i]]
722 ccf = numpy.average(
722 ccf = numpy.average(
723 self.data_cspc[pairsIndexList[i], :, :], axis=0)
723 self.data_cspc[pairsIndexList[i], :, :], axis=0)
724 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
724 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
725 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
725 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
726 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
726 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
727 if phase:
727 if phase:
728 data = numpy.arctan2(avgcoherenceComplex.imag,
728 data = numpy.arctan2(avgcoherenceComplex.imag,
729 avgcoherenceComplex.real) * 180 / numpy.pi
729 avgcoherenceComplex.real) * 180 / numpy.pi
730 else:
730 else:
731 data = numpy.abs(avgcoherenceComplex)
731 data = numpy.abs(avgcoherenceComplex)
732
732
733 z.append(data)
733 z.append(data)
734
734
735 return numpy.array(z)
735 return numpy.array(z)
736
736
737 def setValue(self, value):
737 def setValue(self, value):
738
738
739 print "This property should not be initialized"
739 print "This property should not be initialized"
740
740
741 return
741 return
742
742
743 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
743 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
744 pairsIndexList = property(
744 pairsIndexList = property(
745 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
745 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
746 normFactor = property(getNormFactor, setValue,
746 normFactor = property(getNormFactor, setValue,
747 "I'm the 'getNormFactor' property.")
747 "I'm the 'getNormFactor' property.")
748 flag_cspc = property(getFlagCspc, setValue)
748 flag_cspc = property(getFlagCspc, setValue)
749 flag_dc = property(getFlagDc, setValue)
749 flag_dc = property(getFlagDc, setValue)
750 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
750 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
751 timeInterval = property(getTimeInterval, setValue,
751 timeInterval = property(getTimeInterval, setValue,
752 "I'm the 'timeInterval' property")
752 "I'm the 'timeInterval' property")
753
753
754
754
755 class SpectraHeis(Spectra):
755 class SpectraHeis(Spectra):
756
756
757 data_spc = None
757 data_spc = None
758
758
759 data_cspc = None
759 data_cspc = None
760
760
761 data_dc = None
761 data_dc = None
762
762
763 nFFTPoints = None
763 nFFTPoints = None
764
764
765 # nPairs = None
765 # nPairs = None
766
766
767 pairsList = None
767 pairsList = None
768
768
769 nCohInt = None
769 nCohInt = None
770
770
771 nIncohInt = None
771 nIncohInt = None
772
772
773 def __init__(self):
773 def __init__(self):
774
774
775 self.radarControllerHeaderObj = RadarControllerHeader()
775 self.radarControllerHeaderObj = RadarControllerHeader()
776
776
777 self.systemHeaderObj = SystemHeader()
777 self.systemHeaderObj = SystemHeader()
778
778
779 self.type = "SpectraHeis"
779 self.type = "SpectraHeis"
780
780
781 # self.dtype = None
781 # self.dtype = None
782
782
783 # self.nChannels = 0
783 # self.nChannels = 0
784
784
785 # self.nHeights = 0
785 # self.nHeights = 0
786
786
787 self.nProfiles = None
787 self.nProfiles = None
788
788
789 self.heightList = None
789 self.heightList = None
790
790
791 self.channelList = None
791 self.channelList = None
792
792
793 # self.channelIndexList = None
793 # self.channelIndexList = None
794
794
795 self.flagNoData = True
795 self.flagNoData = True
796
796
797 self.flagDiscontinuousBlock = False
797 self.flagDiscontinuousBlock = False
798
798
799 # self.nPairs = 0
799 # self.nPairs = 0
800
800
801 self.utctime = None
801 self.utctime = None
802
802
803 self.blocksize = None
803 self.blocksize = None
804
804
805 self.profileIndex = 0
805 self.profileIndex = 0
806
806
807 self.nCohInt = 1
807 self.nCohInt = 1
808
808
809 self.nIncohInt = 1
809 self.nIncohInt = 1
810
810
811 def getNormFactor(self):
811 def getNormFactor(self):
812 pwcode = 1
812 pwcode = 1
813 if self.flagDecodeData:
813 if self.flagDecodeData:
814 pwcode = numpy.sum(self.code[0]**2)
814 pwcode = numpy.sum(self.code[0]**2)
815
815
816 normFactor = self.nIncohInt * self.nCohInt * pwcode
816 normFactor = self.nIncohInt * self.nCohInt * pwcode
817
817
818 return normFactor
818 return normFactor
819
819
820 def getTimeInterval(self):
820 def getTimeInterval(self):
821
821
822 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
822 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
823
823
824 return timeInterval
824 return timeInterval
825
825
826 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
826 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
827 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
827 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
828
828
829
829
830 class Fits(JROData):
830 class Fits(JROData):
831
831
832 heightList = None
832 heightList = None
833
833
834 channelList = None
834 channelList = None
835
835
836 flagNoData = True
836 flagNoData = True
837
837
838 flagDiscontinuousBlock = False
838 flagDiscontinuousBlock = False
839
839
840 useLocalTime = False
840 useLocalTime = False
841
841
842 utctime = None
842 utctime = None
843
843
844 timeZone = None
844 timeZone = None
845
845
846 # ippSeconds = None
846 # ippSeconds = None
847
847
848 # timeInterval = None
848 # timeInterval = None
849
849
850 nCohInt = None
850 nCohInt = None
851
851
852 nIncohInt = None
852 nIncohInt = None
853
853
854 noise = None
854 noise = None
855
855
856 windowOfFilter = 1
856 windowOfFilter = 1
857
857
858 # Speed of ligth
858 # Speed of ligth
859 C = 3e8
859 C = 3e8
860
860
861 frequency = 49.92e6
861 frequency = 49.92e6
862
862
863 realtime = False
863 realtime = False
864
864
865 def __init__(self):
865 def __init__(self):
866
866
867 self.type = "Fits"
867 self.type = "Fits"
868
868
869 self.nProfiles = None
869 self.nProfiles = None
870
870
871 self.heightList = None
871 self.heightList = None
872
872
873 self.channelList = None
873 self.channelList = None
874
874
875 # self.channelIndexList = None
875 # self.channelIndexList = None
876
876
877 self.flagNoData = True
877 self.flagNoData = True
878
878
879 self.utctime = None
879 self.utctime = None
880
880
881 self.nCohInt = 1
881 self.nCohInt = 1
882
882
883 self.nIncohInt = 1
883 self.nIncohInt = 1
884
884
885 self.useLocalTime = True
885 self.useLocalTime = True
886
886
887 self.profileIndex = 0
887 self.profileIndex = 0
888
888
889 # self.utctime = None
889 # self.utctime = None
890 # self.timeZone = None
890 # self.timeZone = None
891 # self.ltctime = None
891 # self.ltctime = None
892 # self.timeInterval = None
892 # self.timeInterval = None
893 # self.header = None
893 # self.header = None
894 # self.data_header = None
894 # self.data_header = None
895 # self.data = None
895 # self.data = None
896 # self.datatime = None
896 # self.datatime = None
897 # self.flagNoData = False
897 # self.flagNoData = False
898 # self.expName = ''
898 # self.expName = ''
899 # self.nChannels = None
899 # self.nChannels = None
900 # self.nSamples = None
900 # self.nSamples = None
901 # self.dataBlocksPerFile = None
901 # self.dataBlocksPerFile = None
902 # self.comments = ''
902 # self.comments = ''
903 #
903 #
904
904
905 def getltctime(self):
905 def getltctime(self):
906
906
907 if self.useLocalTime:
907 if self.useLocalTime:
908 return self.utctime - self.timeZone * 60
908 return self.utctime - self.timeZone * 60
909
909
910 return self.utctime
910 return self.utctime
911
911
912 def getDatatime(self):
912 def getDatatime(self):
913
913
914 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
914 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
915 return datatime
915 return datatime
916
916
917 def getTimeRange(self):
917 def getTimeRange(self):
918
918
919 datatime = []
919 datatime = []
920
920
921 datatime.append(self.ltctime)
921 datatime.append(self.ltctime)
922 datatime.append(self.ltctime + self.timeInterval)
922 datatime.append(self.ltctime + self.timeInterval)
923
923
924 datatime = numpy.array(datatime)
924 datatime = numpy.array(datatime)
925
925
926 return datatime
926 return datatime
927
927
928 def getHeiRange(self):
928 def getHeiRange(self):
929
929
930 heis = self.heightList
930 heis = self.heightList
931
931
932 return heis
932 return heis
933
933
934 def getNHeights(self):
934 def getNHeights(self):
935
935
936 return len(self.heightList)
936 return len(self.heightList)
937
937
938 def getNChannels(self):
938 def getNChannels(self):
939
939
940 return len(self.channelList)
940 return len(self.channelList)
941
941
942 def getChannelIndexList(self):
942 def getChannelIndexList(self):
943
943
944 return range(self.nChannels)
944 return range(self.nChannels)
945
945
946 def getNoise(self, type=1):
946 def getNoise(self, type=1):
947
947
948 #noise = numpy.zeros(self.nChannels)
948 #noise = numpy.zeros(self.nChannels)
949
949
950 if type == 1:
950 if type == 1:
951 noise = self.getNoisebyHildebrand()
951 noise = self.getNoisebyHildebrand()
952
952
953 if type == 2:
953 if type == 2:
954 noise = self.getNoisebySort()
954 noise = self.getNoisebySort()
955
955
956 if type == 3:
956 if type == 3:
957 noise = self.getNoisebyWindow()
957 noise = self.getNoisebyWindow()
958
958
959 return noise
959 return noise
960
960
961 def getTimeInterval(self):
961 def getTimeInterval(self):
962
962
963 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
963 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
964
964
965 return timeInterval
965 return timeInterval
966
966
967 datatime = property(getDatatime, "I'm the 'datatime' property")
967 datatime = property(getDatatime, "I'm the 'datatime' property")
968 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
968 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
969 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
969 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
970 channelIndexList = property(
970 channelIndexList = property(
971 getChannelIndexList, "I'm the 'channelIndexList' property.")
971 getChannelIndexList, "I'm the 'channelIndexList' property.")
972 noise = property(getNoise, "I'm the 'nHeights' property.")
972 noise = property(getNoise, "I'm the 'nHeights' property.")
973
973
974 ltctime = property(getltctime, "I'm the 'ltctime' property")
974 ltctime = property(getltctime, "I'm the 'ltctime' property")
975 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
975 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
976
976
977
977
978 class Correlation(JROData):
978 class Correlation(JROData):
979
979
980 noise = None
980 noise = None
981
981
982 SNR = None
982 SNR = None
983
983
984 #--------------------------------------------------
984 #--------------------------------------------------
985
985
986 mode = None
986 mode = None
987
987
988 split = False
988 split = False
989
989
990 data_cf = None
990 data_cf = None
991
991
992 lags = None
992 lags = None
993
993
994 lagRange = None
994 lagRange = None
995
995
996 pairsList = None
996 pairsList = None
997
997
998 normFactor = None
998 normFactor = None
999
999
1000 #--------------------------------------------------
1000 #--------------------------------------------------
1001
1001
1002 # calculateVelocity = None
1002 # calculateVelocity = None
1003
1003
1004 nLags = None
1004 nLags = None
1005
1005
1006 nPairs = None
1006 nPairs = None
1007
1007
1008 nAvg = None
1008 nAvg = None
1009
1009
1010 def __init__(self):
1010 def __init__(self):
1011 '''
1011 '''
1012 Constructor
1012 Constructor
1013 '''
1013 '''
1014 self.radarControllerHeaderObj = RadarControllerHeader()
1014 self.radarControllerHeaderObj = RadarControllerHeader()
1015
1015
1016 self.systemHeaderObj = SystemHeader()
1016 self.systemHeaderObj = SystemHeader()
1017
1017
1018 self.type = "Correlation"
1018 self.type = "Correlation"
1019
1019
1020 self.data = None
1020 self.data = None
1021
1021
1022 self.dtype = None
1022 self.dtype = None
1023
1023
1024 self.nProfiles = None
1024 self.nProfiles = None
1025
1025
1026 self.heightList = None
1026 self.heightList = None
1027
1027
1028 self.channelList = None
1028 self.channelList = None
1029
1029
1030 self.flagNoData = True
1030 self.flagNoData = True
1031
1031
1032 self.flagDiscontinuousBlock = False
1032 self.flagDiscontinuousBlock = False
1033
1033
1034 self.utctime = None
1034 self.utctime = None
1035
1035
1036 self.timeZone = None
1036 self.timeZone = None
1037
1037
1038 self.dstFlag = None
1038 self.dstFlag = None
1039
1039
1040 self.errorCount = None
1040 self.errorCount = None
1041
1041
1042 self.blocksize = None
1042 self.blocksize = None
1043
1043
1044 self.flagDecodeData = False # asumo q la data no esta decodificada
1044 self.flagDecodeData = False # asumo q la data no esta decodificada
1045
1045
1046 self.flagDeflipData = False # asumo q la data no esta sin flip
1046 self.flagDeflipData = False # asumo q la data no esta sin flip
1047
1047
1048 self.pairsList = None
1048 self.pairsList = None
1049
1049
1050 self.nPoints = None
1050 self.nPoints = None
1051
1051
1052 def getPairsList(self):
1052 def getPairsList(self):
1053
1053
1054 return self.pairsList
1054 return self.pairsList
1055
1055
1056 def getNoise(self, mode=2):
1056 def getNoise(self, mode=2):
1057
1057
1058 indR = numpy.where(self.lagR == 0)[0][0]
1058 indR = numpy.where(self.lagR == 0)[0][0]
1059 indT = numpy.where(self.lagT == 0)[0][0]
1059 indT = numpy.where(self.lagT == 0)[0][0]
1060
1060
1061 jspectra0 = self.data_corr[:, :, indR, :]
1061 jspectra0 = self.data_corr[:, :, indR, :]
1062 jspectra = copy.copy(jspectra0)
1062 jspectra = copy.copy(jspectra0)
1063
1063
1064 num_chan = jspectra.shape[0]
1064 num_chan = jspectra.shape[0]
1065 num_hei = jspectra.shape[2]
1065 num_hei = jspectra.shape[2]
1066
1066
1067 freq_dc = jspectra.shape[1] / 2
1067 freq_dc = jspectra.shape[1] / 2
1068 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
1068 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
1069
1069
1070 if ind_vel[0] < 0:
1070 if ind_vel[0] < 0:
1071 ind_vel[range(0, 1)] = ind_vel[range(0, 1)] + self.num_prof
1071 ind_vel[range(0, 1)] = ind_vel[range(0, 1)] + self.num_prof
1072
1072
1073 if mode == 1:
1073 if mode == 1:
1074 jspectra[:, freq_dc, :] = (
1074 jspectra[:, freq_dc, :] = (
1075 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
1075 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
1076
1076
1077 if mode == 2:
1077 if mode == 2:
1078
1078
1079 vel = numpy.array([-2, -1, 1, 2])
1079 vel = numpy.array([-2, -1, 1, 2])
1080 xx = numpy.zeros([4, 4])
1080 xx = numpy.zeros([4, 4])
1081
1081
1082 for fil in range(4):
1082 for fil in range(4):
1083 xx[fil, :] = vel[fil]**numpy.asarray(range(4))
1083 xx[fil, :] = vel[fil]**numpy.asarray(range(4))
1084
1084
1085 xx_inv = numpy.linalg.inv(xx)
1085 xx_inv = numpy.linalg.inv(xx)
1086 xx_aux = xx_inv[0, :]
1086 xx_aux = xx_inv[0, :]
1087
1087
1088 for ich in range(num_chan):
1088 for ich in range(num_chan):
1089 yy = jspectra[ich, ind_vel, :]
1089 yy = jspectra[ich, ind_vel, :]
1090 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
1090 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
1091
1091
1092 junkid = jspectra[ich, freq_dc, :] <= 0
1092 junkid = jspectra[ich, freq_dc, :] <= 0
1093 cjunkid = sum(junkid)
1093 cjunkid = sum(junkid)
1094
1094
1095 if cjunkid.any():
1095 if cjunkid.any():
1096 jspectra[ich, freq_dc, junkid.nonzero()] = (
1096 jspectra[ich, freq_dc, junkid.nonzero()] = (
1097 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
1097 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
1098
1098
1099 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
1099 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
1100
1100
1101 return noise
1101 return noise
1102
1102
1103 def getTimeInterval(self):
1103 def getTimeInterval(self):
1104
1104
1105 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1105 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1106
1106
1107 return timeInterval
1107 return timeInterval
1108
1108
1109 def splitFunctions(self):
1109 def splitFunctions(self):
1110
1110
1111 pairsList = self.pairsList
1111 pairsList = self.pairsList
1112 ccf_pairs = []
1112 ccf_pairs = []
1113 acf_pairs = []
1113 acf_pairs = []
1114 ccf_ind = []
1114 ccf_ind = []
1115 acf_ind = []
1115 acf_ind = []
1116 for l in range(len(pairsList)):
1116 for l in range(len(pairsList)):
1117 chan0 = pairsList[l][0]
1117 chan0 = pairsList[l][0]
1118 chan1 = pairsList[l][1]
1118 chan1 = pairsList[l][1]
1119
1119
1120 # Obteniendo pares de Autocorrelacion
1120 # Obteniendo pares de Autocorrelacion
1121 if chan0 == chan1:
1121 if chan0 == chan1:
1122 acf_pairs.append(chan0)
1122 acf_pairs.append(chan0)
1123 acf_ind.append(l)
1123 acf_ind.append(l)
1124 else:
1124 else:
1125 ccf_pairs.append(pairsList[l])
1125 ccf_pairs.append(pairsList[l])
1126 ccf_ind.append(l)
1126 ccf_ind.append(l)
1127
1127
1128 data_acf = self.data_cf[acf_ind]
1128 data_acf = self.data_cf[acf_ind]
1129 data_ccf = self.data_cf[ccf_ind]
1129 data_ccf = self.data_cf[ccf_ind]
1130
1130
1131 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1131 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1132
1132
1133 def getNormFactor(self):
1133 def getNormFactor(self):
1134 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1134 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1135 acf_pairs = numpy.array(acf_pairs)
1135 acf_pairs = numpy.array(acf_pairs)
1136 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1136 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1137
1137
1138 for p in range(self.nPairs):
1138 for p in range(self.nPairs):
1139 pair = self.pairsList[p]
1139 pair = self.pairsList[p]
1140
1140
1141 ch0 = pair[0]
1141 ch0 = pair[0]
1142 ch1 = pair[1]
1142 ch1 = pair[1]
1143
1143
1144 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1144 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1145 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1145 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1146 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1146 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1147
1147
1148 return normFactor
1148 return normFactor
1149
1149
1150 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1150 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1151 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1151 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1152
1152
1153
1153
1154 class Parameters(Spectra):
1154 class Parameters(Spectra):
1155
1155
1156 experimentInfo = None # Information about the experiment
1156 experimentInfo = None # Information about the experiment
1157
1157
1158 # Information from previous data
1158 # Information from previous data
1159
1159
1160 inputUnit = None # Type of data to be processed
1160 inputUnit = None # Type of data to be processed
1161
1161
1162 operation = None # Type of operation to parametrize
1162 operation = None # Type of operation to parametrize
1163
1163
1164 # normFactor = None #Normalization Factor
1164 # normFactor = None #Normalization Factor
1165
1165
1166 groupList = None # List of Pairs, Groups, etc
1166 groupList = None # List of Pairs, Groups, etc
1167
1167
1168 # Parameters
1168 # Parameters
1169
1169
1170 data_param = None # Parameters obtained
1170 data_param = None # Parameters obtained
1171
1171
1172 data_pre = None # Data Pre Parametrization
1172 data_pre = None # Data Pre Parametrization
1173
1173
1174 data_SNR = None # Signal to Noise Ratio
1174 data_SNR = None # Signal to Noise Ratio
1175
1175
1176 # heightRange = None #Heights
1176 # heightRange = None #Heights
1177
1177
1178 abscissaList = None # Abscissa, can be velocities, lags or time
1178 abscissaList = None # Abscissa, can be velocities, lags or time
1179
1179
1180 # noise = None #Noise Potency
1180 # noise = None #Noise Potency
1181
1181
1182 utctimeInit = None # Initial UTC time
1182 utctimeInit = None # Initial UTC time
1183
1183
1184 paramInterval = None # Time interval to calculate Parameters in seconds
1184 paramInterval = None # Time interval to calculate Parameters in seconds
1185
1185
1186 useLocalTime = True
1186 useLocalTime = True
1187
1187
1188 # Fitting
1188 # Fitting
1189
1189
1190 data_error = None # Error of the estimation
1190 data_error = None # Error of the estimation
1191
1191
1192 constants = None
1192 constants = None
1193
1193
1194 library = None
1194 library = None
1195
1195
1196 # Output signal
1196 # Output signal
1197
1197
1198 outputInterval = None # Time interval to calculate output signal in seconds
1198 outputInterval = None # Time interval to calculate output signal in seconds
1199
1199
1200 data_output = None # Out signal
1200 data_output = None # Out signal
1201
1201
1202 nAvg = None
1202 nAvg = None
1203
1203
1204 noise_estimation = None
1204 noise_estimation = None
1205
1205
1206 GauSPC = None # Fit gaussian SPC
1206 GauSPC = None # Fit gaussian SPC
1207
1207
1208 def __init__(self):
1208 def __init__(self):
1209 '''
1209 '''
1210 Constructor
1210 Constructor
1211 '''
1211 '''
1212 self.radarControllerHeaderObj = RadarControllerHeader()
1212 self.radarControllerHeaderObj = RadarControllerHeader()
1213
1213
1214 self.systemHeaderObj = SystemHeader()
1214 self.systemHeaderObj = SystemHeader()
1215
1215
1216 self.type = "Parameters"
1216 self.type = "Parameters"
1217
1217
1218 def getTimeRange1(self, interval):
1218 def getTimeRange1(self, interval):
1219
1219
1220 datatime = []
1220 datatime = []
1221
1221
1222 if self.useLocalTime:
1222 if self.useLocalTime:
1223 time1 = self.utctimeInit - self.timeZone * 60
1223 time1 = self.utctimeInit - self.timeZone * 60
1224 else:
1224 else:
1225 time1 = self.utctimeInit
1225 time1 = self.utctimeInit
1226
1226
1227 datatime.append(time1)
1227 datatime.append(time1)
1228 datatime.append(time1 + interval)
1228 datatime.append(time1 + interval)
1229 datatime = numpy.array(datatime)
1229 datatime = numpy.array(datatime)
1230
1230
1231 return datatime
1231 return datatime
1232
1232
1233 def getTimeInterval(self):
1233 def getTimeInterval(self):
1234
1234
1235 if hasattr(self, 'timeInterval1'):
1235 if hasattr(self, 'timeInterval1'):
1236 return self.timeInterval1
1236 return self.timeInterval1
1237 else:
1237 else:
1238 return self.paramInterval
1238 return self.paramInterval
1239
1239
1240 def setValue(self, value):
1240 def setValue(self, value):
1241
1241
1242 print "This property should not be initialized"
1242 print "This property should not be initialized"
1243
1243
1244 return
1244 return
1245
1245
1246 def getNoise(self):
1246 def getNoise(self):
1247
1247
1248 return self.spc_noise
1248 return self.spc_noise
1249
1249
1250 timeInterval = property(getTimeInterval)
1250 timeInterval = property(getTimeInterval)
1251 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1251 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
General Comments 0
You need to be logged in to leave comments. Login now