schain Commit - r1314:d42d076a7255 · Jicamarca Repository

modifiacion de la clase pulse pair, correcion de errores del Noise con removeDC y los 4 momentos

avaldez -

r1314:d42d076a7255

parent child

Context file:

r1314:d42d076a7255

schainpy/model/proc/jroproc_voltage.py +59 -36

             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
                                                        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):
                     '''
                     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)
+                    #················· Calculo de Frecuencia y Velocidad doppler········
-                    data_intensity   = lag_0/(self.n*self.nCohInt)#*self.nCohInt)
                     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_freq        = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
-                    data_velocity    = (-1.0*self.lambda_/(4*math.pi*self.ippSec))*numpy.angle(lag_1)#self.ippSec*self.nCohInt
+                    data_velocity    = (self.lambda_/2.0)*data_freq
-                    self.noise  = numpy.zeros([self.__nch,self.__nHeis])
+                    #················ Potencia promedio estimada de la Senal···········
-                    for i in range(self.__nch):
+                    lag_0            = lag_0/self.n
-                        self.noise[i]=dataOut.getNoise(channel=i)
+                    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*self.nCohInt))*tmp*numpy.sign(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
                     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])
+                        data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n   = self.pushData(dataOut=dataOut)
-                        #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)
                         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_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
-                    data_intensity, data_velocity,data_snrPP,data_specwidth = self.pulsePairbyProfiles(dataOut)
                     self.__lastdatatime           = datatime
-                    if data_intensity is None:
+                    if data_power is None:
-                        return None, None,None,None,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.dataPP_POW      = data_intensity # S
-                        dataOut.data_velocity   = data_velocity  #valor para velocidad
+                        dataOut.dataPP_POWER    = data_power     # P
-                        dataOut.data_snrPP      = data_snrPP     # valor para snr
+                        dataOut.dataPP_DOP      = data_velocity
-                        dataOut.data_specwidth  = data_specwidth
+                        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
             #

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