##// END OF EJS Templates
schain-jc
RSS
Se corrige error en la operacion removeDC modo=1 para espectros y cross-espectros....
Daniel Valdez -
r447:26a17bc7ac67
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,1917 +1,1917
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 jspectra[:freq_dc,:] = (jspectra[:ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
989 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
990 990
991 991 if jcspectraExist:
992 jcspectra[:freq_dc,:] = (jcspectra[:ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
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 1198
1199 1199 class IncohInt(Operation):
1200 1200
1201 1201
1202 1202 __profIndex = 0
1203 1203 __withOverapping = False
1204 1204
1205 1205 __byTime = False
1206 1206 __initime = None
1207 1207 __lastdatatime = None
1208 1208 __integrationtime = None
1209 1209
1210 1210 __buffer_spc = None
1211 1211 __buffer_cspc = None
1212 1212 __buffer_dc = None
1213 1213
1214 1214 __dataReady = False
1215 1215
1216 1216 __timeInterval = None
1217 1217
1218 1218 n = None
1219 1219
1220 1220
1221 1221
1222 1222 def __init__(self):
1223 1223
1224 1224 self.__isConfig = False
1225 1225
1226 1226 def setup(self, n=None, timeInterval=None, overlapping=False):
1227 1227 """
1228 1228 Set the parameters of the integration class.
1229 1229
1230 1230 Inputs:
1231 1231
1232 1232 n : Number of coherent integrations
1233 1233 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1234 1234 overlapping :
1235 1235
1236 1236 """
1237 1237
1238 1238 self.__initime = None
1239 1239 self.__lastdatatime = 0
1240 1240 self.__buffer_spc = None
1241 1241 self.__buffer_cspc = None
1242 1242 self.__buffer_dc = None
1243 1243 self.__dataReady = False
1244 1244
1245 1245
1246 1246 if n == None and timeInterval == None:
1247 1247 raise ValueError, "n or timeInterval should be specified ..."
1248 1248
1249 1249 if n != None:
1250 1250 self.n = n
1251 1251 self.__byTime = False
1252 1252 else:
1253 1253 self.__integrationtime = timeInterval #if (type(timeInterval)!=integer) -> change this line
1254 1254 self.n = 9999
1255 1255 self.__byTime = True
1256 1256
1257 1257 if overlapping:
1258 1258 self.__withOverapping = True
1259 1259 else:
1260 1260 self.__withOverapping = False
1261 1261 self.__buffer_spc = 0
1262 1262 self.__buffer_cspc = 0
1263 1263 self.__buffer_dc = 0
1264 1264
1265 1265 self.__profIndex = 0
1266 1266
1267 1267 def putData(self, data_spc, data_cspc, data_dc):
1268 1268
1269 1269 """
1270 1270 Add a profile to the __buffer_spc and increase in one the __profileIndex
1271 1271
1272 1272 """
1273 1273
1274 1274 if not self.__withOverapping:
1275 1275 self.__buffer_spc += data_spc
1276 1276
1277 1277 if data_cspc == None:
1278 1278 self.__buffer_cspc = None
1279 1279 else:
1280 1280 self.__buffer_cspc += data_cspc
1281 1281
1282 1282 if data_dc == None:
1283 1283 self.__buffer_dc = None
1284 1284 else:
1285 1285 self.__buffer_dc += data_dc
1286 1286
1287 1287 self.__profIndex += 1
1288 1288 return
1289 1289
1290 1290 #Overlapping data
1291 1291 nChannels, nFFTPoints, nHeis = data_spc.shape
1292 1292 data_spc = numpy.reshape(data_spc, (1, nChannels, nFFTPoints, nHeis))
1293 1293 if data_cspc != None:
1294 1294 data_cspc = numpy.reshape(data_cspc, (1, -1, nFFTPoints, nHeis))
1295 1295 if data_dc != None:
1296 1296 data_dc = numpy.reshape(data_dc, (1, -1, nHeis))
1297 1297
1298 1298 #If the buffer is empty then it takes the data value
1299 1299 if self.__buffer_spc == None:
1300 1300 self.__buffer_spc = data_spc
1301 1301
1302 1302 if data_cspc == None:
1303 1303 self.__buffer_cspc = None
1304 1304 else:
1305 1305 self.__buffer_cspc += data_cspc
1306 1306
1307 1307 if data_dc == None:
1308 1308 self.__buffer_dc = None
1309 1309 else:
1310 1310 self.__buffer_dc += data_dc
1311 1311
1312 1312 self.__profIndex += 1
1313 1313 return
1314 1314
1315 1315 #If the buffer length is lower than n then stakcing the data value
1316 1316 if self.__profIndex < self.n:
1317 1317 self.__buffer_spc = numpy.vstack((self.__buffer_spc, data_spc))
1318 1318
1319 1319 if data_cspc != None:
1320 1320 self.__buffer_cspc = numpy.vstack((self.__buffer_cspc, data_cspc))
1321 1321
1322 1322 if data_dc != None:
1323 1323 self.__buffer_dc = numpy.vstack((self.__buffer_dc, data_dc))
1324 1324
1325 1325 self.__profIndex += 1
1326 1326 return
1327 1327
1328 1328 #If the buffer length is equal to n then replacing the last buffer value with the data value
1329 1329 self.__buffer_spc = numpy.roll(self.__buffer_spc, -1, axis=0)
1330 1330 self.__buffer_spc[self.n-1] = data_spc
1331 1331
1332 1332 if data_cspc != None:
1333 1333 self.__buffer_cspc = numpy.roll(self.__buffer_cspc, -1, axis=0)
1334 1334 self.__buffer_cspc[self.n-1] = data_cspc
1335 1335
1336 1336 if data_dc != None:
1337 1337 self.__buffer_dc = numpy.roll(self.__buffer_dc, -1, axis=0)
1338 1338 self.__buffer_dc[self.n-1] = data_dc
1339 1339
1340 1340 self.__profIndex = self.n
1341 1341 return
1342 1342
1343 1343
1344 1344 def pushData(self):
1345 1345 """
1346 1346 Return the sum of the last profiles and the profiles used in the sum.
1347 1347
1348 1348 Affected:
1349 1349
1350 1350 self.__profileIndex
1351 1351
1352 1352 """
1353 1353 data_spc = None
1354 1354 data_cspc = None
1355 1355 data_dc = None
1356 1356
1357 1357 if not self.__withOverapping:
1358 1358 data_spc = self.__buffer_spc
1359 1359 data_cspc = self.__buffer_cspc
1360 1360 data_dc = self.__buffer_dc
1361 1361
1362 1362 n = self.__profIndex
1363 1363
1364 1364 self.__buffer_spc = 0
1365 1365 self.__buffer_cspc = 0
1366 1366 self.__buffer_dc = 0
1367 1367 self.__profIndex = 0
1368 1368
1369 1369 return data_spc, data_cspc, data_dc, n
1370 1370
1371 1371 #Integration with Overlapping
1372 1372 data_spc = numpy.sum(self.__buffer_spc, axis=0)
1373 1373
1374 1374 if self.__buffer_cspc != None:
1375 1375 data_cspc = numpy.sum(self.__buffer_cspc, axis=0)
1376 1376
1377 1377 if self.__buffer_dc != None:
1378 1378 data_dc = numpy.sum(self.__buffer_dc, axis=0)
1379 1379
1380 1380 n = self.__profIndex
1381 1381
1382 1382 return data_spc, data_cspc, data_dc, n
1383 1383
1384 1384 def byProfiles(self, *args):
1385 1385
1386 1386 self.__dataReady = False
1387 1387 avgdata_spc = None
1388 1388 avgdata_cspc = None
1389 1389 avgdata_dc = None
1390 1390 n = None
1391 1391
1392 1392 self.putData(*args)
1393 1393
1394 1394 if self.__profIndex == self.n:
1395 1395
1396 1396 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1397 1397 self.__dataReady = True
1398 1398
1399 1399 return avgdata_spc, avgdata_cspc, avgdata_dc
1400 1400
1401 1401 def byTime(self, datatime, *args):
1402 1402
1403 1403 self.__dataReady = False
1404 1404 avgdata_spc = None
1405 1405 avgdata_cspc = None
1406 1406 avgdata_dc = None
1407 1407 n = None
1408 1408
1409 1409 self.putData(*args)
1410 1410
1411 1411 if (datatime - self.__initime) >= self.__integrationtime:
1412 1412 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1413 1413 self.n = n
1414 1414 self.__dataReady = True
1415 1415
1416 1416 return avgdata_spc, avgdata_cspc, avgdata_dc
1417 1417
1418 1418 def integrate(self, datatime, *args):
1419 1419
1420 1420 if self.__initime == None:
1421 1421 self.__initime = datatime
1422 1422
1423 1423 if self.__byTime:
1424 1424 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
1425 1425 else:
1426 1426 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1427 1427
1428 1428 self.__lastdatatime = datatime
1429 1429
1430 1430 if avgdata_spc == None:
1431 1431 return None, None, None, None
1432 1432
1433 1433 avgdatatime = self.__initime
1434 1434 try:
1435 1435 self.__timeInterval = (self.__lastdatatime - self.__initime)/(self.n - 1)
1436 1436 except:
1437 1437 self.__timeInterval = self.__lastdatatime - self.__initime
1438 1438
1439 1439 deltatime = datatime -self.__lastdatatime
1440 1440
1441 1441 if not self.__withOverapping:
1442 1442 self.__initime = datatime
1443 1443 else:
1444 1444 self.__initime += deltatime
1445 1445
1446 1446 return avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc
1447 1447
1448 1448 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1449 1449
1450 1450 if n==1:
1451 1451 dataOut.flagNoData = False
1452 1452 return
1453 1453
1454 1454 if not self.__isConfig:
1455 1455 self.setup(n, timeInterval, overlapping)
1456 1456 self.__isConfig = True
1457 1457
1458 1458 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1459 1459 dataOut.data_spc,
1460 1460 dataOut.data_cspc,
1461 1461 dataOut.data_dc)
1462 1462
1463 1463 # dataOut.timeInterval *= n
1464 1464 dataOut.flagNoData = True
1465 1465
1466 1466 if self.__dataReady:
1467 1467
1468 1468 dataOut.data_spc = avgdata_spc
1469 1469 dataOut.data_cspc = avgdata_cspc
1470 1470 dataOut.data_dc = avgdata_dc
1471 1471
1472 1472 dataOut.nIncohInt *= self.n
1473 1473 dataOut.utctime = avgdatatime
1474 1474 #dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt * dataOut.nIncohInt * dataOut.nFFTPoints
1475 1475 dataOut.timeInterval = self.__timeInterval*self.n
1476 1476 dataOut.flagNoData = False
1477 1477
1478 1478 class ProfileConcat(Operation):
1479 1479
1480 1480 __isConfig = False
1481 1481 buffer = None
1482 1482
1483 1483 def __init__(self):
1484 1484
1485 1485 self.profileIndex = 0
1486 1486
1487 1487 def reset(self):
1488 1488 self.buffer = numpy.zeros_like(self.buffer)
1489 1489 self.start_index = 0
1490 1490 self.times = 1
1491 1491
1492 1492 def setup(self, data, m, n=1):
1493 1493 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
1494 1494 self.profiles = data.shape[1]
1495 1495 self.start_index = 0
1496 1496 self.times = 1
1497 1497
1498 1498 def concat(self, data):
1499 1499
1500 1500 self.buffer[:,self.start_index:self.profiles*self.times] = data.copy()
1501 1501 self.start_index = self.start_index + self.profiles
1502 1502
1503 1503 def run(self, dataOut, m):
1504 1504
1505 1505 dataOut.flagNoData = True
1506 1506
1507 1507 if not self.__isConfig:
1508 1508 self.setup(dataOut.data, m, 1)
1509 1509 self.__isConfig = True
1510 1510
1511 1511 self.concat(dataOut.data)
1512 1512 self.times += 1
1513 1513 if self.times > m:
1514 1514 dataOut.data = self.buffer
1515 1515 self.reset()
1516 1516 dataOut.flagNoData = False
1517 1517 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
1518 1518 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1519 1519 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * 5
1520 1520 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
1521 1521
1522 1522
1523 1523
1524 1524 class ProfileSelector(Operation):
1525 1525
1526 1526 profileIndex = None
1527 1527 # Tamanho total de los perfiles
1528 1528 nProfiles = None
1529 1529
1530 1530 def __init__(self):
1531 1531
1532 1532 self.profileIndex = 0
1533 1533
1534 1534 def incIndex(self):
1535 1535 self.profileIndex += 1
1536 1536
1537 1537 if self.profileIndex >= self.nProfiles:
1538 1538 self.profileIndex = 0
1539 1539
1540 1540 def isProfileInRange(self, minIndex, maxIndex):
1541 1541
1542 1542 if self.profileIndex < minIndex:
1543 1543 return False
1544 1544
1545 1545 if self.profileIndex > maxIndex:
1546 1546 return False
1547 1547
1548 1548 return True
1549 1549
1550 1550 def isProfileInList(self, profileList):
1551 1551
1552 1552 if self.profileIndex not in profileList:
1553 1553 return False
1554 1554
1555 1555 return True
1556 1556
1557 1557 def run(self, dataOut, profileList=None, profileRangeList=None):
1558 1558
1559 1559 dataOut.flagNoData = True
1560 1560 self.nProfiles = dataOut.nProfiles
1561 1561
1562 1562 if profileList != None:
1563 1563 if self.isProfileInList(profileList):
1564 1564 dataOut.flagNoData = False
1565 1565
1566 1566 self.incIndex()
1567 1567 return 1
1568 1568
1569 1569
1570 1570 elif profileRangeList != None:
1571 1571 minIndex = profileRangeList[0]
1572 1572 maxIndex = profileRangeList[1]
1573 1573 if self.isProfileInRange(minIndex, maxIndex):
1574 1574 dataOut.flagNoData = False
1575 1575
1576 1576 self.incIndex()
1577 1577 return 1
1578 1578
1579 1579 else:
1580 1580 raise ValueError, "ProfileSelector needs profileList or profileRangeList"
1581 1581
1582 1582 return 0
1583 1583
1584 1584 class SpectraHeisProc(ProcessingUnit):
1585 1585 def __init__(self):
1586 1586 self.objectDict = {}
1587 1587 # self.buffer = None
1588 1588 # self.firstdatatime = None
1589 1589 # self.profIndex = 0
1590 1590 self.dataOut = SpectraHeis()
1591 1591
1592 1592 def __updateObjFromInput(self):
1593 1593 self.dataOut.timeZone = self.dataIn.timeZone
1594 1594 self.dataOut.dstFlag = self.dataIn.dstFlag
1595 1595 self.dataOut.errorCount = self.dataIn.errorCount
1596 1596 self.dataOut.useLocalTime = self.dataIn.useLocalTime
1597 1597
1598 1598 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()#
1599 1599 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()#
1600 1600 self.dataOut.channelList = self.dataIn.channelList
1601 1601 self.dataOut.heightList = self.dataIn.heightList
1602 1602 # self.dataOut.dtype = self.dataIn.dtype
1603 1603 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
1604 1604 # self.dataOut.nHeights = self.dataIn.nHeights
1605 1605 # self.dataOut.nChannels = self.dataIn.nChannels
1606 1606 self.dataOut.nBaud = self.dataIn.nBaud
1607 1607 self.dataOut.nCode = self.dataIn.nCode
1608 1608 self.dataOut.code = self.dataIn.code
1609 1609 # self.dataOut.nProfiles = 1
1610 1610 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
1611 1611 self.dataOut.nFFTPoints = self.dataIn.nHeights
1612 1612 # self.dataOut.channelIndexList = self.dataIn.channelIndexList
1613 1613 # self.dataOut.flagNoData = self.dataIn.flagNoData
1614 1614 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
1615 1615 self.dataOut.utctime = self.dataIn.utctime
1616 1616 # self.dataOut.utctime = self.firstdatatime
1617 1617 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
1618 1618 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
1619 1619 # self.dataOut.flagShiftFFT = self.dataIn.flagShiftFFT
1620 1620 self.dataOut.nCohInt = self.dataIn.nCohInt
1621 1621 self.dataOut.nIncohInt = 1
1622 1622 self.dataOut.ippSeconds= self.dataIn.ippSeconds
1623 1623 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
1624 1624
1625 1625 self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nIncohInt
1626 1626 # self.dataOut.set=self.dataIn.set
1627 1627 # self.dataOut.deltaHeight=self.dataIn.deltaHeight
1628 1628
1629 1629
1630 1630 def __updateObjFromFits(self):
1631 1631 self.dataOut.utctime = self.dataIn.utctime
1632 1632 self.dataOut.channelIndexList = self.dataIn.channelIndexList
1633 1633
1634 1634 self.dataOut.channelList = self.dataIn.channelList
1635 1635 self.dataOut.heightList = self.dataIn.heightList
1636 1636 self.dataOut.data_spc = self.dataIn.data
1637 1637 self.dataOut.timeInterval = self.dataIn.timeInterval
1638 1638 self.dataOut.timeZone = self.dataIn.timeZone
1639 1639 self.dataOut.useLocalTime = True
1640 1640 # self.dataOut.
1641 1641 # self.dataOut.
1642 1642
1643 1643 def __getFft(self):
1644 1644
1645 1645 fft_volt = numpy.fft.fft(self.dataIn.data, axis=1)
1646 1646 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
1647 1647 spc = numpy.abs(fft_volt * numpy.conjugate(fft_volt))/(self.dataOut.nFFTPoints)
1648 1648 self.dataOut.data_spc = spc
1649 1649
1650 1650 def init(self):
1651 1651
1652 1652 self.dataOut.flagNoData = True
1653 1653
1654 1654 if self.dataIn.type == "Fits":
1655 1655 self.__updateObjFromFits()
1656 1656 self.dataOut.flagNoData = False
1657 1657 return
1658 1658
1659 1659 if self.dataIn.type == "SpectraHeis":
1660 1660 self.dataOut.copy(self.dataIn)
1661 1661 return
1662 1662
1663 1663 if self.dataIn.type == "Voltage":
1664 1664 self.__updateObjFromInput()
1665 1665 self.__getFft()
1666 1666 self.dataOut.flagNoData = False
1667 1667
1668 1668 return
1669 1669
1670 1670 raise ValueError, "The type object %s is not valid"%(self.dataIn.type)
1671 1671
1672 1672
1673 1673 def selectChannels(self, channelList):
1674 1674
1675 1675 channelIndexList = []
1676 1676
1677 1677 for channel in channelList:
1678 1678 index = self.dataOut.channelList.index(channel)
1679 1679 channelIndexList.append(index)
1680 1680
1681 1681 self.selectChannelsByIndex(channelIndexList)
1682 1682
1683 1683 def selectChannelsByIndex(self, channelIndexList):
1684 1684 """
1685 1685 Selecciona un bloque de datos en base a canales segun el channelIndexList
1686 1686
1687 1687 Input:
1688 1688 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
1689 1689
1690 1690 Affected:
1691 1691 self.dataOut.data
1692 1692 self.dataOut.channelIndexList
1693 1693 self.dataOut.nChannels
1694 1694 self.dataOut.m_ProcessingHeader.totalSpectra
1695 1695 self.dataOut.systemHeaderObj.numChannels
1696 1696 self.dataOut.m_ProcessingHeader.blockSize
1697 1697
1698 1698 Return:
1699 1699 None
1700 1700 """
1701 1701
1702 1702 for channelIndex in channelIndexList:
1703 1703 if channelIndex not in self.dataOut.channelIndexList:
1704 1704 print channelIndexList
1705 1705 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
1706 1706
1707 1707 nChannels = len(channelIndexList)
1708 1708
1709 1709 data_spc = self.dataOut.data_spc[channelIndexList,:]
1710 1710
1711 1711 self.dataOut.data_spc = data_spc
1712 1712 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
1713 1713
1714 1714 return 1
1715 1715
1716 1716 class IncohInt4SpectraHeis(Operation):
1717 1717
1718 1718 __isConfig = False
1719 1719
1720 1720 __profIndex = 0
1721 1721 __withOverapping = False
1722 1722
1723 1723 __byTime = False
1724 1724 __initime = None
1725 1725 __lastdatatime = None
1726 1726 __integrationtime = None
1727 1727
1728 1728 __buffer = None
1729 1729
1730 1730 __dataReady = False
1731 1731
1732 1732 n = None
1733 1733
1734 1734
1735 1735 def __init__(self):
1736 1736
1737 1737 self.__isConfig = False
1738 1738
1739 1739 def setup(self, n=None, timeInterval=None, overlapping=False):
1740 1740 """
1741 1741 Set the parameters of the integration class.
1742 1742
1743 1743 Inputs:
1744 1744
1745 1745 n : Number of coherent integrations
1746 1746 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1747 1747 overlapping :
1748 1748
1749 1749 """
1750 1750
1751 1751 self.__initime = None
1752 1752 self.__lastdatatime = 0
1753 1753 self.__buffer = None
1754 1754 self.__dataReady = False
1755 1755
1756 1756
1757 1757 if n == None and timeInterval == None:
1758 1758 raise ValueError, "n or timeInterval should be specified ..."
1759 1759
1760 1760 if n != None:
1761 1761 self.n = n
1762 1762 self.__byTime = False
1763 1763 else:
1764 1764 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
1765 1765 self.n = 9999
1766 1766 self.__byTime = True
1767 1767
1768 1768 if overlapping:
1769 1769 self.__withOverapping = True
1770 1770 self.__buffer = None
1771 1771 else:
1772 1772 self.__withOverapping = False
1773 1773 self.__buffer = 0
1774 1774
1775 1775 self.__profIndex = 0
1776 1776
1777 1777 def putData(self, data):
1778 1778
1779 1779 """
1780 1780 Add a profile to the __buffer and increase in one the __profileIndex
1781 1781
1782 1782 """
1783 1783
1784 1784 if not self.__withOverapping:
1785 1785 self.__buffer += data.copy()
1786 1786 self.__profIndex += 1
1787 1787 return
1788 1788
1789 1789 #Overlapping data
1790 1790 nChannels, nHeis = data.shape
1791 1791 data = numpy.reshape(data, (1, nChannels, nHeis))
1792 1792
1793 1793 #If the buffer is empty then it takes the data value
1794 1794 if self.__buffer == None:
1795 1795 self.__buffer = data
1796 1796 self.__profIndex += 1
1797 1797 return
1798 1798
1799 1799 #If the buffer length is lower than n then stakcing the data value
1800 1800 if self.__profIndex < self.n:
1801 1801 self.__buffer = numpy.vstack((self.__buffer, data))
1802 1802 self.__profIndex += 1
1803 1803 return
1804 1804
1805 1805 #If the buffer length is equal to n then replacing the last buffer value with the data value
1806 1806 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
1807 1807 self.__buffer[self.n-1] = data
1808 1808 self.__profIndex = self.n
1809 1809 return
1810 1810
1811 1811
1812 1812 def pushData(self):
1813 1813 """
1814 1814 Return the sum of the last profiles and the profiles used in the sum.
1815 1815
1816 1816 Affected:
1817 1817
1818 1818 self.__profileIndex
1819 1819
1820 1820 """
1821 1821
1822 1822 if not self.__withOverapping:
1823 1823 data = self.__buffer
1824 1824 n = self.__profIndex
1825 1825
1826 1826 self.__buffer = 0
1827 1827 self.__profIndex = 0
1828 1828
1829 1829 return data, n
1830 1830
1831 1831 #Integration with Overlapping
1832 1832 data = numpy.sum(self.__buffer, axis=0)
1833 1833 n = self.__profIndex
1834 1834
1835 1835 return data, n
1836 1836
1837 1837 def byProfiles(self, data):
1838 1838
1839 1839 self.__dataReady = False
1840 1840 avgdata = None
1841 1841 n = None
1842 1842
1843 1843 self.putData(data)
1844 1844
1845 1845 if self.__profIndex == self.n:
1846 1846
1847 1847 avgdata, n = self.pushData()
1848 1848 self.__dataReady = True
1849 1849
1850 1850 return avgdata
1851 1851
1852 1852 def byTime(self, data, datatime):
1853 1853
1854 1854 self.__dataReady = False
1855 1855 avgdata = None
1856 1856 n = None
1857 1857
1858 1858 self.putData(data)
1859 1859
1860 1860 if (datatime - self.__initime) >= self.__integrationtime:
1861 1861 avgdata, n = self.pushData()
1862 1862 self.n = n
1863 1863 self.__dataReady = True
1864 1864
1865 1865 return avgdata
1866 1866
1867 1867 def integrate(self, data, datatime=None):
1868 1868
1869 1869 if self.__initime == None:
1870 1870 self.__initime = datatime
1871 1871
1872 1872 if self.__byTime:
1873 1873 avgdata = self.byTime(data, datatime)
1874 1874 else:
1875 1875 avgdata = self.byProfiles(data)
1876 1876
1877 1877
1878 1878 self.__lastdatatime = datatime
1879 1879
1880 1880 if avgdata == None:
1881 1881 return None, None
1882 1882
1883 1883 avgdatatime = self.__initime
1884 1884
1885 1885 deltatime = datatime -self.__lastdatatime
1886 1886
1887 1887 if not self.__withOverapping:
1888 1888 self.__initime = datatime
1889 1889 else:
1890 1890 self.__initime += deltatime
1891 1891
1892 1892 return avgdata, avgdatatime
1893 1893
1894 1894 def run(self, dataOut, **kwargs):
1895 1895
1896 1896 if not self.__isConfig:
1897 1897 self.setup(**kwargs)
1898 1898 self.__isConfig = True
1899 1899
1900 1900 avgdata, avgdatatime = self.integrate(dataOut.data_spc, dataOut.utctime)
1901 1901
1902 1902 # dataOut.timeInterval *= n
1903 1903 dataOut.flagNoData = True
1904 1904
1905 1905 if self.__dataReady:
1906 1906 dataOut.data_spc = avgdata
1907 1907 dataOut.nIncohInt *= self.n
1908 1908 # dataOut.nCohInt *= self.n
1909 1909 dataOut.utctime = avgdatatime
1910 1910 dataOut.timeInterval = dataOut.ippSeconds * dataOut.nIncohInt
1911 1911 # dataOut.timeInterval = self.__timeInterval*self.n
1912 1912 dataOut.flagNoData = False
1913 1913
1914 1914
1915 1915
1916 1916
1917 1917 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now