##// END OF EJS Templates
schain
RSS
Bug Fixed: Actualizacion del valor estimado del ruido
Daniel Valdez -
r452:b91c52c913af
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,636 +1,636
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 os, sys
8 8 import copy
9 9 import numpy
10 10 import datetime
11 11
12 12 from jroheaderIO import SystemHeader, RadarControllerHeader
13 13
14 14
15 15 def hildebrand_sekhon(data, navg):
16 16
17 17 data = data.copy()
18 18
19 19 sortdata = numpy.sort(data,axis=None)
20 20 lenOfData = len(sortdata)
21 21 nums_min = lenOfData/10
22 22
23 23 if (lenOfData/10) > 2:
24 24 nums_min = lenOfData/10
25 25 else:
26 26 nums_min = 2
27 27
28 28 sump = 0.
29 29
30 30 sumq = 0.
31 31
32 32 j = 0
33 33
34 34 cont = 1
35 35
36 36 while((cont==1)and(j<lenOfData)):
37 37
38 38 sump += sortdata[j]
39 39
40 40 sumq += sortdata[j]**2
41 41
42 42 j += 1
43 43
44 44 if j > nums_min:
45 45 rtest = float(j)/(j-1) + 1.0/navg
46 46 if ((sumq*j) > (rtest*sump**2)):
47 47 j = j - 1
48 48 sump = sump - sortdata[j]
49 49 sumq = sumq - sortdata[j]**2
50 50 cont = 0
51 51
52 52 lnoise = sump /j
53 53 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
54 54 return lnoise
55 55
56 56 class JROData:
57 57
58 58 # m_BasicHeader = BasicHeader()
59 59 # m_ProcessingHeader = ProcessingHeader()
60 60
61 61 systemHeaderObj = SystemHeader()
62 62
63 63 radarControllerHeaderObj = RadarControllerHeader()
64 64
65 65 # data = None
66 66
67 67 type = None
68 68
69 69 dtype = None
70 70
71 71 # nChannels = None
72 72
73 73 # nHeights = None
74 74
75 75 nProfiles = None
76 76
77 77 heightList = None
78 78
79 79 channelList = None
80 80
81 81 flagNoData = True
82 82
83 83 flagTimeBlock = False
84 84
85 85 useLocalTime = False
86 86
87 87 utctime = None
88 88
89 89 timeZone = None
90 90
91 91 dstFlag = None
92 92
93 93 errorCount = None
94 94
95 95 blocksize = None
96 96
97 97 nCode = None
98 98
99 99 nBaud = None
100 100
101 101 code = None
102 102
103 103 flagDecodeData = False #asumo q la data no esta decodificada
104 104
105 105 flagDeflipData = False #asumo q la data no esta sin flip
106 106
107 107 flagShiftFFT = False
108 108
109 109 ippSeconds = None
110 110
111 111 timeInterval = None
112 112
113 113 nCohInt = None
114 114
115 115 noise = None
116 116
117 117 windowOfFilter = 1
118 118
119 119 #Speed of ligth
120 120 C = 3e8
121 121
122 122 frequency = 49.92e6
123 123
124 124 realtime = False
125 125
126 126 def __init__(self):
127 127
128 128 raise ValueError, "This class has not been implemented"
129 129
130 130 def copy(self, inputObj=None):
131 131
132 132 if inputObj == None:
133 133 return copy.deepcopy(self)
134 134
135 135 for key in inputObj.__dict__.keys():
136 136 self.__dict__[key] = inputObj.__dict__[key]
137 137
138 138 def deepcopy(self):
139 139
140 140 return copy.deepcopy(self)
141 141
142 142 def isEmpty(self):
143 143
144 144 return self.flagNoData
145 145
146 146 def getNoise(self):
147 147
148 148 raise ValueError, "Not implemented"
149 149
150 150 def getNChannels(self):
151 151
152 152 return len(self.channelList)
153 153
154 154 def getChannelIndexList(self):
155 155
156 156 return range(self.nChannels)
157 157
158 158 def getNHeights(self):
159 159
160 160 return len(self.heightList)
161 161
162 162 def getHeiRange(self, extrapoints=0):
163 163
164 164 heis = self.heightList
165 165 # deltah = self.heightList[1] - self.heightList[0]
166 166 #
167 167 # heis.append(self.heightList[-1])
168 168
169 169 return heis
170 170
171 171 def getltctime(self):
172 172
173 173 if self.useLocalTime:
174 174 return self.utctime - self.timeZone*60
175 175
176 176 return self.utctime
177 177
178 178 def getDatatime(self):
179 179
180 180 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
181 181 return datatime
182 182
183 183 def getTimeRange(self):
184 184
185 185 datatime = []
186 186
187 187 datatime.append(self.ltctime)
188 188 datatime.append(self.ltctime + self.timeInterval)
189 189
190 190 datatime = numpy.array(datatime)
191 191
192 192 return datatime
193 193
194 194 def getFmax(self):
195 195
196 196 PRF = 1./(self.ippSeconds * self.nCohInt)
197 197
198 198 fmax = PRF/2.
199 199
200 200 return fmax
201 201
202 202 def getVmax(self):
203 203
204 204 _lambda = self.C/self.frequency
205 205
206 206 vmax = self.getFmax() * _lambda
207 207
208 208 return vmax
209 209
210 210 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
211 211 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
212 212 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
213 213 noise = property(getNoise, "I'm the 'nHeights' property.")
214 214 datatime = property(getDatatime, "I'm the 'datatime' property")
215 215 ltctime = property(getltctime, "I'm the 'ltctime' property")
216 216
217 217 class Voltage(JROData):
218 218
219 219 #data es un numpy array de 2 dmensiones (canales, alturas)
220 220 data = None
221 221
222 222 def __init__(self):
223 223 '''
224 224 Constructor
225 225 '''
226 226
227 227 self.radarControllerHeaderObj = RadarControllerHeader()
228 228
229 229 self.systemHeaderObj = SystemHeader()
230 230
231 231 self.type = "Voltage"
232 232
233 233 self.data = None
234 234
235 235 self.dtype = None
236 236
237 237 # self.nChannels = 0
238 238
239 239 # self.nHeights = 0
240 240
241 241 self.nProfiles = None
242 242
243 243 self.heightList = None
244 244
245 245 self.channelList = None
246 246
247 247 # self.channelIndexList = None
248 248
249 249 self.flagNoData = True
250 250
251 251 self.flagTimeBlock = False
252 252
253 253 self.utctime = None
254 254
255 255 self.timeZone = None
256 256
257 257 self.dstFlag = None
258 258
259 259 self.errorCount = None
260 260
261 261 self.nCohInt = None
262 262
263 263 self.blocksize = None
264 264
265 265 self.flagDecodeData = False #asumo q la data no esta decodificada
266 266
267 267 self.flagDeflipData = False #asumo q la data no esta sin flip
268 268
269 269 self.flagShiftFFT = False
270 270
271 271
272 272 def getNoisebyHildebrand(self):
273 273 """
274 274 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
275 275
276 276 Return:
277 277 noiselevel
278 278 """
279 279
280 280 for channel in range(self.nChannels):
281 281 daux = self.data_spc[channel,:,:]
282 282 self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt)
283 283
284 284 return self.noise
285 285
286 286 def getNoise(self, type = 1):
287 287
288 288 self.noise = numpy.zeros(self.nChannels)
289 289
290 290 if type == 1:
291 291 noise = self.getNoisebyHildebrand()
292 292
293 293 return 10*numpy.log10(noise)
294 294
295 295 class Spectra(JROData):
296 296
297 297 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
298 298 data_spc = None
299 299
300 300 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
301 301 data_cspc = None
302 302
303 303 #data es un numpy array de 2 dmensiones (canales, alturas)
304 304 data_dc = None
305 305
306 306 nFFTPoints = None
307 307
308 308 nPairs = None
309 309
310 310 pairsList = None
311 311
312 312 nIncohInt = None
313 313
314 314 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
315 315
316 316 nCohInt = None #se requiere para determinar el valor de timeInterval
317 317
318 318 ippFactor = None
319 319
320 320 def __init__(self):
321 321 '''
322 322 Constructor
323 323 '''
324 324
325 325 self.radarControllerHeaderObj = RadarControllerHeader()
326 326
327 327 self.systemHeaderObj = SystemHeader()
328 328
329 329 self.type = "Spectra"
330 330
331 331 # self.data = None
332 332
333 333 self.dtype = None
334 334
335 335 # self.nChannels = 0
336 336
337 337 # self.nHeights = 0
338 338
339 339 self.nProfiles = None
340 340
341 341 self.heightList = None
342 342
343 343 self.channelList = None
344 344
345 345 # self.channelIndexList = None
346 346
347 347 self.flagNoData = True
348 348
349 349 self.flagTimeBlock = False
350 350
351 351 self.utctime = None
352 352
353 353 self.nCohInt = None
354 354
355 355 self.nIncohInt = None
356 356
357 357 self.blocksize = None
358 358
359 359 self.nFFTPoints = None
360 360
361 361 self.wavelength = None
362 362
363 363 self.flagDecodeData = False #asumo q la data no esta decodificada
364 364
365 365 self.flagDeflipData = False #asumo q la data no esta sin flip
366 366
367 367 self.flagShiftFFT = False
368 368
369 369 self.ippFactor = 1
370 370
371 371 self.noise = None
372 372
373 373 def getNoisebyHildebrand(self):
374 374 """
375 375 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
376 376
377 377 Return:
378 378 noiselevel
379 379 """
380
380 noise = numpy.zeros(self.nChannels)
381 381 for channel in range(self.nChannels):
382 382 daux = self.data_spc[channel,:,:]
383 self.noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
383 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
384 384
385 return self.noise
385 return noise
386 386
387 387 def getNoise(self):
388 if self.noise == None:
389 self.noise = numpy.zeros(self.nChannels)
390 self.noise = self.getNoisebyHildebrand()
391
388 if self.noise != None:
392 389 return self.noise
390 else:
391 noise = self.getNoisebyHildebrand()
392 return noise
393 393
394 394
395 395 def getFreqRange(self, extrapoints=0):
396 396
397 397 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
398 398 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
399 399
400 400 return freqrange
401 401
402 402 def getVelRange(self, extrapoints=0):
403 403
404 404 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
405 405 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
406 406
407 407 return velrange
408 408
409 409 def getNPairs(self):
410 410
411 411 return len(self.pairsList)
412 412
413 413 def getPairsIndexList(self):
414 414
415 415 return range(self.nPairs)
416 416
417 417 def getNormFactor(self):
418 418 pwcode = 1
419 419 if self.flagDecodeData:
420 420 pwcode = numpy.sum(self.code[0]**2)
421 421 normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode
422 422
423 423 return normFactor
424 424
425 425 def getFlagCspc(self):
426 426
427 427 if self.data_cspc == None:
428 428 return True
429 429
430 430 return False
431 431
432 432 def getFlagDc(self):
433 433
434 434 if self.data_dc == None:
435 435 return True
436 436
437 437 return False
438 438
439 439 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
440 440 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
441 441 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
442 442 flag_cspc = property(getFlagCspc)
443 443 flag_dc = property(getFlagDc)
444 444
445 445 class SpectraHeis(JROData):
446 446
447 447 data_spc = None
448 448
449 449 data_cspc = None
450 450
451 451 data_dc = None
452 452
453 453 nFFTPoints = None
454 454
455 455 nPairs = None
456 456
457 457 pairsList = None
458 458
459 459 nIncohInt = None
460 460
461 461 def __init__(self):
462 462
463 463 self.radarControllerHeaderObj = RadarControllerHeader()
464 464
465 465 self.systemHeaderObj = SystemHeader()
466 466
467 467 self.type = "SpectraHeis"
468 468
469 469 self.dtype = None
470 470
471 471 # self.nChannels = 0
472 472
473 473 # self.nHeights = 0
474 474
475 475 self.nProfiles = None
476 476
477 477 self.heightList = None
478 478
479 479 self.channelList = None
480 480
481 481 # self.channelIndexList = None
482 482
483 483 self.flagNoData = True
484 484
485 485 self.flagTimeBlock = False
486 486
487 487 self.nPairs = 0
488 488
489 489 self.utctime = None
490 490
491 491 self.blocksize = None
492 492
493 493 class Fits:
494 494
495 495 heightList = None
496 496
497 497 channelList = None
498 498
499 499 flagNoData = True
500 500
501 501 flagTimeBlock = False
502 502
503 503 useLocalTime = False
504 504
505 505 utctime = None
506 506
507 507 timeZone = None
508 508
509 509 ippSeconds = None
510 510
511 511 timeInterval = None
512 512
513 513 nCohInt = None
514 514
515 515 nIncohInt = None
516 516
517 517 noise = None
518 518
519 519 windowOfFilter = 1
520 520
521 521 #Speed of ligth
522 522 C = 3e8
523 523
524 524 frequency = 49.92e6
525 525
526 526 realtime = False
527 527
528 528
529 529 def __init__(self):
530 530
531 531 self.type = "Fits"
532 532
533 533 self.nProfiles = None
534 534
535 535 self.heightList = None
536 536
537 537 self.channelList = None
538 538
539 539 # self.channelIndexList = None
540 540
541 541 self.flagNoData = True
542 542
543 543 self.utctime = None
544 544
545 545 self.nCohInt = None
546 546
547 547 self.nIncohInt = None
548 548
549 549 self.useLocalTime = True
550 550
551 551 # self.utctime = None
552 552 # self.timeZone = None
553 553 # self.ltctime = None
554 554 # self.timeInterval = None
555 555 # self.header = None
556 556 # self.data_header = None
557 557 # self.data = None
558 558 # self.datatime = None
559 559 # self.flagNoData = False
560 560 # self.expName = ''
561 561 # self.nChannels = None
562 562 # self.nSamples = None
563 563 # self.dataBlocksPerFile = None
564 564 # self.comments = ''
565 565 #
566 566
567 567
568 568 def getltctime(self):
569 569
570 570 if self.useLocalTime:
571 571 return self.utctime - self.timeZone*60
572 572
573 573 return self.utctime
574 574
575 575 def getDatatime(self):
576 576
577 577 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
578 578 return datatime
579 579
580 580 def getTimeRange(self):
581 581
582 582 datatime = []
583 583
584 584 datatime.append(self.ltctime)
585 585 datatime.append(self.ltctime + self.timeInterval)
586 586
587 587 datatime = numpy.array(datatime)
588 588
589 589 return datatime
590 590
591 591 def getHeiRange(self):
592 592
593 593 heis = self.heightList
594 594
595 595 return heis
596 596
597 597 def isEmpty(self):
598 598
599 599 return self.flagNoData
600 600
601 601 def getNHeights(self):
602 602
603 603 return len(self.heightList)
604 604
605 605 def getNChannels(self):
606 606
607 607 return len(self.channelList)
608 608
609 609 def getChannelIndexList(self):
610 610
611 611 return range(self.nChannels)
612 612
613 613 def getNoise(self, type = 1):
614 614
615 615 self.noise = numpy.zeros(self.nChannels)
616 616
617 617 if type == 1:
618 618 noise = self.getNoisebyHildebrand()
619 619
620 620 if type == 2:
621 621 noise = self.getNoisebySort()
622 622
623 623 if type == 3:
624 624 noise = self.getNoisebyWindow()
625 625
626 626 return noise
627 627
628 628 datatime = property(getDatatime, "I'm the 'datatime' property")
629 629 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
630 630 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
631 631 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
632 632 noise = property(getNoise, "I'm the 'nHeights' property.")
633 633 datatime = property(getDatatime, "I'm the 'datatime' property")
634 634 ltctime = property(getltctime, "I'm the 'ltctime' property")
635 635
636 636 ltctime = property(getltctime, "I'm the 'ltctime' property") No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now