| @@ -0,0 +1,108 | |||
|
|
1 | import numpy | |
|
|
2 | ||
|
|
3 | from Model.JROHeader import * | |
|
|
4 | from Model.Voltage import Voltage | |
|
|
5 | ||
|
|
6 | def hildebrand_sekhon(Data, navg=1 ): | |
|
|
7 | """ | |
|
|
8 | 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 | |
|
|
10 | densities is equal to the mean spectral density for white Gaussian noise | |
|
|
11 | ||
|
|
12 | Inputs: | |
|
|
13 | Data : heights | |
|
|
14 | navg : numbers of averages | |
|
|
15 | ||
|
|
16 | Return: | |
|
|
17 | -1 : any error | |
|
|
18 | anoise : noise's level | |
|
|
19 | """ | |
|
|
20 | divisor = 8 | |
|
|
21 | ratio = 7 / divisor | |
|
|
22 | data = Data.reshape(-1) | |
|
|
23 | npts = data.size #numbers of points of the data | |
|
|
24 | ||
|
|
25 | if npts < 32: | |
|
|
26 | print "error in noise - requires at least 32 points" | |
|
|
27 | return -1.0 | |
|
|
28 | ||
|
|
29 | # data sorted in ascending order | |
|
|
30 | nmin = int(npts/divisor + ratio); | |
|
|
31 | s = 0.0 | |
|
|
32 | s2 = 0.0 | |
|
|
33 | data2 = data[:npts] | |
|
|
34 | data2.sort() | |
|
|
35 | ||
|
|
36 | for i in range(nmin): | |
|
|
37 | s += data2[i] | |
|
|
38 | s2 += data2[i]**2; | |
|
|
39 | ||
|
|
40 | icount = nmin | |
|
|
41 | iflag = 0 | |
|
|
42 | ||
|
|
43 | for i in range(nmin, npts): | |
|
|
44 | s += data2[i]; | |
|
|
45 | s2 += data2[i]**2 | |
|
|
46 | icount=icount+1; | |
|
|
47 | p = s / float(icount); | |
|
|
48 | p2 = p**2; | |
|
|
49 | q = s2 / float(icount) - p2; | |
|
|
50 | leftc = p2; | |
|
|
51 | rightc = q * float(navg); | |
|
|
52 | ||
|
|
53 | if leftc > rightc: | |
|
|
54 | iflag = 1; #No weather signal | |
|
|
55 | # Signal detect: R2 < 1 (R2 = leftc/rightc) | |
|
|
56 | if(leftc < rightc): | |
|
|
57 | if iflag: | |
|
|
58 | break | |
|
|
59 | ||
|
|
60 | anoise = 0.0; | |
|
|
61 | for j in range(i): | |
|
|
62 | anoise += data2[j]; | |
|
|
63 | ||
|
|
64 | anoise = anoise / float(i); | |
|
|
65 | ||
|
|
66 | return anoise; | |
|
|
67 | ||
|
|
68 | ||
|
|
69 | class Noise(): | |
|
|
70 | """ | |
|
|
71 | Clase que implementa los metodos necesarios para deternimar el nivel de ruido en un Spectro Doppler | |
|
|
72 | """ | |
|
|
73 | m_DataObj = None | |
|
|
74 | ||
|
|
75 | ||
|
|
76 | def __init__(self, m_Voltage=None): | |
|
|
77 | """ | |
|
|
78 | Inicializador de la clase Noise para la la determinacion del nivel de ruido en un Spectro Doppler. | |
|
|
79 | ||
|
|
80 | Affected: | |
|
|
81 | self.m_DataObj | |
|
|
82 | ||
|
|
83 | Return: | |
|
|
84 | None | |
|
|
85 | """ | |
|
|
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" | |
|
|
91 | ||
|
|
92 | self.m_DataObj = m_Voltage | |
|
|
93 | ||
|
|
94 | ||
|
|
95 | def getNoiseLevelByHildebrandSekhon(self): | |
|
|
96 | """ | |
|
|
97 | Determino el nivel de ruido usando el metodo Hildebrand-Sekhon | |
|
|
98 | ||
|
|
99 | Return: | |
|
|
100 | noise level | |
|
|
101 | """ | |
|
|
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,:] | |
|
|
105 | ||
|
|
106 | noiselevel = hildebrand_sekhon(data) | |
|
|
107 | ||
|
|
108 | print noiselevel | |
General Comments 0
You need to be logged in to leave comments.
Login now