schain Commit - r794:e2c8bf5aa032 · Jicamarca Repository

spectra_lags and spectra_acf added (Kudeki)

Miguel Valdez -

r794:e2c8bf5aa032

parent child

Context file:

r794:e2c8bf5aa032

Collapse all files

schainpy/model/proc/jroproc_spectra_acf.py created 644 +736 0

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

schainpy/model/proc/__init__.py +2 -1

             '''
             $Author: murco $
             $Id: Processor.py 1 2012-11-12 18:56:07Z murco $
             '''
             from jroproc_voltage import *
             from jroproc_spectra import *
             from jroproc_heispectra import *
             from jroproc_amisr import *
             from jroproc_correlation import *
             from jroproc_parameters import *
-            from jroproc_spectra_lags import *
  No newline at end of file
+            from jroproc_spectra_lags import *
+            from jroproc_spectra_acf import *
  No newline at end of file

schainpy/model/proc/jroproc_spectra_lags.py +4 -4

             import numpy
             from jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jrodata import Spectra
             from schainpy.model.data.jrodata import hildebrand_sekhon
             class SpectraLagsProc(ProcessingUnit):
                 def __init__(self):
                     ProcessingUnit.__init__(self)
                     self.__input_buffer = None
                     self.firstdatatime = None
                     self.profIndex = 0
                     self.dataOut = Spectra()
                     self.id_min = None
                     self.id_max = None
                     self.__codeIndex = 0
                     self.__lags_buffer = None
                 def __updateSpecFromVoltage(self):
                     self.dataOut.plotting = "spectra_lags"
                     self.dataOut.timeZone = self.dataIn.timeZone
                     self.dataOut.dstFlag = self.dataIn.dstFlag
                     self.dataOut.errorCount = self.dataIn.errorCount
                     self.dataOut.useLocalTime = self.dataIn.useLocalTime
                     self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
                     self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
                     self.dataOut.ippSeconds = self.dataIn.getDeltaH()*(10**-6)/0.15
                     self.dataOut.channelList = self.dataIn.channelList
                     self.dataOut.heightList = self.dataIn.heightList
                     self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
                     self.dataOut.nBaud = self.dataIn.nBaud
                     self.dataOut.nCode = self.dataIn.nCode
                     self.dataOut.code = self.dataIn.code
             #         self.dataOut.nProfiles = self.dataOut.nFFTPoints
                     self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
                     self.dataOut.utctime = self.firstdatatime
                     self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
                     self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
                     self.dataOut.flagShiftFFT = False
                     self.dataOut.nCohInt = self.dataIn.nCohInt
                     self.dataOut.nIncohInt = 1
                     self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
                     self.dataOut.frequency = self.dataIn.frequency
                     self.dataOut.realtime = self.dataIn.realtime
                     self.dataOut.azimuth = self.dataIn.azimuth
                     self.dataOut.zenith = self.dataIn.zenith
                     self.dataOut.beam.codeList = self.dataIn.beam.codeList
                     self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
                     self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
                 def __createLagsBlock(self, voltages):
                     if self.__lags_buffer is None:
                         self.__lags_buffer = numpy.zeros((self.dataOut.nChannels, self.dataOut.nProfiles, self.dataOut.nHeights), dtype='complex')
                     nsegments = self.dataOut.nHeights - self.dataOut.nProfiles
             #         codes = numpy.conjugate(self.__input_buffer[:,9:169])/10000
                     for i in range(nsegments):
                         self.__lags_buffer[:,:,i] = voltages[:,i:i+self.dataOut.nProfiles]#*codes
                     return self.__lags_buffer
                 def __decodeData(self, volt_buffer, pulseIndex=None):
                     if pulseIndex is None:
                         return volt_buffer
                     codes = numpy.conjugate(self.__input_buffer[:,pulseIndex[0]:pulseIndex[1]])/10000
                     nsegments = self.dataOut.nHeights - self.dataOut.nProfiles
                     for i in range(nsegments):
                         volt_buffer[:,:,i] = volt_buffer[:,:,i]*codes
                     return volt_buffer
                 def __getFft(self, datablock):
                     """
                     Convierte valores de Voltaje a Spectra
                     Affected:
                         self.dataOut.data_spc
                         self.dataOut.data_cspc
                         self.dataOut.data_dc
                         self.dataOut.heightList
                         self.profIndex
                         self.__input_buffer
                         self.dataOut.flagNoData
                     """
                     fft_volt = numpy.fft.fft(datablock, n=self.dataOut.nFFTPoints, axis=1)
