##// END OF EJS Templates
schain
RSS
update numpy.complex -> numpy.complex_ uso de libreria numpy mayor a 1.23
Alexander Valdez -
r1722:aebcee300801
parent child
Show More
Side by Side Unified
Context file:
r1722:aebcee300801 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1625 +1,1625
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 7 from time import time
8 8 # voltage proc master
9 9
10 10 class VoltageProc(ProcessingUnit):
11 11
12 12 def __init__(self):
13 13
14 14 ProcessingUnit.__init__(self)
15 15
16 16 self.dataOut = Voltage()
17 17 self.flip = 1
18 18 self.setupReq = False
19 19
20 20 def run(self, runNextUnit = 0):
21 21
22 22 if self.dataIn.type == 'AMISR':
23 23 self.__updateObjFromAmisrInput()
24 24
25 25 if self.dataIn.type == 'Voltage':
26 26 self.dataOut.copy(self.dataIn)
27 27 self.dataOut.runNextUnit = runNextUnit
28 28
29 29
30 30 def __updateObjFromAmisrInput(self):
31 31
32 32 self.dataOut.timeZone = self.dataIn.timeZone
33 33 self.dataOut.dstFlag = self.dataIn.dstFlag
34 34 self.dataOut.errorCount = self.dataIn.errorCount
35 35 self.dataOut.useLocalTime = self.dataIn.useLocalTime
36 36
37 37 self.dataOut.flagNoData = self.dataIn.flagNoData
38 38 self.dataOut.data = self.dataIn.data
39 39 self.dataOut.utctime = self.dataIn.utctime
40 40 self.dataOut.channelList = self.dataIn.channelList
41 41 #self.dataOut.timeInterval = self.dataIn.timeInterval
42 42 self.dataOut.heightList = self.dataIn.heightList
43 43 self.dataOut.nProfiles = self.dataIn.nProfiles
44 44
45 45 self.dataOut.nCohInt = self.dataIn.nCohInt
46 46 self.dataOut.ippSeconds = self.dataIn.ippSeconds
47 47 self.dataOut.frequency = self.dataIn.frequency
48 48
49 49 self.dataOut.azimuth = self.dataIn.azimuth
50 50 self.dataOut.zenith = self.dataIn.zenith
51 51
52 52 self.dataOut.beam.codeList = self.dataIn.beam.codeList
53 53 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
54 54 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
55 55
56 56
57 57 class selectChannels(Operation):
58 58
59 59 def run(self, dataOut, channelList):
60 60
61 61 channelIndexList = []
62 62 self.dataOut = dataOut
63 63 for channel in channelList:
64 64 if channel not in self.dataOut.channelList:
65 65 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
66 66
67 67 index = self.dataOut.channelList.index(channel)
68 68 channelIndexList.append(index)
69 69 self.selectChannelsByIndex(channelIndexList)
70 70 return self.dataOut
71 71
72 72 def selectChannelsByIndex(self, channelIndexList):
73 73 """
74 74 Selecciona un bloque de datos en base a canales segun el channelIndexList
75 75
76 76 Input:
77 77 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
78 78
79 79 Affected:
80 80 self.dataOut.data
81 81 self.dataOut.channelIndexList
82 82 self.dataOut.nChannels
83 83 self.dataOut.m_ProcessingHeader.totalSpectra
84 84 self.dataOut.systemHeaderObj.numChannels
85 85 self.dataOut.m_ProcessingHeader.blockSize
86 86
87 87 Return:
88 88 None
89 89 """
90 90
91 91 for channelIndex in channelIndexList:
92 92 if channelIndex not in self.dataOut.channelIndexList:
93 93 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
94 94
95 95 if self.dataOut.type == 'Voltage':
96 96 if self.dataOut.flagDataAsBlock:
97 97 """
98 98 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
99 99 """
100 100 data = self.dataOut.data[channelIndexList,:,:]
101 101 else:
102 102 data = self.dataOut.data[channelIndexList,:]
103 103
104 104 self.dataOut.data = data
105 105 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
106 106 self.dataOut.channelList = range(len(channelIndexList))
107 107
108 108 elif self.dataOut.type == 'Spectra':
109 109 data_spc = self.dataOut.data_spc[channelIndexList, :]
110 110 data_dc = self.dataOut.data_dc[channelIndexList, :]
111 111
112 112 self.dataOut.data_spc = data_spc
113 113 self.dataOut.data_dc = data_dc
114 114
115 115 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
116 116 self.dataOut.channelList = range(len(channelIndexList))
117 117 self.__selectPairsByChannel(channelIndexList)
118 118
119 119 return 1
120 120
121 121 def __selectPairsByChannel(self, channelList=None):
122 122
123 123 if channelList == None:
124 124 return
125 125
126 126 pairsIndexListSelected = []
127 127 for pairIndex in self.dataOut.pairsIndexList:
128 128 # First pair
129 129 if self.dataOut.pairsList[pairIndex][0] not in channelList:
130 130 continue
131 131 # Second pair
132 132 if self.dataOut.pairsList[pairIndex][1] not in channelList:
133 133 continue
134 134
135 135 pairsIndexListSelected.append(pairIndex)
136 136
137 137 if not pairsIndexListSelected:
138 138 self.dataOut.data_cspc = None
139 139 self.dataOut.pairsList = []
140 140 return
141 141
142 142 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
143 143 self.dataOut.pairsList = [self.dataOut.pairsList[i]
144 144 for i in pairsIndexListSelected]
145 145
146 146 return
147 147
148 148 class selectHeights(Operation):
149 149
150 150 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
151 151 """
152 152 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
153 153 minHei <= height <= maxHei
154 154
155 155 Input:
156 156 minHei : valor minimo de altura a considerar
157 157 maxHei : valor maximo de altura a considerar
158 158
159 159 Affected:
160 160 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
161 161
162 162 Return:
163 163 1 si el metodo se ejecuto con exito caso contrario devuelve 0
164 164 """
165 165
166 166 self.dataOut = dataOut
167 167
168 168 if type(minHei) == int or type(minHei) == float:
169 169 v_minHei= True
170 170 else:
171 171 v_minHei= False
172 172
173 173 if v_minHei and maxHei:
174 174 if (minHei < self.dataOut.heightList[0]):
175 175 minHei = self.dataOut.heightList[0]
176 176
177 177 if (maxHei > self.dataOut.heightList[-1]):
178 178 maxHei = self.dataOut.heightList[-1]
179 179
180 180 minIndex = 0
181 181 maxIndex = 0
182 182 heights = self.dataOut.heightList
183 183 inda = numpy.where(heights >= minHei)
184 184 indb = numpy.where(heights <= maxHei)
185 185
186 186 try:
187 187 minIndex = inda[0][0]
188 188 except:
189 189 minIndex = 0
190 190
191 191 try:
192 192 maxIndex = indb[0][-1]
193 193 except:
194 194 maxIndex = len(heights)
195 195 self.selectHeightsByIndex(minIndex, maxIndex)
196 196
197 197 return self.dataOut
198 198
199 199 def selectHeightsByIndex(self, minIndex, maxIndex):
200 200 """
201 201 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
202 202 minIndex <= index <= maxIndex
203 203
204 204 Input:
205 205 minIndex : valor de indice minimo de altura a considerar
206 206 maxIndex : valor de indice maximo de altura a considerar
207 207
208 208 Affected:
209 209 self.dataOut.data
210 210 self.dataOut.heightList
211 211
212 212 Return:
213 213 1 si el metodo se ejecuto con exito caso contrario devuelve 0
214 214 """
215 215
216 216 if self.dataOut.type == 'Voltage':
217 217 if (minIndex < 0) or (minIndex > maxIndex):
218 218 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
219 219
220 220 if (maxIndex >= self.dataOut.nHeights):
221 221 maxIndex = self.dataOut.nHeights
222 222
223 223 #voltage
224 224 if self.dataOut.flagDataAsBlock:
225 225 """
226 226 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
227 227 """
228 228 data = self.dataOut.data[:,:, minIndex:maxIndex]
229 229 else:
230 230 data = self.dataOut.data[:, minIndex:maxIndex]
231 231
232 232 # firstHeight = self.dataOut.heightList[minIndex]
233 233
234 234 self.dataOut.data = data
235 235 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
236 236
237 237 if self.dataOut.nHeights <= 1:
238 238 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
239 239 elif self.dataOut.type == 'Spectra':
240 240 if (minIndex < 0) or (minIndex > maxIndex):
241 241 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
242 242 minIndex, maxIndex))
243 243
244 244 if (maxIndex >= self.dataOut.nHeights):
245 245 maxIndex = self.dataOut.nHeights - 1
246 246
247 247 # Spectra
248 248 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
249 249
250 250 data_cspc = None
251 251 if self.dataOut.data_cspc is not None:
252 252 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
253 253
254 254 data_dc = None
255 255 if self.dataOut.data_dc is not None:
256 256 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
257 257
258 258 self.dataOut.data_spc = data_spc
259 259 self.dataOut.data_cspc = data_cspc
260 260 self.dataOut.data_dc = data_dc
261 261
262 262 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
263 263
264 264 return 1
265 265
266 266
267 267 class filterByHeights(Operation):
268 268
269 269 def run(self, dataOut, window):
270 270
271 271 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
272 272
273 273 if window == None:
274 274 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
275 275
276 276 newdelta = deltaHeight * window
277 277 r = dataOut.nHeights % window
278 278 newheights = (dataOut.nHeights-r)/window
279 279
280 280 if newheights <= 1:
281 281 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
282 282
283 283 if dataOut.flagDataAsBlock:
284 284 """
285 285 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
286 286 """
287 287 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
288 288 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
289 289 buffer = numpy.sum(buffer,3)
290 290
291 291 else:
292 292 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
293 293 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
294 294 buffer = numpy.sum(buffer,2)
295 295
296 296 dataOut.data = buffer
297 297 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
298 298 dataOut.windowOfFilter = window
299 299
300 300 return dataOut
301 301
302 302
303 303 class setH0(Operation):
304 304
305 305 def run(self, dataOut, h0, deltaHeight = None):
306 306
307 307 if not deltaHeight:
308 308 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
309 309
310 310 nHeights = dataOut.nHeights
311 311
312 312 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
313 313
314 314 dataOut.heightList = newHeiRange
315 315
316 316 return dataOut
317 317
318 318
319 319 class deFlip(Operation):
320 320
321 321 def run(self, dataOut, channelList = []):
322 322
323 323 data = dataOut.data.copy()
324 324
325 325 if dataOut.flagDataAsBlock:
326 326 flip = self.flip
327 327 profileList = list(range(dataOut.nProfiles))
328 328
329 329 if not channelList:
330 330 for thisProfile in profileList:
331 331 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
332 332 flip *= -1.0
333 333 else:
334 334 for thisChannel in channelList:
335 335 if thisChannel not in dataOut.channelList:
336 336 continue
337 337
338 338 for thisProfile in profileList:
339 339 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
340 340 flip *= -1.0
341 341
342 342 self.flip = flip
343 343
344 344 else:
345 345 if not channelList:
346 346 data[:,:] = data[:,:]*self.flip
347 347 else:
348 348 for thisChannel in channelList:
349 349 if thisChannel not in dataOut.channelList:
350 350 continue
351 351
352 352 data[thisChannel,:] = data[thisChannel,:]*self.flip
353 353
354 354 self.flip *= -1.
355 355
356 356 dataOut.data = data
357 357
358 358 return dataOut
359 359
360 360
361 361 class setAttribute(Operation):
362 362 '''
363 363 Set an arbitrary attribute(s) to dataOut
364 364 '''
365 365
366 366 def __init__(self):
367 367
368 368 Operation.__init__(self)
369 369 self._ready = False
370 370
371 371 def run(self, dataOut, **kwargs):
372 372
373 373 for key, value in kwargs.items():
374 374 setattr(dataOut, key, value)
375 375
376 376 return dataOut
377 377
378 378
379 379 @MPDecorator
380 380 class printAttribute(Operation):
381 381 '''
382 382 Print an arbitrary attribute of dataOut
383 383 '''
384 384
385 385 def __init__(self):
386 386
387 387 Operation.__init__(self)
388 388
389 389 def run(self, dataOut, attributes):
390 390
391 391 if isinstance(attributes, str):
392 392 attributes = [attributes]
393 393 for attr in attributes:
394 394 if hasattr(dataOut, attr):
395 395 log.log(getattr(dataOut, attr), attr)
396 396
397 397
398 398 class interpolateHeights(Operation):
399 399
400 400 def run(self, dataOut, topLim, botLim):
401 401 #69 al 72 para julia
402 402 #82-84 para meteoros
403 403 if len(numpy.shape(dataOut.data))==2:
404 404 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
405 405 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
406 406 #dataOut.data[:,botLim:limSup+1] = sampInterp
407 407 dataOut.data[:,botLim:topLim+1] = sampInterp
408 408 else:
409 409 nHeights = dataOut.data.shape[2]
410 410 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
411 411 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
412 412 f = interpolate.interp1d(x, y, axis = 2)
413 413 xnew = numpy.arange(botLim,topLim+1)
414 414 ynew = f(xnew)
415 415 dataOut.data[:,:,botLim:topLim+1] = ynew
416 416
417 417 return dataOut
418 418
419 419
420 420 class CohInt(Operation):
421 421
422 422 isConfig = False
423 423 __profIndex = 0
424 424 __byTime = False
425 425 __initime = None
426 426 __lastdatatime = None
427 427 __integrationtime = None
428 428 __buffer = None
429 429 __bufferStride = []
430 430 __dataReady = False
431 431 __profIndexStride = 0
432 432 __dataToPutStride = False
433 433 n = None
434 434
435 435 def __init__(self, **kwargs):
436 436
437 437 Operation.__init__(self, **kwargs)
438 438
439 439 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
440 440 """
441 441 Set the parameters of the integration class.
442 442
443 443 Inputs:
444 444
445 445 n : Number of coherent integrations
446 446 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
447 447 overlapping :
448 448 """
449 449
450 450 self.__initime = None
451 451 self.__lastdatatime = 0
452 452 self.__buffer = None
453 453 self.__dataReady = False
454 454 self.byblock = byblock
455 455 self.stride = stride
456 456
457 457 if n == None and timeInterval == None:
458 458 raise ValueError("n or timeInterval should be specified ...")
459 459
460 460 if n != None:
461 461 self.n = n
462 462 self.__byTime = False
463 463 else:
464 464 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
465 465 self.n = 9999
466 466 self.__byTime = True
467 467
468 468 if overlapping:
469 469 self.__withOverlapping = True
470 470 self.__buffer = None
471 471 else:
472 472 self.__withOverlapping = False
473 473 self.__buffer = 0
474 474
475 475 self.__profIndex = 0
476 476
477 477 def putData(self, data):
478 478
479 479 """
480 480 Add a profile to the __buffer and increase in one the __profileIndex
481 481
482 482 """
483 483
484 484 if not self.__withOverlapping:
485 485 self.__buffer += data.copy()
486 486 self.__profIndex += 1
487 487 return
488 488
489 489 #Overlapping data
490 490 nChannels, nHeis = data.shape
491 491 data = numpy.reshape(data, (1, nChannels, nHeis))
492 492
493 493 #If the buffer is empty then it takes the data value
494 494 if self.__buffer is None:
495 495 self.__buffer = data
496 496 self.__profIndex += 1
497 497 return
498 498
499 499 #If the buffer length is lower than n then stakcing the data value
500 500 if self.__profIndex < self.n:
501 501 self.__buffer = numpy.vstack((self.__buffer, data))
502 502 self.__profIndex += 1
503 503 return
504 504
505 505 #If the buffer length is equal to n then replacing the last buffer value with the data value
506 506 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
507 507 self.__buffer[self.n-1] = data
508 508 self.__profIndex = self.n
509 509 return
510 510
511 511
512 512 def pushData(self):
513 513 """
514 514 Return the sum of the last profiles and the profiles used in the sum.
515 515
516 516 Affected:
517 517
518 518 self.__profileIndex
519 519
520 520 """
521 521
522 522 if not self.__withOverlapping:
523 523 data = self.__buffer
524 524 n = self.__profIndex
525 525
526 526 self.__buffer = 0
527 527 self.__profIndex = 0
528 528
529 529 return data, n
530 530
531 531 #Integration with Overlapping
532 532 data = numpy.sum(self.__buffer, axis=0)
533 533 # print data
534 534 # raise
535 535 n = self.__profIndex
536 536
537 537 return data, n
538 538
539 539 def byProfiles(self, data):
540 540
541 541 self.__dataReady = False
542 542 avgdata = None
543 543 # n = None
544 544 # print data
545 545 # raise
546 546 self.putData(data)
547 547
548 548 if self.__profIndex == self.n:
549 549 avgdata, n = self.pushData()
550 550 self.__dataReady = True
551 551
552 552 return avgdata
553 553
554 554 def byTime(self, data, datatime):
555 555
556 556 self.__dataReady = False
557 557 avgdata = None
558 558 n = None
559 559
560 560 self.putData(data)
561 561
562 562 if (datatime - self.__initime) >= self.__integrationtime:
563 563 avgdata, n = self.pushData()
564 564 self.n = n
565 565 self.__dataReady = True
566 566
567 567 return avgdata
568 568
569 569 def integrateByStride(self, data, datatime):
570 570 # print data
571 571 if self.__profIndex == 0:
572 572 self.__buffer = [[data.copy(), datatime]]
573 573 else:
574 574 self.__buffer.append([data.copy(),datatime])
575 575 self.__profIndex += 1
576 576 self.__dataReady = False
577 577
578 578 if self.__profIndex == self.n * self.stride :
579 579 self.__dataToPutStride = True
580 580 self.__profIndexStride = 0
581 581 self.__profIndex = 0
582 582 self.__bufferStride = []
583 583 for i in range(self.stride):
584 584 current = self.__buffer[i::self.stride]
585 585 data = numpy.sum([t[0] for t in current], axis=0)
586 586 avgdatatime = numpy.average([t[1] for t in current])
587 587 # print data
588 588 self.__bufferStride.append((data, avgdatatime))
589 589
590 590 if self.__dataToPutStride:
591 591 self.__dataReady = True
592 592 self.__profIndexStride += 1
593 593 if self.__profIndexStride == self.stride:
594 594 self.__dataToPutStride = False
595 595 # print self.__bufferStride[self.__profIndexStride - 1]
596 596 # raise
597 597 return self.__bufferStride[self.__profIndexStride - 1]
598 598
599 599
600 600 return None, None
601 601
602 602 def integrate(self, data, datatime=None):
603 603
604 604 if self.__initime == None:
605 605 self.__initime = datatime
606 606
607 607 if self.__byTime:
608 608 avgdata = self.byTime(data, datatime)
609 609 else:
610 610 avgdata = self.byProfiles(data)
611 611
612 612
613 613 self.__lastdatatime = datatime
614 614
615 615 if avgdata is None:
616 616 return None, None
617 617
618 618 avgdatatime = self.__initime
619 619
620 620 deltatime = datatime - self.__lastdatatime
621 621
622 622 if not self.__withOverlapping:
623 623 self.__initime = datatime
624 624 else:
625 625 self.__initime += deltatime
626 626
627 627 return avgdata, avgdatatime
628 628
629 629 def integrateByBlock(self, dataOut):
630 630
631 631 times = int(dataOut.data.shape[1]/self.n)
632 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
632 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex_)
633 633
634 634 id_min = 0
635 635 id_max = self.n
636 636
637 637 for i in range(times):
638 638 junk = dataOut.data[:,id_min:id_max,:]
639 639 avgdata[:,i,:] = junk.sum(axis=1)
640 640 id_min += self.n
641 641 id_max += self.n
642 642
643 643 timeInterval = dataOut.ippSeconds*self.n
644 644 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
645 645 self.__dataReady = True
646 646 return avgdata, avgdatatime
647 647
648 648 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
649 649
650 650 if not self.isConfig:
651 651 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
652 652 self.isConfig = True
653 653 if dataOut.flagDataAsBlock:
654 654 """
655 655 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
656 656 """
657 657 avgdata, avgdatatime = self.integrateByBlock(dataOut)
658 658 dataOut.nProfiles /= self.n
659 659 else:
660 660 if stride is None:
661 661 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
662 662 else:
663 663 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
664 664
665 665
666 666 # dataOut.timeInterval *= n
667 667 dataOut.flagNoData = True
668 668
669 669 if self.__dataReady:
670 670 dataOut.data = avgdata
671 671 if not dataOut.flagCohInt:
672 672 dataOut.nCohInt *= self.n
673 673 dataOut.flagCohInt = True
674 674 dataOut.utctime = avgdatatime
675 675 # print avgdata, avgdatatime
676 676 # raise
677 677 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
678 678 dataOut.flagNoData = False
679 679 return dataOut
680 680
681 681 class Decoder(Operation):
682 682
683 683 isConfig = False
684 684 __profIndex = 0
685 685
686 686 code = None
687 687
688 688 nCode = None
689 689 nBaud = None
690 690
691 691 def __init__(self, **kwargs):
692 692
693 693 Operation.__init__(self, **kwargs)
694 694
695 695 self.times = None
696 696 self.osamp = None
697 697 # self.__setValues = False
698 698 self.isConfig = False
699 699 self.setupReq = False
700 700 def setup(self, code, osamp, dataOut):
701 701
702 702 self.__profIndex = 0
703 703
704 704 self.code = code
705 705
706 706 self.nCode = len(code)
707 707 self.nBaud = len(code[0])
708 708
709 709 if (osamp != None) and (osamp >1):
710 710 self.osamp = osamp
711 711 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
712 712 self.nBaud = self.nBaud*self.osamp
713 713
714 714 self.__nChannels = dataOut.nChannels
715 715 self.__nProfiles = dataOut.nProfiles
716 716 self.__nHeis = dataOut.nHeights
717 717
718 718 if self.__nHeis < self.nBaud:
719 719 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
720 720
721 721 #Frequency
722 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
722 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex_)
723 723
724 724 __codeBuffer[:,0:self.nBaud] = self.code
725 725
726 726 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
727 727
728 728 if dataOut.flagDataAsBlock:
729 729
730 730 self.ndatadec = self.__nHeis #- self.nBaud + 1
731 731
732 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
732 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex_)
733 733
734 734 else:
735 735
736 736 #Time
737 737 self.ndatadec = self.__nHeis #- self.nBaud + 1
738 738
739 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
739 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex_)
740 740
741 741 def __convolutionInFreq(self, data):
742 742
743 743 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
744 744
745 745 fft_data = numpy.fft.fft(data, axis=1)
746 746
747 747 conv = fft_data*fft_code
748 748
749 749 data = numpy.fft.ifft(conv,axis=1)
750 750
751 751 return data
752 752
753 753 def __convolutionInFreqOpt(self, data):
754 754
755 755 raise NotImplementedError
756 756
757 757 def __convolutionInTime(self, data):
758 758
759 759 code = self.code[self.__profIndex]
760 760 for i in range(self.__nChannels):
761 761 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
762 762
763 763 return self.datadecTime
764 764
765 765 def __convolutionByBlockInTime(self, data):
766 766
767 767 repetitions = int(self.__nProfiles / self.nCode)
768 768 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
769 769 junk = junk.flatten()
770 770 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
771 771 profilesList = range(self.__nProfiles)
772 772
773 773 for i in range(self.__nChannels):
774 774 for j in profilesList:
775 775 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
776 776 return self.datadecTime
777 777
778 778 def __convolutionByBlockInFreq(self, data):
779 779
780 780 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
781 781
782 782
783 783 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
784 784
785 785 fft_data = numpy.fft.fft(data, axis=2)
786 786
787 787 conv = fft_data*fft_code
788 788
789 789 data = numpy.fft.ifft(conv,axis=2)
790 790
791 791 return data
792 792
793 793
794 794 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
795 795
796 796 if dataOut.flagDecodeData:
797 797 print("This data is already decoded, recoding again ...")
798 798
799 799 if not self.isConfig:
800 800
801 801 if code is None:
802 802 if dataOut.code is None:
803 803 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
804 804
805 805 code = dataOut.code
806 806 else:
807 807 code = numpy.array(code).reshape(nCode,nBaud)
808 808 self.setup(code, osamp, dataOut)
809 809
810 810 self.isConfig = True
811 811
812 812 if mode == 3:
813 813 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
814 814
815 815 if times != None:
816 816 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
817 817
818 818 if self.code is None:
819 819 print("Fail decoding: Code is not defined.")
820 820 return
821 821
822 822 self.__nProfiles = dataOut.nProfiles
823 823 datadec = None
824 824
825 825 if mode == 3:
826 826 mode = 0
827 827
828 828 if dataOut.flagDataAsBlock:
829 829 """
830 830 Decoding when data have been read as block,
831 831 """
832 832
833 833 if mode == 0:
834 834 datadec = self.__convolutionByBlockInTime(dataOut.data)
835 835 if mode == 1:
836 836 datadec = self.__convolutionByBlockInFreq(dataOut.data)
837 837 else:
838 838 """
839 839 Decoding when data have been read profile by profile
840 840 """
841 841 if mode == 0:
842 842 datadec = self.__convolutionInTime(dataOut.data)
843 843
844 844 if mode == 1:
845 845 datadec = self.__convolutionInFreq(dataOut.data)
846 846
847 847 if mode == 2:
848 848 datadec = self.__convolutionInFreqOpt(dataOut.data)
849 849
850 850 if datadec is None:
851 851 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
852 852
853 853 dataOut.code = self.code
854 854 dataOut.nCode = self.nCode
855 855 dataOut.nBaud = self.nBaud
856 856
857 857 dataOut.data = datadec
858 858 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
859 859
860 860 dataOut.flagDecodeData = True #asumo q la data esta decodificada
861 861
862 862 if self.__profIndex == self.nCode-1:
863 863 self.__profIndex = 0
864 864 return dataOut
865 865
866 866 self.__profIndex += 1
867 867
868 868 return dataOut
869 869
870 870
871 871 class ProfileConcat(Operation):
872 872
873 873 isConfig = False
874 874 buffer = None
875 875
876 876 def __init__(self, **kwargs):
877 877
878 878 Operation.__init__(self, **kwargs)
879 879 self.profileIndex = 0
880 880
881 881 def reset(self):
882 882 self.buffer = numpy.zeros_like(self.buffer)
883 883 self.start_index = 0
884 884 self.times = 1
885 885
886 886 def setup(self, data, m, n=1):
887 887 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
888 888 self.nHeights = data.shape[1]#.nHeights
889 889 self.start_index = 0
890 890 self.times = 1
891 891
892 892 def concat(self, data):
893 893
894 894 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
895 895 self.start_index = self.start_index + self.nHeights
896 896
897 897 def run(self, dataOut, m):
898 898 dataOut.flagNoData = True
899 899
900 900 if not self.isConfig:
901 901 self.setup(dataOut.data, m, 1)
902 902 self.isConfig = True
903 903
904 904 if dataOut.flagDataAsBlock:
905 905 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
906 906
907 907 else:
908 908 self.concat(dataOut.data)
909 909 self.times += 1
910 910 if self.times > m:
911 911 dataOut.data = self.buffer
912 912 self.reset()
913 913 dataOut.flagNoData = False
914 914 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
915 915 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
916 916 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
917 917 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
918 918 dataOut.ippSeconds *= m
919 919 return dataOut
920 920
921 921 class ProfileSelector(Operation):
922 922
923 923 profileIndex = None
924 924 # Tamanho total de los perfiles
925 925 nProfiles = None
926 926
927 927 def __init__(self, **kwargs):
928 928
929 929 Operation.__init__(self, **kwargs)
930 930 self.profileIndex = 0
931 931
932 932 def incProfileIndex(self):
933 933
934 934 self.profileIndex += 1
935 935
936 936 if self.profileIndex >= self.nProfiles:
937 937 self.profileIndex = 0
938 938
939 939 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
940 940
941 941 if profileIndex < minIndex:
942 942 return False
943 943
944 944 if profileIndex > maxIndex:
945 945 return False
946 946
947 947 return True
948 948
949 949 def isThisProfileInList(self, profileIndex, profileList):
950 950
951 951 if profileIndex not in profileList:
952 952 return False
953 953
954 954 return True
955 955
956 956 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
957 957
958 958 """
959 959 ProfileSelector:
960 960
961 961 Inputs:
962 962 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
963 963
964 964 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
965 965
966 966 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
967 967
968 968 """
969 969
970 970 if rangeList is not None:
971 971 if type(rangeList[0]) not in (tuple, list):
972 972 rangeList = [rangeList]
973 973
974 974 dataOut.flagNoData = True
975 975
976 976 if dataOut.flagDataAsBlock:
977 977 """
978 978 data dimension = [nChannels, nProfiles, nHeis]
979 979 """
980 980 if profileList != None:
981 981 dataOut.data = dataOut.data[:,profileList,:]
982 982
983 983 if profileRangeList != None:
984 984 minIndex = profileRangeList[0]
985 985 maxIndex = profileRangeList[1]
986 986 profileList = list(range(minIndex, maxIndex+1))
987 987
988 988 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
989 989
990 990 if rangeList != None:
991 991
992 992 profileList = []
993 993
994 994 for thisRange in rangeList:
995 995 minIndex = thisRange[0]
996 996 maxIndex = thisRange[1]
997 997
998 998 profileList.extend(list(range(minIndex, maxIndex+1)))
999 999
1000 1000 dataOut.data = dataOut.data[:,profileList,:]
1001 1001
1002 1002 dataOut.nProfiles = len(profileList)
1003 1003 dataOut.profileIndex = dataOut.nProfiles - 1
1004 1004 dataOut.flagNoData = False
1005 1005
1006 1006 return dataOut
1007 1007
1008 1008 """
1009 1009 data dimension = [nChannels, nHeis]
1010 1010 """
1011 1011
1012 1012 if profileList != None:
1013 1013
1014 1014 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1015 1015
1016 1016 self.nProfiles = len(profileList)
1017 1017 dataOut.nProfiles = self.nProfiles
1018 1018 dataOut.profileIndex = self.profileIndex
1019 1019 dataOut.flagNoData = False
1020 1020
1021 1021 self.incProfileIndex()
1022 1022 return dataOut
1023 1023
1024 1024 if profileRangeList != None:
1025 1025
1026 1026 minIndex = profileRangeList[0]
1027 1027 maxIndex = profileRangeList[1]
1028 1028
1029 1029 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1030 1030
1031 1031 self.nProfiles = maxIndex - minIndex + 1
1032 1032 dataOut.nProfiles = self.nProfiles
1033 1033 dataOut.profileIndex = self.profileIndex
1034 1034 dataOut.flagNoData = False
1035 1035
1036 1036 self.incProfileIndex()
1037 1037 return dataOut
1038 1038
1039 1039 if rangeList != None:
1040 1040
1041 1041 nProfiles = 0
1042 1042
1043 1043 for thisRange in rangeList:
1044 1044 minIndex = thisRange[0]
1045 1045 maxIndex = thisRange[1]
1046 1046
1047 1047 nProfiles += maxIndex - minIndex + 1
1048 1048
1049 1049 for thisRange in rangeList:
1050 1050
1051 1051 minIndex = thisRange[0]
1052 1052 maxIndex = thisRange[1]
1053 1053
1054 1054 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1055 1055
1056 1056 self.nProfiles = nProfiles
1057 1057 dataOut.nProfiles = self.nProfiles
1058 1058 dataOut.profileIndex = self.profileIndex
1059 1059 dataOut.flagNoData = False
1060 1060
1061 1061 self.incProfileIndex()
1062 1062
1063 1063 break
1064 1064
1065 1065 return dataOut
1066 1066
1067 1067
1068 1068 if beam != None: #beam is only for AMISR data
1069 1069 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1070 1070 dataOut.flagNoData = False
1071 1071 dataOut.profileIndex = self.profileIndex
1072 1072
1073 1073 self.incProfileIndex()
1074 1074
1075 1075 return dataOut
1076 1076
1077 1077 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1078 1078
1079 1079
1080 1080 class Reshaper(Operation):
1081 1081
1082 1082 def __init__(self, **kwargs):
1083 1083
1084 1084 Operation.__init__(self, **kwargs)
1085 1085
1086 1086 self.__buffer = None
1087 1087 self.__nitems = 0
1088 1088
1089 1089 def __appendProfile(self, dataOut, nTxs):
1090 1090
1091 1091 if self.__buffer is None:
1092 1092 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1093 1093 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1094 1094
1095 1095 ini = dataOut.nHeights * self.__nitems
1096 1096 end = ini + dataOut.nHeights
1097 1097
1098 1098 self.__buffer[:, ini:end] = dataOut.data
1099 1099
1100 1100 self.__nitems += 1
1101 1101
1102 1102 return int(self.__nitems*nTxs)
1103 1103
1104 1104 def __getBuffer(self):
1105 1105
1106 1106 if self.__nitems == int(1./self.__nTxs):
1107 1107
1108 1108 self.__nitems = 0
1109 1109
1110 1110 return self.__buffer.copy()
1111 1111
1112 1112 return None
1113 1113
1114 1114 def __checkInputs(self, dataOut, shape, nTxs):
1115 1115
1116 1116 if shape is None and nTxs is None:
1117 1117 raise ValueError("Reshaper: shape of factor should be defined")
1118 1118
1119 1119 if nTxs:
1120 1120 if nTxs < 0:
1121 1121 raise ValueError("nTxs should be greater than 0")
1122 1122
1123 1123 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1124 1124 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1125 1125
1126 1126 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1127 1127
1128 1128 return shape, nTxs
1129 1129
1130 1130 if len(shape) != 2 and len(shape) != 3:
1131 1131 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1132 1132
1133 1133 if len(shape) == 2:
1134 1134 shape_tuple = [dataOut.nChannels]
1135 1135 shape_tuple.extend(shape)
1136 1136 else:
1137 1137 shape_tuple = list(shape)
1138 1138
1139 1139 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1140 1140
1141 1141 return shape_tuple, nTxs
1142 1142
1143 1143 def run(self, dataOut, shape=None, nTxs=None):
1144 1144
1145 1145 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1146 1146
1147 1147 dataOut.flagNoData = True
1148 1148 profileIndex = None
1149 1149
1150 1150 if dataOut.flagDataAsBlock:
1151 1151
1152 1152 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1153 1153 dataOut.flagNoData = False
1154 1154
1155 1155 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1156 1156
1157 1157 else:
1158 1158
1159 1159 if self.__nTxs < 1:
1160 1160
1161 1161 self.__appendProfile(dataOut, self.__nTxs)
1162 1162 new_data = self.__getBuffer()
1163 1163
1164 1164 if new_data is not None:
1165 1165 dataOut.data = new_data
1166 1166 dataOut.flagNoData = False
1167 1167
1168 1168 profileIndex = dataOut.profileIndex*nTxs
1169 1169
1170 1170 else:
1171 1171 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1172 1172
1173 1173 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1174 1174
1175 1175 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1176 1176
1177 1177 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1178 1178
1179 1179 dataOut.profileIndex = profileIndex
1180 1180
1181 1181 dataOut.ippSeconds /= self.__nTxs
1182 1182
1183 1183 return dataOut
1184 1184
1185 1185 class SplitProfiles(Operation):
1186 1186
1187 1187 def __init__(self, **kwargs):
1188 1188
1189 1189 Operation.__init__(self, **kwargs)
1190 1190
1191 1191 def run(self, dataOut, n):
1192 1192
1193 1193 dataOut.flagNoData = True
1194 1194 profileIndex = None
1195 1195
1196 1196 if dataOut.flagDataAsBlock:
1197 1197
1198 1198 #nchannels, nprofiles, nsamples
1199 1199 shape = dataOut.data.shape
1200 1200
1201 1201 if shape[2] % n != 0:
1202 1202 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1203 1203
1204 1204 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1205 1205
1206 1206 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1207 1207 dataOut.flagNoData = False
1208 1208
1209 1209 profileIndex = int(dataOut.nProfiles/n) - 1
1210 1210
1211 1211 else:
1212 1212
1213 1213 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1214 1214
1215 1215 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1216 1216
1217 1217 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1218 1218
1219 1219 dataOut.nProfiles = int(dataOut.nProfiles*n)
1220 1220
1221 1221 dataOut.profileIndex = profileIndex
1222 1222
1223 1223 dataOut.ippSeconds /= n
1224 1224
1225 1225 return dataOut
1226 1226
1227 1227 class CombineProfiles(Operation):
1228 1228 def __init__(self, **kwargs):
1229 1229
1230 1230 Operation.__init__(self, **kwargs)
1231 1231
1232 1232 self.__remData = None
1233 1233 self.__profileIndex = 0
1234 1234
1235 1235 def run(self, dataOut, n):
1236 1236
1237 1237 dataOut.flagNoData = True
1238 1238 profileIndex = None
1239 1239
1240 1240 if dataOut.flagDataAsBlock:
1241 1241
1242 1242 #nchannels, nprofiles, nsamples
1243 1243 shape = dataOut.data.shape
1244 1244 new_shape = shape[0], shape[1]/n, shape[2]*n
1245 1245
1246 1246 if shape[1] % n != 0:
1247 1247 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1248 1248
1249 1249 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1250 1250 dataOut.flagNoData = False
1251 1251
1252 1252 profileIndex = int(dataOut.nProfiles*n) - 1
1253 1253
1254 1254 else:
1255 1255
1256 1256 #nchannels, nsamples
1257 1257 if self.__remData is None:
1258 1258 newData = dataOut.data
1259 1259 else:
1260 1260 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1261 1261
1262 1262 self.__profileIndex += 1
1263 1263
1264 1264 if self.__profileIndex < n:
1265 1265 self.__remData = newData
1266 1266 #continue
1267 1267 return
1268 1268
1269 1269 self.__profileIndex = 0
1270 1270 self.__remData = None
1271 1271
1272 1272 dataOut.data = newData
1273 1273 dataOut.flagNoData = False
1274 1274
1275 1275 profileIndex = dataOut.profileIndex/n
1276 1276
1277 1277
1278 1278 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1279 1279
1280 1280 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1281 1281
1282 1282 dataOut.nProfiles = int(dataOut.nProfiles/n)
1283 1283
1284 1284 dataOut.profileIndex = profileIndex
1285 1285
1286 1286 dataOut.ippSeconds *= n
1287 1287
1288 1288 return dataOut
1289 1289
1290 1290 class PulsePairVoltage(Operation):
1291 1291 '''
1292 1292 Function PulsePair(Signal Power, Velocity)
1293 1293 The real component of Lag[0] provides Intensity Information
1294 1294 The imag component of Lag[1] Phase provides Velocity Information
1295 1295
1296 1296 Configuration Parameters:
1297 1297 nPRF = Number of Several PRF
1298 1298 theta = Degree Azimuth angel Boundaries
1299 1299
1300 1300 Input:
1301 1301 self.dataOut
1302 1302 lag[N]
1303 1303 Affected:
1304 1304 self.dataOut.spc
1305 1305 '''
1306 1306 isConfig = False
1307 1307 __profIndex = 0
1308 1308 __initime = None
1309 1309 __lastdatatime = None
1310 1310 __buffer = None
1311 1311 noise = None
1312 1312 __dataReady = False
1313 1313 n = None
1314 1314 __nch = 0
1315 1315 __nHeis = 0
1316 1316 removeDC = False
1317 1317 ipp = None
1318 1318 lambda_ = 0
1319 1319
1320 1320 def __init__(self,**kwargs):
1321 1321 Operation.__init__(self,**kwargs)
1322 1322
1323 1323 def setup(self, dataOut, n = None, removeDC=False):
1324 1324 '''
1325 1325 n= Numero de PRF's de entrada
1326 1326 '''
1327 1327 self.__initime = None
1328 1328 self.__lastdatatime = 0
1329 1329 self.__dataReady = False
1330 1330 self.__buffer = 0
1331 1331 self.__profIndex = 0
1332 1332 self.noise = None
1333 1333 self.__nch = dataOut.nChannels
1334 1334 self.__nHeis = dataOut.nHeights
1335 1335 self.removeDC = removeDC
1336 1336 self.lambda_ = 3.0e8/(9345.0e6)
1337 1337 self.ippSec = dataOut.ippSeconds
1338 1338 self.nCohInt = dataOut.nCohInt
1339 1339 print("IPPseconds",dataOut.ippSeconds)
1340 1340
1341 1341 print("ELVALOR DE n es:", n)
1342 1342 if n == None:
1343 1343 raise ValueError("n should be specified.")
1344 1344
1345 1345 if n != None:
1346 1346 if n<2:
1347 1347 raise ValueError("n should be greater than 2")
1348 1348
1349 1349 self.n = n
1350 1350 self.__nProf = n
1351 1351
1352 1352 self.__buffer = numpy.zeros((dataOut.nChannels,
1353 1353 n,
1354 1354 dataOut.nHeights),
1355 1355 dtype='complex')
1356 1356
1357 1357 def putData(self,data):
1358 1358 '''
1359 1359 Add a profile to he __buffer and increase in one the __profiel Index
1360 1360 '''
1361 1361 self.__buffer[:,self.__profIndex,:]= data
1362 1362 self.__profIndex += 1
1363 1363 return
1364 1364
1365 1365 def pushData(self,dataOut):
1366 1366 '''
1367 1367 Return the PULSEPAIR and the profiles used in the operation
1368 1368 Affected : self.__profileIndex
1369 1369 '''
1370 1370 #----------------- Remove DC-----------------------------------
1371 1371 if self.removeDC==True:
1372 1372 mean = numpy.mean(self.__buffer,1)
1373 1373 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1374 1374 dc= numpy.tile(tmp,[1,self.__nProf,1])
1375 1375 self.__buffer = self.__buffer - dc
1376 1376 #------------------Calculo de Potencia ------------------------
1377 1377 pair0 = self.__buffer*numpy.conj(self.__buffer)
1378 1378 pair0 = pair0.real
1379 1379 lag_0 = numpy.sum(pair0,1)
1380 1380 #------------------Calculo de Ruido x canal--------------------
1381 1381 self.noise = numpy.zeros(self.__nch)
1382 1382 for i in range(self.__nch):
1383 1383 daux = numpy.sort(pair0[i,:,:],axis= None)
1384 1384 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1385 1385
1386 1386 self.noise = self.noise.reshape(self.__nch,1)
1387 1387 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1388 1388 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1389 1389 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1390 1390 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1391 1391 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1392 1392 #-------------------- Power --------------------------------------------------
1393 1393 data_power = lag_0/(self.n*self.nCohInt)
1394 1394 #------------------ Senal ---------------------------------------------------
1395 1395 data_intensity = pair0 - noise_buffer
1396 1396 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1397 1397 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1398 1398 for i in range(self.__nch):
1399 1399 for j in range(self.__nHeis):
1400 1400 if data_intensity[i][j] < 0:
1401 1401 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1402 1402
1403 1403 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1404 1404 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1405 1405 lag_1 = numpy.sum(pair1,1)
1406 1406 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1407 1407 data_velocity = (self.lambda_/2.0)*data_freq
1408 1408
1409 1409 #---------------- Potencia promedio estimada de la Senal-----------
1410 1410 lag_0 = lag_0/self.n
1411 1411 S = lag_0-self.noise
1412 1412
1413 1413 #---------------- Frecuencia Doppler promedio ---------------------
1414 1414 lag_1 = lag_1/(self.n-1)
1415 1415 R1 = numpy.abs(lag_1)
1416 1416
1417 1417 #---------------- Calculo del SNR----------------------------------
1418 1418 data_snrPP = S/self.noise
1419 1419 for i in range(self.__nch):
1420 1420 for j in range(self.__nHeis):
1421 1421 if data_snrPP[i][j] < 1.e-20:
1422 1422 data_snrPP[i][j] = 1.e-20
1423 1423
1424 1424 #----------------- Calculo del ancho espectral ----------------------
1425 1425 L = S/R1
1426 1426 L = numpy.where(L<0,1,L)
1427 1427 L = numpy.log(L)
1428 1428 tmp = numpy.sqrt(numpy.absolute(L))
1429 1429 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1430 1430 n = self.__profIndex
1431 1431
1432 1432 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1433 1433 self.__profIndex = 0
1434 1434 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
1435 1435
1436 1436
1437 1437 def pulsePairbyProfiles(self,dataOut):
1438 1438
1439 1439 self.__dataReady = False
1440 1440 data_power = None
1441 1441 data_intensity = None
1442 1442 data_velocity = None
1443 1443 data_specwidth = None
1444 1444 data_snrPP = None
1445 1445 self.putData(data=dataOut.data)
1446 1446 if self.__profIndex == self.n:
1447 1447 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1448 1448 self.__dataReady = True
1449 1449
1450 1450 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
1451 1451
1452 1452
1453 1453 def pulsePairOp(self, dataOut, datatime= None):
1454 1454
1455 1455 if self.__initime == None:
1456 1456 self.__initime = datatime
1457 1457 data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
1458 1458 self.__lastdatatime = datatime
1459 1459
1460 1460 if data_power is None:
1461 1461 return None, None, None,None,None,None
1462 1462
1463 1463 avgdatatime = self.__initime
1464 1464 deltatime = datatime - self.__lastdatatime
1465 1465 self.__initime = datatime
1466 1466
1467 1467 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
1468 1468
1469 1469 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1470 1470
1471 1471 if not self.isConfig:
1472 1472 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1473 1473 self.isConfig = True
1474 1474 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1475 1475 dataOut.flagNoData = True
1476 1476
1477 1477 if self.__dataReady:
1478 1478 dataOut.nCohInt *= self.n
1479 1479 dataOut.dataPP_POW = data_intensity # S
1480 1480 dataOut.dataPP_POWER = data_power # P
1481 1481 dataOut.dataPP_DOP = data_velocity
1482 1482 dataOut.dataPP_SNR = data_snrPP
1483 1483 dataOut.dataPP_WIDTH = data_specwidth
1484 1484 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1485 1485 dataOut.utctime = avgdatatime
1486 1486 dataOut.flagNoData = False
1487 1487 return dataOut
1488 1488
1489 1489
1490 1490
1491 1491 # import collections
1492 1492 # from scipy.stats import mode
1493 1493 #
1494 1494 # class Synchronize(Operation):
1495 1495 #
1496 1496 # isConfig = False
1497 1497 # __profIndex = 0
1498 1498 #
1499 1499 # def __init__(self, **kwargs):
1500 1500 #
1501 1501 # Operation.__init__(self, **kwargs)
1502 1502 # # self.isConfig = False
1503 1503 # self.__powBuffer = None
1504 1504 # self.__startIndex = 0
1505 1505 # self.__pulseFound = False
1506 1506 #
1507 1507 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1508 1508 #
1509 1509 # #Read data
1510 1510 #
1511 1511 # powerdB = dataOut.getPower(channel = channel)
1512 1512 # noisedB = dataOut.getNoise(channel = channel)[0]
1513 1513 #
1514 1514 # self.__powBuffer.extend(powerdB.flatten())
1515 1515 #
1516 1516 # dataArray = numpy.array(self.__powBuffer)
1517 1517 #
1518 1518 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1519 1519 #
1520 1520 # maxValue = numpy.nanmax(filteredPower)
1521 1521 #
1522 1522 # if maxValue < noisedB + 10:
1523 1523 # #No se encuentra ningun pulso de transmision
1524 1524 # return None
1525 1525 #
1526 1526 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1527 1527 #
1528 1528 # if len(maxValuesIndex) < 2:
1529 1529 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1530 1530 # return None
1531 1531 #
1532 1532 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1533 1533 #
1534 1534 # #Seleccionar solo valores con un espaciamiento de nSamples
1535 1535 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1536 1536 #
1537 1537 # if len(pulseIndex) < 2:
1538 1538 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1539 1539 # return None
1540 1540 #
1541 1541 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1542 1542 #
1543 1543 # #remover senales que se distancien menos de 10 unidades o muestras
1544 1544 # #(No deberian existir IPP menor a 10 unidades)
1545 1545 #
1546 1546 # realIndex = numpy.where(spacing > 10 )[0]
1547 1547 #
1548 1548 # if len(realIndex) < 2:
1549 1549 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1550 1550 # return None
1551 1551 #
1552 1552 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1553 1553 # realPulseIndex = pulseIndex[realIndex]
1554 1554 #
1555 1555 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1556 1556 #
1557 1557 # print "IPP = %d samples" %period
1558 1558 #
1559 1559 # self.__newNSamples = dataOut.nHeights #int(period)
1560 1560 # self.__startIndex = int(realPulseIndex[0])
1561 1561 #
1562 1562 # return 1
1563 1563 #
1564 1564 #
1565 1565 # def setup(self, nSamples, nChannels, buffer_size = 4):
1566 1566 #
1567 1567 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1568 1568 # maxlen = buffer_size*nSamples)
1569 1569 #
1570 1570 # bufferList = []
1571 1571 #
1572 1572 # for i in range(nChannels):
1573 1573 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1574 1574 # maxlen = buffer_size*nSamples)
1575 1575 #
1576 1576 # bufferList.append(bufferByChannel)
1577 1577 #
1578 1578 # self.__nSamples = nSamples
1579 1579 # self.__nChannels = nChannels
1580 1580 # self.__bufferList = bufferList
1581 1581 #
1582 1582 # def run(self, dataOut, channel = 0):
1583 1583 #
1584 1584 # if not self.isConfig:
1585 1585 # nSamples = dataOut.nHeights
1586 1586 # nChannels = dataOut.nChannels
1587 1587 # self.setup(nSamples, nChannels)
1588 1588 # self.isConfig = True
1589 1589 #
1590 1590 # #Append new data to internal buffer
1591 1591 # for thisChannel in range(self.__nChannels):
1592 1592 # bufferByChannel = self.__bufferList[thisChannel]
1593 1593 # bufferByChannel.extend(dataOut.data[thisChannel])
1594 1594 #
1595 1595 # if self.__pulseFound:
1596 1596 # self.__startIndex -= self.__nSamples
1597 1597 #
1598 1598 # #Finding Tx Pulse
1599 1599 # if not self.__pulseFound:
1600 1600 # indexFound = self.__findTxPulse(dataOut, channel)
1601 1601 #
1602 1602 # if indexFound == None:
1603 1603 # dataOut.flagNoData = True
1604 1604 # return
1605 1605 #
1606 1606 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1607 1607 # self.__pulseFound = True
1608 1608 # self.__startIndex = indexFound
1609 1609 #
1610 1610 # #If pulse was found ...
1611 1611 # for thisChannel in range(self.__nChannels):
1612 1612 # bufferByChannel = self.__bufferList[thisChannel]
1613 1613 # #print self.__startIndex
1614 1614 # x = numpy.array(bufferByChannel)
1615 1615 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1616 1616 #
1617 1617 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1618 1618 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1619 1619 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1620 1620 #
1621 1621 # dataOut.data = self.__arrayBuffer
1622 1622 #
1623 1623 # self.__startIndex += self.__newNSamples
1624 1624 #
1625 1625 # return No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now