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