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