##// END OF EJS Templates
schain
RSS
Fix bug in noise estimation
jespinoza -
r1364:ec95f8e6ce79
parent child
Show More
Side by Side Unified
Context file:
r1364:ec95f8e6ce79 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,82 +1,82
1 #include <Python.h>
1 #include <Python.h>
2 #include <numpy/arrayobject.h>
2 #include <numpy/arrayobject.h>
3 #include <math.h>
3 #include <math.h>
4
4
5
5
6 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) {
6 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) {
7 double navg;
7 double navg;
8 PyObject *data_obj, *data_array;
8 PyObject *data_obj, *data_array;
9
9
10 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) {
10 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) {
11 return NULL;
11 return NULL;
12 }
12 }
13
13
14 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
14 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
15
15
16 if (data_array == NULL) {
16 if (data_array == NULL) {
17 Py_XDECREF(data_array);
17 Py_XDECREF(data_array);
18 Py_XDECREF(data_obj);
18 Py_XDECREF(data_obj);
19 return NULL;
19 return NULL;
20 }
20 }
21 double *sortdata = (double*)PyArray_DATA(data_array);
21 double *sortdata = (double*)PyArray_DATA(data_array);
22 int lenOfData = (int)PyArray_SIZE(data_array) ;
22 int lenOfData = (int)PyArray_SIZE(data_array) ;
23 double nums_min = lenOfData*0.2;
23 double nums_min = lenOfData*0.2;
24 if (nums_min <= 5) nums_min = 5;
24 if (nums_min <= 5) nums_min = 5;
25 double sump = 0;
25 double sump = 0;
26 double sumq = 0;
26 double sumq = 0;
27 int j = 0;
27 long j = 0;
28 int cont = 1;
28 int cont = 1;
29 double rtest = 0;
29 double rtest = 0;
30 while ((cont == 1) && (j < lenOfData)) {
30 while ((cont == 1) && (j < lenOfData)) {
31 sump = sump + sortdata[j];
31 sump = sump + sortdata[j];
32 sumq = sumq + pow(sortdata[j], 2);
32 sumq = sumq + pow(sortdata[j], 2);
33 if (j > nums_min) {
33 if (j > nums_min) {
34 rtest = (double)j/(j-1) + 1/navg;
34 rtest = (double)j/(j-1) + 1/navg;
35 if ((sumq*j) > (rtest*pow(sump, 2))) {
35 if ((sumq*j) > (rtest*pow(sump, 2))) {
36 j = j - 1;
36 j = j - 1;
37 sump = sump - sortdata[j];
37 sump = sump - sortdata[j];
38 sumq = sumq - pow(sortdata[j],2);
38 sumq = sumq - pow(sortdata[j],2);
39 cont = 0;
39 cont = 0;
40 }
40 }
41 }
41 }
42 j = j + 1;
42 j = j + 1;
43 }
43 }
44
44
45 double lnoise = sump / j;
45 double lnoise = sump / j;
46
46
47 Py_DECREF(data_array);
47 Py_DECREF(data_array);
48
48
49 return PyLong_FromLong(lnoise);
49 // return PyLong_FromLong(lnoise);
50 //return Py_BuildValue("d", lnoise);
50 return PyFloat_FromDouble(lnoise);
51 }
51 }
52
52
53
53
54 static PyMethodDef noiseMethods[] = {
54 static PyMethodDef noiseMethods[] = {
55 { "hildebrand_sekhon", hildebrand_sekhon, METH_VARARGS, "Get noise with hildebrand_sekhon algorithm" },
55 { "hildebrand_sekhon", hildebrand_sekhon, METH_VARARGS, "Get noise with hildebrand_sekhon algorithm" },
56 { NULL, NULL, 0, NULL }
56 { NULL, NULL, 0, NULL }
57 };
57 };
58
58
59 #if PY_MAJOR_VERSION >= 3
59 #if PY_MAJOR_VERSION >= 3
60
60
61 static struct PyModuleDef noisemodule = {
61 static struct PyModuleDef noisemodule = {
62 PyModuleDef_HEAD_INIT,
62 PyModuleDef_HEAD_INIT,
63 "_noise",
63 "_noise",
64 "Get noise with hildebrand_sekhon algorithm",
64 "Get noise with hildebrand_sekhon algorithm",
65 -1,
65 -1,
66 noiseMethods
66 noiseMethods
67 };
67 };
68
68
69 #endif
69 #endif
70
70
71 #if PY_MAJOR_VERSION >= 3
71 #if PY_MAJOR_VERSION >= 3
72 PyMODINIT_FUNC PyInit__noise(void) {
72 PyMODINIT_FUNC PyInit__noise(void) {
73 Py_Initialize();
73 Py_Initialize();
74 import_array();
74 import_array();
75 return PyModule_Create(&noisemodule);
75 return PyModule_Create(&noisemodule);
76 }
76 }
77 #else
77 #else
78 PyMODINIT_FUNC init_noise() {
78 PyMODINIT_FUNC init_noise() {
79 Py_InitModule("_noise", noiseMethods);
79 Py_InitModule("_noise", noiseMethods);
80 import_array();
80 import_array();
81 }
81 }
82 #endif
82 #endif
Show file before | Show file after | Hide comments
@@ -1,1066 +1,1069
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 # All rights reserved.
2 # All rights reserved.
3 #
3 #
4 # Distributed under the terms of the BSD 3-clause license.
4 # Distributed under the terms of the BSD 3-clause license.
5 """Definition of diferent Data objects for different types of data
5 """Definition of diferent Data objects for different types of data
6
6
7 Here you will find the diferent data objects for the different types
7 Here you will find the diferent data objects for the different types
8 of data, this data objects must be used as dataIn or dataOut objects in
8 of data, this data objects must be used as dataIn or dataOut objects in
9 processing units and operations. Currently the supported data objects are:
9 processing units and operations. Currently the supported data objects are:
10 Voltage, Spectra, SpectraHeis, Fits, Correlation and Parameters
10 Voltage, Spectra, SpectraHeis, Fits, Correlation and Parameters
11 """
11 """
12
12
13 import copy
13 import copy
14 import numpy
14 import numpy
15 import datetime
15 import datetime
16 import json
16 import json
17
17
18 import schainpy.admin
18 import schainpy.admin
19 from schainpy.utils import log
19 from schainpy.utils import log
20 from .jroheaderIO import SystemHeader, RadarControllerHeader
20 from .jroheaderIO import SystemHeader, RadarControllerHeader
21 from schainpy.model.data import _noise
21 from schainpy.model.data import _noise
22
22
23
23
24 def getNumpyDtype(dataTypeCode):
24 def getNumpyDtype(dataTypeCode):
25
25
26 if dataTypeCode == 0:
26 if dataTypeCode == 0:
27 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
27 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
28 elif dataTypeCode == 1:
28 elif dataTypeCode == 1:
29 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
29 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
30 elif dataTypeCode == 2:
30 elif dataTypeCode == 2:
31 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
31 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
32 elif dataTypeCode == 3:
32 elif dataTypeCode == 3:
33 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
33 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
34 elif dataTypeCode == 4:
34 elif dataTypeCode == 4:
35 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
35 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
36 elif dataTypeCode == 5:
36 elif dataTypeCode == 5:
37 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
37 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
38 else:
38 else:
39 raise ValueError('dataTypeCode was not defined')
39 raise ValueError('dataTypeCode was not defined')
40
40
41 return numpyDtype
41 return numpyDtype
42
42
43
43
44 def getDataTypeCode(numpyDtype):
44 def getDataTypeCode(numpyDtype):
45
45
46 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
46 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
47 datatype = 0
47 datatype = 0
48 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
48 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
49 datatype = 1
49 datatype = 1
50 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
50 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
51 datatype = 2
51 datatype = 2
52 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
52 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
53 datatype = 3
53 datatype = 3
54 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
54 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
55 datatype = 4
55 datatype = 4
56 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
56 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
57 datatype = 5
57 datatype = 5
58 else:
58 else:
59 datatype = None
59 datatype = None
60
60
61 return datatype
61 return datatype
62
62
63
63
64 def hildebrand_sekhon(data, navg):
64 def hildebrand_sekhon(data, navg):
65 """
65 """
66 This method is for the objective determination of the noise level in Doppler spectra. This
66 This method is for the objective determination of the noise level in Doppler spectra. This
67 implementation technique is based on the fact that the standard deviation of the spectral
67 implementation technique is based on the fact that the standard deviation of the spectral
68 densities is equal to the mean spectral density for white Gaussian noise
68 densities is equal to the mean spectral density for white Gaussian noise
69
69
70 Inputs:
70 Inputs:
71 Data : heights
71 Data : heights
72 navg : numbers of averages
72 navg : numbers of averages
73
73
74 Return:
74 Return:
75 mean : noise's level
75 mean : noise's level
76 """
76 """
77
77
78 sortdata = numpy.sort(data, axis=None)
78 sortdata = numpy.sort(data, axis=None)
79 '''
79 '''
80 lenOfData = len(sortdata)
80 lenOfData = len(sortdata)
81 nums_min = lenOfData*0.2
81 nums_min = lenOfData*0.2
82
82
83 if nums_min <= 5:
83 if nums_min <= 5:
84
84
85 nums_min = 5
85 nums_min = 5
86
86
87 sump = 0.
87 sump = 0.
88 sumq = 0.
88 sumq = 0.
89
89
90 j = 0
90 j = 0
91 cont = 1
91 cont = 1
92
92
93 while((cont == 1)and(j < lenOfData)):
93 while((cont == 1)and(j < lenOfData)):
94
94
95 sump += sortdata[j]
95 sump += sortdata[j]
96 sumq += sortdata[j]**2
96 sumq += sortdata[j]**2
97
97
98 if j > nums_min:
98 if j > nums_min:
99 rtest = float(j)/(j-1) + 1.0/navg
99 rtest = float(j)/(j-1) + 1.0/navg
100 if ((sumq*j) > (rtest*sump**2)):
100 if ((sumq*j) > (rtest*sump**2)):
101 j = j - 1
101 j = j - 1
102 sump = sump - sortdata[j]
102 sump = sump - sortdata[j]
103 sumq = sumq - sortdata[j]**2
103 sumq = sumq - sortdata[j]**2
104 cont = 0
104 cont = 0
105
105
106 j += 1
106 j += 1
107
107
108 lnoise = sump / j
108 lnoise = sump / j
109 '''
109 '''
110 return _noise.hildebrand_sekhon(sortdata, navg)
110 return _noise.hildebrand_sekhon(sortdata, navg)
111
111
112
112
113 class Beam:
113 class Beam:
114
114
115 def __init__(self):
115 def __init__(self):
116 self.codeList = []
116 self.codeList = []
117 self.azimuthList = []
117 self.azimuthList = []
118 self.zenithList = []
118 self.zenithList = []
119
119
120
120
121 class GenericData(object):
121 class GenericData(object):
122
122
123 flagNoData = True
123 flagNoData = True
124
124
125 def copy(self, inputObj=None):
125 def copy(self, inputObj=None):
126
126
127 if inputObj == None:
127 if inputObj == None:
128 return copy.deepcopy(self)
128 return copy.deepcopy(self)
129
129
130 for key in list(inputObj.__dict__.keys()):
130 for key in list(inputObj.__dict__.keys()):
131
131
132 attribute = inputObj.__dict__[key]
132 attribute = inputObj.__dict__[key]
133
133
134 # If this attribute is a tuple or list
134 # If this attribute is a tuple or list
135 if type(inputObj.__dict__[key]) in (tuple, list):
135 if type(inputObj.__dict__[key]) in (tuple, list):
136 self.__dict__[key] = attribute[:]
136 self.__dict__[key] = attribute[:]
137 continue
137 continue
138
138
139 # If this attribute is another object or instance
139 # If this attribute is another object or instance
140 if hasattr(attribute, '__dict__'):
140 if hasattr(attribute, '__dict__'):
141 self.__dict__[key] = attribute.copy()
141 self.__dict__[key] = attribute.copy()
142 continue
142 continue
143
143
144 self.__dict__[key] = inputObj.__dict__[key]
144 self.__dict__[key] = inputObj.__dict__[key]
145
145
146 def deepcopy(self):
146 def deepcopy(self):
147
147
148 return copy.deepcopy(self)
148 return copy.deepcopy(self)
149
149
150 def isEmpty(self):
150 def isEmpty(self):
151
151
152 return self.flagNoData
152 return self.flagNoData
153
153
154 def isReady(self):
154 def isReady(self):
155
155
156 return not self.flagNoData
156 return not self.flagNoData
157
157
158
158
159 class JROData(GenericData):
159 class JROData(GenericData):
160
160
161 systemHeaderObj = SystemHeader()
161 systemHeaderObj = SystemHeader()
162 radarControllerHeaderObj = RadarControllerHeader()
162 radarControllerHeaderObj = RadarControllerHeader()
163 type = None
163 type = None
164 datatype = None # dtype but in string
164 datatype = None # dtype but in string
165 nProfiles = None
165 nProfiles = None
166 heightList = None
166 heightList = None
167 channelList = None
167 channelList = None
168 flagDiscontinuousBlock = False
168 flagDiscontinuousBlock = False
169 useLocalTime = False
169 useLocalTime = False
170 utctime = None
170 utctime = None
171 timeZone = None
171 timeZone = None
172 dstFlag = None
172 dstFlag = None
173 errorCount = None
173 errorCount = None
174 blocksize = None
174 blocksize = None
175 flagDecodeData = False # asumo q la data no esta decodificada
175 flagDecodeData = False # asumo q la data no esta decodificada
176 flagDeflipData = False # asumo q la data no esta sin flip
176 flagDeflipData = False # asumo q la data no esta sin flip
177 flagShiftFFT = False
177 flagShiftFFT = False
178 nCohInt = None
178 nCohInt = None
179 windowOfFilter = 1
179 windowOfFilter = 1
180 C = 3e8
180 C = 3e8
181 frequency = 49.92e6
181 frequency = 49.92e6
182 realtime = False
182 realtime = False
183 beacon_heiIndexList = None
183 beacon_heiIndexList = None
184 last_block = None
184 last_block = None
185 blocknow = None
185 blocknow = None
186 azimuth = None
186 azimuth = None
187 zenith = None
187 zenith = None
188 beam = Beam()
188 beam = Beam()
189 profileIndex = None
189 profileIndex = None
190 error = None
190 error = None
191 data = None
191 data = None
192 nmodes = None
192 nmodes = None
193 metadata_list = ['heightList', 'timeZone', 'type']
193 metadata_list = ['heightList', 'timeZone', 'type']
194
194
195 def __str__(self):
195 def __str__(self):
196
196
197 return '{} - {}'.format(self.type, self.datatime())
197 return '{} - {}'.format(self.type, self.datatime())
198
198
199 def getNoise(self):
199 def getNoise(self):
200
200
201 raise NotImplementedError
201 raise NotImplementedError
202
202
203 @property
203 @property
204 def nChannels(self):
204 def nChannels(self):
205
205
206 return len(self.channelList)
206 return len(self.channelList)
207
207
208 @property
208 @property
209 def channelIndexList(self):
209 def channelIndexList(self):
210
210
211 return list(range(self.nChannels))
211 return list(range(self.nChannels))
212
212
213 @property
213 @property
214 def nHeights(self):
214 def nHeights(self):
215
215
216 return len(self.heightList)
216 return len(self.heightList)
217
217
218 def getDeltaH(self):
218 def getDeltaH(self):
219
219
220 return self.heightList[1] - self.heightList[0]
220 return self.heightList[1] - self.heightList[0]
221
221
222 @property
222 @property
223 def ltctime(self):
223 def ltctime(self):
224
224
225 if self.useLocalTime:
225 if self.useLocalTime:
226 return self.utctime - self.timeZone * 60
226 return self.utctime - self.timeZone * 60
227
227
228 return self.utctime
228 return self.utctime
229
229
230 @property
230 @property
231 def datatime(self):
231 def datatime(self):
232
232
233 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
233 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
234 return datatimeValue
234 return datatimeValue
235
235
236 def getTimeRange(self):
236 def getTimeRange(self):
237
237
238 datatime = []
238 datatime = []
239
239
240 datatime.append(self.ltctime)
240 datatime.append(self.ltctime)
241 datatime.append(self.ltctime + self.timeInterval + 1)
241 datatime.append(self.ltctime + self.timeInterval + 1)
242
242
243 datatime = numpy.array(datatime)
243 datatime = numpy.array(datatime)
244
244
245 return datatime
245 return datatime
246
246
247 def getFmaxTimeResponse(self):
247 def getFmaxTimeResponse(self):
248
248
249 period = (10**-6) * self.getDeltaH() / (0.15)
249 period = (10**-6) * self.getDeltaH() / (0.15)
250
250
251 PRF = 1. / (period * self.nCohInt)
251 PRF = 1. / (period * self.nCohInt)
252
252
253 fmax = PRF
253 fmax = PRF
254
254
255 return fmax
255 return fmax
256
256
257 def getFmax(self):
257 def getFmax(self):
258 PRF = 1. / (self.ippSeconds * self.nCohInt)
258 PRF = 1. / (self.ippSeconds * self.nCohInt)
259
259
260 fmax = PRF
260 fmax = PRF
261 return fmax
261 return fmax
262
262
263 def getVmax(self):
263 def getVmax(self):
264
264
265 _lambda = self.C / self.frequency
265 _lambda = self.C / self.frequency
266
266
267 vmax = self.getFmax() * _lambda / 2
267 vmax = self.getFmax() * _lambda / 2
268
268
269 return vmax
269 return vmax
270
270
271 @property
271 @property
272 def ippSeconds(self):
272 def ippSeconds(self):
273 '''
273 '''
274 '''
274 '''
275 return self.radarControllerHeaderObj.ippSeconds
275 return self.radarControllerHeaderObj.ippSeconds
276
276
277 @ippSeconds.setter
277 @ippSeconds.setter
278 def ippSeconds(self, ippSeconds):
278 def ippSeconds(self, ippSeconds):
279 '''
279 '''
280 '''
280 '''
281 self.radarControllerHeaderObj.ippSeconds = ippSeconds
281 self.radarControllerHeaderObj.ippSeconds = ippSeconds
282
282
283 @property
283 @property
284 def code(self):
284 def code(self):
285 '''
285 '''
286 '''
286 '''
287 return self.radarControllerHeaderObj.code
287 return self.radarControllerHeaderObj.code
288
288
289 @code.setter
289 @code.setter
290 def code(self, code):
290 def code(self, code):
291 '''
291 '''
292 '''
292 '''
293 self.radarControllerHeaderObj.code = code
293 self.radarControllerHeaderObj.code = code
294
294
295 @property
295 @property
296 def nCode(self):
296 def nCode(self):
297 '''
297 '''
298 '''
298 '''
299 return self.radarControllerHeaderObj.nCode
299 return self.radarControllerHeaderObj.nCode
300
300
301 @nCode.setter
301 @nCode.setter
302 def nCode(self, ncode):
302 def nCode(self, ncode):
303 '''
303 '''
304 '''
304 '''
305 self.radarControllerHeaderObj.nCode = ncode
305 self.radarControllerHeaderObj.nCode = ncode
306
306
307 @property
307 @property
308 def nBaud(self):
308 def nBaud(self):
309 '''
309 '''
310 '''
310 '''
311 return self.radarControllerHeaderObj.nBaud
311 return self.radarControllerHeaderObj.nBaud
312
312
313 @nBaud.setter
313 @nBaud.setter
314 def nBaud(self, nbaud):
314 def nBaud(self, nbaud):
315 '''
315 '''
316 '''
316 '''
317 self.radarControllerHeaderObj.nBaud = nbaud
317 self.radarControllerHeaderObj.nBaud = nbaud
318
318
319 @property
319 @property
320 def ipp(self):
320 def ipp(self):
321 '''
321 '''
322 '''
322 '''
323 return self.radarControllerHeaderObj.ipp
323 return self.radarControllerHeaderObj.ipp
324
324
325 @ipp.setter
325 @ipp.setter
326 def ipp(self, ipp):
326 def ipp(self, ipp):
327 '''
327 '''
328 '''
328 '''
329 self.radarControllerHeaderObj.ipp = ipp
329 self.radarControllerHeaderObj.ipp = ipp
330
330
331 @property
331 @property
332 def metadata(self):
332 def metadata(self):
333 '''
333 '''
334 '''
334 '''
335
335
336 return {attr: getattr(self, attr) for attr in self.metadata_list}
336 return {attr: getattr(self, attr) for attr in self.metadata_list}
337
337
338
338
339 class Voltage(JROData):
339 class Voltage(JROData):
340
340
341 dataPP_POW = None
341 dataPP_POW = None
342 dataPP_DOP = None
342 dataPP_DOP = None
343 dataPP_WIDTH = None
343 dataPP_WIDTH = None
344 dataPP_SNR = None
344 dataPP_SNR = None
345
345
346 def __init__(self):
346 def __init__(self):
347 '''
347 '''
348 Constructor
348 Constructor
349 '''
349 '''
350
350
351 self.useLocalTime = True
351 self.useLocalTime = True
352 self.radarControllerHeaderObj = RadarControllerHeader()
352 self.radarControllerHeaderObj = RadarControllerHeader()
353 self.systemHeaderObj = SystemHeader()
353 self.systemHeaderObj = SystemHeader()
354 self.type = "Voltage"
354 self.type = "Voltage"
355 self.data = None
355 self.data = None
356 self.nProfiles = None
356 self.nProfiles = None
357 self.heightList = None
357 self.heightList = None
358 self.channelList = None
358 self.channelList = None
359 self.flagNoData = True
359 self.flagNoData = True
360 self.flagDiscontinuousBlock = False
360 self.flagDiscontinuousBlock = False
361 self.utctime = None
361 self.utctime = None
362 self.timeZone = 0
362 self.timeZone = 0
363 self.dstFlag = None
363 self.dstFlag = None
364 self.errorCount = None
364 self.errorCount = None
365 self.nCohInt = None
365 self.nCohInt = None
366 self.blocksize = None
366 self.blocksize = None
367 self.flagCohInt = False
367 self.flagCohInt = False
368 self.flagDecodeData = False # asumo q la data no esta decodificada
368 self.flagDecodeData = False # asumo q la data no esta decodificada
369 self.flagDeflipData = False # asumo q la data no esta sin flip
369 self.flagDeflipData = False # asumo q la data no esta sin flip
370 self.flagShiftFFT = False
370 self.flagShiftFFT = False
371 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
371 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
372 self.profileIndex = 0
372 self.profileIndex = 0
373 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
373 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
374 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
374 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
375
375
376 def getNoisebyHildebrand(self, channel=None):
376 def getNoisebyHildebrand(self, channel=None):
377 """
377 """
378 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
378 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
379
379
380 Return:
380 Return:
381 noiselevel
381 noiselevel
382 """
382 """
383
383
384 if channel != None:
384 if channel != None:
385 data = self.data[channel]
385 data = self.data[channel]
386 nChannels = 1
386 nChannels = 1
387 else:
387 else:
388 data = self.data
388 data = self.data
389 nChannels = self.nChannels
389 nChannels = self.nChannels
390
390
391 noise = numpy.zeros(nChannels)
391 noise = numpy.zeros(nChannels)
392 power = data * numpy.conjugate(data)
392 power = data * numpy.conjugate(data)
393
393
394 for thisChannel in range(nChannels):
394 for thisChannel in range(nChannels):
395 if nChannels == 1:
395 if nChannels == 1:
396 daux = power[:].real
396 daux = power[:].real
397 else:
397 else:
398 daux = power[thisChannel, :].real
398 daux = power[thisChannel, :].real
399 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
399 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
400
400
401 return noise
401 return noise
402
402
403 def getNoise(self, type=1, channel=None):
403 def getNoise(self, type=1, channel=None):
404
404
405 if type == 1:
405 if type == 1:
406 noise = self.getNoisebyHildebrand(channel)
406 noise = self.getNoisebyHildebrand(channel)
407
407
408 return noise
408 return noise
409
409
410 def getPower(self, channel=None):
410 def getPower(self, channel=None):
411
411
412 if channel != None:
412 if channel != None:
413 data = self.data[channel]
413 data = self.data[channel]
414 else:
414 else:
415 data = self.data
415 data = self.data
416
416
417 power = data * numpy.conjugate(data)
417 power = data * numpy.conjugate(data)
418 powerdB = 10 * numpy.log10(power.real)
418 powerdB = 10 * numpy.log10(power.real)
419 powerdB = numpy.squeeze(powerdB)
419 powerdB = numpy.squeeze(powerdB)
420
420
421 return powerdB
421 return powerdB
422
422
423 @property
423 @property
424 def timeInterval(self):
424 def timeInterval(self):
425
425
426 return self.ippSeconds * self.nCohInt
426 return self.ippSeconds * self.nCohInt
427
427
428 noise = property(getNoise, "I'm the 'nHeights' property.")
428 noise = property(getNoise, "I'm the 'nHeights' property.")
429
429
430
430
431 class Spectra(JROData):
431 class Spectra(JROData):
432
432
433 def __init__(self):
433 def __init__(self):
434 '''
434 '''
435 Constructor
435 Constructor
436 '''
436 '''
437
437
438 self.data_dc = None
439 self.data_spc = None
440 self.data_cspc = None
438 self.useLocalTime = True
441 self.useLocalTime = True
439 self.radarControllerHeaderObj = RadarControllerHeader()
442 self.radarControllerHeaderObj = RadarControllerHeader()
440 self.systemHeaderObj = SystemHeader()
443 self.systemHeaderObj = SystemHeader()
441 self.type = "Spectra"
444 self.type = "Spectra"
442 self.timeZone = 0
445 self.timeZone = 0
443 self.nProfiles = None
446 self.nProfiles = None
444 self.heightList = None
447 self.heightList = None
445 self.channelList = None
448 self.channelList = None
446 self.pairsList = None
449 self.pairsList = None
447 self.flagNoData = True
450 self.flagNoData = True
448 self.flagDiscontinuousBlock = False
451 self.flagDiscontinuousBlock = False
449 self.utctime = None
452 self.utctime = None
450 self.nCohInt = None
453 self.nCohInt = None
451 self.nIncohInt = None
454 self.nIncohInt = None
452 self.blocksize = None
455 self.blocksize = None
453 self.nFFTPoints = None
456 self.nFFTPoints = None
454 self.wavelength = None
457 self.wavelength = None
455 self.flagDecodeData = False # asumo q la data no esta decodificada
458 self.flagDecodeData = False # asumo q la data no esta decodificada
456 self.flagDeflipData = False # asumo q la data no esta sin flip
459 self.flagDeflipData = False # asumo q la data no esta sin flip
457 self.flagShiftFFT = False
460 self.flagShiftFFT = False
458 self.ippFactor = 1
461 self.ippFactor = 1
459 self.beacon_heiIndexList = []
462 self.beacon_heiIndexList = []
460 self.noise_estimation = None
463 self.noise_estimation = None
461 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
464 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
462 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
465 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
463
466
464 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
467 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
465 """
468 """
466 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
469 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
467
470
468 Return:
471 Return:
469 noiselevel
472 noiselevel
470 """
473 """
471
474
472 noise = numpy.zeros(self.nChannels)
475 noise = numpy.zeros(self.nChannels)
473
476
474 for channel in range(self.nChannels):
477 for channel in range(self.nChannels):
475 daux = self.data_spc[channel,
478 daux = self.data_spc[channel,
476 xmin_index:xmax_index, ymin_index:ymax_index]
479 xmin_index:xmax_index, ymin_index:ymax_index]
477 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
480 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
478
481
479 return noise
482 return noise
480
483
481 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
484 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
482
485
483 if self.noise_estimation is not None:
486 if self.noise_estimation is not None:
484 # this was estimated by getNoise Operation defined in jroproc_spectra.py
487 # this was estimated by getNoise Operation defined in jroproc_spectra.py
485 return self.noise_estimation
488 return self.noise_estimation
486 else:
489 else:
487 noise = self.getNoisebyHildebrand(
490 noise = self.getNoisebyHildebrand(
488 xmin_index, xmax_index, ymin_index, ymax_index)
491 xmin_index, xmax_index, ymin_index, ymax_index)
489 return noise
492 return noise
490
493
491 def getFreqRangeTimeResponse(self, extrapoints=0):
494 def getFreqRangeTimeResponse(self, extrapoints=0):
492
495
493 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
496 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
494 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
497 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
495
498
496 return freqrange
499 return freqrange
497
500
498 def getAcfRange(self, extrapoints=0):
501 def getAcfRange(self, extrapoints=0):
499
502
500 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
503 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
501 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
504 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
502
505
503 return freqrange
506 return freqrange
504
507
505 def getFreqRange(self, extrapoints=0):
508 def getFreqRange(self, extrapoints=0):
506
509
507 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
510 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
508 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
511 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
509
512
510 return freqrange
513 return freqrange
511
514
512 def getVelRange(self, extrapoints=0):
515 def getVelRange(self, extrapoints=0):
513
516
514 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
517 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
515 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
518 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
516
519
517 if self.nmodes:
520 if self.nmodes:
518 return velrange/self.nmodes
521 return velrange/self.nmodes
519 else:
522 else:
520 return velrange
523 return velrange
521
524
522 @property
525 @property
523 def nPairs(self):
526 def nPairs(self):
524
527
525 return len(self.pairsList)
528 return len(self.pairsList)
526
529
527 @property
530 @property
528 def pairsIndexList(self):
531 def pairsIndexList(self):
529
532
530 return list(range(self.nPairs))
533 return list(range(self.nPairs))
531
534
532 @property
535 @property
533 def normFactor(self):
536 def normFactor(self):
534
537
535 pwcode = 1
538 pwcode = 1
536
539
537 if self.flagDecodeData:
540 if self.flagDecodeData:
538 pwcode = numpy.sum(self.code[0]**2)
541 pwcode = numpy.sum(self.code[0]**2)
539 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
542 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
540 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
543 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
541
544
542 return normFactor
545 return normFactor
543
546
544 @property
547 @property
545 def flag_cspc(self):
548 def flag_cspc(self):
546
549
547 if self.data_cspc is None:
550 if self.data_cspc is None:
548 return True
551 return True
549
552
550 return False
553 return False
551
554
552 @property
555 @property
553 def flag_dc(self):
556 def flag_dc(self):
554
557
555 if self.data_dc is None:
558 if self.data_dc is None:
556 return True
559 return True
557
560
558 return False
561 return False
559
562
560 @property
563 @property
561 def timeInterval(self):
564 def timeInterval(self):
562
565
563 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
566 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
564 if self.nmodes:
567 if self.nmodes:
565 return self.nmodes*timeInterval
568 return self.nmodes*timeInterval
566 else:
569 else:
567 return timeInterval
570 return timeInterval
568
571
569 def getPower(self):
572 def getPower(self):
570
573
571 factor = self.normFactor
574 factor = self.normFactor
572 z = self.data_spc / factor
575 z = self.data_spc / factor
573 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
576 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
574 avg = numpy.average(z, axis=1)
577 avg = numpy.average(z, axis=1)
575
578
576 return 10 * numpy.log10(avg)
579 return 10 * numpy.log10(avg)
577
580
578 def getCoherence(self, pairsList=None, phase=False):
581 def getCoherence(self, pairsList=None, phase=False):
579
582
580 z = []
583 z = []
581 if pairsList is None:
584 if pairsList is None:
582 pairsIndexList = self.pairsIndexList
585 pairsIndexList = self.pairsIndexList
583 else:
586 else:
584 pairsIndexList = []
587 pairsIndexList = []
585 for pair in pairsList:
588 for pair in pairsList:
586 if pair not in self.pairsList:
589 if pair not in self.pairsList:
587 raise ValueError("Pair %s is not in dataOut.pairsList" % (
590 raise ValueError("Pair %s is not in dataOut.pairsList" % (
588 pair))
591 pair))
589 pairsIndexList.append(self.pairsList.index(pair))
592 pairsIndexList.append(self.pairsList.index(pair))
590 for i in range(len(pairsIndexList)):
593 for i in range(len(pairsIndexList)):
591 pair = self.pairsList[pairsIndexList[i]]
594 pair = self.pairsList[pairsIndexList[i]]
592 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
595 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
593 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
596 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
594 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
597 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
595 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
598 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
596 if phase:
599 if phase:
597 data = numpy.arctan2(avgcoherenceComplex.imag,
600 data = numpy.arctan2(avgcoherenceComplex.imag,
598 avgcoherenceComplex.real) * 180 / numpy.pi
601 avgcoherenceComplex.real) * 180 / numpy.pi
599 else:
602 else:
600 data = numpy.abs(avgcoherenceComplex)
603 data = numpy.abs(avgcoherenceComplex)
601
604
602 z.append(data)
605 z.append(data)
603
606
604 return numpy.array(z)
607 return numpy.array(z)
605
608
606 def setValue(self, value):
609 def setValue(self, value):
607
610
608 print("This property should not be initialized")
611 print("This property should not be initialized")
609
612
610 return
613 return
611
614
612 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
615 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
613
616
614
617
615 class SpectraHeis(Spectra):
618 class SpectraHeis(Spectra):
616
619
617 def __init__(self):
620 def __init__(self):
618
621
619 self.radarControllerHeaderObj = RadarControllerHeader()
622 self.radarControllerHeaderObj = RadarControllerHeader()
620 self.systemHeaderObj = SystemHeader()
623 self.systemHeaderObj = SystemHeader()
621 self.type = "SpectraHeis"
624 self.type = "SpectraHeis"
622 self.nProfiles = None
625 self.nProfiles = None
623 self.heightList = None
626 self.heightList = None
624 self.channelList = None
627 self.channelList = None
625 self.flagNoData = True
628 self.flagNoData = True
626 self.flagDiscontinuousBlock = False
629 self.flagDiscontinuousBlock = False
627 self.utctime = None
630 self.utctime = None
628 self.blocksize = None
631 self.blocksize = None
629 self.profileIndex = 0
632 self.profileIndex = 0
630 self.nCohInt = 1
633 self.nCohInt = 1
631 self.nIncohInt = 1
634 self.nIncohInt = 1
632
635
633 @property
636 @property
634 def normFactor(self):
637 def normFactor(self):
635 pwcode = 1
638 pwcode = 1
636 if self.flagDecodeData:
639 if self.flagDecodeData:
637 pwcode = numpy.sum(self.code[0]**2)
640 pwcode = numpy.sum(self.code[0]**2)
638
641
639 normFactor = self.nIncohInt * self.nCohInt * pwcode
642 normFactor = self.nIncohInt * self.nCohInt * pwcode
640
643
641 return normFactor
644 return normFactor
642
645
643 @property
646 @property
644 def timeInterval(self):
647 def timeInterval(self):
645
648
646 return self.ippSeconds * self.nCohInt * self.nIncohInt
649 return self.ippSeconds * self.nCohInt * self.nIncohInt
647
650
648
651
649 class Fits(JROData):
652 class Fits(JROData):
650
653
651 def __init__(self):
654 def __init__(self):
652
655
653 self.type = "Fits"
656 self.type = "Fits"
654 self.nProfiles = None
657 self.nProfiles = None
655 self.heightList = None
658 self.heightList = None
656 self.channelList = None
659 self.channelList = None
657 self.flagNoData = True
660 self.flagNoData = True
658 self.utctime = None
661 self.utctime = None
659 self.nCohInt = 1
662 self.nCohInt = 1
660 self.nIncohInt = 1
663 self.nIncohInt = 1
661 self.useLocalTime = True
664 self.useLocalTime = True
662 self.profileIndex = 0
665 self.profileIndex = 0
663 self.timeZone = 0
666 self.timeZone = 0
664
667
665 def getTimeRange(self):
668 def getTimeRange(self):
666
669
667 datatime = []
670 datatime = []
668
671
669 datatime.append(self.ltctime)
672 datatime.append(self.ltctime)
670 datatime.append(self.ltctime + self.timeInterval)
673 datatime.append(self.ltctime + self.timeInterval)
671
674
672 datatime = numpy.array(datatime)
675 datatime = numpy.array(datatime)
673
676
674 return datatime
677 return datatime
675
678
676 def getChannelIndexList(self):
679 def getChannelIndexList(self):
677
680
678 return list(range(self.nChannels))
681 return list(range(self.nChannels))
679
682
680 def getNoise(self, type=1):
683 def getNoise(self, type=1):
681
684
682
685
683 if type == 1:
686 if type == 1:
684 noise = self.getNoisebyHildebrand()
687 noise = self.getNoisebyHildebrand()
685
688
686 if type == 2:
689 if type == 2:
687 noise = self.getNoisebySort()
690 noise = self.getNoisebySort()
688
691
689 if type == 3:
692 if type == 3:
690 noise = self.getNoisebyWindow()
693 noise = self.getNoisebyWindow()
691
694
692 return noise
695 return noise
693
696
694 @property
697 @property
695 def timeInterval(self):
698 def timeInterval(self):
696
699
697 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
700 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
698
701
699 return timeInterval
702 return timeInterval
700
703
701 @property
704 @property
702 def ippSeconds(self):
705 def ippSeconds(self):
703 '''
706 '''
704 '''
707 '''
705 return self.ipp_sec
708 return self.ipp_sec
706
709
707 noise = property(getNoise, "I'm the 'nHeights' property.")
710 noise = property(getNoise, "I'm the 'nHeights' property.")
708
711
709
712
710 class Correlation(JROData):
713 class Correlation(JROData):
711
714
712 def __init__(self):
715 def __init__(self):
713 '''
716 '''
714 Constructor
717 Constructor
715 '''
718 '''
716 self.radarControllerHeaderObj = RadarControllerHeader()
719 self.radarControllerHeaderObj = RadarControllerHeader()
717 self.systemHeaderObj = SystemHeader()
720 self.systemHeaderObj = SystemHeader()
718 self.type = "Correlation"
721 self.type = "Correlation"
719 self.data = None
722 self.data = None
720 self.dtype = None
723 self.dtype = None
721 self.nProfiles = None
724 self.nProfiles = None
722 self.heightList = None
725 self.heightList = None
723 self.channelList = None
726 self.channelList = None
724 self.flagNoData = True
727 self.flagNoData = True
725 self.flagDiscontinuousBlock = False
728 self.flagDiscontinuousBlock = False
726 self.utctime = None
729 self.utctime = None
727 self.timeZone = 0
730 self.timeZone = 0
728 self.dstFlag = None
731 self.dstFlag = None
729 self.errorCount = None
732 self.errorCount = None
730 self.blocksize = None
733 self.blocksize = None
731 self.flagDecodeData = False # asumo q la data no esta decodificada
734 self.flagDecodeData = False # asumo q la data no esta decodificada
732 self.flagDeflipData = False # asumo q la data no esta sin flip
735 self.flagDeflipData = False # asumo q la data no esta sin flip
733 self.pairsList = None
736 self.pairsList = None
734 self.nPoints = None
737 self.nPoints = None
735
738
736 def getPairsList(self):
739 def getPairsList(self):
737
740
738 return self.pairsList
741 return self.pairsList
739
742
740 def getNoise(self, mode=2):
743 def getNoise(self, mode=2):
741
744
742 indR = numpy.where(self.lagR == 0)[0][0]
745 indR = numpy.where(self.lagR == 0)[0][0]
743 indT = numpy.where(self.lagT == 0)[0][0]
746 indT = numpy.where(self.lagT == 0)[0][0]
744
747
745 jspectra0 = self.data_corr[:, :, indR, :]
748 jspectra0 = self.data_corr[:, :, indR, :]
746 jspectra = copy.copy(jspectra0)
749 jspectra = copy.copy(jspectra0)
747
750
748 num_chan = jspectra.shape[0]
751 num_chan = jspectra.shape[0]
749 num_hei = jspectra.shape[2]
752 num_hei = jspectra.shape[2]
750
753
751 freq_dc = jspectra.shape[1] / 2
754 freq_dc = jspectra.shape[1] / 2
752 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
755 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
753
756
754 if ind_vel[0] < 0:
757 if ind_vel[0] < 0:
755 ind_vel[list(range(0, 1))] = ind_vel[list(
758 ind_vel[list(range(0, 1))] = ind_vel[list(
756 range(0, 1))] + self.num_prof
759 range(0, 1))] + self.num_prof
757
760
758 if mode == 1:
761 if mode == 1:
759 jspectra[:, freq_dc, :] = (
762 jspectra[:, freq_dc, :] = (
760 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
763 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
761
764
762 if mode == 2:
765 if mode == 2:
763
766
764 vel = numpy.array([-2, -1, 1, 2])
767 vel = numpy.array([-2, -1, 1, 2])
765 xx = numpy.zeros([4, 4])
768 xx = numpy.zeros([4, 4])
766
769
767 for fil in range(4):
770 for fil in range(4):
768 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
771 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
769
772
770 xx_inv = numpy.linalg.inv(xx)
773 xx_inv = numpy.linalg.inv(xx)
771 xx_aux = xx_inv[0, :]
774 xx_aux = xx_inv[0, :]
772
775
773 for ich in range(num_chan):
776 for ich in range(num_chan):
774 yy = jspectra[ich, ind_vel, :]
777 yy = jspectra[ich, ind_vel, :]
775 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
778 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
776
779
777 junkid = jspectra[ich, freq_dc, :] <= 0
780 junkid = jspectra[ich, freq_dc, :] <= 0
778 cjunkid = sum(junkid)
781 cjunkid = sum(junkid)
779
782
780 if cjunkid.any():
783 if cjunkid.any():
781 jspectra[ich, freq_dc, junkid.nonzero()] = (
784 jspectra[ich, freq_dc, junkid.nonzero()] = (
782 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
785 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
783
786
784 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
787 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
785
788
786 return noise
789 return noise
787
790
788 @property
791 @property
789 def timeInterval(self):
792 def timeInterval(self):
790
793
791 return self.ippSeconds * self.nCohInt * self.nProfiles
794 return self.ippSeconds * self.nCohInt * self.nProfiles
792
795
793 def splitFunctions(self):
796 def splitFunctions(self):
794
797
795 pairsList = self.pairsList
798 pairsList = self.pairsList
796 ccf_pairs = []
799 ccf_pairs = []
797 acf_pairs = []
800 acf_pairs = []
798 ccf_ind = []
801 ccf_ind = []
799 acf_ind = []
802 acf_ind = []
800 for l in range(len(pairsList)):
803 for l in range(len(pairsList)):
801 chan0 = pairsList[l][0]
804 chan0 = pairsList[l][0]
802 chan1 = pairsList[l][1]
805 chan1 = pairsList[l][1]
803
806
804 # Obteniendo pares de Autocorrelacion
807 # Obteniendo pares de Autocorrelacion
805 if chan0 == chan1:
808 if chan0 == chan1:
806 acf_pairs.append(chan0)
809 acf_pairs.append(chan0)
807 acf_ind.append(l)
810 acf_ind.append(l)
808 else:
811 else:
809 ccf_pairs.append(pairsList[l])
812 ccf_pairs.append(pairsList[l])
810 ccf_ind.append(l)
813 ccf_ind.append(l)
811
814
812 data_acf = self.data_cf[acf_ind]
815 data_acf = self.data_cf[acf_ind]
813 data_ccf = self.data_cf[ccf_ind]
816 data_ccf = self.data_cf[ccf_ind]
814
817
815 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
818 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
816
819
817 @property
820 @property
818 def normFactor(self):
821 def normFactor(self):
819 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
822 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
820 acf_pairs = numpy.array(acf_pairs)
823 acf_pairs = numpy.array(acf_pairs)
821 normFactor = numpy.zeros((self.nPairs, self.nHeights))
824 normFactor = numpy.zeros((self.nPairs, self.nHeights))
822
825
823 for p in range(self.nPairs):
826 for p in range(self.nPairs):
824 pair = self.pairsList[p]
827 pair = self.pairsList[p]
825
828
826 ch0 = pair[0]
829 ch0 = pair[0]
827 ch1 = pair[1]
830 ch1 = pair[1]
828
831
829 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
832 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
830 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
833 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
831 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
834 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
832
835
833 return normFactor
836 return normFactor
834
837
835
838
836 class Parameters(Spectra):
839 class Parameters(Spectra):
837
840
838 groupList = None # List of Pairs, Groups, etc
841 groupList = None # List of Pairs, Groups, etc
839 data_param = None # Parameters obtained
842 data_param = None # Parameters obtained
840 data_pre = None # Data Pre Parametrization
843 data_pre = None # Data Pre Parametrization
841 data_SNR = None # Signal to Noise Ratio
844 data_SNR = None # Signal to Noise Ratio
842 abscissaList = None # Abscissa, can be velocities, lags or time
845 abscissaList = None # Abscissa, can be velocities, lags or time
843 utctimeInit = None # Initial UTC time
846 utctimeInit = None # Initial UTC time
844 paramInterval = None # Time interval to calculate Parameters in seconds
847 paramInterval = None # Time interval to calculate Parameters in seconds
845 useLocalTime = True
848 useLocalTime = True
846 # Fitting
849 # Fitting
847 data_error = None # Error of the estimation
850 data_error = None # Error of the estimation
848 constants = None
851 constants = None
849 library = None
852 library = None
850 # Output signal
853 # Output signal
851 outputInterval = None # Time interval to calculate output signal in seconds
854 outputInterval = None # Time interval to calculate output signal in seconds
852 data_output = None # Out signal
855 data_output = None # Out signal
853 nAvg = None
856 nAvg = None
854 noise_estimation = None
857 noise_estimation = None
855 GauSPC = None # Fit gaussian SPC
858 GauSPC = None # Fit gaussian SPC
856
859
857 def __init__(self):
860 def __init__(self):
858 '''
861 '''
859 Constructor
862 Constructor
860 '''
863 '''
861 self.radarControllerHeaderObj = RadarControllerHeader()
864 self.radarControllerHeaderObj = RadarControllerHeader()
862 self.systemHeaderObj = SystemHeader()
865 self.systemHeaderObj = SystemHeader()
863 self.type = "Parameters"
866 self.type = "Parameters"
864 self.timeZone = 0
867 self.timeZone = 0
865
868
866 def getTimeRange1(self, interval):
869 def getTimeRange1(self, interval):
867
870
868 datatime = []
871 datatime = []
869
872
870 if self.useLocalTime:
873 if self.useLocalTime:
871 time1 = self.utctimeInit - self.timeZone * 60
874 time1 = self.utctimeInit - self.timeZone * 60
872 else:
875 else:
873 time1 = self.utctimeInit
876 time1 = self.utctimeInit
874
877
875 datatime.append(time1)
878 datatime.append(time1)
876 datatime.append(time1 + interval)
879 datatime.append(time1 + interval)
877 datatime = numpy.array(datatime)
880 datatime = numpy.array(datatime)
878
881
879 return datatime
882 return datatime
880
883
881 @property
884 @property
882 def timeInterval(self):
885 def timeInterval(self):
883
886
884 if hasattr(self, 'timeInterval1'):
887 if hasattr(self, 'timeInterval1'):
885 return self.timeInterval1
888 return self.timeInterval1
886 else:
889 else:
887 return self.paramInterval
890 return self.paramInterval
888
891
889 def setValue(self, value):
892 def setValue(self, value):
890
893
891 print("This property should not be initialized")
894 print("This property should not be initialized")
892
895
893 return
896 return
894
897
895 def getNoise(self):
898 def getNoise(self):
896
899
897 return self.spc_noise
900 return self.spc_noise
898
901
899 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
902 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
900
903
901
904
902 class PlotterData(object):
905 class PlotterData(object):
903 '''
906 '''
904 Object to hold data to be plotted
907 Object to hold data to be plotted
905 '''
908 '''
906
909
907 MAXNUMX = 200
910 MAXNUMX = 200
908 MAXNUMY = 200
911 MAXNUMY = 200
909
912
910 def __init__(self, code, exp_code, localtime=True):
913 def __init__(self, code, exp_code, localtime=True):
911
914
912 self.key = code
915 self.key = code
913 self.exp_code = exp_code
916 self.exp_code = exp_code
914 self.ready = False
917 self.ready = False
915 self.flagNoData = False
918 self.flagNoData = False
916 self.localtime = localtime
919 self.localtime = localtime
917 self.data = {}
920 self.data = {}
918 self.meta = {}
921 self.meta = {}
919 self.__heights = []
922 self.__heights = []
920
923
921 def __str__(self):
924 def __str__(self):
922 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
925 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
923 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
926 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
924
927
925 def __len__(self):
928 def __len__(self):
926 return len(self.data)
929 return len(self.data)
927
930
928 def __getitem__(self, key):
931 def __getitem__(self, key):
929 if isinstance(key, int):
932 if isinstance(key, int):
930 return self.data[self.times[key]]
933 return self.data[self.times[key]]
931 elif isinstance(key, str):
934 elif isinstance(key, str):
932 ret = numpy.array([self.data[x][key] for x in self.times])
935 ret = numpy.array([self.data[x][key] for x in self.times])
933 if ret.ndim > 1:
936 if ret.ndim > 1:
934 ret = numpy.swapaxes(ret, 0, 1)
937 ret = numpy.swapaxes(ret, 0, 1)
935 return ret
938 return ret
936
939
937 def __contains__(self, key):
940 def __contains__(self, key):
938 return key in self.data[self.min_time]
941 return key in self.data[self.min_time]
939
942
940 def setup(self):
943 def setup(self):
941 '''
944 '''
942 Configure object
945 Configure object
943 '''
946 '''
944 self.type = ''
947 self.type = ''
945 self.ready = False
948 self.ready = False
946 del self.data
949 del self.data
947 self.data = {}
950 self.data = {}
948 self.__heights = []
951 self.__heights = []
949 self.__all_heights = set()
952 self.__all_heights = set()
950
953
951 def shape(self, key):
954 def shape(self, key):
952 '''
955 '''
953 Get the shape of the one-element data for the given key
956 Get the shape of the one-element data for the given key
954 '''
957 '''
955
958
956 if len(self.data[self.min_time][key]):
959 if len(self.data[self.min_time][key]):
957 return self.data[self.min_time][key].shape
960 return self.data[self.min_time][key].shape
958 return (0,)
961 return (0,)
959
962
960 def update(self, data, tm, meta={}):
963 def update(self, data, tm, meta={}):
961 '''
964 '''
962 Update data object with new dataOut
965 Update data object with new dataOut
963 '''
966 '''
964
967
965 self.data[tm] = data
968 self.data[tm] = data
966
969
967 for key, value in meta.items():
970 for key, value in meta.items():
968 setattr(self, key, value)
971 setattr(self, key, value)
969
972
970 def normalize_heights(self):
973 def normalize_heights(self):
971 '''
974 '''
972 Ensure same-dimension of the data for different heighList
975 Ensure same-dimension of the data for different heighList
973 '''
976 '''
974
977
975 H = numpy.array(list(self.__all_heights))
978 H = numpy.array(list(self.__all_heights))
976 H.sort()
979 H.sort()
977 for key in self.data:
980 for key in self.data:
978 shape = self.shape(key)[:-1] + H.shape
981 shape = self.shape(key)[:-1] + H.shape
979 for tm, obj in list(self.data[key].items()):
982 for tm, obj in list(self.data[key].items()):
980 h = self.__heights[self.times.tolist().index(tm)]
983 h = self.__heights[self.times.tolist().index(tm)]
981 if H.size == h.size:
984 if H.size == h.size:
982 continue
985 continue
983 index = numpy.where(numpy.in1d(H, h))[0]
986 index = numpy.where(numpy.in1d(H, h))[0]
984 dummy = numpy.zeros(shape) + numpy.nan
987 dummy = numpy.zeros(shape) + numpy.nan
985 if len(shape) == 2:
988 if len(shape) == 2:
986 dummy[:, index] = obj
989 dummy[:, index] = obj
987 else:
990 else:
988 dummy[index] = obj
991 dummy[index] = obj
989 self.data[key][tm] = dummy
992 self.data[key][tm] = dummy
990
993
991 self.__heights = [H for tm in self.times]
994 self.__heights = [H for tm in self.times]
992
995
993 def jsonify(self, tm, plot_name, plot_type, decimate=False):
996 def jsonify(self, tm, plot_name, plot_type, decimate=False):
994 '''
997 '''
995 Convert data to json
998 Convert data to json
996 '''
999 '''
997
1000
998 meta = {}
1001 meta = {}
999 meta['xrange'] = []
1002 meta['xrange'] = []
1000 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1003 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1001 tmp = self.data[tm][self.key]
1004 tmp = self.data[tm][self.key]
1002 shape = tmp.shape
1005 shape = tmp.shape
1003 if len(shape) == 2:
1006 if len(shape) == 2:
1004 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1007 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1005 elif len(shape) == 3:
1008 elif len(shape) == 3:
1006 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1009 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1007 data = self.roundFloats(
1010 data = self.roundFloats(
1008 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1011 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1009 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1012 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1010 else:
1013 else:
1011 data = self.roundFloats(self.data[tm][self.key].tolist())
1014 data = self.roundFloats(self.data[tm][self.key].tolist())
1012
1015
1013 ret = {
1016 ret = {
1014 'plot': plot_name,
1017 'plot': plot_name,
1015 'code': self.exp_code,
1018 'code': self.exp_code,
1016 'time': float(tm),
1019 'time': float(tm),
1017 'data': data,
1020 'data': data,
1018 }
1021 }
1019 meta['type'] = plot_type
1022 meta['type'] = plot_type
1020 meta['interval'] = float(self.interval)
1023 meta['interval'] = float(self.interval)
1021 meta['localtime'] = self.localtime
1024 meta['localtime'] = self.localtime
1022 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1025 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1023 meta.update(self.meta)
1026 meta.update(self.meta)
1024 ret['metadata'] = meta
1027 ret['metadata'] = meta
1025 return json.dumps(ret)
1028 return json.dumps(ret)
1026
1029
1027 @property
1030 @property
1028 def times(self):
1031 def times(self):
1029 '''
1032 '''
1030 Return the list of times of the current data
1033 Return the list of times of the current data
1031 '''
1034 '''
1032
1035
1033 ret = [t for t in self.data]
1036 ret = [t for t in self.data]
1034 ret.sort()
1037 ret.sort()
1035 return numpy.array(ret)
1038 return numpy.array(ret)
1036
1039
1037 @property
1040 @property
1038 def min_time(self):
1041 def min_time(self):
1039 '''
1042 '''
1040 Return the minimun time value
1043 Return the minimun time value
1041 '''
1044 '''
1042
1045
1043 return self.times[0]
1046 return self.times[0]
1044
1047
1045 @property
1048 @property
1046 def max_time(self):
1049 def max_time(self):
1047 '''
1050 '''
1048 Return the maximun time value
1051 Return the maximun time value
1049 '''
1052 '''
1050
1053
1051 return self.times[-1]
1054 return self.times[-1]
1052
1055
1053 # @property
1056 # @property
1054 # def heights(self):
1057 # def heights(self):
1055 # '''
1058 # '''
1056 # Return the list of heights of the current data
1059 # Return the list of heights of the current data
1057 # '''
1060 # '''
1058
1061
1059 # return numpy.array(self.__heights[-1])
1062 # return numpy.array(self.__heights[-1])
1060
1063
1061 @staticmethod
1064 @staticmethod
1062 def roundFloats(obj):
1065 def roundFloats(obj):
1063 if isinstance(obj, list):
1066 if isinstance(obj, list):
1064 return list(map(PlotterData.roundFloats, obj))
1067 return list(map(PlotterData.roundFloats, obj))
1065 elif isinstance(obj, float):
1068 elif isinstance(obj, float):
1066 return round(obj, 2)
1069 return round(obj, 2)
General Comments 0
You need to be logged in to leave comments. Login now