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