##// END OF EJS Templates
schain-jc
RSS
Fix saving shifted pdata, add shift_fft parameter in SpectraProc to correct shifted pdata
jespinoza -
r1132:b323193ace39
parent child
Show More
Side by Side Unified
Context file:
r1132:b323193ace39 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,679 +1,679
1 1 '''
2 2 Created on Jul 2, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import numpy
7 7
8 8 from jroIO_base import LOCALTIME, JRODataReader, JRODataWriter
9 9 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
10 10 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader
11 11 from schainpy.model.data.jrodata import Spectra
12 12
13 13 class SpectraReader(JRODataReader, ProcessingUnit):
14 14 """
15 15 Esta clase permite leer datos de espectros desde archivos procesados (.pdata). La lectura
16 16 de los datos siempre se realiza por bloques. Los datos leidos (array de 3 dimensiones)
17 17 son almacenados en tres buffer's para el Self Spectra, el Cross Spectra y el DC Channel.
18 18
19 19 paresCanalesIguales * alturas * perfiles (Self Spectra)
20 20 paresCanalesDiferentes * alturas * perfiles (Cross Spectra)
21 21 canales * alturas (DC Channels)
22 22
23 23 Esta clase contiene instancias (objetos) de las clases BasicHeader, SystemHeader,
24 24 RadarControllerHeader y Spectra. Los tres primeros se usan para almacenar informacion de la
25 25 cabecera de datos (metadata), y el cuarto (Spectra) para obtener y almacenar un bloque de
26 26 datos desde el "buffer" cada vez que se ejecute el metodo "getData".
27 27
28 28 Example:
29 29 dpath = "/home/myuser/data"
30 30
31 31 startTime = datetime.datetime(2010,1,20,0,0,0,0,0,0)
32 32
33 33 endTime = datetime.datetime(2010,1,21,23,59,59,0,0,0)
34 34
35 35 readerObj = SpectraReader()
36 36
37 37 readerObj.setup(dpath, startTime, endTime)
38 38
39 39 while(True):
40 40
41 41 readerObj.getData()
42 42
43 43 print readerObj.data_spc
44 44
45 45 print readerObj.data_cspc
46 46
47 47 print readerObj.data_dc
48 48
49 49 if readerObj.flagNoMoreFiles:
50 50 break
51 51
52 52 """
53 53
54 54 pts2read_SelfSpectra = 0
55 55
56 56 pts2read_CrossSpectra = 0
57 57
58 58 pts2read_DCchannels = 0
59 59
60 60 ext = ".pdata"
61 61
62 62 optchar = "P"
63 63
64 64 dataOut = None
65 65
66 66 nRdChannels = None
67 67
68 68 nRdPairs = None
69 69
70 70 rdPairList = []
71 71
72 72 def __init__(self, **kwargs):
73 73 """
74 74 Inicializador de la clase SpectraReader para la lectura de datos de espectros.
75 75
76 76 Inputs:
77 77 dataOut : Objeto de la clase Spectra. Este objeto sera utilizado para
78 78 almacenar un perfil de datos cada vez que se haga un requerimiento
79 79 (getData). El perfil sera obtenido a partir del buffer de datos,
80 80 si el buffer esta vacio se hara un nuevo proceso de lectura de un
81 81 bloque de datos.
82 82 Si este parametro no es pasado se creara uno internamente.
83 83
84 84 Affected:
85 85 self.dataOut
86 86
87 87 Return : None
88 88 """
89 89
90 90 #Eliminar de la base la herencia
91 91 ProcessingUnit.__init__(self, **kwargs)
92 92
93 93 # self.isConfig = False
94 94
95 95 self.pts2read_SelfSpectra = 0
96 96
97 97 self.pts2read_CrossSpectra = 0
98 98
99 99 self.pts2read_DCchannels = 0
100 100
101 101 self.datablock = None
102 102
103 103 self.utc = None
104 104
105 105 self.ext = ".pdata"
106 106
107 107 self.optchar = "P"
108 108
109 109 self.basicHeaderObj = BasicHeader(LOCALTIME)
110 110
111 111 self.systemHeaderObj = SystemHeader()
112 112
113 113 self.radarControllerHeaderObj = RadarControllerHeader()
114 114
115 115 self.processingHeaderObj = ProcessingHeader()
116 116
117 117 self.online = 0
118 118
119 119 self.fp = None
120 120
121 121 self.idFile = None
122 122
123 123 self.dtype = None
124 124
125 125 self.fileSizeByHeader = None
126 126
127 127 self.filenameList = []
128 128
129 129 self.filename = None
130 130
131 131 self.fileSize = None
132 132
133 133 self.firstHeaderSize = 0
134 134
135 135 self.basicHeaderSize = 24
136 136
137 137 self.pathList = []
138 138
139 139 self.lastUTTime = 0
140 140
141 141 self.maxTimeStep = 30
142 142
143 143 self.flagNoMoreFiles = 0
144 144
145 145 self.set = 0
146 146
147 147 self.path = None
148 148
149 149 self.delay = 60 #seconds
150 150
151 151 self.nTries = 3 #quantity tries
152 152
153 153 self.nFiles = 3 #number of files for searching
154 154
155 155 self.nReadBlocks = 0
156 156
157 157 self.flagIsNewFile = 1
158 158
159 159 self.__isFirstTimeOnline = 1
160 160
161 161 # self.ippSeconds = 0
162 162
163 163 self.flagDiscontinuousBlock = 0
164 164
165 165 self.flagIsNewBlock = 0
166 166
167 167 self.nTotalBlocks = 0
168 168
169 169 self.blocksize = 0
170 170
171 171 self.dataOut = self.createObjByDefault()
172 172
173 173 self.profileIndex = 1 #Always
174 174
175 175
176 176 def createObjByDefault(self):
177 177
178 178 dataObj = Spectra()
179 179
180 180 return dataObj
181 181
182 182 def __hasNotDataInBuffer(self):
183 183 return 1
184 184
185 185
186 186 def getBlockDimension(self):
187 187 """
188 188 Obtiene la cantidad de puntos a leer por cada bloque de datos
189 189
190 190 Affected:
191 191 self.nRdChannels
192 192 self.nRdPairs
193 193 self.pts2read_SelfSpectra
194 194 self.pts2read_CrossSpectra
195 195 self.pts2read_DCchannels
196 196 self.blocksize
197 197 self.dataOut.nChannels
198 198 self.dataOut.nPairs
199 199
200 200 Return:
201 201 None
202 202 """
203 203 self.nRdChannels = 0
204 204 self.nRdPairs = 0
205 205 self.rdPairList = []
206 206
207 207 for i in range(0, self.processingHeaderObj.totalSpectra*2, 2):
208 208 if self.processingHeaderObj.spectraComb[i] == self.processingHeaderObj.spectraComb[i+1]:
209 209 self.nRdChannels = self.nRdChannels + 1 #par de canales iguales
210 210 else:
211 211 self.nRdPairs = self.nRdPairs + 1 #par de canales diferentes
212 212 self.rdPairList.append((self.processingHeaderObj.spectraComb[i], self.processingHeaderObj.spectraComb[i+1]))
213 213
214 214 pts2read = self.processingHeaderObj.nHeights * self.processingHeaderObj.profilesPerBlock
215 215
216 216 self.pts2read_SelfSpectra = int(self.nRdChannels * pts2read)
217 217 self.blocksize = self.pts2read_SelfSpectra
218 218
219 219 if self.processingHeaderObj.flag_cspc:
220 220 self.pts2read_CrossSpectra = int(self.nRdPairs * pts2read)
221 221 self.blocksize += self.pts2read_CrossSpectra
222 222
223 223 if self.processingHeaderObj.flag_dc:
224 224 self.pts2read_DCchannels = int(self.systemHeaderObj.nChannels * self.processingHeaderObj.nHeights)
225 225 self.blocksize += self.pts2read_DCchannels
226 226
227 227 # self.blocksize = self.pts2read_SelfSpectra + self.pts2read_CrossSpectra + self.pts2read_DCchannels
228 228
229 229
230 230 def readBlock(self):
231 231 """
232 232 Lee el bloque de datos desde la posicion actual del puntero del archivo
233 233 (self.fp) y actualiza todos los parametros relacionados al bloque de datos
234 234 (metadata + data). La data leida es almacenada en el buffer y el contador del buffer
235 235 es seteado a 0
236 236
237 237 Return: None
238 238
239 239 Variables afectadas:
240 240
241 241 self.flagIsNewFile
242 242 self.flagIsNewBlock
243 243 self.nTotalBlocks
244 244 self.data_spc
245 245 self.data_cspc
246 246 self.data_dc
247 247
248 248 Exceptions:
249 249 Si un bloque leido no es un bloque valido
250 250 """
251 251 blockOk_flag = False
252 252 fpointer = self.fp.tell()
253 253
254 254 spc = numpy.fromfile( self.fp, self.dtype[0], self.pts2read_SelfSpectra )
255 255 spc = spc.reshape( (self.nRdChannels, self.processingHeaderObj.nHeights, self.processingHeaderObj.profilesPerBlock) ) #transforma a un arreglo 3D
256 256
257 257 if self.processingHeaderObj.flag_cspc:
258 258 cspc = numpy.fromfile( self.fp, self.dtype, self.pts2read_CrossSpectra )
259 259 cspc = cspc.reshape( (self.nRdPairs, self.processingHeaderObj.nHeights, self.processingHeaderObj.profilesPerBlock) ) #transforma a un arreglo 3D
260 260
261 261 if self.processingHeaderObj.flag_dc:
262 262 dc = numpy.fromfile( self.fp, self.dtype, self.pts2read_DCchannels ) #int(self.processingHeaderObj.nHeights*self.systemHeaderObj.nChannels) )
263 263 dc = dc.reshape( (self.systemHeaderObj.nChannels, self.processingHeaderObj.nHeights) ) #transforma a un arreglo 2D
264 264
265 265
266 if not(self.processingHeaderObj.shif_fft):
266 if self.processingHeaderObj.shif_fft:
267 267 #desplaza a la derecha en el eje 2 determinadas posiciones
268 268 shift = int(self.processingHeaderObj.profilesPerBlock/2)
269 269 spc = numpy.roll( spc, shift , axis=2 )
270 270
271 271 if self.processingHeaderObj.flag_cspc:
272 272 #desplaza a la derecha en el eje 2 determinadas posiciones
273 273 cspc = numpy.roll( cspc, shift, axis=2 )
274 274
275 275 #Dimensions : nChannels, nProfiles, nSamples
276 276 spc = numpy.transpose( spc, (0,2,1) )
277 277 self.data_spc = spc
278 278
279 279 if self.processingHeaderObj.flag_cspc:
280 280 cspc = numpy.transpose( cspc, (0,2,1) )
281 281 self.data_cspc = cspc['real'] + cspc['imag']*1j
282 282 else:
283 283 self.data_cspc = None
284 284
285 285 if self.processingHeaderObj.flag_dc:
286 286 self.data_dc = dc['real'] + dc['imag']*1j
287 287 else:
288 288 self.data_dc = None
289 289
290 290 self.flagIsNewFile = 0
291 291 self.flagIsNewBlock = 1
292 292
293 293 self.nTotalBlocks += 1
294 294 self.nReadBlocks += 1
295 295
296 296 return 1
297 297
298 298 def getFirstHeader(self):
299 299
300 300 self.getBasicHeader()
301 301
302 302 self.dataOut.systemHeaderObj = self.systemHeaderObj.copy()
303 303
304 304 self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy()
305 305
306 306 # self.dataOut.ippSeconds = self.ippSeconds
307 307
308 308 # self.dataOut.timeInterval = self.radarControllerHeaderObj.ippSeconds * self.processingHeaderObj.nCohInt * self.processingHeaderObj.nIncohInt * self.processingHeaderObj.profilesPerBlock
309 309
310 310 self.dataOut.dtype = self.dtype
311 311
312 312 # self.dataOut.nPairs = self.nPairs
313 313
314 314 self.dataOut.pairsList = self.rdPairList
315 315
316 316 self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock
317 317
318 318 self.dataOut.nFFTPoints = self.processingHeaderObj.profilesPerBlock
319 319
320 320 self.dataOut.nCohInt = self.processingHeaderObj.nCohInt
321 321
322 322 self.dataOut.nIncohInt = self.processingHeaderObj.nIncohInt
323 323
324 324 xf = self.processingHeaderObj.firstHeight + self.processingHeaderObj.nHeights*self.processingHeaderObj.deltaHeight
325 325
326 326 self.dataOut.heightList = numpy.arange(self.processingHeaderObj.firstHeight, xf, self.processingHeaderObj.deltaHeight)
327 327
328 328 self.dataOut.channelList = range(self.systemHeaderObj.nChannels)
329 329
330 330 self.dataOut.flagShiftFFT = True #Data is always shifted
331 331
332 332 self.dataOut.flagDecodeData = self.processingHeaderObj.flag_decode #asumo q la data no esta decodificada
333 333
334 334 self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip #asumo q la data esta sin flip
335 335
336 336 def getData(self):
337 337 """
338 338 First method to execute before "RUN" is called.
339 339
340 340 Copia el buffer de lectura a la clase "Spectra",
341 341 con todos los parametros asociados a este (metadata). cuando no hay datos en el buffer de
342 342 lectura es necesario hacer una nueva lectura de los bloques de datos usando "readNextBlock"
343 343
344 344 Return:
345 345 0 : Si no hay mas archivos disponibles
346 346 1 : Si hizo una buena copia del buffer
347 347
348 348 Affected:
349 349 self.dataOut
350 350
351 351 self.flagDiscontinuousBlock
352 352 self.flagIsNewBlock
353 353 """
354 354
355 355 if self.flagNoMoreFiles:
356 356 self.dataOut.flagNoData = True
357 357 print 'Process finished'
358 358 return 0
359 359
360 360 self.flagDiscontinuousBlock = 0
361 361 self.flagIsNewBlock = 0
362 362
363 363 if self.__hasNotDataInBuffer():
364 364
365 365 if not( self.readNextBlock() ):
366 366 self.dataOut.flagNoData = True
367 367 return 0
368 368
369 369 #data es un numpy array de 3 dmensiones (perfiles, alturas y canales)
370 370
371 371 if self.data_spc is None:
372 372 self.dataOut.flagNoData = True
373 373 return 0
374 374
375 375 self.getBasicHeader()
376 376
377 377 self.getFirstHeader()
378 378
379 379 self.dataOut.data_spc = self.data_spc
380 380
381 381 self.dataOut.data_cspc = self.data_cspc
382 382
383 383 self.dataOut.data_dc = self.data_dc
384 384
385 385 self.dataOut.flagNoData = False
386 386
387 387 self.dataOut.realtime = self.online
388 388
389 389 return self.dataOut.data_spc
390 390
391 391 class SpectraWriter(JRODataWriter, Operation):
392 392
393 393 """
394 394 Esta clase permite escribir datos de espectros a archivos procesados (.pdata). La escritura
395 395 de los datos siempre se realiza por bloques.
396 396 """
397 397
398 398 ext = ".pdata"
399 399
400 400 optchar = "P"
401 401
402 402 shape_spc_Buffer = None
403 403
404 404 shape_cspc_Buffer = None
405 405
406 406 shape_dc_Buffer = None
407 407
408 408 data_spc = None
409 409
410 410 data_cspc = None
411 411
412 412 data_dc = None
413 413
414 414 # dataOut = None
415 415
416 416 def __init__(self, **kwargs):
417 417 """
418 418 Inicializador de la clase SpectraWriter para la escritura de datos de espectros.
419 419
420 420 Affected:
421 421 self.dataOut
422 422 self.basicHeaderObj
423 423 self.systemHeaderObj
424 424 self.radarControllerHeaderObj
425 425 self.processingHeaderObj
426 426
427 427 Return: None
428 428 """
429 429
430 430 Operation.__init__(self, **kwargs)
431 431
432 432 self.isConfig = False
433 433
434 434 self.nTotalBlocks = 0
435 435
436 436 self.data_spc = None
437 437
438 438 self.data_cspc = None
439 439
440 440 self.data_dc = None
441 441
442 442 self.fp = None
443 443
444 444 self.flagIsNewFile = 1
445 445
446 446 self.nTotalBlocks = 0
447 447
448 448 self.flagIsNewBlock = 0
449 449
450 450 self.setFile = None
451 451
452 452 self.dtype = None
453 453
454 454 self.path = None
455 455
456 456 self.noMoreFiles = 0
457 457
458 458 self.filename = None
459 459
460 460 self.basicHeaderObj = BasicHeader(LOCALTIME)
461 461
462 462 self.systemHeaderObj = SystemHeader()
463 463
464 464 self.radarControllerHeaderObj = RadarControllerHeader()
465 465
466 466 self.processingHeaderObj = ProcessingHeader()
467 467
468 468
469 469 def hasAllDataInBuffer(self):
470 470 return 1
471 471
472 472
473 473 def setBlockDimension(self):
474 474 """
475 475 Obtiene las formas dimensionales del los subbloques de datos que componen un bloque
476 476
477 477 Affected:
478 478 self.shape_spc_Buffer
479 479 self.shape_cspc_Buffer
480 480 self.shape_dc_Buffer
481 481
482 482 Return: None
483 483 """
484 484 self.shape_spc_Buffer = (self.dataOut.nChannels,
485 485 self.processingHeaderObj.nHeights,
486 486 self.processingHeaderObj.profilesPerBlock)
487 487
488 488 self.shape_cspc_Buffer = (self.dataOut.nPairs,
489 489 self.processingHeaderObj.nHeights,
490 490 self.processingHeaderObj.profilesPerBlock)
491 491
492 492 self.shape_dc_Buffer = (self.dataOut.nChannels,
493 493 self.processingHeaderObj.nHeights)
494 494
495 495
496 496 def writeBlock(self):
497 497 """
498 498 Escribe el buffer en el file designado
499 499
500 500 Affected:
501 501 self.data_spc
502 502 self.data_cspc
503 503 self.data_dc
504 504 self.flagIsNewFile
505 505 self.flagIsNewBlock
506 506 self.nTotalBlocks
507 507 self.nWriteBlocks
508 508
509 509 Return: None
510 510 """
511 511
512 512 spc = numpy.transpose( self.data_spc, (0,2,1) )
513 if not( self.processingHeaderObj.shif_fft ):
513 if self.processingHeaderObj.shif_fft:
514 514 spc = numpy.roll( spc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
515 515 data = spc.reshape((-1))
516 516 data = data.astype(self.dtype[0])
517 517 data.tofile(self.fp)
518 518
519 519 if self.data_cspc is not None:
520 520 data = numpy.zeros( self.shape_cspc_Buffer, self.dtype )
521 521 cspc = numpy.transpose( self.data_cspc, (0,2,1) )
522 if not( self.processingHeaderObj.shif_fft ):
522 if self.processingHeaderObj.shif_fft:
523 523 cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
524 524 data['real'] = cspc.real
525 525 data['imag'] = cspc.imag
526 526 data = data.reshape((-1))
527 527 data.tofile(self.fp)
528 528
529 529 if self.data_dc is not None:
530 530 data = numpy.zeros( self.shape_dc_Buffer, self.dtype )
531 531 dc = self.data_dc
532 532 data['real'] = dc.real
533 533 data['imag'] = dc.imag
534 534 data = data.reshape((-1))
535 535 data.tofile(self.fp)
536 536
537 537 # self.data_spc.fill(0)
538 538 #
539 539 # if self.data_dc is not None:
540 540 # self.data_dc.fill(0)
541 541 #
542 542 # if self.data_cspc is not None:
543 543 # self.data_cspc.fill(0)
544 544
545 545 self.flagIsNewFile = 0
546 546 self.flagIsNewBlock = 1
547 547 self.nTotalBlocks += 1
548 548 self.nWriteBlocks += 1
549 549 self.blockIndex += 1
550 550
551 551 # print "[Writing] Block = %d04" %self.blockIndex
552 552
553 553 def putData(self):
554 554 """
555 555 Setea un bloque de datos y luego los escribe en un file
556 556
557 557 Affected:
558 558 self.data_spc
559 559 self.data_cspc
560 560 self.data_dc
561 561
562 562 Return:
563 563 0 : Si no hay data o no hay mas files que puedan escribirse
564 564 1 : Si se escribio la data de un bloque en un file
565 565 """
566 566
567 567 if self.dataOut.flagNoData:
568 568 return 0
569 569
570 570 self.flagIsNewBlock = 0
571 571
572 572 if self.dataOut.flagDiscontinuousBlock:
573 573 self.data_spc.fill(0)
574 574 if self.dataOut.data_cspc is not None:
575 575 self.data_cspc.fill(0)
576 576 if self.dataOut.data_dc is not None:
577 577 self.data_dc.fill(0)
578 578 self.setNextFile()
579 579
580 580 if self.flagIsNewFile == 0:
581 581 self.setBasicHeader()
582 582
583 583 self.data_spc = self.dataOut.data_spc.copy()
584 584
585 585 if self.dataOut.data_cspc is not None:
586 586 self.data_cspc = self.dataOut.data_cspc.copy()
587 587
588 588 if self.dataOut.data_dc is not None:
589 589 self.data_dc = self.dataOut.data_dc.copy()
590 590
591 591 # #self.processingHeaderObj.dataBlocksPerFile)
592 592 if self.hasAllDataInBuffer():
593 593 # self.setFirstHeader()
594 594 self.writeNextBlock()
595 595
596 596 return 1
597 597
598 598 def __getBlockSize(self):
599 599 '''
600 600 Este metodos determina el cantidad de bytes para un bloque de datos de tipo Spectra
601 601 '''
602 602
603 603 dtype_width = self.getDtypeWidth()
604 604
605 605 pts2write = self.dataOut.nHeights * self.dataOut.nFFTPoints
606 606
607 607 pts2write_SelfSpectra = int(self.dataOut.nChannels * pts2write)
608 608 blocksize = (pts2write_SelfSpectra*dtype_width)
609 609
610 610 if self.dataOut.data_cspc is not None:
611 611 pts2write_CrossSpectra = int(self.dataOut.nPairs * pts2write)
612 612 blocksize += (pts2write_CrossSpectra*dtype_width*2)
613 613
614 614 if self.dataOut.data_dc is not None:
615 615 pts2write_DCchannels = int(self.dataOut.nChannels * self.dataOut.nHeights)
616 616 blocksize += (pts2write_DCchannels*dtype_width*2)
617 617
618 618 # blocksize = blocksize #* datatypeValue * 2 #CORREGIR ESTO
619 619
620 620 return blocksize
621 621
622 622 def setFirstHeader(self):
623 623
624 624 """
625 625 Obtiene una copia del First Header
626 626
627 627 Affected:
628 628 self.systemHeaderObj
629 629 self.radarControllerHeaderObj
630 630 self.dtype
631 631
632 632 Return:
633 633 None
634 634 """
635 635
636 636 self.systemHeaderObj = self.dataOut.systemHeaderObj.copy()
637 637 self.systemHeaderObj.nChannels = self.dataOut.nChannels
638 638 self.radarControllerHeaderObj = self.dataOut.radarControllerHeaderObj.copy()
639 639
640 640 self.processingHeaderObj.dtype = 1 # Spectra
641 641 self.processingHeaderObj.blockSize = self.__getBlockSize()
642 642 self.processingHeaderObj.profilesPerBlock = self.dataOut.nFFTPoints
643 643 self.processingHeaderObj.dataBlocksPerFile = self.blocksPerFile
644 644 self.processingHeaderObj.nWindows = 1 #podria ser 1 o self.dataOut.processingHeaderObj.nWindows
645 645 self.processingHeaderObj.nCohInt = self.dataOut.nCohInt# Se requiere para determinar el valor de timeInterval
646 646 self.processingHeaderObj.nIncohInt = self.dataOut.nIncohInt
647 647 self.processingHeaderObj.totalSpectra = self.dataOut.nPairs + self.dataOut.nChannels
648 648 self.processingHeaderObj.shif_fft = self.dataOut.flagShiftFFT
649 649
650 650 if self.processingHeaderObj.totalSpectra > 0:
651 651 channelList = []
652 652 for channel in range(self.dataOut.nChannels):
653 653 channelList.append(channel)
654 654 channelList.append(channel)
655 655
656 656 pairsList = []
657 657 if self.dataOut.nPairs > 0:
658 658 for pair in self.dataOut.pairsList:
659 659 pairsList.append(pair[0])
660 660 pairsList.append(pair[1])
661 661
662 662 spectraComb = channelList + pairsList
663 663 spectraComb = numpy.array(spectraComb, dtype="u1")
664 664 self.processingHeaderObj.spectraComb = spectraComb
665 665
666 666 if self.dataOut.code is not None:
667 667 self.processingHeaderObj.code = self.dataOut.code
668 668 self.processingHeaderObj.nCode = self.dataOut.nCode
669 669 self.processingHeaderObj.nBaud = self.dataOut.nBaud
670 670
671 671 if self.processingHeaderObj.nWindows != 0:
672 672 self.processingHeaderObj.firstHeight = self.dataOut.heightList[0]
673 673 self.processingHeaderObj.deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
674 674 self.processingHeaderObj.nHeights = self.dataOut.nHeights
675 675 self.processingHeaderObj.samplesWin = self.dataOut.nHeights
676 676
677 677 self.processingHeaderObj.processFlags = self.getProcessFlags()
678 678
679 679 self.setBasicHeader()
Show file before | Show file after | Hide comments
@@ -1,943 +1,952
1 1 import itertools
2 2
3 3 import numpy
4 4
5 5 from jroproc_base import ProcessingUnit, Operation
6 6 from schainpy.model.data.jrodata import Spectra
7 7 from schainpy.model.data.jrodata import hildebrand_sekhon
8 8
9 9
10 10 class SpectraProc(ProcessingUnit):
11 11
12 12 def __init__(self, **kwargs):
13 13
14 14 ProcessingUnit.__init__(self, **kwargs)
15 15
16 16 self.buffer = None
17 17 self.firstdatatime = None
18 18 self.profIndex = 0
19 19 self.dataOut = Spectra()
20 20 self.id_min = None
21 21 self.id_max = None
22 22
23 23 def __updateSpecFromVoltage(self):
24 24
25 25 self.dataOut.timeZone = self.dataIn.timeZone
26 26 self.dataOut.dstFlag = self.dataIn.dstFlag
27 27 self.dataOut.errorCount = self.dataIn.errorCount
28 28 self.dataOut.useLocalTime = self.dataIn.useLocalTime
29 29 try:
30 30 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
31 31 except:
32 32 pass
33 33 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
34 34 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
35 35 self.dataOut.channelList = self.dataIn.channelList
36 36 self.dataOut.heightList = self.dataIn.heightList
37 37 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
38 38
39 39 self.dataOut.nBaud = self.dataIn.nBaud
40 40 self.dataOut.nCode = self.dataIn.nCode
41 41 self.dataOut.code = self.dataIn.code
42 42 self.dataOut.nProfiles = self.dataOut.nFFTPoints
43 43
44 44 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
45 45 self.dataOut.utctime = self.firstdatatime
46 46 # asumo q la data esta decodificada
47 47 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
48 48 # asumo q la data esta sin flip
49 49 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
50 50 self.dataOut.flagShiftFFT = False
51 51
52 52 self.dataOut.nCohInt = self.dataIn.nCohInt
53 53 self.dataOut.nIncohInt = 1
54 54
55 55 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
56 56
57 57 self.dataOut.frequency = self.dataIn.frequency
58 58 self.dataOut.realtime = self.dataIn.realtime
59 59
60 60 self.dataOut.azimuth = self.dataIn.azimuth
61 61 self.dataOut.zenith = self.dataIn.zenith
62 62
63 63 self.dataOut.beam.codeList = self.dataIn.beam.codeList
64 64 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
65 65 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
66 66
67 67 def __getFft(self):
68 68 """
69 69 Convierte valores de Voltaje a Spectra
70 70
71 71 Affected:
72 72 self.dataOut.data_spc
73 73 self.dataOut.data_cspc
74 74 self.dataOut.data_dc
75 75 self.dataOut.heightList
76 76 self.profIndex
77 77 self.buffer
78 78 self.dataOut.flagNoData
79 79 """
80 80 fft_volt = numpy.fft.fft(
81 81 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
82 82 fft_volt = fft_volt.astype(numpy.dtype('complex'))
83 83 dc = fft_volt[:, 0, :]
84 84
85 85 # calculo de self-spectra
86 86 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
87 87 spc = fft_volt * numpy.conjugate(fft_volt)
88 88 spc = spc.real
89 89
90 90 blocksize = 0
91 91 blocksize += dc.size
92 92 blocksize += spc.size
93 93
94 94 cspc = None
95 95 pairIndex = 0
96 96 if self.dataOut.pairsList != None:
97 97 # calculo de cross-spectra
98 98 cspc = numpy.zeros(
99 99 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
100 100 for pair in self.dataOut.pairsList:
101 101 if pair[0] not in self.dataOut.channelList:
102 102 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
103 103 str(pair), str(self.dataOut.channelList))
104 104 if pair[1] not in self.dataOut.channelList:
105 105 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
106 106 str(pair), str(self.dataOut.channelList))
107 107
108 108 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
109 109 numpy.conjugate(fft_volt[pair[1], :, :])
110 110 pairIndex += 1
111 111 blocksize += cspc.size
112 112
113 113 self.dataOut.data_spc = spc
114 114 self.dataOut.data_cspc = cspc
115 115 self.dataOut.data_dc = dc
116 116 self.dataOut.blockSize = blocksize
117 117 self.dataOut.flagShiftFFT = True
118 118
119 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None):
119 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None, shift_fft=False):
120 120
121 121 self.dataOut.flagNoData = True
122 122
123 123 if self.dataIn.type == "Spectra":
124 124 self.dataOut.copy(self.dataIn)
125 125 if not pairsList:
126 126 pairsList = itertools.combinations(self.dataOut.channelList, 2)
127 127 if self.dataOut.data_cspc is not None:
128 128 self.__selectPairs(pairsList)
129 if shift_fft:
130 #desplaza a la derecha en el eje 2 determinadas posiciones
131 shift = int(self.dataOut.nFFTPoints/2)
132 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
133
134 if self.dataOut.data_cspc is not None:
135 #desplaza a la derecha en el eje 2 determinadas posiciones
136 self.dataOut.data_cspc = numpy.roll(self.dataOut.cspc, shift, axis=1)
137
129 138 return True
130 139
131 140 if self.dataIn.type == "Voltage":
132 141
133 142 if nFFTPoints == None:
134 143 raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
135 144
136 145 if nProfiles == None:
137 146 nProfiles = nFFTPoints
138 147
139 148 if ippFactor == None:
140 149 ippFactor = 1
141 150
142 151 self.dataOut.ippFactor = ippFactor
143 152
144 153 self.dataOut.nFFTPoints = nFFTPoints
145 154 self.dataOut.pairsList = pairsList
146 155
147 156 if self.buffer is None:
148 157 self.buffer = numpy.zeros((self.dataIn.nChannels,
149 158 nProfiles,
150 159 self.dataIn.nHeights),
151 160 dtype='complex')
152 161
153 162 if self.dataIn.flagDataAsBlock:
154 163 # data dimension: [nChannels, nProfiles, nSamples]
155 164 nVoltProfiles = self.dataIn.data.shape[1]
156 165 # nVoltProfiles = self.dataIn.nProfiles
157 166
158 167 if nVoltProfiles == nProfiles:
159 168 self.buffer = self.dataIn.data.copy()
160 169 self.profIndex = nVoltProfiles
161 170
162 171 elif nVoltProfiles < nProfiles:
163 172
164 173 if self.profIndex == 0:
165 174 self.id_min = 0
166 175 self.id_max = nVoltProfiles
167 176
168 177 self.buffer[:, self.id_min:self.id_max,
169 178 :] = self.dataIn.data
170 179 self.profIndex += nVoltProfiles
171 180 self.id_min += nVoltProfiles
172 181 self.id_max += nVoltProfiles
173 182 else:
174 183 raise ValueError, "The type object %s has %d profiles, it should just has %d profiles" % (
175 184 self.dataIn.type, self.dataIn.data.shape[1], nProfiles)
176 185 self.dataOut.flagNoData = True
177 186 return 0
178 187 else:
179 188 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
180 189 self.profIndex += 1
181 190
182 191 if self.firstdatatime == None:
183 192 self.firstdatatime = self.dataIn.utctime
184 193
185 194 if self.profIndex == nProfiles:
186 195 self.__updateSpecFromVoltage()
187 196 self.__getFft()
188 197
189 198 self.dataOut.flagNoData = False
190 199 self.firstdatatime = None
191 200 self.profIndex = 0
192 201
193 202 return True
194 203
195 204 raise ValueError, "The type of input object '%s' is not valid" % (
196 205 self.dataIn.type)
197 206
198 207 def __selectPairs(self, pairsList):
199 208
200 209 if not pairsList:
201 210 return
202 211
203 212 pairs = []
204 213 pairsIndex = []
205 214
206 215 for pair in pairsList:
207 216 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
208 217 continue
209 218 pairs.append(pair)
210 219 pairsIndex.append(pairs.index(pair))
211 220
212 221 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
213 222 self.dataOut.pairsList = pairs
214 223
215 224 return
216 225
217 226 def __selectPairsByChannel(self, channelList=None):
218 227
219 228 if channelList == None:
220 229 return
221 230
222 231 pairsIndexListSelected = []
223 232 for pairIndex in self.dataOut.pairsIndexList:
224 233 # First pair
225 234 if self.dataOut.pairsList[pairIndex][0] not in channelList:
226 235 continue
227 236 # Second pair
228 237 if self.dataOut.pairsList[pairIndex][1] not in channelList:
229 238 continue
230 239
231 240 pairsIndexListSelected.append(pairIndex)
232 241
233 242 if not pairsIndexListSelected:
234 243 self.dataOut.data_cspc = None
235 244 self.dataOut.pairsList = []
236 245 return
237 246
238 247 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
239 248 self.dataOut.pairsList = [self.dataOut.pairsList[i]
240 249 for i in pairsIndexListSelected]
241 250
242 251 return
243 252
244 253 def selectChannels(self, channelList):
245 254
246 255 channelIndexList = []
247 256
248 257 for channel in channelList:
249 258 if channel not in self.dataOut.channelList:
250 259 raise ValueError, "Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" % (
251 260 channel, str(self.dataOut.channelList))
252 261
253 262 index = self.dataOut.channelList.index(channel)
254 263 channelIndexList.append(index)
255 264
256 265 self.selectChannelsByIndex(channelIndexList)
257 266
258 267 def selectChannelsByIndex(self, channelIndexList):
259 268 """
260 269 Selecciona un bloque de datos en base a canales segun el channelIndexList
261 270
262 271 Input:
263 272 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
264 273
265 274 Affected:
266 275 self.dataOut.data_spc
267 276 self.dataOut.channelIndexList
268 277 self.dataOut.nChannels
269 278
270 279 Return:
271 280 None
272 281 """
273 282
274 283 for channelIndex in channelIndexList:
275 284 if channelIndex not in self.dataOut.channelIndexList:
276 285 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " % (
277 286 channelIndex, self.dataOut.channelIndexList)
278 287
279 288 # nChannels = len(channelIndexList)
280 289
281 290 data_spc = self.dataOut.data_spc[channelIndexList, :]
282 291 data_dc = self.dataOut.data_dc[channelIndexList, :]
283 292
284 293 self.dataOut.data_spc = data_spc
285 294 self.dataOut.data_dc = data_dc
286 295
287 296 self.dataOut.channelList = [
288 297 self.dataOut.channelList[i] for i in channelIndexList]
289 298 # self.dataOut.nChannels = nChannels
290 299
291 300 self.__selectPairsByChannel(self.dataOut.channelList)
292 301
293 302 return 1
294 303
295 304 def selectHeights(self, minHei, maxHei):
296 305 """
297 306 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
298 307 minHei <= height <= maxHei
299 308
300 309 Input:
301 310 minHei : valor minimo de altura a considerar
302 311 maxHei : valor maximo de altura a considerar
303 312
304 313 Affected:
305 314 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
306 315
307 316 Return:
308 317 1 si el metodo se ejecuto con exito caso contrario devuelve 0
309 318 """
310 319
311 320 if (minHei > maxHei):
312 321 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (
313 322 minHei, maxHei)
314 323
315 324 if (minHei < self.dataOut.heightList[0]):
316 325 minHei = self.dataOut.heightList[0]
317 326
318 327 if (maxHei > self.dataOut.heightList[-1]):
319 328 maxHei = self.dataOut.heightList[-1]
320 329
321 330 minIndex = 0
322 331 maxIndex = 0
323 332 heights = self.dataOut.heightList
324 333
325 334 inda = numpy.where(heights >= minHei)
326 335 indb = numpy.where(heights <= maxHei)
327 336
328 337 try:
329 338 minIndex = inda[0][0]
330 339 except:
331 340 minIndex = 0
332 341
333 342 try:
334 343 maxIndex = indb[0][-1]
335 344 except:
336 345 maxIndex = len(heights)
337 346
338 347 self.selectHeightsByIndex(minIndex, maxIndex)
339 348
340 349 return 1
341 350
342 351 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
343 352 newheis = numpy.where(
344 353 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
345 354
346 355 if hei_ref != None:
347 356 newheis = numpy.where(self.dataOut.heightList > hei_ref)
348 357
349 358 minIndex = min(newheis[0])
350 359 maxIndex = max(newheis[0])
351 360 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
352 361 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
353 362
354 363 # determina indices
355 364 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
356 365 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
357 366 avg_dB = 10 * \
358 367 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
359 368 beacon_dB = numpy.sort(avg_dB)[-nheis:]
360 369 beacon_heiIndexList = []
361 370 for val in avg_dB.tolist():
362 371 if val >= beacon_dB[0]:
363 372 beacon_heiIndexList.append(avg_dB.tolist().index(val))
364 373
365 374 #data_spc = data_spc[:,:,beacon_heiIndexList]
366 375 data_cspc = None
367 376 if self.dataOut.data_cspc is not None:
368 377 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
369 378 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
370 379
371 380 data_dc = None
372 381 if self.dataOut.data_dc is not None:
373 382 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
374 383 #data_dc = data_dc[:,beacon_heiIndexList]
375 384
376 385 self.dataOut.data_spc = data_spc
377 386 self.dataOut.data_cspc = data_cspc
378 387 self.dataOut.data_dc = data_dc
379 388 self.dataOut.heightList = heightList
380 389 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
381 390
382 391 return 1
383 392
384 393 def selectHeightsByIndex(self, minIndex, maxIndex):
385 394 """
386 395 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
387 396 minIndex <= index <= maxIndex
388 397
389 398 Input:
390 399 minIndex : valor de indice minimo de altura a considerar
391 400 maxIndex : valor de indice maximo de altura a considerar
392 401
393 402 Affected:
394 403 self.dataOut.data_spc
395 404 self.dataOut.data_cspc
396 405 self.dataOut.data_dc
397 406 self.dataOut.heightList
398 407
399 408 Return:
400 409 1 si el metodo se ejecuto con exito caso contrario devuelve 0
401 410 """
402 411
403 412 if (minIndex < 0) or (minIndex > maxIndex):
404 413 raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (
405 414 minIndex, maxIndex)
406 415
407 416 if (maxIndex >= self.dataOut.nHeights):
408 417 maxIndex = self.dataOut.nHeights - 1
409 418
410 419 # Spectra
411 420 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
412 421
413 422 data_cspc = None
414 423 if self.dataOut.data_cspc is not None:
415 424 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
416 425
417 426 data_dc = None
418 427 if self.dataOut.data_dc is not None:
419 428 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
420 429
421 430 self.dataOut.data_spc = data_spc
422 431 self.dataOut.data_cspc = data_cspc
423 432 self.dataOut.data_dc = data_dc
424 433
425 434 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
426 435
427 436 return 1
428 437
429 438 def removeDC(self, mode=2):
430 439 jspectra = self.dataOut.data_spc
431 440 jcspectra = self.dataOut.data_cspc
432 441
433 442 num_chan = jspectra.shape[0]
434 443 num_hei = jspectra.shape[2]
435 444
436 445 if jcspectra is not None:
437 446 jcspectraExist = True
438 447 num_pairs = jcspectra.shape[0]
439 448 else:
440 449 jcspectraExist = False
441 450
442 451 freq_dc = jspectra.shape[1] / 2
443 452 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
444 453
445 454 if ind_vel[0] < 0:
446 455 ind_vel[range(0, 1)] = ind_vel[range(0, 1)] + self.num_prof
447 456
448 457 if mode == 1:
449 458 jspectra[:, freq_dc, :] = (
450 459 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
451 460
452 461 if jcspectraExist:
453 462 jcspectra[:, freq_dc, :] = (
454 463 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
455 464
456 465 if mode == 2:
457 466
458 467 vel = numpy.array([-2, -1, 1, 2])
459 468 xx = numpy.zeros([4, 4])
460 469
461 470 for fil in range(4):
462 471 xx[fil, :] = vel[fil]**numpy.asarray(range(4))
463 472
464 473 xx_inv = numpy.linalg.inv(xx)
465 474 xx_aux = xx_inv[0, :]
466 475
467 476 for ich in range(num_chan):
468 477 yy = jspectra[ich, ind_vel, :]
469 478 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
470 479
471 480 junkid = jspectra[ich, freq_dc, :] <= 0
472 481 cjunkid = sum(junkid)
473 482
474 483 if cjunkid.any():
475 484 jspectra[ich, freq_dc, junkid.nonzero()] = (
476 485 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
477 486
478 487 if jcspectraExist:
479 488 for ip in range(num_pairs):
480 489 yy = jcspectra[ip, ind_vel, :]
481 490 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
482 491
483 492 self.dataOut.data_spc = jspectra
484 493 self.dataOut.data_cspc = jcspectra
485 494
486 495 return 1
487 496
488 497 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
489 498
490 499 jspectra = self.dataOut.data_spc
491 500 jcspectra = self.dataOut.data_cspc
492 501 jnoise = self.dataOut.getNoise()
493 502 num_incoh = self.dataOut.nIncohInt
494 503
495 504 num_channel = jspectra.shape[0]
496 505 num_prof = jspectra.shape[1]
497 506 num_hei = jspectra.shape[2]
498 507
499 508 # hei_interf
500 509 if hei_interf is None:
501 510 count_hei = num_hei / 2 # Como es entero no importa
502 511 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
503 512 hei_interf = numpy.asarray(hei_interf)[0]
504 513 # nhei_interf
505 514 if (nhei_interf == None):
506 515 nhei_interf = 5
507 516 if (nhei_interf < 1):
508 517 nhei_interf = 1
509 518 if (nhei_interf > count_hei):
510 519 nhei_interf = count_hei
511 520 if (offhei_interf == None):
512 521 offhei_interf = 0
513 522
514 523 ind_hei = range(num_hei)
515 524 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
516 525 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
517 526 mask_prof = numpy.asarray(range(num_prof))
518 527 num_mask_prof = mask_prof.size
519 528 comp_mask_prof = [0, num_prof / 2]
520 529
521 530 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
522 531 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
523 532 jnoise = numpy.nan
524 533 noise_exist = jnoise[0] < numpy.Inf
525 534
526 535 # Subrutina de Remocion de la Interferencia
527 536 for ich in range(num_channel):
528 537 # Se ordena los espectros segun su potencia (menor a mayor)
529 538 power = jspectra[ich, mask_prof, :]
530 539 power = power[:, hei_interf]
531 540 power = power.sum(axis=0)
532 541 psort = power.ravel().argsort()
533 542
534 543 # Se estima la interferencia promedio en los Espectros de Potencia empleando
535 544 junkspc_interf = jspectra[ich, :, hei_interf[psort[range(
536 545 offhei_interf, nhei_interf + offhei_interf)]]]
537 546
538 547 if noise_exist:
539 548 # tmp_noise = jnoise[ich] / num_prof
540 549 tmp_noise = jnoise[ich]
541 550 junkspc_interf = junkspc_interf - tmp_noise
542 551 #junkspc_interf[:,comp_mask_prof] = 0
543 552
544 553 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
545 554 jspc_interf = jspc_interf.transpose()
546 555 # Calculando el espectro de interferencia promedio
547 556 noiseid = numpy.where(
548 557 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
549 558 noiseid = noiseid[0]
550 559 cnoiseid = noiseid.size
551 560 interfid = numpy.where(
552 561 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
553 562 interfid = interfid[0]
554 563 cinterfid = interfid.size
555 564
556 565 if (cnoiseid > 0):
557 566 jspc_interf[noiseid] = 0
558 567
559 568 # Expandiendo los perfiles a limpiar
560 569 if (cinterfid > 0):
561 570 new_interfid = (
562 571 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
563 572 new_interfid = numpy.asarray(new_interfid)
564 573 new_interfid = {x for x in new_interfid}
565 574 new_interfid = numpy.array(list(new_interfid))
566 575 new_cinterfid = new_interfid.size
567 576 else:
568 577 new_cinterfid = 0
569 578
570 579 for ip in range(new_cinterfid):
571 580 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
572 581 jspc_interf[new_interfid[ip]
573 582 ] = junkspc_interf[ind[nhei_interf / 2], new_interfid[ip]]
574 583
575 584 jspectra[ich, :, ind_hei] = jspectra[ich, :,
576 585 ind_hei] - jspc_interf # Corregir indices
577 586
578 587 # Removiendo la interferencia del punto de mayor interferencia
579 588 ListAux = jspc_interf[mask_prof].tolist()
580 589 maxid = ListAux.index(max(ListAux))
581 590
582 591 if cinterfid > 0:
583 592 for ip in range(cinterfid * (interf == 2) - 1):
584 593 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
585 594 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
586 595 cind = len(ind)
587 596
588 597 if (cind > 0):
589 598 jspectra[ich, interfid[ip], ind] = tmp_noise * \
590 599 (1 + (numpy.random.uniform(cind) - 0.5) /
591 600 numpy.sqrt(num_incoh))
592 601
593 602 ind = numpy.array([-2, -1, 1, 2])
594 603 xx = numpy.zeros([4, 4])
595 604
596 605 for id1 in range(4):
597 606 xx[:, id1] = ind[id1]**numpy.asarray(range(4))
598 607
599 608 xx_inv = numpy.linalg.inv(xx)
600 609 xx = xx_inv[:, 0]
601 610 ind = (ind + maxid + num_mask_prof) % num_mask_prof
602 611 yy = jspectra[ich, mask_prof[ind], :]
603 612 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
604 613 yy.transpose(), xx)
605 614
606 615 indAux = (jspectra[ich, :, :] < tmp_noise *
607 616 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
608 617 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
609 618 (1 - 1 / numpy.sqrt(num_incoh))
610 619
611 620 # Remocion de Interferencia en el Cross Spectra
612 621 if jcspectra is None:
613 622 return jspectra, jcspectra
614 623 num_pairs = jcspectra.size / (num_prof * num_hei)
615 624 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
616 625
617 626 for ip in range(num_pairs):
618 627
619 628 #-------------------------------------------
620 629
621 630 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
622 631 cspower = cspower[:, hei_interf]
623 632 cspower = cspower.sum(axis=0)
624 633
625 634 cspsort = cspower.ravel().argsort()
626 635 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[range(
627 636 offhei_interf, nhei_interf + offhei_interf)]]]
628 637 junkcspc_interf = junkcspc_interf.transpose()
629 638 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
630 639
631 640 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
632 641
633 642 median_real = numpy.median(numpy.real(
634 643 junkcspc_interf[mask_prof[ind[range(3 * num_prof / 4)]], :]))
635 644 median_imag = numpy.median(numpy.imag(
636 645 junkcspc_interf[mask_prof[ind[range(3 * num_prof / 4)]], :]))
637 646 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
638 647 median_real, median_imag)
639 648
640 649 for iprof in range(num_prof):
641 650 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
642 651 jcspc_interf[iprof] = junkcspc_interf[iprof,
643 652 ind[nhei_interf / 2]]
644 653
645 654 # Removiendo la Interferencia
646 655 jcspectra[ip, :, ind_hei] = jcspectra[ip,
647 656 :, ind_hei] - jcspc_interf
648 657
649 658 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
650 659 maxid = ListAux.index(max(ListAux))
651 660
652 661 ind = numpy.array([-2, -1, 1, 2])
653 662 xx = numpy.zeros([4, 4])
654 663
655 664 for id1 in range(4):
656 665 xx[:, id1] = ind[id1]**numpy.asarray(range(4))
657 666
658 667 xx_inv = numpy.linalg.inv(xx)
659 668 xx = xx_inv[:, 0]
660 669
661 670 ind = (ind + maxid + num_mask_prof) % num_mask_prof
662 671 yy = jcspectra[ip, mask_prof[ind], :]
663 672 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
664 673
665 674 # Guardar Resultados
666 675 self.dataOut.data_spc = jspectra
667 676 self.dataOut.data_cspc = jcspectra
668 677
669 678 return 1
670 679
671 680 def setRadarFrequency(self, frequency=None):
672 681
673 682 if frequency != None:
674 683 self.dataOut.frequency = frequency
675 684
676 685 return 1
677 686
678 687 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
679 688 # validacion de rango
680 689 if minHei == None:
681 690 minHei = self.dataOut.heightList[0]
682 691
683 692 if maxHei == None:
684 693 maxHei = self.dataOut.heightList[-1]
685 694
686 695 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
687 696 print 'minHei: %.2f is out of the heights range' % (minHei)
688 697 print 'minHei is setting to %.2f' % (self.dataOut.heightList[0])
689 698 minHei = self.dataOut.heightList[0]
690 699
691 700 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
692 701 print 'maxHei: %.2f is out of the heights range' % (maxHei)
693 702 print 'maxHei is setting to %.2f' % (self.dataOut.heightList[-1])
694 703 maxHei = self.dataOut.heightList[-1]
695 704
696 705 # validacion de velocidades
697 706 velrange = self.dataOut.getVelRange(1)
698 707
699 708 if minVel == None:
700 709 minVel = velrange[0]
701 710
702 711 if maxVel == None:
703 712 maxVel = velrange[-1]
704 713
705 714 if (minVel < velrange[0]) or (minVel > maxVel):
706 715 print 'minVel: %.2f is out of the velocity range' % (minVel)
707 716 print 'minVel is setting to %.2f' % (velrange[0])
708 717 minVel = velrange[0]
709 718
710 719 if (maxVel > velrange[-1]) or (maxVel < minVel):
711 720 print 'maxVel: %.2f is out of the velocity range' % (maxVel)
712 721 print 'maxVel is setting to %.2f' % (velrange[-1])
713 722 maxVel = velrange[-1]
714 723
715 724 # seleccion de indices para rango
716 725 minIndex = 0
717 726 maxIndex = 0
718 727 heights = self.dataOut.heightList
719 728
720 729 inda = numpy.where(heights >= minHei)
721 730 indb = numpy.where(heights <= maxHei)
722 731
723 732 try:
724 733 minIndex = inda[0][0]
725 734 except:
726 735 minIndex = 0
727 736
728 737 try:
729 738 maxIndex = indb[0][-1]
730 739 except:
731 740 maxIndex = len(heights)
732 741
733 742 if (minIndex < 0) or (minIndex > maxIndex):
734 743 raise ValueError, "some value in (%d,%d) is not valid" % (
735 744 minIndex, maxIndex)
736 745
737 746 if (maxIndex >= self.dataOut.nHeights):
738 747 maxIndex = self.dataOut.nHeights - 1
739 748
740 749 # seleccion de indices para velocidades
741 750 indminvel = numpy.where(velrange >= minVel)
742 751 indmaxvel = numpy.where(velrange <= maxVel)
743 752 try:
744 753 minIndexVel = indminvel[0][0]
745 754 except:
746 755 minIndexVel = 0
747 756
748 757 try:
749 758 maxIndexVel = indmaxvel[0][-1]
750 759 except:
751 760 maxIndexVel = len(velrange)
752 761
753 762 # seleccion del espectro
754 763 data_spc = self.dataOut.data_spc[:,
755 764 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
756 765 # estimacion de ruido
757 766 noise = numpy.zeros(self.dataOut.nChannels)
758 767
759 768 for channel in range(self.dataOut.nChannels):
760 769 daux = data_spc[channel, :, :]
761 770 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
762 771
763 772 self.dataOut.noise_estimation = noise.copy()
764 773
765 774 return 1
766 775
767 776
768 777 class IncohInt(Operation):
769 778
770 779 __profIndex = 0
771 780 __withOverapping = False
772 781
773 782 __byTime = False
774 783 __initime = None
775 784 __lastdatatime = None
776 785 __integrationtime = None
777 786
778 787 __buffer_spc = None
779 788 __buffer_cspc = None
780 789 __buffer_dc = None
781 790
782 791 __dataReady = False
783 792
784 793 __timeInterval = None
785 794
786 795 n = None
787 796
788 797 def __init__(self, **kwargs):
789 798
790 799 Operation.__init__(self, **kwargs)
791 800 # self.isConfig = False
792 801
793 802 def setup(self, n=None, timeInterval=None, overlapping=False):
794 803 """
795 804 Set the parameters of the integration class.
796 805
797 806 Inputs:
798 807
799 808 n : Number of coherent integrations
800 809 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
801 810 overlapping :
802 811
803 812 """
804 813
805 814 self.__initime = None
806 815 self.__lastdatatime = 0
807 816
808 817 self.__buffer_spc = 0
809 818 self.__buffer_cspc = 0
810 819 self.__buffer_dc = 0
811 820
812 821 self.__profIndex = 0
813 822 self.__dataReady = False
814 823 self.__byTime = False
815 824
816 825 if n is None and timeInterval is None:
817 826 raise ValueError, "n or timeInterval should be specified ..."
818 827
819 828 if n is not None:
820 829 self.n = int(n)
821 830 else:
822 831 # if (type(timeInterval)!=integer) -> change this line
823 832 self.__integrationtime = int(timeInterval)
824 833 self.n = None
825 834 self.__byTime = True
826 835
827 836 def putData(self, data_spc, data_cspc, data_dc):
828 837 """
829 838 Add a profile to the __buffer_spc and increase in one the __profileIndex
830 839
831 840 """
832 841
833 842 self.__buffer_spc += data_spc
834 843
835 844 if data_cspc is None:
836 845 self.__buffer_cspc = None
837 846 else:
838 847 self.__buffer_cspc += data_cspc
839 848
840 849 if data_dc is None:
841 850 self.__buffer_dc = None
842 851 else:
843 852 self.__buffer_dc += data_dc
844 853
845 854 self.__profIndex += 1
846 855
847 856 return
848 857
849 858 def pushData(self):
850 859 """
851 860 Return the sum of the last profiles and the profiles used in the sum.
852 861
853 862 Affected:
854 863
855 864 self.__profileIndex
856 865
857 866 """
858 867
859 868 data_spc = self.__buffer_spc
860 869 data_cspc = self.__buffer_cspc
861 870 data_dc = self.__buffer_dc
862 871 n = self.__profIndex
863 872
864 873 self.__buffer_spc = 0
865 874 self.__buffer_cspc = 0
866 875 self.__buffer_dc = 0
867 876 self.__profIndex = 0
868 877
869 878 return data_spc, data_cspc, data_dc, n
870 879
871 880 def byProfiles(self, *args):
872 881
873 882 self.__dataReady = False
874 883 avgdata_spc = None
875 884 avgdata_cspc = None
876 885 avgdata_dc = None
877 886
878 887 self.putData(*args)
879 888
880 889 if self.__profIndex == self.n:
881 890
882 891 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
883 892 self.n = n
884 893 self.__dataReady = True
885 894
886 895 return avgdata_spc, avgdata_cspc, avgdata_dc
887 896
888 897 def byTime(self, datatime, *args):
889 898
890 899 self.__dataReady = False
891 900 avgdata_spc = None
892 901 avgdata_cspc = None
893 902 avgdata_dc = None
894 903
895 904 self.putData(*args)
896 905
897 906 if (datatime - self.__initime) >= self.__integrationtime:
898 907 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
899 908 self.n = n
900 909 self.__dataReady = True
901 910
902 911 return avgdata_spc, avgdata_cspc, avgdata_dc
903 912
904 913 def integrate(self, datatime, *args):
905 914
906 915 if self.__profIndex == 0:
907 916 self.__initime = datatime
908 917
909 918 if self.__byTime:
910 919 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
911 920 datatime, *args)
912 921 else:
913 922 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
914 923
915 924 if not self.__dataReady:
916 925 return None, None, None, None
917 926
918 927 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
919 928
920 929 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
921 930 if n == 1:
922 931 return
923 932
924 933 dataOut.flagNoData = True
925 934
926 935 if not self.isConfig:
927 936 self.setup(n, timeInterval, overlapping)
928 937 self.isConfig = True
929 938
930 939 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
931 940 dataOut.data_spc,
932 941 dataOut.data_cspc,
933 942 dataOut.data_dc)
934 943
935 944 if self.__dataReady:
936 945
937 946 dataOut.data_spc = avgdata_spc
938 947 dataOut.data_cspc = avgdata_cspc
939 948 dataOut.data_dc = avgdata_dc
940 949
941 950 dataOut.nIncohInt *= self.n
942 951 dataOut.utctime = avgdatatime
943 952 dataOut.flagNoData = False
General Comments 0
You need to be logged in to leave comments. Login now