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',' 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: 1 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