-            #         dc = fft_volt[:,0,:]
+                    dc = fft_volt[:,0,:]
                     #calculo de self-spectra
                     fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
                     spc = fft_volt * numpy.conjugate(fft_volt)
                     spc = spc.real
                     blocksize = 0
-            #         blocksize += dc.size
+                    blocksize += dc.size
                     blocksize += spc.size
                     cspc = None
                     pairIndex = 0
-                    if self.dataOut.pairsList != None:
+                    if self.dataOut.pairsList != []:
                         #calculo de cross-spectra
                         cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
                         for pair in self.dataOut.pairsList:
                             if pair[0] not in self.dataOut.channelList:
                                 raise ValueError, "Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
                             if pair[1] not in self.dataOut.channelList:
                                 raise ValueError, "Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList))
                             chan_index0 = self.dataOut.channelList.index(pair[0])
                             chan_index1 = self.dataOut.channelList.index(pair[1])
                             cspc[pairIndex,:,:] = fft_volt[chan_index0,:,:] * numpy.conjugate(fft_volt[chan_index1,:,:])
                             pairIndex += 1
                         blocksize += cspc.size
                     self.dataOut.data_spc = spc
                     self.dataOut.data_cspc = cspc
-            #         self.dataOut.data_dc = dc
+                    self.dataOut.data_dc = dc
                     self.dataOut.blockSize = blocksize
                     self.dataOut.flagShiftFFT = True
                 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], code=None, nCode=None, nBaud=None, codeFromHeader=False, pulseIndex=None):
                     self.dataOut.flagNoData = True
                     self.code = None
                     if codeFromHeader:
                         if self.dataIn.code is not None:
                             self.code = self.dataIn.code
                     if code is not None:
                         self.code = numpy.array(code).reshape(nCode,nBaud)
                     if self.dataIn.type == "Voltage":
                         if nFFTPoints == None:
                             raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
                         if nProfiles == None:
                             nProfiles = nFFTPoints
                         self.profIndex == nProfiles
                         self.firstdatatime = self.dataIn.utctime
                         self.dataOut.ippFactor = 1
                         self.dataOut.nFFTPoints = nFFTPoints
                         self.dataOut.nProfiles = nProfiles
                         self.dataOut.pairsList = pairsList
                         self.__updateSpecFromVoltage()
                         if not self.dataIn.flagDataAsBlock:
                             self.__input_buffer = self.dataIn.data.copy()
                             lags_block = self.__createLagsBlock(self.__input_buffer)
                             lags_block = self.__decodeData(lags_block, pulseIndex)
                         else:
                             self.__input_buffer = self.dataIn.data.copy()
                         self.__getFft(lags_block)
                         self.dataOut.flagNoData = False
                         return True
                     raise ValueError, "The type of input object '%s' is not valid"%(self.dataIn.type)
                 def __selectPairs(self, pairsList):
                     if channelList == None:
                         return
                     pairsIndexListSelected = []
                     for thisPair in pairsList:
                         if thisPair not in self.dataOut.pairsList:
                             continue
                         pairIndex = self.dataOut.pairsList.index(thisPair)
                         pairsIndexListSelected.append(pairIndex)
                     if not pairsIndexListSelected:
                         self.dataOut.data_cspc = None
                         self.dataOut.pairsList = []
                         return
                     self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
                     self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
                     return
                 def __selectPairsByChannel(self, channelList=None):
                     if channelList == None:
                         return
                     pairsIndexListSelected = []
                     for pairIndex in self.dataOut.pairsIndexList:
                         #First pair
                         if self.dataOut.pairsList[pairIndex][0] not in channelList:
                             continue
                         #Second pair
                         if self.dataOut.pairsList[pairIndex][1] not in channelList:
                             continue
                         pairsIndexListSelected.append(pairIndex)
                     if not pairsIndexListSelected:
                         self.dataOut.data_cspc = None
                         self.dataOut.pairsList = []
                         return
                     self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
                     self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
                     return
                 def selectChannels(self, channelList):
                     channelIndexList = []
                     for channel in channelList:
                         if channel not in self.dataOut.channelList:
                             raise ValueError, "Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList))
                         index = self.dataOut.channelList.index(channel)
                         channelIndexList.append(index)
                     self.selectChannelsByIndex(channelIndexList)
                 def selectChannelsByIndex(self, channelIndexList):
                     """
                     Selecciona un bloque de datos en base a canales segun el channelIndexList
                     Input:
                         channelIndexList    :    lista sencilla de canales a seleccionar por ej. [2,3,7]
                     Affected:
                         self.dataOut.data_spc
                         self.dataOut.channelIndexList
                         self.dataOut.nChannels
                     Return:
                         None
                     """
                     for channelIndex in channelIndexList:
                         if channelIndex not in self.dataOut.channelIndexList:
                             raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList)
             #         nChannels = len(channelIndexList)
                     data_spc = self.dataOut.data_spc[channelIndexList,:]
                     data_dc = self.dataOut.data_dc[channelIndexList,:]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.data_dc = data_dc
                     self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
             #        self.dataOut.nChannels = nChannels
                     self.__selectPairsByChannel(self.dataOut.channelList)
                     return 1
                 def selectHeights(self, minHei, maxHei):
                     """
                     Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
                     minHei <= height <= maxHei
                     Input:
                         minHei    :    valor minimo de altura a considerar
                         maxHei    :    valor maximo de altura a considerar
                     Affected:
                         Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
                     Return:
 si el metodo se ejecuto con exito caso contrario devuelve 0
                     """
                     if (minHei > maxHei):
                         raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
                     if (minHei < self.dataOut.heightList[0]):
                         minHei = self.dataOut.heightList[0]
                     if (maxHei > self.dataOut.heightList[-1]):
                         maxHei = self.dataOut.heightList[-1]
                     minIndex = 0
                     maxIndex = 0
                     heights = self.dataOut.heightList
                     inda = numpy.where(heights >= minHei)
                     indb = numpy.where(heights <= maxHei)
                     try:
                         minIndex = inda[0][0]
                     except:
                         minIndex = 0
                     try:
                         maxIndex = indb[0][-1]
                     except:
                         maxIndex = len(heights)
                     self.selectHeightsByIndex(minIndex, maxIndex)
                     return 1
                 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
                     newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
                     if hei_ref != None:
                         newheis = numpy.where(self.dataOut.heightList>hei_ref)
                     minIndex = min(newheis[0])
                     maxIndex = max(newheis[0])
                     data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
                     heightList = self.dataOut.heightList[minIndex:maxIndex+1]
                     # determina indices
                     nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
                     avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
                     beacon_dB = numpy.sort(avg_dB)[-nheis:]
                     beacon_heiIndexList = []
                     for val in avg_dB.tolist():
                         if val >= beacon_dB[0]:
                             beacon_heiIndexList.append(avg_dB.tolist().index(val))
                     #data_spc = data_spc[:,:,beacon_heiIndexList]
                     data_cspc = None
                     if self.dataOut.data_cspc is not None:
                         data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
                         #data_cspc = data_cspc[:,:,beacon_heiIndexList]
                     data_dc = None
                     if self.dataOut.data_dc is not None:
                         data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
                         #data_dc = data_dc[:,beacon_heiIndexList]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.data_cspc = data_cspc
                     self.dataOut.data_dc = data_dc
                     self.dataOut.heightList = heightList
                     self.dataOut.beacon_heiIndexList = beacon_heiIndexList
                     return 1
                 def selectHeightsByIndex(self, minIndex, maxIndex):
                     """
                     Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
                     minIndex <= index <= maxIndex
                     Input:
                         minIndex    :    valor de indice minimo de altura a considerar
                         maxIndex    :    valor de indice maximo de altura a considerar
                     Affected:
                         self.dataOut.data_spc
                         self.dataOut.data_cspc
                         self.dataOut.data_dc
                         self.dataOut.heightList
                     Return:
 si el metodo se ejecuto con exito caso contrario devuelve 0
                     """
                     if (minIndex < 0) or (minIndex > maxIndex):
                         raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex)
                     if (maxIndex >= self.dataOut.nHeights):
                         maxIndex = self.dataOut.nHeights-1
                     #Spectra
                     data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
                     data_cspc = None
                     if self.dataOut.data_cspc is not None:
                         data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
                     data_dc = None
                     if self.dataOut.data_dc is not None:
                         data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.data_cspc = data_cspc
                     self.dataOut.data_dc = data_dc
                     self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
                     return 1
                 def removeDC(self, mode = 2):
                     jspectra = self.dataOut.data_spc
                     jcspectra = self.dataOut.data_cspc
                     num_chan = jspectra.shape[0]
                     num_hei = jspectra.shape[2]
                     if jcspectra is not None:
                         jcspectraExist = True
                         num_pairs = jcspectra.shape[0]
                     else:   jcspectraExist = False
                     freq_dc = jspectra.shape[1]/2
                     ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
                     if ind_vel[0]<0:
                         ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
                     if mode == 1:
                         jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
                         if jcspectraExist:
                             jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
                     if mode == 2:
                         vel = numpy.array([-2,-1,1,2])
                         xx = numpy.zeros([4,4])
                         for fil in range(4):
                             xx[fil,:] = vel[fil]**numpy.asarray(range(4))
                         xx_inv = numpy.linalg.inv(xx)
                         xx_aux = xx_inv[0,:]
                         for ich in range(num_chan):
                             yy = jspectra[ich,ind_vel,:]
                             jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
                             junkid = jspectra[ich,freq_dc,:]<=0
                             cjunkid = sum(junkid)
                             if cjunkid.any():
                                 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
                         if jcspectraExist:
                             for ip in range(num_pairs):
                                 yy = jcspectra[ip,ind_vel,:]
                                 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
                     self.dataOut.data_spc = jspectra
                     self.dataOut.data_cspc = jcspectra
                     return 1
                 def removeInterference(self,  interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
                     jspectra = self.dataOut.data_spc
                     jcspectra = self.dataOut.data_cspc
                     jnoise = self.dataOut.getNoise()
                     num_incoh = self.dataOut.nIncohInt
                     num_channel  = jspectra.shape[0]
                     num_prof  = jspectra.shape[1]
                     num_hei   = jspectra.shape[2]
                     #hei_interf
                     if hei_interf is None:
                         count_hei = num_hei/2   #Como es entero no importa
                         hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
                         hei_interf = numpy.asarray(hei_interf)[0]
                     #nhei_interf
                     if (nhei_interf == None):
                         nhei_interf = 5
                     if (nhei_interf < 1):
                         nhei_interf = 1
                     if (nhei_interf > count_hei):
                         nhei_interf = count_hei
                     if (offhei_interf == None):
                         offhei_interf = 0
                     ind_hei = range(num_hei)
             #         mask_prof = numpy.asarray(range(num_prof - 2)) + 1
             #         mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
                     mask_prof = numpy.asarray(range(num_prof))
                     num_mask_prof = mask_prof.size
                     comp_mask_prof = [0, num_prof/2]
                     #noise_exist:    Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
                     if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
                         jnoise = numpy.nan
                     noise_exist = jnoise[0] < numpy.Inf
                     #Subrutina de Remocion de la Interferencia
                     for ich in range(num_channel):
                         #Se ordena los espectros segun su potencia (menor a mayor)
                         power = jspectra[ich,mask_prof,:]
                         power = power[:,hei_interf]
                         power = power.sum(axis = 0)
                         psort = power.ravel().argsort()
                         #Se estima la interferencia promedio en los Espectros de Potencia empleando
                         junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
                         if noise_exist:
                         #    tmp_noise = jnoise[ich] / num_prof
                             tmp_noise = jnoise[ich]
                         junkspc_interf = junkspc_interf - tmp_noise
                         #junkspc_interf[:,comp_mask_prof] = 0
                         jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
                         jspc_interf = jspc_interf.transpose()
                         #Calculando el espectro de interferencia promedio
                         noiseid =  numpy.where(jspc_interf <= tmp_noise/ numpy.sqrt(num_incoh))
                         noiseid = noiseid[0]
                         cnoiseid = noiseid.size
                         interfid = numpy.where(jspc_interf > tmp_noise/ numpy.sqrt(num_incoh))
                         interfid = interfid[0]
                         cinterfid = interfid.size
                         if (cnoiseid > 0):   jspc_interf[noiseid] = 0
                         #Expandiendo los perfiles a limpiar
                         if (cinterfid > 0):
                             new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
                             new_interfid = numpy.asarray(new_interfid)
                             new_interfid = {x for x in new_interfid}
                             new_interfid = numpy.array(list(new_interfid))
                             new_cinterfid = new_interfid.size
                         else: new_cinterfid = 0
                         for ip in range(new_cinterfid):
                             ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
                             jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
                         jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
                         #Removiendo la interferencia del punto de mayor interferencia
                         ListAux = jspc_interf[mask_prof].tolist()
                         maxid = ListAux.index(max(ListAux))
                         if cinterfid > 0:
                             for ip in range(cinterfid*(interf == 2) - 1):
                                 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/numpy.sqrt(num_incoh))).nonzero()
                                 cind = len(ind)
                                 if (cind > 0):
                                     jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/numpy.sqrt(num_incoh))
                             ind = numpy.array([-2,-1,1,2])
                             xx = numpy.zeros([4,4])
                             for id1 in range(4):
                                 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
                             xx_inv = numpy.linalg.inv(xx)
                             xx = xx_inv[:,0]
                             ind = (ind + maxid + num_mask_prof)%num_mask_prof
                             yy = jspectra[ich,mask_prof[ind],:]
                             jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
                         indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/numpy.sqrt(num_incoh))).nonzero()
                         jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/numpy.sqrt(num_incoh))
                     #Remocion de Interferencia en el Cross Spectra
                     if jcspectra is None: return jspectra, jcspectra
                     num_pairs = jcspectra.size/(num_prof*num_hei)
                     jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
                     for ip in range(num_pairs):
                         #-------------------------------------------
                         cspower = numpy.abs(jcspectra[ip,mask_prof,:])
                         cspower = cspower[:,hei_interf]
                         cspower = cspower.sum(axis = 0)
                         cspsort = cspower.ravel().argsort()
                         junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
                         junkcspc_interf = junkcspc_interf.transpose()
                         jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
                         ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
                         median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
                         median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
                         junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
                         for iprof in range(num_prof):
                             ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
                             jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
                         #Removiendo la Interferencia
                         jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
                         ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
                         maxid = ListAux.index(max(ListAux))
                         ind = numpy.array([-2,-1,1,2])
                         xx = numpy.zeros([4,4])
                         for id1 in range(4):
                             xx[:,id1] = ind[id1]**numpy.asarray(range(4))
                         xx_inv = numpy.linalg.inv(xx)
                         xx = xx_inv[:,0]
                         ind = (ind + maxid + num_mask_prof)%num_mask_prof
                         yy = jcspectra[ip,mask_prof[ind],:]
                         jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
                     #Guardar Resultados
                     self.dataOut.data_spc = jspectra
                     self.dataOut.data_cspc = jcspectra
                     return 1
                 def setRadarFrequency(self, frequency=None):
                     if frequency != None:
                         self.dataOut.frequency = frequency
                     return 1
                 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
                     #validacion de rango
                     if minHei == None:
                         minHei = self.dataOut.heightList[0]
                     if maxHei == None:
                         maxHei = self.dataOut.heightList[-1]
                     if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
                         print 'minHei: %.2f is out of the heights range'%(minHei)
                         print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
                         minHei = self.dataOut.heightList[0]
                     if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
                         print 'maxHei: %.2f is out of the heights range'%(maxHei)
                         print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
                         maxHei = self.dataOut.heightList[-1]
                     # validacion de velocidades
                     velrange = self.dataOut.getVelRange(1)
                     if minVel == None:
                         minVel = velrange[0]
                     if maxVel == None:
                         maxVel = velrange[-1]
                     if (minVel < velrange[0]) or (minVel > maxVel):
                         print 'minVel: %.2f is out of the velocity range'%(minVel)
                         print 'minVel is setting to %.2f'%(velrange[0])
                         minVel = velrange[0]
                     if (maxVel > velrange[-1]) or (maxVel < minVel):
                         print 'maxVel: %.2f is out of the velocity range'%(maxVel)
                         print 'maxVel is setting to %.2f'%(velrange[-1])
                         maxVel = velrange[-1]
                     # seleccion de indices para rango
                     minIndex = 0
                     maxIndex = 0
                     heights = self.dataOut.heightList
                     inda = numpy.where(heights >= minHei)
                     indb = numpy.where(heights <= maxHei)
                     try:
                         minIndex = inda[0][0]
                     except:
                         minIndex = 0
                     try:
                         maxIndex = indb[0][-1]
                     except:
                         maxIndex = len(heights)
                     if (minIndex < 0) or (minIndex > maxIndex):
                         raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
                     if (maxIndex >= self.dataOut.nHeights):
                         maxIndex = self.dataOut.nHeights-1
                     # seleccion de indices para velocidades
                     indminvel = numpy.where(velrange >= minVel)
                     indmaxvel = numpy.where(velrange <= maxVel)
                     try:
                         minIndexVel = indminvel[0][0]
                     except:
                         minIndexVel = 0
                     try:
                         maxIndexVel = indmaxvel[0][-1]
                     except:
                         maxIndexVel = len(velrange)
                     #seleccion del espectro
                     data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
                     #estimacion de ruido
                     noise = numpy.zeros(self.dataOut.nChannels)
                     for channel in range(self.dataOut.nChannels):
                         daux = data_spc[channel,:,:]
                         noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
                     self.dataOut.noise_estimation = noise.copy()
                     return 1

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages