##// END OF EJS Templates
schain-jc
RSS
formatting y raise cuando digitalrf recibe spectra
José Chávez -
r1120:4cb1c6729a0e
parent child
Show More
Side by Side Unified
Context file:
r1120:4cb1c6729a0e Loading diff...:
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (659 lines changed) Show them Hide them
@@ -13,75 +13,77 SPEED_OF_LIGHT = 299792458
13 13 SPEED_OF_LIGHT = 3e8
14 14
15 15 BASIC_STRUCTURE = numpy.dtype([
16 ('nSize','<u4'),
17 ('nVersion','<u2'),
18 ('nDataBlockId','<u4'),
19 ('nUtime','<u4'),
20 ('nMilsec','<u2'),
21 ('nTimezone','<i2'),
22 ('nDstflag','<i2'),
23 ('nErrorCount','<u4')
24 ])
16 ('nSize', '<u4'),
17 ('nVersion', '<u2'),
18 ('nDataBlockId', '<u4'),
19 ('nUtime', '<u4'),
20 ('nMilsec', '<u2'),
21 ('nTimezone', '<i2'),
22 ('nDstflag', '<i2'),
23 ('nErrorCount', '<u4')
24 ])
25 25
26 26 SYSTEM_STRUCTURE = numpy.dtype([
27 ('nSize','<u4'),
28 ('nNumSamples','<u4'),
29 ('nNumProfiles','<u4'),
30 ('nNumChannels','<u4'),
31 ('nADCResolution','<u4'),
32 ('nPCDIOBusWidth','<u4'),
33 ])
27 ('nSize', '<u4'),
28 ('nNumSamples', '<u4'),
29 ('nNumProfiles', '<u4'),
30 ('nNumChannels', '<u4'),
31 ('nADCResolution', '<u4'),
32 ('nPCDIOBusWidth', '<u4'),
33 ])
34 34
35 35 RADAR_STRUCTURE = numpy.dtype([
36 ('nSize','<u4'),
37 ('nExpType','<u4'),
38 ('nNTx','<u4'),
39 ('fIpp','<f4'),
40 ('fTxA','<f4'),
41 ('fTxB','<f4'),
42 ('nNumWindows','<u4'),
43 ('nNumTaus','<u4'),
44 ('nCodeType','<u4'),
45 ('nLine6Function','<u4'),
46 ('nLine5Function','<u4'),
47 ('fClock','<f4'),
48 ('nPrePulseBefore','<u4'),
49 ('nPrePulseAfter','<u4'),
50 ('sRangeIPP','<a20'),
51 ('sRangeTxA','<a20'),
52 ('sRangeTxB','<a20'),
53 ])
54
55 SAMPLING_STRUCTURE = numpy.dtype([('h0','<f4'),('dh','<f4'),('nsa','<u4')])
36 ('nSize', '<u4'),
37 ('nExpType', '<u4'),
38 ('nNTx', '<u4'),
39 ('fIpp', '<f4'),
40 ('fTxA', '<f4'),
41 ('fTxB', '<f4'),
42 ('nNumWindows', '<u4'),
43 ('nNumTaus', '<u4'),
44 ('nCodeType', '<u4'),
45 ('nLine6Function', '<u4'),
46 ('nLine5Function', '<u4'),
47 ('fClock', '<f4'),
48 ('nPrePulseBefore', '<u4'),
49 ('nPrePulseAfter', '<u4'),
50 ('sRangeIPP', '<a20'),
51 ('sRangeTxA', '<a20'),
52 ('sRangeTxB', '<a20'),
53 ])
54
55 SAMPLING_STRUCTURE = numpy.dtype(
56 [('h0', '<f4'), ('dh', '<f4'), ('nsa', '<u4')])
56 57
57 58
58 59 PROCESSING_STRUCTURE = numpy.dtype([
59 ('nSize','<u4'),
60 ('nDataType','<u4'),
61 ('nSizeOfDataBlock','<u4'),
62 ('nProfilesperBlock','<u4'),
63 ('nDataBlocksperFile','<u4'),
64 ('nNumWindows','<u4'),
65 ('nProcessFlags','<u4'),
66 ('nCoherentIntegrations','<u4'),
67 ('nIncoherentIntegrations','<u4'),
68 ('nTotalSpectra','<u4')
69 ])
60 ('nSize', '<u4'),
61 ('nDataType', '<u4'),
62 ('nSizeOfDataBlock', '<u4'),
63 ('nProfilesperBlock', '<u4'),
64 ('nDataBlocksperFile', '<u4'),
65 ('nNumWindows', '<u4'),
66 ('nProcessFlags', '<u4'),
67 ('nCoherentIntegrations', '<u4'),
68 ('nIncoherentIntegrations', '<u4'),
69 ('nTotalSpectra', '<u4')
70 ])
71
70 72
71 73 class Header(object):
72
74
73 75 def __init__(self):
74 76 raise NotImplementedError
75
77
76 78 def copy(self):
77 79 return copy.deepcopy(self)
78
80
79 81 def read(self):
80
82
81 83 raise NotImplementedError
82
84
83 85 def write(self):
84
86
85 87 raise NotImplementedError
86 88
87 89 def getAllowedArgs(self):
@@ -98,32 +100,33 class Header(object):
98 100 for x in args:
99 101 asDict[x] = self[x]
100 102 return asDict
101
103
102 104 def __getitem__(self, name):
103 105 return getattr(self, name)
104 106
105 107 def printInfo(self):
106
107 message = "#"*50 + "\n"
108
109 message = "#" * 50 + "\n"
108 110 message += self.__class__.__name__.upper() + "\n"
109 message += "#"*50 + "\n"
110
111 message += "#" * 50 + "\n"
112
111 113 keyList = self.__dict__.keys()
112 114 keyList.sort()
113
115
114 116 for key in keyList:
115 message += "%s = %s" %(key, self.__dict__[key]) + "\n"
116
117 message += "%s = %s" % (key, self.__dict__[key]) + "\n"
118
117 119 if "size" not in keyList:
118 120 attr = getattr(self, "size")
119
121
120 122 if attr:
121 message += "%s = %s" %("size", attr) + "\n"
122
123 message += "%s = %s" % ("size", attr) + "\n"
124
123 125 print message
124 126
127
125 128 class BasicHeader(Header):
126
129
127 130 size = None
128 131 version = None
129 132 dataBlock = None
@@ -136,9 +139,9 class BasicHeader(Header):
136 139 datatime = None
137 140 structure = BASIC_STRUCTURE
138 141 __LOCALTIME = None
139
142
140 143 def __init__(self, useLocalTime=True):
141
144
142 145 self.size = 24
143 146 self.version = 0
144 147 self.dataBlock = 0
@@ -147,17 +150,17 class BasicHeader(Header):
147 150 self.timeZone = 0
148 151 self.dstFlag = 0
149 152 self.errorCount = 0
150
153
151 154 self.useLocalTime = useLocalTime
152
155
153 156 def read(self, fp):
154
157
155 158 self.length = 0
156 159 try:
157 160 if hasattr(fp, 'read'):
158 header = numpy.fromfile(fp, BASIC_STRUCTURE,1)
161 header = numpy.fromfile(fp, BASIC_STRUCTURE, 1)
159 162 else:
160 header = numpy.fromstring(fp, BASIC_STRUCTURE,1)
163 header = numpy.fromstring(fp, BASIC_STRUCTURE, 1)
161 164 except Exception, e:
162 165 print "BasicHeader: "
163 166 print e
@@ -171,38 +174,40 class BasicHeader(Header):
171 174 self.timeZone = int(header['nTimezone'][0])
172 175 self.dstFlag = int(header['nDstflag'][0])
173 176 self.errorCount = int(header['nErrorCount'][0])
174
177
175 178 if self.size < 24:
176 179 return 0
177 180
178 181 self.length = header.nbytes
179 182 return 1
180
183
181 184 def write(self, fp):
182
183 headerTuple = (self.size,self.version,self.dataBlock,self.utc,self.miliSecond,self.timeZone,self.dstFlag,self.errorCount)
184 header = numpy.array(headerTuple, BASIC_STRUCTURE)
185
186 headerTuple = (self.size, self.version, self.dataBlock, self.utc,
187 self.miliSecond, self.timeZone, self.dstFlag, self.errorCount)
188 header = numpy.array(headerTuple, BASIC_STRUCTURE)
185 189 header.tofile(fp)
186
190
187 191 return 1
188
192
189 193 def get_ltc(self):
190
191 return self.utc - self.timeZone*60
192
194
195 return self.utc - self.timeZone * 60
196
193 197 def set_ltc(self, value):
194
195 self.utc = value + self.timeZone*60
196
198
199 self.utc = value + self.timeZone * 60
200
197 201 def get_datatime(self):
198
202
199 203 return datetime.datetime.utcfromtimestamp(self.ltc)
200
204
201 205 ltc = property(get_ltc, set_ltc)
202 206 datatime = property(get_datatime)
203 207
208
204 209 class SystemHeader(Header):
205
210
206 211 size = None
207 212 nSamples = None
208 213 nProfiles = None
@@ -212,8 +217,8 class SystemHeader(Header):
212 217 structure = SYSTEM_STRUCTURE
213 218
214 219 def __init__(self, nSamples=0, nProfiles=0, nChannels=0, adcResolution=14, pciDioBusWidth=0):
215
216 self.size = 24
220
221 self.size = 24
217 222 self.nSamples = nSamples
218 223 self.nProfiles = nProfiles
219 224 self.nChannels = nChannels
@@ -230,45 +235,48 class SystemHeader(Header):
230 235
231 236 try:
232 237 if hasattr(fp, 'read'):
233 header = numpy.fromfile(fp, SYSTEM_STRUCTURE,1)
238 header = numpy.fromfile(fp, SYSTEM_STRUCTURE, 1)
234 239 else:
235 header = numpy.fromstring(fp, SYSTEM_STRUCTURE,1)
240 header = numpy.fromstring(fp, SYSTEM_STRUCTURE, 1)
236 241 except Exception, e:
237 242 print "System Header: " + str(e)
238 243 return 0
239
244
240 245 self.size = header['nSize'][0]
241 246 self.nSamples = header['nNumSamples'][0]
242 247 self.nProfiles = header['nNumProfiles'][0]
243 248 self.nChannels = header['nNumChannels'][0]
244 249 self.adcResolution = header['nADCResolution'][0]
245 250 self.pciDioBusWidth = header['nPCDIOBusWidth'][0]
246
247
251
248 252 if startFp is not None:
249 253 endFp = self.size + startFp
250
254
251 255 if fp.tell() > endFp:
252 sys.stderr.write("Warning %s: Size value read from System Header is lower than it has to be\n" %fp.name)
256 sys.stderr.write(
257 "Warning %s: Size value read from System Header is lower than it has to be\n" % fp.name)
253 258 return 0
254
259
255 260 if fp.tell() < endFp:
256 sys.stderr.write("Warning %s: Size value read from System Header size is greater than it has to be\n" %fp.name)
261 sys.stderr.write(
262 "Warning %s: Size value read from System Header size is greater than it has to be\n" % fp.name)
257 263 return 0
258
264
259 265 self.length = header.nbytes
260 266 return 1
261
267
262 268 def write(self, fp):
263
264 headerTuple = (self.size,self.nSamples,self.nProfiles,self.nChannels,self.adcResolution,self.pciDioBusWidth)
265 header = numpy.array(headerTuple,SYSTEM_STRUCTURE)
269
270 headerTuple = (self.size, self.nSamples, self.nProfiles,
271 self.nChannels, self.adcResolution, self.pciDioBusWidth)
272 header = numpy.array(headerTuple, SYSTEM_STRUCTURE)
266 273 header.tofile(fp)
267
274
268 275 return 1
269 276
277
270 278 class RadarControllerHeader(Header):
271
279
272 280 expType = None
273 281 nTx = None
274 282 ipp = None
@@ -287,7 +295,7 class RadarControllerHeader(Header):
287 295 rangeTxB = None
288 296 structure = RADAR_STRUCTURE
289 297 __size = None
290
298
291 299 def __init__(self, expType=2, nTx=1,
292 300 ipp=None, txA=0, txB=0,
293 301 nWindows=None, nHeights=None, firstHeight=None, deltaHeight=None,
@@ -295,8 +303,8 class RadarControllerHeader(Header):
295 303 prePulseBefore=0, prePulseAfter=0,
296 304 codeType=0, nCode=0, nBaud=0, code=None,
297 305 flip1=0, flip2=0):
298
299 # self.size = 116
306
307 # self.size = 116
300 308 self.expType = expType
301 309 self.nTx = nTx
302 310 self.ipp = ipp
@@ -305,7 +313,7 class RadarControllerHeader(Header):
305 313 self.rangeIpp = ipp
306 314 self.rangeTxA = txA
307 315 self.rangeTxB = txB
308
316
309 317 self.nWindows = nWindows
310 318 self.numTaus = numTaus
311 319 self.codeType = codeType
@@ -314,23 +322,23 class RadarControllerHeader(Header):
314 322 self.fClock = fClock
315 323 self.prePulseBefore = prePulseBefore
316 324 self.prePulseAfter = prePulseAfter
317
325
318 326 self.nHeights = nHeights
319 327 self.firstHeight = firstHeight
320 328 self.deltaHeight = deltaHeight
321 329 self.samplesWin = nHeights
322
330
323 331 self.nCode = nCode
324 332 self.nBaud = nBaud
325 333 self.code = code
326 334 self.flip1 = flip1
327 335 self.flip2 = flip2
328
329 self.code_size = int(numpy.ceil(self.nBaud/32.))*self.nCode*4
336
337 self.code_size = int(numpy.ceil(self.nBaud / 32.)) * self.nCode * 4
330 338 # self.dynamic = numpy.array([],numpy.dtype('byte'))
331
339
332 340 if self.fClock is None and self.deltaHeight is not None:
333 self.fClock = 0.15/(deltaHeight*1e-6) #0.15Km / (height * 1u)
341 self.fClock = 0.15 / (deltaHeight * 1e-6) # 0.15Km / (height * 1u)
334 342
335 343 def read(self, fp):
336 344 self.length = 0
@@ -342,14 +350,14 class RadarControllerHeader(Header):
342 350
343 351 try:
344 352 if hasattr(fp, 'read'):
345 header = numpy.fromfile(fp, RADAR_STRUCTURE,1)
353 header = numpy.fromfile(fp, RADAR_STRUCTURE, 1)
346 354 else:
347 header = numpy.fromstring(fp, RADAR_STRUCTURE,1)
355 header = numpy.fromstring(fp, RADAR_STRUCTURE, 1)
348 356 self.length += header.nbytes
349 357 except Exception, e:
350 358 print "RadarControllerHeader: " + str(e)
351 359 return 0
352
360
353 361 size = int(header['nSize'][0])
354 362 self.expType = int(header['nExpType'][0])
355 363 self.nTx = int(header['nNTx'][0])
@@ -367,12 +375,14 class RadarControllerHeader(Header):
367 375 self.rangeIpp = header['sRangeIPP'][0]
368 376 self.rangeTxA = header['sRangeTxA'][0]
369 377 self.rangeTxB = header['sRangeTxB'][0]
370
378
371 379 try:
372 380 if hasattr(fp, 'read'):
373 samplingWindow = numpy.fromfile(fp, SAMPLING_STRUCTURE, self.nWindows)
381 samplingWindow = numpy.fromfile(
382 fp, SAMPLING_STRUCTURE, self.nWindows)
374 383 else:
375 samplingWindow = numpy.fromstring(fp[self.length:], SAMPLING_STRUCTURE, self.nWindows)
384 samplingWindow = numpy.fromstring(
385 fp[self.length:], SAMPLING_STRUCTURE, self.nWindows)
376 386 self.length += samplingWindow.nbytes
377 387 except Exception, e:
378 388 print "RadarControllerHeader: " + str(e)
@@ -381,24 +391,21 class RadarControllerHeader(Header):
381 391 self.firstHeight = samplingWindow['h0']
382 392 self.deltaHeight = samplingWindow['dh']
383 393 self.samplesWin = samplingWindow['nsa']
384
385
386 394
387 395 try:
388 396 if hasattr(fp, 'read'):
389 397 self.Taus = numpy.fromfile(fp, '<f4', self.numTaus)
390 398 else:
391 self.Taus = numpy.fromstring(fp[self.length:], '<f4', self.numTaus)
399 self.Taus = numpy.fromstring(
400 fp[self.length:], '<f4', self.numTaus)
392 401 self.length += self.Taus.nbytes
393 402 except Exception, e:
394 403 print "RadarControllerHeader: " + str(e)
395 404 return 0
396
397
398
405
399 406 self.code_size = 0
400 407 if self.codeType != 0:
401
408
402 409 try:
403 410 if hasattr(fp, 'read'):
404 411 self.nCode = numpy.fromfile(fp, '<u4', 1)[0]
@@ -406,58 +413,63 class RadarControllerHeader(Header):
406 413 self.nBaud = numpy.fromfile(fp, '<u4', 1)[0]
407 414 self.length += self.nBaud.nbytes
408 415 else:
409 self.nCode = numpy.fromstring(fp[self.length:], '<u4', 1)[0]
416 self.nCode = numpy.fromstring(
417 fp[self.length:], '<u4', 1)[0]
410 418 self.length += self.nCode.nbytes
411 self.nBaud = numpy.fromstring(fp[self.length:], '<u4', 1)[0]
419 self.nBaud = numpy.fromstring(
420 fp[self.length:], '<u4', 1)[0]
412 421 self.length += self.nBaud.nbytes
413 422 except Exception, e:
414 423 print "RadarControllerHeader: " + str(e)
415 return 0
416 code = numpy.empty([self.nCode,self.nBaud],dtype='i1')
417
424 return 0
425 code = numpy.empty([self.nCode, self.nBaud], dtype='i1')
426
418 427 for ic in range(self.nCode):
419 428 try:
420 429 if hasattr(fp, 'read'):
421 temp = numpy.fromfile(fp,'u4', int(numpy.ceil(self.nBaud/32.)))
430 temp = numpy.fromfile(fp, 'u4', int(
431 numpy.ceil(self.nBaud / 32.)))
422 432 else:
423 temp = numpy.fromstring(fp,'u4', int(numpy.ceil(self.nBaud/32.)))
433 temp = numpy.fromstring(
434 fp, 'u4', int(numpy.ceil(self.nBaud / 32.)))
424 435 self.length += temp.nbytes
425 436 except Exception, e:
426 437 print "RadarControllerHeader: " + str(e)
427 438 return 0
428 439
429 for ib in range(self.nBaud-1,-1,-1):
430 code[ic,ib] = temp[ib/32]%2
431 temp[ib/32] = temp[ib/32]/2
432
433 self.code = 2.0*code - 1.0
434 self.code_size = int(numpy.ceil(self.nBaud/32.))*self.nCode*4
435
440 for ib in range(self.nBaud - 1, -1, -1):
441 code[ic, ib] = temp[ib / 32] % 2
442 temp[ib / 32] = temp[ib / 32] / 2
443
444 self.code = 2.0 * code - 1.0
445 self.code_size = int(numpy.ceil(self.nBaud / 32.)) * self.nCode * 4
446
436 447 # if self.line5Function == RCfunction.FLIP:
437 448 # self.flip1 = numpy.fromfile(fp,'<u4',1)
438 #
449 #
439 450 # if self.line6Function == RCfunction.FLIP:
440 451 # self.flip2 = numpy.fromfile(fp,'<u4',1)
441 452 if startFp is not None:
442 453 endFp = size + startFp
443
454
444 455 if fp.tell() != endFp:
445 # fp.seek(endFp)
446 print "%s: Radar Controller Header size is not consistent: from data [%d] != from header field [%d]" %(fp.name, fp.tell()-startFp, size)
456 # fp.seek(endFp)
457 print "%s: Radar Controller Header size is not consistent: from data [%d] != from header field [%d]" % (fp.name, fp.tell() - startFp, size)
447 458 # return 0
448 459
449 460 if fp.tell() > endFp:
450 sys.stderr.write("Warning %s: Size value read from Radar Controller header is lower than it has to be\n" %fp.name)
461 sys.stderr.write(
462 "Warning %s: Size value read from Radar Controller header is lower than it has to be\n" % fp.name)
451 463 # return 0
452
464
453 465 if fp.tell() < endFp:
454 sys.stderr.write("Warning %s: Size value read from Radar Controller header is greater than it has to be\n" %fp.name)
455
456
466 sys.stderr.write(
467 "Warning %s: Size value read from Radar Controller header is greater than it has to be\n" % fp.name)
468
457 469 return 1
458
470
459 471 def write(self, fp):
460
472
461 473 headerTuple = (self.size,
462 474 self.expType,
463 475 self.nTx,
@@ -475,83 +487,87 class RadarControllerHeader(Header):
475 487 self.rangeIpp,
476 488 self.rangeTxA,
477 489 self.rangeTxB)
478
479 header = numpy.array(headerTuple,RADAR_STRUCTURE)
490
491 header = numpy.array(headerTuple, RADAR_STRUCTURE)
480 492 header.tofile(fp)
481
482 sampleWindowTuple = (self.firstHeight,self.deltaHeight,self.samplesWin)
483 samplingWindow = numpy.array(sampleWindowTuple,SAMPLING_STRUCTURE)
493
494 sampleWindowTuple = (
495 self.firstHeight, self.deltaHeight, self.samplesWin)
496 samplingWindow = numpy.array(sampleWindowTuple, SAMPLING_STRUCTURE)
484 497 samplingWindow.tofile(fp)
485
498
486 499 if self.numTaus > 0:
487 500 self.Taus.tofile(fp)
488
489 if self.codeType !=0:
501
502 if self.codeType != 0:
490 503 nCode = numpy.array(self.nCode, '<u4')
491 504 nCode.tofile(fp)
492 505 nBaud = numpy.array(self.nBaud, '<u4')
493 506 nBaud.tofile(fp)
494 code1 = (self.code + 1.0)/2.
495
507 code1 = (self.code + 1.0) / 2.
508
496 509 for ic in range(self.nCode):
497 tempx = numpy.zeros(numpy.ceil(self.nBaud/32.))
510 tempx = numpy.zeros(numpy.ceil(self.nBaud / 32.))
498 511 start = 0
499 512 end = 32
500 513 for i in range(len(tempx)):
501 code_selected = code1[ic,start:end]
502 for j in range(len(code_selected)-1,-1,-1):
514 code_selected = code1[ic, start:end]
515 for j in range(len(code_selected) - 1, -1, -1):
503 516 if code_selected[j] == 1:
504 tempx[i] = tempx[i] + 2**(len(code_selected)-1-j)
517 tempx[i] = tempx[i] + \
518 2**(len(code_selected) - 1 - j)
505 519 start = start + 32
506 520 end = end + 32
507
521
508 522 tempx = tempx.astype('u4')
509 523 tempx.tofile(fp)
510
524
511 525 # if self.line5Function == RCfunction.FLIP:
512 526 # self.flip1.tofile(fp)
513 #
527 #
514 528 # if self.line6Function == RCfunction.FLIP:
515 529 # self.flip2.tofile(fp)
516
530
517 531 return 1
518 532
519 533 def get_ippSeconds(self):
520 534 '''
521 535 '''
522 536 ippSeconds = 2.0 * 1000 * self.ipp / SPEED_OF_LIGHT
523
537
524 538 return ippSeconds
525
539
526 540 def set_ippSeconds(self, ippSeconds):
527 541 '''
528 542 '''
529
530 self.ipp = ippSeconds * SPEED_OF_LIGHT / (2.0*1000)
531
543
544 self.ipp = ippSeconds * SPEED_OF_LIGHT / (2.0 * 1000)
545
532 546 return
533
547
534 548 def get_size(self):
535
536 self.__size = 116 + 12*self.nWindows + 4*self.numTaus
537
549
550 self.__size = 116 + 12 * self.nWindows + 4 * self.numTaus
551
538 552 if self.codeType != 0:
539 self.__size += 4 + 4 + 4*self.nCode*numpy.ceil(self.nBaud/32.)
540
553 self.__size += 4 + 4 + 4 * self.nCode * \
554 numpy.ceil(self.nBaud / 32.)
555
541 556 return self.__size
542
557
543 558 def set_size(self, value):
544
559
545 560 raise IOError, "size is a property and it cannot be set, just read"
546
561
547 562 return
548
563
549 564 ippSeconds = property(get_ippSeconds, set_ippSeconds)
550 565 size = property(get_size, set_size)
551 566
567
552 568 class ProcessingHeader(Header):
553
554 # size = None
569
570 # size = None
555 571 dtype = None
556 572 blockSize = None
557 573 profilesPerBlock = None
@@ -564,13 +580,13 class ProcessingHeader(Header):
564 580 structure = PROCESSING_STRUCTURE
565 581 flag_dc = None
566 582 flag_cspc = None
567
568 def __init__(self, dtype=0, blockSize=0, profilesPerBlock=0, dataBlocksPerFile=0, nWindows=0,processFlags=0, nCohInt=0,
569 nIncohInt=0, totalSpectra=0, nHeights=0, firstHeight=0, deltaHeight=0, samplesWin=0, spectraComb=0, nCode=0,
570 code=0, nBaud=None, shif_fft=False, flag_dc=False, flag_cspc=False, flag_decode=False, flag_deflip=False
571 ):
572
573 # self.size = 0
583
584 def __init__(self, dtype=0, blockSize=0, profilesPerBlock=0, dataBlocksPerFile=0, nWindows=0, processFlags=0, nCohInt=0,
585 nIncohInt=0, totalSpectra=0, nHeights=0, firstHeight=0, deltaHeight=0, samplesWin=0, spectraComb=0, nCode=0,
586 code=0, nBaud=None, shif_fft=False, flag_dc=False, flag_cspc=False, flag_decode=False, flag_deflip=False
587 ):
588
589 # self.size = 0
574 590 self.dtype = dtype
575 591 self.blockSize = blockSize
576 592 self.profilesPerBlock = 0
@@ -580,7 +596,7 class ProcessingHeader(Header):
580 596 self.nCohInt = 0
581 597 self.nIncohInt = 0
582 598 self.totalSpectra = 0
583
599
584 600 self.nHeights = 0
585 601 self.firstHeight = 0
586 602 self.deltaHeight = 0
@@ -589,7 +605,7 class ProcessingHeader(Header):
589 605 self.nCode = None
590 606 self.code = None
591 607 self.nBaud = None
592
608
593 609 self.shif_fft = False
594 610 self.flag_dc = False
595 611 self.flag_cspc = False
@@ -604,7 +620,7 class ProcessingHeader(Header):
604 620 except Exception, e:
605 621 startFp = None
606 622 pass
607
623
608 624 try:
609 625 if hasattr(fp, 'read'):
610 626 header = numpy.fromfile(fp, PROCESSING_STRUCTURE, 1)
@@ -614,7 +630,7 class ProcessingHeader(Header):
614 630 except Exception, e:
615 631 print "ProcessingHeader: " + str(e)
616 632 return 0
617
633
618 634 size = int(header['nSize'][0])
619 635 self.dtype = int(header['nDataType'][0])
620 636 self.blockSize = int(header['nSizeOfDataBlock'][0])
@@ -625,94 +641,98 class ProcessingHeader(Header):
625 641 self.nCohInt = int(header['nCoherentIntegrations'][0])
626 642 self.nIncohInt = int(header['nIncoherentIntegrations'][0])
627 643 self.totalSpectra = int(header['nTotalSpectra'][0])
628
644
629 645 try:
630 646 if hasattr(fp, 'read'):
631 samplingWindow = numpy.fromfile(fp, SAMPLING_STRUCTURE, self.nWindows)
647 samplingWindow = numpy.fromfile(
648 fp, SAMPLING_STRUCTURE, self.nWindows)
632 649 else:
633 samplingWindow = numpy.fromstring(fp[self.length:], SAMPLING_STRUCTURE, self.nWindows)
650 samplingWindow = numpy.fromstring(
651 fp[self.length:], SAMPLING_STRUCTURE, self.nWindows)
634 652 self.length += samplingWindow.nbytes
635 653 except Exception, e:
636 654 print "ProcessingHeader: " + str(e)
637 655 return 0
638
656
639 657 self.nHeights = int(numpy.sum(samplingWindow['nsa']))
640 658 self.firstHeight = float(samplingWindow['h0'][0])
641 659 self.deltaHeight = float(samplingWindow['dh'][0])
642 660 self.samplesWin = samplingWindow['nsa'][0]
643
644 661
645 662 try:
646 663 if hasattr(fp, 'read'):
647 self.spectraComb = numpy.fromfile(fp, 'u1', 2*self.totalSpectra)
664 self.spectraComb = numpy.fromfile(
665 fp, 'u1', 2 * self.totalSpectra)
648 666 else:
649 self.spectraComb = numpy.fromstring(fp[self.length:], 'u1', 2*self.totalSpectra)
667 self.spectraComb = numpy.fromstring(
668 fp[self.length:], 'u1', 2 * self.totalSpectra)
650 669 self.length += self.spectraComb.nbytes
651 670 except Exception, e:
652 671 print "ProcessingHeader: " + str(e)
653 672 return 0
654
673
655 674 if ((self.processFlags & PROCFLAG.DEFINE_PROCESS_CODE) == PROCFLAG.DEFINE_PROCESS_CODE):
656 self.nCode = int(numpy.fromfile(fp,'<u4',1))
657 self.nBaud = int(numpy.fromfile(fp,'<u4',1))
658 self.code = numpy.fromfile(fp,'<f4',self.nCode*self.nBaud).reshape(self.nCode,self.nBaud)
659
675 self.nCode = int(numpy.fromfile(fp, '<u4', 1))
676 self.nBaud = int(numpy.fromfile(fp, '<u4', 1))
677 self.code = numpy.fromfile(
678 fp, '<f4', self.nCode * self.nBaud).reshape(self.nCode, self.nBaud)
679
660 680 if ((self.processFlags & PROCFLAG.EXP_NAME_ESP) == PROCFLAG.EXP_NAME_ESP):
661 exp_name_len = int(numpy.fromfile(fp,'<u4',1))
662 exp_name = numpy.fromfile(fp,'u1',exp_name_len+1)
663
681 exp_name_len = int(numpy.fromfile(fp, '<u4', 1))
682 exp_name = numpy.fromfile(fp, 'u1', exp_name_len + 1)
683
664 684 if ((self.processFlags & PROCFLAG.SHIFT_FFT_DATA) == PROCFLAG.SHIFT_FFT_DATA):
665 685 self.shif_fft = True
666 686 else:
667 687 self.shif_fft = False
668
688
669 689 if ((self.processFlags & PROCFLAG.SAVE_CHANNELS_DC) == PROCFLAG.SAVE_CHANNELS_DC):
670 690 self.flag_dc = True
671 691 else:
672 692 self.flag_dc = False
673
693
674 694 if ((self.processFlags & PROCFLAG.DECODE_DATA) == PROCFLAG.DECODE_DATA):
675 695 self.flag_decode = True
676 696 else:
677 697 self.flag_decode = False
678
698
679 699 if ((self.processFlags & PROCFLAG.DEFLIP_DATA) == PROCFLAG.DEFLIP_DATA):
680 700 self.flag_deflip = True
681 701 else:
682 702 self.flag_deflip = False
683
703
684 704 nChannels = 0
685 705 nPairs = 0
686 706 pairList = []
687
688 for i in range( 0, self.totalSpectra*2, 2 ):
689 if self.spectraComb[i] == self.spectraComb[i+1]:
690 nChannels = nChannels + 1 #par de canales iguales
707
708 for i in range(0, self.totalSpectra * 2, 2):
709 if self.spectraComb[i] == self.spectraComb[i + 1]:
710 nChannels = nChannels + 1 # par de canales iguales
691 711 else:
692 nPairs = nPairs + 1 #par de canales diferentes
693 pairList.append( (self.spectraComb[i], self.spectraComb[i+1]) )
694
712 nPairs = nPairs + 1 # par de canales diferentes
713 pairList.append((self.spectraComb[i], self.spectraComb[i + 1]))
714
695 715 self.flag_cspc = False
696 716 if nPairs > 0:
697 717 self.flag_cspc = True
698
699
700
718
701 719 if startFp is not None:
702 720 endFp = size + startFp
703 721 if fp.tell() > endFp:
704 sys.stderr.write("Warning: Processing header size is lower than it has to be")
722 sys.stderr.write(
723 "Warning: Processing header size is lower than it has to be")
705 724 return 0
706
725
707 726 if fp.tell() < endFp:
708 sys.stderr.write("Warning: Processing header size is greater than it is considered")
709
727 sys.stderr.write(
728 "Warning: Processing header size is greater than it is considered")
729
710 730 return 1
711
731
712 732 def write(self, fp):
713 #Clear DEFINE_PROCESS_CODE
733 # Clear DEFINE_PROCESS_CODE
714 734 self.processFlags = self.processFlags & (~PROCFLAG.DEFINE_PROCESS_CODE)
715
735
716 736 headerTuple = (self.size,
717 737 self.dtype,
718 738 self.blockSize,
@@ -723,154 +743,163 class ProcessingHeader(Header):
723 743 self.nCohInt,
724 744 self.nIncohInt,
725 745 self.totalSpectra)
726
727 header = numpy.array(headerTuple,PROCESSING_STRUCTURE)
746
747 header = numpy.array(headerTuple, PROCESSING_STRUCTURE)
728 748 header.tofile(fp)
729
749
730 750 if self.nWindows != 0:
731 sampleWindowTuple = (self.firstHeight,self.deltaHeight,self.samplesWin)
732 samplingWindow = numpy.array(sampleWindowTuple,SAMPLING_STRUCTURE)
751 sampleWindowTuple = (
752 self.firstHeight, self.deltaHeight, self.samplesWin)
753 samplingWindow = numpy.array(sampleWindowTuple, SAMPLING_STRUCTURE)
733 754 samplingWindow.tofile(fp)
734
755
735 756 if self.totalSpectra != 0:
736 # spectraComb = numpy.array([],numpy.dtype('u1'))
757 # spectraComb = numpy.array([],numpy.dtype('u1'))
737 758 spectraComb = self.spectraComb
738 759 spectraComb.tofile(fp)
739
760
740 761 # if self.processFlags & PROCFLAG.DEFINE_PROCESS_CODE == PROCFLAG.DEFINE_PROCESS_CODE:
741 762 # nCode = numpy.array([self.nCode], numpy.dtype('u4')) #Probar con un dato que almacene codigo, hasta el momento no se hizo la prueba
742 763 # nCode.tofile(fp)
743 #
764 #
744 765 # nBaud = numpy.array([self.nBaud], numpy.dtype('u4'))
745 766 # nBaud.tofile(fp)
746 #
767 #
747 768 # code = self.code.reshape(self.nCode*self.nBaud)
748 769 # code = code.astype(numpy.dtype('<f4'))
749 770 # code.tofile(fp)
750
771
751 772 return 1
752 773
753 774 def get_size(self):
754
755 self.__size = 40 + 12*self.nWindows + 2*self.totalSpectra
756
775
776 self.__size = 40 + 12 * self.nWindows + 2 * self.totalSpectra
777
757 778 # if self.processFlags & PROCFLAG.DEFINE_PROCESS_CODE == PROCFLAG.DEFINE_PROCESS_CODE:
758 779 # self.__size += 4 + 4 + 4*self.nCode*numpy.ceil(self.nBaud/32.)
759 780 # self.__size += 4 + 4 + 4 * self.nCode * self.nBaud
760
781
761 782 return self.__size
762
783
763 784 def set_size(self, value):
764
785
765 786 raise IOError, "size is a property and it cannot be set, just read"
766
787
767 788 return
768
789
769 790 size = property(get_size, set_size)
770
791
792
771 793 class RCfunction:
772 NONE=0
773 FLIP=1
774 CODE=2
775 SAMPLING=3
776 LIN6DIV256=4
777 SYNCHRO=5
794 NONE = 0
795 FLIP = 1
796 CODE = 2
797 SAMPLING = 3
798 LIN6DIV256 = 4
799 SYNCHRO = 5
800
778 801
779 802 class nCodeType:
780 NONE=0
781 USERDEFINE=1
782 BARKER2=2
783 BARKER3=3
784 BARKER4=4
785 BARKER5=5
786 BARKER7=6
787 BARKER11=7
788 BARKER13=8
789 AC128=9
790 COMPLEMENTARYCODE2=10
791 COMPLEMENTARYCODE4=11
792 COMPLEMENTARYCODE8=12
793 COMPLEMENTARYCODE16=13
794 COMPLEMENTARYCODE32=14
795 COMPLEMENTARYCODE64=15
796 COMPLEMENTARYCODE128=16
797 CODE_BINARY28=17
798
799 class PROCFLAG:
800
803 NONE = 0
804 USERDEFINE = 1
805 BARKER2 = 2
806 BARKER3 = 3
807 BARKER4 = 4
808 BARKER5 = 5
809 BARKER7 = 6
810 BARKER11 = 7
811 BARKER13 = 8
812 AC128 = 9
813 COMPLEMENTARYCODE2 = 10
814 COMPLEMENTARYCODE4 = 11
815 COMPLEMENTARYCODE8 = 12
816 COMPLEMENTARYCODE16 = 13
817 COMPLEMENTARYCODE32 = 14
818 COMPLEMENTARYCODE64 = 15
819 COMPLEMENTARYCODE128 = 16
820 CODE_BINARY28 = 17
821
822
823 class PROCFLAG:
824
801 825 COHERENT_INTEGRATION = numpy.uint32(0x00000001)
802 DECODE_DATA = numpy.uint32(0x00000002)
826 DECODE_DATA = numpy.uint32(0x00000002)
803 827 SPECTRA_CALC = numpy.uint32(0x00000004)
804 INCOHERENT_INTEGRATION = numpy.uint32(0x00000008)
828 INCOHERENT_INTEGRATION = numpy.uint32(0x00000008)
805 829 POST_COHERENT_INTEGRATION = numpy.uint32(0x00000010)
806 830 SHIFT_FFT_DATA = numpy.uint32(0x00000020)
807
831
808 832 DATATYPE_CHAR = numpy.uint32(0x00000040)
809 833 DATATYPE_SHORT = numpy.uint32(0x00000080)
810 834 DATATYPE_LONG = numpy.uint32(0x00000100)
811 835 DATATYPE_INT64 = numpy.uint32(0x00000200)
812 836 DATATYPE_FLOAT = numpy.uint32(0x00000400)
813 837 DATATYPE_DOUBLE = numpy.uint32(0x00000800)
814
815 DATAARRANGE_CONTIGUOUS_CH = numpy.uint32(0x00001000)
816 DATAARRANGE_CONTIGUOUS_H = numpy.uint32(0x00002000)
817 DATAARRANGE_CONTIGUOUS_P = numpy.uint32(0x00004000)
818
838
839 DATAARRANGE_CONTIGUOUS_CH = numpy.uint32(0x00001000)
840 DATAARRANGE_CONTIGUOUS_H = numpy.uint32(0x00002000)
841 DATAARRANGE_CONTIGUOUS_P = numpy.uint32(0x00004000)
842
819 843 SAVE_CHANNELS_DC = numpy.uint32(0x00008000)
820 DEFLIP_DATA = numpy.uint32(0x00010000)
821 DEFINE_PROCESS_CODE = numpy.uint32(0x00020000)
822
844 DEFLIP_DATA = numpy.uint32(0x00010000)
845 DEFINE_PROCESS_CODE = numpy.uint32(0x00020000)
846
823 847 ACQ_SYS_NATALIA = numpy.uint32(0x00040000)
824 848 ACQ_SYS_ECHOTEK = numpy.uint32(0x00080000)
825 849 ACQ_SYS_ADRXD = numpy.uint32(0x000C0000)
826 850 ACQ_SYS_JULIA = numpy.uint32(0x00100000)
827 851 ACQ_SYS_XXXXXX = numpy.uint32(0x00140000)
828
852
829 853 EXP_NAME_ESP = numpy.uint32(0x00200000)
830 854 CHANNEL_NAMES_ESP = numpy.uint32(0x00400000)
831
855
832 856 OPERATION_MASK = numpy.uint32(0x0000003F)
833 857 DATATYPE_MASK = numpy.uint32(0x00000FC0)
834 858 DATAARRANGE_MASK = numpy.uint32(0x00007000)
835 859 ACQ_SYS_MASK = numpy.uint32(0x001C0000)
836
837 dtype0 = numpy.dtype([('real','<i1'),('imag','<i1')])
838 dtype1 = numpy.dtype([('real','<i2'),('imag','<i2')])
839 dtype2 = numpy.dtype([('real','<i4'),('imag','<i4')])
840 dtype3 = numpy.dtype([('real','<i8'),('imag','<i8')])
841 dtype4 = numpy.dtype([('real','<f4'),('imag','<f4')])
842 dtype5 = numpy.dtype([('real','<f8'),('imag','<f8')])
860
861
862 dtype0 = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
863 dtype1 = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
864 dtype2 = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
865 dtype3 = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
866 dtype4 = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
867 dtype5 = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
843 868
844 869 NUMPY_DTYPE_LIST = [dtype0, dtype1, dtype2, dtype3, dtype4, dtype5]
845 870
846 PROCFLAG_DTYPE_LIST = [PROCFLAG.DATATYPE_CHAR,
847 PROCFLAG.DATATYPE_SHORT,
848 PROCFLAG.DATATYPE_LONG,
849 PROCFLAG.DATATYPE_INT64,
850 PROCFLAG.DATATYPE_FLOAT,
871 PROCFLAG_DTYPE_LIST = [PROCFLAG.DATATYPE_CHAR,
872 PROCFLAG.DATATYPE_SHORT,
873 PROCFLAG.DATATYPE_LONG,
874 PROCFLAG.DATATYPE_INT64,
875 PROCFLAG.DATATYPE_FLOAT,
851 876 PROCFLAG.DATATYPE_DOUBLE]
852 877
853 878 DTYPE_WIDTH = [1, 2, 4, 8, 4, 8]
854 879
880
855 881 def get_dtype_index(numpy_dtype):
856
882
857 883 index = None
858
884
859 885 for i in range(len(NUMPY_DTYPE_LIST)):
860 886 if numpy_dtype == NUMPY_DTYPE_LIST[i]:
861 887 index = i
862 888 break
863
889
864 890 return index
865 891
892
866 893 def get_numpy_dtype(index):
867
894
868 895 return NUMPY_DTYPE_LIST[index]
869 896
897
870 898 def get_procflag_dtype(index):
871
899
872 900 return PROCFLAG_DTYPE_LIST[index]
873 901
902
874 903 def get_dtype_width(index):
875
876 return DTYPE_WIDTH[index] No newline at end of file
904
905 return DTYPE_WIDTH[index]
Show file before | Show file after | Hide comments
@@ -657,11 +657,11 class DigitalRFWriter(Operation):
657 657 self.metadata_dict['flagDecodeData'] = self.dataOut.flagDecodeData
658 658 self.metadata_dict['flagDeflipData'] = self.dataOut.flagDeflipData
659 659 self.metadata_dict['flagShiftFFT'] = self.dataOut.flagShiftFFT
660 self.metadata_dict['flagDataAsBlock'] = self.dataOut.flagDataAsBlock
661 660 self.metadata_dict['useLocalTime'] = self.dataOut.useLocalTime
662 661 self.metadata_dict['nCohInt'] = self.dataOut.nCohInt
663
664 return
662 self.metadata_dict['type'] = self.dataOut.type
663 self.metadata_dict['flagDataAsBlock'] = getattr(
664 self.dataOut, 'flagDataAsBlock', None) # chequear
665 665
666 666 def setup(self, dataOut, path, frequency, fileCadence, dirCadence, metadataCadence, set=0, metadataFile='metadata', ext='.h5'):
667 667 '''
@@ -678,10 +678,14 class DigitalRFWriter(Operation):
678 678 self.__dtype = dataOut.dtype[0]
679 679 self.__nSamples = dataOut.systemHeaderObj.nSamples
680 680 self.__nProfiles = dataOut.nProfiles
681 self.__blocks_per_file = dataOut.processingHeaderObj.dataBlocksPerFile
682 681
683 self.arr_data = arr_data = numpy.ones((self.__nSamples, len(
684 self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
682 if self.dataOut.type != 'Voltage':
683 raise 'Digital RF cannot be used with this data type'
684 self.arr_data = numpy.ones((1, dataOut.nFFTPoints * len(
685 self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
686 else:
687 self.arr_data = numpy.ones((self.__nSamples, len(
688 self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
685 689
686 690 file_cadence_millisecs = 1000
687 691
@@ -702,14 +706,11 class DigitalRFWriter(Operation):
702 706 fileCadence, start_global_index,
703 707 sample_rate_numerator, sample_rate_denominator, uuid, compression_level, checksum,
704 708 is_complex, num_subchannels, is_continuous, marching_periods)
705
706 709 metadata_dir = os.path.join(path, 'metadata')
707 710 os.system('mkdir %s' % (metadata_dir))
708
709 711 self.digitalMetadataWriteObj = digital_rf.DigitalMetadataWriter(metadata_dir, dirCadence, 1, # 236, file_cadence_millisecs / 1000
710 712 sample_rate_numerator, sample_rate_denominator,
711 713 metadataFile)
712
713 714 self.isConfig = True
714 715 self.currentSample = 0
715 716 self.oldAverage = 0
@@ -717,7 +718,6 class DigitalRFWriter(Operation):
717 718 return
718 719
719 720 def writeMetadata(self):
720 print '[Writing] - Writing metadata'
721 721 start_idx = self.__sample_rate * self.dataOut.utctime
722 722
723 723 self.metadata_dict['processingHeader'] = self.dataOut.processingHeaderObj.getAsDict(
@@ -741,10 +741,19 class DigitalRFWriter(Operation):
741 741 return
742 742
743 743 def writeData(self):
744 for i in range(self.dataOut.systemHeaderObj.nSamples):
744 if self.dataOut.type != 'Voltage':
745 raise 'Digital RF cannot be used with this data type'
745 746 for channel in self.dataOut.channelList:
746 self.arr_data[i][channel]['r'] = self.dataOut.data[channel][i].real
747 self.arr_data[i][channel]['i'] = self.dataOut.data[channel][i].imag
747 for i in range(self.dataOut.nFFTPoints):
748 self.arr_data[1][channel * self.dataOut.nFFTPoints +
749 i]['r'] = self.dataOut.data[channel][i].real
750 self.arr_data[1][channel * self.dataOut.nFFTPoints +
751 i]['i'] = self.dataOut.data[channel][i].imag
752 else:
753 for i in range(self.dataOut.systemHeaderObj.nSamples):
754 for channel in self.dataOut.channelList:
755 self.arr_data[i][channel]['r'] = self.dataOut.data[channel][i].real
756 self.arr_data[i][channel]['i'] = self.dataOut.data[channel][i].imag
748 757
749 758 def f(): return self.digitalWriteObj.rf_write(self.arr_data)
750 759 self.timeit(f)
Show file before | Show file after | Hide comments
@@ -6,6 +6,7 from jroproc_base import ProcessingUnit, Operation
6 6 from schainpy.model.data.jrodata import Spectra
7 7 from schainpy.model.data.jrodata import hildebrand_sekhon
8 8
9
9 10 class SpectraProc(ProcessingUnit):
10 11
11 12 def __init__(self, **kwargs):
@@ -25,12 +26,15 class SpectraProc(ProcessingUnit):
25 26 self.dataOut.dstFlag = self.dataIn.dstFlag
26 27 self.dataOut.errorCount = self.dataIn.errorCount
27 28 self.dataOut.useLocalTime = self.dataIn.useLocalTime
28
29 try:
30 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
31 except:
32 pass
29 33 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
30 34 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
31 35 self.dataOut.channelList = self.dataIn.channelList
32 36 self.dataOut.heightList = self.dataIn.heightList
33 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
37 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
34 38
35 39 self.dataOut.nBaud = self.dataIn.nBaud
36 40 self.dataOut.nCode = self.dataIn.nCode
@@ -39,8 +43,10 class SpectraProc(ProcessingUnit):
39 43
40 44 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
41 45 self.dataOut.utctime = self.firstdatatime
42 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
43 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
46 # asumo q la data esta decodificada
47 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
48 # asumo q la data esta sin flip
49 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
44 50 self.dataOut.flagShiftFFT = False
45 51
46 52 self.dataOut.nCohInt = self.dataIn.nCohInt
@@ -71,12 +77,13 class SpectraProc(ProcessingUnit):
71 77 self.buffer
72 78 self.dataOut.flagNoData
73 79 """
74 fft_volt = numpy.fft.fft(self.buffer,n=self.dataOut.nFFTPoints,axis=1)
80 fft_volt = numpy.fft.fft(
81 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
75 82 fft_volt = fft_volt.astype(numpy.dtype('complex'))
76 dc = fft_volt[:,0,:]
83 dc = fft_volt[:, 0, :]
77 84
78 #calculo de self-spectra
79 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
85 # calculo de self-spectra
86 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
80 87 spc = fft_volt * numpy.conjugate(fft_volt)
81 88 spc = spc.real
82 89
@@ -87,15 +94,19 class SpectraProc(ProcessingUnit):
87 94 cspc = None
88 95 pairIndex = 0
89 96 if self.dataOut.pairsList != None:
90 #calculo de cross-spectra
91 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
97 # calculo de cross-spectra
98 cspc = numpy.zeros(
99 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
92 100 for pair in self.dataOut.pairsList:
93 101 if pair[0] not in self.dataOut.channelList:
94 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
102 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
103 str(pair), str(self.dataOut.channelList))
95 104 if pair[1] not in self.dataOut.channelList:
96 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
105 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
106 str(pair), str(self.dataOut.channelList))
97 107
98 cspc[pairIndex,:,:] = fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:])
108 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
109 numpy.conjugate(fft_volt[pair[1], :, :])
99 110 pairIndex += 1
100 111 blocksize += cspc.size
101 112
@@ -112,7 +123,7 class SpectraProc(ProcessingUnit):
112 123 if self.dataIn.type == "Spectra":
113 124 self.dataOut.copy(self.dataIn)
114 125 if not pairsList:
115 pairsList = itertools.combinations(self.dataOut.channelList, 2)
126 pairsList = itertools.combinations(self.dataOut.channelList, 2)
116 127 if self.dataOut.data_cspc is not None:
117 128 self.__selectPairs(pairsList)
118 129 return True
@@ -134,13 +145,13 class SpectraProc(ProcessingUnit):
134 145 self.dataOut.pairsList = pairsList
135 146
136 147 if self.buffer is None:
137 self.buffer = numpy.zeros( (self.dataIn.nChannels,
138 nProfiles,
139 self.dataIn.nHeights),
148 self.buffer = numpy.zeros((self.dataIn.nChannels,
149 nProfiles,
150 self.dataIn.nHeights),
140 151 dtype='complex')
141 152
142 153 if self.dataIn.flagDataAsBlock:
143 #data dimension: [nChannels, nProfiles, nSamples]
154 # data dimension: [nChannels, nProfiles, nSamples]
144 155 nVoltProfiles = self.dataIn.data.shape[1]
145 156 # nVoltProfiles = self.dataIn.nProfiles
146 157
@@ -154,16 +165,18 class SpectraProc(ProcessingUnit):
154 165 self.id_min = 0
155 166 self.id_max = nVoltProfiles
156 167
157 self.buffer[:,self.id_min:self.id_max,:] = self.dataIn.data
168 self.buffer[:, self.id_min:self.id_max,
169 :] = self.dataIn.data
158 170 self.profIndex += nVoltProfiles
159 171 self.id_min += nVoltProfiles
160 172 self.id_max += nVoltProfiles
161 173 else:
162 raise ValueError, "The type object %s has %d profiles, it should just has %d profiles"%(self.dataIn.type,self.dataIn.data.shape[1],nProfiles)
174 raise ValueError, "The type object %s has %d profiles, it should just has %d profiles" % (
175 self.dataIn.type, self.dataIn.data.shape[1], nProfiles)
163 176 self.dataOut.flagNoData = True
164 177 return 0
165 178 else:
166 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
179 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
167 180 self.profIndex += 1
168 181
169 182 if self.firstdatatime == None:
@@ -179,11 +192,12 class SpectraProc(ProcessingUnit):
179 192
180 193 return True
181 194
182 raise ValueError, "The type of input object '%s' is not valid"%(self.dataIn.type)
195 raise ValueError, "The type of input object '%s' is not valid" % (
196 self.dataIn.type)
183 197
184 198 def __selectPairs(self, pairsList):
185 199
186 if not pairsList:
200 if not pairsList:
187 201 return
188 202
189 203 pairs = []
@@ -194,7 +208,7 class SpectraProc(ProcessingUnit):
194 208 continue
195 209 pairs.append(pair)
196 210 pairsIndex.append(pairs.index(pair))
197
211
198 212 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
199 213 self.dataOut.pairsList = pairs
200 214
@@ -207,10 +221,10 class SpectraProc(ProcessingUnit):
207 221
208 222 pairsIndexListSelected = []
209 223 for pairIndex in self.dataOut.pairsIndexList:
210 #First pair
224 # First pair
211 225 if self.dataOut.pairsList[pairIndex][0] not in channelList:
212 226 continue
213 #Second pair
227 # Second pair
214 228 if self.dataOut.pairsList[pairIndex][1] not in channelList:
215 229 continue
216 230
@@ -222,7 +236,8 class SpectraProc(ProcessingUnit):
222 236 return
223 237
224 238 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
225 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
239 self.dataOut.pairsList = [self.dataOut.pairsList[i]
240 for i in pairsIndexListSelected]
226 241
227 242 return
228 243
@@ -232,7 +247,8 class SpectraProc(ProcessingUnit):
232 247
233 248 for channel in channelList:
234 249 if channel not in self.dataOut.channelList:
235 raise ValueError, "Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList))
250 raise ValueError, "Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" % (
251 channel, str(self.dataOut.channelList))
236 252
237 253 index = self.dataOut.channelList.index(channel)
238 254 channelIndexList.append(index)
@@ -257,17 +273,19 class SpectraProc(ProcessingUnit):
257 273
258 274 for channelIndex in channelIndexList:
259 275 if channelIndex not in self.dataOut.channelIndexList:
260 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList)
276 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " % (
277 channelIndex, self.dataOut.channelIndexList)
261 278
262 279 # nChannels = len(channelIndexList)
263 280
264 data_spc = self.dataOut.data_spc[channelIndexList,:]
265 data_dc = self.dataOut.data_dc[channelIndexList,:]
281 data_spc = self.dataOut.data_spc[channelIndexList, :]
282 data_dc = self.dataOut.data_dc[channelIndexList, :]
266 283
267 284 self.dataOut.data_spc = data_spc
268 285 self.dataOut.data_dc = data_dc
269 286
270 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
287 self.dataOut.channelList = [
288 self.dataOut.channelList[i] for i in channelIndexList]
271 289 # self.dataOut.nChannels = nChannels
272 290
273 291 self.__selectPairsByChannel(self.dataOut.channelList)
@@ -291,7 +309,8 class SpectraProc(ProcessingUnit):
291 309 """
292 310
293 311 if (minHei > maxHei):
294 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
312 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (
313 minHei, maxHei)
295 314
296 315 if (minHei < self.dataOut.heightList[0]):
297 316 minHei = self.dataOut.heightList[0]
@@ -320,20 +339,23 class SpectraProc(ProcessingUnit):
320 339
321 340 return 1
322 341
323 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
324 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
342 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
343 newheis = numpy.where(
344 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
325 345
326 346 if hei_ref != None:
327 newheis = numpy.where(self.dataOut.heightList>hei_ref)
347 newheis = numpy.where(self.dataOut.heightList > hei_ref)
328 348
329 349 minIndex = min(newheis[0])
330 350 maxIndex = max(newheis[0])
331 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
332 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
351 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
352 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
333 353
334 354 # determina indices
335 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
336 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
355 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
356 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
357 avg_dB = 10 * \
358 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
337 359 beacon_dB = numpy.sort(avg_dB)[-nheis:]
338 360 beacon_heiIndexList = []
339 361 for val in avg_dB.tolist():
@@ -343,12 +365,12 class SpectraProc(ProcessingUnit):
343 365 #data_spc = data_spc[:,:,beacon_heiIndexList]
344 366 data_cspc = None
345 367 if self.dataOut.data_cspc is not None:
346 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
368 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
347 369 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
348 370
349 371 data_dc = None
350 372 if self.dataOut.data_dc is not None:
351 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
373 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
352 374 #data_dc = data_dc[:,beacon_heiIndexList]
353 375
354 376 self.dataOut.data_spc = data_spc
@@ -359,7 +381,6 class SpectraProc(ProcessingUnit):
359 381
360 382 return 1
361 383
362
363 384 def selectHeightsByIndex(self, minIndex, maxIndex):
364 385 """
365 386 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
@@ -380,104 +401,107 class SpectraProc(ProcessingUnit):
380 401 """
381 402
382 403 if (minIndex < 0) or (minIndex > maxIndex):
383 raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex)
404 raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (
405 minIndex, maxIndex)
384 406
385 407 if (maxIndex >= self.dataOut.nHeights):
386 maxIndex = self.dataOut.nHeights-1
408 maxIndex = self.dataOut.nHeights - 1
387 409
388 #Spectra
389 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
410 # Spectra
411 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
390 412
391 413 data_cspc = None
392 414 if self.dataOut.data_cspc is not None:
393 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
415 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
394 416
395 417 data_dc = None
396 418 if self.dataOut.data_dc is not None:
397 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
419 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
398 420
399 421 self.dataOut.data_spc = data_spc
400 422 self.dataOut.data_cspc = data_cspc
401 423 self.dataOut.data_dc = data_dc
402 424
403 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
425 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
404 426
405 427 return 1
406 428
407 def removeDC(self, mode = 2):
429 def removeDC(self, mode=2):
408 430 jspectra = self.dataOut.data_spc
409 431 jcspectra = self.dataOut.data_cspc
410 432
411
412 433 num_chan = jspectra.shape[0]
413 434 num_hei = jspectra.shape[2]
414 435
415 436 if jcspectra is not None:
416 437 jcspectraExist = True
417 438 num_pairs = jcspectra.shape[0]
418 else: jcspectraExist = False
439 else:
440 jcspectraExist = False
419 441
420 freq_dc = jspectra.shape[1]/2
421 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
442 freq_dc = jspectra.shape[1] / 2
443 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
422 444
423 if ind_vel[0]<0:
424 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
445 if ind_vel[0] < 0:
446 ind_vel[range(0, 1)] = ind_vel[range(0, 1)] + self.num_prof
425 447
426 448 if mode == 1:
427 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
449 jspectra[:, freq_dc, :] = (
450 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
428 451
429 452 if jcspectraExist:
430 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
453 jcspectra[:, freq_dc, :] = (
454 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
431 455
432 456 if mode == 2:
433 457
434 vel = numpy.array([-2,-1,1,2])
435 xx = numpy.zeros([4,4])
458 vel = numpy.array([-2, -1, 1, 2])
459 xx = numpy.zeros([4, 4])
436 460
437 461 for fil in range(4):
438 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
462 xx[fil, :] = vel[fil]**numpy.asarray(range(4))
439 463
440 464 xx_inv = numpy.linalg.inv(xx)
441 xx_aux = xx_inv[0,:]
465 xx_aux = xx_inv[0, :]
442 466
443 467 for ich in range(num_chan):
444 yy = jspectra[ich,ind_vel,:]
445 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
468 yy = jspectra[ich, ind_vel, :]
469 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
446 470
447 junkid = jspectra[ich,freq_dc,:]<=0
471 junkid = jspectra[ich, freq_dc, :] <= 0
448 472 cjunkid = sum(junkid)
449 473
450 474 if cjunkid.any():
451 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
475 jspectra[ich, freq_dc, junkid.nonzero()] = (
476 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
452 477
453 478 if jcspectraExist:
454 479 for ip in range(num_pairs):
455 yy = jcspectra[ip,ind_vel,:]
456 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
457
480 yy = jcspectra[ip, ind_vel, :]
481 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
458 482
459 483 self.dataOut.data_spc = jspectra
460 484 self.dataOut.data_cspc = jcspectra
461 485
462 486 return 1
463 487
464 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
488 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
465 489
466 490 jspectra = self.dataOut.data_spc
467 491 jcspectra = self.dataOut.data_cspc
468 492 jnoise = self.dataOut.getNoise()
469 493 num_incoh = self.dataOut.nIncohInt
470 494
471 num_channel = jspectra.shape[0]
472 num_prof = jspectra.shape[1]
473 num_hei = jspectra.shape[2]
495 num_channel = jspectra.shape[0]
496 num_prof = jspectra.shape[1]
497 num_hei = jspectra.shape[2]
474 498
475 #hei_interf
499 # hei_interf
476 500 if hei_interf is None:
477 count_hei = num_hei/2 #Como es entero no importa
501 count_hei = num_hei / 2 # Como es entero no importa
478 502 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
479 503 hei_interf = numpy.asarray(hei_interf)[0]
480 #nhei_interf
504 # nhei_interf
481 505 if (nhei_interf == None):
482 506 nhei_interf = 5
483 507 if (nhei_interf < 1):
@@ -492,136 +516,153 class SpectraProc(ProcessingUnit):
492 516 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
493 517 mask_prof = numpy.asarray(range(num_prof))
494 518 num_mask_prof = mask_prof.size
495 comp_mask_prof = [0, num_prof/2]
496
519 comp_mask_prof = [0, num_prof / 2]
497 520
498 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
521 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
499 522 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
500 523 jnoise = numpy.nan
501 524 noise_exist = jnoise[0] < numpy.Inf
502 525
503 #Subrutina de Remocion de la Interferencia
526 # Subrutina de Remocion de la Interferencia
504 527 for ich in range(num_channel):
505 #Se ordena los espectros segun su potencia (menor a mayor)
506 power = jspectra[ich,mask_prof,:]
507 power = power[:,hei_interf]
508 power = power.sum(axis = 0)
528 # Se ordena los espectros segun su potencia (menor a mayor)
529 power = jspectra[ich, mask_prof, :]
530 power = power[:, hei_interf]
531 power = power.sum(axis=0)
509 532 psort = power.ravel().argsort()
510 533
511 #Se estima la interferencia promedio en los Espectros de Potencia empleando
512 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
534 # Se estima la interferencia promedio en los Espectros de Potencia empleando
535 junkspc_interf = jspectra[ich, :, hei_interf[psort[range(
536 offhei_interf, nhei_interf + offhei_interf)]]]
513 537
514 538 if noise_exist:
515 # tmp_noise = jnoise[ich] / num_prof
539 # tmp_noise = jnoise[ich] / num_prof
516 540 tmp_noise = jnoise[ich]
517 541 junkspc_interf = junkspc_interf - tmp_noise
518 542 #junkspc_interf[:,comp_mask_prof] = 0
519 543
520 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
544 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
521 545 jspc_interf = jspc_interf.transpose()
522 #Calculando el espectro de interferencia promedio
523 noiseid = numpy.where(jspc_interf <= tmp_noise/ numpy.sqrt(num_incoh))
546 # Calculando el espectro de interferencia promedio
547 noiseid = numpy.where(
548 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
524 549 noiseid = noiseid[0]
525 550 cnoiseid = noiseid.size
526 interfid = numpy.where(jspc_interf > tmp_noise/ numpy.sqrt(num_incoh))
551 interfid = numpy.where(
552 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
527 553 interfid = interfid[0]
528 554 cinterfid = interfid.size
529 555
530 if (cnoiseid > 0): jspc_interf[noiseid] = 0
556 if (cnoiseid > 0):
557 jspc_interf[noiseid] = 0
531 558
532 #Expandiendo los perfiles a limpiar
559 # Expandiendo los perfiles a limpiar
533 560 if (cinterfid > 0):
534 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
561 new_interfid = (
562 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
535 563 new_interfid = numpy.asarray(new_interfid)
536 564 new_interfid = {x for x in new_interfid}
537 565 new_interfid = numpy.array(list(new_interfid))
538 566 new_cinterfid = new_interfid.size
539 else: new_cinterfid = 0
567 else:
568 new_cinterfid = 0
540 569
541 570 for ip in range(new_cinterfid):
542 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
543 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
544
571 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
572 jspc_interf[new_interfid[ip]
573 ] = junkspc_interf[ind[nhei_interf / 2], new_interfid[ip]]
545 574
546 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
575 jspectra[ich, :, ind_hei] = jspectra[ich, :,
576 ind_hei] - jspc_interf # Corregir indices
547 577
548 #Removiendo la interferencia del punto de mayor interferencia
578 # Removiendo la interferencia del punto de mayor interferencia
549 579 ListAux = jspc_interf[mask_prof].tolist()
550 580 maxid = ListAux.index(max(ListAux))
551 581
552
553 582 if cinterfid > 0:
554 for ip in range(cinterfid*(interf == 2) - 1):
555 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/numpy.sqrt(num_incoh))).nonzero()
583 for ip in range(cinterfid * (interf == 2) - 1):
584 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
585 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
556 586 cind = len(ind)
557 587
558 588 if (cind > 0):
559 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/numpy.sqrt(num_incoh))
589 jspectra[ich, interfid[ip], ind] = tmp_noise * \
590 (1 + (numpy.random.uniform(cind) - 0.5) /
591 numpy.sqrt(num_incoh))
560 592
561 ind = numpy.array([-2,-1,1,2])
562 xx = numpy.zeros([4,4])
593 ind = numpy.array([-2, -1, 1, 2])
594 xx = numpy.zeros([4, 4])
563 595
564 596 for id1 in range(4):
565 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
597 xx[:, id1] = ind[id1]**numpy.asarray(range(4))
566 598
567 599 xx_inv = numpy.linalg.inv(xx)
568 xx = xx_inv[:,0]
569 ind = (ind + maxid + num_mask_prof)%num_mask_prof
570 yy = jspectra[ich,mask_prof[ind],:]
571 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
572
573
574 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/numpy.sqrt(num_incoh))).nonzero()
575 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/numpy.sqrt(num_incoh))
576
577 #Remocion de Interferencia en el Cross Spectra
578 if jcspectra is None: return jspectra, jcspectra
579 num_pairs = jcspectra.size/(num_prof*num_hei)
600 xx = xx_inv[:, 0]
601 ind = (ind + maxid + num_mask_prof) % num_mask_prof
602 yy = jspectra[ich, mask_prof[ind], :]
603 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
604 yy.transpose(), xx)
605
606 indAux = (jspectra[ich, :, :] < tmp_noise *
607 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
608 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
609 (1 - 1 / numpy.sqrt(num_incoh))
610
611 # Remocion de Interferencia en el Cross Spectra
612 if jcspectra is None:
613 return jspectra, jcspectra
614 num_pairs = jcspectra.size / (num_prof * num_hei)
580 615 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
581 616
582 617 for ip in range(num_pairs):
583 618
584 619 #-------------------------------------------
585 620
586 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
587 cspower = cspower[:,hei_interf]
588 cspower = cspower.sum(axis = 0)
621 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
622 cspower = cspower[:, hei_interf]
623 cspower = cspower.sum(axis=0)
589 624
590 625 cspsort = cspower.ravel().argsort()
591 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
626 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[range(
627 offhei_interf, nhei_interf + offhei_interf)]]]
592 628 junkcspc_interf = junkcspc_interf.transpose()
593 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
629 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
594 630
595 631 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
596 632
597 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
598 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
599 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
633 median_real = numpy.median(numpy.real(
634 junkcspc_interf[mask_prof[ind[range(3 * num_prof / 4)]], :]))
635 median_imag = numpy.median(numpy.imag(
636 junkcspc_interf[mask_prof[ind[range(3 * num_prof / 4)]], :]))
637 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
638 median_real, median_imag)
600 639
601 640 for iprof in range(num_prof):
602 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
603 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
641 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
642 jcspc_interf[iprof] = junkcspc_interf[iprof,
643 ind[nhei_interf / 2]]
604 644
605 #Removiendo la Interferencia
606 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
645 # Removiendo la Interferencia
646 jcspectra[ip, :, ind_hei] = jcspectra[ip,
647 :, ind_hei] - jcspc_interf
607 648
608 649 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
609 650 maxid = ListAux.index(max(ListAux))
610 651
611 ind = numpy.array([-2,-1,1,2])
612 xx = numpy.zeros([4,4])
652 ind = numpy.array([-2, -1, 1, 2])
653 xx = numpy.zeros([4, 4])
613 654
614 655 for id1 in range(4):
615 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
656 xx[:, id1] = ind[id1]**numpy.asarray(range(4))
616 657
617 658 xx_inv = numpy.linalg.inv(xx)
618 xx = xx_inv[:,0]
659 xx = xx_inv[:, 0]
619 660
620 ind = (ind + maxid + num_mask_prof)%num_mask_prof
621 yy = jcspectra[ip,mask_prof[ind],:]
622 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
661 ind = (ind + maxid + num_mask_prof) % num_mask_prof
662 yy = jcspectra[ip, mask_prof[ind], :]
663 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
623 664
624 #Guardar Resultados
665 # Guardar Resultados
625 666 self.dataOut.data_spc = jspectra
626 667 self.dataOut.data_cspc = jcspectra
627 668
@@ -635,7 +676,7 class SpectraProc(ProcessingUnit):
635 676 return 1
636 677
637 678 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
638 #validacion de rango
679 # validacion de rango
639 680 if minHei == None:
640 681 minHei = self.dataOut.heightList[0]
641 682
@@ -643,13 +684,13 class SpectraProc(ProcessingUnit):
643 684 maxHei = self.dataOut.heightList[-1]
644 685
645 686 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
646 print 'minHei: %.2f is out of the heights range'%(minHei)
647 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
687 print 'minHei: %.2f is out of the heights range' % (minHei)
688 print 'minHei is setting to %.2f' % (self.dataOut.heightList[0])
648 689 minHei = self.dataOut.heightList[0]
649 690
650 691 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
651 print 'maxHei: %.2f is out of the heights range'%(maxHei)
652 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
692 print 'maxHei: %.2f is out of the heights range' % (maxHei)
693 print 'maxHei is setting to %.2f' % (self.dataOut.heightList[-1])
653 694 maxHei = self.dataOut.heightList[-1]
654 695
655 696 # validacion de velocidades
@@ -662,13 +703,13 class SpectraProc(ProcessingUnit):
662 703 maxVel = velrange[-1]
663 704
664 705 if (minVel < velrange[0]) or (minVel > maxVel):
665 print 'minVel: %.2f is out of the velocity range'%(minVel)
666 print 'minVel is setting to %.2f'%(velrange[0])
706 print 'minVel: %.2f is out of the velocity range' % (minVel)
707 print 'minVel is setting to %.2f' % (velrange[0])
667 708 minVel = velrange[0]
668 709
669 710 if (maxVel > velrange[-1]) or (maxVel < minVel):
670 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
671 print 'maxVel is setting to %.2f'%(velrange[-1])
711 print 'maxVel: %.2f is out of the velocity range' % (maxVel)
712 print 'maxVel is setting to %.2f' % (velrange[-1])
672 713 maxVel = velrange[-1]
673 714
674 715 # seleccion de indices para rango
@@ -690,10 +731,11 class SpectraProc(ProcessingUnit):
690 731 maxIndex = len(heights)
691 732
692 733 if (minIndex < 0) or (minIndex > maxIndex):
693 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
734 raise ValueError, "some value in (%d,%d) is not valid" % (
735 minIndex, maxIndex)
694 736
695 737 if (maxIndex >= self.dataOut.nHeights):
696 maxIndex = self.dataOut.nHeights-1
738 maxIndex = self.dataOut.nHeights - 1
697 739
698 740 # seleccion de indices para velocidades
699 741 indminvel = numpy.where(velrange >= minVel)
@@ -708,24 +750,25 class SpectraProc(ProcessingUnit):
708 750 except:
709 751 maxIndexVel = len(velrange)
710 752
711 #seleccion del espectro
712 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
713 #estimacion de ruido
753 # seleccion del espectro
754 data_spc = self.dataOut.data_spc[:,
755 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
756 # estimacion de ruido
714 757 noise = numpy.zeros(self.dataOut.nChannels)
715 758
716 759 for channel in range(self.dataOut.nChannels):
717 daux = data_spc[channel,:,:]
760 daux = data_spc[channel, :, :]
718 761 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
719 762
720 763 self.dataOut.noise_estimation = noise.copy()
721 764
722 765 return 1
723 766
724 class IncohInt(Operation):
725 767
768 class IncohInt(Operation):
726 769
727 770 __profIndex = 0
728 __withOverapping = False
771 __withOverapping = False
729 772
730 773 __byTime = False
731 774 __initime = None
@@ -742,8 +785,6 class IncohInt(Operation):
742 785
743 786 n = None
744 787
745
746
747 788 def __init__(self, **kwargs):
748 789
749 790 Operation.__init__(self, **kwargs)
@@ -778,12 +819,12 class IncohInt(Operation):
778 819 if n is not None:
779 820 self.n = int(n)
780 821 else:
781 self.__integrationtime = int(timeInterval) #if (type(timeInterval)!=integer) -> change this line
822 # if (type(timeInterval)!=integer) -> change this line
823 self.__integrationtime = int(timeInterval)
782 824 self.n = None
783 825 self.__byTime = True
784 826
785 827 def putData(self, data_spc, data_cspc, data_dc):
786
787 828 """
788 829 Add a profile to the __buffer_spc and increase in one the __profileIndex
789 830
@@ -866,7 +907,8 class IncohInt(Operation):
866 907 self.__initime = datatime
867 908
868 909 if self.__byTime:
869 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
910 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
911 datatime, *args)
870 912 else:
871 913 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
872 914
@@ -876,7 +918,7 class IncohInt(Operation):
876 918 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
877 919
878 920 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
879 if n==1:
921 if n == 1:
880 922 return
881 923
882 924 dataOut.flagNoData = True
General Comments 0
You need to be logged in to leave comments. Login now