##// END OF EJS Templates
schain
RSS
Noise.py:...
Victor Sarmiento -
r93:1e5ff7f4bb15
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,108 +1,112
1 import numpy
1 import numpy
2
2 from Model.Spectra import Spectra
3 from Model.JROHeader import *
4 from Model.Voltage import Voltage
5
3
6 def hildebrand_sekhon(Data, navg=1 ):
4 def hildebrand_sekhon(Data, navg=1):
7 """
5 """
8 This method is for the objective determination of de noise level in Doppler spectra. This
6 This method is for the objective determination of de noise level in Doppler spectra. This
9 implementation technique is based on the fact that the standard deviation of the spectral
7 implementation technique is based on the fact that the standard deviation of the spectral
10 densities is equal to the mean spectral density for white Gaussian noise
8 densities is equal to the mean spectral density for white Gaussian noise
11
9
12 Inputs:
10 Inputs:
13 Data : heights
11 Data : heights
14 navg : numbers of averages
12 navg : numbers of averages
15
13
16 Return:
14 Return:
17 -1 : any error
15 -1 : any error
18 anoise : noise's level
16 anoise : noise's level
19 """
17 """
20 divisor = 8
18 divisor = 8
21 ratio = 7 / divisor
19 ratio = 7 / divisor
22 data = Data.reshape(-1)
20 data = Data.reshape(-1)
23 npts = data.size #numbers of points of the data
21 npts = data.size #numbers of points of the data
24
22
25 if npts < 32:
23 if npts < 32:
26 print "error in noise - requires at least 32 points"
24 print "error in noise - requires at least 32 points"
27 return -1.0
25 return -1.0
28
26
29 # data sorted in ascending order
27 # data sorted in ascending order
30 nmin = int(npts/divisor + ratio);
28 nmin = int(npts/divisor + ratio);
31 s = 0.0
29 s = 0.0
32 s2 = 0.0
30 s2 = 0.0
33 data2 = data[:npts]
31 data2 = data[:npts]
34 data2.sort()
32 data2.sort()
35
33
36 for i in range(nmin):
34 for i in range(nmin):
37 s += data2[i]
35 s += data2[i]
38 s2 += data2[i]**2;
36 s2 += data2[i]**2;
39
37
40 icount = nmin
38 icount = nmin
41 iflag = 0
39 iflag = 0
42
40
43 for i in range(nmin, npts):
41 for i in range(nmin, npts):
44 s += data2[i];
42 s += data2[i];
45 s2 += data2[i]**2
43 s2 += data2[i]**2
46 icount=icount+1;
44 icount=icount+1;
47 p = s / float(icount);
45 p = s / float(icount);
48 p2 = p**2;
46 p2 = p**2;
49 q = s2 / float(icount) - p2;
47 q = s2 / float(icount) - p2;
50 leftc = p2;
48 leftc = p2;
51 rightc = q * float(navg);
49 rightc = q * float(navg);
52
50
53 if leftc > rightc:
51 if leftc > rightc:
54 iflag = 1; #No weather signal
52 iflag = 1; #No weather signal
55 # Signal detect: R2 < 1 (R2 = leftc/rightc)
53 # Signal detect: R2 < 1 (R2 = leftc/rightc)
56 if(leftc < rightc):
54 if(leftc < rightc):
57 if iflag:
55 if iflag:
58 break
56 break
59
57
60 anoise = 0.0;
58 anoise = 0.0;
61 for j in range(i):
59 for j in range(i):
62 anoise += data2[j];
60 anoise += data2[j];
63
61
64 anoise = anoise / float(i);
62 anoise = anoise / float(i);
65
63
66 return anoise;
64 return anoise;
67
65
68
66
69 class Noise():
67 class Noise():
70 """
68 """
71 Clase que implementa los metodos necesarios para deternimar el nivel de ruido en un Spectro Doppler
69 Clase que implementa los metodos necesarios para deternimar el nivel de ruido en un Spectro Doppler
72 """
70 """
73 m_DataObj = None
71 m_DataObj = None
74
72
75
73
76 def __init__(self, m_Voltage=None):
74 def __init__(self, m_Spectra=None):
77 """
75 """
78 Inicializador de la clase Noise para la la determinacion del nivel de ruido en un Spectro Doppler.
76 Inicializador de la clase Noise para la la determinacion del nivel de ruido en un Spectro Doppler.
79
77
80 Affected:
78 Affected:
81 self.m_DataObj
79 self.m_DataObj
82
80
83 Return:
81 Return:
84 None
82 None
85 """
83 """
86 # if m_Voltage == None:
84 if m_Spectra == None:
87 # m_Voltage = Voltage()
85 m_Spectra = Spectra()
88 #
86
89 # if not(isinstance(m_Voltage, Voltage)):
87 if not(isinstance(m_Spectra, Spectra)):
90 # raise ValueError, "in Noise class, m_Voltage must be an Voltage class object"
88 raise ValueError, "in Noise class, m_Spectra must be an Spectra class object"
91
89
92 self.m_DataObj = m_Voltage
90 self.m_DataObj = m_Spectra
93
91
94
92
95 def getNoiseLevelByHildebrandSekhon(self):
93 def getNoiseLevelByHildebrandSekhon(self):
96 """
94 """
97 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
95 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
98
96
99 Return:
97 Return:
100 noise level
98 noise level
101 """
99 """
102 data = self.m_DataObj.data # heights x perfiles
100 data = self.m_DataObj.data_spc
103 #heights = numpy.transpose( Data, (2,0,1) ) # channel x profile x height
101 daux = None
104 #data = Data[0,0,:]
102
103 for channel in range(self.m_DataObj.nChannels):
104 daux = data[channel,:,:]
105 noiselevel = hildebrand_sekhon(daux)
106 print noiselevel
105
107
106 noiselevel = hildebrand_sekhon(data)
107
108
109 for pair in range(self.m_DataObj.nPairs):
110 daux = data[pair,:,:]
111 noiselevel = hildebrand_sekhon(daux)
108 print noiselevel
112 print noiselevel
General Comments 0
You need to be logged in to leave comments. Login now