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