diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index 4f4748a..a54b500 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -2,7 +2,7 @@ import sys import numpy,math from scipy import interpolate from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator -from schainpy.model.data.jrodata import Voltage +from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon from schainpy.utils import log from time import time @@ -1355,9 +1355,6 @@ class PulsePairVoltage(Operation): n, dataOut.nHeights), dtype='complex') - #self.noise = numpy.zeros([self.__nch,self.__nHeis]) - #for i in range(self.__nch): - # self.noise[i]=dataOut.getNoise(channel=i) def putData(self,data): ''' @@ -1372,101 +1369,127 @@ class PulsePairVoltage(Operation): Return the PULSEPAIR and the profiles used in the operation Affected : self.__profileIndex ''' + #················· Remove DC··································· if self.removeDC==True: mean = numpy.mean(self.__buffer,1) tmp = mean.reshape(self.__nch,1,self.__nHeis) dc= numpy.tile(tmp,[1,self.__nProf,1]) self.__buffer = self.__buffer - dc + #··················Calculo de Potencia ························ + pair0 = self.__buffer*numpy.conj(self.__buffer) + pair0 = pair0.real + lag_0 = numpy.sum(pair0,1) + #··················Calculo de Ruido x canal···················· + self.noise = numpy.zeros(self.__nch) + for i in range(self.__nch): + daux = numpy.sort(pair0[i,:,:],axis= None) + self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt) + + self.noise = self.noise.reshape(self.__nch,1) + self.noise = numpy.tile(self.noise,[1,self.__nHeis]) + noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis) + noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1]) + #·················· Potencia recibida= P , Potencia senal = S , Ruido= N·· + #·················· P= S+N ,P=lag_0/N ································· + #···················· Power ·················································· + data_power = lag_0/(self.n*self.nCohInt) + #------------------ Senal ··················································· + data_intensity = pair0 - noise_buffer + data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt) + #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt) + for i in range(self.__nch): + for j in range(self.__nHeis): + if data_intensity[i][j] < 0: + data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j])) - lag_0 = numpy.sum(self.__buffer*numpy.conj(self.__buffer),1) - data_intensity = lag_0/(self.n*self.nCohInt)#*self.nCohInt) - + #················· Calculo de Frecuencia y Velocidad doppler········ pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:]) lag_1 = numpy.sum(pair1,1) - #angle = numpy.angle(numpy.sum(pair1,1))*180/(math.pi) - data_velocity = (-1.0*self.lambda_/(4*math.pi*self.ippSec))*numpy.angle(lag_1)#self.ippSec*self.nCohInt + data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1) + data_velocity = (self.lambda_/2.0)*data_freq - self.noise = numpy.zeros([self.__nch,self.__nHeis]) - for i in range(self.__nch): - self.noise[i]=dataOut.getNoise(channel=i) + #················ Potencia promedio estimada de la Senal··········· + lag_0 = lag_0/self.n + S = lag_0-self.noise - lag_0 = lag_0.real/(self.n) + #················ Frecuencia Doppler promedio ····················· lag_1 = lag_1/(self.n-1) R1 = numpy.abs(lag_1) - S = (lag_0-self.noise) + #················ Calculo del SNR·································· data_snrPP = S/self.noise - data_snrPP = numpy.where(data_snrPP<0,1,data_snrPP) + for i in range(self.__nch): + for j in range(self.__nHeis): + if data_snrPP[i][j] < 1.e-20: + data_snrPP[i][j] = 1.e-20 + #················· Calculo del ancho espectral ······················ L = S/R1 L = numpy.where(L<0,1,L) L = numpy.log(L) - tmp = numpy.sqrt(numpy.absolute(L)) - - data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec))*tmp*numpy.sign(L) - #data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec))*k + data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L) n = self.__profIndex self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex') self.__profIndex = 0 - return data_intensity,data_velocity,data_snrPP,data_specwidth,n + return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n + def pulsePairbyProfiles(self,dataOut): self.__dataReady = False + data_power = None data_intensity = None data_velocity = None data_specwidth = None data_snrPP = None self.putData(data=dataOut.data) if self.__profIndex == self.n: - #self.noise = numpy.zeros([self.__nch,self.__nHeis]) - #for i in range(self.__nch): - # self.noise[i]=data.getNoise(channel=i) - #print(self.noise.shape) - data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut) + data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut) self.__dataReady = True - return data_intensity, data_velocity,data_snrPP,data_specwidth + return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth + def pulsePairOp(self, dataOut, datatime= None): if self.__initime == None: self.__initime = datatime - #print("hola") - data_intensity, data_velocity,data_snrPP,data_specwidth = self.pulsePairbyProfiles(dataOut) + data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut) self.__lastdatatime = datatime - if data_intensity is None: - return None, None,None,None,None + if data_power is None: + return None, None, None,None,None,None avgdatatime = self.__initime deltatime = datatime - self.__lastdatatime self.__initime = datatime - return data_intensity, data_velocity,data_snrPP,data_specwidth,avgdatatime + return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs): if not self.isConfig: self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs) self.isConfig = True - data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime) + data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime) dataOut.flagNoData = True if self.__dataReady: dataOut.nCohInt *= self.n - dataOut.data_intensity = data_intensity #valor para intensidad - dataOut.data_velocity = data_velocity #valor para velocidad - dataOut.data_snrPP = data_snrPP # valor para snr - dataOut.data_specwidth = data_specwidth + dataOut.dataPP_POW = data_intensity # S + dataOut.dataPP_POWER = data_power # P + dataOut.dataPP_DOP = data_velocity + dataOut.dataPP_SNR = data_snrPP + dataOut.dataPP_WIDTH = data_specwidth dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo. dataOut.utctime = avgdatatime dataOut.flagNoData = False return dataOut + # import collections # from scipy.stats import mode #