|
|
import numpy
|
|
|
from Model.Spectra import Spectra
|
|
|
|
|
|
def hildebrand_sekhon(Data, navg=1):
|
|
|
"""
|
|
|
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
|
|
|
"""
|
|
|
divisor = 8
|
|
|
ratio = 7 / divisor
|
|
|
data = Data.reshape(-1)
|
|
|
npts = data.size #numbers of points of the data
|
|
|
|
|
|
if npts < 32:
|
|
|
print "error in noise - requires at least 32 points"
|
|
|
return -1.0
|
|
|
|
|
|
# data sorted in ascending order
|
|
|
nmin = int(npts/divisor + ratio);
|
|
|
s = 0.0
|
|
|
s2 = 0.0
|
|
|
data2 = data[:npts]
|
|
|
data2.sort()
|
|
|
|
|
|
for i in range(nmin):
|
|
|
s += data2[i]
|
|
|
s2 += data2[i]**2;
|
|
|
|
|
|
icount = nmin
|
|
|
iflag = 0
|
|
|
|
|
|
for i in range(nmin, npts):
|
|
|
s += data2[i];
|
|
|
s2 += data2[i]**2
|
|
|
icount=icount+1;
|
|
|
p = s / float(icount);
|
|
|
p2 = p**2;
|
|
|
q = s2 / float(icount) - p2;
|
|
|
leftc = p2;
|
|
|
rightc = q * float(navg);
|
|
|
|
|
|
if leftc > rightc:
|
|
|
iflag = 1; #No weather signal
|
|
|
# Signal detect: R2 < 1 (R2 = leftc/rightc)
|
|
|
if(leftc < rightc):
|
|
|
if iflag:
|
|
|
break
|
|
|
|
|
|
anoise = 0.0;
|
|
|
for j in range(i):
|
|
|
anoise += data2[j];
|
|
|
|
|
|
anoise = anoise / float(i);
|
|
|
|
|
|
return anoise;
|
|
|
|
|
|
|
|
|
class Noise:
|
|
|
"""
|
|
|
Clase que implementa los metodos necesarios para deternimar el nivel de ruido en un Spectro Doppler
|
|
|
"""
|
|
|
m_DataObj = None
|
|
|
|
|
|
|
|
|
def __init__(self, m_Spectra=None):
|
|
|
"""
|
|
|
Inicializador de la clase Noise para la la determinacion del nivel de ruido en un Spectro Doppler.
|
|
|
|
|
|
Affected:
|
|
|
self.m_DataObj
|
|
|
|
|
|
Return:
|
|
|
None
|
|
|
"""
|
|
|
if m_Spectra == None:
|
|
|
m_Spectra = Spectra()
|
|
|
|
|
|
if not(isinstance(m_Spectra, Spectra)):
|
|
|
raise ValueError, "in Noise class, m_Spectra must be an Spectra class object"
|
|
|
|
|
|
self.m_DataObj = m_Spectra
|
|
|
|
|
|
|
|
|
def getNoiseLevelByHildebrandSekhon(self):
|
|
|
"""
|
|
|
Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
|
|
|
|
|
|
Return:
|
|
|
noise level
|
|
|
"""
|
|
|
data = self.m_DataObj.data_spc
|
|
|
daux = None
|
|
|
|
|
|
for channel in range(self.m_DataObj.nChannels):
|
|
|
daux = data[channel,:,:]
|
|
|
noiselevel = hildebrand_sekhon(daux)
|
|
|
print noiselevel
|
|
|
|
|
|
|
|
|
for pair in range(self.m_DataObj.nPairs):
|
|
|
daux = data[pair,:,:]
|
|
|
noiselevel = hildebrand_sekhon(daux)
|
|
|
print noiselevel
|
|
|
|
