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