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