##// END OF EJS Templates
schain-jc
RSS
Add C extension for hildebrand_sekhon noise estimation and update requirements versions in setup
Juan C. Valdez -
r878:0ef791f9f681
parent child
Show More
Show file before | Show file after | Hide comments
@@ -0,0 +1,55
1 #include <Python.h>
2 #include <numpy/arrayobject.h>
3 #include <math.h>
4
5 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args);
6
7 static PyMethodDef extensionsMethods[] = {
8 { "hildebrand_sekhon", (PyCFunction)hildebrand_sekhon, METH_VARARGS, "get noise with" },
9 { NULL, NULL, 0, NULL }
10 };
11
12 PyMODINIT_FUNC initcSchain() {
13 Py_InitModule("cSchain", extensionsMethods);
14 import_array();
15 }
16
17 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) {
18 /* Do your stuff here. */
19 double navg;
20 PyObject *data_obj, *data_array;
21
22 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) return NULL;
23 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
24 if (data_array == NULL) {
25 Py_XDECREF(data_array);
26 Py_XDECREF(data_obj);
27 return NULL;
28 }
29 double *sortdata = (double*)PyArray_DATA(data_array);
30 int lenOfData = (int)PyArray_SIZE(data_array) ;
31 double nums_min = lenOfData*0.2;
32 if (nums_min <= 5) nums_min = 5;
33 double sump = 0;
34 double sumq = 0;
35 int j = 0;
36 int cont = 1;
37 double rtest = 0;
38 while ((cont == 1) && (j < lenOfData)) {
39 sump = sump + sortdata[j];
40 sumq = sumq + pow(sortdata[j], 2);
41 if (j > nums_min) {
42 rtest = (double)j/(j-1) + 1/navg;
43 if ((sumq*j) > (rtest*pow(sump, 2))) {
44 j = j - 1;
45 sump = sump - sortdata[j];
46 sumq = sumq - pow(sortdata[j],2);
47 cont = 0;
48 }
49 }
50 j = j + 1;
51 }
52
53 double lnoise = sump / j;
54 return Py_BuildValue("d", lnoise);
55 }
Show file before | Show file after | Hide comments
@@ -1,1183 +1,1189
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10
11 11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 from schainpy import cSchain
13
12 14
13 15 def getNumpyDtype(dataTypeCode):
14 16
15 17 if dataTypeCode == 0:
16 18 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
17 19 elif dataTypeCode == 1:
18 20 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
19 21 elif dataTypeCode == 2:
20 22 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
21 23 elif dataTypeCode == 3:
22 24 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
23 25 elif dataTypeCode == 4:
24 26 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
25 27 elif dataTypeCode == 5:
26 28 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
27 29 else:
28 30 raise ValueError, 'dataTypeCode was not defined'
29 31
30 32 return numpyDtype
31 33
32 34 def getDataTypeCode(numpyDtype):
33 35
34 36 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
35 37 datatype = 0
36 38 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
37 39 datatype = 1
38 40 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
39 41 datatype = 2
40 42 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
41 43 datatype = 3
42 44 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
43 45 datatype = 4
44 46 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
45 47 datatype = 5
46 48 else:
47 49 datatype = None
48 50
49 51 return datatype
50 52
51 53 def hildebrand_sekhon(data, navg):
52 54 """
53 55 This method is for the objective determination of the noise level in Doppler spectra. This
54 56 implementation technique is based on the fact that the standard deviation of the spectral
55 57 densities is equal to the mean spectral density for white Gaussian noise
56 58
57 59 Inputs:
58 60 Data : heights
59 61 navg : numbers of averages
60 62
61 63 Return:
62 64 -1 : any error
63 65 anoise : noise's level
64 66 """
65 67
66 68 sortdata = numpy.sort(data,axis=None)
67 lenOfData = len(sortdata)
68 nums_min = lenOfData*0.2
69
70 if nums_min <= 5:
71 nums_min = 5
69 # lenOfData = len(sortdata)
70 # nums_min = lenOfData*0.2
71 #
72 # if nums_min <= 5:
73 # nums_min = 5
74 #
75 # sump = 0.
76 #
77 # sumq = 0.
78 #
79 # j = 0
80 #
81 # cont = 1
82 #
83 # while((cont==1)and(j<lenOfData)):
84 #
85 # sump += sortdata[j]
86 #
87 # sumq += sortdata[j]**2
88 #
89 # if j > nums_min:
90 # rtest = float(j)/(j-1) + 1.0/navg
91 # if ((sumq*j) > (rtest*sump**2)):
92 # j = j - 1
93 # sump = sump - sortdata[j]
94 # sumq = sumq - sortdata[j]**2
95 # cont = 0
96 #
97 # j += 1
98 #
99 # lnoise = sump /j
100 #
101 # return lnoise
102
103 return cSchain.hildebrand_sekhon(sortdata, navg)
72 104
73 sump = 0.
74
75 sumq = 0.
76
77 j = 0
78
79 cont = 1
80
81 while((cont==1)and(j<lenOfData)):
82
83 sump += sortdata[j]
84
85 sumq += sortdata[j]**2
86
87 if j > nums_min:
88 rtest = float(j)/(j-1) + 1.0/navg
89 if ((sumq*j) > (rtest*sump**2)):
90 j = j - 1
91 sump = sump - sortdata[j]
92 sumq = sumq - sortdata[j]**2
93 cont = 0
94
95 j += 1
96
97 lnoise = sump /j
98 # stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
99 return lnoise
100 105
101 106 class Beam:
107
102 108 def __init__(self):
103 109 self.codeList = []
104 110 self.azimuthList = []
105 111 self.zenithList = []
106 112
107 113 class GenericData(object):
108 114
109 115 flagNoData = True
110 116
111 117 def __init__(self):
112 118
113 119 raise NotImplementedError
114 120
115 121 def copy(self, inputObj=None):
116 122
117 123 if inputObj == None:
118 124 return copy.deepcopy(self)
119 125
120 126 for key in inputObj.__dict__.keys():
121 127
122 128 attribute = inputObj.__dict__[key]
123 129
124 130 #If this attribute is a tuple or list
125 131 if type(inputObj.__dict__[key]) in (tuple, list):
126 132 self.__dict__[key] = attribute[:]
127 133 continue
128 134
129 135 #If this attribute is another object or instance
130 136 if hasattr(attribute, '__dict__'):
131 137 self.__dict__[key] = attribute.copy()
132 138 continue
133 139
134 140 self.__dict__[key] = inputObj.__dict__[key]
135 141
136 142 def deepcopy(self):
137 143
138 144 return copy.deepcopy(self)
139 145
140 146 def isEmpty(self):
141 147
142 148 return self.flagNoData
143 149
144 150 class JROData(GenericData):
145 151
146 152 # m_BasicHeader = BasicHeader()
147 153 # m_ProcessingHeader = ProcessingHeader()
148 154
149 155 systemHeaderObj = SystemHeader()
150 156
151 157 radarControllerHeaderObj = RadarControllerHeader()
152 158
153 159 # data = None
154 160
155 161 type = None
156 162
157 163 datatype = None #dtype but in string
158 164
159 165 # dtype = None
160 166
161 167 # nChannels = None
162 168
163 169 # nHeights = None
164 170
165 171 nProfiles = None
166 172
167 173 heightList = None
168 174
169 175 channelList = None
170 176
171 177 flagDiscontinuousBlock = False
172 178
173 179 useLocalTime = False
174 180
175 181 utctime = None
176 182
177 183 timeZone = None
178 184
179 185 dstFlag = None
180 186
181 187 errorCount = None
182 188
183 189 blocksize = None
184 190
185 191 # nCode = None
186 192 #
187 193 # nBaud = None
188 194 #
189 195 # code = None
190 196
191 197 flagDecodeData = False #asumo q la data no esta decodificada
192 198
193 199 flagDeflipData = False #asumo q la data no esta sin flip
194 200
195 201 flagShiftFFT = False
196 202
197 203 # ippSeconds = None
198 204
199 205 # timeInterval = None
200 206
201 207 nCohInt = None
202 208
203 209 # noise = None
204 210
205 211 windowOfFilter = 1
206 212
207 213 #Speed of ligth
208 214 C = 3e8
209 215
210 216 frequency = 49.92e6
211 217
212 218 realtime = False
213 219
214 220 beacon_heiIndexList = None
215 221
216 222 last_block = None
217 223
218 224 blocknow = None
219 225
220 226 azimuth = None
221 227
222 228 zenith = None
223 229
224 230 beam = Beam()
225 231
226 232 profileIndex = None
227 233
228 234 def __init__(self):
229 235
230 236 raise NotImplementedError
231 237
232 238 def getNoise(self):
233 239
234 240 raise NotImplementedError
235 241
236 242 def getNChannels(self):
237 243
238 244 return len(self.channelList)
239 245
240 246 def getChannelIndexList(self):
241 247
242 248 return range(self.nChannels)
243 249
244 250 def getNHeights(self):
245 251
246 252 return len(self.heightList)
247 253
248 254 def getHeiRange(self, extrapoints=0):
249 255
250 256 heis = self.heightList
251 257 # deltah = self.heightList[1] - self.heightList[0]
252 258 #
253 259 # heis.append(self.heightList[-1])
254 260
255 261 return heis
256 262
257 263 def getDeltaH(self):
258 264
259 265 delta = self.heightList[1] - self.heightList[0]
260 266
261 267 return delta
262 268
263 269 def getltctime(self):
264 270
265 271 if self.useLocalTime:
266 272 return self.utctime - self.timeZone*60
267 273
268 274 return self.utctime
269 275
270 276 def getDatatime(self):
271 277
272 278 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
273 279 return datatimeValue
274 280
275 281 def getTimeRange(self):
276 282
277 283 datatime = []
278 284
279 285 datatime.append(self.ltctime)
280 286 datatime.append(self.ltctime + self.timeInterval+1)
281 287
282 288 datatime = numpy.array(datatime)
283 289
284 290 return datatime
285 291
286 292 def getFmaxTimeResponse(self):
287 293
288 294 period = (10**-6)*self.getDeltaH()/(0.15)
289 295
290 296 PRF = 1./(period * self.nCohInt)
291 297
292 298 fmax = PRF
293 299
294 300 return fmax
295 301
296 302 def getFmax(self):
297 303
298 304 PRF = 1./(self.ippSeconds * self.nCohInt)
299 305
300 306 fmax = PRF
301 307
302 308 return fmax
303 309
304 310 def getVmax(self):
305 311
306 312 _lambda = self.C/self.frequency
307 313
308 314 vmax = self.getFmax() * _lambda/2
309 315
310 316 return vmax
311 317
312 318 def get_ippSeconds(self):
313 319 '''
314 320 '''
315 321 return self.radarControllerHeaderObj.ippSeconds
316 322
317 323 def set_ippSeconds(self, ippSeconds):
318 324 '''
319 325 '''
320 326
321 327 self.radarControllerHeaderObj.ippSeconds = ippSeconds
322 328
323 329 return
324 330
325 331 def get_dtype(self):
326 332 '''
327 333 '''
328 334 return getNumpyDtype(self.datatype)
329 335
330 336 def set_dtype(self, numpyDtype):
331 337 '''
332 338 '''
333 339
334 340 self.datatype = getDataTypeCode(numpyDtype)
335 341
336 342 def get_code(self):
337 343 '''
338 344 '''
339 345 return self.radarControllerHeaderObj.code
340 346
341 347 def set_code(self, code):
342 348 '''
343 349 '''
344 350 self.radarControllerHeaderObj.code = code
345 351
346 352 return
347 353
348 354 def get_ncode(self):
349 355 '''
350 356 '''
351 357 return self.radarControllerHeaderObj.nCode
352 358
353 359 def set_ncode(self, nCode):
354 360 '''
355 361 '''
356 362 self.radarControllerHeaderObj.nCode = nCode
357 363
358 364 return
359 365
360 366 def get_nbaud(self):
361 367 '''
362 368 '''
363 369 return self.radarControllerHeaderObj.nBaud
364 370
365 371 def set_nbaud(self, nBaud):
366 372 '''
367 373 '''
368 374 self.radarControllerHeaderObj.nBaud = nBaud
369 375
370 376 return
371 377
372 378 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
373 379 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
374 380 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
375 381 #noise = property(getNoise, "I'm the 'nHeights' property.")
376 382 datatime = property(getDatatime, "I'm the 'datatime' property")
377 383 ltctime = property(getltctime, "I'm the 'ltctime' property")
378 384 ippSeconds = property(get_ippSeconds, set_ippSeconds)
379 385 dtype = property(get_dtype, set_dtype)
380 386 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
381 387 code = property(get_code, set_code)
382 388 nCode = property(get_ncode, set_ncode)
383 389 nBaud = property(get_nbaud, set_nbaud)
384 390
385 391 class Voltage(JROData):
386 392
387 393 #data es un numpy array de 2 dmensiones (canales, alturas)
388 394 data = None
389 395
390 396 def __init__(self):
391 397 '''
392 398 Constructor
393 399 '''
394 400
395 401 self.useLocalTime = True
396 402
397 403 self.radarControllerHeaderObj = RadarControllerHeader()
398 404
399 405 self.systemHeaderObj = SystemHeader()
400 406
401 407 self.type = "Voltage"
402 408
403 409 self.data = None
404 410
405 411 # self.dtype = None
406 412
407 413 # self.nChannels = 0
408 414
409 415 # self.nHeights = 0
410 416
411 417 self.nProfiles = None
412 418
413 419 self.heightList = None
414 420
415 421 self.channelList = None
416 422
417 423 # self.channelIndexList = None
418 424
419 425 self.flagNoData = True
420 426
421 427 self.flagDiscontinuousBlock = False
422 428
423 429 self.utctime = None
424 430
425 431 self.timeZone = None
426 432
427 433 self.dstFlag = None
428 434
429 435 self.errorCount = None
430 436
431 437 self.nCohInt = None
432 438
433 439 self.blocksize = None
434 440
435 441 self.flagDecodeData = False #asumo q la data no esta decodificada
436 442
437 443 self.flagDeflipData = False #asumo q la data no esta sin flip
438 444
439 445 self.flagShiftFFT = False
440 446
441 447 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
442 448
443 449 self.profileIndex = 0
444 450
445 451 def getNoisebyHildebrand(self, channel = None):
446 452 """
447 453 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
448 454
449 455 Return:
450 456 noiselevel
451 457 """
452 458
453 459 if channel != None:
454 460 data = self.data[channel]
455 461 nChannels = 1
456 462 else:
457 463 data = self.data
458 464 nChannels = self.nChannels
459 465
460 466 noise = numpy.zeros(nChannels)
461 467 power = data * numpy.conjugate(data)
462 468
463 469 for thisChannel in range(nChannels):
464 470 if nChannels == 1:
465 471 daux = power[:].real
466 472 else:
467 473 daux = power[thisChannel,:].real
468 474 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
469 475
470 476 return noise
471 477
472 478 def getNoise(self, type = 1, channel = None):
473 479
474 480 if type == 1:
475 481 noise = self.getNoisebyHildebrand(channel)
476 482
477 483 return noise
478 484
479 485 def getPower(self, channel = None):
480 486
481 487 if channel != None:
482 488 data = self.data[channel]
483 489 else:
484 490 data = self.data
485 491
486 492 power = data * numpy.conjugate(data)
487 493 powerdB = 10*numpy.log10(power.real)
488 494 powerdB = numpy.squeeze(powerdB)
489 495
490 496 return powerdB
491 497
492 498 def getTimeInterval(self):
493 499
494 500 timeInterval = self.ippSeconds * self.nCohInt
495 501
496 502 return timeInterval
497 503
498 504 noise = property(getNoise, "I'm the 'nHeights' property.")
499 505 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
500 506
501 507 class Spectra(JROData):
502 508
503 509 #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
504 510 data_spc = None
505 511
506 512 #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
507 513 data_cspc = None
508 514
509 515 #data dc es un numpy array de 2 dmensiones (canales, alturas)
510 516 data_dc = None
511 517
512 518 #data power
513 519 data_pwr = None
514 520
515 521 nFFTPoints = None
516 522
517 523 # nPairs = None
518 524
519 525 pairsList = None
520 526
521 527 nIncohInt = None
522 528
523 529 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
524 530
525 531 nCohInt = None #se requiere para determinar el valor de timeInterval
526 532
527 533 ippFactor = None
528 534
529 535 profileIndex = 0
530 536
531 537 plotting = "spectra"
532 538
533 539 def __init__(self):
534 540 '''
535 541 Constructor
536 542 '''
537 543
538 544 self.useLocalTime = True
539 545
540 546 self.radarControllerHeaderObj = RadarControllerHeader()
541 547
542 548 self.systemHeaderObj = SystemHeader()
543 549
544 550 self.type = "Spectra"
545 551
546 552 # self.data = None
547 553
548 554 # self.dtype = None
549 555
550 556 # self.nChannels = 0
551 557
552 558 # self.nHeights = 0
553 559
554 560 self.nProfiles = None
555 561
556 562 self.heightList = None
557 563
558 564 self.channelList = None
559 565
560 566 # self.channelIndexList = None
561 567
562 568 self.pairsList = None
563 569
564 570 self.flagNoData = True
565 571
566 572 self.flagDiscontinuousBlock = False
567 573
568 574 self.utctime = None
569 575
570 576 self.nCohInt = None
571 577
572 578 self.nIncohInt = None
573 579
574 580 self.blocksize = None
575 581
576 582 self.nFFTPoints = None
577 583
578 584 self.wavelength = None
579 585
580 586 self.flagDecodeData = False #asumo q la data no esta decodificada
581 587
582 588 self.flagDeflipData = False #asumo q la data no esta sin flip
583 589
584 590 self.flagShiftFFT = False
585 591
586 592 self.ippFactor = 1
587 593
588 594 #self.noise = None
589 595
590 596 self.beacon_heiIndexList = []
591 597
592 598 self.noise_estimation = None
593 599
594 600
595 601 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
596 602 """
597 603 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
598 604
599 605 Return:
600 606 noiselevel
601 607 """
602 608
603 609 noise = numpy.zeros(self.nChannels)
604 610
605 611 for channel in range(self.nChannels):
606 612 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
607 613 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
608 614
609 615 return noise
610 616
611 617 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
612 618
613 619 if self.noise_estimation is not None:
614 620 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
615 621 else:
616 622 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
617 623 return noise
618 624
619 625 def getFreqRangeTimeResponse(self, extrapoints=0):
620 626
621 627 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor)
622 628 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
623 629
624 630 return freqrange
625 631
626 632 def getAcfRange(self, extrapoints=0):
627 633
628 634 deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor))
629 635 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
630 636
631 637 return freqrange
632 638
633 639 def getFreqRange(self, extrapoints=0):
634 640
635 641 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
636 642 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
637 643
638 644 return freqrange
639 645
640 646 def getVelRange(self, extrapoints=0):
641 647
642 648 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
643 649 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2
644 650
645 651 return velrange
646 652
647 653 def getNPairs(self):
648 654
649 655 return len(self.pairsList)
650 656
651 657 def getPairsIndexList(self):
652 658
653 659 return range(self.nPairs)
654 660
655 661 def getNormFactor(self):
656 662
657 663 pwcode = 1
658 664
659 665 if self.flagDecodeData:
660 666 pwcode = numpy.sum(self.code[0]**2)
661 667 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
662 668 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
663 669
664 670 return normFactor
665 671
666 672 def getFlagCspc(self):
667 673
668 674 if self.data_cspc is None:
669 675 return True
670 676
671 677 return False
672 678
673 679 def getFlagDc(self):
674 680
675 681 if self.data_dc is None:
676 682 return True
677 683
678 684 return False
679 685
680 686 def getTimeInterval(self):
681 687
682 688 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
683 689
684 690 return timeInterval
685 691
686 692 def getPower(self):
687 693
688 694 factor = self.normFactor
689 695 z = self.data_spc/factor
690 696 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
691 697 avg = numpy.average(z, axis=1)
692 698
693 699 return 10*numpy.log10(avg)
694 700
695 701 def setValue(self, value):
696 702
697 703 print "This property should not be initialized"
698 704
699 705 return
700 706
701 707 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
702 708 pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
703 709 normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.")
704 710 flag_cspc = property(getFlagCspc, setValue)
705 711 flag_dc = property(getFlagDc, setValue)
706 712 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
707 713 timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property")
708 714
709 715 class SpectraHeis(Spectra):
710 716
711 717 data_spc = None
712 718
713 719 data_cspc = None
714 720
715 721 data_dc = None
716 722
717 723 nFFTPoints = None
718 724
719 725 # nPairs = None
720 726
721 727 pairsList = None
722 728
723 729 nCohInt = None
724 730
725 731 nIncohInt = None
726 732
727 733 def __init__(self):
728 734
729 735 self.radarControllerHeaderObj = RadarControllerHeader()
730 736
731 737 self.systemHeaderObj = SystemHeader()
732 738
733 739 self.type = "SpectraHeis"
734 740
735 741 # self.dtype = None
736 742
737 743 # self.nChannels = 0
738 744
739 745 # self.nHeights = 0
740 746
741 747 self.nProfiles = None
742 748
743 749 self.heightList = None
744 750
745 751 self.channelList = None
746 752
747 753 # self.channelIndexList = None
748 754
749 755 self.flagNoData = True
750 756
751 757 self.flagDiscontinuousBlock = False
752 758
753 759 # self.nPairs = 0
754 760
755 761 self.utctime = None
756 762
757 763 self.blocksize = None
758 764
759 765 self.profileIndex = 0
760 766
761 767 self.nCohInt = 1
762 768
763 769 self.nIncohInt = 1
764 770
765 771 def getNormFactor(self):
766 772 pwcode = 1
767 773 if self.flagDecodeData:
768 774 pwcode = numpy.sum(self.code[0]**2)
769 775
770 776 normFactor = self.nIncohInt*self.nCohInt*pwcode
771 777
772 778 return normFactor
773 779
774 780 def getTimeInterval(self):
775 781
776 782 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
777 783
778 784 return timeInterval
779 785
780 786 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
781 787 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
782 788
783 789 class Fits(JROData):
784 790
785 791 heightList = None
786 792
787 793 channelList = None
788 794
789 795 flagNoData = True
790 796
791 797 flagDiscontinuousBlock = False
792 798
793 799 useLocalTime = False
794 800
795 801 utctime = None
796 802
797 803 timeZone = None
798 804
799 805 # ippSeconds = None
800 806
801 807 # timeInterval = None
802 808
803 809 nCohInt = None
804 810
805 811 nIncohInt = None
806 812
807 813 noise = None
808 814
809 815 windowOfFilter = 1
810 816
811 817 #Speed of ligth
812 818 C = 3e8
813 819
814 820 frequency = 49.92e6
815 821
816 822 realtime = False
817 823
818 824
819 825 def __init__(self):
820 826
821 827 self.type = "Fits"
822 828
823 829 self.nProfiles = None
824 830
825 831 self.heightList = None
826 832
827 833 self.channelList = None
828 834
829 835 # self.channelIndexList = None
830 836
831 837 self.flagNoData = True
832 838
833 839 self.utctime = None
834 840
835 841 self.nCohInt = 1
836 842
837 843 self.nIncohInt = 1
838 844
839 845 self.useLocalTime = True
840 846
841 847 self.profileIndex = 0
842 848
843 849 # self.utctime = None
844 850 # self.timeZone = None
845 851 # self.ltctime = None
846 852 # self.timeInterval = None
847 853 # self.header = None
848 854 # self.data_header = None
849 855 # self.data = None
850 856 # self.datatime = None
851 857 # self.flagNoData = False
852 858 # self.expName = ''
853 859 # self.nChannels = None
854 860 # self.nSamples = None
855 861 # self.dataBlocksPerFile = None
856 862 # self.comments = ''
857 863 #
858 864
859 865
860 866 def getltctime(self):
861 867
862 868 if self.useLocalTime:
863 869 return self.utctime - self.timeZone*60
864 870
865 871 return self.utctime
866 872
867 873 def getDatatime(self):
868 874
869 875 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
870 876 return datatime
871 877
872 878 def getTimeRange(self):
873 879
874 880 datatime = []
875 881
876 882 datatime.append(self.ltctime)
877 883 datatime.append(self.ltctime + self.timeInterval)
878 884
879 885 datatime = numpy.array(datatime)
880 886
881 887 return datatime
882 888
883 889 def getHeiRange(self):
884 890
885 891 heis = self.heightList
886 892
887 893 return heis
888 894
889 895 def getNHeights(self):
890 896
891 897 return len(self.heightList)
892 898
893 899 def getNChannels(self):
894 900
895 901 return len(self.channelList)
896 902
897 903 def getChannelIndexList(self):
898 904
899 905 return range(self.nChannels)
900 906
901 907 def getNoise(self, type = 1):
902 908
903 909 #noise = numpy.zeros(self.nChannels)
904 910
905 911 if type == 1:
906 912 noise = self.getNoisebyHildebrand()
907 913
908 914 if type == 2:
909 915 noise = self.getNoisebySort()
910 916
911 917 if type == 3:
912 918 noise = self.getNoisebyWindow()
913 919
914 920 return noise
915 921
916 922 def getTimeInterval(self):
917 923
918 924 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
919 925
920 926 return timeInterval
921 927
922 928 datatime = property(getDatatime, "I'm the 'datatime' property")
923 929 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
924 930 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
925 931 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
926 932 noise = property(getNoise, "I'm the 'nHeights' property.")
927 933
928 934 ltctime = property(getltctime, "I'm the 'ltctime' property")
929 935 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
930 936
931 937
932 938 class Correlation(JROData):
933 939
934 940 noise = None
935 941
936 942 SNR = None
937 943
938 944 #--------------------------------------------------
939 945
940 946 mode = None
941 947
942 948 split = False
943 949
944 950 data_cf = None
945 951
946 952 lags = None
947 953
948 954 lagRange = None
949 955
950 956 pairsList = None
951 957
952 958 normFactor = None
953 959
954 960 #--------------------------------------------------
955 961
956 962 # calculateVelocity = None
957 963
958 964 nLags = None
959 965
960 966 nPairs = None
961 967
962 968 nAvg = None
963 969
964 970
965 971 def __init__(self):
966 972 '''
967 973 Constructor
968 974 '''
969 975 self.radarControllerHeaderObj = RadarControllerHeader()
970 976
971 977 self.systemHeaderObj = SystemHeader()
972 978
973 979 self.type = "Correlation"
974 980
975 981 self.data = None
976 982
977 983 self.dtype = None
978 984
979 985 self.nProfiles = None
980 986
981 987 self.heightList = None
982 988
983 989 self.channelList = None
984 990
985 991 self.flagNoData = True
986 992
987 993 self.flagDiscontinuousBlock = False
988 994
989 995 self.utctime = None
990 996
991 997 self.timeZone = None
992 998
993 999 self.dstFlag = None
994 1000
995 1001 self.errorCount = None
996 1002
997 1003 self.blocksize = None
998 1004
999 1005 self.flagDecodeData = False #asumo q la data no esta decodificada
1000 1006
1001 1007 self.flagDeflipData = False #asumo q la data no esta sin flip
1002 1008
1003 1009 self.pairsList = None
1004 1010
1005 1011 self.nPoints = None
1006 1012
1007 1013 def getPairsList(self):
1008 1014
1009 1015 return self.pairsList
1010 1016
1011 1017 def getNoise(self, mode = 2):
1012 1018
1013 1019 indR = numpy.where(self.lagR == 0)[0][0]
1014 1020 indT = numpy.where(self.lagT == 0)[0][0]
1015 1021
1016 1022 jspectra0 = self.data_corr[:,:,indR,:]
1017 1023 jspectra = copy.copy(jspectra0)
1018 1024
1019 1025 num_chan = jspectra.shape[0]
1020 1026 num_hei = jspectra.shape[2]
1021 1027
1022 1028 freq_dc = jspectra.shape[1]/2
1023 1029 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1024 1030
1025 1031 if ind_vel[0]<0:
1026 1032 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1027 1033
1028 1034 if mode == 1:
1029 1035 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1030 1036
1031 1037 if mode == 2:
1032 1038
1033 1039 vel = numpy.array([-2,-1,1,2])
1034 1040 xx = numpy.zeros([4,4])
1035 1041
1036 1042 for fil in range(4):
1037 1043 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1038 1044
1039 1045 xx_inv = numpy.linalg.inv(xx)
1040 1046 xx_aux = xx_inv[0,:]
1041 1047
1042 1048 for ich in range(num_chan):
1043 1049 yy = jspectra[ich,ind_vel,:]
1044 1050 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1045 1051
1046 1052 junkid = jspectra[ich,freq_dc,:]<=0
1047 1053 cjunkid = sum(junkid)
1048 1054
1049 1055 if cjunkid.any():
1050 1056 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1051 1057
1052 1058 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1053 1059
1054 1060 return noise
1055 1061
1056 1062 def getTimeInterval(self):
1057 1063
1058 1064 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
1059 1065
1060 1066 return timeInterval
1061 1067
1062 1068 def splitFunctions(self):
1063 1069
1064 1070 pairsList = self.pairsList
1065 1071 ccf_pairs = []
1066 1072 acf_pairs = []
1067 1073 ccf_ind = []
1068 1074 acf_ind = []
1069 1075 for l in range(len(pairsList)):
1070 1076 chan0 = pairsList[l][0]
1071 1077 chan1 = pairsList[l][1]
1072 1078
1073 1079 #Obteniendo pares de Autocorrelacion
1074 1080 if chan0 == chan1:
1075 1081 acf_pairs.append(chan0)
1076 1082 acf_ind.append(l)
1077 1083 else:
1078 1084 ccf_pairs.append(pairsList[l])
1079 1085 ccf_ind.append(l)
1080 1086
1081 1087 data_acf = self.data_cf[acf_ind]
1082 1088 data_ccf = self.data_cf[ccf_ind]
1083 1089
1084 1090 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1085 1091
1086 1092 def getNormFactor(self):
1087 1093 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1088 1094 acf_pairs = numpy.array(acf_pairs)
1089 1095 normFactor = numpy.zeros((self.nPairs,self.nHeights))
1090 1096
1091 1097 for p in range(self.nPairs):
1092 1098 pair = self.pairsList[p]
1093 1099
1094 1100 ch0 = pair[0]
1095 1101 ch1 = pair[1]
1096 1102
1097 1103 ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1)
1098 1104 ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1)
1099 1105 normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max)
1100 1106
1101 1107 return normFactor
1102 1108
1103 1109 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1104 1110 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1105 1111
1106 1112 class Parameters(JROData):
1107 1113
1108 1114 experimentInfo = None #Information about the experiment
1109 1115
1110 1116 #Information from previous data
1111 1117
1112 1118 inputUnit = None #Type of data to be processed
1113 1119
1114 1120 operation = None #Type of operation to parametrize
1115 1121
1116 1122 normFactor = None #Normalization Factor
1117 1123
1118 1124 groupList = None #List of Pairs, Groups, etc
1119 1125
1120 1126 #Parameters
1121 1127
1122 1128 data_param = None #Parameters obtained
1123 1129
1124 1130 data_pre = None #Data Pre Parametrization
1125 1131
1126 1132 data_SNR = None #Signal to Noise Ratio
1127 1133
1128 1134 # heightRange = None #Heights
1129 1135
1130 1136 abscissaList = None #Abscissa, can be velocities, lags or time
1131 1137
1132 1138 noise = None #Noise Potency
1133 1139
1134 1140 utctimeInit = None #Initial UTC time
1135 1141
1136 1142 paramInterval = None #Time interval to calculate Parameters in seconds
1137 1143
1138 1144 useLocalTime = True
1139 1145
1140 1146 #Fitting
1141 1147
1142 1148 data_error = None #Error of the estimation
1143 1149
1144 1150 constants = None
1145 1151
1146 1152 library = None
1147 1153
1148 1154 #Output signal
1149 1155
1150 1156 outputInterval = None #Time interval to calculate output signal in seconds
1151 1157
1152 1158 data_output = None #Out signal
1153 1159
1154 1160 nAvg = None
1155 1161
1156 1162
1157 1163 def __init__(self):
1158 1164 '''
1159 1165 Constructor
1160 1166 '''
1161 1167 self.radarControllerHeaderObj = RadarControllerHeader()
1162 1168
1163 1169 self.systemHeaderObj = SystemHeader()
1164 1170
1165 1171 self.type = "Parameters"
1166 1172
1167 1173 def getTimeRange1(self, interval):
1168 1174
1169 1175 datatime = []
1170 1176
1171 1177 if self.useLocalTime:
1172 1178 time1 = self.utctimeInit - self.timeZone*60
1173 1179 else:
1174 1180 time1 = self.utctimeInit
1175 1181
1176 1182 # datatime.append(self.utctimeInit)
1177 1183 # datatime.append(self.utctimeInit + self.outputInterval)
1178 1184 datatime.append(time1)
1179 1185 datatime.append(time1 + interval)
1180 1186
1181 1187 datatime = numpy.array(datatime)
1182 1188
1183 1189 return datatime
Show file before | Show file after | Hide comments
@@ -1,45 +1,48
1 1 '''
2 2 Created on Jul 16, 2014
3 3
4 @author: @author: Miguel Urco
4 @author: Miguel Urco
5 5 '''
6 6
7 7 from schainpy import __version__
8 8 from setuptools import setup, Extension
9 9
10 10 setup(name="schainpy",
11 11 version=__version__,
12 12 description="Python tools to read, write and process Jicamarca data",
13 13 author="Miguel Urco",
14 14 author_email="miguel.urco@jro.igp.gob.pe",
15 15 url="http://jro.igp.gob.pe",
16 16 packages = {'schainpy',
17 17 'schainpy.model',
18 18 'schainpy.model.data',
19 19 'schainpy.model.graphics',
20 20 'schainpy.model.io',
21 21 'schainpy.model.proc',
22 22 'schainpy.model.serializer',
23 23 'schainpy.model.utils',
24 24 'schainpy.gui',
25 25 'schainpy.gui.figures',
26 26 'schainpy.gui.viewcontroller',
27 27 'schainpy.gui.viewer',
28 28 'schainpy.gui.viewer.windows'},
29 ext_package='schainpy',
29 30 py_modules=[''],
30 31 package_data={'': ['schain.conf.template'],
31 32 'schainpy.gui.figures': ['*.png','*.jpg'],
32 33 },
33 34 include_package_data=False,
34 35 scripts =['schainpy/gui/schainGUI',
35 36 'schainpy/scripts/schain'],
37 ext_modules=[Extension("cSchain", ["schainpy/model/proc/extensions.c"])],
36 38 install_requires=[
37 "scipy >= 0.9.0",
38 "h5py >= 2.0.1",
39 "matplotlib >= 1.0.0",
40 "pyfits >= 2.0.0",
41 "numpy >= 1.6.0",
42 "paramiko",
43 "paho-mqtt"
39 "scipy >= 0.14.0",
40 "h5py >= 2.2.1",
41 "matplotlib >= 1.4.2",
42 "pyfits >= 3.4",
43 "numpy >= 1.11.2",
44 "paramiko >= 2.1.2",
45 "paho-mqtt >= 1.2",
46 "zmq",
44 47 ],
45 48 ) No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now