##// 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
@@ -30,7 +30,7 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
Show file before | Show file after | Hide comments
@@ -15,41 +15,43 from schainpy import cSchain
15 15 def getNumpyDtype(dataTypeCode):
16 16
17 17 if dataTypeCode == 0:
18 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
18 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
19 19 elif dataTypeCode == 1:
20 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
20 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
21 21 elif dataTypeCode == 2:
22 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
22 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
23 23 elif dataTypeCode == 3:
24 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
24 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
25 25 elif dataTypeCode == 4:
26 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
26 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
27 27 elif dataTypeCode == 5:
28 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
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 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
37 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
37 38 datatype = 0
38 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
39 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
39 40 datatype = 1
40 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
41 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
41 42 datatype = 2
42 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
43 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
43 44 datatype = 3
44 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
45 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
45 46 datatype = 4
46 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
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
@@ -110,6 +112,7 class Beam:
110 112 self.azimuthList = []
111 113 self.zenithList = []
112 114
115
113 116 class GenericData(object):
114 117
115 118 flagNoData = True
@@ -123,12 +126,12 class GenericData(object):
123 126
124 127 attribute = inputObj.__dict__[key]
125 128
126 #If this attribute is a tuple or list
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 #If this attribute is another object or instance
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
@@ -143,10 +146,11 class GenericData(object):
143 146
144 147 return self.flagNoData
145 148
149
146 150 class JROData(GenericData):
147 151
148 # m_BasicHeader = BasicHeader()
149 # m_ProcessingHeader = ProcessingHeader()
152 # m_BasicHeader = BasicHeader()
153 # m_ProcessingHeader = ProcessingHeader()
150 154
151 155 systemHeaderObj = SystemHeader()
152 156
@@ -156,7 +160,7 class JROData(GenericData):
156 160
157 161 type = None
158 162
159 datatype = None #dtype but in string
163 datatype = None # dtype but in string
160 164
161 165 # dtype = None
162 166
@@ -190,9 +194,9 class JROData(GenericData):
190 194 #
191 195 # code = None
192 196
193 flagDecodeData = False #asumo q la data no esta decodificada
197 flagDecodeData = False # asumo q la data no esta decodificada
194 198
195 flagDeflipData = False #asumo q la data no esta sin flip
199 flagDeflipData = False # asumo q la data no esta sin flip
196 200
197 201 flagShiftFFT = False
198 202
@@ -206,7 +210,7 class JROData(GenericData):
206 210
207 211 windowOfFilter = 1
208 212
209 #Speed of ligth
213 # Speed of ligth
210 214 C = 3e8
211 215
212 216 frequency = 49.92e6
@@ -261,7 +265,7 class JROData(GenericData):
261 265 def getltctime(self):
262 266
263 267 if self.useLocalTime:
264 return self.utctime - self.timeZone*60
268 return self.utctime - self.timeZone * 60
265 269
266 270 return self.utctime
267 271
@@ -275,7 +279,7 class JROData(GenericData):
275 279 datatime = []
276 280
277 281 datatime.append(self.ltctime)
278 datatime.append(self.ltctime + self.timeInterval+1)
282 datatime.append(self.ltctime + self.timeInterval + 1)
279 283
280 284 datatime = numpy.array(datatime)
281 285
@@ -283,25 +287,25 class JROData(GenericData):
283 287
284 288 def getFmaxTimeResponse(self):
285 289
286 period = (10**-6)*self.getDeltaH()/(0.15)
290 period = (10**-6) * self.getDeltaH() / (0.15)
287 291
288 PRF = 1./(period * self.nCohInt)
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 PRF = 1./(self.ippSeconds * self.nCohInt)
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 _lambda = self.C/self.frequency
306 _lambda = self.C / self.frequency
303 307
304 vmax = self.getFmax() * _lambda/2
308 vmax = self.getFmax() * _lambda / 2
305 309
306 310 return vmax
307 311
@@ -366,7 +370,8 class JROData(GenericData):
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")
@@ -378,9 +383,10 class JROData(GenericData):
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 #data es un numpy array de 2 dmensiones (canales, alturas)
389 # data es un numpy array de 2 dmensiones (canales, alturas)
384 390 data = None
385 391
386 392 def __init__(self):
@@ -428,17 +434,17 class Voltage(JROData):
428 434
429 435 self.blocksize = None
430 436
431 self.flagDecodeData = False #asumo q la data no esta decodificada
437 self.flagDecodeData = False # asumo q la data no esta decodificada
432 438
433 self.flagDeflipData = False #asumo q la data no esta sin flip
439 self.flagDeflipData = False # asumo q la data no esta sin flip
434 440
435 441 self.flagShiftFFT = False
436 442
437 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
443 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
438 444
439 445 self.profileIndex = 0
440 446
441 def getNoisebyHildebrand(self, channel = None):
447 def getNoisebyHildebrand(self, channel=None):
442 448 """
443 449 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
444 450
@@ -460,19 +466,19 class Voltage(JROData):
460 466 if nChannels == 1:
461 467 daux = power[:].real
462 468 else:
463 daux = power[thisChannel,:].real
469 daux = power[thisChannel, :].real
464 470 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
465 471
466 472 return noise
467 473
468 def getNoise(self, type = 1, channel = None):
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 def getPower(self, channel = None):
481 def getPower(self, channel=None):
476 482
477 483 if channel != None:
478 484 data = self.data[channel]
@@ -480,7 +486,7 class Voltage(JROData):
480 486 data = self.data
481 487
482 488 power = data * numpy.conjugate(data)
483 powerdB = 10*numpy.log10(power.real)
489 powerdB = 10 * numpy.log10(power.real)
484 490 powerdB = numpy.squeeze(powerdB)
485 491
486 492 return powerdB
@@ -494,18 +500,19 class Voltage(JROData):
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 #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
506 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
500 507 data_spc = None
501 508
502 #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
509 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
503 510 data_cspc = None
504 511
505 #data dc es un numpy array de 2 dmensiones (canales, alturas)
512 # data dc es un numpy array de 2 dmensiones (canales, alturas)
506 513 data_dc = None
507 514
508 #data power
515 # data power
509 516 data_pwr = None
510 517
511 518 nFFTPoints = None
@@ -516,9 +523,9 class Spectra(JROData):
516 523
517 524 nIncohInt = None
518 525
519 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
526 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
520 527
521 nCohInt = None #se requiere para determinar el valor de timeInterval
528 nCohInt = None # se requiere para determinar el valor de timeInterval
522 529
523 530 ippFactor = None
524 531
@@ -573,9 +580,9 class Spectra(JROData):
573 580
574 581 self.wavelength = None
575 582
576 self.flagDecodeData = False #asumo q la data no esta decodificada
583 self.flagDecodeData = False # asumo q la data no esta decodificada
577 584
578 self.flagDeflipData = False #asumo q la data no esta sin flip
585 self.flagDeflipData = False # asumo q la data no esta sin flip
579 586
580 587 self.flagShiftFFT = False
581 588
@@ -587,7 +594,6 class Spectra(JROData):
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
@@ -599,7 +605,8 class Spectra(JROData):
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
@@ -607,36 +614,45 class Spectra(JROData):
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 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
618 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
626 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
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 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
625 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
635 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
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 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
632 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
644 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
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 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
639 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2
653 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
654 velrange = deltav * (numpy.arange(self.nFFTPoints +
655 extrapoints) - self.nFFTPoints / 2.) # - deltav/2
640 656
641 657 return velrange
642 658
@@ -655,7 +671,8 class Spectra(JROData):
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
@@ -682,11 +699,11 class Spectra(JROData):
682 699 def getPower(self):
683 700
684 701 factor = self.normFactor
685 z = self.data_spc/factor
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 return 10*numpy.log10(avg)
706 return 10 * numpy.log10(avg)
690 707
691 708 def getCoherence(self, pairsList=None, phase=False):
692 709
@@ -697,17 +714,19 class Spectra(JROData):
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)
701 pairsIndexList.append(self.pairsList.index(pair))
717 raise ValueError, "Pair %s is not in dataOut.pairsList" % (
718 pair)
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 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
726 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
708 727 if phase:
709 728 data = numpy.arctan2(avgcoherenceComplex.imag,
710 avgcoherenceComplex.real)*180/numpy.pi
729 avgcoherenceComplex.real) * 180 / numpy.pi
711 730 else:
712 731 data = numpy.abs(avgcoherenceComplex)
713 732
@@ -722,12 +741,16 class Spectra(JROData):
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
@@ -790,7 +813,7 class SpectraHeis(Spectra):
790 813 if self.flagDecodeData:
791 814 pwcode = numpy.sum(self.code[0]**2)
792 815
793 normFactor = self.nIncohInt*self.nCohInt*pwcode
816 normFactor = self.nIncohInt * self.nCohInt * pwcode
794 817
795 818 return normFactor
796 819
@@ -803,6 +826,7 class SpectraHeis(Spectra):
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
@@ -831,14 +855,13 class Fits(JROData):
831 855
832 856 windowOfFilter = 1
833 857
834 #Speed of ligth
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"
@@ -879,11 +902,10 class Fits(JROData):
879 902 # self.comments = ''
880 903 #
881 904
882
883 905 def getltctime(self):
884 906
885 907 if self.useLocalTime:
886 return self.utctime - self.timeZone*60
908 return self.utctime - self.timeZone * 60
887 909
888 910 return self.utctime
889 911
@@ -921,7 +943,7 class Fits(JROData):
921 943
922 944 return range(self.nChannels)
923 945
924 def getNoise(self, type = 1):
946 def getNoise(self, type=1):
925 947
926 948 #noise = numpy.zeros(self.nChannels)
927 949
@@ -945,7 +967,8 class Fits(JROData):
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")
@@ -984,7 +1007,6 class Correlation(JROData):
984 1007
985 1008 nAvg = None
986 1009
987
988 1010 def __init__(self):
989 1011 '''
990 1012 Constructor
@@ -1019,9 +1041,9 class Correlation(JROData):
1019 1041
1020 1042 self.blocksize = None
1021 1043
1022 self.flagDecodeData = False #asumo q la data no esta decodificada
1044 self.flagDecodeData = False # asumo q la data no esta decodificada
1023 1045
1024 self.flagDeflipData = False #asumo q la data no esta sin flip
1046 self.flagDeflipData = False # asumo q la data no esta sin flip
1025 1047
1026 1048 self.pairsList = None
1027 1049
@@ -1031,48 +1053,50 class Correlation(JROData):
1031 1053
1032 1054 return self.pairsList
1033 1055
1034 def getNoise(self, mode = 2):
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 jspectra0 = self.data_corr[:,:,indR,:]
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 freq_dc = jspectra.shape[1]/2
1046 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1067 freq_dc = jspectra.shape[1] / 2
1068 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
1047 1069
1048 if ind_vel[0]<0:
1049 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1070 if ind_vel[0] < 0:
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 vel = numpy.array([-2,-1,1,2])
1057 xx = numpy.zeros([4,4])
1079 vel = numpy.array([-2, -1, 1, 2])
1080 xx = numpy.zeros([4, 4])
1058 1081
1059 1082 for fil in range(4):
1060 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1083 xx[fil, :] = vel[fil]**numpy.asarray(range(4))
1061 1084
1062 1085 xx_inv = numpy.linalg.inv(xx)
1063 xx_aux = xx_inv[0,:]
1086 xx_aux = xx_inv[0, :]
1064 1087
1065 1088 for ich in range(num_chan):
1066 yy = jspectra[ich,ind_vel,:]
1067 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1089 yy = jspectra[ich, ind_vel, :]
1090 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
1068 1091
1069 junkid = jspectra[ich,freq_dc,:]<=0
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 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1099 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
1076 1100
1077 1101 return noise
1078 1102
@@ -1093,7 +1117,7 class Correlation(JROData):
1093 1117 chan0 = pairsList[l][0]
1094 1118 chan1 = pairsList[l][1]
1095 1119
1096 #Obteniendo pares de Autocorrelacion
1120 # Obteniendo pares de Autocorrelacion
1097 1121 if chan0 == chan1:
1098 1122 acf_pairs.append(chan0)
1099 1123 acf_ind.append(l)
@@ -1109,7 +1133,7 class Correlation(JROData):
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 normFactor = numpy.zeros((self.nPairs,self.nHeights))
1136 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1113 1137
1114 1138 for p in range(self.nPairs):
1115 1139 pair = self.pairsList[p]
@@ -1117,69 +1141,69 class Correlation(JROData):
1117 1141 ch0 = pair[0]
1118 1142 ch1 = pair[1]
1119 1143
1120 ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1)
1121 ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1)
1122 normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max)
1144 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1145 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
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 experimentInfo = None #Information about the experiment
1156 experimentInfo = None # Information about the experiment
1132 1157
1133 #Information from previous data
1158 # Information from previous data
1134 1159
1135 inputUnit = None #Type of data to be processed
1160 inputUnit = None # Type of data to be processed
1136 1161
1137 operation = None #Type of operation to parametrize
1162 operation = None # Type of operation to parametrize
1138 1163
1139 #normFactor = None #Normalization Factor
1164 # normFactor = None #Normalization Factor
1140 1165
1141 groupList = None #List of Pairs, Groups, etc
1166 groupList = None # List of Pairs, Groups, etc
1142 1167
1143 #Parameters
1168 # Parameters
1144 1169
1145 data_param = None #Parameters obtained
1170 data_param = None # Parameters obtained
1146 1171
1147 data_pre = None #Data Pre Parametrization
1172 data_pre = None # Data Pre Parametrization
1148 1173
1149 data_SNR = None #Signal to Noise Ratio
1174 data_SNR = None # Signal to Noise Ratio
1150 1175
1151 1176 # heightRange = None #Heights
1152 1177
1153 abscissaList = None #Abscissa, can be velocities, lags or time
1178 abscissaList = None # Abscissa, can be velocities, lags or time
1154 1179
1155 1180 # noise = None #Noise Potency
1156 1181
1157 utctimeInit = None #Initial UTC time
1182 utctimeInit = None # Initial UTC time
1158 1183
1159 paramInterval = None #Time interval to calculate Parameters in seconds
1184 paramInterval = None # Time interval to calculate Parameters in seconds
1160 1185
1161 1186 useLocalTime = True
1162 1187
1163 #Fitting
1188 # Fitting
1164 1189
1165 data_error = None #Error of the estimation
1190 data_error = None # Error of the estimation
1166 1191
1167 1192 constants = None
1168 1193
1169 1194 library = None
1170 1195
1171 #Output signal
1196 # Output signal
1172 1197
1173 outputInterval = None #Time interval to calculate output signal in seconds
1198 outputInterval = None # Time interval to calculate output signal in seconds
1174 1199
1175 data_output = None #Out signal
1200 data_output = None # Out signal
1176 1201
1177 1202 nAvg = None
1178 1203
1179 1204 noise_estimation = None
1180
1181 GauSPC = None #Fit gaussian SPC
1182 1205
1206 GauSPC = None # Fit gaussian SPC
1183 1207
1184 1208 def __init__(self):
1185 1209 '''
@@ -1196,7 +1220,7 class Parameters(Spectra):
1196 1220 datatime = []
1197 1221
1198 1222 if self.useLocalTime:
1199 time1 = self.utctimeInit - self.timeZone*60
1223 time1 = self.utctimeInit - self.timeZone * 60
1200 1224 else:
1201 1225 time1 = self.utctimeInit
1202 1226
Show file before | Show file after | Hide comments
@@ -15,7 +15,6 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:
@@ -34,6 +33,7 try:
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
@@ -63,7 +63,7 class DigitalRFReader(ProcessingUnit):
63 63
64 64 def __getCurrentSecond(self):
65 65
66 return self.__thisUnixSample/self.__sample_rate
66 return self.__thisUnixSample / self.__sample_rate
67 67
68 68 thisSecond = property(__getCurrentSecond, "I'm the 'thisSecond' property.")
69 69
@@ -71,33 +71,35 class DigitalRFReader(ProcessingUnit):
71 71 '''
72 72 In this method will be initialized every parameter of dataOut object (header, no data)
73 73 '''
74 ippSeconds = 1.0*self.__nSamples/self.__sample_rate
74 ippSeconds = 1.0 * self.__nSamples / self.__sample_rate
75
76 nProfiles = 1.0 / ippSeconds # Number of profiles in one second
75 77
76 nProfiles = 1.0/ippSeconds # Number of profiles in one second
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 txA=0,
83 txB=0,
84 nWindows=1,
85 nHeights=self.__nSamples,
86 firstHeight=self.__firstHeigth,
87 deltaHeight=self.__deltaHeigth,
88 codeType=self.__codeType,
89 nCode=self.__nCode, nBaud=self.__nBaud,
90 code = self.__code)
91
83 txA=0,
84 txB=0,
85 nWindows=1,
86 nHeights=self.__nSamples,
87 firstHeight=self.__firstHeigth,
88 deltaHeight=self.__deltaHeigth,
89 codeType=self.__codeType,
90 nCode=self.__nCode, nBaud=self.__nBaud,
91 code=self.__code)
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
@@ -108,7 +110,9 class DigitalRFReader(ProcessingUnit):
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
@@ -124,25 +128,29 class DigitalRFReader(ProcessingUnit):
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
146 154
147 155 self.dataOut.ippSeconds = ippSeconds
148 156
@@ -159,7 +167,8 class DigitalRFReader(ProcessingUnit):
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
@@ -179,7 +188,8 class DigitalRFReader(ProcessingUnit):
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)
@@ -194,7 +204,7 class DigitalRFReader(ProcessingUnit):
194 204
195 205 thisDatetime = startDatetime
196 206
197 while(thisDatetime<=endDatatime):
207 while(thisDatetime <= endDatatime):
198 208
199 209 thisDate = thisDatetime.date()
200 210
@@ -209,18 +219,23 class DigitalRFReader(ProcessingUnit):
209 219
210 220 return dateList
211 221
212 def setup(self, path = None,
213 startDate = None,
214 endDate = None,
215 startTime = datetime.time(0,0,0),
216 endTime = datetime.time(23,59,59),
217 channelList = None,
218 nSamples = None,
219 online = False,
220 delay = 60,
221 buffer_size = 1024,
222 ippKm=None,
223 **kwargs):
222 def setup(self, path=None,
223 startDate=None,
224 endDate=None,
225 startTime=datetime.time(0, 0, 0),
226 endTime=datetime.time(23, 59, 59),
227 channelList=None,
228 nSamples=None,
229 online=False,
230 delay=60,
231 buffer_size=1024,
232 ippKm=None,
233 nCohInt=1,
234 nCode=1,
235 nBaud=1,
236 flagDecodeData=False,
237 code=numpy.ones((1, 1), dtype=numpy.int),
238 **kwargs):
224 239 '''
225 240 In this method we should set all initial parameters.
226 241
@@ -236,37 +251,43 class DigitalRFReader(ProcessingUnit):
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 raise ValueError, "[Reading] Directory %s does not exist" %path
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 raise ValueError, "[Reading] Directory %s does not have any files" %path
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']
@@ -275,7 +296,6 class DigitalRFReader(ProcessingUnit):
275 296 self.dtype = cPickle.loads(self.fixed_metadata_dict['dtype'])
276 297 except:
277 298 pass
278
279 299
280 300 self.__frequency = None
281 301
@@ -283,12 +303,11 class DigitalRFReader(ProcessingUnit):
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:
@@ -296,10 +315,6 class DigitalRFReader(ProcessingUnit):
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']
@@ -307,8 +322,7 class DigitalRFReader(ProcessingUnit):
307 322 code = self.__radarControllerHeader['code']
308 323 except:
309 324 pass
310
311
325
312 326 if not ippKm:
313 327 try:
314 328 # seconds to km
@@ -322,42 +336,46 class DigitalRFReader(ProcessingUnit):
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 startUTCSecond = start_index/self.__sample_rate
351 startUTCSecond = start_index / self.__sample_rate
335 352
336 if start_index > startUTCSecond*self.__sample_rate:
337 startUTCSecond = start_index/self.__sample_rate
353 if start_index > startUTCSecond * self.__sample_rate:
354 startUTCSecond = start_index / self.__sample_rate
338 355
339 356 if not endUTCSecond:
340 endUTCSecond = end_index/self.__sample_rate
357 endUTCSecond = end_index / self.__sample_rate
341 358
342 if end_index < endUTCSecond*self.__sample_rate:
343 endUTCSecond = end_index/self.__sample_rate
359 if end_index < endUTCSecond * self.__sample_rate:
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 nSamples = int(ippKm / (1e6*0.15/self.__sample_rate))
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 thisChannelName = channelNameList[thisIndexChannel]
353 start_index, end_index = self.digitalReadObj.get_bounds(thisChannelName)
369 thisChannelName = channelNameList[thisIndexChannel]
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 self.i= 0
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
@@ -369,36 +387,44 class DigitalRFReader(ProcessingUnit):
369 387 self.__channelNameList = channelNameListFiltered
370 388 self.__channelBoundList = channelBoundList
371 389 self.__nSamples = nSamples
372 self.__samples_to_read = long(nSamples) # FIJO: AHORA 40
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 # self.__thisUnixSample = int(endUTCSecond*self.__sample_rate - 4*self.__samples_to_read)
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 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)
394 )
395
396 print "[Reading] Starting process from %s to %s" %(datetime.datetime.utcfromtimestamp(startUTCSecond - self.__timezone),
397 datetime.datetime.utcfromtimestamp(endUTCSecond - self.__timezone)
398 )
413 print "[Reading] Digital RF Data was found from %s to %s " % (
414 datetime.datetime.utcfromtimestamp(
415 self.__startUTCSecond - self.__timezone),
416 datetime.datetime.utcfromtimestamp(
417 self.__endUTCSecond - self.__timezone)
418 )
419
420 print "[Reading] Starting process from %s to %s" % (datetime.datetime.utcfromtimestamp(startUTCSecond - self.__timezone),
421 datetime.datetime.utcfromtimestamp(
422 endUTCSecond - self.__timezone)
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]" %(
@@ -413,18 +439,21 class DigitalRFReader(ProcessingUnit):
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 if start_index > self.__startUTCSecond*self.__sample_rate:
419 self.__startUTCSecond = 1.0*start_index/self.__sample_rate
445 if start_index > self.__startUTCSecond * self.__sample_rate:
446 self.__startUTCSecond = 1.0 * start_index / self.__sample_rate
420 447
421 if end_index > self.__endUTCSecond*self.__sample_rate:
422 self.__endUTCSecond = 1.0*end_index/self.__sample_rate
448 if end_index > self.__endUTCSecond * self.__sample_rate:
449 self.__endUTCSecond = 1.0 * end_index / self.__sample_rate
423 450 print
424 print "[Reading] New timerange found [%s, %s] " %(
425 datetime.datetime.utcfromtimestamp(self.__startUTCSecond - self.__timezone),
426 datetime.datetime.utcfromtimestamp(self.__endUTCSecond - self.__timezone)
427 )
451 print "[Reading] New timerange found [%s, %s] " % (
452 datetime.datetime.utcfromtimestamp(
453 self.__startUTCSecond - self.__timezone),
454 datetime.datetime.utcfromtimestamp(
455 self.__endUTCSecond - self.__timezone)
456 )
428 457
429 458 return True
430 459
@@ -434,12 +463,14 class DigitalRFReader(ProcessingUnit):
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 def __readNextBlock(self, seconds=30, volt_scale = 1):
473 def __readNextBlock(self, seconds=30, volt_scale=1):
443 474 '''
444 475 '''
445 476
@@ -447,21 +478,21 class DigitalRFReader(ProcessingUnit):
447 478 self.__flagDiscontinuousBlock = False
448 479 self.__thisUnixSample += self.__samples_to_read
449 480
450 if self.__thisUnixSample + 2*self.__samples_to_read > self.__endUTCSecond*self.__sample_rate:
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 if self.__thisUnixSample + 2*self.__samples_to_read > self.__endUTCSecond*self.__sample_rate:
488 if self.__thisUnixSample + 2 * self.__samples_to_read > self.__endUTCSecond * self.__sample_rate:
458 489 return False
459 self.__thisUnixSample -= self.__samples_to_read
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()
@@ -469,47 +500,48 class DigitalRFReader(ProcessingUnit):
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 #read next profile
510 # read next profile
478 511 self.__flagDiscontinuousBlock = True
479 print "[Reading] %s" %datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone), e
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 print "[Reading] %s: Too few samples were found, just %d/%d samples" %(datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
485 result.shape[0],
486 self.__samples_to_read)
517 print "[Reading] %s: Too few samples were found, just %d/%d samples" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
518 result.shape[0],
519 self.__samples_to_read)
487 520 break
488 521
489 self.__data_buffer[indexSubchannel,:] = result*volt_scale
522 self.__data_buffer[indexSubchannel, :] = result * volt_scale
490 523
491 524 indexChannel += 1
492 525
493 526 dataOk = True
494
495 self.__utctime = self.__thisUnixSample/self.__sample_rate
527
528 self.__utctime = self.__thisUnixSample / self.__sample_rate
496 529
497 530 if not dataOk:
498 531 return False
499 532
500 print "[Reading] %s: %d samples <> %f sec" %(datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
501 self.__samples_to_read,
502 self.__timeInterval)
533 print "[Reading] %s: %d samples <> %f sec" % (datetime.datetime.utcfromtimestamp(self.thisSecond - self.__timezone),
534 self.__samples_to_read,
535 self.__timeInterval)
503 536
504 537 self.__bufferIndex = 0
505 538
506 539 return True
507 540
508 541 def __isBufferEmpty(self):
509 return self.__bufferIndex > self.__samples_to_read - self.__nSamples #40960 - 40
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
@@ -535,7 +567,7 class DigitalRFReader(ProcessingUnit):
535 567 while True:
536 568 if self.__readNextBlock():
537 569 break
538 if self.__thisUnixSample > self.__endUTCSecond*self.__sample_rate:
570 if self.__thisUnixSample > self.__endUTCSecond * self.__sample_rate:
539 571 return False
540 572
541 573 if self.__flagDiscontinuousBlock:
@@ -549,11 +581,13 class DigitalRFReader(ProcessingUnit):
549 581 if err_counter > nTries:
550 582 return False
551 583
552 print '[Reading] waiting %d seconds to read a new block' %seconds
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
@@ -583,12 +617,11 class DigitalRFReader(ProcessingUnit):
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
@@ -596,6 +629,7 class DigitalRFReader(ProcessingUnit):
596 629
597 630 return
598 631
632
599 633 class DigitalRFWriter(Operation):
600 634 '''
601 635 classdocs
@@ -607,7 +641,7 class DigitalRFWriter(Operation):
607 641 '''
608 642 Operation.__init__(self, **kwargs)
609 643 self.metadata_dict = {}
610 self.dataOut = None
644 self.dataOut = None
611 645 self.dtype = None
612 646
613 647 def setHeader(self):
@@ -624,7 +658,7 class DigitalRFWriter(Operation):
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'):
@@ -636,7 +670,7 class DigitalRFWriter(Operation):
636 670 self.setHeader()
637 671 self.__ippSeconds = dataOut.ippSeconds
638 672 self.__deltaH = dataOut.getDeltaH()
639 self.__sample_rate = 1e6*0.15/self.__deltaH
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]
@@ -644,16 +678,18 class DigitalRFWriter(Operation):
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
@@ -663,45 +699,47 class DigitalRFWriter(Operation):
663 699 marching_periods = False
664 700
665 701 self.digitalWriteObj = digital_rf.DigitalRFWriter(path, self.__dtype, dirCadence,
666 fileCadence, start_global_index,
667 sample_rate_numerator, sample_rate_denominator, uuid, compression_level, checksum,
668 is_complex, num_subchannels, is_continuous, marching_periods)
669
702 fileCadence, start_global_index,
703 sample_rate_numerator, sample_rate_denominator, uuid, compression_level, checksum,
704 is_complex, num_subchannels, is_continuous, marching_periods)
705
670 706 metadata_dir = os.path.join(path, 'metadata')
671 707 os.system('mkdir %s' % (metadata_dir))
672
673 self.digitalMetadataWriteObj = digital_rf.DigitalMetadataWriter(metadata_dir, dirCadence, 1, ##236, file_cadence_millisecs / 1000
674 sample_rate_numerator, sample_rate_denominator,
675 metadataFile)
676 708
709 self.digitalMetadataWriteObj = digital_rf.DigitalMetadataWriter(metadata_dir, dirCadence, 1, # 236, file_cadence_millisecs / 1000
710 sample_rate_numerator, sample_rate_denominator,
711 metadataFile)
677 712
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
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()
722
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:
@@ -710,9 +748,9 class DigitalRFWriter(Operation):
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
@@ -722,14 +760,15 class DigitalRFWriter(Operation):
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):
@@ -739,7 +778,7 class DigitalRFWriter(Operation):
739 778 self.digitalWriteObj.close()
740 779 except:
741 780 pass
742
781
743 782 # raise
744 783 if __name__ == '__main__':
745 784
@@ -748,4 +787,4 if __name__ == '__main__':
748 787 while True:
749 788 readObj.run(path='/home/jchavez/jicamarca/mocked_data/')
750 789 # readObj.printInfo()
751 # readObj.printNumberOfBlock()
790 # readObj.printNumberOfBlock()
Show file before | Show file after | Hide comments
@@ -616,7 +616,6 class Decoder(Operation):
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
@@ -666,7 +665,6 class Decoder(Operation):
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
General Comments 0
You need to be logged in to leave comments. Login now