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)
-                     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
              #

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