##// END OF EJS Templates
schain-jc
RSS
remove profilehooks imports
Juan C. Espinoza -
r1050:d6a89c93744f
parent child
Show More
Side by Side Unified
Context file:
r1050:d6a89c93744f Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1306 +1,1283
1 1 import sys
2 2 import numpy
3 from profilehooks import profile
4 3 from scipy import interpolate
5 4 from schainpy import cSchain
6 5 from jroproc_base import ProcessingUnit, Operation
7 6 from schainpy.model.data.jrodata import Voltage
8 7 from time import time
9 8
10 9 class VoltageProc(ProcessingUnit):
11 10
12 11
13 12 def __init__(self, **kwargs):
14 13
15 14 ProcessingUnit.__init__(self, **kwargs)
16 15
17 16 # self.objectDict = {}
18 17 self.dataOut = Voltage()
19 18 self.flip = 1
20 19
21 20 def run(self):
22 21 if self.dataIn.type == 'AMISR':
23 22 self.__updateObjFromAmisrInput()
24 23
25 24 if self.dataIn.type == 'Voltage':
26 25 self.dataOut.copy(self.dataIn)
27 26
28 27 # self.dataOut.copy(self.dataIn)
29 28
30 29 def __updateObjFromAmisrInput(self):
31 30
32 31 self.dataOut.timeZone = self.dataIn.timeZone
33 32 self.dataOut.dstFlag = self.dataIn.dstFlag
34 33 self.dataOut.errorCount = self.dataIn.errorCount
35 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
36 35
37 36 self.dataOut.flagNoData = self.dataIn.flagNoData
38 37 self.dataOut.data = self.dataIn.data
39 38 self.dataOut.utctime = self.dataIn.utctime
40 39 self.dataOut.channelList = self.dataIn.channelList
41 40 # self.dataOut.timeInterval = self.dataIn.timeInterval
42 41 self.dataOut.heightList = self.dataIn.heightList
43 42 self.dataOut.nProfiles = self.dataIn.nProfiles
44 43
45 44 self.dataOut.nCohInt = self.dataIn.nCohInt
46 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
47 46 self.dataOut.frequency = self.dataIn.frequency
48 47
49 48 self.dataOut.azimuth = self.dataIn.azimuth
50 49 self.dataOut.zenith = self.dataIn.zenith
51 50
52 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
53 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
54 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
55 54 #
56 55 # pass#
57 56 #
58 57 # def init(self):
59 58 #
60 59 #
61 60 # if self.dataIn.type == 'AMISR':
62 61 # self.__updateObjFromAmisrInput()
63 62 #
64 63 # if self.dataIn.type == 'Voltage':
65 64 # self.dataOut.copy(self.dataIn)
66 65 # # No necesita copiar en cada init() los atributos de dataIn
67 66 # # la copia deberia hacerse por cada nuevo bloque de datos
68 67
69 68 def selectChannels(self, channelList):
70 69
71 70 channelIndexList = []
72 71
73 72 for channel in channelList:
74 73 if channel not in self.dataOut.channelList:
75 74 raise ValueError, "Channel %d is not in %s" %(channel, str(self.dataOut.channelList))
76 75
77 76 index = self.dataOut.channelList.index(channel)
78 77 channelIndexList.append(index)
79 78
80 79 self.selectChannelsByIndex(channelIndexList)
81 80
82 81 def selectChannelsByIndex(self, channelIndexList):
83 82 """
84 83 Selecciona un bloque de datos en base a canales segun el channelIndexList
85 84
86 85 Input:
87 86 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
88 87
89 88 Affected:
90 89 self.dataOut.data
91 90 self.dataOut.channelIndexList
92 91 self.dataOut.nChannels
93 92 self.dataOut.m_ProcessingHeader.totalSpectra
94 93 self.dataOut.systemHeaderObj.numChannels
95 94 self.dataOut.m_ProcessingHeader.blockSize
96 95
97 96 Return:
98 97 None
99 98 """
100 99
101 100 for channelIndex in channelIndexList:
102 101 if channelIndex not in self.dataOut.channelIndexList:
103 102 print channelIndexList
104 103 raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
105 104
106 105 if self.dataOut.flagDataAsBlock:
107 106 """
108 107 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
109 108 """
110 109 data = self.dataOut.data[channelIndexList,:,:]
111 110 else:
112 111 data = self.dataOut.data[channelIndexList,:]
113 112
114 113 self.dataOut.data = data
115 114 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
116 115 # self.dataOut.nChannels = nChannels
117 116
118 117 return 1
119 118
120 119 def selectHeights(self, minHei=None, maxHei=None):
121 120 """
122 121 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
123 122 minHei <= height <= maxHei
124 123
125 124 Input:
126 125 minHei : valor minimo de altura a considerar
127 126 maxHei : valor maximo de altura a considerar
128 127
129 128 Affected:
130 129 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
131 130
132 131 Return:
133 132 1 si el metodo se ejecuto con exito caso contrario devuelve 0
134 133 """
135 134
136 135 if minHei == None:
137 136 minHei = self.dataOut.heightList[0]
138 137
139 138 if maxHei == None:
140 139 maxHei = self.dataOut.heightList[-1]
141 140
142 141 if (minHei < self.dataOut.heightList[0]):
143 142 minHei = self.dataOut.heightList[0]
144 143
145 144 if (maxHei > self.dataOut.heightList[-1]):
146 145 maxHei = self.dataOut.heightList[-1]
147 146
148 147 minIndex = 0
149 148 maxIndex = 0
150 149 heights = self.dataOut.heightList
151 150
152 151 inda = numpy.where(heights >= minHei)
153 152 indb = numpy.where(heights <= maxHei)
154 153
155 154 try:
156 155 minIndex = inda[0][0]
157 156 except:
158 157 minIndex = 0
159 158
160 159 try:
161 160 maxIndex = indb[0][-1]
162 161 except:
163 162 maxIndex = len(heights)
164 163
165 164 self.selectHeightsByIndex(minIndex, maxIndex)
166 165
167 166 return 1
168 167
169 168
170 169 def selectHeightsByIndex(self, minIndex, maxIndex):
171 170 """
172 171 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
173 172 minIndex <= index <= maxIndex
174 173
175 174 Input:
176 175 minIndex : valor de indice minimo de altura a considerar
177 176 maxIndex : valor de indice maximo de altura a considerar
178 177
179 178 Affected:
180 179 self.dataOut.data
181 180 self.dataOut.heightList
182 181
183 182 Return:
184 183 1 si el metodo se ejecuto con exito caso contrario devuelve 0
185 184 """
186 185
187 186 if (minIndex < 0) or (minIndex > maxIndex):
188 187 raise ValueError, "Height index range (%d,%d) is not valid" % (minIndex, maxIndex)
189 188
190 189 if (maxIndex >= self.dataOut.nHeights):
191 190 maxIndex = self.dataOut.nHeights
192 191
193 192 #voltage
194 193 if self.dataOut.flagDataAsBlock:
195 194 """
196 195 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
197 196 """
198 197 data = self.dataOut.data[:,:, minIndex:maxIndex]
199 198 else:
200 199 data = self.dataOut.data[:, minIndex:maxIndex]
201 200
202 201 # firstHeight = self.dataOut.heightList[minIndex]
203 202
204 203 self.dataOut.data = data
205 204 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
206 205
207 206 if self.dataOut.nHeights <= 1:
208 207 raise ValueError, "selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights)
209 208
210 209 return 1
211 210
212 211
213 212 def filterByHeights(self, window):
214 213
215 214 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
216 215
217 216 if window == None:
218 217 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
219 218
220 219 newdelta = deltaHeight * window
221 220 r = self.dataOut.nHeights % window
222 221 newheights = (self.dataOut.nHeights-r)/window
223 222
224 223 if newheights <= 1:
225 224 raise ValueError, "filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(self.dataOut.nHeights, window)
226 225
227 226 if self.dataOut.flagDataAsBlock:
228 227 """
229 228 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
230 229 """
231 230 buffer = self.dataOut.data[:, :, 0:self.dataOut.nHeights-r]
232 231 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nProfiles,self.dataOut.nHeights/window,window)
233 232 buffer = numpy.sum(buffer,3)
234 233
235 234 else:
236 235 buffer = self.dataOut.data[:,0:self.dataOut.nHeights-r]
237 236 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nHeights/window,window)
238 237 buffer = numpy.sum(buffer,2)
239 238
240 239 self.dataOut.data = buffer
241 240 self.dataOut.heightList = self.dataOut.heightList[0] + numpy.arange( newheights )*newdelta
242 241 self.dataOut.windowOfFilter = window
243 242
244 243 def setH0(self, h0, deltaHeight = None):
245 244
246 245 if not deltaHeight:
247 246 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
248 247
249 248 nHeights = self.dataOut.nHeights
250 249
251 250 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
252 251
253 252 self.dataOut.heightList = newHeiRange
254 253
255 254 def deFlip(self, channelList = []):
256 255
257 256 data = self.dataOut.data.copy()
258 257
259 258 if self.dataOut.flagDataAsBlock:
260 259 flip = self.flip
261 260 profileList = range(self.dataOut.nProfiles)
262 261
263 262 if not channelList:
264 263 for thisProfile in profileList:
265 264 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
266 265 flip *= -1.0
267 266 else:
268 267 for thisChannel in channelList:
269 268 if thisChannel not in self.dataOut.channelList:
270 269 continue
271 270
272 271 for thisProfile in profileList:
273 272 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
274 273 flip *= -1.0
275 274
276 275 self.flip = flip
277 276
278 277 else:
279 278 if not channelList:
280 279 data[:,:] = data[:,:]*self.flip
281 280 else:
282 281 for thisChannel in channelList:
283 282 if thisChannel not in self.dataOut.channelList:
284 283 continue
285 284
286 285 data[thisChannel,:] = data[thisChannel,:]*self.flip
287 286
288 287 self.flip *= -1.
289 288
290 289 self.dataOut.data = data
291 290
292 291 def setRadarFrequency(self, frequency=None):
293 292
294 293 if frequency != None:
295 294 self.dataOut.frequency = frequency
296 295
297 296 return 1
298 297
299 298 def interpolateHeights(self, topLim, botLim):
300 299 #69 al 72 para julia
301 300 #82-84 para meteoros
302 301 if len(numpy.shape(self.dataOut.data))==2:
303 302 sampInterp = (self.dataOut.data[:,botLim-1] + self.dataOut.data[:,topLim+1])/2
304 303 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
305 304 #self.dataOut.data[:,botLim:limSup+1] = sampInterp
306 305 self.dataOut.data[:,botLim:topLim+1] = sampInterp
307 306 else:
308 307 nHeights = self.dataOut.data.shape[2]
309 308 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
310 309 y = self.dataOut.data[:,:,range(botLim)+range(topLim+1,nHeights)]
311 310 f = interpolate.interp1d(x, y, axis = 2)
312 311 xnew = numpy.arange(botLim,topLim+1)
313 312 ynew = f(xnew)
314 313
315 314 self.dataOut.data[:,:,botLim:topLim+1] = ynew
316 315
317 316 # import collections
318 317
319 318 class CohInt(Operation):
320 319
321 320 isConfig = False
322 321
323 322 __profIndex = 0
324 323 __withOverapping = False
325 324
326 325 __byTime = False
327 326 __initime = None
328 327 __lastdatatime = None
329 328 __integrationtime = None
330 329
331 330 __buffer = None
332 331
333 332 __dataReady = False
334 333
335 334 n = None
336 335
337 336 def __init__(self, **kwargs):
338 337
339 338 Operation.__init__(self, **kwargs)
340 339
341 340 # self.isConfig = False
342 341
343 342 def setup(self, n=None, timeInterval=None, overlapping=False, byblock=False):
344 343 """
345 344 Set the parameters of the integration class.
346 345
347 346 Inputs:
348 347
349 348 n : Number of coherent integrations
350 349 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
351 350 overlapping :
352 351 """
353 352
354 353 self.__initime = None
355 354 self.__lastdatatime = 0
356 355 self.__buffer = None
357 356 self.__dataReady = False
358 357 self.byblock = byblock
359 358
360 359 if n == None and timeInterval == None:
361 360 raise ValueError, "n or timeInterval should be specified ..."
362 361
363 362 if n != None:
364 363 self.n = n
365 364 self.__byTime = False
366 365 else:
367 366 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
368 367 self.n = 9999
369 368 self.__byTime = True
370 369
371 370 if overlapping:
372 371 self.__withOverapping = True
373 372 self.__buffer = None
374 373 else:
375 374 self.__withOverapping = False
376 375 self.__buffer = 0
377 376
378 377 self.__profIndex = 0
379 378
380 379 def putData(self, data):
381 380
382 381 """
383 382 Add a profile to the __buffer and increase in one the __profileIndex
384 383
385 384 """
386 385
387 386 if not self.__withOverapping:
388 387 self.__buffer += data.copy()
389 388 self.__profIndex += 1
390 389 return
391 390
392 391 #Overlapping data
393 392 nChannels, nHeis = data.shape
394 393 data = numpy.reshape(data, (1, nChannels, nHeis))
395 394
396 395 #If the buffer is empty then it takes the data value
397 396 if self.__buffer is None:
398 397 self.__buffer = data
399 398 self.__profIndex += 1
400 399 return
401 400
402 401 #If the buffer length is lower than n then stakcing the data value
403 402 if self.__profIndex < self.n:
404 403 self.__buffer = numpy.vstack((self.__buffer, data))
405 404 self.__profIndex += 1
406 405 return
407 406
408 407 #If the buffer length is equal to n then replacing the last buffer value with the data value
409 408 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
410 409 self.__buffer[self.n-1] = data
411 410 self.__profIndex = self.n
412 411 return
413 412
414 413
415 414 def pushData(self):
416 415 """
417 416 Return the sum of the last profiles and the profiles used in the sum.
418 417
419 418 Affected:
420 419
421 420 self.__profileIndex
422 421
423 422 """
424 423
425 424 if not self.__withOverapping:
426 425 data = self.__buffer
427 426 n = self.__profIndex
428 427
429 428 self.__buffer = 0
430 429 self.__profIndex = 0
431 430
432 431 return data, n
433 432
434 433 #Integration with Overlapping
435 434 data = numpy.sum(self.__buffer, axis=0)
436 435 n = self.__profIndex
437 436
438 437 return data, n
439 438
440 439 def byProfiles(self, data):
441 440
442 441 self.__dataReady = False
443 442 avgdata = None
444 443 # n = None
445 444
446 445 self.putData(data)
447 446
448 447 if self.__profIndex == self.n:
449 448
450 449 avgdata, n = self.pushData()
451 450 self.__dataReady = True
452 451
453 452 return avgdata
454 453
455 454 def byTime(self, data, datatime):
456 455
457 456 self.__dataReady = False
458 457 avgdata = None
459 458 n = None
460 459
461 460 self.putData(data)
462 461
463 462 if (datatime - self.__initime) >= self.__integrationtime:
464 463 avgdata, n = self.pushData()
465 464 self.n = n
466 465 self.__dataReady = True
467 466
468 467 return avgdata
469 468
470 469 def integrate(self, data, datatime=None):
471 470
472 471 if self.__initime == None:
473 472 self.__initime = datatime
474 473
475 474 if self.__byTime:
476 475 avgdata = self.byTime(data, datatime)
477 476 else:
478 477 avgdata = self.byProfiles(data)
479 478
480 479
481 480 self.__lastdatatime = datatime
482 481
483 482 if avgdata is None:
484 483 return None, None
485 484
486 485 avgdatatime = self.__initime
487 486
488 487 deltatime = datatime -self.__lastdatatime
489 488
490 489 if not self.__withOverapping:
491 490 self.__initime = datatime
492 491 else:
493 492 self.__initime += deltatime
494 493
495 494 return avgdata, avgdatatime
496 495
497 496 def integrateByBlock(self, dataOut):
498 497
499 498 times = int(dataOut.data.shape[1]/self.n)
500 499 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
501 500
502 501 id_min = 0
503 502 id_max = self.n
504 503
505 504 for i in range(times):
506 505 junk = dataOut.data[:,id_min:id_max,:]
507 506 avgdata[:,i,:] = junk.sum(axis=1)
508 507 id_min += self.n
509 508 id_max += self.n
510 509
511 510 timeInterval = dataOut.ippSeconds*self.n
512 511 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
513 512 self.__dataReady = True
514 513 return avgdata, avgdatatime
515 514
516 515
517 516 def run(self, dataOut, n=None, timeInterval=None, overlapping=False, byblock=False, **kwargs):
518 517 if not self.isConfig:
519 518 self.setup(n=n, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
520 519 self.isConfig = True
521 520
522 521 if dataOut.flagDataAsBlock:
523 522 """
524 523 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
525 524 """
526 525 avgdata, avgdatatime = self.integrateByBlock(dataOut)
527 526 dataOut.nProfiles /= self.n
528 527 else:
529 528 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
530 529
531 530 # dataOut.timeInterval *= n
532 531 dataOut.flagNoData = True
533 532
534 533 if self.__dataReady:
535 534 dataOut.data = avgdata
536 535 dataOut.nCohInt *= self.n
537 536 dataOut.utctime = avgdatatime
538 537 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
539 538 dataOut.flagNoData = False
540 539
541 540 class Decoder(Operation):
542 541
543 542 isConfig = False
544 543 __profIndex = 0
545 544
546 545 code = None
547 546
548 547 nCode = None
549 548 nBaud = None
550 549
551 550 def __init__(self, **kwargs):
552 551
553 552 Operation.__init__(self, **kwargs)
554 553
555 554 self.times = None
556 555 self.osamp = None
557 556 # self.__setValues = False
558 557 self.isConfig = False
559 558
560 559 def setup(self, code, osamp, dataOut):
561 560
562 561 self.__profIndex = 0
563 562
564 563 self.code = code
565 564
566 565 self.nCode = len(code)
567 566 self.nBaud = len(code[0])
568 567
569 568 if (osamp != None) and (osamp >1):
570 569 self.osamp = osamp
571 570 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
572 571 self.nBaud = self.nBaud*self.osamp
573 572
574 573 self.__nChannels = dataOut.nChannels
575 574 self.__nProfiles = dataOut.nProfiles
576 575 self.__nHeis = dataOut.nHeights
577 576
578 577 if self.__nHeis < self.nBaud:
579 578 raise ValueError, 'Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud)
580 579
581 580 #Frequency
582 581 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
583 582
584 583 __codeBuffer[:,0:self.nBaud] = self.code
585 584
586 585 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
587 586
588 587 if dataOut.flagDataAsBlock:
589 588
590 589 self.ndatadec = self.__nHeis #- self.nBaud + 1
591 590
592 591 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
593 592
594 593 else:
595 594
596 595 #Time
597 596 self.ndatadec = self.__nHeis #- self.nBaud + 1
598 597
599 598 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
600 599
601 600 def __convolutionInFreq(self, data):
602 601
603 602 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
604 603
605 604 fft_data = numpy.fft.fft(data, axis=1)
606 605
607 606 conv = fft_data*fft_code
608 607
609 608 data = numpy.fft.ifft(conv,axis=1)
610 609
611 610 return data
612 611
613 612 def __convolutionInFreqOpt(self, data):
614 613
615 614 raise NotImplementedError
616 615
617 616 def __convolutionInTime(self, data):
618 617
619 618 code = self.code[self.__profIndex]
620 619
621 620 for i in range(self.__nChannels):
622 621 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
623 622
624 623 return self.datadecTime
625 624
626 #@profile
627 def oldCorrelate(self, i, data, code_block):
628 profilesList = xrange(self.__nProfiles)
629 for j in profilesList:
630 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
631
632 #@profile
633 625 def __convolutionByBlockInTime(self, data):
634 626
635 627 repetitions = self.__nProfiles / self.nCode
636 628
637 629 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
638 630 junk = junk.flatten()
639 631 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
640 632 profilesList = xrange(self.__nProfiles)
641 633
642 # def toVectorize(a,b):
643 # return numpy.correlate(a,b, mode='full')
644 # vectorized = numpy.vectorize(toVectorize, signature='(n),(m)->(k)')
645 for i in range(self.__nChannels):
646 # self.datadecTime[i,:,:] = numpy.array([numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:] for j in profilesList ])
647 # def func(i, j):
648 # self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
649 # map(lambda j: func(i, j), range(self.__nProfiles))
650 #print data[i,:,:].shape
651 # self.datadecTime[i,:,:] = vectorized(data[i,:,:], code_block[:,:])[:,self.nBaud-1:]
652 for j in profilesList:
653 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
654 # print data[i,:,:]
655 # print cSchain.correlateByBlock(data[i,:,:], code_block, 2)
656 # self.datadecTime[i,:,:] = cSchain.correlateByBlock(data[i,:,:], code_block, 2)
657 # print self.datadecTime[i,:,:]
658 #print self.datadecTime[i,:,:].shape
659 return self.datadecTime
660
634 for i in range(self.__nChannels):
635 for j in profilesList:
636 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
637 return self.datadecTime
661 638
662 639 def __convolutionByBlockInFreq(self, data):
663 640
664 641 raise NotImplementedError, "Decoder by frequency fro Blocks not implemented"
665 642
666 643
667 644 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
668 645
669 646 fft_data = numpy.fft.fft(data, axis=2)
670 647
671 648 conv = fft_data*fft_code
672 649
673 650 data = numpy.fft.ifft(conv,axis=2)
674 651
675 652 return data
676 653
677 654
678 655 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
679 656
680 657 if dataOut.flagDecodeData:
681 658 print "This data is already decoded, recoding again ..."
682 659
683 660 if not self.isConfig:
684 661
685 662 if code is None:
686 663 if dataOut.code is None:
687 664 raise ValueError, "Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type
688 665
689 666 code = dataOut.code
690 667 else:
691 668 code = numpy.array(code).reshape(nCode,nBaud)
692 669
693 670 self.setup(code, osamp, dataOut)
694 671
695 672 self.isConfig = True
696 673
697 674 if mode == 3:
698 675 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
699 676
700 677 if times != None:
701 678 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
702 679
703 680 if self.code is None:
704 681 print "Fail decoding: Code is not defined."
705 682 return
706 683
707 684 self.__nProfiles = dataOut.nProfiles
708 685 datadec = None
709 686
710 687 if mode == 3:
711 688 mode = 0
712 689
713 690 if dataOut.flagDataAsBlock:
714 691 """
715 692 Decoding when data have been read as block,
716 693 """
717 694
718 695 if mode == 0:
719 696 datadec = self.__convolutionByBlockInTime(dataOut.data)
720 697 if mode == 1:
721 698 datadec = self.__convolutionByBlockInFreq(dataOut.data)
722 699 else:
723 700 """
724 701 Decoding when data have been read profile by profile
725 702 """
726 703 if mode == 0:
727 704 datadec = self.__convolutionInTime(dataOut.data)
728 705
729 706 if mode == 1:
730 707 datadec = self.__convolutionInFreq(dataOut.data)
731 708
732 709 if mode == 2:
733 710 datadec = self.__convolutionInFreqOpt(dataOut.data)
734 711
735 712 if datadec is None:
736 713 raise ValueError, "Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode
737 714
738 715 dataOut.code = self.code
739 716 dataOut.nCode = self.nCode
740 717 dataOut.nBaud = self.nBaud
741 718
742 719 dataOut.data = datadec
743 720
744 721 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
745 722
746 723 dataOut.flagDecodeData = True #asumo q la data esta decodificada
747 724
748 725 if self.__profIndex == self.nCode-1:
749 726 self.__profIndex = 0
750 727 return 1
751 728
752 729 self.__profIndex += 1
753 730
754 731 return 1
755 732 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
756 733
757 734
758 735 class ProfileConcat(Operation):
759 736
760 737 isConfig = False
761 738 buffer = None
762 739
763 740 def __init__(self, **kwargs):
764 741
765 742 Operation.__init__(self, **kwargs)
766 743 self.profileIndex = 0
767 744
768 745 def reset(self):
769 746 self.buffer = numpy.zeros_like(self.buffer)
770 747 self.start_index = 0
771 748 self.times = 1
772 749
773 750 def setup(self, data, m, n=1):
774 751 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
775 752 self.nHeights = data.shape[1]#.nHeights
776 753 self.start_index = 0
777 754 self.times = 1
778 755
779 756 def concat(self, data):
780 757
781 758 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
782 759 self.start_index = self.start_index + self.nHeights
783 760
784 761 def run(self, dataOut, m):
785 762
786 763 dataOut.flagNoData = True
787 764
788 765 if not self.isConfig:
789 766 self.setup(dataOut.data, m, 1)
790 767 self.isConfig = True
791 768
792 769 if dataOut.flagDataAsBlock:
793 770 raise ValueError, "ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False"
794 771
795 772 else:
796 773 self.concat(dataOut.data)
797 774 self.times += 1
798 775 if self.times > m:
799 776 dataOut.data = self.buffer
800 777 self.reset()
801 778 dataOut.flagNoData = False
802 779 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
803 780 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
804 781 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
805 782 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
806 783 dataOut.ippSeconds *= m
807 784
808 785 class ProfileSelector(Operation):
809 786
810 787 profileIndex = None
811 788 # Tamanho total de los perfiles
812 789 nProfiles = None
813 790
814 791 def __init__(self, **kwargs):
815 792
816 793 Operation.__init__(self, **kwargs)
817 794 self.profileIndex = 0
818 795
819 796 def incProfileIndex(self):
820 797
821 798 self.profileIndex += 1
822 799
823 800 if self.profileIndex >= self.nProfiles:
824 801 self.profileIndex = 0
825 802
826 803 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
827 804
828 805 if profileIndex < minIndex:
829 806 return False
830 807
831 808 if profileIndex > maxIndex:
832 809 return False
833 810
834 811 return True
835 812
836 813 def isThisProfileInList(self, profileIndex, profileList):
837 814
838 815 if profileIndex not in profileList:
839 816 return False
840 817
841 818 return True
842 819
843 820 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
844 821
845 822 """
846 823 ProfileSelector:
847 824
848 825 Inputs:
849 826 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
850 827
851 828 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
852 829
853 830 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
854 831
855 832 """
856 833
857 834 if rangeList is not None:
858 835 if type(rangeList[0]) not in (tuple, list):
859 836 rangeList = [rangeList]
860 837
861 838 dataOut.flagNoData = True
862 839
863 840 if dataOut.flagDataAsBlock:
864 841 """
865 842 data dimension = [nChannels, nProfiles, nHeis]
866 843 """
867 844 if profileList != None:
868 845 dataOut.data = dataOut.data[:,profileList,:]
869 846
870 847 if profileRangeList != None:
871 848 minIndex = profileRangeList[0]
872 849 maxIndex = profileRangeList[1]
873 850 profileList = range(minIndex, maxIndex+1)
874 851
875 852 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
876 853
877 854 if rangeList != None:
878 855
879 856 profileList = []
880 857
881 858 for thisRange in rangeList:
882 859 minIndex = thisRange[0]
883 860 maxIndex = thisRange[1]
884 861
885 862 profileList.extend(range(minIndex, maxIndex+1))
886 863
887 864 dataOut.data = dataOut.data[:,profileList,:]
888 865
889 866 dataOut.nProfiles = len(profileList)
890 867 dataOut.profileIndex = dataOut.nProfiles - 1
891 868 dataOut.flagNoData = False
892 869
893 870 return True
894 871
895 872 """
896 873 data dimension = [nChannels, nHeis]
897 874 """
898 875
899 876 if profileList != None:
900 877
901 878 if self.isThisProfileInList(dataOut.profileIndex, profileList):
902 879
903 880 self.nProfiles = len(profileList)
904 881 dataOut.nProfiles = self.nProfiles
905 882 dataOut.profileIndex = self.profileIndex
906 883 dataOut.flagNoData = False
907 884
908 885 self.incProfileIndex()
909 886 return True
910 887
911 888 if profileRangeList != None:
912 889
913 890 minIndex = profileRangeList[0]
914 891 maxIndex = profileRangeList[1]
915 892
916 893 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
917 894
918 895 self.nProfiles = maxIndex - minIndex + 1
919 896 dataOut.nProfiles = self.nProfiles
920 897 dataOut.profileIndex = self.profileIndex
921 898 dataOut.flagNoData = False
922 899
923 900 self.incProfileIndex()
924 901 return True
925 902
926 903 if rangeList != None:
927 904
928 905 nProfiles = 0
929 906
930 907 for thisRange in rangeList:
931 908 minIndex = thisRange[0]
932 909 maxIndex = thisRange[1]
933 910
934 911 nProfiles += maxIndex - minIndex + 1
935 912
936 913 for thisRange in rangeList:
937 914
938 915 minIndex = thisRange[0]
939 916 maxIndex = thisRange[1]
940 917
941 918 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
942 919
943 920 self.nProfiles = nProfiles
944 921 dataOut.nProfiles = self.nProfiles
945 922 dataOut.profileIndex = self.profileIndex
946 923 dataOut.flagNoData = False
947 924
948 925 self.incProfileIndex()
949 926
950 927 break
951 928
952 929 return True
953 930
954 931
955 932 if beam != None: #beam is only for AMISR data
956 933 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
957 934 dataOut.flagNoData = False
958 935 dataOut.profileIndex = self.profileIndex
959 936
960 937 self.incProfileIndex()
961 938
962 939 return True
963 940
964 941 raise ValueError, "ProfileSelector needs profileList, profileRangeList or rangeList parameter"
965 942
966 943 return False
967 944
968 945 class Reshaper(Operation):
969 946
970 947 def __init__(self, **kwargs):
971 948
972 949 Operation.__init__(self, **kwargs)
973 950
974 951 self.__buffer = None
975 952 self.__nitems = 0
976 953
977 954 def __appendProfile(self, dataOut, nTxs):
978 955
979 956 if self.__buffer is None:
980 957 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
981 958 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
982 959
983 960 ini = dataOut.nHeights * self.__nitems
984 961 end = ini + dataOut.nHeights
985 962
986 963 self.__buffer[:, ini:end] = dataOut.data
987 964
988 965 self.__nitems += 1
989 966
990 967 return int(self.__nitems*nTxs)
991 968
992 969 def __getBuffer(self):
993 970
994 971 if self.__nitems == int(1./self.__nTxs):
995 972
996 973 self.__nitems = 0
997 974
998 975 return self.__buffer.copy()
999 976
1000 977 return None
1001 978
1002 979 def __checkInputs(self, dataOut, shape, nTxs):
1003 980
1004 981 if shape is None and nTxs is None:
1005 982 raise ValueError, "Reshaper: shape of factor should be defined"
1006 983
1007 984 if nTxs:
1008 985 if nTxs < 0:
1009 986 raise ValueError, "nTxs should be greater than 0"
1010 987
1011 988 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1012 989 raise ValueError, "nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs))
1013 990
1014 991 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1015 992
1016 993 return shape, nTxs
1017 994
1018 995 if len(shape) != 2 and len(shape) != 3:
1019 996 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)
1020 997
1021 998 if len(shape) == 2:
1022 999 shape_tuple = [dataOut.nChannels]
1023 1000 shape_tuple.extend(shape)
1024 1001 else:
1025 1002 shape_tuple = list(shape)
1026 1003
1027 1004 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1028 1005
1029 1006 return shape_tuple, nTxs
1030 1007
1031 1008 def run(self, dataOut, shape=None, nTxs=None):
1032 1009
1033 1010 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1034 1011
1035 1012 dataOut.flagNoData = True
1036 1013 profileIndex = None
1037 1014
1038 1015 if dataOut.flagDataAsBlock:
1039 1016
1040 1017 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1041 1018 dataOut.flagNoData = False
1042 1019
1043 1020 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1044 1021
1045 1022 else:
1046 1023
1047 1024 if self.__nTxs < 1:
1048 1025
1049 1026 self.__appendProfile(dataOut, self.__nTxs)
1050 1027 new_data = self.__getBuffer()
1051 1028
1052 1029 if new_data is not None:
1053 1030 dataOut.data = new_data
1054 1031 dataOut.flagNoData = False
1055 1032
1056 1033 profileIndex = dataOut.profileIndex*nTxs
1057 1034
1058 1035 else:
1059 1036 raise ValueError, "nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)"
1060 1037
1061 1038 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1062 1039
1063 1040 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1064 1041
1065 1042 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1066 1043
1067 1044 dataOut.profileIndex = profileIndex
1068 1045
1069 1046 dataOut.ippSeconds /= self.__nTxs
1070 1047
1071 1048 class SplitProfiles(Operation):
1072 1049
1073 1050 def __init__(self, **kwargs):
1074 1051
1075 1052 Operation.__init__(self, **kwargs)
1076 1053
1077 1054 def run(self, dataOut, n):
1078 1055
1079 1056 dataOut.flagNoData = True
1080 1057 profileIndex = None
1081 1058
1082 1059 if dataOut.flagDataAsBlock:
1083 1060
1084 1061 #nchannels, nprofiles, nsamples
1085 1062 shape = dataOut.data.shape
1086 1063
1087 1064 if shape[2] % n != 0:
1088 1065 raise ValueError, "Could not split the data, n=%d has to be multiple of %d" %(n, shape[2])
1089 1066
1090 1067 new_shape = shape[0], shape[1]*n, shape[2]/n
1091 1068
1092 1069 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1093 1070 dataOut.flagNoData = False
1094 1071
1095 1072 profileIndex = int(dataOut.nProfiles/n) - 1
1096 1073
1097 1074 else:
1098 1075
1099 1076 raise ValueError, "Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)"
1100 1077
1101 1078 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1102 1079
1103 1080 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1104 1081
1105 1082 dataOut.nProfiles = int(dataOut.nProfiles*n)
1106 1083
1107 1084 dataOut.profileIndex = profileIndex
1108 1085
1109 1086 dataOut.ippSeconds /= n
1110 1087
1111 1088 class CombineProfiles(Operation):
1112 1089 def __init__(self, **kwargs):
1113 1090
1114 1091 Operation.__init__(self, **kwargs)
1115 1092
1116 1093 self.__remData = None
1117 1094 self.__profileIndex = 0
1118 1095
1119 1096 def run(self, dataOut, n):
1120 1097
1121 1098 dataOut.flagNoData = True
1122 1099 profileIndex = None
1123 1100
1124 1101 if dataOut.flagDataAsBlock:
1125 1102
1126 1103 #nchannels, nprofiles, nsamples
1127 1104 shape = dataOut.data.shape
1128 1105 new_shape = shape[0], shape[1]/n, shape[2]*n
1129 1106
1130 1107 if shape[1] % n != 0:
1131 1108 raise ValueError, "Could not split the data, n=%d has to be multiple of %d" %(n, shape[1])
1132 1109
1133 1110 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1134 1111 dataOut.flagNoData = False
1135 1112
1136 1113 profileIndex = int(dataOut.nProfiles*n) - 1
1137 1114
1138 1115 else:
1139 1116
1140 1117 #nchannels, nsamples
1141 1118 if self.__remData is None:
1142 1119 newData = dataOut.data
1143 1120 else:
1144 1121 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1145 1122
1146 1123 self.__profileIndex += 1
1147 1124
1148 1125 if self.__profileIndex < n:
1149 1126 self.__remData = newData
1150 1127 #continue
1151 1128 return
1152 1129
1153 1130 self.__profileIndex = 0
1154 1131 self.__remData = None
1155 1132
1156 1133 dataOut.data = newData
1157 1134 dataOut.flagNoData = False
1158 1135
1159 1136 profileIndex = dataOut.profileIndex/n
1160 1137
1161 1138
1162 1139 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1163 1140
1164 1141 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1165 1142
1166 1143 dataOut.nProfiles = int(dataOut.nProfiles/n)
1167 1144
1168 1145 dataOut.profileIndex = profileIndex
1169 1146
1170 1147 dataOut.ippSeconds *= n
1171 1148
1172 1149 # import collections
1173 1150 # from scipy.stats import mode
1174 1151 #
1175 1152 # class Synchronize(Operation):
1176 1153 #
1177 1154 # isConfig = False
1178 1155 # __profIndex = 0
1179 1156 #
1180 1157 # def __init__(self, **kwargs):
1181 1158 #
1182 1159 # Operation.__init__(self, **kwargs)
1183 1160 # # self.isConfig = False
1184 1161 # self.__powBuffer = None
1185 1162 # self.__startIndex = 0
1186 1163 # self.__pulseFound = False
1187 1164 #
1188 1165 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1189 1166 #
1190 1167 # #Read data
1191 1168 #
1192 1169 # powerdB = dataOut.getPower(channel = channel)
1193 1170 # noisedB = dataOut.getNoise(channel = channel)[0]
1194 1171 #
1195 1172 # self.__powBuffer.extend(powerdB.flatten())
1196 1173 #
1197 1174 # dataArray = numpy.array(self.__powBuffer)
1198 1175 #
1199 1176 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1200 1177 #
1201 1178 # maxValue = numpy.nanmax(filteredPower)
1202 1179 #
1203 1180 # if maxValue < noisedB + 10:
1204 1181 # #No se encuentra ningun pulso de transmision
1205 1182 # return None
1206 1183 #
1207 1184 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1208 1185 #
1209 1186 # if len(maxValuesIndex) < 2:
1210 1187 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1211 1188 # return None
1212 1189 #
1213 1190 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1214 1191 #
1215 1192 # #Seleccionar solo valores con un espaciamiento de nSamples
1216 1193 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1217 1194 #
1218 1195 # if len(pulseIndex) < 2:
1219 1196 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1220 1197 # return None
1221 1198 #
1222 1199 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1223 1200 #
1224 1201 # #remover senales que se distancien menos de 10 unidades o muestras
1225 1202 # #(No deberian existir IPP menor a 10 unidades)
1226 1203 #
1227 1204 # realIndex = numpy.where(spacing > 10 )[0]
1228 1205 #
1229 1206 # if len(realIndex) < 2:
1230 1207 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1231 1208 # return None
1232 1209 #
1233 1210 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1234 1211 # realPulseIndex = pulseIndex[realIndex]
1235 1212 #
1236 1213 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1237 1214 #
1238 1215 # print "IPP = %d samples" %period
1239 1216 #
1240 1217 # self.__newNSamples = dataOut.nHeights #int(period)
1241 1218 # self.__startIndex = int(realPulseIndex[0])
1242 1219 #
1243 1220 # return 1
1244 1221 #
1245 1222 #
1246 1223 # def setup(self, nSamples, nChannels, buffer_size = 4):
1247 1224 #
1248 1225 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1249 1226 # maxlen = buffer_size*nSamples)
1250 1227 #
1251 1228 # bufferList = []
1252 1229 #
1253 1230 # for i in range(nChannels):
1254 1231 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1255 1232 # maxlen = buffer_size*nSamples)
1256 1233 #
1257 1234 # bufferList.append(bufferByChannel)
1258 1235 #
1259 1236 # self.__nSamples = nSamples
1260 1237 # self.__nChannels = nChannels
1261 1238 # self.__bufferList = bufferList
1262 1239 #
1263 1240 # def run(self, dataOut, channel = 0):
1264 1241 #
1265 1242 # if not self.isConfig:
1266 1243 # nSamples = dataOut.nHeights
1267 1244 # nChannels = dataOut.nChannels
1268 1245 # self.setup(nSamples, nChannels)
1269 1246 # self.isConfig = True
1270 1247 #
1271 1248 # #Append new data to internal buffer
1272 1249 # for thisChannel in range(self.__nChannels):
1273 1250 # bufferByChannel = self.__bufferList[thisChannel]
1274 1251 # bufferByChannel.extend(dataOut.data[thisChannel])
1275 1252 #
1276 1253 # if self.__pulseFound:
1277 1254 # self.__startIndex -= self.__nSamples
1278 1255 #
1279 1256 # #Finding Tx Pulse
1280 1257 # if not self.__pulseFound:
1281 1258 # indexFound = self.__findTxPulse(dataOut, channel)
1282 1259 #
1283 1260 # if indexFound == None:
1284 1261 # dataOut.flagNoData = True
1285 1262 # return
1286 1263 #
1287 1264 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1288 1265 # self.__pulseFound = True
1289 1266 # self.__startIndex = indexFound
1290 1267 #
1291 1268 # #If pulse was found ...
1292 1269 # for thisChannel in range(self.__nChannels):
1293 1270 # bufferByChannel = self.__bufferList[thisChannel]
1294 1271 # #print self.__startIndex
1295 1272 # x = numpy.array(bufferByChannel)
1296 1273 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1297 1274 #
1298 1275 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1299 1276 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1300 1277 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1301 1278 #
1302 1279 # dataOut.data = self.__arrayBuffer
1303 1280 #
1304 1281 # self.__startIndex += self.__newNSamples
1305 1282 #
1306 1283 # return
Show file before | Show file after | Hide comments
@@ -1,502 +1,501
1 1 '''
2 2 @author: Juan C. Espinoza
3 3 '''
4 4
5 5 import time
6 6 import json
7 7 import numpy
8 8 import paho.mqtt.client as mqtt
9 9 import zmq
10 from profilehooks import profile
11 10 import datetime
12 11 from zmq.utils.monitor import recv_monitor_message
13 12 from functools import wraps
14 13 from threading import Thread
15 14 from multiprocessing import Process
16 15
17 16 from schainpy.model.proc.jroproc_base import Operation, ProcessingUnit
18 17 from schainpy.model.data.jrodata import JROData
19 18
20 19 MAXNUMX = 100
21 20 MAXNUMY = 100
22 21
23 22 class PrettyFloat(float):
24 23 def __repr__(self):
25 24 return '%.2f' % self
26 25
27 26 def roundFloats(obj):
28 27 if isinstance(obj, list):
29 28 return map(roundFloats, obj)
30 29 elif isinstance(obj, float):
31 30 return round(obj, 2)
32 31
33 32 def decimate(z, MAXNUMY):
34 33 # dx = int(len(self.x)/self.__MAXNUMX) + 1
35 34
36 35 dy = int(len(z[0])/MAXNUMY) + 1
37 36
38 37 return z[::, ::dy]
39 38
40 39 class throttle(object):
41 40 """Decorator that prevents a function from being called more than once every
42 41 time period.
43 42 To create a function that cannot be called more than once a minute, but
44 43 will sleep until it can be called:
45 44 @throttle(minutes=1)
46 45 def foo():
47 46 pass
48 47
49 48 for i in range(10):
50 49 foo()
51 50 print "This function has run %s times." % i
52 51 """
53 52
54 53 def __init__(self, seconds=0, minutes=0, hours=0):
55 54 self.throttle_period = datetime.timedelta(
56 55 seconds=seconds, minutes=minutes, hours=hours
57 56 )
58 57
59 58 self.time_of_last_call = datetime.datetime.min
60 59
61 60 def __call__(self, fn):
62 61 @wraps(fn)
63 62 def wrapper(*args, **kwargs):
64 63 now = datetime.datetime.now()
65 64 time_since_last_call = now - self.time_of_last_call
66 65 time_left = self.throttle_period - time_since_last_call
67 66
68 67 if time_left > datetime.timedelta(seconds=0):
69 68 return
70 69
71 70 self.time_of_last_call = datetime.datetime.now()
72 71 return fn(*args, **kwargs)
73 72
74 73 return wrapper
75 74
76 75
77 76 class PublishData(Operation):
78 77 """Clase publish."""
79 78
80 79 def __init__(self, **kwargs):
81 80 """Inicio."""
82 81 Operation.__init__(self, **kwargs)
83 82 self.isConfig = False
84 83 self.client = None
85 84 self.zeromq = None
86 85 self.mqtt = None
87 86
88 87 def on_disconnect(self, client, userdata, rc):
89 88 if rc != 0:
90 89 print("Unexpected disconnection.")
91 90 self.connect()
92 91
93 92 def connect(self):
94 93 print 'trying to connect'
95 94 try:
96 95 self.client.connect(
97 96 host=self.host,
98 97 port=self.port,
99 98 keepalive=60*10,
100 99 bind_address='')
101 100 self.client.loop_start()
102 101 # self.client.publish(
103 102 # self.topic + 'SETUP',
104 103 # json.dumps(setup),
105 104 # retain=True
106 105 # )
107 106 except:
108 107 print "MQTT Conection error."
109 108 self.client = False
110 109
111 110 def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, verbose=True, **kwargs):
112 111 self.counter = 0
113 112 self.topic = kwargs.get('topic', 'schain')
114 113 self.delay = kwargs.get('delay', 0)
115 114 self.plottype = kwargs.get('plottype', 'spectra')
116 115 self.host = kwargs.get('host', "10.10.10.82")
117 116 self.port = kwargs.get('port', 3000)
118 117 self.clientId = clientId
119 118 self.cnt = 0
120 119 self.zeromq = zeromq
121 120 self.mqtt = kwargs.get('plottype', 0)
122 121 self.client = None
123 122 self.verbose = verbose
124 123 self.dataOut.firstdata = True
125 124 setup = []
126 125 if mqtt is 1:
127 126 self.client = mqtt.Client(
128 127 client_id=self.clientId + self.topic + 'SCHAIN',
129 128 clean_session=True)
130 129 self.client.on_disconnect = self.on_disconnect
131 130 self.connect()
132 131 for plot in self.plottype:
133 132 setup.append({
134 133 'plot': plot,
135 134 'topic': self.topic + plot,
136 135 'title': getattr(self, plot + '_' + 'title', False),
137 136 'xlabel': getattr(self, plot + '_' + 'xlabel', False),
138 137 'ylabel': getattr(self, plot + '_' + 'ylabel', False),
139 138 'xrange': getattr(self, plot + '_' + 'xrange', False),
140 139 'yrange': getattr(self, plot + '_' + 'yrange', False),
141 140 'zrange': getattr(self, plot + '_' + 'zrange', False),
142 141 })
143 142 if zeromq is 1:
144 143 context = zmq.Context()
145 144 self.zmq_socket = context.socket(zmq.PUSH)
146 145 server = kwargs.get('server', 'zmq.pipe')
147 146
148 147 if 'tcp://' in server:
149 148 address = server
150 149 else:
151 150 address = 'ipc:///tmp/%s' % server
152 151
153 152 self.zmq_socket.connect(address)
154 153 time.sleep(1)
155 154
156 155
157 156 def publish_data(self):
158 157 self.dataOut.finished = False
159 158 if self.mqtt is 1:
160 159 yData = self.dataOut.heightList[:2].tolist()
161 160 if self.plottype == 'spectra':
162 161 data = getattr(self.dataOut, 'data_spc')
163 162 z = data/self.dataOut.normFactor
164 163 zdB = 10*numpy.log10(z)
165 164 xlen, ylen = zdB[0].shape
166 165 dx = int(xlen/MAXNUMX) + 1
167 166 dy = int(ylen/MAXNUMY) + 1
168 167 Z = [0 for i in self.dataOut.channelList]
169 168 for i in self.dataOut.channelList:
170 169 Z[i] = zdB[i][::dx, ::dy].tolist()
171 170 payload = {
172 171 'timestamp': self.dataOut.utctime,
173 172 'data': roundFloats(Z),
174 173 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
175 174 'interval': self.dataOut.getTimeInterval(),
176 175 'type': self.plottype,
177 176 'yData': yData
178 177 }
179 178 # print payload
180 179
181 180 elif self.plottype in ('rti', 'power'):
182 181 data = getattr(self.dataOut, 'data_spc')
183 182 z = data/self.dataOut.normFactor
184 183 avg = numpy.average(z, axis=1)
185 184 avgdB = 10*numpy.log10(avg)
186 185 xlen, ylen = z[0].shape
187 186 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
188 187 AVG = [0 for i in self.dataOut.channelList]
189 188 for i in self.dataOut.channelList:
190 189 AVG[i] = avgdB[i][::dy].tolist()
191 190 payload = {
192 191 'timestamp': self.dataOut.utctime,
193 192 'data': roundFloats(AVG),
194 193 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
195 194 'interval': self.dataOut.getTimeInterval(),
196 195 'type': self.plottype,
197 196 'yData': yData
198 197 }
199 198 elif self.plottype == 'noise':
200 199 noise = self.dataOut.getNoise()/self.dataOut.normFactor
201 200 noisedB = 10*numpy.log10(noise)
202 201 payload = {
203 202 'timestamp': self.dataOut.utctime,
204 203 'data': roundFloats(noisedB.reshape(-1, 1).tolist()),
205 204 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
206 205 'interval': self.dataOut.getTimeInterval(),
207 206 'type': self.plottype,
208 207 'yData': yData
209 208 }
210 209 elif self.plottype == 'snr':
211 210 data = getattr(self.dataOut, 'data_SNR')
212 211 avgdB = 10*numpy.log10(data)
213 212
214 213 ylen = data[0].size
215 214 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
216 215 AVG = [0 for i in self.dataOut.channelList]
217 216 for i in self.dataOut.channelList:
218 217 AVG[i] = avgdB[i][::dy].tolist()
219 218 payload = {
220 219 'timestamp': self.dataOut.utctime,
221 220 'data': roundFloats(AVG),
222 221 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
223 222 'type': self.plottype,
224 223 'yData': yData
225 224 }
226 225 else:
227 226 print "Tipo de grafico invalido"
228 227 payload = {
229 228 'data': 'None',
230 229 'timestamp': 'None',
231 230 'type': None
232 231 }
233 232 # print 'Publishing data to {}'.format(self.host)
234 233 self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0)
235 234
236 235 if self.zeromq is 1:
237 236 if self.verbose:
238 237 print '[Sending] {} - {}'.format(self.dataOut.type, self.dataOut.datatime)
239 238 self.zmq_socket.send_pyobj(self.dataOut)
240 239 self.dataOut.firstdata = False
241 240
242 241
243 242 def run(self, dataOut, **kwargs):
244 243 self.dataOut = dataOut
245 244 if not self.isConfig:
246 245 self.setup(**kwargs)
247 246 self.isConfig = True
248 247
249 248 self.publish_data()
250 249 time.sleep(self.delay)
251 250
252 251 def close(self):
253 252 if self.zeromq is 1:
254 253 self.dataOut.finished = True
255 254 self.zmq_socket.send_pyobj(self.dataOut)
256 255 self.zmq_socket.close()
257 256 if self.client:
258 257 self.client.loop_stop()
259 258 self.client.disconnect()
260 259
261 260
262 261 class ReceiverData(ProcessingUnit):
263 262
264 263 def __init__(self, **kwargs):
265 264
266 265 ProcessingUnit.__init__(self, **kwargs)
267 266
268 267 self.isConfig = False
269 268 server = kwargs.get('server', 'zmq.pipe')
270 269 if 'tcp://' in server:
271 270 address = server
272 271 else:
273 272 address = 'ipc:///tmp/%s' % server
274 273
275 274 self.address = address
276 275 self.dataOut = JROData()
277 276
278 277 def setup(self):
279 278
280 279 self.context = zmq.Context()
281 280 self.receiver = self.context.socket(zmq.PULL)
282 281 self.receiver.bind(self.address)
283 282 time.sleep(0.5)
284 283 print '[Starting] ReceiverData from {}'.format(self.address)
285 284
286 285
287 286 def run(self):
288 287
289 288 if not self.isConfig:
290 289 self.setup()
291 290 self.isConfig = True
292 291
293 292 self.dataOut = self.receiver.recv_pyobj()
294 293 print '[Receiving] {} - {}'.format(self.dataOut.type,
295 294 self.dataOut.datatime.ctime())
296 295
297 296
298 297 class PlotterReceiver(ProcessingUnit, Process):
299 298
300 299 throttle_value = 5
301 300
302 301 def __init__(self, **kwargs):
303 302
304 303 ProcessingUnit.__init__(self, **kwargs)
305 304 Process.__init__(self)
306 305 self.mp = False
307 306 self.isConfig = False
308 307 self.isWebConfig = False
309 308 self.plottypes = []
310 309 self.connections = 0
311 310 server = kwargs.get('server', 'zmq.pipe')
312 311 plot_server = kwargs.get('plot_server', 'zmq.web')
313 312 if 'tcp://' in server:
314 313 address = server
315 314 else:
316 315 address = 'ipc:///tmp/%s' % server
317 316
318 317 if 'tcp://' in plot_server:
319 318 plot_address = plot_server
320 319 else:
321 320 plot_address = 'ipc:///tmp/%s' % plot_server
322 321
323 322 self.address = address
324 323 self.plot_address = plot_address
325 324 self.plottypes = [s.strip() for s in kwargs.get('plottypes', 'rti').split(',')]
326 325 self.realtime = kwargs.get('realtime', False)
327 326 self.throttle_value = kwargs.get('throttle', 5)
328 327 self.sendData = self.initThrottle(self.throttle_value)
329 328 self.setup()
330 329
331 330 def setup(self):
332 331
333 332 self.data = {}
334 333 self.data['times'] = []
335 334 for plottype in self.plottypes:
336 335 self.data[plottype] = {}
337 336 self.data['noise'] = {}
338 337 self.data['throttle'] = self.throttle_value
339 338 self.data['ENDED'] = False
340 339 self.isConfig = True
341 340 self.data_web = {}
342 341
343 342 def event_monitor(self, monitor):
344 343
345 344 events = {}
346 345
347 346 for name in dir(zmq):
348 347 if name.startswith('EVENT_'):
349 348 value = getattr(zmq, name)
350 349 events[value] = name
351 350
352 351 while monitor.poll():
353 352 evt = recv_monitor_message(monitor)
354 353 if evt['event'] == 32:
355 354 self.connections += 1
356 355 if evt['event'] == 512:
357 356 pass
358 357 if self.connections == 0 and self.started is True:
359 358 self.ended = True
360 359
361 360 evt.update({'description': events[evt['event']]})
362 361
363 362 if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
364 363 break
365 364 monitor.close()
366 365 print("event monitor thread done!")
367 366
368 367 def initThrottle(self, throttle_value):
369 368
370 369 @throttle(seconds=throttle_value)
371 370 def sendDataThrottled(fn_sender, data):
372 371 fn_sender(data)
373 372
374 373 return sendDataThrottled
375 374
376 375
377 376 def send(self, data):
378 377 # print '[sending] data=%s size=%s' % (data.keys(), len(data['times']))
379 378 self.sender.send_pyobj(data)
380 379
381 380
382 381 def update(self):
383 382 t = self.dataOut.utctime
384 383
385 384 if t in self.data['times']:
386 385 return
387 386
388 387 self.data['times'].append(t)
389 388 self.data['dataOut'] = self.dataOut
390 389
391 390 for plottype in self.plottypes:
392 391 if plottype == 'spc':
393 392 z = self.dataOut.data_spc/self.dataOut.normFactor
394 393 self.data[plottype] = 10*numpy.log10(z)
395 394 self.data['noise'][t] = 10*numpy.log10(self.dataOut.getNoise()/self.dataOut.normFactor)
396 395 if plottype == 'cspc':
397 396 jcoherence = self.dataOut.data_cspc/numpy.sqrt(self.dataOut.data_spc*self.dataOut.data_spc)
398 397 self.data['cspc_coh'] = numpy.abs(jcoherence)
399 398 self.data['cspc_phase'] = numpy.arctan2(jcoherence.imag, jcoherence.real)*180/numpy.pi
400 399 if plottype == 'rti':
401 400 self.data[plottype][t] = self.dataOut.getPower()
402 401 if plottype == 'snr':
403 402 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_SNR)
404 403 if plottype == 'dop':
405 404 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_DOP)
406 405 if plottype == 'mean':
407 406 self.data[plottype][t] = self.dataOut.data_MEAN
408 407 if plottype == 'std':
409 408 self.data[plottype][t] = self.dataOut.data_STD
410 409 if plottype == 'coh':
411 410 self.data[plottype][t] = self.dataOut.getCoherence()
412 411 if plottype == 'phase':
413 412 self.data[plottype][t] = self.dataOut.getCoherence(phase=True)
414 413 if plottype == 'output':
415 414 self.data[plottype][t] = self.dataOut.data_output
416 415 if plottype == 'param':
417 416 self.data[plottype][t] = self.dataOut.data_param
418 417 if self.realtime:
419 418 self.data_web['timestamp'] = t
420 419 if plottype == 'spc':
421 420 self.data_web[plottype] = roundFloats(decimate(self.data[plottype]).tolist())
422 421 elif plottype == 'cspc':
423 422 self.data_web['cspc_coh'] = roundFloats(decimate(self.data['cspc_coh']).tolist())
424 423 self.data_web['cspc_phase'] = roundFloats(decimate(self.data['cspc_phase']).tolist())
425 424 elif plottype == 'noise':
426 425 self.data_web['noise'] = roundFloats(self.data['noise'][t].tolist())
427 426 else:
428 427 self.data_web[plottype] = roundFloats(decimate(self.data[plottype][t]).tolist())
429 428 self.data_web['interval'] = self.dataOut.getTimeInterval()
430 429 self.data_web['type'] = plottype
431 430
432 431 def run(self):
433 432
434 433 print '[Starting] {} from {}'.format(self.name, self.address)
435 434
436 435 self.context = zmq.Context()
437 436 self.receiver = self.context.socket(zmq.PULL)
438 437 self.receiver.bind(self.address)
439 438 monitor = self.receiver.get_monitor_socket()
440 439 self.sender = self.context.socket(zmq.PUB)
441 440 if self.realtime:
442 441 self.sender_web = self.context.socket(zmq.PUB)
443 442 self.sender_web.connect(self.plot_address)
444 443 time.sleep(1)
445 444
446 445 if 'server' in self.kwargs:
447 446 self.sender.bind("ipc:///tmp/{}.plots".format(self.kwargs['server']))
448 447 else:
449 448 self.sender.bind("ipc:///tmp/zmq.plots")
450 449
451 450 time.sleep(3)
452 451
453 452 t = Thread(target=self.event_monitor, args=(monitor,))
454 453 t.start()
455 454
456 455 while True:
457 456 self.dataOut = self.receiver.recv_pyobj()
458 457 # print '[Receiving] {} - {}'.format(self.dataOut.type,
459 458 # self.dataOut.datatime.ctime())
460 459
461 460 self.update()
462 461
463 462 if self.dataOut.firstdata is True:
464 463 self.data['STARTED'] = True
465 464
466 465 if self.dataOut.finished is True:
467 466 self.send(self.data)
468 467 self.connections -= 1
469 468 if self.connections == 0 and self.started:
470 469 self.ended = True
471 470 self.data['ENDED'] = True
472 471 self.send(self.data)
473 472 self.setup()
474 473 self.started = False
475 474 else:
476 475 if self.realtime:
477 476 self.send(self.data)
478 477 self.sender_web.send_string(json.dumps(self.data_web))
479 478 else:
480 479 self.sendData(self.send, self.data)
481 480 self.started = True
482 481
483 482 self.data['STARTED'] = False
484 483 return
485 484
486 485 def sendToWeb(self):
487 486
488 487 if not self.isWebConfig:
489 488 context = zmq.Context()
490 489 sender_web_config = context.socket(zmq.PUB)
491 490 if 'tcp://' in self.plot_address:
492 491 dum, address, port = self.plot_address.split(':')
493 492 conf_address = '{}:{}:{}'.format(dum, address, int(port)+1)
494 493 else:
495 494 conf_address = self.plot_address + '.config'
496 495 sender_web_config.bind(conf_address)
497 496 time.sleep(1)
498 497 for kwargs in self.operationKwargs.values():
499 498 if 'plot' in kwargs:
500 499 print '[Sending] Config data to web for {}'.format(kwargs['code'].upper())
501 500 sender_web_config.send_string(json.dumps(kwargs))
502 501 self.isWebConfig = True
General Comments 0
You need to be logged in to leave comments. Login now