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