##// END OF EJS Templates
schain
RSS
En jrodata.py se eliminan los metodos sorting_bruce, getNoisebySort, getNoisebyWindow....
Daniel Valdez -
r450:3d2d18302dc5
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,725 +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 def hildebrand_sekhon(data, navg):
15 """
16 This method is for the objective determination of de noise level in Doppler spectra. This
17 implementation technique is based on the fact that the standard deviation of the spectral
18 densities is equal to the mean spectral density for white Gaussian noise
19
20 Inputs:
21 Data : heights
22 navg : numbers of averages
23
24 Return:
25 -1 : any error
26 anoise : noise's level
27 """
28
29 dataflat = data.copy().reshape(-1)
30 dataflat.sort()
31 npts = dataflat.size #numbers of points of the data
32 npts_noise = 0.2*npts
33
34 if npts < 32:
35 print "error in noise - requires at least 32 points"
36 return -1.0
37
38 dataflat2 = numpy.power(dataflat,2)
39
40 cs = numpy.cumsum(dataflat)
41 cs2 = numpy.cumsum(dataflat2)
42
43 # data sorted in ascending order
44 nmin = int((npts + 7.)/8)
45 14
46 for i in range(nmin, npts):
47 s = cs[i]
48 s2 = cs2[i]
49 p = s / float(i);
50 p2 = p**2;
51 q = s2 / float(i) - p2;
52 leftc = p2;
53 rightc = q * float(navg);
54 R2 = leftc/rightc
55
56 # Signal detect: R2 < 1 (R2 = leftc/rightc)
57 if R2 < 1:
58 npts_noise = i
59 break
60
61
62 anoise = numpy.average(dataflat[0:npts_noise])
63
64 return anoise;
65
66 def sorting_bruce(data, navg):
15 def hildebrand_sekhon(data, navg):
67 16
68 17 data = data.copy()
69 18
70 sortdata = numpy.sort(data)
71 lenOfData = len(data)
19 sortdata = numpy.sort(data,axis=None)
20 lenOfData = len(sortdata)
72 21 nums_min = lenOfData/10
73 22
74 if (lenOfData/10) > 0:
23 if (lenOfData/10) > 2:
75 24 nums_min = lenOfData/10
76 25 else:
77 nums_min = 0
78
79 rtest = 1.0 + 1.0/navg
26 nums_min = 2
80 27
81 28 sump = 0.
82 29
83 30 sumq = 0.
84 31
85 32 j = 0
86 33
87 34 cont = 1
88 35
89 36 while((cont==1)and(j<lenOfData)):
90 37
91 38 sump += sortdata[j]
92 39
93 40 sumq += sortdata[j]**2
94 41
95 42 j += 1
96 43
97 44 if j > nums_min:
98 if ((sumq*j) <= (rtest*sump**2)):
99 lnoise = sump / j
100 else:
45 rtest = float(j)/(j-1) + 1.0/navg
46 if ((sumq*j) > (rtest*sump**2)):
101 47 j = j - 1
102 48 sump = sump - sortdata[j]
103 49 sumq = sumq - sortdata[j]**2
104 50 cont = 0
105 51
106 if j == nums_min:
107 52 lnoise = sump /j
108
53 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
109 54 return lnoise
110 55
111 56 class JROData:
112 57
113 58 # m_BasicHeader = BasicHeader()
114 59 # m_ProcessingHeader = ProcessingHeader()
115 60
116 61 systemHeaderObj = SystemHeader()
117 62
118 63 radarControllerHeaderObj = RadarControllerHeader()
119 64
120 65 # data = None
121 66
122 67 type = None
123 68
124 69 dtype = None
125 70
126 71 # nChannels = None
127 72
128 73 # nHeights = None
129 74
130 75 nProfiles = None
131 76
132 77 heightList = None
133 78
134 79 channelList = None
135 80
136 81 flagNoData = True
137 82
138 83 flagTimeBlock = False
139 84
140 85 useLocalTime = False
141 86
142 87 utctime = None
143 88
144 89 timeZone = None
145 90
146 91 dstFlag = None
147 92
148 93 errorCount = None
149 94
150 95 blocksize = None
151 96
152 97 nCode = None
153 98
154 99 nBaud = None
155 100
156 101 code = None
157 102
158 103 flagDecodeData = False #asumo q la data no esta decodificada
159 104
160 105 flagDeflipData = False #asumo q la data no esta sin flip
161 106
162 107 flagShiftFFT = False
163 108
164 109 ippSeconds = None
165 110
166 111 timeInterval = None
167 112
168 113 nCohInt = None
169 114
170 115 noise = None
171 116
172 117 windowOfFilter = 1
173 118
174 119 #Speed of ligth
175 120 C = 3e8
176 121
177 122 frequency = 49.92e6
178 123
179 124 realtime = False
180 125
181 126 def __init__(self):
182 127
183 128 raise ValueError, "This class has not been implemented"
184 129
185 130 def copy(self, inputObj=None):
186 131
187 132 if inputObj == None:
188 133 return copy.deepcopy(self)
189 134
190 135 for key in inputObj.__dict__.keys():
191 136 self.__dict__[key] = inputObj.__dict__[key]
192 137
193 138 def deepcopy(self):
194 139
195 140 return copy.deepcopy(self)
196 141
197 142 def isEmpty(self):
198 143
199 144 return self.flagNoData
200 145
201 146 def getNoise(self):
202 147
203 148 raise ValueError, "Not implemented"
204 149
205 150 def getNChannels(self):
206 151
207 152 return len(self.channelList)
208 153
209 154 def getChannelIndexList(self):
210 155
211 156 return range(self.nChannels)
212 157
213 158 def getNHeights(self):
214 159
215 160 return len(self.heightList)
216 161
217 162 def getHeiRange(self, extrapoints=0):
218 163
219 164 heis = self.heightList
220 165 # deltah = self.heightList[1] - self.heightList[0]
221 166 #
222 167 # heis.append(self.heightList[-1])
223 168
224 169 return heis
225 170
226 171 def getltctime(self):
227 172
228 173 if self.useLocalTime:
229 174 return self.utctime - self.timeZone*60
230 175
231 176 return self.utctime
232 177
233 178 def getDatatime(self):
234 179
235 180 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
236 181 return datatime
237 182
238 183 def getTimeRange(self):
239 184
240 185 datatime = []
241 186
242 187 datatime.append(self.ltctime)
243 188 datatime.append(self.ltctime + self.timeInterval)
244 189
245 190 datatime = numpy.array(datatime)
246 191
247 192 return datatime
248 193
249 194 def getFmax(self):
250 195
251 196 PRF = 1./(self.ippSeconds * self.nCohInt)
252 197
253 198 fmax = PRF/2.
254 199
255 200 return fmax
256 201
257 202 def getVmax(self):
258 203
259 204 _lambda = self.C/self.frequency
260 205
261 206 vmax = self.getFmax() * _lambda
262 207
263 208 return vmax
264 209
265 210 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
266 211 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
267 212 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
268 213 noise = property(getNoise, "I'm the 'nHeights' property.")
269 214 datatime = property(getDatatime, "I'm the 'datatime' property")
270 215 ltctime = property(getltctime, "I'm the 'ltctime' property")
271 216
272 217 class Voltage(JROData):
273 218
274 219 #data es un numpy array de 2 dmensiones (canales, alturas)
275 220 data = None
276 221
277 222 def __init__(self):
278 223 '''
279 224 Constructor
280 225 '''
281 226
282 227 self.radarControllerHeaderObj = RadarControllerHeader()
283 228
284 229 self.systemHeaderObj = SystemHeader()
285 230
286 231 self.type = "Voltage"
287 232
288 233 self.data = None
289 234
290 235 self.dtype = None
291 236
292 237 # self.nChannels = 0
293 238
294 239 # self.nHeights = 0
295 240
296 241 self.nProfiles = None
297 242
298 243 self.heightList = None
299 244
300 245 self.channelList = None
301 246
302 247 # self.channelIndexList = None
303 248
304 249 self.flagNoData = True
305 250
306 251 self.flagTimeBlock = False
307 252
308 253 self.utctime = None
309 254
310 255 self.timeZone = None
311 256
312 257 self.dstFlag = None
313 258
314 259 self.errorCount = None
315 260
316 261 self.nCohInt = None
317 262
318 263 self.blocksize = None
319 264
320 265 self.flagDecodeData = False #asumo q la data no esta decodificada
321 266
322 267 self.flagDeflipData = False #asumo q la data no esta sin flip
323 268
324 269 self.flagShiftFFT = False
325 270
326 271
327 272 def getNoisebyHildebrand(self):
328 273 """
329 274 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
330 275
331 276 Return:
332 277 noiselevel
333 278 """
334 279
335 280 for channel in range(self.nChannels):
336 281 daux = self.data_spc[channel,:,:]
337 282 self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt)
338 283
339 284 return self.noise
340 285
341 286 def getNoise(self, type = 1):
342 287
343 288 self.noise = numpy.zeros(self.nChannels)
344 289
345 290 if type == 1:
346 291 noise = self.getNoisebyHildebrand()
347 292
348 293 return 10*numpy.log10(noise)
349 294
350 295 class Spectra(JROData):
351 296
352 297 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
353 298 data_spc = None
354 299
355 300 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
356 301 data_cspc = None
357 302
358 303 #data es un numpy array de 2 dmensiones (canales, alturas)
359 304 data_dc = None
360 305
361 306 nFFTPoints = None
362 307
363 308 nPairs = None
364 309
365 310 pairsList = None
366 311
367 312 nIncohInt = None
368 313
369 314 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
370 315
371 316 nCohInt = None #se requiere para determinar el valor de timeInterval
372 317
373 318 ippFactor = None
374 319
375 320 def __init__(self):
376 321 '''
377 322 Constructor
378 323 '''
379 324
380 325 self.radarControllerHeaderObj = RadarControllerHeader()
381 326
382 327 self.systemHeaderObj = SystemHeader()
383 328
384 329 self.type = "Spectra"
385 330
386 331 # self.data = None
387 332
388 333 self.dtype = None
389 334
390 335 # self.nChannels = 0
391 336
392 337 # self.nHeights = 0
393 338
394 339 self.nProfiles = None
395 340
396 341 self.heightList = None
397 342
398 343 self.channelList = None
399 344
400 345 # self.channelIndexList = None
401 346
402 347 self.flagNoData = True
403 348
404 349 self.flagTimeBlock = False
405 350
406 351 self.utctime = None
407 352
408 353 self.nCohInt = None
409 354
410 355 self.nIncohInt = None
411 356
412 357 self.blocksize = None
413 358
414 359 self.nFFTPoints = None
415 360
416 361 self.wavelength = None
417 362
418 363 self.flagDecodeData = False #asumo q la data no esta decodificada
419 364
420 365 self.flagDeflipData = False #asumo q la data no esta sin flip
421 366
422 367 self.flagShiftFFT = False
423 368
424 369 self.ippFactor = 1
425 370
371 self.noise = None
372
426 373 def getNoisebyHildebrand(self):
427 374 """
428 375 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
429 376
430 377 Return:
431 378 noiselevel
432 379 """
433 380
434 381 for channel in range(self.nChannels):
435 382 daux = self.data_spc[channel,:,:]
436 383 self.noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
437 384
438 385 return self.noise
439 386
440 def getNoisebyWindow(self, heiIndexMin=0, heiIndexMax=-1, freqIndexMin=0, freqIndexMax=-1):
441 """
442 Determina el ruido del canal utilizando la ventana indicada con las coordenadas:
443 (heiIndexMIn, freqIndexMin) hasta (heiIndexMax, freqIndexMAx)
444
445 Inputs:
446 heiIndexMin: Limite inferior del eje de alturas
447 heiIndexMax: Limite superior del eje de alturas
448 freqIndexMin: Limite inferior del eje de frecuencia
449 freqIndexMax: Limite supoerior del eje de frecuencia
450 """
451
452 data = self.data_spc[:, heiIndexMin:heiIndexMax, freqIndexMin:freqIndexMax]
453
454 for channel in range(self.nChannels):
455 daux = data[channel,:,:]
456 self.noise[channel] = numpy.average(daux)
457
458 return self.noise
459
460 def getNoisebySort(self):
461
462 for channel in range(self.nChannels):
463 daux = self.data_spc[channel,:,:]
464 self.noise[channel] = sorting_bruce(daux, self.nIncohInt)
465
466 return self.noise
467
468 387 def getNoise(self, type = 1):
469 388 if self.noise == None:
470 389 self.noise = numpy.zeros(self.nChannels)
471
472 if type == 1:
473 390 self.noise = self.getNoisebyHildebrand()
474 391
475 if type == 2:
476 self.noise = self.getNoisebySort()
477
478 if type == 3:
479 self.noise = self.getNoisebyWindow()
480
481 392 return self.noise
482 393
483 394
484 395 def getFreqRange(self, extrapoints=0):
485 396
486 397 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
487 398 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
488 399
489 400 return freqrange
490 401
491 402 def getVelRange(self, extrapoints=0):
492 403
493 404 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
494 405 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
495 406
496 407 return velrange
497 408
498 409 def getNPairs(self):
499 410
500 411 return len(self.pairsList)
501 412
502 413 def getPairsIndexList(self):
503 414
504 415 return range(self.nPairs)
505 416
506 417 def getNormFactor(self):
507 418 pwcode = 1
508 419 if self.flagDecodeData:
509 420 pwcode = numpy.sum(self.code[0]**2)
510 421 normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode
511 422
512 423 return normFactor
513 424
514 425 def getFlagCspc(self):
515 426
516 427 if self.data_cspc == None:
517 428 return True
518 429
519 430 return False
520 431
521 432 def getFlagDc(self):
522 433
523 434 if self.data_dc == None:
524 435 return True
525 436
526 437 return False
527 438
528 439 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
529 440 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
530 441 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
531 442 flag_cspc = property(getFlagCspc)
532 443 flag_dc = property(getFlagDc)
533 444
534 445 class SpectraHeis(JROData):
535 446
536 447 data_spc = None
537 448
538 449 data_cspc = None
539 450
540 451 data_dc = None
541 452
542 453 nFFTPoints = None
543 454
544 455 nPairs = None
545 456
546 457 pairsList = None
547 458
548 459 nIncohInt = None
549 460
550 461 def __init__(self):
551 462
552 463 self.radarControllerHeaderObj = RadarControllerHeader()
553 464
554 465 self.systemHeaderObj = SystemHeader()
555 466
556 467 self.type = "SpectraHeis"
557 468
558 469 self.dtype = None
559 470
560 471 # self.nChannels = 0
561 472
562 473 # self.nHeights = 0
563 474
564 475 self.nProfiles = None
565 476
566 477 self.heightList = None
567 478
568 479 self.channelList = None
569 480
570 481 # self.channelIndexList = None
571 482
572 483 self.flagNoData = True
573 484
574 485 self.flagTimeBlock = False
575 486
576 487 self.nPairs = 0
577 488
578 489 self.utctime = None
579 490
580 491 self.blocksize = None
581 492
582 493 class Fits:
583 494
584 495 heightList = None
585 496
586 497 channelList = None
587 498
588 499 flagNoData = True
589 500
590 501 flagTimeBlock = False
591 502
592 503 useLocalTime = False
593 504
594 505 utctime = None
595 506
596 507 timeZone = None
597 508
598 509 ippSeconds = None
599 510
600 511 timeInterval = None
601 512
602 513 nCohInt = None
603 514
604 515 nIncohInt = None
605 516
606 517 noise = None
607 518
608 519 windowOfFilter = 1
609 520
610 521 #Speed of ligth
611 522 C = 3e8
612 523
613 524 frequency = 49.92e6
614 525
615 526 realtime = False
616 527
617 528
618 529 def __init__(self):
619 530
620 531 self.type = "Fits"
621 532
622 533 self.nProfiles = None
623 534
624 535 self.heightList = None
625 536
626 537 self.channelList = None
627 538
628 539 # self.channelIndexList = None
629 540
630 541 self.flagNoData = True
631 542
632 543 self.utctime = None
633 544
634 545 self.nCohInt = None
635 546
636 547 self.nIncohInt = None
637 548
638 549 self.useLocalTime = True
639 550
640 551 # self.utctime = None
641 552 # self.timeZone = None
642 553 # self.ltctime = None
643 554 # self.timeInterval = None
644 555 # self.header = None
645 556 # self.data_header = None
646 557 # self.data = None
647 558 # self.datatime = None
648 559 # self.flagNoData = False
649 560 # self.expName = ''
650 561 # self.nChannels = None
651 562 # self.nSamples = None
652 563 # self.dataBlocksPerFile = None
653 564 # self.comments = ''
654 565 #
655 566
656 567
657 568 def getltctime(self):
658 569
659 570 if self.useLocalTime:
660 571 return self.utctime - self.timeZone*60
661 572
662 573 return self.utctime
663 574
664 575 def getDatatime(self):
665 576
666 577 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
667 578 return datatime
668 579
669 580 def getTimeRange(self):
670 581
671 582 datatime = []
672 583
673 584 datatime.append(self.ltctime)
674 585 datatime.append(self.ltctime + self.timeInterval)
675 586
676 587 datatime = numpy.array(datatime)
677 588
678 589 return datatime
679 590
680 591 def getHeiRange(self):
681 592
682 593 heis = self.heightList
683 594
684 595 return heis
685 596
686 597 def isEmpty(self):
687 598
688 599 return self.flagNoData
689 600
690 601 def getNHeights(self):
691 602
692 603 return len(self.heightList)
693 604
694 605 def getNChannels(self):
695 606
696 607 return len(self.channelList)
697 608
698 609 def getChannelIndexList(self):
699 610
700 611 return range(self.nChannels)
701 612
702 613 def getNoise(self, type = 1):
703 614
704 615 self.noise = numpy.zeros(self.nChannels)
705 616
706 617 if type == 1:
707 618 noise = self.getNoisebyHildebrand()
708 619
709 620 if type == 2:
710 621 noise = self.getNoisebySort()
711 622
712 623 if type == 3:
713 624 noise = self.getNoisebyWindow()
714 625
715 626 return noise
716 627
717 628 datatime = property(getDatatime, "I'm the 'datatime' property")
718 629 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
719 630 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
720 631 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
721 632 noise = property(getNoise, "I'm the 'nHeights' property.")
722 633 datatime = property(getDatatime, "I'm the 'datatime' property")
723 634 ltctime = property(getltctime, "I'm the 'ltctime' property")
724 635
725 636 ltctime = property(getltctime, "I'm the 'ltctime' property") No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1960 +1,2004
1 1 '''
2 2
3 3 $Author: dsuarez $
4 4 $Id: Processor.py 1 2012-11-12 18:56:07Z dsuarez $
5 5 '''
6 6 import os
7 7 import numpy
8 8 import datetime
9 9 import time
10 10 import math
11 11 from jrodata import *
12 12 from jrodataIO import *
13 13 from jroplot import *
14 14
15 15 try:
16 16 import cfunctions
17 17 except:
18 18 pass
19 19
20 20 class ProcessingUnit:
21 21
22 22 """
23 23 Esta es la clase base para el procesamiento de datos.
24 24
25 25 Contiene el metodo "call" para llamar operaciones. Las operaciones pueden ser:
26 26 - Metodos internos (callMethod)
27 27 - Objetos del tipo Operation (callObject). Antes de ser llamados, estos objetos
28 28 tienen que ser agreagados con el metodo "add".
29 29
30 30 """
31 31 # objeto de datos de entrada (Voltage, Spectra o Correlation)
32 32 dataIn = None
33 33
34 34 # objeto de datos de entrada (Voltage, Spectra o Correlation)
35 35 dataOut = None
36 36
37 37
38 38 objectDict = None
39 39
40 40 def __init__(self):
41 41
42 42 self.objectDict = {}
43 43
44 44 def init(self):
45 45
46 46 raise ValueError, "Not implemented"
47 47
48 48 def addOperation(self, object, objId):
49 49
50 50 """
51 51 Agrega el objeto "object" a la lista de objetos "self.objectList" y retorna el
52 52 identificador asociado a este objeto.
53 53
54 54 Input:
55 55
56 56 object : objeto de la clase "Operation"
57 57
58 58 Return:
59 59
60 60 objId : identificador del objeto, necesario para ejecutar la operacion
61 61 """
62 62
63 63 self.objectDict[objId] = object
64 64
65 65 return objId
66 66
67 67 def operation(self, **kwargs):
68 68
69 69 """
70 70 Operacion directa sobre la data (dataOut.data). Es necesario actualizar los valores de los
71 71 atributos del objeto dataOut
72 72
73 73 Input:
74 74
75 75 **kwargs : Diccionario de argumentos de la funcion a ejecutar
76 76 """
77 77
78 78 raise ValueError, "ImplementedError"
79 79
80 80 def callMethod(self, name, **kwargs):
81 81
82 82 """
83 83 Ejecuta el metodo con el nombre "name" y con argumentos **kwargs de la propia clase.
84 84
85 85 Input:
86 86 name : nombre del metodo a ejecutar
87 87
88 88 **kwargs : diccionario con los nombres y valores de la funcion a ejecutar.
89 89
90 90 """
91 91 if name != 'run':
92 92
93 93 if name == 'init' and self.dataIn.isEmpty():
94 94 self.dataOut.flagNoData = True
95 95 return False
96 96
97 97 if name != 'init' and self.dataOut.isEmpty():
98 98 return False
99 99
100 100 methodToCall = getattr(self, name)
101 101
102 102 methodToCall(**kwargs)
103 103
104 104 if name != 'run':
105 105 return True
106 106
107 107 if self.dataOut.isEmpty():
108 108 return False
109 109
110 110 return True
111 111
112 112 def callObject(self, objId, **kwargs):
113 113
114 114 """
115 115 Ejecuta la operacion asociada al identificador del objeto "objId"
116 116
117 117 Input:
118 118
119 119 objId : identificador del objeto a ejecutar
120 120
121 121 **kwargs : diccionario con los nombres y valores de la funcion a ejecutar.
122 122
123 123 Return:
124 124
125 125 None
126 126 """
127 127
128 128 if self.dataOut.isEmpty():
129 129 return False
130 130
131 131 object = self.objectDict[objId]
132 132
133 133 object.run(self.dataOut, **kwargs)
134 134
135 135 return True
136 136
137 137 def call(self, operationConf, **kwargs):
138 138
139 139 """
140 140 Return True si ejecuta la operacion "operationConf.name" con los
141 141 argumentos "**kwargs". False si la operacion no se ha ejecutado.
142 142 La operacion puede ser de dos tipos:
143 143
144 144 1. Un metodo propio de esta clase:
145 145
146 146 operation.type = "self"
147 147
148 148 2. El metodo "run" de un objeto del tipo Operation o de un derivado de ella:
149 149 operation.type = "other".
150 150
151 151 Este objeto de tipo Operation debe de haber sido agregado antes con el metodo:
152 152 "addOperation" e identificado con el operation.id
153 153
154 154
155 155 con el id de la operacion.
156 156
157 157 Input:
158 158
159 159 Operation : Objeto del tipo operacion con los atributos: name, type y id.
160 160
161 161 """
162 162
163 163 if operationConf.type == 'self':
164 164 sts = self.callMethod(operationConf.name, **kwargs)
165 165
166 166 if operationConf.type == 'other':
167 167 sts = self.callObject(operationConf.id, **kwargs)
168 168
169 169 return sts
170 170
171 171 def setInput(self, dataIn):
172 172
173 173 self.dataIn = dataIn
174 174
175 175 def getOutput(self):
176 176
177 177 return self.dataOut
178 178
179 179 class Operation():
180 180
181 181 """
182 182 Clase base para definir las operaciones adicionales que se pueden agregar a la clase ProcessingUnit
183 183 y necesiten acumular informacion previa de los datos a procesar. De preferencia usar un buffer de
184 184 acumulacion dentro de esta clase
185 185
186 186 Ejemplo: Integraciones coherentes, necesita la informacion previa de los n perfiles anteriores (bufffer)
187 187
188 188 """
189 189
190 190 __buffer = None
191 191 __isConfig = False
192 192
193 193 def __init__(self):
194 194
195 195 pass
196 196
197 197 def run(self, dataIn, **kwargs):
198 198
199 199 """
200 200 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los atributos del objeto dataIn.
201 201
202 202 Input:
203 203
204 204 dataIn : objeto del tipo JROData
205 205
206 206 Return:
207 207
208 208 None
209 209
210 210 Affected:
211 211 __buffer : buffer de recepcion de datos.
212 212
213 213 """
214 214
215 215 raise ValueError, "ImplementedError"
216 216
217 217 class VoltageProc(ProcessingUnit):
218 218
219 219
220 220 def __init__(self):
221 221
222 222 self.objectDict = {}
223 223 self.dataOut = Voltage()
224 224 self.flip = 1
225 225
226 226 def init(self):
227 227
228 228 self.dataOut.copy(self.dataIn)
229 229 # No necesita copiar en cada init() los atributos de dataIn
230 230 # la copia deberia hacerse por cada nuevo bloque de datos
231 231
232 232 def selectChannels(self, channelList):
233 233
234 234 channelIndexList = []
235 235
236 236 for channel in channelList:
237 237 index = self.dataOut.channelList.index(channel)
238 238 channelIndexList.append(index)
239 239
240 240 self.selectChannelsByIndex(channelIndexList)
241 241
242 242 def selectChannelsByIndex(self, channelIndexList):
243 243 """
244 244 Selecciona un bloque de datos en base a canales segun el channelIndexList
245 245
246 246 Input:
247 247 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
248 248
249 249 Affected:
250 250 self.dataOut.data
251 251 self.dataOut.channelIndexList
252 252 self.dataOut.nChannels
253 253 self.dataOut.m_ProcessingHeader.totalSpectra
254 254 self.dataOut.systemHeaderObj.numChannels
255 255 self.dataOut.m_ProcessingHeader.blockSize
256 256
257 257 Return:
258 258 None
259 259 """
260 260
261 261 for channelIndex in channelIndexList:
262 262 if channelIndex not in self.dataOut.channelIndexList:
263 263 print channelIndexList
264 264 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
265 265
266 266 nChannels = len(channelIndexList)
267 267
268 268 data = self.dataOut.data[channelIndexList,:]
269 269
270 270 self.dataOut.data = data
271 271 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
272 272 # self.dataOut.nChannels = nChannels
273 273
274 274 return 1
275 275
276 276 def selectHeights(self, minHei=None, maxHei=None):
277 277 """
278 278 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
279 279 minHei <= height <= maxHei
280 280
281 281 Input:
282 282 minHei : valor minimo de altura a considerar
283 283 maxHei : valor maximo de altura a considerar
284 284
285 285 Affected:
286 286 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
287 287
288 288 Return:
289 289 1 si el metodo se ejecuto con exito caso contrario devuelve 0
290 290 """
291 291
292 292 if minHei == None:
293 293 minHei = self.dataOut.heightList[0]
294 294
295 295 if maxHei == None:
296 296 maxHei = self.dataOut.heightList[-1]
297 297
298 298 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
299 299 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
300 300
301 301
302 302 if (maxHei > self.dataOut.heightList[-1]):
303 303 maxHei = self.dataOut.heightList[-1]
304 304 # raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
305 305
306 306 minIndex = 0
307 307 maxIndex = 0
308 308 heights = self.dataOut.heightList
309 309
310 310 inda = numpy.where(heights >= minHei)
311 311 indb = numpy.where(heights <= maxHei)
312 312
313 313 try:
314 314 minIndex = inda[0][0]
315 315 except:
316 316 minIndex = 0
317 317
318 318 try:
319 319 maxIndex = indb[0][-1]
320 320 except:
321 321 maxIndex = len(heights)
322 322
323 323 self.selectHeightsByIndex(minIndex, maxIndex)
324 324
325 325 return 1
326 326
327 327
328 328 def selectHeightsByIndex(self, minIndex, maxIndex):
329 329 """
330 330 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
331 331 minIndex <= index <= maxIndex
332 332
333 333 Input:
334 334 minIndex : valor de indice minimo de altura a considerar
335 335 maxIndex : valor de indice maximo de altura a considerar
336 336
337 337 Affected:
338 338 self.dataOut.data
339 339 self.dataOut.heightList
340 340
341 341 Return:
342 342 1 si el metodo se ejecuto con exito caso contrario devuelve 0
343 343 """
344 344
345 345 if (minIndex < 0) or (minIndex > maxIndex):
346 346 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
347 347
348 348 if (maxIndex >= self.dataOut.nHeights):
349 349 maxIndex = self.dataOut.nHeights-1
350 350 # raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
351 351
352 352 nHeights = maxIndex - minIndex + 1
353 353
354 354 #voltage
355 355 data = self.dataOut.data[:,minIndex:maxIndex+1]
356 356
357 357 firstHeight = self.dataOut.heightList[minIndex]
358 358
359 359 self.dataOut.data = data
360 360 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
361 361
362 362 return 1
363 363
364 364
365 365 def filterByHeights(self, window):
366 366 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
367 367
368 368 if window == None:
369 369 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
370 370
371 371 newdelta = deltaHeight * window
372 372 r = self.dataOut.data.shape[1] % window
373 373 buffer = self.dataOut.data[:,0:self.dataOut.data.shape[1]-r]
374 374 buffer = buffer.reshape(self.dataOut.data.shape[0],self.dataOut.data.shape[1]/window,window)
375 375 buffer = numpy.sum(buffer,2)
376 376 self.dataOut.data = buffer
377 377 self.dataOut.heightList = numpy.arange(self.dataOut.heightList[0],newdelta*(self.dataOut.nHeights-r)/window,newdelta)
378 378 self.dataOut.windowOfFilter = window
379 379
380 380 def deFlip(self):
381 381 self.dataOut.data *= self.flip
382 382 self.flip *= -1.
383 383
384 384 def setRadarFrequency(self, frequency=None):
385 385 if frequency != None:
386 386 self.dataOut.frequency = frequency
387 387
388 388 return 1
389 389
390 390 class CohInt(Operation):
391 391
392 392 __isConfig = False
393 393
394 394 __profIndex = 0
395 395 __withOverapping = False
396 396
397 397 __byTime = False
398 398 __initime = None
399 399 __lastdatatime = None
400 400 __integrationtime = None
401 401
402 402 __buffer = None
403 403
404 404 __dataReady = False
405 405
406 406 n = None
407 407
408 408
409 409 def __init__(self):
410 410
411 411 self.__isConfig = False
412 412
413 413 def setup(self, n=None, timeInterval=None, overlapping=False):
414 414 """
415 415 Set the parameters of the integration class.
416 416
417 417 Inputs:
418 418
419 419 n : Number of coherent integrations
420 420 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
421 421 overlapping :
422 422
423 423 """
424 424
425 425 self.__initime = None
426 426 self.__lastdatatime = 0
427 427 self.__buffer = None
428 428 self.__dataReady = False
429 429
430 430
431 431 if n == None and timeInterval == None:
432 432 raise ValueError, "n or timeInterval should be specified ..."
433 433
434 434 if n != None:
435 435 self.n = n
436 436 self.__byTime = False
437 437 else:
438 438 self.__integrationtime = timeInterval * 60. #if (type(timeInterval)!=integer) -> change this line
439 439 self.n = 9999
440 440 self.__byTime = True
441 441
442 442 if overlapping:
443 443 self.__withOverapping = True
444 444 self.__buffer = None
445 445 else:
446 446 self.__withOverapping = False
447 447 self.__buffer = 0
448 448
449 449 self.__profIndex = 0
450 450
451 451 def putData(self, data):
452 452
453 453 """
454 454 Add a profile to the __buffer and increase in one the __profileIndex
455 455
456 456 """
457 457
458 458 if not self.__withOverapping:
459 459 self.__buffer += data.copy()
460 460 self.__profIndex += 1
461 461 return
462 462
463 463 #Overlapping data
464 464 nChannels, nHeis = data.shape
465 465 data = numpy.reshape(data, (1, nChannels, nHeis))
466 466
467 467 #If the buffer is empty then it takes the data value
468 468 if self.__buffer == None:
469 469 self.__buffer = data
470 470 self.__profIndex += 1
471 471 return
472 472
473 473 #If the buffer length is lower than n then stakcing the data value
474 474 if self.__profIndex < self.n:
475 475 self.__buffer = numpy.vstack((self.__buffer, data))
476 476 self.__profIndex += 1
477 477 return
478 478
479 479 #If the buffer length is equal to n then replacing the last buffer value with the data value
480 480 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
481 481 self.__buffer[self.n-1] = data
482 482 self.__profIndex = self.n
483 483 return
484 484
485 485
486 486 def pushData(self):
487 487 """
488 488 Return the sum of the last profiles and the profiles used in the sum.
489 489
490 490 Affected:
491 491
492 492 self.__profileIndex
493 493
494 494 """
495 495
496 496 if not self.__withOverapping:
497 497 data = self.__buffer
498 498 n = self.__profIndex
499 499
500 500 self.__buffer = 0
501 501 self.__profIndex = 0
502 502
503 503 return data, n
504 504
505 505 #Integration with Overlapping
506 506 data = numpy.sum(self.__buffer, axis=0)
507 507 n = self.__profIndex
508 508
509 509 return data, n
510 510
511 511 def byProfiles(self, data):
512 512
513 513 self.__dataReady = False
514 514 avgdata = None
515 515 n = None
516 516
517 517 self.putData(data)
518 518
519 519 if self.__profIndex == self.n:
520 520
521 521 avgdata, n = self.pushData()
522 522 self.__dataReady = True
523 523
524 524 return avgdata
525 525
526 526 def byTime(self, data, datatime):
527 527
528 528 self.__dataReady = False
529 529 avgdata = None
530 530 n = None
531 531
532 532 self.putData(data)
533 533
534 534 if (datatime - self.__initime) >= self.__integrationtime:
535 535 avgdata, n = self.pushData()
536 536 self.n = n
537 537 self.__dataReady = True
538 538
539 539 return avgdata
540 540
541 541 def integrate(self, data, datatime=None):
542 542
543 543 if self.__initime == None:
544 544 self.__initime = datatime
545 545
546 546 if self.__byTime:
547 547 avgdata = self.byTime(data, datatime)
548 548 else:
549 549 avgdata = self.byProfiles(data)
550 550
551 551
552 552 self.__lastdatatime = datatime
553 553
554 554 if avgdata == None:
555 555 return None, None
556 556
557 557 avgdatatime = self.__initime
558 558
559 559 deltatime = datatime -self.__lastdatatime
560 560
561 561 if not self.__withOverapping:
562 562 self.__initime = datatime
563 563 else:
564 564 self.__initime += deltatime
565 565
566 566 return avgdata, avgdatatime
567 567
568 568 def run(self, dataOut, **kwargs):
569 569
570 570 if not self.__isConfig:
571 571 self.setup(**kwargs)
572 572 self.__isConfig = True
573 573
574 574 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
575 575
576 576 # dataOut.timeInterval *= n
577 577 dataOut.flagNoData = True
578 578
579 579 if self.__dataReady:
580 580 dataOut.data = avgdata
581 581 dataOut.nCohInt *= self.n
582 582 dataOut.utctime = avgdatatime
583 583 dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
584 584 dataOut.flagNoData = False
585 585
586 586
587 587 class Decoder(Operation):
588 588
589 589 __isConfig = False
590 590 __profIndex = 0
591 591
592 592 code = None
593 593
594 594 nCode = None
595 595 nBaud = None
596 596
597 597 def __init__(self):
598 598
599 599 self.__isConfig = False
600 600
601 601 def setup(self, code, shape):
602 602
603 603 self.__profIndex = 0
604 604
605 605 self.code = code
606 606
607 607 self.nCode = len(code)
608 608 self.nBaud = len(code[0])
609 609
610 610 self.__nChannels, self.__nHeis = shape
611 611
612 612 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
613 613
614 614 __codeBuffer[:,0:self.nBaud] = self.code
615 615
616 616 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
617 617
618 618 self.ndatadec = self.__nHeis - self.nBaud + 1
619 619
620 620 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
621 621
622 622 def convolutionInFreq(self, data):
623 623
624 624 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
625 625
626 626 fft_data = numpy.fft.fft(data, axis=1)
627 627
628 628 conv = fft_data*fft_code
629 629
630 630 data = numpy.fft.ifft(conv,axis=1)
631 631
632 632 datadec = data[:,:-self.nBaud+1]
633 633
634 634 return datadec
635 635
636 636 def convolutionInFreqOpt(self, data):
637 637
638 638 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
639 639
640 640 data = cfunctions.decoder(fft_code, data)
641 641
642 642 datadec = data[:,:-self.nBaud+1]
643 643
644 644 return datadec
645 645
646 646 def convolutionInTime(self, data):
647 647
648 648 code = self.code[self.__profIndex]
649 649
650 650 for i in range(self.__nChannels):
651 651 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='valid')
652 652
653 653 return self.datadecTime
654 654
655 655 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0):
656 656
657 657 if not self.__isConfig:
658 658
659 659 if code == None:
660 660 code = dataOut.code
661 661 else:
662 662 code = numpy.array(code).reshape(nCode,nBaud)
663 663 dataOut.code = code
664 664 dataOut.nCode = nCode
665 665 dataOut.nBaud = nBaud
666 666
667 667 if code == None:
668 668 return 1
669 669
670 670 self.setup(code, dataOut.data.shape)
671 671 self.__isConfig = True
672 672
673 673 if mode == 0:
674 674 datadec = self.convolutionInTime(dataOut.data)
675 675
676 676 if mode == 1:
677 677 datadec = self.convolutionInFreq(dataOut.data)
678 678
679 679 if mode == 2:
680 680 datadec = self.convolutionInFreqOpt(dataOut.data)
681 681
682 682 dataOut.data = datadec
683 683
684 684 dataOut.heightList = dataOut.heightList[0:self.ndatadec]
685 685
686 686 dataOut.flagDecodeData = True #asumo q la data no esta decodificada
687 687
688 688 if self.__profIndex == self.nCode-1:
689 689 self.__profIndex = 0
690 690 return 1
691 691
692 692 self.__profIndex += 1
693 693
694 694 return 1
695 695 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
696 696
697 697
698 698
699 699 class SpectraProc(ProcessingUnit):
700 700
701 701 def __init__(self):
702 702
703 703 self.objectDict = {}
704 704 self.buffer = None
705 705 self.firstdatatime = None
706 706 self.profIndex = 0
707 707 self.dataOut = Spectra()
708 708
709 709 def __updateObjFromInput(self):
710 710
711 711 self.dataOut.timeZone = self.dataIn.timeZone
712 712 self.dataOut.dstFlag = self.dataIn.dstFlag
713 713 self.dataOut.errorCount = self.dataIn.errorCount
714 714 self.dataOut.useLocalTime = self.dataIn.useLocalTime
715 715
716 716 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
717 717 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
718 718 self.dataOut.channelList = self.dataIn.channelList
719 719 self.dataOut.heightList = self.dataIn.heightList
720 720 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
721 721 # self.dataOut.nHeights = self.dataIn.nHeights
722 722 # self.dataOut.nChannels = self.dataIn.nChannels
723 723 self.dataOut.nBaud = self.dataIn.nBaud
724 724 self.dataOut.nCode = self.dataIn.nCode
725 725 self.dataOut.code = self.dataIn.code
726 726 self.dataOut.nProfiles = self.dataOut.nFFTPoints
727 727 # self.dataOut.channelIndexList = self.dataIn.channelIndexList
728 728 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
729 729 self.dataOut.utctime = self.firstdatatime
730 730 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
731 731 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
732 732 # self.dataOut.flagShiftFFT = self.dataIn.flagShiftFFT
733 733 self.dataOut.nCohInt = self.dataIn.nCohInt
734 734 self.dataOut.nIncohInt = 1
735 735 self.dataOut.ippSeconds = self.dataIn.ippSeconds
736 736 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
737 737
738 738 self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nFFTPoints*self.dataOut.nIncohInt
739 739 self.dataOut.frequency = self.dataIn.frequency
740 740 self.dataOut.realtime = self.dataIn.realtime
741 741
742 742 def __getFft(self):
743 743 """
744 744 Convierte valores de Voltaje a Spectra
745 745
746 746 Affected:
747 747 self.dataOut.data_spc
748 748 self.dataOut.data_cspc
749 749 self.dataOut.data_dc
750 750 self.dataOut.heightList
751 751 self.profIndex
752 752 self.buffer
753 753 self.dataOut.flagNoData
754 754 """
755 755 fft_volt = numpy.fft.fft(self.buffer,n=self.dataOut.nFFTPoints,axis=1)
756 756 fft_volt = fft_volt.astype(numpy.dtype('complex'))
757 757 dc = fft_volt[:,0,:]
758 758
759 759 #calculo de self-spectra
760 760 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
761 761 spc = fft_volt * numpy.conjugate(fft_volt)
762 762 spc = spc.real
763 763
764 764 blocksize = 0
765 765 blocksize += dc.size
766 766 blocksize += spc.size
767 767
768 768 cspc = None
769 769 pairIndex = 0
770 770 if self.dataOut.pairsList != None:
771 771 #calculo de cross-spectra
772 772 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
773 773 for pair in self.dataOut.pairsList:
774 774 cspc[pairIndex,:,:] = fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:])
775 775 pairIndex += 1
776 776 blocksize += cspc.size
777 777
778 778 self.dataOut.data_spc = spc
779 779 self.dataOut.data_cspc = cspc
780 780 self.dataOut.data_dc = dc
781 781 self.dataOut.blockSize = blocksize
782 782 self.dataOut.flagShiftFFT = False
783 783
784 784 def init(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None):
785 785
786 786 self.dataOut.flagNoData = True
787 787
788 788 if self.dataIn.type == "Spectra":
789 789 self.dataOut.copy(self.dataIn)
790 790 return
791 791
792 792 if self.dataIn.type == "Voltage":
793 793
794 794 if nFFTPoints == None:
795 795 raise ValueError, "This SpectraProc.init() need nFFTPoints input variable"
796 796
797 797 if pairsList == None:
798 798 nPairs = 0
799 799 else:
800 800 nPairs = len(pairsList)
801 801
802 802 if ippFactor == None:
803 803 ippFactor = 1
804 804 self.dataOut.ippFactor = ippFactor
805 805
806 806 self.dataOut.nFFTPoints = nFFTPoints
807 807 self.dataOut.pairsList = pairsList
808 808 self.dataOut.nPairs = nPairs
809 809
810 810 if self.buffer == None:
811 811 self.buffer = numpy.zeros((self.dataIn.nChannels,
812 812 nProfiles,
813 813 self.dataIn.nHeights),
814 814 dtype='complex')
815 815
816 816
817 817 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
818 818 self.profIndex += 1
819 819
820 820 if self.firstdatatime == None:
821 821 self.firstdatatime = self.dataIn.utctime
822 822
823 823 if self.profIndex == nProfiles:
824 824 self.__updateObjFromInput()
825 825 self.__getFft()
826 826
827 827 self.dataOut.flagNoData = False
828 828
829 829 self.buffer = None
830 830 self.firstdatatime = None
831 831 self.profIndex = 0
832 832
833 833 return
834 834
835 835 raise ValueError, "The type object %s is not valid"%(self.dataIn.type)
836 836
837 837 def selectChannels(self, channelList):
838 838
839 839 channelIndexList = []
840 840
841 841 for channel in channelList:
842 842 index = self.dataOut.channelList.index(channel)
843 843 channelIndexList.append(index)
844 844
845 845 self.selectChannelsByIndex(channelIndexList)
846 846
847 847 def selectChannelsByIndex(self, channelIndexList):
848 848 """
849 849 Selecciona un bloque de datos en base a canales segun el channelIndexList
850 850
851 851 Input:
852 852 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
853 853
854 854 Affected:
855 855 self.dataOut.data_spc
856 856 self.dataOut.channelIndexList
857 857 self.dataOut.nChannels
858 858
859 859 Return:
860 860 None
861 861 """
862 862
863 863 for channelIndex in channelIndexList:
864 864 if channelIndex not in self.dataOut.channelIndexList:
865 865 print channelIndexList
866 866 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
867 867
868 868 nChannels = len(channelIndexList)
869 869
870 870 data_spc = self.dataOut.data_spc[channelIndexList,:]
871 871
872 872 self.dataOut.data_spc = data_spc
873 873 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
874 874 # self.dataOut.nChannels = nChannels
875 875
876 876 return 1
877 877
878 878 def selectHeights(self, minHei, maxHei):
879 879 """
880 880 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
881 881 minHei <= height <= maxHei
882 882
883 883 Input:
884 884 minHei : valor minimo de altura a considerar
885 885 maxHei : valor maximo de altura a considerar
886 886
887 887 Affected:
888 888 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
889 889
890 890 Return:
891 891 1 si el metodo se ejecuto con exito caso contrario devuelve 0
892 892 """
893 893 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
894 894 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
895 895
896 896 if (maxHei > self.dataOut.heightList[-1]):
897 897 maxHei = self.dataOut.heightList[-1]
898 898 # raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
899 899
900 900 minIndex = 0
901 901 maxIndex = 0
902 902 heights = self.dataOut.heightList
903 903
904 904 inda = numpy.where(heights >= minHei)
905 905 indb = numpy.where(heights <= maxHei)
906 906
907 907 try:
908 908 minIndex = inda[0][0]
909 909 except:
910 910 minIndex = 0
911 911
912 912 try:
913 913 maxIndex = indb[0][-1]
914 914 except:
915 915 maxIndex = len(heights)
916 916
917 917 self.selectHeightsByIndex(minIndex, maxIndex)
918 918
919 919 return 1
920 920
921 921
922 922 def selectHeightsByIndex(self, minIndex, maxIndex):
923 923 """
924 924 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
925 925 minIndex <= index <= maxIndex
926 926
927 927 Input:
928 928 minIndex : valor de indice minimo de altura a considerar
929 929 maxIndex : valor de indice maximo de altura a considerar
930 930
931 931 Affected:
932 932 self.dataOut.data_spc
933 933 self.dataOut.data_cspc
934 934 self.dataOut.data_dc
935 935 self.dataOut.heightList
936 936
937 937 Return:
938 938 1 si el metodo se ejecuto con exito caso contrario devuelve 0
939 939 """
940 940
941 941 if (minIndex < 0) or (minIndex > maxIndex):
942 942 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
943 943
944 944 if (maxIndex >= self.dataOut.nHeights):
945 945 maxIndex = self.dataOut.nHeights-1
946 946 # raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
947 947
948 948 nHeights = maxIndex - minIndex + 1
949 949
950 950 #Spectra
951 951 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
952 952
953 953 data_cspc = None
954 954 if self.dataOut.data_cspc != None:
955 955 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
956 956
957 957 data_dc = None
958 958 if self.dataOut.data_dc != None:
959 959 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
960 960
961 961 self.dataOut.data_spc = data_spc
962 962 self.dataOut.data_cspc = data_cspc
963 963 self.dataOut.data_dc = data_dc
964 964
965 965 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
966 966
967 967 return 1
968 968
969 969 def removeDC(self, mode = 2):
970 970 jspectra = self.dataOut.data_spc
971 971 jcspectra = self.dataOut.data_cspc
972 972
973 973
974 974 num_chan = jspectra.shape[0]
975 975 num_hei = jspectra.shape[2]
976 976
977 977 if jcspectra != None:
978 978 jcspectraExist = True
979 979 num_pairs = jcspectra.shape[0]
980 980 else: jcspectraExist = False
981 981
982 982 freq_dc = jspectra.shape[1]/2
983 983 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
984 984
985 985 if ind_vel[0]<0:
986 986 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
987 987
988 988 if mode == 1:
989 989 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
990 990
991 991 if jcspectraExist:
992 992 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
993 993
994 994 if mode == 2:
995 995
996 996 vel = numpy.array([-2,-1,1,2])
997 997 xx = numpy.zeros([4,4])
998 998
999 999 for fil in range(4):
1000 1000 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1001 1001
1002 1002 xx_inv = numpy.linalg.inv(xx)
1003 1003 xx_aux = xx_inv[0,:]
1004 1004
1005 1005 for ich in range(num_chan):
1006 1006 yy = jspectra[ich,ind_vel,:]
1007 1007 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1008 1008
1009 1009 junkid = jspectra[ich,freq_dc,:]<=0
1010 1010 cjunkid = sum(junkid)
1011 1011
1012 1012 if cjunkid.any():
1013 1013 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1014 1014
1015 1015 if jcspectraExist:
1016 1016 for ip in range(num_pairs):
1017 1017 yy = jcspectra[ip,ind_vel,:]
1018 1018 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
1019 1019
1020 1020
1021 1021 self.dataOut.data_spc = jspectra
1022 1022 self.dataOut.data_cspc = jcspectra
1023 1023
1024 1024 return 1
1025 1025
1026 1026 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
1027 1027
1028 1028 jspectra = self.dataOut.data_spc
1029 1029 jcspectra = self.dataOut.data_cspc
1030 1030 jnoise = self.dataOut.getNoise()
1031 1031 num_incoh = self.dataOut.nIncohInt
1032 1032
1033 1033 num_channel = jspectra.shape[0]
1034 1034 num_prof = jspectra.shape[1]
1035 1035 num_hei = jspectra.shape[2]
1036 1036
1037 1037 #hei_interf
1038 1038 if hei_interf == None:
1039 1039 count_hei = num_hei/2 #Como es entero no importa
1040 1040 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
1041 1041 hei_interf = numpy.asarray(hei_interf)[0]
1042 1042 #nhei_interf
1043 1043 if (nhei_interf == None):
1044 1044 nhei_interf = 5
1045 1045 if (nhei_interf < 1):
1046 1046 nhei_interf = 1
1047 1047 if (nhei_interf > count_hei):
1048 1048 nhei_interf = count_hei
1049 1049 if (offhei_interf == None):
1050 1050 offhei_interf = 0
1051 1051
1052 1052 ind_hei = range(num_hei)
1053 1053 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1054 1054 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1055 1055 mask_prof = numpy.asarray(range(num_prof))
1056 1056 num_mask_prof = mask_prof.size
1057 1057 comp_mask_prof = [0, num_prof/2]
1058 1058
1059 1059
1060 1060 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1061 1061 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1062 1062 jnoise = numpy.nan
1063 1063 noise_exist = jnoise[0] < numpy.Inf
1064 1064
1065 1065 #Subrutina de Remocion de la Interferencia
1066 1066 for ich in range(num_channel):
1067 1067 #Se ordena los espectros segun su potencia (menor a mayor)
1068 1068 power = jspectra[ich,mask_prof,:]
1069 1069 power = power[:,hei_interf]
1070 1070 power = power.sum(axis = 0)
1071 1071 psort = power.ravel().argsort()
1072 1072
1073 1073 #Se estima la interferencia promedio en los Espectros de Potencia empleando
1074 1074 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
1075 1075
1076 1076 if noise_exist:
1077 1077 # tmp_noise = jnoise[ich] / num_prof
1078 1078 tmp_noise = jnoise[ich]
1079 1079 junkspc_interf = junkspc_interf - tmp_noise
1080 1080 #junkspc_interf[:,comp_mask_prof] = 0
1081 1081
1082 1082 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
1083 1083 jspc_interf = jspc_interf.transpose()
1084 1084 #Calculando el espectro de interferencia promedio
1085 1085 noiseid = numpy.where(jspc_interf <= tmp_noise/ math.sqrt(num_incoh))
1086 1086 noiseid = noiseid[0]
1087 1087 cnoiseid = noiseid.size
1088 1088 interfid = numpy.where(jspc_interf > tmp_noise/ math.sqrt(num_incoh))
1089 1089 interfid = interfid[0]
1090 1090 cinterfid = interfid.size
1091 1091
1092 1092 if (cnoiseid > 0): jspc_interf[noiseid] = 0
1093 1093
1094 1094 #Expandiendo los perfiles a limpiar
1095 1095 if (cinterfid > 0):
1096 1096 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
1097 1097 new_interfid = numpy.asarray(new_interfid)
1098 1098 new_interfid = {x for x in new_interfid}
1099 1099 new_interfid = numpy.array(list(new_interfid))
1100 1100 new_cinterfid = new_interfid.size
1101 1101 else: new_cinterfid = 0
1102 1102
1103 1103 for ip in range(new_cinterfid):
1104 1104 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
1105 1105 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
1106 1106
1107 1107
1108 1108 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
1109 1109
1110 1110 #Removiendo la interferencia del punto de mayor interferencia
1111 1111 ListAux = jspc_interf[mask_prof].tolist()
1112 1112 maxid = ListAux.index(max(ListAux))
1113 1113
1114 1114
1115 1115 if cinterfid > 0:
1116 1116 for ip in range(cinterfid*(interf == 2) - 1):
1117 1117 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/math.sqrt(num_incoh))).nonzero()
1118 1118 cind = len(ind)
1119 1119
1120 1120 if (cind > 0):
1121 1121 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/math.sqrt(num_incoh))
1122 1122
1123 1123 ind = numpy.array([-2,-1,1,2])
1124 1124 xx = numpy.zeros([4,4])
1125 1125
1126 1126 for id1 in range(4):
1127 1127 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
1128 1128
1129 1129 xx_inv = numpy.linalg.inv(xx)
1130 1130 xx = xx_inv[:,0]
1131 1131 ind = (ind + maxid + num_mask_prof)%num_mask_prof
1132 1132 yy = jspectra[ich,mask_prof[ind],:]
1133 1133 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
1134 1134
1135 1135
1136 1136 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/math.sqrt(num_incoh))).nonzero()
1137 1137 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/math.sqrt(num_incoh))
1138 1138
1139 1139 #Remocion de Interferencia en el Cross Spectra
1140 1140 if jcspectra == None: return jspectra, jcspectra
1141 1141 num_pairs = jcspectra.size/(num_prof*num_hei)
1142 1142 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1143 1143
1144 1144 for ip in range(num_pairs):
1145 1145
1146 1146 #-------------------------------------------
1147 1147
1148 1148 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
1149 1149 cspower = cspower[:,hei_interf]
1150 1150 cspower = cspower.sum(axis = 0)
1151 1151
1152 1152 cspsort = cspower.ravel().argsort()
1153 1153 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
1154 1154 junkcspc_interf = junkcspc_interf.transpose()
1155 1155 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
1156 1156
1157 1157 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1158 1158
1159 1159 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
1160 1160 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
1161 1161 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
1162 1162
1163 1163 for iprof in range(num_prof):
1164 1164 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
1165 1165 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
1166 1166
1167 1167 #Removiendo la Interferencia
1168 1168 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
1169 1169
1170 1170 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1171 1171 maxid = ListAux.index(max(ListAux))
1172 1172
1173 1173 ind = numpy.array([-2,-1,1,2])
1174 1174 xx = numpy.zeros([4,4])
1175 1175
1176 1176 for id1 in range(4):
1177 1177 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
1178 1178
1179 1179 xx_inv = numpy.linalg.inv(xx)
1180 1180 xx = xx_inv[:,0]
1181 1181
1182 1182 ind = (ind + maxid + num_mask_prof)%num_mask_prof
1183 1183 yy = jcspectra[ip,mask_prof[ind],:]
1184 1184 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
1185 1185
1186 1186 #Guardar Resultados
1187 1187 self.dataOut.data_spc = jspectra
1188 1188 self.dataOut.data_cspc = jcspectra
1189 1189
1190 1190 return 1
1191 1191
1192 1192 def setRadarFrequency(self, frequency=None):
1193 1193 if frequency != None:
1194 1194 self.dataOut.frequency = frequency
1195 1195
1196 1196 return 1
1197 1197
1198 def getNoise(self, minHei, maxHei):
1198 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
1199 #validacion de rango
1200 if minHei == None:
1201 minHei = self.dataOut.heightList[0]
1202
1203 if maxHei == None:
1204 maxHei = self.dataOut.heightList[-1]
1199 1205
1200 1206 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
1201 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
1207 print 'minHei: %.2f is out of the heights range'%(minHei)
1208 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
1209 minHei = self.dataOut.heightList[0]
1202 1210
1203 if (maxHei > self.dataOut.heightList[-1]):
1211 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
1212 print 'maxHei: %.2f is out of the heights range'%(maxHei)
1213 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
1204 1214 maxHei = self.dataOut.heightList[-1]
1205 1215
1216 # validacion de velocidades
1217 velrange = self.dataOut.getVelRange(1)
1218
1219 if minVel == None:
1220 minVel = velrange[0]
1221
1222 if maxVel == None:
1223 maxVel = velrange[-1]
1224
1225 if (minVel < velrange[0]) or (minVel > maxVel):
1226 print 'minVel: %.2f is out of the velocity range'%(minVel)
1227 print 'minVel is setting to %.2f'%(velrange[0])
1228 minVel = velrange[0]
1229
1230 if (maxVel > velrange[-1]) or (maxVel < minVel):
1231 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
1232 print 'maxVel is setting to %.2f'%(velrange[-1])
1233 maxVel = velrange[-1]
1234
1235 # seleccion de indices para rango
1206 1236 minIndex = 0
1207 1237 maxIndex = 0
1208 1238 heights = self.dataOut.heightList
1209 1239
1210 1240 inda = numpy.where(heights >= minHei)
1211 1241 indb = numpy.where(heights <= maxHei)
1212 1242
1213 1243 try:
1214 1244 minIndex = inda[0][0]
1215 1245 except:
1216 1246 minIndex = 0
1217 1247
1218 1248 try:
1219 1249 maxIndex = indb[0][-1]
1220 1250 except:
1221 1251 maxIndex = len(heights)
1222 1252
1223 1253 if (minIndex < 0) or (minIndex > maxIndex):
1224 1254 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
1225 1255
1226 1256 if (maxIndex >= self.dataOut.nHeights):
1227 1257 maxIndex = self.dataOut.nHeights-1
1228 1258
1229 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
1259 # seleccion de indices para velocidades
1260 indminvel = numpy.where(velrange >= minVel)
1261 indmaxvel = numpy.where(velrange <= maxVel)
1262 try:
1263 minIndexVel = indminvel[0][0]
1264 except:
1265 minIndexVel = 0
1266
1267 try:
1268 maxIndexVel = indmaxvel[0][-1]
1269 except:
1270 maxIndexVel = len(velrange)
1230 1271
1272 #seleccion del espectro
1273 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
1274 #estimacion de ruido
1231 1275 noise = numpy.zeros(self.dataOut.nChannels)
1232 1276
1233 1277 for channel in range(self.dataOut.nChannels):
1234 1278 daux = data_spc[channel,:,:]
1235 1279 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
1236 1280
1237 1281 self.dataOut.noise = noise.copy()
1238 1282
1239 1283 return 1
1240 1284
1241 1285
1242 1286 class IncohInt(Operation):
1243 1287
1244 1288
1245 1289 __profIndex = 0
1246 1290 __withOverapping = False
1247 1291
1248 1292 __byTime = False
1249 1293 __initime = None
1250 1294 __lastdatatime = None
1251 1295 __integrationtime = None
1252 1296
1253 1297 __buffer_spc = None
1254 1298 __buffer_cspc = None
1255 1299 __buffer_dc = None
1256 1300
1257 1301 __dataReady = False
1258 1302
1259 1303 __timeInterval = None
1260 1304
1261 1305 n = None
1262 1306
1263 1307
1264 1308
1265 1309 def __init__(self):
1266 1310
1267 1311 self.__isConfig = False
1268 1312
1269 1313 def setup(self, n=None, timeInterval=None, overlapping=False):
1270 1314 """
1271 1315 Set the parameters of the integration class.
1272 1316
1273 1317 Inputs:
1274 1318
1275 1319 n : Number of coherent integrations
1276 1320 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1277 1321 overlapping :
1278 1322
1279 1323 """
1280 1324
1281 1325 self.__initime = None
1282 1326 self.__lastdatatime = 0
1283 1327 self.__buffer_spc = None
1284 1328 self.__buffer_cspc = None
1285 1329 self.__buffer_dc = None
1286 1330 self.__dataReady = False
1287 1331
1288 1332
1289 1333 if n == None and timeInterval == None:
1290 1334 raise ValueError, "n or timeInterval should be specified ..."
1291 1335
1292 1336 if n != None:
1293 1337 self.n = n
1294 1338 self.__byTime = False
1295 1339 else:
1296 1340 self.__integrationtime = timeInterval #if (type(timeInterval)!=integer) -> change this line
1297 1341 self.n = 9999
1298 1342 self.__byTime = True
1299 1343
1300 1344 if overlapping:
1301 1345 self.__withOverapping = True
1302 1346 else:
1303 1347 self.__withOverapping = False
1304 1348 self.__buffer_spc = 0
1305 1349 self.__buffer_cspc = 0
1306 1350 self.__buffer_dc = 0
1307 1351
1308 1352 self.__profIndex = 0
1309 1353
1310 1354 def putData(self, data_spc, data_cspc, data_dc):
1311 1355
1312 1356 """
1313 1357 Add a profile to the __buffer_spc and increase in one the __profileIndex
1314 1358
1315 1359 """
1316 1360
1317 1361 if not self.__withOverapping:
1318 1362 self.__buffer_spc += data_spc
1319 1363
1320 1364 if data_cspc == None:
1321 1365 self.__buffer_cspc = None
1322 1366 else:
1323 1367 self.__buffer_cspc += data_cspc
1324 1368
1325 1369 if data_dc == None:
1326 1370 self.__buffer_dc = None
1327 1371 else:
1328 1372 self.__buffer_dc += data_dc
1329 1373
1330 1374 self.__profIndex += 1
1331 1375 return
1332 1376
1333 1377 #Overlapping data
1334 1378 nChannels, nFFTPoints, nHeis = data_spc.shape
1335 1379 data_spc = numpy.reshape(data_spc, (1, nChannels, nFFTPoints, nHeis))
1336 1380 if data_cspc != None:
1337 1381 data_cspc = numpy.reshape(data_cspc, (1, -1, nFFTPoints, nHeis))
1338 1382 if data_dc != None:
1339 1383 data_dc = numpy.reshape(data_dc, (1, -1, nHeis))
1340 1384
1341 1385 #If the buffer is empty then it takes the data value
1342 1386 if self.__buffer_spc == None:
1343 1387 self.__buffer_spc = data_spc
1344 1388
1345 1389 if data_cspc == None:
1346 1390 self.__buffer_cspc = None
1347 1391 else:
1348 1392 self.__buffer_cspc += data_cspc
1349 1393
1350 1394 if data_dc == None:
1351 1395 self.__buffer_dc = None
1352 1396 else:
1353 1397 self.__buffer_dc += data_dc
1354 1398
1355 1399 self.__profIndex += 1
1356 1400 return
1357 1401
1358 1402 #If the buffer length is lower than n then stakcing the data value
1359 1403 if self.__profIndex < self.n:
1360 1404 self.__buffer_spc = numpy.vstack((self.__buffer_spc, data_spc))
1361 1405
1362 1406 if data_cspc != None:
1363 1407 self.__buffer_cspc = numpy.vstack((self.__buffer_cspc, data_cspc))
1364 1408
1365 1409 if data_dc != None:
1366 1410 self.__buffer_dc = numpy.vstack((self.__buffer_dc, data_dc))
1367 1411
1368 1412 self.__profIndex += 1
1369 1413 return
1370 1414
1371 1415 #If the buffer length is equal to n then replacing the last buffer value with the data value
1372 1416 self.__buffer_spc = numpy.roll(self.__buffer_spc, -1, axis=0)
1373 1417 self.__buffer_spc[self.n-1] = data_spc
1374 1418
1375 1419 if data_cspc != None:
1376 1420 self.__buffer_cspc = numpy.roll(self.__buffer_cspc, -1, axis=0)
1377 1421 self.__buffer_cspc[self.n-1] = data_cspc
1378 1422
1379 1423 if data_dc != None:
1380 1424 self.__buffer_dc = numpy.roll(self.__buffer_dc, -1, axis=0)
1381 1425 self.__buffer_dc[self.n-1] = data_dc
1382 1426
1383 1427 self.__profIndex = self.n
1384 1428 return
1385 1429
1386 1430
1387 1431 def pushData(self):
1388 1432 """
1389 1433 Return the sum of the last profiles and the profiles used in the sum.
1390 1434
1391 1435 Affected:
1392 1436
1393 1437 self.__profileIndex
1394 1438
1395 1439 """
1396 1440 data_spc = None
1397 1441 data_cspc = None
1398 1442 data_dc = None
1399 1443
1400 1444 if not self.__withOverapping:
1401 1445 data_spc = self.__buffer_spc
1402 1446 data_cspc = self.__buffer_cspc
1403 1447 data_dc = self.__buffer_dc
1404 1448
1405 1449 n = self.__profIndex
1406 1450
1407 1451 self.__buffer_spc = 0
1408 1452 self.__buffer_cspc = 0
1409 1453 self.__buffer_dc = 0
1410 1454 self.__profIndex = 0
1411 1455
1412 1456 return data_spc, data_cspc, data_dc, n
1413 1457
1414 1458 #Integration with Overlapping
1415 1459 data_spc = numpy.sum(self.__buffer_spc, axis=0)
1416 1460
1417 1461 if self.__buffer_cspc != None:
1418 1462 data_cspc = numpy.sum(self.__buffer_cspc, axis=0)
1419 1463
1420 1464 if self.__buffer_dc != None:
1421 1465 data_dc = numpy.sum(self.__buffer_dc, axis=0)
1422 1466
1423 1467 n = self.__profIndex
1424 1468
1425 1469 return data_spc, data_cspc, data_dc, n
1426 1470
1427 1471 def byProfiles(self, *args):
1428 1472
1429 1473 self.__dataReady = False
1430 1474 avgdata_spc = None
1431 1475 avgdata_cspc = None
1432 1476 avgdata_dc = None
1433 1477 n = None
1434 1478
1435 1479 self.putData(*args)
1436 1480
1437 1481 if self.__profIndex == self.n:
1438 1482
1439 1483 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1440 1484 self.__dataReady = True
1441 1485
1442 1486 return avgdata_spc, avgdata_cspc, avgdata_dc
1443 1487
1444 1488 def byTime(self, datatime, *args):
1445 1489
1446 1490 self.__dataReady = False
1447 1491 avgdata_spc = None
1448 1492 avgdata_cspc = None
1449 1493 avgdata_dc = None
1450 1494 n = None
1451 1495
1452 1496 self.putData(*args)
1453 1497
1454 1498 if (datatime - self.__initime) >= self.__integrationtime:
1455 1499 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1456 1500 self.n = n
1457 1501 self.__dataReady = True
1458 1502
1459 1503 return avgdata_spc, avgdata_cspc, avgdata_dc
1460 1504
1461 1505 def integrate(self, datatime, *args):
1462 1506
1463 1507 if self.__initime == None:
1464 1508 self.__initime = datatime
1465 1509
1466 1510 if self.__byTime:
1467 1511 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
1468 1512 else:
1469 1513 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1470 1514
1471 1515 self.__lastdatatime = datatime
1472 1516
1473 1517 if avgdata_spc == None:
1474 1518 return None, None, None, None
1475 1519
1476 1520 avgdatatime = self.__initime
1477 1521 try:
1478 1522 self.__timeInterval = (self.__lastdatatime - self.__initime)/(self.n - 1)
1479 1523 except:
1480 1524 self.__timeInterval = self.__lastdatatime - self.__initime
1481 1525
1482 1526 deltatime = datatime -self.__lastdatatime
1483 1527
1484 1528 if not self.__withOverapping:
1485 1529 self.__initime = datatime
1486 1530 else:
1487 1531 self.__initime += deltatime
1488 1532
1489 1533 return avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc
1490 1534
1491 1535 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1492 1536
1493 1537 if n==1:
1494 1538 dataOut.flagNoData = False
1495 1539 return
1496 1540
1497 1541 if not self.__isConfig:
1498 1542 self.setup(n, timeInterval, overlapping)
1499 1543 self.__isConfig = True
1500 1544
1501 1545 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1502 1546 dataOut.data_spc,
1503 1547 dataOut.data_cspc,
1504 1548 dataOut.data_dc)
1505 1549
1506 1550 # dataOut.timeInterval *= n
1507 1551 dataOut.flagNoData = True
1508 1552
1509 1553 if self.__dataReady:
1510 1554
1511 1555 dataOut.data_spc = avgdata_spc
1512 1556 dataOut.data_cspc = avgdata_cspc
1513 1557 dataOut.data_dc = avgdata_dc
1514 1558
1515 1559 dataOut.nIncohInt *= self.n
1516 1560 dataOut.utctime = avgdatatime
1517 1561 #dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt * dataOut.nIncohInt * dataOut.nFFTPoints
1518 1562 dataOut.timeInterval = self.__timeInterval*self.n
1519 1563 dataOut.flagNoData = False
1520 1564
1521 1565 class ProfileConcat(Operation):
1522 1566
1523 1567 __isConfig = False
1524 1568 buffer = None
1525 1569
1526 1570 def __init__(self):
1527 1571
1528 1572 self.profileIndex = 0
1529 1573
1530 1574 def reset(self):
1531 1575 self.buffer = numpy.zeros_like(self.buffer)
1532 1576 self.start_index = 0
1533 1577 self.times = 1
1534 1578
1535 1579 def setup(self, data, m, n=1):
1536 1580 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
1537 1581 self.profiles = data.shape[1]
1538 1582 self.start_index = 0
1539 1583 self.times = 1
1540 1584
1541 1585 def concat(self, data):
1542 1586
1543 1587 self.buffer[:,self.start_index:self.profiles*self.times] = data.copy()
1544 1588 self.start_index = self.start_index + self.profiles
1545 1589
1546 1590 def run(self, dataOut, m):
1547 1591
1548 1592 dataOut.flagNoData = True
1549 1593
1550 1594 if not self.__isConfig:
1551 1595 self.setup(dataOut.data, m, 1)
1552 1596 self.__isConfig = True
1553 1597
1554 1598 self.concat(dataOut.data)
1555 1599 self.times += 1
1556 1600 if self.times > m:
1557 1601 dataOut.data = self.buffer
1558 1602 self.reset()
1559 1603 dataOut.flagNoData = False
1560 1604 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
1561 1605 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1562 1606 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * 5
1563 1607 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
1564 1608
1565 1609
1566 1610
1567 1611 class ProfileSelector(Operation):
1568 1612
1569 1613 profileIndex = None
1570 1614 # Tamanho total de los perfiles
1571 1615 nProfiles = None
1572 1616
1573 1617 def __init__(self):
1574 1618
1575 1619 self.profileIndex = 0
1576 1620
1577 1621 def incIndex(self):
1578 1622 self.profileIndex += 1
1579 1623
1580 1624 if self.profileIndex >= self.nProfiles:
1581 1625 self.profileIndex = 0
1582 1626
1583 1627 def isProfileInRange(self, minIndex, maxIndex):
1584 1628
1585 1629 if self.profileIndex < minIndex:
1586 1630 return False
1587 1631
1588 1632 if self.profileIndex > maxIndex:
1589 1633 return False
1590 1634
1591 1635 return True
1592 1636
1593 1637 def isProfileInList(self, profileList):
1594 1638
1595 1639 if self.profileIndex not in profileList:
1596 1640 return False
1597 1641
1598 1642 return True
1599 1643
1600 1644 def run(self, dataOut, profileList=None, profileRangeList=None):
1601 1645
1602 1646 dataOut.flagNoData = True
1603 1647 self.nProfiles = dataOut.nProfiles
1604 1648
1605 1649 if profileList != None:
1606 1650 if self.isProfileInList(profileList):
1607 1651 dataOut.flagNoData = False
1608 1652
1609 1653 self.incIndex()
1610 1654 return 1
1611 1655
1612 1656
1613 1657 elif profileRangeList != None:
1614 1658 minIndex = profileRangeList[0]
1615 1659 maxIndex = profileRangeList[1]
1616 1660 if self.isProfileInRange(minIndex, maxIndex):
1617 1661 dataOut.flagNoData = False
1618 1662
1619 1663 self.incIndex()
1620 1664 return 1
1621 1665
1622 1666 else:
1623 1667 raise ValueError, "ProfileSelector needs profileList or profileRangeList"
1624 1668
1625 1669 return 0
1626 1670
1627 1671 class SpectraHeisProc(ProcessingUnit):
1628 1672 def __init__(self):
1629 1673 self.objectDict = {}
1630 1674 # self.buffer = None
1631 1675 # self.firstdatatime = None
1632 1676 # self.profIndex = 0
1633 1677 self.dataOut = SpectraHeis()
1634 1678
1635 1679 def __updateObjFromInput(self):
1636 1680 self.dataOut.timeZone = self.dataIn.timeZone
1637 1681 self.dataOut.dstFlag = self.dataIn.dstFlag
1638 1682 self.dataOut.errorCount = self.dataIn.errorCount
1639 1683 self.dataOut.useLocalTime = self.dataIn.useLocalTime
1640 1684
1641 1685 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()#
1642 1686 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()#
1643 1687 self.dataOut.channelList = self.dataIn.channelList
1644 1688 self.dataOut.heightList = self.dataIn.heightList
1645 1689 # self.dataOut.dtype = self.dataIn.dtype
1646 1690 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
1647 1691 # self.dataOut.nHeights = self.dataIn.nHeights
1648 1692 # self.dataOut.nChannels = self.dataIn.nChannels
1649 1693 self.dataOut.nBaud = self.dataIn.nBaud
1650 1694 self.dataOut.nCode = self.dataIn.nCode
1651 1695 self.dataOut.code = self.dataIn.code
1652 1696 # self.dataOut.nProfiles = 1
1653 1697 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
1654 1698 self.dataOut.nFFTPoints = self.dataIn.nHeights
1655 1699 # self.dataOut.channelIndexList = self.dataIn.channelIndexList
1656 1700 # self.dataOut.flagNoData = self.dataIn.flagNoData
1657 1701 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
1658 1702 self.dataOut.utctime = self.dataIn.utctime
1659 1703 # self.dataOut.utctime = self.firstdatatime
1660 1704 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
1661 1705 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
1662 1706 # self.dataOut.flagShiftFFT = self.dataIn.flagShiftFFT
1663 1707 self.dataOut.nCohInt = self.dataIn.nCohInt
1664 1708 self.dataOut.nIncohInt = 1
1665 1709 self.dataOut.ippSeconds= self.dataIn.ippSeconds
1666 1710 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
1667 1711
1668 1712 self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nIncohInt
1669 1713 # self.dataOut.set=self.dataIn.set
1670 1714 # self.dataOut.deltaHeight=self.dataIn.deltaHeight
1671 1715
1672 1716
1673 1717 def __updateObjFromFits(self):
1674 1718 self.dataOut.utctime = self.dataIn.utctime
1675 1719 self.dataOut.channelIndexList = self.dataIn.channelIndexList
1676 1720
1677 1721 self.dataOut.channelList = self.dataIn.channelList
1678 1722 self.dataOut.heightList = self.dataIn.heightList
1679 1723 self.dataOut.data_spc = self.dataIn.data
1680 1724 self.dataOut.timeInterval = self.dataIn.timeInterval
1681 1725 self.dataOut.timeZone = self.dataIn.timeZone
1682 1726 self.dataOut.useLocalTime = True
1683 1727 # self.dataOut.
1684 1728 # self.dataOut.
1685 1729
1686 1730 def __getFft(self):
1687 1731
1688 1732 fft_volt = numpy.fft.fft(self.dataIn.data, axis=1)
1689 1733 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
1690 1734 spc = numpy.abs(fft_volt * numpy.conjugate(fft_volt))/(self.dataOut.nFFTPoints)
1691 1735 self.dataOut.data_spc = spc
1692 1736
1693 1737 def init(self):
1694 1738
1695 1739 self.dataOut.flagNoData = True
1696 1740
1697 1741 if self.dataIn.type == "Fits":
1698 1742 self.__updateObjFromFits()
1699 1743 self.dataOut.flagNoData = False
1700 1744 return
1701 1745
1702 1746 if self.dataIn.type == "SpectraHeis":
1703 1747 self.dataOut.copy(self.dataIn)
1704 1748 return
1705 1749
1706 1750 if self.dataIn.type == "Voltage":
1707 1751 self.__updateObjFromInput()
1708 1752 self.__getFft()
1709 1753 self.dataOut.flagNoData = False
1710 1754
1711 1755 return
1712 1756
1713 1757 raise ValueError, "The type object %s is not valid"%(self.dataIn.type)
1714 1758
1715 1759
1716 1760 def selectChannels(self, channelList):
1717 1761
1718 1762 channelIndexList = []
1719 1763
1720 1764 for channel in channelList:
1721 1765 index = self.dataOut.channelList.index(channel)
1722 1766 channelIndexList.append(index)
1723 1767
1724 1768 self.selectChannelsByIndex(channelIndexList)
1725 1769
1726 1770 def selectChannelsByIndex(self, channelIndexList):
1727 1771 """
1728 1772 Selecciona un bloque de datos en base a canales segun el channelIndexList
1729 1773
1730 1774 Input:
1731 1775 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
1732 1776
1733 1777 Affected:
1734 1778 self.dataOut.data
1735 1779 self.dataOut.channelIndexList
1736 1780 self.dataOut.nChannels
1737 1781 self.dataOut.m_ProcessingHeader.totalSpectra
1738 1782 self.dataOut.systemHeaderObj.numChannels
1739 1783 self.dataOut.m_ProcessingHeader.blockSize
1740 1784
1741 1785 Return:
1742 1786 None
1743 1787 """
1744 1788
1745 1789 for channelIndex in channelIndexList:
1746 1790 if channelIndex not in self.dataOut.channelIndexList:
1747 1791 print channelIndexList
1748 1792 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
1749 1793
1750 1794 nChannels = len(channelIndexList)
1751 1795
1752 1796 data_spc = self.dataOut.data_spc[channelIndexList,:]
1753 1797
1754 1798 self.dataOut.data_spc = data_spc
1755 1799 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
1756 1800
1757 1801 return 1
1758 1802
1759 1803 class IncohInt4SpectraHeis(Operation):
1760 1804
1761 1805 __isConfig = False
1762 1806
1763 1807 __profIndex = 0
1764 1808 __withOverapping = False
1765 1809
1766 1810 __byTime = False
1767 1811 __initime = None
1768 1812 __lastdatatime = None
1769 1813 __integrationtime = None
1770 1814
1771 1815 __buffer = None
1772 1816
1773 1817 __dataReady = False
1774 1818
1775 1819 n = None
1776 1820
1777 1821
1778 1822 def __init__(self):
1779 1823
1780 1824 self.__isConfig = False
1781 1825
1782 1826 def setup(self, n=None, timeInterval=None, overlapping=False):
1783 1827 """
1784 1828 Set the parameters of the integration class.
1785 1829
1786 1830 Inputs:
1787 1831
1788 1832 n : Number of coherent integrations
1789 1833 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1790 1834 overlapping :
1791 1835
1792 1836 """
1793 1837
1794 1838 self.__initime = None
1795 1839 self.__lastdatatime = 0
1796 1840 self.__buffer = None
1797 1841 self.__dataReady = False
1798 1842
1799 1843
1800 1844 if n == None and timeInterval == None:
1801 1845 raise ValueError, "n or timeInterval should be specified ..."
1802 1846
1803 1847 if n != None:
1804 1848 self.n = n
1805 1849 self.__byTime = False
1806 1850 else:
1807 1851 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
1808 1852 self.n = 9999
1809 1853 self.__byTime = True
1810 1854
1811 1855 if overlapping:
1812 1856 self.__withOverapping = True
1813 1857 self.__buffer = None
1814 1858 else:
1815 1859 self.__withOverapping = False
1816 1860 self.__buffer = 0
1817 1861
1818 1862 self.__profIndex = 0
1819 1863
1820 1864 def putData(self, data):
1821 1865
1822 1866 """
1823 1867 Add a profile to the __buffer and increase in one the __profileIndex
1824 1868
1825 1869 """
1826 1870
1827 1871 if not self.__withOverapping:
1828 1872 self.__buffer += data.copy()
1829 1873 self.__profIndex += 1
1830 1874 return
1831 1875
1832 1876 #Overlapping data
1833 1877 nChannels, nHeis = data.shape
1834 1878 data = numpy.reshape(data, (1, nChannels, nHeis))
1835 1879
1836 1880 #If the buffer is empty then it takes the data value
1837 1881 if self.__buffer == None:
1838 1882 self.__buffer = data
1839 1883 self.__profIndex += 1
1840 1884 return
1841 1885
1842 1886 #If the buffer length is lower than n then stakcing the data value
1843 1887 if self.__profIndex < self.n:
1844 1888 self.__buffer = numpy.vstack((self.__buffer, data))
1845 1889 self.__profIndex += 1
1846 1890 return
1847 1891
1848 1892 #If the buffer length is equal to n then replacing the last buffer value with the data value
1849 1893 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
1850 1894 self.__buffer[self.n-1] = data
1851 1895 self.__profIndex = self.n
1852 1896 return
1853 1897
1854 1898
1855 1899 def pushData(self):
1856 1900 """
1857 1901 Return the sum of the last profiles and the profiles used in the sum.
1858 1902
1859 1903 Affected:
1860 1904
1861 1905 self.__profileIndex
1862 1906
1863 1907 """
1864 1908
1865 1909 if not self.__withOverapping:
1866 1910 data = self.__buffer
1867 1911 n = self.__profIndex
1868 1912
1869 1913 self.__buffer = 0
1870 1914 self.__profIndex = 0
1871 1915
1872 1916 return data, n
1873 1917
1874 1918 #Integration with Overlapping
1875 1919 data = numpy.sum(self.__buffer, axis=0)
1876 1920 n = self.__profIndex
1877 1921
1878 1922 return data, n
1879 1923
1880 1924 def byProfiles(self, data):
1881 1925
1882 1926 self.__dataReady = False
1883 1927 avgdata = None
1884 1928 n = None
1885 1929
1886 1930 self.putData(data)
1887 1931
1888 1932 if self.__profIndex == self.n:
1889 1933
1890 1934 avgdata, n = self.pushData()
1891 1935 self.__dataReady = True
1892 1936
1893 1937 return avgdata
1894 1938
1895 1939 def byTime(self, data, datatime):
1896 1940
1897 1941 self.__dataReady = False
1898 1942 avgdata = None
1899 1943 n = None
1900 1944
1901 1945 self.putData(data)
1902 1946
1903 1947 if (datatime - self.__initime) >= self.__integrationtime:
1904 1948 avgdata, n = self.pushData()
1905 1949 self.n = n
1906 1950 self.__dataReady = True
1907 1951
1908 1952 return avgdata
1909 1953
1910 1954 def integrate(self, data, datatime=None):
1911 1955
1912 1956 if self.__initime == None:
1913 1957 self.__initime = datatime
1914 1958
1915 1959 if self.__byTime:
1916 1960 avgdata = self.byTime(data, datatime)
1917 1961 else:
1918 1962 avgdata = self.byProfiles(data)
1919 1963
1920 1964
1921 1965 self.__lastdatatime = datatime
1922 1966
1923 1967 if avgdata == None:
1924 1968 return None, None
1925 1969
1926 1970 avgdatatime = self.__initime
1927 1971
1928 1972 deltatime = datatime -self.__lastdatatime
1929 1973
1930 1974 if not self.__withOverapping:
1931 1975 self.__initime = datatime
1932 1976 else:
1933 1977 self.__initime += deltatime
1934 1978
1935 1979 return avgdata, avgdatatime
1936 1980
1937 1981 def run(self, dataOut, **kwargs):
1938 1982
1939 1983 if not self.__isConfig:
1940 1984 self.setup(**kwargs)
1941 1985 self.__isConfig = True
1942 1986
1943 1987 avgdata, avgdatatime = self.integrate(dataOut.data_spc, dataOut.utctime)
1944 1988
1945 1989 # dataOut.timeInterval *= n
1946 1990 dataOut.flagNoData = True
1947 1991
1948 1992 if self.__dataReady:
1949 1993 dataOut.data_spc = avgdata
1950 1994 dataOut.nIncohInt *= self.n
1951 1995 # dataOut.nCohInt *= self.n
1952 1996 dataOut.utctime = avgdatatime
1953 1997 dataOut.timeInterval = dataOut.ippSeconds * dataOut.nIncohInt
1954 1998 # dataOut.timeInterval = self.__timeInterval*self.n
1955 1999 dataOut.flagNoData = False
1956 2000
1957 2001
1958 2002
1959 2003
1960 2004 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now