diff --git a/schainpy/model/jrodata.py b/schainpy/model/jrodata.py index 05737ba..17c11e9 100644 --- a/schainpy/model/jrodata.py +++ b/schainpy/model/jrodata.py @@ -11,73 +11,20 @@ import datetime from jroheaderIO import SystemHeader, RadarControllerHeader -def hildebrand_sekhon(data, navg): - """ - This method is for the objective determination of de noise level in Doppler spectra. This - implementation technique is based on the fact that the standard deviation of the spectral - densities is equal to the mean spectral density for white Gaussian noise - - Inputs: - Data : heights - navg : numbers of averages - - Return: - -1 : any error - anoise : noise's level - """ - - dataflat = data.copy().reshape(-1) - dataflat.sort() - npts = dataflat.size #numbers of points of the data - npts_noise = 0.2*npts - - if npts < 32: - print "error in noise - requires at least 32 points" - return -1.0 - - dataflat2 = numpy.power(dataflat,2) - - cs = numpy.cumsum(dataflat) - cs2 = numpy.cumsum(dataflat2) - - # data sorted in ascending order - nmin = int((npts + 7.)/8) - - for i in range(nmin, npts): - s = cs[i] - s2 = cs2[i] - p = s / float(i); - p2 = p**2; - q = s2 / float(i) - p2; - leftc = p2; - rightc = q * float(navg); - R2 = leftc/rightc - - # Signal detect: R2 < 1 (R2 = leftc/rightc) - if R2 < 1: - npts_noise = i - break - - - anoise = numpy.average(dataflat[0:npts_noise]) - - return anoise; -def sorting_bruce(data, navg): +def hildebrand_sekhon(data, navg): data = data.copy() - sortdata = numpy.sort(data) - lenOfData = len(data) + sortdata = numpy.sort(data,axis=None) + lenOfData = len(sortdata) nums_min = lenOfData/10 - if (lenOfData/10) > 0: + if (lenOfData/10) > 2: nums_min = lenOfData/10 else: - nums_min = 0 - - rtest = 1.0 + 1.0/navg - + nums_min = 2 + sump = 0. sumq = 0. @@ -95,17 +42,15 @@ def sorting_bruce(data, navg): j += 1 if j > nums_min: - if ((sumq*j) <= (rtest*sump**2)): - lnoise = sump / j - else: + rtest = float(j)/(j-1) + 1.0/navg + if ((sumq*j) > (rtest*sump**2)): j = j - 1 sump = sump - sortdata[j] sumq = sumq - sortdata[j]**2 cont = 0 - - if j == nums_min: - lnoise = sump /j - + + lnoise = sump /j + stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1)) return lnoise class JROData: @@ -422,6 +367,8 @@ class Spectra(JROData): self.flagShiftFFT = False self.ippFactor = 1 + + self.noise = None def getNoisebyHildebrand(self): """ @@ -436,48 +383,12 @@ class Spectra(JROData): self.noise[channel] = hildebrand_sekhon(daux, self.nIncohInt) return self.noise - - def getNoisebyWindow(self, heiIndexMin=0, heiIndexMax=-1, freqIndexMin=0, freqIndexMax=-1): - """ - Determina el ruido del canal utilizando la ventana indicada con las coordenadas: - (heiIndexMIn, freqIndexMin) hasta (heiIndexMax, freqIndexMAx) - - Inputs: - heiIndexMin: Limite inferior del eje de alturas - heiIndexMax: Limite superior del eje de alturas - freqIndexMin: Limite inferior del eje de frecuencia - freqIndexMax: Limite supoerior del eje de frecuencia - """ - - data = self.data_spc[:, heiIndexMin:heiIndexMax, freqIndexMin:freqIndexMax] - for channel in range(self.nChannels): - daux = data[channel,:,:] - self.noise[channel] = numpy.average(daux) - - return self.noise - - def getNoisebySort(self): - - for channel in range(self.nChannels): - daux = self.data_spc[channel,:,:] - self.noise[channel] = sorting_bruce(daux, self.nIncohInt) - - return self.noise - def getNoise(self, type = 1): if self.noise == None: self.noise = numpy.zeros(self.nChannels) - - if type == 1: - self.noise = self.getNoisebyHildebrand() - - if type == 2: - self.noise = self.getNoisebySort() - - if type == 3: - self.noise = self.getNoisebyWindow() - + self.noise = self.getNoisebyHildebrand() + return self.noise diff --git a/schainpy/model/jroprocessing.py b/schainpy/model/jroprocessing.py index 169bfb3..e3e99c0 100644 --- a/schainpy/model/jroprocessing.py +++ b/schainpy/model/jroprocessing.py @@ -1195,14 +1195,44 @@ class SpectraProc(ProcessingUnit): return 1 - def getNoise(self, minHei, maxHei): + 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): - raise ValueError, "some value in (%d,%d) is not valid" % (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]): + 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 @@ -1226,8 +1256,22 @@ class SpectraProc(ProcessingUnit): if (maxIndex >= self.dataOut.nHeights): maxIndex = self.dataOut.nHeights-1 - data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+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):