##// END OF EJS Templates
schain
RSS
spectra_lags: nProfiles and nFFTPoints can be different
Miguel Valdez -
r770:c5f19c5ee4fb
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,909 +1,910
1 1 import numpy
2 2
3 3 from jroproc_base import ProcessingUnit, Operation
4 4 from schainpy.model.data.jrodata import Spectra
5 5 from schainpy.model.data.jrodata import hildebrand_sekhon
6 6
7 7 class SpectraLagsProc(ProcessingUnit):
8 8
9 9 def __init__(self):
10 10
11 11 ProcessingUnit.__init__(self)
12 12
13 13 self.buffer = None
14 14 self.firstdatatime = None
15 15 self.profIndex = 0
16 16 self.dataOut = Spectra()
17 17 self.id_min = None
18 18 self.id_max = None
19 19
20 20 def __updateSpecFromVoltage(self):
21 21
22 22 self.dataOut.timeZone = self.dataIn.timeZone
23 23 self.dataOut.dstFlag = self.dataIn.dstFlag
24 24 self.dataOut.errorCount = self.dataIn.errorCount
25 25 self.dataOut.useLocalTime = self.dataIn.useLocalTime
26 26
27 27 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
28 28 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
29 29 self.dataOut.ippSeconds = self.dataIn.getDeltaH()*(10**-6)/0.15
30 30
31 31 self.dataOut.channelList = self.dataIn.channelList
32 32 self.dataOut.heightList = self.dataIn.heightList
33 33 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
34 34
35 35 self.dataOut.nBaud = self.dataIn.nBaud
36 36 self.dataOut.nCode = self.dataIn.nCode
37 37 self.dataOut.code = self.dataIn.code
38 self.dataOut.nProfiles = self.dataOut.nFFTPoints
38 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
39 39
40 40 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
41 41 self.dataOut.utctime = self.firstdatatime
42 42 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
43 43 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
44 44 self.dataOut.flagShiftFFT = False
45 45
46 46 self.dataOut.nCohInt = self.dataIn.nCohInt
47 47 self.dataOut.nIncohInt = 1
48 48
49 49 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
50 50
51 51 self.dataOut.frequency = self.dataIn.frequency
52 52 self.dataOut.realtime = self.dataIn.realtime
53 53
54 54 self.dataOut.azimuth = self.dataIn.azimuth
55 55 self.dataOut.zenith = self.dataIn.zenith
56 56
57 57 self.dataOut.beam.codeList = self.dataIn.beam.codeList
58 58 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
59 59 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
60 60
61 61 def __decodeData(self, nProfiles, code):
62 62
63 63 if code is None:
64 64 return
65 65
66 66 for i in range(nProfiles):
67 67 self.buffer[:,i,:] = self.buffer[:,i,:]*code[0][i]
68 68
69 69 def __getFft(self):
70 70 """
71 71 Convierte valores de Voltaje a Spectra
72 72
73 73 Affected:
74 74 self.dataOut.data_spc
75 75 self.dataOut.data_cspc
76 76 self.dataOut.data_dc
77 77 self.dataOut.heightList
78 78 self.profIndex
79 79 self.buffer
80 80 self.dataOut.flagNoData
81 81 """
82 82 nsegments = self.dataOut.nHeights
83 83
84 84 _fft_buffer = numpy.zeros((self.dataOut.nChannels, self.dataOut.nProfiles, nsegments), dtype='complex')
85 85
86 86 for i in range(nsegments):
87 87 try:
88 88 _fft_buffer[:,:,i] = self.buffer[:,i:i+self.dataOut.nProfiles]
89 89
90 90 if self.code is not None:
91 91 _fft_buffer[:,:,i] = _fft_buffer[:,:,i]*self.code[0]
92 92 except:
93 93 pass
94 94
95 95 fft_volt = numpy.fft.fft(_fft_buffer, n=self.dataOut.nFFTPoints, axis=1)
96 96 fft_volt = fft_volt.astype(numpy.dtype('complex'))
97 97 dc = fft_volt[:,0,:]
98 98
99 99 #calculo de self-spectra
100 100 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
101 101 spc = fft_volt * numpy.conjugate(fft_volt)
102 102 spc = spc.real
103 103
104 104 blocksize = 0
105 105 blocksize += dc.size
106 106 blocksize += spc.size
107 107
108 108 cspc = None
109 109 pairIndex = 0
110 110
111 111 if self.dataOut.pairsList != None:
112 112 #calculo de cross-spectra
113 113 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
114 114 for pair in self.dataOut.pairsList:
115 115 if pair[0] not in self.dataOut.channelList:
116 116 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
117 117 if pair[1] not in self.dataOut.channelList:
118 118 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
119 119
120 120 chan_index0 = self.dataOut.channelList.index(pair[0])
121 121 chan_index1 = self.dataOut.channelList.index(pair[1])
122 122
123 123 cspc[pairIndex,:,:] = fft_volt[chan_index0,:,:] * numpy.conjugate(fft_volt[chan_index1,:,:])
124 124 pairIndex += 1
125 125 blocksize += cspc.size
126 126
127 127 self.dataOut.data_spc = spc
128 128 self.dataOut.data_cspc = cspc
129 129 self.dataOut.data_dc = dc
130 130 self.dataOut.blockSize = blocksize
131 131 self.dataOut.flagShiftFFT = True
132 132
133 133 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], code=None, nCode=1, nBaud=1):
134 134
135 135 self.dataOut.flagNoData = True
136 136
137 137 if code is not None:
138 138 self.code = numpy.array(code).reshape(nCode,nBaud)
139 139 else:
140 140 self.code = None
141 141
142 142 if self.dataIn.type == "Voltage":
143 143
144 144 if nFFTPoints == None:
145 145 raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
146 146
147 147 if nProfiles == None:
148 148 nProfiles = nFFTPoints
149 149
150 150 self.dataOut.ippFactor = 1
151 151
152 152 self.dataOut.nFFTPoints = nFFTPoints
153 self.dataOut.nProfiles = nProfiles
153 154 self.dataOut.pairsList = pairsList
154 155
155 156 # if self.buffer is None:
156 157 # self.buffer = numpy.zeros( (self.dataIn.nChannels, nProfiles, self.dataIn.nHeights),
157 158 # dtype='complex')
158 159
159 160 if not self.dataIn.flagDataAsBlock:
160 161 self.buffer = self.dataIn.data.copy()
161 162
162 163 # for i in range(self.dataIn.nHeights):
163 164 # self.buffer[:, self.profIndex, self.profIndex:] = voltage_data[:,:self.dataIn.nHeights - self.profIndex]
164 165 #
165 166 # self.profIndex += 1
166 167
167 168 else:
168 169 raise ValueError, ""
169 170
170 171 self.firstdatatime = self.dataIn.utctime
171 172
172 173 self.profIndex == nProfiles
173 174
174 175 self.__updateSpecFromVoltage()
175 176
176 177 self.__getFft()
177 178
178 179 self.dataOut.flagNoData = False
179 180
180 181 return True
181 182
182 183 raise ValueError, "The type of input object '%s' is not valid"%(self.dataIn.type)
183 184
184 185 def __selectPairs(self, pairsList):
185 186
186 187 if channelList == None:
187 188 return
188 189
189 190 pairsIndexListSelected = []
190 191
191 192 for thisPair in pairsList:
192 193
193 194 if thisPair not in self.dataOut.pairsList:
194 195 continue
195 196
196 197 pairIndex = self.dataOut.pairsList.index(thisPair)
197 198
198 199 pairsIndexListSelected.append(pairIndex)
199 200
200 201 if not pairsIndexListSelected:
201 202 self.dataOut.data_cspc = None
202 203 self.dataOut.pairsList = []
203 204 return
204 205
205 206 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
206 207 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
207 208
208 209 return
209 210
210 211 def __selectPairsByChannel(self, channelList=None):
211 212
212 213 if channelList == None:
213 214 return
214 215
215 216 pairsIndexListSelected = []
216 217 for pairIndex in self.dataOut.pairsIndexList:
217 218 #First pair
218 219 if self.dataOut.pairsList[pairIndex][0] not in channelList:
219 220 continue
220 221 #Second pair
221 222 if self.dataOut.pairsList[pairIndex][1] not in channelList:
222 223 continue
223 224
224 225 pairsIndexListSelected.append(pairIndex)
225 226
226 227 if not pairsIndexListSelected:
227 228 self.dataOut.data_cspc = None
228 229 self.dataOut.pairsList = []
229 230 return
230 231
231 232 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
232 233 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
233 234
234 235 return
235 236
236 237 def selectChannels(self, channelList):
237 238
238 239 channelIndexList = []
239 240
240 241 for channel in channelList:
241 242 if channel not in self.dataOut.channelList:
242 243 raise ValueError, "Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList))
243 244
244 245 index = self.dataOut.channelList.index(channel)
245 246 channelIndexList.append(index)
246 247
247 248 self.selectChannelsByIndex(channelIndexList)
248 249
249 250 def selectChannelsByIndex(self, channelIndexList):
250 251 """
251 252 Selecciona un bloque de datos en base a canales segun el channelIndexList
252 253
253 254 Input:
254 255 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
255 256
256 257 Affected:
257 258 self.dataOut.data_spc
258 259 self.dataOut.channelIndexList
259 260 self.dataOut.nChannels
260 261
261 262 Return:
262 263 None
263 264 """
264 265
265 266 for channelIndex in channelIndexList:
266 267 if channelIndex not in self.dataOut.channelIndexList:
267 268 raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList)
268 269
269 270 # nChannels = len(channelIndexList)
270 271
271 272 data_spc = self.dataOut.data_spc[channelIndexList,:]
272 273 data_dc = self.dataOut.data_dc[channelIndexList,:]
273 274
274 275 self.dataOut.data_spc = data_spc
275 276 self.dataOut.data_dc = data_dc
276 277
277 278 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
278 279 # self.dataOut.nChannels = nChannels
279 280
280 281 self.__selectPairsByChannel(self.dataOut.channelList)
281 282
282 283 return 1
283 284
284 285 def selectHeights(self, minHei, maxHei):
285 286 """
286 287 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
287 288 minHei <= height <= maxHei
288 289
289 290 Input:
290 291 minHei : valor minimo de altura a considerar
291 292 maxHei : valor maximo de altura a considerar
292 293
293 294 Affected:
294 295 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
295 296
296 297 Return:
297 298 1 si el metodo se ejecuto con exito caso contrario devuelve 0
298 299 """
299 300
300 301 if (minHei > maxHei):
301 302 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
302 303
303 304 if (minHei < self.dataOut.heightList[0]):
304 305 minHei = self.dataOut.heightList[0]
305 306
306 307 if (maxHei > self.dataOut.heightList[-1]):
307 308 maxHei = self.dataOut.heightList[-1]
308 309
309 310 minIndex = 0
310 311 maxIndex = 0
311 312 heights = self.dataOut.heightList
312 313
313 314 inda = numpy.where(heights >= minHei)
314 315 indb = numpy.where(heights <= maxHei)
315 316
316 317 try:
317 318 minIndex = inda[0][0]
318 319 except:
319 320 minIndex = 0
320 321
321 322 try:
322 323 maxIndex = indb[0][-1]
323 324 except:
324 325 maxIndex = len(heights)
325 326
326 327 self.selectHeightsByIndex(minIndex, maxIndex)
327 328
328 329 return 1
329 330
330 331 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
331 332 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
332 333
333 334 if hei_ref != None:
334 335 newheis = numpy.where(self.dataOut.heightList>hei_ref)
335 336
336 337 minIndex = min(newheis[0])
337 338 maxIndex = max(newheis[0])
338 339 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
339 340 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
340 341
341 342 # determina indices
342 343 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
343 344 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
344 345 beacon_dB = numpy.sort(avg_dB)[-nheis:]
345 346 beacon_heiIndexList = []
346 347 for val in avg_dB.tolist():
347 348 if val >= beacon_dB[0]:
348 349 beacon_heiIndexList.append(avg_dB.tolist().index(val))
349 350
350 351 #data_spc = data_spc[:,:,beacon_heiIndexList]
351 352 data_cspc = None
352 353 if self.dataOut.data_cspc is not None:
353 354 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
354 355 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
355 356
356 357 data_dc = None
357 358 if self.dataOut.data_dc is not None:
358 359 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
359 360 #data_dc = data_dc[:,beacon_heiIndexList]
360 361
361 362 self.dataOut.data_spc = data_spc
362 363 self.dataOut.data_cspc = data_cspc
363 364 self.dataOut.data_dc = data_dc
364 365 self.dataOut.heightList = heightList
365 366 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
366 367
367 368 return 1
368 369
369 370
370 371 def selectHeightsByIndex(self, minIndex, maxIndex):
371 372 """
372 373 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
373 374 minIndex <= index <= maxIndex
374 375
375 376 Input:
376 377 minIndex : valor de indice minimo de altura a considerar
377 378 maxIndex : valor de indice maximo de altura a considerar
378 379
379 380 Affected:
380 381 self.dataOut.data_spc
381 382 self.dataOut.data_cspc
382 383 self.dataOut.data_dc
383 384 self.dataOut.heightList
384 385
385 386 Return:
386 387 1 si el metodo se ejecuto con exito caso contrario devuelve 0
387 388 """
388 389
389 390 if (minIndex < 0) or (minIndex > maxIndex):
390 391 raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex)
391 392
392 393 if (maxIndex >= self.dataOut.nHeights):
393 394 maxIndex = self.dataOut.nHeights-1
394 395
395 396 #Spectra
396 397 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
397 398
398 399 data_cspc = None
399 400 if self.dataOut.data_cspc is not None:
400 401 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
401 402
402 403 data_dc = None
403 404 if self.dataOut.data_dc is not None:
404 405 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
405 406
406 407 self.dataOut.data_spc = data_spc
407 408 self.dataOut.data_cspc = data_cspc
408 409 self.dataOut.data_dc = data_dc
409 410
410 411 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
411 412
412 413 return 1
413 414
414 415 def removeDC(self, mode = 2):
415 416 jspectra = self.dataOut.data_spc
416 417 jcspectra = self.dataOut.data_cspc
417 418
418 419
419 420 num_chan = jspectra.shape[0]
420 421 num_hei = jspectra.shape[2]
421 422
422 423 if jcspectra is not None:
423 424 jcspectraExist = True
424 425 num_pairs = jcspectra.shape[0]
425 426 else: jcspectraExist = False
426 427
427 428 freq_dc = jspectra.shape[1]/2
428 429 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
429 430
430 431 if ind_vel[0]<0:
431 432 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
432 433
433 434 if mode == 1:
434 435 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
435 436
436 437 if jcspectraExist:
437 438 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
438 439
439 440 if mode == 2:
440 441
441 442 vel = numpy.array([-2,-1,1,2])
442 443 xx = numpy.zeros([4,4])
443 444
444 445 for fil in range(4):
445 446 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
446 447
447 448 xx_inv = numpy.linalg.inv(xx)
448 449 xx_aux = xx_inv[0,:]
449 450
450 451 for ich in range(num_chan):
451 452 yy = jspectra[ich,ind_vel,:]
452 453 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
453 454
454 455 junkid = jspectra[ich,freq_dc,:]<=0
455 456 cjunkid = sum(junkid)
456 457
457 458 if cjunkid.any():
458 459 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
459 460
460 461 if jcspectraExist:
461 462 for ip in range(num_pairs):
462 463 yy = jcspectra[ip,ind_vel,:]
463 464 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
464 465
465 466
466 467 self.dataOut.data_spc = jspectra
467 468 self.dataOut.data_cspc = jcspectra
468 469
469 470 return 1
470 471
471 472 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
472 473
473 474 jspectra = self.dataOut.data_spc
474 475 jcspectra = self.dataOut.data_cspc
475 476 jnoise = self.dataOut.getNoise()
476 477 num_incoh = self.dataOut.nIncohInt
477 478
478 479 num_channel = jspectra.shape[0]
479 480 num_prof = jspectra.shape[1]
480 481 num_hei = jspectra.shape[2]
481 482
482 483 #hei_interf
483 484 if hei_interf is None:
484 485 count_hei = num_hei/2 #Como es entero no importa
485 486 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
486 487 hei_interf = numpy.asarray(hei_interf)[0]
487 488 #nhei_interf
488 489 if (nhei_interf == None):
489 490 nhei_interf = 5
490 491 if (nhei_interf < 1):
491 492 nhei_interf = 1
492 493 if (nhei_interf > count_hei):
493 494 nhei_interf = count_hei
494 495 if (offhei_interf == None):
495 496 offhei_interf = 0
496 497
497 498 ind_hei = range(num_hei)
498 499 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
499 500 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
500 501 mask_prof = numpy.asarray(range(num_prof))
501 502 num_mask_prof = mask_prof.size
502 503 comp_mask_prof = [0, num_prof/2]
503 504
504 505
505 506 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
506 507 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
507 508 jnoise = numpy.nan
508 509 noise_exist = jnoise[0] < numpy.Inf
509 510
510 511 #Subrutina de Remocion de la Interferencia
511 512 for ich in range(num_channel):
512 513 #Se ordena los espectros segun su potencia (menor a mayor)
513 514 power = jspectra[ich,mask_prof,:]
514 515 power = power[:,hei_interf]
515 516 power = power.sum(axis = 0)
516 517 psort = power.ravel().argsort()
517 518
518 519 #Se estima la interferencia promedio en los Espectros de Potencia empleando
519 520 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
520 521
521 522 if noise_exist:
522 523 # tmp_noise = jnoise[ich] / num_prof
523 524 tmp_noise = jnoise[ich]
524 525 junkspc_interf = junkspc_interf - tmp_noise
525 526 #junkspc_interf[:,comp_mask_prof] = 0
526 527
527 528 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
528 529 jspc_interf = jspc_interf.transpose()
529 530 #Calculando el espectro de interferencia promedio
530 531 noiseid = numpy.where(jspc_interf <= tmp_noise/ numpy.sqrt(num_incoh))
531 532 noiseid = noiseid[0]
532 533 cnoiseid = noiseid.size
533 534 interfid = numpy.where(jspc_interf > tmp_noise/ numpy.sqrt(num_incoh))
534 535 interfid = interfid[0]
535 536 cinterfid = interfid.size
536 537
537 538 if (cnoiseid > 0): jspc_interf[noiseid] = 0
538 539
539 540 #Expandiendo los perfiles a limpiar
540 541 if (cinterfid > 0):
541 542 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
542 543 new_interfid = numpy.asarray(new_interfid)
543 544 new_interfid = {x for x in new_interfid}
544 545 new_interfid = numpy.array(list(new_interfid))
545 546 new_cinterfid = new_interfid.size
546 547 else: new_cinterfid = 0
547 548
548 549 for ip in range(new_cinterfid):
549 550 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
550 551 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
551 552
552 553
553 554 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
554 555
555 556 #Removiendo la interferencia del punto de mayor interferencia
556 557 ListAux = jspc_interf[mask_prof].tolist()
557 558 maxid = ListAux.index(max(ListAux))
558 559
559 560
560 561 if cinterfid > 0:
561 562 for ip in range(cinterfid*(interf == 2) - 1):
562 563 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/numpy.sqrt(num_incoh))).nonzero()
563 564 cind = len(ind)
564 565
565 566 if (cind > 0):
566 567 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/numpy.sqrt(num_incoh))
567 568
568 569 ind = numpy.array([-2,-1,1,2])
569 570 xx = numpy.zeros([4,4])
570 571
571 572 for id1 in range(4):
572 573 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
573 574
574 575 xx_inv = numpy.linalg.inv(xx)
575 576 xx = xx_inv[:,0]
576 577 ind = (ind + maxid + num_mask_prof)%num_mask_prof
577 578 yy = jspectra[ich,mask_prof[ind],:]
578 579 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
579 580
580 581
581 582 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/numpy.sqrt(num_incoh))).nonzero()
582 583 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/numpy.sqrt(num_incoh))
583 584
584 585 #Remocion de Interferencia en el Cross Spectra
585 586 if jcspectra is None: return jspectra, jcspectra
586 587 num_pairs = jcspectra.size/(num_prof*num_hei)
587 588 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
588 589
589 590 for ip in range(num_pairs):
590 591
591 592 #-------------------------------------------
592 593
593 594 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
594 595 cspower = cspower[:,hei_interf]
595 596 cspower = cspower.sum(axis = 0)
596 597
597 598 cspsort = cspower.ravel().argsort()
598 599 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
599 600 junkcspc_interf = junkcspc_interf.transpose()
600 601 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
601 602
602 603 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
603 604
604 605 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
605 606 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
606 607 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
607 608
608 609 for iprof in range(num_prof):
609 610 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
610 611 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
611 612
612 613 #Removiendo la Interferencia
613 614 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
614 615
615 616 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
616 617 maxid = ListAux.index(max(ListAux))
617 618
618 619 ind = numpy.array([-2,-1,1,2])
619 620 xx = numpy.zeros([4,4])
620 621
621 622 for id1 in range(4):
622 623 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
623 624
624 625 xx_inv = numpy.linalg.inv(xx)
625 626 xx = xx_inv[:,0]
626 627
627 628 ind = (ind + maxid + num_mask_prof)%num_mask_prof
628 629 yy = jcspectra[ip,mask_prof[ind],:]
629 630 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
630 631
631 632 #Guardar Resultados
632 633 self.dataOut.data_spc = jspectra
633 634 self.dataOut.data_cspc = jcspectra
634 635
635 636 return 1
636 637
637 638 def setRadarFrequency(self, frequency=None):
638 639
639 640 if frequency != None:
640 641 self.dataOut.frequency = frequency
641 642
642 643 return 1
643 644
644 645 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
645 646 #validacion de rango
646 647 if minHei == None:
647 648 minHei = self.dataOut.heightList[0]
648 649
649 650 if maxHei == None:
650 651 maxHei = self.dataOut.heightList[-1]
651 652
652 653 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
653 654 print 'minHei: %.2f is out of the heights range'%(minHei)
654 655 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
655 656 minHei = self.dataOut.heightList[0]
656 657
657 658 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
658 659 print 'maxHei: %.2f is out of the heights range'%(maxHei)
659 660 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
660 661 maxHei = self.dataOut.heightList[-1]
661 662
662 663 # validacion de velocidades
663 664 velrange = self.dataOut.getVelRange(1)
664 665
665 666 if minVel == None:
666 667 minVel = velrange[0]
667 668
668 669 if maxVel == None:
669 670 maxVel = velrange[-1]
670 671
671 672 if (minVel < velrange[0]) or (minVel > maxVel):
672 673 print 'minVel: %.2f is out of the velocity range'%(minVel)
673 674 print 'minVel is setting to %.2f'%(velrange[0])
674 675 minVel = velrange[0]
675 676
676 677 if (maxVel > velrange[-1]) or (maxVel < minVel):
677 678 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
678 679 print 'maxVel is setting to %.2f'%(velrange[-1])
679 680 maxVel = velrange[-1]
680 681
681 682 # seleccion de indices para rango
682 683 minIndex = 0
683 684 maxIndex = 0
684 685 heights = self.dataOut.heightList
685 686
686 687 inda = numpy.where(heights >= minHei)
687 688 indb = numpy.where(heights <= maxHei)
688 689
689 690 try:
690 691 minIndex = inda[0][0]
691 692 except:
692 693 minIndex = 0
693 694
694 695 try:
695 696 maxIndex = indb[0][-1]
696 697 except:
697 698 maxIndex = len(heights)
698 699
699 700 if (minIndex < 0) or (minIndex > maxIndex):
700 701 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
701 702
702 703 if (maxIndex >= self.dataOut.nHeights):
703 704 maxIndex = self.dataOut.nHeights-1
704 705
705 706 # seleccion de indices para velocidades
706 707 indminvel = numpy.where(velrange >= minVel)
707 708 indmaxvel = numpy.where(velrange <= maxVel)
708 709 try:
709 710 minIndexVel = indminvel[0][0]
710 711 except:
711 712 minIndexVel = 0
712 713
713 714 try:
714 715 maxIndexVel = indmaxvel[0][-1]
715 716 except:
716 717 maxIndexVel = len(velrange)
717 718
718 719 #seleccion del espectro
719 720 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
720 721 #estimacion de ruido
721 722 noise = numpy.zeros(self.dataOut.nChannels)
722 723
723 724 for channel in range(self.dataOut.nChannels):
724 725 daux = data_spc[channel,:,:]
725 726 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
726 727
727 728 self.dataOut.noise_estimation = noise.copy()
728 729
729 730 return 1
730 731
731 732 class IncohIntLags(Operation):
732 733
733 734
734 735 __profIndex = 0
735 736 __withOverapping = False
736 737
737 738 __byTime = False
738 739 __initime = None
739 740 __lastdatatime = None
740 741 __integrationtime = None
741 742
742 743 __buffer_spc = None
743 744 __buffer_cspc = None
744 745 __buffer_dc = None
745 746
746 747 __dataReady = False
747 748
748 749 __timeInterval = None
749 750
750 751 n = None
751 752
752 753
753 754
754 755 def __init__(self):
755 756
756 757 Operation.__init__(self)
757 758 # self.isConfig = False
758 759
759 760 def setup(self, n=None, timeInterval=None, overlapping=False):
760 761 """
761 762 Set the parameters of the integration class.
762 763
763 764 Inputs:
764 765
765 766 n : Number of coherent integrations
766 767 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
767 768 overlapping :
768 769
769 770 """
770 771
771 772 self.__initime = None
772 773 self.__lastdatatime = 0
773 774
774 775 self.__buffer_spc = 0
775 776 self.__buffer_cspc = 0
776 777 self.__buffer_dc = 0
777 778
778 779 self.__profIndex = 0
779 780 self.__dataReady = False
780 781 self.__byTime = False
781 782
782 783 if n is None and timeInterval is None:
783 784 raise ValueError, "n or timeInterval should be specified ..."
784 785
785 786 if n is not None:
786 787 self.n = int(n)
787 788 else:
788 789 self.__integrationtime = int(timeInterval) #if (type(timeInterval)!=integer) -> change this line
789 790 self.n = None
790 791 self.__byTime = True
791 792
792 793 def putData(self, data_spc, data_cspc, data_dc):
793 794
794 795 """
795 796 Add a profile to the __buffer_spc and increase in one the __profileIndex
796 797
797 798 """
798 799
799 800 self.__buffer_spc += data_spc
800 801
801 802 if data_cspc is None:
802 803 self.__buffer_cspc = None
803 804 else:
804 805 self.__buffer_cspc += data_cspc
805 806
806 807 if data_dc is None:
807 808 self.__buffer_dc = None
808 809 else:
809 810 self.__buffer_dc += data_dc
810 811
811 812 self.__profIndex += 1
812 813
813 814 return
814 815
815 816 def pushData(self):
816 817 """
817 818 Return the sum of the last profiles and the profiles used in the sum.
818 819
819 820 Affected:
820 821
821 822 self.__profileIndex
822 823
823 824 """
824 825
825 826 data_spc = self.__buffer_spc
826 827 data_cspc = self.__buffer_cspc
827 828 data_dc = self.__buffer_dc
828 829 n = self.__profIndex
829 830
830 831 self.__buffer_spc = 0
831 832 self.__buffer_cspc = 0
832 833 self.__buffer_dc = 0
833 834 self.__profIndex = 0
834 835
835 836 return data_spc, data_cspc, data_dc, n
836 837
837 838 def byProfiles(self, *args):
838 839
839 840 self.__dataReady = False
840 841 avgdata_spc = None
841 842 avgdata_cspc = None
842 843 avgdata_dc = None
843 844
844 845 self.putData(*args)
845 846
846 847 if self.__profIndex == self.n:
847 848
848 849 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
849 850 self.n = n
850 851 self.__dataReady = True
851 852
852 853 return avgdata_spc, avgdata_cspc, avgdata_dc
853 854
854 855 def byTime(self, datatime, *args):
855 856
856 857 self.__dataReady = False
857 858 avgdata_spc = None
858 859 avgdata_cspc = None
859 860 avgdata_dc = None
860 861
861 862 self.putData(*args)
862 863
863 864 if (datatime - self.__initime) >= self.__integrationtime:
864 865 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
865 866 self.n = n
866 867 self.__dataReady = True
867 868
868 869 return avgdata_spc, avgdata_cspc, avgdata_dc
869 870
870 871 def integrate(self, datatime, *args):
871 872
872 873 if self.__profIndex == 0:
873 874 self.__initime = datatime
874 875
875 876 if self.__byTime:
876 877 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
877 878 else:
878 879 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
879 880
880 881 if not self.__dataReady:
881 882 return None, None, None, None
882 883
883 884 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
884 885
885 886 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
886 887
887 888 if n==1:
888 889 return
889 890
890 891 dataOut.flagNoData = True
891 892
892 893 if not self.isConfig:
893 894 self.setup(n, timeInterval, overlapping)
894 895 self.isConfig = True
895 896
896 897 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
897 898 dataOut.data_spc,
898 899 dataOut.data_cspc,
899 900 dataOut.data_dc)
900 901
901 902 if self.__dataReady:
902 903
903 904 dataOut.data_spc = avgdata_spc
904 905 dataOut.data_cspc = avgdata_cspc
905 906 dataOut.data_dc = avgdata_dc
906 907
907 908 dataOut.nIncohInt *= self.n
908 909 dataOut.utctime = avgdatatime
909 910 dataOut.flagNoData = False
General Comments 0
You need to be logged in to leave comments. Login now