##// END OF EJS Templates
schain-jc
RSS
formatting, template actualizado, decimation a 300
José Chávez -
r1092:57cb558a168b
parent child
Show More
Side by Side Unified
Context file:
r1092:57cb558a168b Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,90 +1,90
1 1 basic = '''from schainpy.controller import Project
2 2
3 3 desc = "{desc}"
4 4 project = Project()
5 5 project.setup(id='200', name="{name}", description=desc)
6 6
7 7 voltage_reader = project.addReadUnit(datatype='VoltageReader',
8 8 path="{path}",
9 9 startDate="{startDate}",
10 10 endDate="{endDate}",
11 11 startTime="{startHour}",
12 12 endTime="{endHour}",
13 13 online=0,
14 14 verbose=1,
15 15 walk=1,
16 16 )
17 17
18 18 voltage_proc = project.addProcUnit(datatype='VoltageProc', inputId=voltage_reader.getId())
19 19
20 20 profile = voltage_proc.addOperation(name='ProfileSelector', optype='other')
21 21 profile.addParameter(name='profileRangeList', value='120,183', format='intlist')
22 22
23 23 rti = voltage_proc.addOperation(name='RTIPlot', optype='other')
24 24 rti.addParameter(name='wintitle', value='Jicamarca Radio Observatory', format='str')
25 25 rti.addParameter(name='showprofile', value='0', format='int')
26 26 rti.addParameter(name='xmin', value='0', format='int')
27 27 rti.addParameter(name='xmax', value='24', format='int')
28 28 rti.addParameter(name='figpath', value="{figpath}", format='str')
29 29 rti.addParameter(name='wr_period', value='5', format='int')
30 30 rti.addParameter(name='exp_code', value='22', format='int')
31 31
32 32
33 controller.start()
33 project.start()
34 34 '''
35 35
36 36 multiprocess = '''from schainpy.controller import Project, MPProject
37 37 from time import sleep
38 38 desc = "{desc}"
39 39
40 40 ####################
41 41 # PLOTTER RECEIVER #
42 42 ####################
43 43 plotter = Project()
44 44 plotter.setup(id='100', name='receiver', description=desc)
45 45
46 46 receiver_plot = plotter.addProcUnit(name='PlotterReceiver')
47 47 receiver_plot.addParameter(name='throttle', value=20, format='int')
48 48 receiver_plot.addParameter(name='plottypes', value='rti', format='str')
49 49
50 50 rti = receiver_plot.addOperation(name='PlotRTIData', optype='other')
51 51 rti.addParameter(name='zmin', value='-40.0', format='float')
52 52 rti.addParameter(name='zmax', value='100.0', format='float')
53 53 rti.addParameter(name='decimation', value='200', format='int')
54 54 rti.addParameter(name='xmin', value='0.0', format='int')
55 55 rti.addParameter(name='colormap', value='jet', format='str')
56 56
57 57 plotter.start()
58 58
59 59 sleep(2)
60 60
61 61 ################
62 62 # DATA EMITTER #
63 63 ################
64 64 project = Project()
65 65 project.setup(id='200', name="{name}", description=desc)
66 66
67 67 spectra_reader = project.addReadUnit(datatype='SpectraReader',
68 68 path="{path}",
69 69 startDate={startDate},
70 70 endDate={endDate},
71 71 startTime="{startHour}",
72 72 endTime="{endHour}",
73 73 online=0,
74 74 verbose=1,
75 75 walk=1,
76 76 )
77 77
78 78 spectra_proc = project.addProcUnit(datatype='Spectra', inputId=spectra_reader.getId())
79 79
80 80 parameters_proc = project.addProcUnit(datatype='ParametersProc', inputId=spectra_proc.getId())
81 81 moments = parameters_proc.addOperation(name='SpectralMoments', optype='other')
82 82
83 83 publish = parameters_proc.addOperation(name='PublishData', optype='other')
84 84 publish.addParameter(name='zeromq', value=1, format='int')
85 85 publish.addParameter(name='verbose', value=0, format='bool')
86 86
87 87 MPProject(project, 16)
88 88
89 89
90 90 '''
Show file before | Show file after | Hide comments
@@ -1,1227 +1,1251
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10
11 11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 12 from schainpy import cSchain
13 13
14 14
15 15 def getNumpyDtype(dataTypeCode):
16 16
17 17 if dataTypeCode == 0:
18 18 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
19 19 elif dataTypeCode == 1:
20 20 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
21 21 elif dataTypeCode == 2:
22 22 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
23 23 elif dataTypeCode == 3:
24 24 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
25 25 elif dataTypeCode == 4:
26 26 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
27 27 elif dataTypeCode == 5:
28 28 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
29 29 else:
30 30 raise ValueError, 'dataTypeCode was not defined'
31 31
32 32 return numpyDtype
33 33
34
34 35 def getDataTypeCode(numpyDtype):
35 36
36 37 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
37 38 datatype = 0
38 39 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
39 40 datatype = 1
40 41 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
41 42 datatype = 2
42 43 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
43 44 datatype = 3
44 45 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
45 46 datatype = 4
46 47 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
47 48 datatype = 5
48 49 else:
49 50 datatype = None
50 51
51 52 return datatype
52 53
54
53 55 def hildebrand_sekhon(data, navg):
54 56 """
55 57 This method is for the objective determination of the noise level in Doppler spectra. This
56 58 implementation technique is based on the fact that the standard deviation of the spectral
57 59 densities is equal to the mean spectral density for white Gaussian noise
58 60
59 61 Inputs:
60 62 Data : heights
61 63 navg : numbers of averages
62 64
63 65 Return:
64 66 -1 : any error
65 67 anoise : noise's level
66 68 """
67 69
68 70 sortdata = numpy.sort(data, axis=None)
69 71 # lenOfData = len(sortdata)
70 72 # nums_min = lenOfData*0.2
71 73 #
72 74 # if nums_min <= 5:
73 75 # nums_min = 5
74 76 #
75 77 # sump = 0.
76 78 #
77 79 # sumq = 0.
78 80 #
79 81 # j = 0
80 82 #
81 83 # cont = 1
82 84 #
83 85 # while((cont==1)and(j<lenOfData)):
84 86 #
85 87 # sump += sortdata[j]
86 88 #
87 89 # sumq += sortdata[j]**2
88 90 #
89 91 # if j > nums_min:
90 92 # rtest = float(j)/(j-1) + 1.0/navg
91 93 # if ((sumq*j) > (rtest*sump**2)):
92 94 # j = j - 1
93 95 # sump = sump - sortdata[j]
94 96 # sumq = sumq - sortdata[j]**2
95 97 # cont = 0
96 98 #
97 99 # j += 1
98 100 #
99 101 # lnoise = sump /j
100 102 #
101 103 # return lnoise
102 104
103 105 return cSchain.hildebrand_sekhon(sortdata, navg)
104 106
105 107
106 108 class Beam:
107 109
108 110 def __init__(self):
109 111 self.codeList = []
110 112 self.azimuthList = []
111 113 self.zenithList = []
112 114
115
113 116 class GenericData(object):
114 117
115 118 flagNoData = True
116 119
117 120 def copy(self, inputObj=None):
118 121
119 122 if inputObj == None:
120 123 return copy.deepcopy(self)
121 124
122 125 for key in inputObj.__dict__.keys():
123 126
124 127 attribute = inputObj.__dict__[key]
125 128
126 129 #If this attribute is a tuple or list
127 130 if type(inputObj.__dict__[key]) in (tuple, list):
128 131 self.__dict__[key] = attribute[:]
129 132 continue
130 133
131 134 #If this attribute is another object or instance
132 135 if hasattr(attribute, '__dict__'):
133 136 self.__dict__[key] = attribute.copy()
134 137 continue
135 138
136 139 self.__dict__[key] = inputObj.__dict__[key]
137 140
138 141 def deepcopy(self):
139 142
140 143 return copy.deepcopy(self)
141 144
142 145 def isEmpty(self):
143 146
144 147 return self.flagNoData
145 148
149
146 150 class JROData(GenericData):
147 151
148 152 # m_BasicHeader = BasicHeader()
149 153 # m_ProcessingHeader = ProcessingHeader()
150 154
151 155 systemHeaderObj = SystemHeader()
152 156
153 157 radarControllerHeaderObj = RadarControllerHeader()
154 158
155 159 # data = None
156 160
157 161 type = None
158 162
159 163 datatype = None #dtype but in string
160 164
161 165 # dtype = None
162 166
163 167 # nChannels = None
164 168
165 169 # nHeights = None
166 170
167 171 nProfiles = None
168 172
169 173 heightList = None
170 174
171 175 channelList = None
172 176
173 177 flagDiscontinuousBlock = False
174 178
175 179 useLocalTime = False
176 180
177 181 utctime = None
178 182
179 183 timeZone = None
180 184
181 185 dstFlag = None
182 186
183 187 errorCount = None
184 188
185 189 blocksize = None
186 190
187 191 # nCode = None
188 192 #
189 193 # nBaud = None
190 194 #
191 195 # code = None
192 196
193 197 flagDecodeData = False #asumo q la data no esta decodificada
194 198
195 199 flagDeflipData = False #asumo q la data no esta sin flip
196 200
197 201 flagShiftFFT = False
198 202
199 203 # ippSeconds = None
200 204
201 205 # timeInterval = None
202 206
203 207 nCohInt = None
204 208
205 209 # noise = None
206 210
207 211 windowOfFilter = 1
208 212
209 213 #Speed of ligth
210 214 C = 3e8
211 215
212 216 frequency = 49.92e6
213 217
214 218 realtime = False
215 219
216 220 beacon_heiIndexList = None
217 221
218 222 last_block = None
219 223
220 224 blocknow = None
221 225
222 226 azimuth = None
223 227
224 228 zenith = None
225 229
226 230 beam = Beam()
227 231
228 232 profileIndex = None
229 233
230 234 def getNoise(self):
231 235
232 236 raise NotImplementedError
233 237
234 238 def getNChannels(self):
235 239
236 240 return len(self.channelList)
237 241
238 242 def getChannelIndexList(self):
239 243
240 244 return range(self.nChannels)
241 245
242 246 def getNHeights(self):
243 247
244 248 return len(self.heightList)
245 249
246 250 def getHeiRange(self, extrapoints=0):
247 251
248 252 heis = self.heightList
249 253 # deltah = self.heightList[1] - self.heightList[0]
250 254 #
251 255 # heis.append(self.heightList[-1])
252 256
253 257 return heis
254 258
255 259 def getDeltaH(self):
256 260
257 261 delta = self.heightList[1] - self.heightList[0]
258 262
259 263 return delta
260 264
261 265 def getltctime(self):
262 266
263 267 if self.useLocalTime:
264 268 return self.utctime - self.timeZone*60
265 269
266 270 return self.utctime
267 271
268 272 def getDatatime(self):
269 273
270 274 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
271 275 return datatimeValue
272 276
273 277 def getTimeRange(self):
274 278
275 279 datatime = []
276 280
277 281 datatime.append(self.ltctime)
278 282 datatime.append(self.ltctime + self.timeInterval+1)
279 283
280 284 datatime = numpy.array(datatime)
281 285
282 286 return datatime
283 287
284 288 def getFmaxTimeResponse(self):
285 289
286 290 period = (10**-6)*self.getDeltaH()/(0.15)
287 291
288 292 PRF = 1./(period * self.nCohInt)
289 293
290 294 fmax = PRF
291 295
292 296 return fmax
293 297
294 298 def getFmax(self):
295 299 PRF = 1./(self.ippSeconds * self.nCohInt)
296 300
297 301 fmax = PRF
298 302 return fmax
299 303
300 304 def getVmax(self):
301 305
302 306 _lambda = self.C/self.frequency
303 307
304 308 vmax = self.getFmax() * _lambda/2
305 309
306 310 return vmax
307 311
308 312 def get_ippSeconds(self):
309 313 '''
310 314 '''
311 315 return self.radarControllerHeaderObj.ippSeconds
312 316
313 317 def set_ippSeconds(self, ippSeconds):
314 318 '''
315 319 '''
316 320
317 321 self.radarControllerHeaderObj.ippSeconds = ippSeconds
318 322
319 323 return
320 324
321 325 def get_dtype(self):
322 326 '''
323 327 '''
324 328 return getNumpyDtype(self.datatype)
325 329
326 330 def set_dtype(self, numpyDtype):
327 331 '''
328 332 '''
329 333
330 334 self.datatype = getDataTypeCode(numpyDtype)
331 335
332 336 def get_code(self):
333 337 '''
334 338 '''
335 339 return self.radarControllerHeaderObj.code
336 340
337 341 def set_code(self, code):
338 342 '''
339 343 '''
340 344 self.radarControllerHeaderObj.code = code
341 345
342 346 return
343 347
344 348 def get_ncode(self):
345 349 '''
346 350 '''
347 351 return self.radarControllerHeaderObj.nCode
348 352
349 353 def set_ncode(self, nCode):
350 354 '''
351 355 '''
352 356 self.radarControllerHeaderObj.nCode = nCode
353 357
354 358 return
355 359
356 360 def get_nbaud(self):
357 361 '''
358 362 '''
359 363 return self.radarControllerHeaderObj.nBaud
360 364
361 365 def set_nbaud(self, nBaud):
362 366 '''
363 367 '''
364 368 self.radarControllerHeaderObj.nBaud = nBaud
365 369
366 370 return
367 371
368 372 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
369 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
373 channelIndexList = property(
374 getChannelIndexList, "I'm the 'channelIndexList' property.")
370 375 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
371 376 #noise = property(getNoise, "I'm the 'nHeights' property.")
372 377 datatime = property(getDatatime, "I'm the 'datatime' property")
373 378 ltctime = property(getltctime, "I'm the 'ltctime' property")
374 379 ippSeconds = property(get_ippSeconds, set_ippSeconds)
375 380 dtype = property(get_dtype, set_dtype)
376 381 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
377 382 code = property(get_code, set_code)
378 383 nCode = property(get_ncode, set_ncode)
379 384 nBaud = property(get_nbaud, set_nbaud)
380 385
386
381 387 class Voltage(JROData):
382 388
383 389 #data es un numpy array de 2 dmensiones (canales, alturas)
384 390 data = None
385 391
386 392 def __init__(self):
387 393 '''
388 394 Constructor
389 395 '''
390 396
391 397 self.useLocalTime = True
392 398
393 399 self.radarControllerHeaderObj = RadarControllerHeader()
394 400
395 401 self.systemHeaderObj = SystemHeader()
396 402
397 403 self.type = "Voltage"
398 404
399 405 self.data = None
400 406
401 407 # self.dtype = None
402 408
403 409 # self.nChannels = 0
404 410
405 411 # self.nHeights = 0
406 412
407 413 self.nProfiles = None
408 414
409 415 self.heightList = None
410 416
411 417 self.channelList = None
412 418
413 419 # self.channelIndexList = None
414 420
415 421 self.flagNoData = True
416 422
417 423 self.flagDiscontinuousBlock = False
418 424
419 425 self.utctime = None
420 426
421 427 self.timeZone = None
422 428
423 429 self.dstFlag = None
424 430
425 431 self.errorCount = None
426 432
427 433 self.nCohInt = None
428 434
429 435 self.blocksize = None
430 436
431 437 self.flagDecodeData = False #asumo q la data no esta decodificada
432 438
433 439 self.flagDeflipData = False #asumo q la data no esta sin flip
434 440
435 441 self.flagShiftFFT = False
436 442
437 443 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
438 444
439 445 self.profileIndex = 0
440 446
441 447 def getNoisebyHildebrand(self, channel = None):
442 448 """
443 449 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
444 450
445 451 Return:
446 452 noiselevel
447 453 """
448 454
449 455 if channel != None:
450 456 data = self.data[channel]
451 457 nChannels = 1
452 458 else:
453 459 data = self.data
454 460 nChannels = self.nChannels
455 461
456 462 noise = numpy.zeros(nChannels)
457 463 power = data * numpy.conjugate(data)
458 464
459 465 for thisChannel in range(nChannels):
460 466 if nChannels == 1:
461 467 daux = power[:].real
462 468 else:
463 469 daux = power[thisChannel,:].real
464 470 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
465 471
466 472 return noise
467 473
468 474 def getNoise(self, type = 1, channel = None):
469 475
470 476 if type == 1:
471 477 noise = self.getNoisebyHildebrand(channel)
472 478
473 479 return noise
474 480
475 481 def getPower(self, channel = None):
476 482
477 483 if channel != None:
478 484 data = self.data[channel]
479 485 else:
480 486 data = self.data
481 487
482 488 power = data * numpy.conjugate(data)
483 489 powerdB = 10*numpy.log10(power.real)
484 490 powerdB = numpy.squeeze(powerdB)
485 491
486 492 return powerdB
487 493
488 494 def getTimeInterval(self):
489 495
490 496 timeInterval = self.ippSeconds * self.nCohInt
491 497
492 498 return timeInterval
493 499
494 500 noise = property(getNoise, "I'm the 'nHeights' property.")
495 501 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
496 502
503
497 504 class Spectra(JROData):
498 505
499 506 #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
500 507 data_spc = None
501 508
502 509 #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
503 510 data_cspc = None
504 511
505 512 #data dc es un numpy array de 2 dmensiones (canales, alturas)
506 513 data_dc = None
507 514
508 515 #data power
509 516 data_pwr = None
510 517
511 518 nFFTPoints = None
512 519
513 520 # nPairs = None
514 521
515 522 pairsList = None
516 523
517 524 nIncohInt = None
518 525
519 526 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
520 527
521 528 nCohInt = None #se requiere para determinar el valor de timeInterval
522 529
523 530 ippFactor = None
524 531
525 532 profileIndex = 0
526 533
527 534 plotting = "spectra"
528 535
529 536 def __init__(self):
530 537 '''
531 538 Constructor
532 539 '''
533 540
534 541 self.useLocalTime = True
535 542
536 543 self.radarControllerHeaderObj = RadarControllerHeader()
537 544
538 545 self.systemHeaderObj = SystemHeader()
539 546
540 547 self.type = "Spectra"
541 548
542 549 # self.data = None
543 550
544 551 # self.dtype = None
545 552
546 553 # self.nChannels = 0
547 554
548 555 # self.nHeights = 0
549 556
550 557 self.nProfiles = None
551 558
552 559 self.heightList = None
553 560
554 561 self.channelList = None
555 562
556 563 # self.channelIndexList = None
557 564
558 565 self.pairsList = None
559 566
560 567 self.flagNoData = True
561 568
562 569 self.flagDiscontinuousBlock = False
563 570
564 571 self.utctime = None
565 572
566 573 self.nCohInt = None
567 574
568 575 self.nIncohInt = None
569 576
570 577 self.blocksize = None
571 578
572 579 self.nFFTPoints = None
573 580
574 581 self.wavelength = None
575 582
576 583 self.flagDecodeData = False #asumo q la data no esta decodificada
577 584
578 585 self.flagDeflipData = False #asumo q la data no esta sin flip
579 586
580 587 self.flagShiftFFT = False
581 588
582 589 self.ippFactor = 1
583 590
584 591 #self.noise = None
585 592
586 593 self.beacon_heiIndexList = []
587 594
588 595 self.noise_estimation = None
589 596
590
591 597 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
592 598 """
593 599 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
594 600
595 601 Return:
596 602 noiselevel
597 603 """
598 604
599 605 noise = numpy.zeros(self.nChannels)
600 606
601 607 for channel in range(self.nChannels):
602 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
608 daux = self.data_spc[channel,
609 xmin_index:xmax_index, ymin_index:ymax_index]
603 610 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
604 611
605 612 return noise
606 613
607 614 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
608 615
609 616 if self.noise_estimation is not None:
610 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
617 # this was estimated by getNoise Operation defined in jroproc_spectra.py
618 return self.noise_estimation
611 619 else:
612 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
620 noise = self.getNoisebyHildebrand(
621 xmin_index, xmax_index, ymin_index, ymax_index)
613 622 return noise
614 623
615 624 def getFreqRangeTimeResponse(self, extrapoints=0):
616 625
617 626 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
618 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
627 freqrange = deltafreq * \
628 (numpy.arange(self.nFFTPoints + extrapoints) -
629 self.nFFTPoints / 2.) - deltafreq / 2
619 630
620 631 return freqrange
621 632
622 633 def getAcfRange(self, extrapoints=0):
623 634
624 635 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
625 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
636 freqrange = deltafreq * \
637 (numpy.arange(self.nFFTPoints + extrapoints) -
638 self.nFFTPoints / 2.) - deltafreq / 2
626 639
627 640 return freqrange
628 641
629 642 def getFreqRange(self, extrapoints=0):
630 643
631 644 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
632 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
645 freqrange = deltafreq * \
646 (numpy.arange(self.nFFTPoints + extrapoints) -
647 self.nFFTPoints / 2.) - deltafreq / 2
633 648
634 649 return freqrange
635 650
636 651 def getVelRange(self, extrapoints=0):
637 652
638 653 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
639 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2
654 velrange = deltav * (numpy.arange(self.nFFTPoints +
655 extrapoints) - self.nFFTPoints / 2.) # - deltav/2
640 656
641 657 return velrange
642 658
643 659 def getNPairs(self):
644 660
645 661 return len(self.pairsList)
646 662
647 663 def getPairsIndexList(self):
648 664
649 665 return range(self.nPairs)
650 666
651 667 def getNormFactor(self):
652 668
653 669 pwcode = 1
654 670
655 671 if self.flagDecodeData:
656 672 pwcode = numpy.sum(self.code[0]**2)
657 673 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
658 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
674 normFactor = self.nProfiles * self.nIncohInt * \
675 self.nCohInt * pwcode * self.windowOfFilter
659 676
660 677 return normFactor
661 678
662 679 def getFlagCspc(self):
663 680
664 681 if self.data_cspc is None:
665 682 return True
666 683
667 684 return False
668 685
669 686 def getFlagDc(self):
670 687
671 688 if self.data_dc is None:
672 689 return True
673 690
674 691 return False
675 692
676 693 def getTimeInterval(self):
677 694
678 695 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
679 696
680 697 return timeInterval
681 698
682 699 def getPower(self):
683 700
684 701 factor = self.normFactor
685 702 z = self.data_spc/factor
686 703 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
687 704 avg = numpy.average(z, axis=1)
688 705
689 706 return 10*numpy.log10(avg)
690 707
691 708 def getCoherence(self, pairsList=None, phase=False):
692 709
693 710 z = []
694 711 if pairsList is None:
695 712 pairsIndexList = self.pairsIndexList
696 713 else:
697 714 pairsIndexList = []
698 715 for pair in pairsList:
699 716 if pair not in self.pairsList:
700 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
717 raise ValueError, "Pair %s is not in dataOut.pairsList" % (
718 pair)
701 719 pairsIndexList.append(self.pairsList.index(pair))
702 720 for i in range(len(pairsIndexList)):
703 721 pair = self.pairsList[pairsIndexList[i]]
704 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
722 ccf = numpy.average(
723 self.data_cspc[pairsIndexList[i], :, :], axis=0)
705 724 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
706 725 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
707 726 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
708 727 if phase:
709 728 data = numpy.arctan2(avgcoherenceComplex.imag,
710 729 avgcoherenceComplex.real)*180/numpy.pi
711 730 else:
712 731 data = numpy.abs(avgcoherenceComplex)
713 732
714 733 z.append(data)
715 734
716 735 return numpy.array(z)
717 736
718 737 def setValue(self, value):
719 738
720 739 print "This property should not be initialized"
721 740
722 741 return
723 742
724 743 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
725 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
726 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
744 pairsIndexList = property(
745 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
746 normFactor = property(getNormFactor, setValue,
747 "I'm the 'getNormFactor' property.")
727 748 flag_cspc = property(getFlagCspc, setValue)
728 749 flag_dc = property(getFlagDc, setValue)
729 750 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
730 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
751 timeInterval = property(getTimeInterval, setValue,
752 "I'm the 'timeInterval' property")
753
731 754
732 755 class SpectraHeis(Spectra):
733 756
734 757 data_spc = None
735 758
736 759 data_cspc = None
737 760
738 761 data_dc = None
739 762
740 763 nFFTPoints = None
741 764
742 765 # nPairs = None
743 766
744 767 pairsList = None
745 768
746 769 nCohInt = None
747 770
748 771 nIncohInt = None
749 772
750 773 def __init__(self):
751 774
752 775 self.radarControllerHeaderObj = RadarControllerHeader()
753 776
754 777 self.systemHeaderObj = SystemHeader()
755 778
756 779 self.type = "SpectraHeis"
757 780
758 781 # self.dtype = None
759 782
760 783 # self.nChannels = 0
761 784
762 785 # self.nHeights = 0
763 786
764 787 self.nProfiles = None
765 788
766 789 self.heightList = None
767 790
768 791 self.channelList = None
769 792
770 793 # self.channelIndexList = None
771 794
772 795 self.flagNoData = True
773 796
774 797 self.flagDiscontinuousBlock = False
775 798
776 799 # self.nPairs = 0
777 800
778 801 self.utctime = None
779 802
780 803 self.blocksize = None
781 804
782 805 self.profileIndex = 0
783 806
784 807 self.nCohInt = 1
785 808
786 809 self.nIncohInt = 1
787 810
788 811 def getNormFactor(self):
789 812 pwcode = 1
790 813 if self.flagDecodeData:
791 814 pwcode = numpy.sum(self.code[0]**2)
792 815
793 816 normFactor = self.nIncohInt*self.nCohInt*pwcode
794 817
795 818 return normFactor
796 819
797 820 def getTimeInterval(self):
798 821
799 822 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
800 823
801 824 return timeInterval
802 825
803 826 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
804 827 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
805 828
829
806 830 class Fits(JROData):
807 831
808 832 heightList = None
809 833
810 834 channelList = None
811 835
812 836 flagNoData = True
813 837
814 838 flagDiscontinuousBlock = False
815 839
816 840 useLocalTime = False
817 841
818 842 utctime = None
819 843
820 844 timeZone = None
821 845
822 846 # ippSeconds = None
823 847
824 848 # timeInterval = None
825 849
826 850 nCohInt = None
827 851
828 852 nIncohInt = None
829 853
830 854 noise = None
831 855
832 856 windowOfFilter = 1
833 857
834 858 #Speed of ligth
835 859 C = 3e8
836 860
837 861 frequency = 49.92e6
838 862
839 863 realtime = False
840 864
841
842 865 def __init__(self):
843 866
844 867 self.type = "Fits"
845 868
846 869 self.nProfiles = None
847 870
848 871 self.heightList = None
849 872
850 873 self.channelList = None
851 874
852 875 # self.channelIndexList = None
853 876
854 877 self.flagNoData = True
855 878
856 879 self.utctime = None
857 880
858 881 self.nCohInt = 1
859 882
860 883 self.nIncohInt = 1
861 884
862 885 self.useLocalTime = True
863 886
864 887 self.profileIndex = 0
865 888
866 889 # self.utctime = None
867 890 # self.timeZone = None
868 891 # self.ltctime = None
869 892 # self.timeInterval = None
870 893 # self.header = None
871 894 # self.data_header = None
872 895 # self.data = None
873 896 # self.datatime = None
874 897 # self.flagNoData = False
875 898 # self.expName = ''
876 899 # self.nChannels = None
877 900 # self.nSamples = None
878 901 # self.dataBlocksPerFile = None
879 902 # self.comments = ''
880 903 #
881 904
882
883 905 def getltctime(self):
884 906
885 907 if self.useLocalTime:
886 908 return self.utctime - self.timeZone*60
887 909
888 910 return self.utctime
889 911
890 912 def getDatatime(self):
891 913
892 914 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
893 915 return datatime
894 916
895 917 def getTimeRange(self):
896 918
897 919 datatime = []
898 920
899 921 datatime.append(self.ltctime)
900 922 datatime.append(self.ltctime + self.timeInterval)
901 923
902 924 datatime = numpy.array(datatime)
903 925
904 926 return datatime
905 927
906 928 def getHeiRange(self):
907 929
908 930 heis = self.heightList
909 931
910 932 return heis
911 933
912 934 def getNHeights(self):
913 935
914 936 return len(self.heightList)
915 937
916 938 def getNChannels(self):
917 939
918 940 return len(self.channelList)
919 941
920 942 def getChannelIndexList(self):
921 943
922 944 return range(self.nChannels)
923 945
924 946 def getNoise(self, type = 1):
925 947
926 948 #noise = numpy.zeros(self.nChannels)
927 949
928 950 if type == 1:
929 951 noise = self.getNoisebyHildebrand()
930 952
931 953 if type == 2:
932 954 noise = self.getNoisebySort()
933 955
934 956 if type == 3:
935 957 noise = self.getNoisebyWindow()
936 958
937 959 return noise
938 960
939 961 def getTimeInterval(self):
940 962
941 963 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
942 964
943 965 return timeInterval
944 966
945 967 datatime = property(getDatatime, "I'm the 'datatime' property")
946 968 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
947 969 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
948 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
970 channelIndexList = property(
971 getChannelIndexList, "I'm the 'channelIndexList' property.")
949 972 noise = property(getNoise, "I'm the 'nHeights' property.")
950 973
951 974 ltctime = property(getltctime, "I'm the 'ltctime' property")
952 975 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
953 976
954 977
955 978 class Correlation(JROData):
956 979
957 980 noise = None
958 981
959 982 SNR = None
960 983
961 984 #--------------------------------------------------
962 985
963 986 mode = None
964 987
965 988 split = False
966 989
967 990 data_cf = None
968 991
969 992 lags = None
970 993
971 994 lagRange = None
972 995
973 996 pairsList = None
974 997
975 998 normFactor = None
976 999
977 1000 #--------------------------------------------------
978 1001
979 1002 # calculateVelocity = None
980 1003
981 1004 nLags = None
982 1005
983 1006 nPairs = None
984 1007
985 1008 nAvg = None
986 1009
987
988 1010 def __init__(self):
989 1011 '''
990 1012 Constructor
991 1013 '''
992 1014 self.radarControllerHeaderObj = RadarControllerHeader()
993 1015
994 1016 self.systemHeaderObj = SystemHeader()
995 1017
996 1018 self.type = "Correlation"
997 1019
998 1020 self.data = None
999 1021
1000 1022 self.dtype = None
1001 1023
1002 1024 self.nProfiles = None
1003 1025
1004 1026 self.heightList = None
1005 1027
1006 1028 self.channelList = None
1007 1029
1008 1030 self.flagNoData = True
1009 1031
1010 1032 self.flagDiscontinuousBlock = False
1011 1033
1012 1034 self.utctime = None
1013 1035
1014 1036 self.timeZone = None
1015 1037
1016 1038 self.dstFlag = None
1017 1039
1018 1040 self.errorCount = None
1019 1041
1020 1042 self.blocksize = None
1021 1043
1022 1044 self.flagDecodeData = False #asumo q la data no esta decodificada
1023 1045
1024 1046 self.flagDeflipData = False #asumo q la data no esta sin flip
1025 1047
1026 1048 self.pairsList = None
1027 1049
1028 1050 self.nPoints = None
1029 1051
1030 1052 def getPairsList(self):
1031 1053
1032 1054 return self.pairsList
1033 1055
1034 1056 def getNoise(self, mode = 2):
1035 1057
1036 1058 indR = numpy.where(self.lagR == 0)[0][0]
1037 1059 indT = numpy.where(self.lagT == 0)[0][0]
1038 1060
1039 1061 jspectra0 = self.data_corr[:,:,indR,:]
1040 1062 jspectra = copy.copy(jspectra0)
1041 1063
1042 1064 num_chan = jspectra.shape[0]
1043 1065 num_hei = jspectra.shape[2]
1044 1066
1045 1067 freq_dc = jspectra.shape[1]/2
1046 1068 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1047 1069
1048 1070 if ind_vel[0]<0:
1049 1071 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1050 1072
1051 1073 if mode == 1:
1052 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1074 jspectra[:, freq_dc, :] = (
1075 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
1053 1076
1054 1077 if mode == 2:
1055 1078
1056 1079 vel = numpy.array([-2,-1,1,2])
1057 1080 xx = numpy.zeros([4,4])
1058 1081
1059 1082 for fil in range(4):
1060 1083 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1061 1084
1062 1085 xx_inv = numpy.linalg.inv(xx)
1063 1086 xx_aux = xx_inv[0,:]
1064 1087
1065 1088 for ich in range(num_chan):
1066 1089 yy = jspectra[ich,ind_vel,:]
1067 1090 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1068 1091
1069 1092 junkid = jspectra[ich,freq_dc,:]<=0
1070 1093 cjunkid = sum(junkid)
1071 1094
1072 1095 if cjunkid.any():
1073 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1096 jspectra[ich, freq_dc, junkid.nonzero()] = (
1097 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
1074 1098
1075 1099 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1076 1100
1077 1101 return noise
1078 1102
1079 1103 def getTimeInterval(self):
1080 1104
1081 1105 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1082 1106
1083 1107 return timeInterval
1084 1108
1085 1109 def splitFunctions(self):
1086 1110
1087 1111 pairsList = self.pairsList
1088 1112 ccf_pairs = []
1089 1113 acf_pairs = []
1090 1114 ccf_ind = []
1091 1115 acf_ind = []
1092 1116 for l in range(len(pairsList)):
1093 1117 chan0 = pairsList[l][0]
1094 1118 chan1 = pairsList[l][1]
1095 1119
1096 1120 #Obteniendo pares de Autocorrelacion
1097 1121 if chan0 == chan1:
1098 1122 acf_pairs.append(chan0)
1099 1123 acf_ind.append(l)
1100 1124 else:
1101 1125 ccf_pairs.append(pairsList[l])
1102 1126 ccf_ind.append(l)
1103 1127
1104 1128 data_acf = self.data_cf[acf_ind]
1105 1129 data_ccf = self.data_cf[ccf_ind]
1106 1130
1107 1131 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1108 1132
1109 1133 def getNormFactor(self):
1110 1134 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1111 1135 acf_pairs = numpy.array(acf_pairs)
1112 1136 normFactor = numpy.zeros((self.nPairs,self.nHeights))
1113 1137
1114 1138 for p in range(self.nPairs):
1115 1139 pair = self.pairsList[p]
1116 1140
1117 1141 ch0 = pair[0]
1118 1142 ch1 = pair[1]
1119 1143
1120 1144 ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1)
1121 1145 ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1)
1122 1146 normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max)
1123 1147
1124 1148 return normFactor
1125 1149
1126 1150 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1127 1151 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1128 1152
1153
1129 1154 class Parameters(Spectra):
1130 1155
1131 1156 experimentInfo = None #Information about the experiment
1132 1157
1133 1158 #Information from previous data
1134 1159
1135 1160 inputUnit = None #Type of data to be processed
1136 1161
1137 1162 operation = None #Type of operation to parametrize
1138 1163
1139 1164 #normFactor = None #Normalization Factor
1140 1165
1141 1166 groupList = None #List of Pairs, Groups, etc
1142 1167
1143 1168 #Parameters
1144 1169
1145 1170 data_param = None #Parameters obtained
1146 1171
1147 1172 data_pre = None #Data Pre Parametrization
1148 1173
1149 1174 data_SNR = None #Signal to Noise Ratio
1150 1175
1151 1176 # heightRange = None #Heights
1152 1177
1153 1178 abscissaList = None #Abscissa, can be velocities, lags or time
1154 1179
1155 1180 # noise = None #Noise Potency
1156 1181
1157 1182 utctimeInit = None #Initial UTC time
1158 1183
1159 1184 paramInterval = None #Time interval to calculate Parameters in seconds
1160 1185
1161 1186 useLocalTime = True
1162 1187
1163 1188 #Fitting
1164 1189
1165 1190 data_error = None #Error of the estimation
1166 1191
1167 1192 constants = None
1168 1193
1169 1194 library = None
1170 1195
1171 1196 #Output signal
1172 1197
1173 1198 outputInterval = None #Time interval to calculate output signal in seconds
1174 1199
1175 1200 data_output = None #Out signal
1176 1201
1177 1202 nAvg = None
1178 1203
1179 1204 noise_estimation = None
1180 1205
1181 1206 GauSPC = None #Fit gaussian SPC
1182 1207
1183
1184 1208 def __init__(self):
1185 1209 '''
1186 1210 Constructor
1187 1211 '''
1188 1212 self.radarControllerHeaderObj = RadarControllerHeader()
1189 1213
1190 1214 self.systemHeaderObj = SystemHeader()
1191 1215
1192 1216 self.type = "Parameters"
1193 1217
1194 1218 def getTimeRange1(self, interval):
1195 1219
1196 1220 datatime = []
1197 1221
1198 1222 if self.useLocalTime:
1199 1223 time1 = self.utctimeInit - self.timeZone*60
1200 1224 else:
1201 1225 time1 = self.utctimeInit
1202 1226
1203 1227 datatime.append(time1)
1204 1228 datatime.append(time1 + interval)
1205 1229 datatime = numpy.array(datatime)
1206 1230
1207 1231 return datatime
1208 1232
1209 1233 def getTimeInterval(self):
1210 1234
1211 1235 if hasattr(self, 'timeInterval1'):
1212 1236 return self.timeInterval1
1213 1237 else:
1214 1238 return self.paramInterval
1215 1239
1216 1240 def setValue(self, value):
1217 1241
1218 1242 print "This property should not be initialized"
1219 1243
1220 1244 return
1221 1245
1222 1246 def getNoise(self):
1223 1247
1224 1248 return self.spc_noise
1225 1249
1226 1250 timeInterval = property(getTimeInterval)
1227 1251 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
Show file before | Show file after | Hide comments
@@ -1,751 +1,790
1 1
2 2 '''
3 3 Created on Jul 3, 2014
4 4
5 5 @author: roj-idl71
6 6 '''
7 7 # SUBCHANNELS EN VEZ DE CHANNELS
8 8 # BENCHMARKS -> PROBLEMAS CON ARCHIVOS GRANDES -> INCONSTANTE EN EL TIEMPO
9 9 # ACTUALIZACION DE VERSION
10 10 # HEADERS
11 11 # MODULO DE ESCRITURA
12 12 # METADATA
13 13
14 14 import os
15 15 import datetime
16 16 import numpy
17 17 import timeit
18 from profilehooks import coverage, profile
19 18 from fractions import Fraction
20 19
21 20 try:
22 21 from gevent import sleep
23 22 except:
24 23 from time import sleep
25 24
26 25 from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader
27 26 from schainpy.model.data.jrodata import Voltage
28 27 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
29 28 from time import time
30 29
31 30 import cPickle
32 31 try:
33 32 import digital_rf
34 33 except:
35 34 print 'You should install "digital_rf" module if you want to read Digital RF data'
36 35
36
37 37 class DigitalRFReader(ProcessingUnit):
38 38 '''
39 39 classdocs
40 40 '''
41 41
42 42 def __init__(self, **kwargs):
43 43 '''
44 44 Constructor
45 45 '''
46 46
47 47 ProcessingUnit.__init__(self, **kwargs)
48 48
49 49 self.dataOut = Voltage()
50 50 self.__printInfo = True
51 51 self.__flagDiscontinuousBlock = False
52 52 self.__bufferIndex = 9999999
53 53 self.__ippKm = None
54 54 self.__codeType = 0
55 55 self.__nCode = None
56 56 self.__nBaud = None
57 57 self.__code = None
58 58 self.dtype = None
59 59
60 60 def close(self):
61 61 print 'Average of writing to digital rf format is ', self.oldAverage * 1000
62 62 return
63 63
64 64 def __getCurrentSecond(self):
65 65
66 66 return self.__thisUnixSample/self.__sample_rate
67 67
68 68 thisSecond = property(__getCurrentSecond, "I'm the 'thisSecond' property.")
69 69
70 70 def __setFileHeader(self):
71 71 '''
72 72 In this method will be initialized every parameter of dataOut object (header, no data)
73 73 '''
74 74 ippSeconds = 1.0*self.__nSamples/self.__sample_rate
75 75
76 76 nProfiles = 1.0/ippSeconds # Number of profiles in one second
77 77
78 78 try:
79 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(self.__radarControllerHeader)
79 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(
80 self.__radarControllerHeader)
80 81 except:
81 82 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(
82 83 txA=0,
83 84 txB=0,
84 85 nWindows=1,
85 86 nHeights=self.__nSamples,
86 87 firstHeight=self.__firstHeigth,
87 88 deltaHeight=self.__deltaHeigth,
88 89 codeType=self.__codeType,
89 90 nCode=self.__nCode, nBaud=self.__nBaud,
90 91 code = self.__code)
91 92
92 93 try:
93 94 self.dataOut.systemHeaderObj = SystemHeader(self.__systemHeader)
94 95 except:
95 96 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
96 97 nProfiles=nProfiles,
97 nChannels=len(self.__channelList),
98 nChannels=len(
99 self.__channelList),
98 100 adcResolution=14)
99 101 self.dataOut.type = "Voltage"
100 102
101 103 self.dataOut.data = None
102 104
103 105 self.dataOut.dtype = self.dtype
104 106
105 107 # self.dataOut.nChannels = 0
106 108
107 109 # self.dataOut.nHeights = 0
108 110
109 111 self.dataOut.nProfiles = int(nProfiles)
110 112
111 self.dataOut.heightList = self.__firstHeigth + numpy.arange(self.__nSamples, dtype = numpy.float)*self.__deltaHeigth
113 self.dataOut.heightList = self.__firstHeigth + \
114 numpy.arange(self.__nSamples, dtype=numpy.float) * \
115 self.__deltaHeigth
112 116
113 117 self.dataOut.channelList = range(self.__num_subchannels)
114 118
115 119 self.dataOut.blocksize = self.dataOut.getNChannels() * self.dataOut.getNHeights()
116 120
117 121 # self.dataOut.channelIndexList = None
118 122
119 123 self.dataOut.flagNoData = True
120 124
121 125 self.dataOut.flagDataAsBlock = False
122 126 # Set to TRUE if the data is discontinuous
123 127 self.dataOut.flagDiscontinuousBlock = False
124 128
125 129 self.dataOut.utctime = None
126 130
127 self.dataOut.timeZone = self.__timezone/60 # timezone like jroheader, difference in minutes between UTC and localtime
131 # timezone like jroheader, difference in minutes between UTC and localtime
132 self.dataOut.timeZone = self.__timezone / 60
128 133
129 134 self.dataOut.dstFlag = 0
130 135
131 136 self.dataOut.errorCount = 0
132 137
133 138 try:
134 self.dataOut.nCohInt = self.fixed_metadata_dict.get('nCohInt', 1)
139 self.dataOut.nCohInt = self.fixed_metadata_dict.get(
140 'nCohInt', self.nCohInt)
135 141
136 self.dataOut.flagDecodeData = self.fixed_metadata_dict['flagDecodeData'] # asumo que la data esta decodificada
142 # asumo que la data esta decodificada
143 self.dataOut.flagDecodeData = self.fixed_metadata_dict.get(
144 'flagDecodeData', self.flagDecodeData)
137 145
138 self.dataOut.flagDeflipData = self.fixed_metadata_dict['flagDeflipData'] # asumo que la data esta sin flip
146 # asumo que la data esta sin flip
147 self.dataOut.flagDeflipData = self.fixed_metadata_dict['flagDeflipData']
139 148
140 149 self.dataOut.flagShiftFFT = self.fixed_metadata_dict['flagShiftFFT']
141 150
142 151 self.dataOut.useLocalTime = self.fixed_metadata_dict['useLocalTime']
143 152 except:
144 153 pass
145 154
146
147 155 self.dataOut.ippSeconds = ippSeconds
148 156
149 157 # Time interval between profiles
150 158 # self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt
151 159
152 160 self.dataOut.frequency = self.__frequency
153 161
154 162 self.dataOut.realtime = self.__online
155 163
156 164 def findDatafiles(self, path, startDate=None, endDate=None):
157 165
158 166 if not os.path.isdir(path):
159 167 return []
160 168
161 169 try:
162 digitalReadObj = digital_rf.DigitalRFReader(path, load_all_metadata=True)
170 digitalReadObj = digital_rf.DigitalRFReader(
171 path, load_all_metadata=True)
163 172 except:
164 173 digitalReadObj = digital_rf.DigitalRFReader(path)
165 174
166 175 channelNameList = digitalReadObj.get_channels()
167 176
168 177 if not channelNameList:
169 178 return []
170 179
171 180 metadata_dict = digitalReadObj.get_rf_file_metadata(channelNameList[0])
172 181
173 182 sample_rate = metadata_dict['sample_rate'][0]
174 183
175 184 this_metadata_file = digitalReadObj.get_metadata(channelNameList[0])
176 185
177 186 try:
178 187 timezone = this_metadata_file['timezone'].value
179 188 except:
180 189 timezone = 0
181 190
182 startUTCSecond, endUTCSecond = digitalReadObj.get_bounds(channelNameList[0])/sample_rate - timezone
191 startUTCSecond, endUTCSecond = digitalReadObj.get_bounds(
192 channelNameList[0]) / sample_rate - timezone
183 193
184 194 startDatetime = datetime.datetime.utcfromtimestamp(startUTCSecond)
185 195 endDatatime = datetime.datetime.utcfromtimestamp(endUTCSecond)
186 196
187 197 if not startDate:
188 198 startDate = startDatetime.date()
189 199
190 200 if not endDate:
191 201 endDate = endDatatime.date()
192 202
193 203 dateList = []
194 204
195 205 thisDatetime = startDatetime
196 206
197 207 while(thisDatetime<=endDatatime):
198 208
199 209 thisDate = thisDatetime.date()
200 210
201 211 if thisDate < startDate:
202 212 continue
203 213
204 214 if thisDate > endDate:
205 215 break
206 216
207 217 dateList.append(thisDate)
208 218 thisDatetime += datetime.timedelta(1)
209 219
210 220 return dateList
211 221
212 222 def setup(self, path = None,
213 223 startDate = None,
214 224 endDate = None,
215 225 startTime = datetime.time(0,0,0),
216 226 endTime = datetime.time(23,59,59),
217 227 channelList = None,
218 228 nSamples = None,
219 229 online = False,
220 230 delay = 60,
221 231 buffer_size = 1024,
222 232 ippKm=None,
233 nCohInt=1,
234 nCode=1,
235 nBaud=1,
236 flagDecodeData=False,
237 code=numpy.ones((1, 1), dtype=numpy.int),
223 238 **kwargs):
224 239 '''
225 240 In this method we should set all initial parameters.
226 241
227 242 Inputs:
228 243 path
229 244 startDate
230 245 endDate
231 246 startTime
232 247 endTime
233 248 set
234 249 expLabel
235 250 ext
236 251 online
237 252 delay
238 253 '''
254 self.nCohInt = nCohInt
255 self.flagDecodeData = flagDecodeData
239 256 self.i = 0
240 257 if not os.path.isdir(path):
241 258 raise ValueError, "[Reading] Directory %s does not exist" %path
242 259
243 260 try:
244 self.digitalReadObj = digital_rf.DigitalRFReader(path, load_all_metadata=True)
261 self.digitalReadObj = digital_rf.DigitalRFReader(
262 path, load_all_metadata=True)
245 263 except:
246 264 self.digitalReadObj = digital_rf.DigitalRFReader(path)
247 265
248 266 channelNameList = self.digitalReadObj.get_channels()
249 267
250 268 if not channelNameList:
251 269 raise ValueError, "[Reading] Directory %s does not have any files" %path
252 270
253 271 if not channelList:
254 272 channelList = range(len(channelNameList))
255 273
256
257 274 ########## Reading metadata ######################
258 275
259 top_properties = self.digitalReadObj.get_properties(channelNameList[channelList[0]])
260
276 top_properties = self.digitalReadObj.get_properties(
277 channelNameList[channelList[0]])
261 278
262 279 self.__num_subchannels = top_properties['num_subchannels']
263 self.__sample_rate = 1.0 * top_properties['sample_rate_numerator'] / top_properties['sample_rate_denominator']
280 self.__sample_rate = 1.0 * \
281 top_properties['sample_rate_numerator'] / \
282 top_properties['sample_rate_denominator']
264 283 # self.__samples_per_file = top_properties['samples_per_file'][0]
265 self.__deltaHeigth = 1e6*0.15/self.__sample_rate ## why 0.15?
284 self.__deltaHeigth = 1e6 * 0.15 / self.__sample_rate # why 0.15?
266 285
267 this_metadata_file = self.digitalReadObj.get_digital_metadata(channelNameList[channelList[0]])
286 this_metadata_file = self.digitalReadObj.get_digital_metadata(
287 channelNameList[channelList[0]])
268 288 metadata_bounds = this_metadata_file.get_bounds()
269 self.fixed_metadata_dict = this_metadata_file.read(metadata_bounds[0])[metadata_bounds[0]] ## GET FIRST HEADER
289 self.fixed_metadata_dict = this_metadata_file.read(
290 metadata_bounds[0])[metadata_bounds[0]] # GET FIRST HEADER
270 291
271 292 try:
272 293 self.__processingHeader = self.fixed_metadata_dict['processingHeader']
273 294 self.__radarControllerHeader = self.fixed_metadata_dict['radarControllerHeader']
274 295 self.__systemHeader = self.fixed_metadata_dict['systemHeader']
275 296 self.dtype = cPickle.loads(self.fixed_metadata_dict['dtype'])
276 297 except:
277 298 pass
278 299
279
280 300 self.__frequency = None
281 301
282 302 self.__frequency = self.fixed_metadata_dict.get('frequency', 1)
283 303
284 304 self.__timezone = self.fixed_metadata_dict.get('timezone', 300)
285 305
286
287 306 try:
288 307 nSamples = self.fixed_metadata_dict['nSamples']
289 308 except:
290 309 nSamples = None
291 310
292 311 self.__firstHeigth = 0
293 312
294 313 try:
295 314 codeType = self.__radarControllerHeader['codeType']
296 315 except:
297 316 codeType = 0
298 317
299 nCode = 1
300 nBaud = 1
301 code = numpy.ones((nCode, nBaud), dtype=numpy.int)
302
303 318 try:
304 319 if codeType:
305 320 nCode = self.__radarControllerHeader['nCode']
306 321 nBaud = self.__radarControllerHeader['nBaud']
307 322 code = self.__radarControllerHeader['code']
308 323 except:
309 324 pass
310 325
311
312 326 if not ippKm:
313 327 try:
314 328 # seconds to km
315 329 ippKm = self.__radarControllerHeader['ipp']
316 330 except:
317 331 ippKm = None
318 332 ####################################################
319 333 self.__ippKm = ippKm
320 334 startUTCSecond = None
321 335 endUTCSecond = None
322 336
323 337 if startDate:
324 338 startDatetime = datetime.datetime.combine(startDate, startTime)
325 startUTCSecond = (startDatetime-datetime.datetime(1970,1,1)).total_seconds() + self.__timezone
339 startUTCSecond = (
340 startDatetime - datetime.datetime(1970, 1, 1)).total_seconds() + self.__timezone
326 341
327 342 if endDate:
328 343 endDatetime = datetime.datetime.combine(endDate, endTime)
329 endUTCSecond = (endDatetime-datetime.datetime(1970,1,1)).total_seconds() + self.__timezone
344 endUTCSecond = (endDatetime - datetime.datetime(1970,
345 1, 1)).total_seconds() + self.__timezone
330 346
331 start_index, end_index = self.digitalReadObj.get_bounds(channelNameList[channelList[0]])
347 start_index, end_index = self.digitalReadObj.get_bounds(
348 channelNameList[channelList[0]])
332 349
333 350 if not startUTCSecond:
334 351 startUTCSecond = start_index/self.__sample_rate
335 352
336 353 if start_index > startUTCSecond*self.__sample_rate:
337 354 startUTCSecond = start_index/self.__sample_rate
338 355
339 356 if not endUTCSecond:
340 357 endUTCSecond = end_index/self.__sample_rate
341 358
342 359 if end_index < endUTCSecond*self.__sample_rate:
343 360 endUTCSecond = end_index/self.__sample_rate
344 361 if not nSamples:
345 362 if not ippKm:
346 363 raise ValueError, "[Reading] nSamples or ippKm should be defined"
347 364 nSamples = int(ippKm / (1e6*0.15/self.__sample_rate))
348 365 channelBoundList = []
349 366 channelNameListFiltered = []
350 367
351 368 for thisIndexChannel in channelList:
352 369 thisChannelName = channelNameList[thisIndexChannel]
353 start_index, end_index = self.digitalReadObj.get_bounds(thisChannelName)
370 start_index, end_index = self.digitalReadObj.get_bounds(
371 thisChannelName)
354 372 channelBoundList.append((start_index, end_index))
355 373 channelNameListFiltered.append(thisChannelName)
356 374
357 375 self.profileIndex = 0
358 376 self.i= 0
359 377 self.__delay = delay
360 378
361 379 self.__codeType = codeType
362 380 self.__nCode = nCode
363 381 self.__nBaud = nBaud
364 382 self.__code = code
365 383
366 384 self.__datapath = path
367 385 self.__online = online
368 386 self.__channelList = channelList
369 387 self.__channelNameList = channelNameListFiltered
370 388 self.__channelBoundList = channelBoundList
371 389 self.__nSamples = nSamples
372 390 self.__samples_to_read = long(nSamples) # FIJO: AHORA 40
373 391 self.__nChannels = len(self.__channelList)
374 392
375 393 self.__startUTCSecond = startUTCSecond
376 394 self.__endUTCSecond = endUTCSecond
377 395
378 self.__timeInterval = 1.0 * self.__samples_to_read/self.__sample_rate # Time interval
396 self.__timeInterval = 1.0 * self.__samples_to_read / \
397 self.__sample_rate # Time interval
379 398
380 399 if online:
381 400 # self.__thisUnixSample = int(endUTCSecond*self.__sample_rate - 4*self.__samples_to_read)
382 401 startUTCSecond = numpy.floor(endUTCSecond)
383 402
384 self.__thisUnixSample = long(startUTCSecond*self.__sample_rate) - self.__samples_to_read ## por que en el otro metodo lo primero q se hace es sumar samplestoread
403 # por que en el otro metodo lo primero q se hace es sumar samplestoread
404 self.__thisUnixSample = long(
405 startUTCSecond * self.__sample_rate) - self.__samples_to_read
385 406
386 self.__data_buffer = numpy.zeros((self.__num_subchannels, self.__samples_to_read), dtype = numpy.complex)
407 self.__data_buffer = numpy.zeros(
408 (self.__num_subchannels, self.__samples_to_read), dtype=numpy.complex)
387 409
388 410 self.__setFileHeader()
389 411 self.isConfig = True
390 412
391 413 print "[Reading] Digital RF Data was found from %s to %s " %(
392 datetime.datetime.utcfromtimestamp(self.__startUTCSecond - self.__timezone),
393 datetime.datetime.utcfromtimestamp(self.__endUTCSecond - self.__timezone)
414 datetime.datetime.utcfromtimestamp(
415 self.__startUTCSecond - self.__timezone),
416 datetime.datetime.utcfromtimestamp(
417 self.__endUTCSecond - self.__timezone)
394 418 )
395 419
396 420 print "[Reading] Starting process from %s to %s" %(datetime.datetime.utcfromtimestamp(startUTCSecond - self.__timezone),
397 datetime.datetime.utcfromtimestamp(endUTCSecond - self.__timezone)
421 datetime.datetime.utcfromtimestamp(
422 endUTCSecond - self.__timezone)
398 423 )
399 424 self.oldAverage = None
400 425 self.count = 0
401 426 self.executionTime = 0
427
402 428 def __reload(self):
403 429 # print
404 430 # print "%s not in range [%s, %s]" %(
405 431 # datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
406 432 # datetime.datetime.utcfromtimestamp(self.__startUTCSecond - self.__timezone),
407 433 # datetime.datetime.utcfromtimestamp(self.__endUTCSecond - self.__timezone)
408 434 # )
409 435 print "[Reading] reloading metadata ..."
410 436
411 437 try:
412 438 self.digitalReadObj.reload(complete_update=True)
413 439 except:
414 440 self.digitalReadObj.reload()
415 441
416 start_index, end_index = self.digitalReadObj.get_bounds(self.__channelNameList[self.__channelList[0]])
442 start_index, end_index = self.digitalReadObj.get_bounds(
443 self.__channelNameList[self.__channelList[0]])
417 444
418 445 if start_index > self.__startUTCSecond*self.__sample_rate:
419 446 self.__startUTCSecond = 1.0*start_index/self.__sample_rate
420 447
421 448 if end_index > self.__endUTCSecond*self.__sample_rate:
422 449 self.__endUTCSecond = 1.0*end_index/self.__sample_rate
423 450 print
424 451 print "[Reading] New timerange found [%s, %s] " %(
425 datetime.datetime.utcfromtimestamp(self.__startUTCSecond - self.__timezone),
426 datetime.datetime.utcfromtimestamp(self.__endUTCSecond - self.__timezone)
452 datetime.datetime.utcfromtimestamp(
453 self.__startUTCSecond - self.__timezone),
454 datetime.datetime.utcfromtimestamp(
455 self.__endUTCSecond - self.__timezone)
427 456 )
428 457
429 458 return True
430 459
431 460 return False
432 461
433 462 def timeit(self, toExecute):
434 463 t0 = time()
435 464 toExecute()
436 465 self.executionTime = time() - t0
437 if self.oldAverage is None: self.oldAverage = self.executionTime
438 self.oldAverage = (self.executionTime + self.count*self.oldAverage) / (self.count + 1.0)
466 if self.oldAverage is None:
467 self.oldAverage = self.executionTime
468 self.oldAverage = (self.executionTime + self.count *
469 self.oldAverage) / (self.count + 1.0)
439 470 self.count = self.count + 1.0
440 471 return
441 472
442 473 def __readNextBlock(self, seconds=30, volt_scale = 1):
443 474 '''
444 475 '''
445 476
446 477 # Set the next data
447 478 self.__flagDiscontinuousBlock = False
448 479 self.__thisUnixSample += self.__samples_to_read
449 480
450 481 if self.__thisUnixSample + 2*self.__samples_to_read > self.__endUTCSecond*self.__sample_rate:
451 482 print "[Reading] There are no more data into selected time-range"
452 483 if self.__online:
453 484 self.__reload()
454 485 else:
455 486 return False
456 487
457 488 if self.__thisUnixSample + 2*self.__samples_to_read > self.__endUTCSecond*self.__sample_rate:
458 489 return False
459 490 self.__thisUnixSample -= self.__samples_to_read
460 491
461 492 indexChannel = 0
462 493
463 494 dataOk = False
464 for thisChannelName in self.__channelNameList: ##TODO VARIOS CHANNELS?
495 for thisChannelName in self.__channelNameList: # TODO VARIOS CHANNELS?
465 496 for indexSubchannel in range(self.__num_subchannels):
466 497 try:
467 498 t0 = time()
468 499 result = self.digitalReadObj.read_vector_c81d(self.__thisUnixSample,
469 500 self.__samples_to_read,
470 501 thisChannelName, sub_channel=indexSubchannel)
471 502 self.executionTime = time() - t0
472 if self.oldAverage is None: self.oldAverage = self.executionTime
473 self.oldAverage = (self.executionTime + self.count*self.oldAverage) / (self.count + 1.0)
503 if self.oldAverage is None:
504 self.oldAverage = self.executionTime
505 self.oldAverage = (
506 self.executionTime + self.count * self.oldAverage) / (self.count + 1.0)
474 507 self.count = self.count + 1.0
475 508
476 509 except IOError, e:
477 510 #read next profile
478 511 self.__flagDiscontinuousBlock = True
479 512 print "[Reading] %s" %datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone), e
480 513 break
481 514
482 515 if result.shape[0] != self.__samples_to_read:
483 516 self.__flagDiscontinuousBlock = True
484 517 print "[Reading] %s: Too few samples were found, just %d/%d samples" %(datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
485 518 result.shape[0],
486 519 self.__samples_to_read)
487 520 break
488 521
489 522 self.__data_buffer[indexSubchannel,:] = result*volt_scale
490 523
491 524 indexChannel += 1
492 525
493 526 dataOk = True
494 527
495 528 self.__utctime = self.__thisUnixSample/self.__sample_rate
496 529
497 530 if not dataOk:
498 531 return False
499 532
500 533 print "[Reading] %s: %d samples <> %f sec" %(datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
501 534 self.__samples_to_read,
502 535 self.__timeInterval)
503 536
504 537 self.__bufferIndex = 0
505 538
506 539 return True
507 540
508 541 def __isBufferEmpty(self):
509 542 return self.__bufferIndex > self.__samples_to_read - self.__nSamples #40960 - 40
510 543
511 544 def getData(self, seconds=30, nTries=5):
512
513 545 '''
514 546 This method gets the data from files and put the data into the dataOut object
515 547
516 548 In addition, increase el the buffer counter in one.
517 549
518 550 Return:
519 551 data : retorna un perfil de voltages (alturas * canales) copiados desde el
520 552 buffer. Si no hay mas archivos a leer retorna None.
521 553
522 554 Affected:
523 555 self.dataOut
524 556 self.profileIndex
525 557 self.flagDiscontinuousBlock
526 558 self.flagIsNewBlock
527 559 '''
528 560
529 561 err_counter = 0
530 562 self.dataOut.flagNoData = True
531 563
532 564 if self.__isBufferEmpty():
533 565 self.__flagDiscontinuousBlock = False
534 566
535 567 while True:
536 568 if self.__readNextBlock():
537 569 break
538 570 if self.__thisUnixSample > self.__endUTCSecond*self.__sample_rate:
539 571 return False
540 572
541 573 if self.__flagDiscontinuousBlock:
542 574 print '[Reading] discontinuous block found ... continue with the next block'
543 575 continue
544 576
545 577 if not self.__online:
546 578 return False
547 579
548 580 err_counter += 1
549 581 if err_counter > nTries:
550 582 return False
551 583
552 584 print '[Reading] waiting %d seconds to read a new block' %seconds
553 585 sleep(seconds)
554 586
555 self.dataOut.data = self.__data_buffer[:,self.__bufferIndex:self.__bufferIndex+self.__nSamples]
556 self.dataOut.utctime = (self.__thisUnixSample + self.__bufferIndex)/self.__sample_rate
587 self.dataOut.data = self.__data_buffer[:,
588 self.__bufferIndex:self.__bufferIndex + self.__nSamples]
589 self.dataOut.utctime = (
590 self.__thisUnixSample + self.__bufferIndex) / self.__sample_rate
557 591 self.dataOut.flagNoData = False
558 592 self.dataOut.flagDiscontinuousBlock = self.__flagDiscontinuousBlock
559 593 self.dataOut.profileIndex = self.profileIndex
560 594
561 595 self.__bufferIndex += self.__nSamples
562 596 self.profileIndex += 1
563 597
564 598 if self.profileIndex == self.dataOut.nProfiles:
565 599 self.profileIndex = 0
566 600
567 601 return True
568 602
569 603 def printInfo(self):
570 604 '''
571 605 '''
572 606 if self.__printInfo == False:
573 607 return
574 608
575 609 # self.systemHeaderObj.printInfo()
576 610 # self.radarControllerHeaderObj.printInfo()
577 611
578 612 self.__printInfo = False
579 613
580 614 def printNumberOfBlock(self):
581 615 '''
582 616 '''
583 617 return
584 618 # print self.profileIndex
585 619
586
587 620 def run(self, **kwargs):
588 621 '''
589 622 This method will be called many times so here you should put all your code
590 623 '''
591 624
592 625 if not self.isConfig:
593 626 self.setup(**kwargs)
594 627 #self.i = self.i+1
595 628 self.getData(seconds=self.__delay)
596 629
597 630 return
598 631
632
599 633 class DigitalRFWriter(Operation):
600 634 '''
601 635 classdocs
602 636 '''
603 637
604 638 def __init__(self, **kwargs):
605 639 '''
606 640 Constructor
607 641 '''
608 642 Operation.__init__(self, **kwargs)
609 643 self.metadata_dict = {}
610 644 self.dataOut = None
611 645 self.dtype = None
612 646
613 647 def setHeader(self):
614 648
615 649 self.metadata_dict['frequency'] = self.dataOut.frequency
616 650 self.metadata_dict['timezone'] = self.dataOut.timeZone
617 651 self.metadata_dict['dtype'] = cPickle.dumps(self.dataOut.dtype)
618 652 self.metadata_dict['nProfiles'] = self.dataOut.nProfiles
619 653 self.metadata_dict['heightList'] = self.dataOut.heightList
620 654 self.metadata_dict['channelList'] = self.dataOut.channelList
621 655 self.metadata_dict['flagDecodeData'] = self.dataOut.flagDecodeData
622 656 self.metadata_dict['flagDeflipData'] = self.dataOut.flagDeflipData
623 657 self.metadata_dict['flagShiftFFT'] = self.dataOut.flagShiftFFT
624 658 self.metadata_dict['flagDataAsBlock'] = self.dataOut.flagDataAsBlock
625 659 self.metadata_dict['useLocalTime'] = self.dataOut.useLocalTime
626 660 self.metadata_dict['nCohInt'] = self.dataOut.nCohInt
627 661
628 662 return
629 663
630 664 def setup(self, dataOut, path, frequency, fileCadence, dirCadence, metadataCadence, set=0, metadataFile='metadata', ext='.h5'):
631 665 '''
632 666 In this method we should set all initial parameters.
633 667 Input:
634 668 dataOut: Input data will also be outputa data
635 669 '''
636 670 self.setHeader()
637 671 self.__ippSeconds = dataOut.ippSeconds
638 672 self.__deltaH = dataOut.getDeltaH()
639 673 self.__sample_rate = 1e6*0.15/self.__deltaH
640 674 self.__dtype = dataOut.dtype
641 675 if len(dataOut.dtype) == 2:
642 676 self.__dtype = dataOut.dtype[0]
643 677 self.__nSamples = dataOut.systemHeaderObj.nSamples
644 678 self.__nProfiles = dataOut.nProfiles
645 679 self.__blocks_per_file = dataOut.processingHeaderObj.dataBlocksPerFile
646 680
647 self.arr_data = arr_data = numpy.ones((self.__nSamples, len(self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
681 self.arr_data = arr_data = numpy.ones((self.__nSamples, len(
682 self.dataOut.channelList)), dtype=[('r', self.__dtype), ('i', self.__dtype)])
648 683
649 file_cadence_millisecs = long(1.0 * self.__blocks_per_file * self.__nProfiles * self.__nSamples / self.__sample_rate) * 1000
684 file_cadence_millisecs = long(
685 1.0 * self.__blocks_per_file * self.__nProfiles * self.__nSamples / self.__sample_rate) * 1000
650 686 sub_cadence_secs = file_cadence_millisecs / 500
651 687
652 688 sample_rate_fraction = Fraction(self.__sample_rate).limit_denominator()
653 689 sample_rate_numerator = long(sample_rate_fraction.numerator)
654 690 sample_rate_denominator = long(sample_rate_fraction.denominator)
655 691 start_global_index = dataOut.utctime * self.__sample_rate
656 692
657 693 uuid = 'prueba'
658 694 compression_level = 1
659 695 checksum = False
660 696 is_complex = True
661 697 num_subchannels = len(dataOut.channelList)
662 698 is_continuous = True
663 699 marching_periods = False
664 700
665 701 self.digitalWriteObj = digital_rf.DigitalRFWriter(path, self.__dtype, dirCadence,
666 702 fileCadence, start_global_index,
667 703 sample_rate_numerator, sample_rate_denominator, uuid, compression_level, checksum,
668 704 is_complex, num_subchannels, is_continuous, marching_periods)
669 705
670 706 metadata_dir = os.path.join(path, 'metadata')
671 707 os.system('mkdir %s' % (metadata_dir))
672 708
673 self.digitalMetadataWriteObj = digital_rf.DigitalMetadataWriter(metadata_dir, dirCadence, 1, ##236, file_cadence_millisecs / 1000
709 self.digitalMetadataWriteObj = digital_rf.DigitalMetadataWriter(metadata_dir, dirCadence, 1, # 236, file_cadence_millisecs / 1000
674 710 sample_rate_numerator, sample_rate_denominator,
675 711 metadataFile)
676 712
677
678 713 self.isConfig = True
679 714 self.currentSample = 0
680 715 self.oldAverage = 0
681 716 self.count = 0
682 717 return
683 718
684 719 def writeMetadata(self):
685 720 print '[Writing] - Writing metadata'
686 721 start_idx = self.__sample_rate * self.dataOut.utctime
687 722
688 self.metadata_dict['processingHeader'] = self.dataOut.processingHeaderObj.getAsDict()
689 self.metadata_dict['radarControllerHeader'] = self.dataOut.radarControllerHeaderObj.getAsDict()
690 self.metadata_dict['systemHeader'] = self.dataOut.systemHeaderObj.getAsDict()
723 self.metadata_dict['processingHeader'] = self.dataOut.processingHeaderObj.getAsDict(
724 )
725 self.metadata_dict['radarControllerHeader'] = self.dataOut.radarControllerHeaderObj.getAsDict(
726 )
727 self.metadata_dict['systemHeader'] = self.dataOut.systemHeaderObj.getAsDict(
728 )
691 729 self.digitalMetadataWriteObj.write(start_idx, self.metadata_dict)
692 730 return
693 731
694
695 732 def timeit(self, toExecute):
696 733 t0 = time()
697 734 toExecute()
698 735 self.executionTime = time() - t0
699 if self.oldAverage is None: self.oldAverage = self.executionTime
700 self.oldAverage = (self.executionTime + self.count*self.oldAverage) / (self.count + 1.0)
736 if self.oldAverage is None:
737 self.oldAverage = self.executionTime
738 self.oldAverage = (self.executionTime + self.count *
739 self.oldAverage) / (self.count + 1.0)
701 740 self.count = self.count + 1.0
702 741 return
703 742
704
705 743 def writeData(self):
706 744 for i in range(self.dataOut.systemHeaderObj.nSamples):
707 745 for channel in self.dataOut.channelList:
708 746 self.arr_data[i][channel]['r'] = self.dataOut.data[channel][i].real
709 747 self.arr_data[i][channel]['i'] = self.dataOut.data[channel][i].imag
710 748
711 749 def f(): return self.digitalWriteObj.rf_write(self.arr_data)
712 750 self.timeit(f)
713 751
714 752 return
715 753
716 754 def run(self, dataOut, frequency=49.92e6, path=None, fileCadence=100, dirCadence=25, metadataCadence=1, **kwargs):
717 755 '''
718 756 This method will be called many times so here you should put all your code
719 757 Inputs:
720 758 dataOut: object with the data
721 759 '''
722 760 # print dataOut.__dict__
723 761 self.dataOut = dataOut
724 762 if not self.isConfig:
725 self.setup(dataOut, path, frequency, fileCadence, dirCadence, metadataCadence, **kwargs)
763 self.setup(dataOut, path, frequency, fileCadence,
764 dirCadence, metadataCadence, **kwargs)
726 765 self.writeMetadata()
727 766
728 767 self.writeData()
729 768
730 769 ## self.currentSample += 1
731 ## if self.dataOut.flagDataAsBlock or self.currentSample == 1:
732 ## self.writeMetadata()
770 # if self.dataOut.flagDataAsBlock or self.currentSample == 1:
771 # self.writeMetadata()
733 772 ## if self.currentSample == self.__nProfiles: self.currentSample = 0
734 773
735 774 def close(self):
736 775 print '[Writing] - Closing files '
737 776 print 'Average of writing to digital rf format is ', self.oldAverage * 1000
738 777 try:
739 778 self.digitalWriteObj.close()
740 779 except:
741 780 pass
742 781
743 782 # raise
744 783 if __name__ == '__main__':
745 784
746 785 readObj = DigitalRFReader()
747 786
748 787 while True:
749 788 readObj.run(path='/home/jchavez/jicamarca/mocked_data/')
750 789 # readObj.printInfo()
751 790 # readObj.printNumberOfBlock()
Show file before | Show file after | Hide comments
@@ -1,1283 +1,1281
1 1 import sys
2 2 import numpy
3 3 from scipy import interpolate
4 4 from schainpy import cSchain
5 5 from jroproc_base import ProcessingUnit, Operation
6 6 from schainpy.model.data.jrodata import Voltage
7 7 from time import time
8 8
9 9 class VoltageProc(ProcessingUnit):
10 10
11 11
12 12 def __init__(self, **kwargs):
13 13
14 14 ProcessingUnit.__init__(self, **kwargs)
15 15
16 16 # self.objectDict = {}
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19
20 20 def run(self):
21 21 if self.dataIn.type == 'AMISR':
22 22 self.__updateObjFromAmisrInput()
23 23
24 24 if self.dataIn.type == 'Voltage':
25 25 self.dataOut.copy(self.dataIn)
26 26
27 27 # self.dataOut.copy(self.dataIn)
28 28
29 29 def __updateObjFromAmisrInput(self):
30 30
31 31 self.dataOut.timeZone = self.dataIn.timeZone
32 32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 33 self.dataOut.errorCount = self.dataIn.errorCount
34 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 35
36 36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 37 self.dataOut.data = self.dataIn.data
38 38 self.dataOut.utctime = self.dataIn.utctime
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 # self.dataOut.timeInterval = self.dataIn.timeInterval
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.nProfiles = self.dataIn.nProfiles
43 43
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 46 self.dataOut.frequency = self.dataIn.frequency
47 47
48 48 self.dataOut.azimuth = self.dataIn.azimuth
49 49 self.dataOut.zenith = self.dataIn.zenith
50 50
51 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54 54 #
55 55 # pass#
56 56 #
57 57 # def init(self):
58 58 #
59 59 #
60 60 # if self.dataIn.type == 'AMISR':
61 61 # self.__updateObjFromAmisrInput()
62 62 #
63 63 # if self.dataIn.type == 'Voltage':
64 64 # self.dataOut.copy(self.dataIn)
65 65 # # No necesita copiar en cada init() los atributos de dataIn
66 66 # # la copia deberia hacerse por cada nuevo bloque de datos
67 67
68 68 def selectChannels(self, channelList):
69 69
70 70 channelIndexList = []
71 71
72 72 for channel in channelList:
73 73 if channel not in self.dataOut.channelList:
74 74 raise ValueError, "Channel %d is not in %s" %(channel, str(self.dataOut.channelList))
75 75
76 76 index = self.dataOut.channelList.index(channel)
77 77 channelIndexList.append(index)
78 78
79 79 self.selectChannelsByIndex(channelIndexList)
80 80
81 81 def selectChannelsByIndex(self, channelIndexList):
82 82 """
83 83 Selecciona un bloque de datos en base a canales segun el channelIndexList
84 84
85 85 Input:
86 86 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
87 87
88 88 Affected:
89 89 self.dataOut.data
90 90 self.dataOut.channelIndexList
91 91 self.dataOut.nChannels
92 92 self.dataOut.m_ProcessingHeader.totalSpectra
93 93 self.dataOut.systemHeaderObj.numChannels
94 94 self.dataOut.m_ProcessingHeader.blockSize
95 95
96 96 Return:
97 97 None
98 98 """
99 99
100 100 for channelIndex in channelIndexList:
101 101 if channelIndex not in self.dataOut.channelIndexList:
102 102 print channelIndexList
103 103 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
104 104
105 105 if self.dataOut.flagDataAsBlock:
106 106 """
107 107 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
108 108 """
109 109 data = self.dataOut.data[channelIndexList,:,:]
110 110 else:
111 111 data = self.dataOut.data[channelIndexList,:]
112 112
113 113 self.dataOut.data = data
114 114 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
115 115 # self.dataOut.nChannels = nChannels
116 116
117 117 return 1
118 118
119 119 def selectHeights(self, minHei=None, maxHei=None):
120 120 """
121 121 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
122 122 minHei <= height <= maxHei
123 123
124 124 Input:
125 125 minHei : valor minimo de altura a considerar
126 126 maxHei : valor maximo de altura a considerar
127 127
128 128 Affected:
129 129 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
130 130
131 131 Return:
132 132 1 si el metodo se ejecuto con exito caso contrario devuelve 0
133 133 """
134 134
135 135 if minHei == None:
136 136 minHei = self.dataOut.heightList[0]
137 137
138 138 if maxHei == None:
139 139 maxHei = self.dataOut.heightList[-1]
140 140
141 141 if (minHei < self.dataOut.heightList[0]):
142 142 minHei = self.dataOut.heightList[0]
143 143
144 144 if (maxHei > self.dataOut.heightList[-1]):
145 145 maxHei = self.dataOut.heightList[-1]
146 146
147 147 minIndex = 0
148 148 maxIndex = 0
149 149 heights = self.dataOut.heightList
150 150
151 151 inda = numpy.where(heights >= minHei)
152 152 indb = numpy.where(heights <= maxHei)
153 153
154 154 try:
155 155 minIndex = inda[0][0]
156 156 except:
157 157 minIndex = 0
158 158
159 159 try:
160 160 maxIndex = indb[0][-1]
161 161 except:
162 162 maxIndex = len(heights)
163 163
164 164 self.selectHeightsByIndex(minIndex, maxIndex)
165 165
166 166 return 1
167 167
168 168
169 169 def selectHeightsByIndex(self, minIndex, maxIndex):
170 170 """
171 171 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
172 172 minIndex <= index <= maxIndex
173 173
174 174 Input:
175 175 minIndex : valor de indice minimo de altura a considerar
176 176 maxIndex : valor de indice maximo de altura a considerar
177 177
178 178 Affected:
179 179 self.dataOut.data
180 180 self.dataOut.heightList
181 181
182 182 Return:
183 183 1 si el metodo se ejecuto con exito caso contrario devuelve 0
184 184 """
185 185
186 186 if (minIndex < 0) or (minIndex > maxIndex):
187 187 raise ValueError, "Height index range (%d,%d) is not valid" % (minIndex, maxIndex)
188 188
189 189 if (maxIndex >= self.dataOut.nHeights):
190 190 maxIndex = self.dataOut.nHeights
191 191
192 192 #voltage
193 193 if self.dataOut.flagDataAsBlock:
194 194 """
195 195 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
196 196 """
197 197 data = self.dataOut.data[:,:, minIndex:maxIndex]
198 198 else:
199 199 data = self.dataOut.data[:, minIndex:maxIndex]
200 200
201 201 # firstHeight = self.dataOut.heightList[minIndex]
202 202
203 203 self.dataOut.data = data
204 204 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
205 205
206 206 if self.dataOut.nHeights <= 1:
207 207 raise ValueError, "selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights)
208 208
209 209 return 1
210 210
211 211
212 212 def filterByHeights(self, window):
213 213
214 214 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
215 215
216 216 if window == None:
217 217 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
218 218
219 219 newdelta = deltaHeight * window
220 220 r = self.dataOut.nHeights % window
221 221 newheights = (self.dataOut.nHeights-r)/window
222 222
223 223 if newheights <= 1:
224 224 raise ValueError, "filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(self.dataOut.nHeights, window)
225 225
226 226 if self.dataOut.flagDataAsBlock:
227 227 """
228 228 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
229 229 """
230 230 buffer = self.dataOut.data[:, :, 0:self.dataOut.nHeights-r]
231 231 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nProfiles,self.dataOut.nHeights/window,window)
232 232 buffer = numpy.sum(buffer,3)
233 233
234 234 else:
235 235 buffer = self.dataOut.data[:,0:self.dataOut.nHeights-r]
236 236 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nHeights/window,window)
237 237 buffer = numpy.sum(buffer,2)
238 238
239 239 self.dataOut.data = buffer
240 240 self.dataOut.heightList = self.dataOut.heightList[0] + numpy.arange( newheights )*newdelta
241 241 self.dataOut.windowOfFilter = window
242 242
243 243 def setH0(self, h0, deltaHeight = None):
244 244
245 245 if not deltaHeight:
246 246 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
247 247
248 248 nHeights = self.dataOut.nHeights
249 249
250 250 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
251 251
252 252 self.dataOut.heightList = newHeiRange
253 253
254 254 def deFlip(self, channelList = []):
255 255
256 256 data = self.dataOut.data.copy()
257 257
258 258 if self.dataOut.flagDataAsBlock:
259 259 flip = self.flip
260 260 profileList = range(self.dataOut.nProfiles)
261 261
262 262 if not channelList:
263 263 for thisProfile in profileList:
264 264 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
265 265 flip *= -1.0
266 266 else:
267 267 for thisChannel in channelList:
268 268 if thisChannel not in self.dataOut.channelList:
269 269 continue
270 270
271 271 for thisProfile in profileList:
272 272 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
273 273 flip *= -1.0
274 274
275 275 self.flip = flip
276 276
277 277 else:
278 278 if not channelList:
279 279 data[:,:] = data[:,:]*self.flip
280 280 else:
281 281 for thisChannel in channelList:
282 282 if thisChannel not in self.dataOut.channelList:
283 283 continue
284 284
285 285 data[thisChannel,:] = data[thisChannel,:]*self.flip
286 286
287 287 self.flip *= -1.
288 288
289 289 self.dataOut.data = data
290 290
291 291 def setRadarFrequency(self, frequency=None):
292 292
293 293 if frequency != None:
294 294 self.dataOut.frequency = frequency
295 295
296 296 return 1
297 297
298 298 def interpolateHeights(self, topLim, botLim):
299 299 #69 al 72 para julia
300 300 #82-84 para meteoros
301 301 if len(numpy.shape(self.dataOut.data))==2:
302 302 sampInterp = (self.dataOut.data[:,botLim-1] + self.dataOut.data[:,topLim+1])/2
303 303 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
304 304 #self.dataOut.data[:,botLim:limSup+1] = sampInterp
305 305 self.dataOut.data[:,botLim:topLim+1] = sampInterp
306 306 else:
307 307 nHeights = self.dataOut.data.shape[2]
308 308 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
309 309 y = self.dataOut.data[:,:,range(botLim)+range(topLim+1,nHeights)]
310 310 f = interpolate.interp1d(x, y, axis = 2)
311 311 xnew = numpy.arange(botLim,topLim+1)
312 312 ynew = f(xnew)
313 313
314 314 self.dataOut.data[:,:,botLim:topLim+1] = ynew
315 315
316 316 # import collections
317 317
318 318 class CohInt(Operation):
319 319
320 320 isConfig = False
321 321
322 322 __profIndex = 0
323 323 __withOverapping = False
324 324
325 325 __byTime = False
326 326 __initime = None
327 327 __lastdatatime = None
328 328 __integrationtime = None
329 329
330 330 __buffer = None
331 331
332 332 __dataReady = False
333 333
334 334 n = None
335 335
336 336 def __init__(self, **kwargs):
337 337
338 338 Operation.__init__(self, **kwargs)
339 339
340 340 # self.isConfig = False
341 341
342 342 def setup(self, n=None, timeInterval=None, overlapping=False, byblock=False):
343 343 """
344 344 Set the parameters of the integration class.
345 345
346 346 Inputs:
347 347
348 348 n : Number of coherent integrations
349 349 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
350 350 overlapping :
351 351 """
352 352
353 353 self.__initime = None
354 354 self.__lastdatatime = 0
355 355 self.__buffer = None
356 356 self.__dataReady = False
357 357 self.byblock = byblock
358 358
359 359 if n == None and timeInterval == None:
360 360 raise ValueError, "n or timeInterval should be specified ..."
361 361
362 362 if n != None:
363 363 self.n = n
364 364 self.__byTime = False
365 365 else:
366 366 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
367 367 self.n = 9999
368 368 self.__byTime = True
369 369
370 370 if overlapping:
371 371 self.__withOverapping = True
372 372 self.__buffer = None
373 373 else:
374 374 self.__withOverapping = False
375 375 self.__buffer = 0
376 376
377 377 self.__profIndex = 0
378 378
379 379 def putData(self, data):
380 380
381 381 """
382 382 Add a profile to the __buffer and increase in one the __profileIndex
383 383
384 384 """
385 385
386 386 if not self.__withOverapping:
387 387 self.__buffer += data.copy()
388 388 self.__profIndex += 1
389 389 return
390 390
391 391 #Overlapping data
392 392 nChannels, nHeis = data.shape
393 393 data = numpy.reshape(data, (1, nChannels, nHeis))
394 394
395 395 #If the buffer is empty then it takes the data value
396 396 if self.__buffer is None:
397 397 self.__buffer = data
398 398 self.__profIndex += 1
399 399 return
400 400
401 401 #If the buffer length is lower than n then stakcing the data value
402 402 if self.__profIndex < self.n:
403 403 self.__buffer = numpy.vstack((self.__buffer, data))
404 404 self.__profIndex += 1
405 405 return
406 406
407 407 #If the buffer length is equal to n then replacing the last buffer value with the data value
408 408 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
409 409 self.__buffer[self.n-1] = data
410 410 self.__profIndex = self.n
411 411 return
412 412
413 413
414 414 def pushData(self):
415 415 """
416 416 Return the sum of the last profiles and the profiles used in the sum.
417 417
418 418 Affected:
419 419
420 420 self.__profileIndex
421 421
422 422 """
423 423
424 424 if not self.__withOverapping:
425 425 data = self.__buffer
426 426 n = self.__profIndex
427 427
428 428 self.__buffer = 0
429 429 self.__profIndex = 0
430 430
431 431 return data, n
432 432
433 433 #Integration with Overlapping
434 434 data = numpy.sum(self.__buffer, axis=0)
435 435 n = self.__profIndex
436 436
437 437 return data, n
438 438
439 439 def byProfiles(self, data):
440 440
441 441 self.__dataReady = False
442 442 avgdata = None
443 443 # n = None
444 444
445 445 self.putData(data)
446 446
447 447 if self.__profIndex == self.n:
448 448
449 449 avgdata, n = self.pushData()
450 450 self.__dataReady = True
451 451
452 452 return avgdata
453 453
454 454 def byTime(self, data, datatime):
455 455
456 456 self.__dataReady = False
457 457 avgdata = None
458 458 n = None
459 459
460 460 self.putData(data)
461 461
462 462 if (datatime - self.__initime) >= self.__integrationtime:
463 463 avgdata, n = self.pushData()
464 464 self.n = n
465 465 self.__dataReady = True
466 466
467 467 return avgdata
468 468
469 469 def integrate(self, data, datatime=None):
470 470
471 471 if self.__initime == None:
472 472 self.__initime = datatime
473 473
474 474 if self.__byTime:
475 475 avgdata = self.byTime(data, datatime)
476 476 else:
477 477 avgdata = self.byProfiles(data)
478 478
479 479
480 480 self.__lastdatatime = datatime
481 481
482 482 if avgdata is None:
483 483 return None, None
484 484
485 485 avgdatatime = self.__initime
486 486
487 487 deltatime = datatime -self.__lastdatatime
488 488
489 489 if not self.__withOverapping:
490 490 self.__initime = datatime
491 491 else:
492 492 self.__initime += deltatime
493 493
494 494 return avgdata, avgdatatime
495 495
496 496 def integrateByBlock(self, dataOut):
497 497
498 498 times = int(dataOut.data.shape[1]/self.n)
499 499 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
500 500
501 501 id_min = 0
502 502 id_max = self.n
503 503
504 504 for i in range(times):
505 505 junk = dataOut.data[:,id_min:id_max,:]
506 506 avgdata[:,i,:] = junk.sum(axis=1)
507 507 id_min += self.n
508 508 id_max += self.n
509 509
510 510 timeInterval = dataOut.ippSeconds*self.n
511 511 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
512 512 self.__dataReady = True
513 513 return avgdata, avgdatatime
514 514
515 515
516 516 def run(self, dataOut, n=None, timeInterval=None, overlapping=False, byblock=False, **kwargs):
517 517 if not self.isConfig:
518 518 self.setup(n=n, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
519 519 self.isConfig = True
520 520
521 521 if dataOut.flagDataAsBlock:
522 522 """
523 523 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
524 524 """
525 525 avgdata, avgdatatime = self.integrateByBlock(dataOut)
526 526 dataOut.nProfiles /= self.n
527 527 else:
528 528 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
529 529
530 530 # dataOut.timeInterval *= n
531 531 dataOut.flagNoData = True
532 532
533 533 if self.__dataReady:
534 534 dataOut.data = avgdata
535 535 dataOut.nCohInt *= self.n
536 536 dataOut.utctime = avgdatatime
537 537 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
538 538 dataOut.flagNoData = False
539 539
540 540 class Decoder(Operation):
541 541
542 542 isConfig = False
543 543 __profIndex = 0
544 544
545 545 code = None
546 546
547 547 nCode = None
548 548 nBaud = None
549 549
550 550 def __init__(self, **kwargs):
551 551
552 552 Operation.__init__(self, **kwargs)
553 553
554 554 self.times = None
555 555 self.osamp = None
556 556 # self.__setValues = False
557 557 self.isConfig = False
558 558
559 559 def setup(self, code, osamp, dataOut):
560 560
561 561 self.__profIndex = 0
562 562
563 563 self.code = code
564 564
565 565 self.nCode = len(code)
566 566 self.nBaud = len(code[0])
567 567
568 568 if (osamp != None) and (osamp >1):
569 569 self.osamp = osamp
570 570 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
571 571 self.nBaud = self.nBaud*self.osamp
572 572
573 573 self.__nChannels = dataOut.nChannels
574 574 self.__nProfiles = dataOut.nProfiles
575 575 self.__nHeis = dataOut.nHeights
576 576
577 577 if self.__nHeis < self.nBaud:
578 578 raise ValueError, 'Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud)
579 579
580 580 #Frequency
581 581 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
582 582
583 583 __codeBuffer[:,0:self.nBaud] = self.code
584 584
585 585 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
586 586
587 587 if dataOut.flagDataAsBlock:
588 588
589 589 self.ndatadec = self.__nHeis #- self.nBaud + 1
590 590
591 591 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
592 592
593 593 else:
594 594
595 595 #Time
596 596 self.ndatadec = self.__nHeis #- self.nBaud + 1
597 597
598 598 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
599 599
600 600 def __convolutionInFreq(self, data):
601 601
602 602 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
603 603
604 604 fft_data = numpy.fft.fft(data, axis=1)
605 605
606 606 conv = fft_data*fft_code
607 607
608 608 data = numpy.fft.ifft(conv,axis=1)
609 609
610 610 return data
611 611
612 612 def __convolutionInFreqOpt(self, data):
613 613
614 614 raise NotImplementedError
615 615
616 616 def __convolutionInTime(self, data):
617 617
618 618 code = self.code[self.__profIndex]
619
620 619 for i in range(self.__nChannels):
621 620 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
622 621
623 622 return self.datadecTime
624 623
625 624 def __convolutionByBlockInTime(self, data):
626 625
627 626 repetitions = self.__nProfiles / self.nCode
628 627
629 628 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
630 629 junk = junk.flatten()
631 630 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
632 631 profilesList = xrange(self.__nProfiles)
633 632
634 633 for i in range(self.__nChannels):
635 634 for j in profilesList:
636 635 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
637 636 return self.datadecTime
638 637
639 638 def __convolutionByBlockInFreq(self, data):
640 639
641 640 raise NotImplementedError, "Decoder by frequency fro Blocks not implemented"
642 641
643 642
644 643 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
645 644
646 645 fft_data = numpy.fft.fft(data, axis=2)
647 646
648 647 conv = fft_data*fft_code
649 648
650 649 data = numpy.fft.ifft(conv,axis=2)
651 650
652 651 return data
653 652
654 653
655 654 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
656 655
657 656 if dataOut.flagDecodeData:
658 657 print "This data is already decoded, recoding again ..."
659 658
660 659 if not self.isConfig:
661 660
662 661 if code is None:
663 662 if dataOut.code is None:
664 663 raise ValueError, "Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type
665 664
666 665 code = dataOut.code
667 666 else:
668 667 code = numpy.array(code).reshape(nCode,nBaud)
669
670 668 self.setup(code, osamp, dataOut)
671 669
672 670 self.isConfig = True
673 671
674 672 if mode == 3:
675 673 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
676 674
677 675 if times != None:
678 676 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
679 677
680 678 if self.code is None:
681 679 print "Fail decoding: Code is not defined."
682 680 return
683 681
684 682 self.__nProfiles = dataOut.nProfiles
685 683 datadec = None
686 684
687 685 if mode == 3:
688 686 mode = 0
689 687
690 688 if dataOut.flagDataAsBlock:
691 689 """
692 690 Decoding when data have been read as block,
693 691 """
694 692
695 693 if mode == 0:
696 694 datadec = self.__convolutionByBlockInTime(dataOut.data)
697 695 if mode == 1:
698 696 datadec = self.__convolutionByBlockInFreq(dataOut.data)
699 697 else:
700 698 """
701 699 Decoding when data have been read profile by profile
702 700 """
703 701 if mode == 0:
704 702 datadec = self.__convolutionInTime(dataOut.data)
705 703
706 704 if mode == 1:
707 705 datadec = self.__convolutionInFreq(dataOut.data)
708 706
709 707 if mode == 2:
710 708 datadec = self.__convolutionInFreqOpt(dataOut.data)
711 709
712 710 if datadec is None:
713 711 raise ValueError, "Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode
714 712
715 713 dataOut.code = self.code
716 714 dataOut.nCode = self.nCode
717 715 dataOut.nBaud = self.nBaud
718 716
719 717 dataOut.data = datadec
720 718
721 719 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
722 720
723 721 dataOut.flagDecodeData = True #asumo q la data esta decodificada
724 722
725 723 if self.__profIndex == self.nCode-1:
726 724 self.__profIndex = 0
727 725 return 1
728 726
729 727 self.__profIndex += 1
730 728
731 729 return 1
732 730 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
733 731
734 732
735 733 class ProfileConcat(Operation):
736 734
737 735 isConfig = False
738 736 buffer = None
739 737
740 738 def __init__(self, **kwargs):
741 739
742 740 Operation.__init__(self, **kwargs)
743 741 self.profileIndex = 0
744 742
745 743 def reset(self):
746 744 self.buffer = numpy.zeros_like(self.buffer)
747 745 self.start_index = 0
748 746 self.times = 1
749 747
750 748 def setup(self, data, m, n=1):
751 749 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
752 750 self.nHeights = data.shape[1]#.nHeights
753 751 self.start_index = 0
754 752 self.times = 1
755 753
756 754 def concat(self, data):
757 755
758 756 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
759 757 self.start_index = self.start_index + self.nHeights
760 758
761 759 def run(self, dataOut, m):
762 760
763 761 dataOut.flagNoData = True
764 762
765 763 if not self.isConfig:
766 764 self.setup(dataOut.data, m, 1)
767 765 self.isConfig = True
768 766
769 767 if dataOut.flagDataAsBlock:
770 768 raise ValueError, "ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False"
771 769
772 770 else:
773 771 self.concat(dataOut.data)
774 772 self.times += 1
775 773 if self.times > m:
776 774 dataOut.data = self.buffer
777 775 self.reset()
778 776 dataOut.flagNoData = False
779 777 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
780 778 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
781 779 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
782 780 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
783 781 dataOut.ippSeconds *= m
784 782
785 783 class ProfileSelector(Operation):
786 784
787 785 profileIndex = None
788 786 # Tamanho total de los perfiles
789 787 nProfiles = None
790 788
791 789 def __init__(self, **kwargs):
792 790
793 791 Operation.__init__(self, **kwargs)
794 792 self.profileIndex = 0
795 793
796 794 def incProfileIndex(self):
797 795
798 796 self.profileIndex += 1
799 797
800 798 if self.profileIndex >= self.nProfiles:
801 799 self.profileIndex = 0
802 800
803 801 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
804 802
805 803 if profileIndex < minIndex:
806 804 return False
807 805
808 806 if profileIndex > maxIndex:
809 807 return False
810 808
811 809 return True
812 810
813 811 def isThisProfileInList(self, profileIndex, profileList):
814 812
815 813 if profileIndex not in profileList:
816 814 return False
817 815
818 816 return True
819 817
820 818 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
821 819
822 820 """
823 821 ProfileSelector:
824 822
825 823 Inputs:
826 824 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
827 825
828 826 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
829 827
830 828 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
831 829
832 830 """
833 831
834 832 if rangeList is not None:
835 833 if type(rangeList[0]) not in (tuple, list):
836 834 rangeList = [rangeList]
837 835
838 836 dataOut.flagNoData = True
839 837
840 838 if dataOut.flagDataAsBlock:
841 839 """
842 840 data dimension = [nChannels, nProfiles, nHeis]
843 841 """
844 842 if profileList != None:
845 843 dataOut.data = dataOut.data[:,profileList,:]
846 844
847 845 if profileRangeList != None:
848 846 minIndex = profileRangeList[0]
849 847 maxIndex = profileRangeList[1]
850 848 profileList = range(minIndex, maxIndex+1)
851 849
852 850 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
853 851
854 852 if rangeList != None:
855 853
856 854 profileList = []
857 855
858 856 for thisRange in rangeList:
859 857 minIndex = thisRange[0]
860 858 maxIndex = thisRange[1]
861 859
862 860 profileList.extend(range(minIndex, maxIndex+1))
863 861
864 862 dataOut.data = dataOut.data[:,profileList,:]
865 863
866 864 dataOut.nProfiles = len(profileList)
867 865 dataOut.profileIndex = dataOut.nProfiles - 1
868 866 dataOut.flagNoData = False
869 867
870 868 return True
871 869
872 870 """
873 871 data dimension = [nChannels, nHeis]
874 872 """
875 873
876 874 if profileList != None:
877 875
878 876 if self.isThisProfileInList(dataOut.profileIndex, profileList):
879 877
880 878 self.nProfiles = len(profileList)
881 879 dataOut.nProfiles = self.nProfiles
882 880 dataOut.profileIndex = self.profileIndex
883 881 dataOut.flagNoData = False
884 882
885 883 self.incProfileIndex()
886 884 return True
887 885
888 886 if profileRangeList != None:
889 887
890 888 minIndex = profileRangeList[0]
891 889 maxIndex = profileRangeList[1]
892 890
893 891 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
894 892
895 893 self.nProfiles = maxIndex - minIndex + 1
896 894 dataOut.nProfiles = self.nProfiles
897 895 dataOut.profileIndex = self.profileIndex
898 896 dataOut.flagNoData = False
899 897
900 898 self.incProfileIndex()
901 899 return True
902 900
903 901 if rangeList != None:
904 902
905 903 nProfiles = 0
906 904
907 905 for thisRange in rangeList:
908 906 minIndex = thisRange[0]
909 907 maxIndex = thisRange[1]
910 908
911 909 nProfiles += maxIndex - minIndex + 1
912 910
913 911 for thisRange in rangeList:
914 912
915 913 minIndex = thisRange[0]
916 914 maxIndex = thisRange[1]
917 915
918 916 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
919 917
920 918 self.nProfiles = nProfiles
921 919 dataOut.nProfiles = self.nProfiles
922 920 dataOut.profileIndex = self.profileIndex
923 921 dataOut.flagNoData = False
924 922
925 923 self.incProfileIndex()
926 924
927 925 break
928 926
929 927 return True
930 928
931 929
932 930 if beam != None: #beam is only for AMISR data
933 931 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
934 932 dataOut.flagNoData = False
935 933 dataOut.profileIndex = self.profileIndex
936 934
937 935 self.incProfileIndex()
938 936
939 937 return True
940 938
941 939 raise ValueError, "ProfileSelector needs profileList, profileRangeList or rangeList parameter"
942 940
943 941 return False
944 942
945 943 class Reshaper(Operation):
946 944
947 945 def __init__(self, **kwargs):
948 946
949 947 Operation.__init__(self, **kwargs)
950 948
951 949 self.__buffer = None
952 950 self.__nitems = 0
953 951
954 952 def __appendProfile(self, dataOut, nTxs):
955 953
956 954 if self.__buffer is None:
957 955 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
958 956 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
959 957
960 958 ini = dataOut.nHeights * self.__nitems
961 959 end = ini + dataOut.nHeights
962 960
963 961 self.__buffer[:, ini:end] = dataOut.data
964 962
965 963 self.__nitems += 1
966 964
967 965 return int(self.__nitems*nTxs)
968 966
969 967 def __getBuffer(self):
970 968
971 969 if self.__nitems == int(1./self.__nTxs):
972 970
973 971 self.__nitems = 0
974 972
975 973 return self.__buffer.copy()
976 974
977 975 return None
978 976
979 977 def __checkInputs(self, dataOut, shape, nTxs):
980 978
981 979 if shape is None and nTxs is None:
982 980 raise ValueError, "Reshaper: shape of factor should be defined"
983 981
984 982 if nTxs:
985 983 if nTxs < 0:
986 984 raise ValueError, "nTxs should be greater than 0"
987 985
988 986 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
989 987 raise ValueError, "nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs))
990 988
991 989 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
992 990
993 991 return shape, nTxs
994 992
995 993 if len(shape) != 2 and len(shape) != 3:
996 994 raise ValueError, "shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights)
997 995
998 996 if len(shape) == 2:
999 997 shape_tuple = [dataOut.nChannels]
1000 998 shape_tuple.extend(shape)
1001 999 else:
1002 1000 shape_tuple = list(shape)
1003 1001
1004 1002 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1005 1003
1006 1004 return shape_tuple, nTxs
1007 1005
1008 1006 def run(self, dataOut, shape=None, nTxs=None):
1009 1007
1010 1008 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1011 1009
1012 1010 dataOut.flagNoData = True
1013 1011 profileIndex = None
1014 1012
1015 1013 if dataOut.flagDataAsBlock:
1016 1014
1017 1015 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1018 1016 dataOut.flagNoData = False
1019 1017
1020 1018 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1021 1019
1022 1020 else:
1023 1021
1024 1022 if self.__nTxs < 1:
1025 1023
1026 1024 self.__appendProfile(dataOut, self.__nTxs)
1027 1025 new_data = self.__getBuffer()
1028 1026
1029 1027 if new_data is not None:
1030 1028 dataOut.data = new_data
1031 1029 dataOut.flagNoData = False
1032 1030
1033 1031 profileIndex = dataOut.profileIndex*nTxs
1034 1032
1035 1033 else:
1036 1034 raise ValueError, "nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)"
1037 1035
1038 1036 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1039 1037
1040 1038 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1041 1039
1042 1040 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1043 1041
1044 1042 dataOut.profileIndex = profileIndex
1045 1043
1046 1044 dataOut.ippSeconds /= self.__nTxs
1047 1045
1048 1046 class SplitProfiles(Operation):
1049 1047
1050 1048 def __init__(self, **kwargs):
1051 1049
1052 1050 Operation.__init__(self, **kwargs)
1053 1051
1054 1052 def run(self, dataOut, n):
1055 1053
1056 1054 dataOut.flagNoData = True
1057 1055 profileIndex = None
1058 1056
1059 1057 if dataOut.flagDataAsBlock:
1060 1058
1061 1059 #nchannels, nprofiles, nsamples
1062 1060 shape = dataOut.data.shape
1063 1061
1064 1062 if shape[2] % n != 0:
1065 1063 raise ValueError, "Could not split the data, n=%d has to be multiple of %d" %(n, shape[2])
1066 1064
1067 1065 new_shape = shape[0], shape[1]*n, shape[2]/n
1068 1066
1069 1067 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1070 1068 dataOut.flagNoData = False
1071 1069
1072 1070 profileIndex = int(dataOut.nProfiles/n) - 1
1073 1071
1074 1072 else:
1075 1073
1076 1074 raise ValueError, "Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)"
1077 1075
1078 1076 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1079 1077
1080 1078 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1081 1079
1082 1080 dataOut.nProfiles = int(dataOut.nProfiles*n)
1083 1081
1084 1082 dataOut.profileIndex = profileIndex
1085 1083
1086 1084 dataOut.ippSeconds /= n
1087 1085
1088 1086 class CombineProfiles(Operation):
1089 1087 def __init__(self, **kwargs):
1090 1088
1091 1089 Operation.__init__(self, **kwargs)
1092 1090
1093 1091 self.__remData = None
1094 1092 self.__profileIndex = 0
1095 1093
1096 1094 def run(self, dataOut, n):
1097 1095
1098 1096 dataOut.flagNoData = True
1099 1097 profileIndex = None
1100 1098
1101 1099 if dataOut.flagDataAsBlock:
1102 1100
1103 1101 #nchannels, nprofiles, nsamples
1104 1102 shape = dataOut.data.shape
1105 1103 new_shape = shape[0], shape[1]/n, shape[2]*n
1106 1104
1107 1105 if shape[1] % n != 0:
1108 1106 raise ValueError, "Could not split the data, n=%d has to be multiple of %d" %(n, shape[1])
1109 1107
1110 1108 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1111 1109 dataOut.flagNoData = False
1112 1110
1113 1111 profileIndex = int(dataOut.nProfiles*n) - 1
1114 1112
1115 1113 else:
1116 1114
1117 1115 #nchannels, nsamples
1118 1116 if self.__remData is None:
1119 1117 newData = dataOut.data
1120 1118 else:
1121 1119 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1122 1120
1123 1121 self.__profileIndex += 1
1124 1122
1125 1123 if self.__profileIndex < n:
1126 1124 self.__remData = newData
1127 1125 #continue
1128 1126 return
1129 1127
1130 1128 self.__profileIndex = 0
1131 1129 self.__remData = None
1132 1130
1133 1131 dataOut.data = newData
1134 1132 dataOut.flagNoData = False
1135 1133
1136 1134 profileIndex = dataOut.profileIndex/n
1137 1135
1138 1136
1139 1137 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1140 1138
1141 1139 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1142 1140
1143 1141 dataOut.nProfiles = int(dataOut.nProfiles/n)
1144 1142
1145 1143 dataOut.profileIndex = profileIndex
1146 1144
1147 1145 dataOut.ippSeconds *= n
1148 1146
1149 1147 # import collections
1150 1148 # from scipy.stats import mode
1151 1149 #
1152 1150 # class Synchronize(Operation):
1153 1151 #
1154 1152 # isConfig = False
1155 1153 # __profIndex = 0
1156 1154 #
1157 1155 # def __init__(self, **kwargs):
1158 1156 #
1159 1157 # Operation.__init__(self, **kwargs)
1160 1158 # # self.isConfig = False
1161 1159 # self.__powBuffer = None
1162 1160 # self.__startIndex = 0
1163 1161 # self.__pulseFound = False
1164 1162 #
1165 1163 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1166 1164 #
1167 1165 # #Read data
1168 1166 #
1169 1167 # powerdB = dataOut.getPower(channel = channel)
1170 1168 # noisedB = dataOut.getNoise(channel = channel)[0]
1171 1169 #
1172 1170 # self.__powBuffer.extend(powerdB.flatten())
1173 1171 #
1174 1172 # dataArray = numpy.array(self.__powBuffer)
1175 1173 #
1176 1174 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1177 1175 #
1178 1176 # maxValue = numpy.nanmax(filteredPower)
1179 1177 #
1180 1178 # if maxValue < noisedB + 10:
1181 1179 # #No se encuentra ningun pulso de transmision
1182 1180 # return None
1183 1181 #
1184 1182 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1185 1183 #
1186 1184 # if len(maxValuesIndex) < 2:
1187 1185 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1188 1186 # return None
1189 1187 #
1190 1188 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1191 1189 #
1192 1190 # #Seleccionar solo valores con un espaciamiento de nSamples
1193 1191 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1194 1192 #
1195 1193 # if len(pulseIndex) < 2:
1196 1194 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1197 1195 # return None
1198 1196 #
1199 1197 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1200 1198 #
1201 1199 # #remover senales que se distancien menos de 10 unidades o muestras
1202 1200 # #(No deberian existir IPP menor a 10 unidades)
1203 1201 #
1204 1202 # realIndex = numpy.where(spacing > 10 )[0]
1205 1203 #
1206 1204 # if len(realIndex) < 2:
1207 1205 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1208 1206 # return None
1209 1207 #
1210 1208 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1211 1209 # realPulseIndex = pulseIndex[realIndex]
1212 1210 #
1213 1211 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1214 1212 #
1215 1213 # print "IPP = %d samples" %period
1216 1214 #
1217 1215 # self.__newNSamples = dataOut.nHeights #int(period)
1218 1216 # self.__startIndex = int(realPulseIndex[0])
1219 1217 #
1220 1218 # return 1
1221 1219 #
1222 1220 #
1223 1221 # def setup(self, nSamples, nChannels, buffer_size = 4):
1224 1222 #
1225 1223 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1226 1224 # maxlen = buffer_size*nSamples)
1227 1225 #
1228 1226 # bufferList = []
1229 1227 #
1230 1228 # for i in range(nChannels):
1231 1229 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1232 1230 # maxlen = buffer_size*nSamples)
1233 1231 #
1234 1232 # bufferList.append(bufferByChannel)
1235 1233 #
1236 1234 # self.__nSamples = nSamples
1237 1235 # self.__nChannels = nChannels
1238 1236 # self.__bufferList = bufferList
1239 1237 #
1240 1238 # def run(self, dataOut, channel = 0):
1241 1239 #
1242 1240 # if not self.isConfig:
1243 1241 # nSamples = dataOut.nHeights
1244 1242 # nChannels = dataOut.nChannels
1245 1243 # self.setup(nSamples, nChannels)
1246 1244 # self.isConfig = True
1247 1245 #
1248 1246 # #Append new data to internal buffer
1249 1247 # for thisChannel in range(self.__nChannels):
1250 1248 # bufferByChannel = self.__bufferList[thisChannel]
1251 1249 # bufferByChannel.extend(dataOut.data[thisChannel])
1252 1250 #
1253 1251 # if self.__pulseFound:
1254 1252 # self.__startIndex -= self.__nSamples
1255 1253 #
1256 1254 # #Finding Tx Pulse
1257 1255 # if not self.__pulseFound:
1258 1256 # indexFound = self.__findTxPulse(dataOut, channel)
1259 1257 #
1260 1258 # if indexFound == None:
1261 1259 # dataOut.flagNoData = True
1262 1260 # return
1263 1261 #
1264 1262 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1265 1263 # self.__pulseFound = True
1266 1264 # self.__startIndex = indexFound
1267 1265 #
1268 1266 # #If pulse was found ...
1269 1267 # for thisChannel in range(self.__nChannels):
1270 1268 # bufferByChannel = self.__bufferList[thisChannel]
1271 1269 # #print self.__startIndex
1272 1270 # x = numpy.array(bufferByChannel)
1273 1271 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1274 1272 #
1275 1273 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1276 1274 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1277 1275 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1278 1276 #
1279 1277 # dataOut.data = self.__arrayBuffer
1280 1278 #
1281 1279 # self.__startIndex += self.__newNSamples
1282 1280 #
1283 1281 # return
General Comments 0
You need to be logged in to leave comments. Login now