##// END OF EJS Templates
schain-jc
RSS
Getnoise se agrega al metodo getFft de SpectraProcessor.
Daniel Valdez -
r111:9c1626ca5c6a
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,963 +1,965
1 1 """
2 2 Created on Feb 7, 2012
3 3
4 4 @autor $Author$
5 5 @version $Id$
6 6
7 7 """
8 8 import os
9 9 import numpy
10 10 import sys
11 11 import time
12 12 import datetime
13 13 import time
14 14 import plplot
15 15
16 16
17 17 def cmap1_init(colormap="gray"):
18 18
19 19 if colormap == None:
20 20 return
21 21
22 22 ncolor = None
23 23 rgb_lvl = None
24 24
25 25 # Routine for defining a specific color map 1 in HLS space.
26 26 # if gray is true, use basic grayscale variation from half-dark to light.
27 27 # otherwise use false color variation from blue (240 deg) to red (360 deg).
28 28
29 29 # Independent variable of control points.
30 30 i = numpy.array((0., 1.))
31 31 if colormap=="gray":
32 32 ncolor = 256
33 33 # Hue for control points. Doesn't matter since saturation is zero.
34 34 h = numpy.array((0., 0.))
35 35 # Lightness ranging from half-dark (for interest) to light.
36 36 l = numpy.array((0.5, 1.))
37 37 # Gray scale has zero saturation
38 38 s = numpy.array((0., 0.))
39 39
40 40 # number of cmap1 colours is 256 in this case.
41 41 plplot.plscmap1n(ncolor)
42 42 # Interpolate between control points to set up cmap1.
43 43 plplot.plscmap1l(0, i, h, l, s)
44 44
45 45 return None
46 46
47 47 if colormap == 'jet':
48 48 ncolor = 256
49 49 pos = numpy.zeros((ncolor))
50 50 r = numpy.zeros((ncolor))
51 51 g = numpy.zeros((ncolor))
52 52 b = numpy.zeros((ncolor))
53 53
54 54 for i in range(ncolor):
55 55 if(i <= 35.0/100*(ncolor-1)): rf = 0.0
56 56 elif (i <= 66.0/100*(ncolor-1)): rf = (100.0/31)*i/(ncolor-1) - 35.0/31
57 57 elif (i <= 89.0/100*(ncolor-1)): rf = 1.0
58 58 else: rf = (-100.0/22)*i/(ncolor-1) + 111.0/22
59 59
60 60 if(i <= 12.0/100*(ncolor-1)): gf = 0.0
61 61 elif(i <= 38.0/100*(ncolor-1)): gf = (100.0/26)*i/(ncolor-1) - 12.0/26
62 62 elif(i <= 64.0/100*(ncolor-1)): gf = 1.0
63 63 elif(i <= 91.0/100*(ncolor-1)): gf = (-100.0/27)*i/(ncolor-1) + 91.0/27
64 64 else: gf = 0.0
65 65
66 66 if(i <= 11.0/100*(ncolor-1)): bf = (50.0/11)*i/(ncolor-1) + 0.5
67 67 elif(i <= 34.0/100*(ncolor-1)): bf = 1.0
68 68 elif(i <= 65.0/100*(ncolor-1)): bf = (-100.0/31)*i/(ncolor-1) + 65.0/31
69 69 else: bf = 0
70 70
71 71 r[i] = rf
72 72 g[i] = gf
73 73 b[i] = bf
74 74
75 75 pos[i] = float(i)/float(ncolor-1)
76 76
77 77
78 78 plplot.plscmap1n(ncolor)
79 79 plplot.plscmap1l(1, pos, r, g, b)
80 80
81 81
82 82
83 83 if colormap=="br_green":
84 84 ncolor = 256
85 85 # Hue ranges from blue (240 deg) to red (0 or 360 deg)
86 86 h = numpy.array((240., 0.))
87 87 # Lightness and saturation are constant (values taken from C example).
88 88 l = numpy.array((0.6, 0.6))
89 89 s = numpy.array((0.8, 0.8))
90 90
91 91 # number of cmap1 colours is 256 in this case.
92 92 plplot.plscmap1n(ncolor)
93 93 # Interpolate between control points to set up cmap1.
94 94 plplot.plscmap1l(0, i, h, l, s)
95 95
96 96 return None
97 97
98 98 if colormap=="tricolor":
99 99 ncolor = 3
100 100 # Hue ranges from blue (240 deg) to red (0 or 360 deg)
101 101 h = numpy.array((240., 0.))
102 102 # Lightness and saturation are constant (values taken from C example).
103 103 l = numpy.array((0.6, 0.6))
104 104 s = numpy.array((0.8, 0.8))
105 105
106 106 # number of cmap1 colours is 256 in this case.
107 107 plplot.plscmap1n(ncolor)
108 108 # Interpolate between control points to set up cmap1.
109 109 plplot.plscmap1l(0, i, h, l, s)
110 110
111 111 return None
112 112
113 113 if colormap == 'rgb' or colormap == 'rgb666':
114 114
115 115 color_sz = 6
116 116 ncolor = color_sz*color_sz*color_sz
117 117 pos = numpy.zeros((ncolor))
118 118 r = numpy.zeros((ncolor))
119 119 g = numpy.zeros((ncolor))
120 120 b = numpy.zeros((ncolor))
121 121 ind = 0
122 122 for ri in range(color_sz):
123 123 for gi in range(color_sz):
124 124 for bi in range(color_sz):
125 125 r[ind] = ri/(color_sz-1.0)
126 126 g[ind] = gi/(color_sz-1.0)
127 127 b[ind] = bi/(color_sz-1.0)
128 128 pos[ind] = ind/(ncolor-1.0)
129 129 ind += 1
130 130 rgb_lvl = [6,6,6] #Levels for RGB colors
131 131
132 132 if colormap == 'rgb676':
133 133 ncolor = 6*7*6
134 134 pos = numpy.zeros((ncolor))
135 135 r = numpy.zeros((ncolor))
136 136 g = numpy.zeros((ncolor))
137 137 b = numpy.zeros((ncolor))
138 138 ind = 0
139 139 for ri in range(8):
140 140 for gi in range(8):
141 141 for bi in range(4):
142 142 r[ind] = ri/(6-1.0)
143 143 g[ind] = gi/(7-1.0)
144 144 b[ind] = bi/(6-1.0)
145 145 pos[ind] = ind/(ncolor-1.0)
146 146 ind += 1
147 147 rgb_lvl = [6,7,6] #Levels for RGB colors
148 148
149 149 if colormap == 'rgb685':
150 150 ncolor = 6*8*5
151 151 pos = numpy.zeros((ncolor))
152 152 r = numpy.zeros((ncolor))
153 153 g = numpy.zeros((ncolor))
154 154 b = numpy.zeros((ncolor))
155 155 ind = 0
156 156 for ri in range(8):
157 157 for gi in range(8):
158 158 for bi in range(4):
159 159 r[ind] = ri/(6-1.0)
160 160 g[ind] = gi/(8-1.0)
161 161 b[ind] = bi/(5-1.0)
162 162 pos[ind] = ind/(ncolor-1.0)
163 163 ind += 1
164 164 rgb_lvl = [6,8,5] #Levels for RGB colors
165 165
166 166 if colormap == 'rgb884':
167 167 ncolor = 8*8*4
168 168 pos = numpy.zeros((ncolor))
169 169 r = numpy.zeros((ncolor))
170 170 g = numpy.zeros((ncolor))
171 171 b = numpy.zeros((ncolor))
172 172 ind = 0
173 173 for ri in range(8):
174 174 for gi in range(8):
175 175 for bi in range(4):
176 176 r[ind] = ri/(8-1.0)
177 177 g[ind] = gi/(8-1.0)
178 178 b[ind] = bi/(4-1.0)
179 179 pos[ind] = ind/(ncolor-1.0)
180 180 ind += 1
181 181 rgb_lvl = [8,8,4] #Levels for RGB colors
182 182
183 183 if ncolor == None:
184 184 raise ValueError, "The colormap selected is not valid"
185 185
186 186 plplot.plscmap1n(ncolor)
187 187 plplot.plscmap1l(1, pos, r, g, b)
188 188
189 189 return rgb_lvl
190 190
191 191 def setColormap(colormap="jet"):
192 192 cmap1_init(colormap)
193 193
194 def savePlplot(indexPlot,filename,ncol,nrow,width,height):
194 def savePlplot(filename,width,height):
195 195 curr_strm = plplot.plgstrm()
196 196 save_strm = plplot.plmkstrm()
197 plplot.plsetopt("geometry", "%dx%d"%(width*ncol,height*nrow))
198 plplot.plsdev("pngcairo")
197 plplot.plsetopt("geometry", "%dx%d"%(width,height))
198 plplot.plsdev("png")
199 199 plplot.plsfnam(filename)
200 200 plplot.plcpstrm(curr_strm,0)
201 201 plplot.plreplot()
202 plplot.plclear()
202 203 plplot.plend1()
203 plplot.plsstrm(indexPlot)
204 print ''
204 plplot.plsstrm(curr_strm)
205 205
206 206
207 207 def initPlplot(indexPlot,ncol,nrow,winTitle,width,height):
208 208 plplot.plsstrm(indexPlot)
209 209 plplot.plparseopts([winTitle],plplot.PL_PARSE_FULL)
210 210 plplot.plsetopt("geometry", "%dx%d"%(width*ncol,height*nrow))
211 211 plplot.plsdev("xwin")
212 212 plplot.plscolbg(255,255,255)
213 213 plplot.plscol0(1,0,0,0)
214 214 plplot.plinit()
215 215 plplot.plspause(False)
216 216 plplot.plssub(ncol,nrow)
217 217
218 218 def clearData(objGraph):
219 219 objGraph.plotBox(objGraph.xrange[0], objGraph.xrange[1], objGraph.yrange[0], objGraph.yrange[1], "bc", "bc")
220 220
221 221 objGraph.setColor(15) #Setting Line Color to White
222 222
223 223 if objGraph.datatype == "complex":
224 224 objGraph.basicXYPlot(objGraph.xdata,objGraph.ydata.real)
225 225 objGraph.basicXYPlot(objGraph.xdata,objGraph.ydata.imag)
226 226
227 227 if objGraph.datatype == "real":
228 228 objGraph.basicXYPlot(objGraph.xdata,objGraph.ydata)
229 229
230 230 objGraph.setColor(1) #Setting Line Color to Black
231 231 # objGraph.setLineStyle(2)
232 232 # objGraph.plotBox(objGraph.xrange[0], objGraph.xrange[1], objGraph.yrange[0], objGraph.yrange[1], "bcntg", "bc")
233 233 # objGraph.setLineStyle(1)
234 234
235 235 def setStrm(indexPlot):
236 236 plplot.plsstrm(indexPlot)
237 237
238 238 def plFlush():
239 239 plplot.plflush()
240 240
241 241 def setPlTitle(pltitle,color):
242 242 setSubpages(1, 0)
243 243 plplot.pladv(0)
244 244 plplot.plvpor(0., 1., 0., 1.)
245 245
246 246 if color == "black":
247 247 plplot.plcol0(1)
248 248 if color == "white":
249 249 plplot.plcol0(15)
250 250
251 251 plplot.plmtex("t",-1., 0.5, 0.5, pltitle)
252 252
253 253 def setSubpages(ncol,nrow):
254 254 plplot.plssub(ncol,nrow)
255 255
256 256 class BaseGraph:
257 257 __name = None
258 258 __xpos = None
259 259 __ypos = None
260 260 __subplot = None
261 261 __xg = None
262 262 __yg = None
263 263 xdata = None
264 264 ydata = None
265 265 getGrid = True
266 266 xaxisIsTime = False
267 267 deltax = None
268 268 xmin = None
269 269 xmax = None
270 270 def __init__(self,name,subplot,xpos,ypos,xlabel,ylabel,title,szchar,xrange,yrange,zrange=None,deltax=1.0):
271 271 self.setName(name)
272 272 self.setScreenPos(xpos, ypos)
273 273 self.setLabels(xlabel,ylabel,title)
274 274 self.setSubPlot(subplot)
275 275 self.setSizeOfChar(szchar)
276 276 self.setXYZrange(xrange,yrange,zrange)
277 277 self.getGrid = True
278 278 self.xaxisIsTime = False
279 279 self.deltax = deltax
280 280
281 281 def setXYZrange(self,xrange,yrange,zrange):
282 282 self.xrange = xrange
283 283 self.yrange = yrange
284 284 self.zrange = zrange
285 285
286 286 def setName(self, name):
287 287 self.__name = name
288 288
289 289 def setScreenPos(self,xpos,ypos):
290 290 self.__xpos = xpos
291 291 self.__ypos = ypos
292 292
293 293 def setXYData(self,xdata=None,ydata=None,datatype="real"):
294 294 if((xdata != None) and (ydata != None)):
295 295 self.xdata = xdata
296 296 self.ydata = ydata
297 297 self.datatype = datatype
298 298
299 299 if((self.xdata == None) and (self.ydata == None)):
300 300 return None
301 301
302 302 return 1
303 303
304 304
305 305 def setLabels(self,xlabel=None,ylabel=None,title=None):
306 306 if xlabel != None: self.xlabel = xlabel
307 307 if ylabel != None: self.ylabel = ylabel
308 308 if title != None: self.title = title
309 309
310 310 def setSubPlot(self,subplot):
311 311 self.__subplot = subplot
312 312
313 313 def setSizeOfChar(self,szchar):
314 314 self.__szchar = szchar
315 315
316 316 def setLineStyle(self,style):
317 317 plplot.pllsty(style)
318 318
319 319 def setColor(self,color):
320 320 plplot.plcol0(color)
321 321
322 322 def setXAxisAsTime(self,value=False):
323 323 self.xaxisIsTime = value
324 324
325 325 def basicLineTimePlot(self, x, y, xmin=None, xmax=None, ymin=None, ymax=None, colline=1):
326 326
327 327 if xmin == None: xmin = x[0]
328 328 if xmax == None: xmax = x[-1]
329 329 if ymin == None: ymin = y[0]
330 330 if ymax == None: ymax = y[-1]
331 331
332 332 plplot.plcol0(colline)
333 333 plplot.plline(x, y)
334 334 plplot.plcol0(1)
335 335
336 336 def basicXYPlot(self, x, y, xmin=None, xmax=None, ymin=None, ymax=None):
337 337
338 338 if xmin == None: xmin = x[0]
339 339 if xmax == None: xmax = x[-1]
340 340 if ymin == None: ymin = y[0]
341 341 if ymax == None: ymax = y[-1]
342 342
343 343 plplot.plline(x, y)
344 344
345 345 def basicPcolorPlot(self, data, x, y, xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None):
346 346 """
347 347 """
348 348 if xmin == None: xmin = x[0]
349 349 if xmax == None: xmax = x[-1]
350 350 if ymin == None: ymin = y[0]
351 351 if ymax == None: ymax = y[-1]
352 352 if zmin == None: zmin = numpy.nanmin(data)
353 353 if zmax == None: zmax = numpy.nanmax(data)
354 354
355 355 plplot.plimage(data,
356 356 float(x[0]),
357 357 float(x[-1]),
358 358 float(y[0]),
359 359 float(y[-1]),
360 360 float(zmin),
361 361 float(zmax),
362 362 float(xmin),
363 363 float(xmax),
364 364 float(ymin),
365 365 float(ymax)
366 366 )
367 367
368 368 def __getBoxpltr(self, x, y, deltax=None, deltay=None):
369 369
370 370 if not(len(x)>0 and len(y)>0):
371 371 raise ValueError, "x axis and y axis are empty"
372 372
373 373 if deltax == None: deltax = x[-1] - x[-2]
374 374 if deltay == None: deltay = y[-1] - y[-2]
375 375
376 376 x1 = numpy.append(x, x[-1] + deltax)
377 377 y1 = numpy.append(y, y[-1] + deltay)
378 378
379 379 xg = (numpy.multiply.outer(x1, numpy.ones(len(y1))))
380 380 yg = (numpy.multiply.outer(numpy.ones(len(x1)), y1))
381 381
382 382 self.__xg = xg
383 383 self.__yg = yg
384 384
385 385 return xg, yg
386 386
387 387
388 388 def advPcolorPlot(self, data, x, y, xmin=None, xmax=None, ymin=None, ymax=None, zmin=0., zmax=0., deltax=1.0, deltay=None, getGrid = True):
389 389 if getGrid:
390 390 xg, yg = self.__getBoxpltr(x, y, deltax, deltay)
391 391 else:
392 392 xg = self.__xg
393 393 yg = self.__yg
394 394
395 395 plplot.plimagefr(data,
396 396 float(xmin),
397 397 float(xmax),
398 398 float(ymin),
399 399 float(ymax),
400 400 0.,
401 401 0.,
402 402 float(zmin),
403 403 float(zmax),
404 404 plplot.pltr2,
405 405 xg,
406 406 yg)
407 407
408 408
409 409 def colorbarPlot(self, xmin=0., xmax=1., ymin=0., ymax=1.):
410 410 data = numpy.arange(256)
411 411 data = numpy.reshape(data, (1,-1))
412 412
413 413 plplot.plimage(data,
414 414 float(xmin),
415 415 float(xmax),
416 416 float(ymin),
417 417 float(ymax),
418 418 0.,
419 419 255.,
420 420 float(xmin),
421 421 float(xmax),
422 422 float(ymin),
423 423 float(ymax))
424 424
425 425 def plotBox(self, xmin, xmax, ymin, ymax, xopt, yopt, nolabels=False):
426 426
427 427 plplot.plschr(0.0,self.__szchar-0.05)
428 428 plplot.pladv(self.__subplot)
429 429 plplot.plvpor(self.__xpos[0], self.__xpos[1], self.__ypos[0], self.__ypos[1])
430 430 plplot.plwind(float(xmin), # self.xrange[0]
431 431 float(xmax), # self.xrange[1]
432 432 float(ymin), # self.yrange[0]
433 433 float(ymax) # self.yrange[1]
434 434 )
435 435
436 436
437 437
438 438 if self.xaxisIsTime:
439 439 plplot.pltimefmt("%H:%M")
440 440 timedelta = (xmax - xmin + 1)/8.
441 441 plplot.plbox(xopt, timedelta, 3, yopt, 0.0, 0)
442 442 else:
443 443 plplot.plbox(xopt, 0.0, 0, yopt, 0.0, 0)
444 444
445 445
446 446 if not(nolabels):
447 447 plplot.pllab(self.xlabel, self.ylabel, self.title)
448 448
449 449
450 450 def delLabels(self):
451 451 self.setColor(15) #Setting Line Color to White
452 452 plplot.pllab(self.xlabel, self.ylabel, self.title)
453 453 self.setColor(1) #Setting Line Color to Black
454 454
455 455
456 456
457 457 def plotImage(self,x,y,z,xrange,yrange,zrange):
458 458 xi = x[0]
459 459 xf = x[-1]
460 460 yi = y[0]
461 461 yf = y[-1]
462 462
463 463 plplot.plimage(z,
464 464 float(xi),
465 465 float(xf),
466 466 float(yi),
467 467 float(yf),
468 468 float(zrange[0]),
469 469 float(zrange[1]),
470 470 float(xi),
471 471 float(xf),
472 472 float(yrange[0]),
473 473 yrange[1])
474 474
475 475 class LinearPlot:
476 476 linearObjDic = {}
477 477 __xpos = None
478 478 __ypos = None
479 479 def __init__(self,indexPlot,nsubplot,winTitle):
480 480 self.width = 700
481 481 self.height = 150
482 482 ncol = 1
483 483 nrow = nsubplot
484 484 initPlplot(indexPlot,ncol,nrow,winTitle,self.width,self.height)
485 485
486 486
487 487 def setFigure(self,indexPlot):
488 488 setStrm(indexPlot)
489 489
490 490 def setPosition(self):
491 491
492 492 xi = 0.07; xf = 0.9 #0.8,0.7,0.5
493 493 yi = 0.15; yf = 0.8
494 494
495 495 xpos = [xi,xf]
496 496 ypos = [yi,yf]
497 497
498 498 self.__xpos = xpos
499 499 self.__ypos = ypos
500 500
501 501 return xpos,ypos
502 502
503 503 def refresh(self):
504 504 plFlush()
505 505
506 506 def setup(self,subplot,xmin,xmax,ymin,ymax,title,xlabel,ylabel):
507 507 szchar = 1.10
508 508 name = "linear"
509 509 key = name + "%d"%subplot
510 510 xrange = [xmin,xmax]
511 511 yrange = [ymin,ymax]
512 512
513 513 xpos,ypos = self.setPosition()
514 514 linearObj = BaseGraph(name,subplot,xpos,ypos,xlabel,ylabel,title,szchar,xrange,yrange)
515 515 linearObj.plotBox(linearObj.xrange[0], linearObj.xrange[1], linearObj.yrange[0], linearObj.yrange[1], "bcnst", "bcnstv")
516 516 self.linearObjDic[key] = linearObj
517 517
518 518 def plot(self,subplot,x,y,type="power"):
519 519 name = "linear"
520 520 key = name + "%d"%subplot
521 521
522 522 linearObj = self.linearObjDic[key]
523 523 linearObj.plotBox(linearObj.xrange[0], linearObj.xrange[1], linearObj.yrange[0], linearObj.yrange[1], "bcst", "bcst")
524 524
525 525 if linearObj.setXYData() != None:
526 526 clearData(linearObj)
527 527
528 528 else:
529 529 if type.lower() == 'power':
530 530 linearObj.setXYData(x,abs(y),"real")
531 531 if type.lower() == 'iq':
532 532 linearObj.setXYData(x,y,"complex")
533 533
534 534 if type.lower() == 'power':
535 535 colline = 9
536 536 linearObj.basicLineTimePlot(x, abs(y), xmin, xmax, ymin, ymax, colline)
537 537 linearObj.setXYData(x,abs(y),"real")
538 538
539 539 if type.lower() == 'iq':
540 540 colline = 9
541 541 linearObj.basicLineTimePlot(x=x, y=y.real, colline=colline)
542 542 colline = 13
543 543 linearObj.basicLineTimePlot(x=x, y=y.imag, colline=colline)
544 544
545 545 linearObj.setXYData(x,y,"complex")
546 546
547 547 linearObj.plotBox(linearObj.xrange[0], linearObj.xrange[1], linearObj.yrange[0], linearObj.yrange[1], "bcst", "bcst")
548 548
549 549
550 550 # linearObj.plotBox(linearObj.xrange[0], linearObj.xrange[1], linearObj.yrange[0], linearObj.yrange[1], "bc", "bc")
551 551 # linearObj.basicXYPlot(data,y)
552 552 # linearObj.setXYData(data,y)
553 553
554 554
555 555
556 556 class SpectraPlot:
557 557 pcolorObjDic = {}
558 558 colorbarObjDic = {}
559 559 pwprofileObjDic = {}
560 560 showColorbar = None
561 561 showPowerProfile = None
562 562 XAxisAsTime = None
563 563 width = None
564 564 height = None
565 565 __spcxpos = None
566 566 __spcypos = None
567 567 __cmapxpos = None
568 568 __cmapypos = None
569 569 __profxpos = None
570 570 __profypos = None
571 571 __lastTitle = None
572 572
573 573 def __init__(self,indexPlot,nsubplot,winTitle,colormap,showColorbar,showPowerProfile,XAxisAsTime):
574 574 self.width = 460
575 575 self.height = 300
576 576 self.showColorbar = showColorbar
577 577 self.showPowerProfile = showPowerProfile
578 578 self.XAxisAsTime = XAxisAsTime
579 579
580 580 nrow = 2
581 581 if (nsubplot%2)==0:
582 582 ncol = nsubplot/nrow
583 583 else:
584 584 ncol = int(nsubplot)/nrow + 1
585 585
586 586 initPlplot(indexPlot,ncol,nrow,winTitle,self.width,self.height)
587 587 setColormap(colormap)
588 588 self.ncol = ncol
589 589 self.nrow = nrow
590 590
591 591 def setFigure(self,indexPlot):
592 592 setStrm(indexPlot)
593 593
594 594 def setSpectraPos(self): #modificar valores de acuerdo al colorbar y pwprofile
595 595 if self.showPowerProfile: xi = 0.09; xf = 0.6 #0.075
596 596 else: xi = 0.15; xf = 0.8 #0.8,0.7,0.5
597 597 yi = 0.15; yf = 0.80
598 598
599 599 xpos = [xi,xf]
600 600 ypos = [yi,yf]
601 601
602 602 self.__spcxpos = xpos
603 603 self.__spcypos = ypos
604 604
605 605 return xpos,ypos
606 606
607 607 def setColorbarScreenPos(self):
608 608
609 609 xi = self.__spcxpos[1] + 0.03; xf = xi + 0.03
610 610 yi = self.__spcypos[0]; yf = self.__spcypos[1]
611 611
612 612 xpos = [xi,xf]
613 613 ypos = [yi,yf]
614 614
615 615 self.__cmapxpos = xpos
616 616 self.__cmapypos = ypos
617 617
618 618 return xpos,ypos
619 619
620 620 def setPowerprofileScreenPos(self):
621 621
622 622 xi = self.__cmapxpos[1] + 0.07; xf = xi + 0.25
623 623 yi = self.__spcypos[0]; yf = self.__spcypos[1]
624 624
625 625 xpos = [xi,xf]
626 626 ypos = [yi,yf]
627 627
628 628 self.__profxpos = [xi,xf]
629 629 self.__profypos = [yi,yf]
630 630
631 631 return xpos,ypos
632 632
633 633 def setup(self,subplot,xmin,xmax,ymin,ymax,zmin,zmax,title,xlabel,ylabel):
634 634 # Config Spectra plot
635 635 szchar = 0.7
636 636 name = "spc"
637 637 key = name + "%d"%subplot
638 638 xrange = [xmin,xmax]
639 639 yrange = [ymin,ymax]
640 640 zrange = [zmin,zmax]
641 641
642 642 xpos,ypos = self.setSpectraPos()
643 643 pcolorObj = BaseGraph(name,subplot,xpos,ypos,xlabel,ylabel,title,szchar,xrange,yrange,zrange)
644 644 pcolorObj.plotBox(pcolorObj.xrange[0], pcolorObj.xrange[1], pcolorObj.yrange[0], pcolorObj.yrange[1], "bcnst", "bcnstv")
645 645 self.pcolorObjDic[key] = pcolorObj
646 646
647 647 # Config Colorbar
648 648 if self.showColorbar:
649 649 szchar = 0.65
650 650 name = "colorbar"
651 651 key = name + "%d"%subplot
652 652
653 653 xpos,ypos = self.setColorbarScreenPos()
654 654 xrange = [0.,1.]
655 655 yrange = [zmin,zmax]
656 656 cmapObj = BaseGraph(name,subplot,xpos,ypos,"","","dB",szchar,xrange,yrange)
657 657 cmapObj.plotBox(cmapObj.xrange[0], cmapObj.xrange[1], cmapObj.yrange[0], cmapObj.yrange[1], "bc", "bcm")
658 658 cmapObj.colorbarPlot(cmapObj.xrange[0], cmapObj.xrange[1], cmapObj.yrange[0], cmapObj.yrange[1])
659 659 cmapObj.plotBox(cmapObj.xrange[0], cmapObj.xrange[1], cmapObj.yrange[0], cmapObj.yrange[1], "bc", "bcmtsv")
660 660 self.colorbarObjDic[key] = cmapObj
661 661
662 662 # Config Power profile
663 663 if self.showPowerProfile:
664 664 szchar = 0.55
665 665 name = "pwprofile"
666 666 key = name + "%d"%subplot
667 667
668 668 xpos,ypos = self.setPowerprofileScreenPos()
669 669 xrange = [zmin,zmax]
670 670 yrange = [ymin,ymax]
671 671 powObj = BaseGraph(name,subplot,xpos,ypos,"dB","","Power Profile",szchar,xrange,yrange)
672 672 powObj.setLineStyle(2)
673 673 powObj.plotBox(powObj.xrange[0], powObj.xrange[1], powObj.yrange[0], powObj.yrange[1], "bcntg", "bc")
674 674 powObj.setLineStyle(1)
675 675 powObj.plotBox(powObj.xrange[0], powObj.xrange[1], powObj.yrange[0], powObj.yrange[1], "bc", "bc")
676 676 self.pwprofileObjDic[key] = powObj
677 677
678 678 def printTitle(self,pltitle):
679 679 if self.__lastTitle != None:
680 680 setPlTitle(self.__lastTitle,"white")
681 681
682 682 self.__lastTitle = pltitle
683 683
684 684 setPlTitle(pltitle,"black")
685 685
686 686 setSubpages(self.ncol,self.nrow)
687 687
688 688 def plot(self,subplot,x,y,z,subtitle):
689 689 # Spectra plot
690 690
691 691 name = "spc"
692 692 key = name + "%d"%subplot
693 693
694 694 # newx = [x[0],x[-1]]
695 695 pcolorObj = self.pcolorObjDic[key]
696 696 pcolorObj.plotBox(pcolorObj.xrange[0], pcolorObj.xrange[1], pcolorObj.yrange[0], pcolorObj.yrange[1], "bcst", "bcst")
697 697 pcolorObj.delLabels()
698 698 pcolorObj.setLabels(title=subtitle)
699 699
700 700 deltax = None; deltay = None
701 701
702 702 pcolorObj.advPcolorPlot(z,
703 703 x,
704 704 y,
705 705 xmin=pcolorObj.xrange[0],
706 706 xmax=pcolorObj.xrange[1],
707 707 ymin=pcolorObj.yrange[0],
708 708 ymax=pcolorObj.yrange[1],
709 709 zmin=pcolorObj.zrange[0],
710 710 zmax=pcolorObj.zrange[1],
711 711 deltax=deltax,
712 712 deltay=deltay,
713 713 getGrid=pcolorObj.getGrid)
714 714
715 715 pcolorObj.getGrid = False
716 716
717 717 pcolorObj.plotBox(pcolorObj.xrange[0], pcolorObj.xrange[1], pcolorObj.yrange[0], pcolorObj.yrange[1], "bcst", "bcst")
718 718
719 719 # Power Profile
720 720 if self.showPowerProfile:
721 721 power = numpy.average(z, axis=0)
722 722 name = "pwprofile"
723 723 key = name + "%d"%subplot
724 724 powObj = self.pwprofileObjDic[key]
725 725
726 726 if powObj.setXYData() != None:
727 727 clearData(powObj)
728 728 powObj.setLineStyle(2)
729 729 powObj.plotBox(powObj.xrange[0], powObj.xrange[1], powObj.yrange[0], powObj.yrange[1], "bcntg", "bc")
730 730 powObj.setLineStyle(1)
731 731 else:
732 732 powObj.setXYData(power,y)
733 733
734 734 powObj.plotBox(powObj.xrange[0], powObj.xrange[1], powObj.yrange[0], powObj.yrange[1], "bc", "bc")
735 735 powObj.basicXYPlot(power,y)
736 736 powObj.setXYData(power,y)
737 737
738 738 def savePlot(self,indexPlot,path):
739 739
740 740 now = datetime.datetime.now().timetuple()
741 741 file = "spc_img%02d_%03d_%02d%02d%02d"%(indexPlot,now[7],now[3],now[4],now[5])
742 742 filename = os.path.join(path,file+".png")
743 savePlplot(indexPlot,filename,self.ncol,self.nrow,self.width,self.height)
743 width = self.width*self.ncol
744 hei = self.height*self.nrow
745 savePlplot(filename,width,hei)
744 746
745 747
746 748
747 749 def refresh(self):
748 750 plFlush()
749 751
750 752 class RtiPlot:
751 753
752 754 pcolorObjDic = {}
753 755 colorbarObjDic = {}
754 756 pwprofileObjDic = {}
755 757 showColorbar = None
756 758 showPowerProfile = None
757 759 XAxisAsTime = None
758 760 widht = None
759 761 height = None
760 762 __rtixpos = None
761 763 __rtiypos = None
762 764 __cmapxpos = None
763 765 __cmapypos = None
764 766 __profxpos = None
765 767 __profypos = None
766 768
767 769 def __init__(self,indexPlot,nsubplot,winTitle,colormap,showColorbar,showPowerProfile,XAxisAsTime):
768 770 self.width = 700
769 771 self.height = 150
770 772 self.showColorbar = showColorbar
771 773 self.showPowerProfile = showPowerProfile
772 774 self.XAxisAsTime = XAxisAsTime
773 775
774 776 ncol = 1
775 777 nrow = nsubplot
776 778 initPlplot(indexPlot,ncol,nrow,winTitle,self.width,self.height)
777 779 setColormap(colormap)
778 780
779 781 def setFigure(self,indexPlot):
780 782 setStrm(indexPlot)
781 783
782 784 def setRtiScreenPos(self):
783 785
784 786 if self.showPowerProfile: xi = 0.07; xf = 0.65
785 787 else: xi = 0.07; xf = 0.9
786 788 yi = 0.15; yf = 0.80
787 789
788 790 xpos = [xi,xf]
789 791 ypos = [yi,yf]
790 792
791 793 self.__rtixpos = xpos
792 794 self.__rtiypos = ypos
793 795
794 796 return xpos,ypos
795 797
796 798 def setColorbarScreenPos(self):
797 799
798 800 xi = self.__rtixpos[1] + 0.03; xf = xi + 0.03
799 801
800 802 yi = self.__rtiypos[0]; yf = self.__rtiypos[1]
801 803
802 804 xpos = [xi,xf]
803 805 ypos = [yi,yf]
804 806
805 807 self.__cmapxpos = xpos
806 808 self.__cmapypos = ypos
807 809
808 810 return xpos,ypos
809 811
810 812 def setPowerprofileScreenPos(self):
811 813
812 814 xi = self.__cmapxpos[1] + 0.05; xf = xi + 0.20
813 815
814 816 yi = self.__rtiypos[0]; yf = self.__rtiypos[1]
815 817
816 818 xpos = [xi,xf]
817 819 ypos = [yi,yf]
818 820
819 821 self.__profxpos = [xi,xf]
820 822 self.__profypos = [yi,yf]
821 823
822 824 return xpos,ypos
823 825
824 826 def setup(self,subplot,xmin,xmax,ymin,ymax,zmin,zmax,title,xlabel,ylabel,timedata,timezone="lt",npoints=100):
825 827 # Config Rti plot
826 828 szchar = 1.10
827 829 name = "rti"
828 830 key = name + "%d"%subplot
829 831
830 832 # xmin, xmax --> minHour, max Hour : valores que definen el ejex x=[horaInicio,horaFinal]
831 833 thisDateTime = datetime.datetime.fromtimestamp(timedata)
832 834 startDateTime = datetime.datetime(thisDateTime.year,thisDateTime.month,thisDateTime.day,xmin,0,0)
833 835 endDateTime = datetime.datetime(thisDateTime.year,thisDateTime.month,thisDateTime.day,xmax,59,59)
834 836 deltaTime = 0
835 837 if timezone == "lt":
836 838 deltaTime = time.timezone
837 839 startTimeInSecs = time.mktime(startDateTime.timetuple()) - deltaTime
838 840 endTimeInSecs = time.mktime(endDateTime.timetuple()) - deltaTime
839 841
840 842 xrange = [startTimeInSecs,endTimeInSecs]
841 843 totalTimeInXrange = endTimeInSecs - startTimeInSecs + 1.
842 844 deltax = totalTimeInXrange / npoints
843 845
844 846 yrange = [ymin,ymax]
845 847 zrange = [zmin,zmax]
846 848
847 849 xpos,ypos = self.setRtiScreenPos()
848 850 pcolorObj = BaseGraph(name,subplot,xpos,ypos,xlabel,ylabel,title,szchar,xrange,yrange,zrange,deltax)
849 851 if self.XAxisAsTime:
850 852 pcolorObj.setXAxisAsTime(self.XAxisAsTime)
851 853 xopt = "bcnstd"
852 854 yopt = "bcnstv"
853 855 else:
854 856 xopt = "bcnst"
855 857 yopt = "bcnstv"
856 858
857 859 pcolorObj.plotBox(pcolorObj.xrange[0], pcolorObj.xrange[1], pcolorObj.yrange[0], pcolorObj.yrange[1], xopt, yopt)
858 860 self.pcolorObjDic[key] = pcolorObj
859 861
860 862
861 863 # Config Colorbar
862 864 if self.showColorbar:
863 865 szchar = 0.9
864 866 name = "colorbar"
865 867 key = name + "%d"%subplot
866 868
867 869 xpos,ypos = self.setColorbarScreenPos()
868 870 xrange = [0.,1.]
869 871 yrange = [zmin,zmax]
870 872 cmapObj = BaseGraph(name,subplot,xpos,ypos,"","","dB",szchar,xrange,yrange)
871 873 cmapObj.plotBox(cmapObj.xrange[0], cmapObj.xrange[1], cmapObj.yrange[0], cmapObj.yrange[1], "bc", "bcm")
872 874 cmapObj.colorbarPlot(cmapObj.xrange[0], cmapObj.xrange[1], cmapObj.yrange[0], cmapObj.yrange[1])
873 875 cmapObj.plotBox(cmapObj.xrange[0], cmapObj.xrange[1], cmapObj.yrange[0], cmapObj.yrange[1], "bc", "bcmtsv")
874 876 self.colorbarObjDic[key] = cmapObj
875 877
876 878
877 879 # Config Power profile
878 880 if self.showPowerProfile:
879 881 szchar = 0.8
880 882 name = "pwprofile"
881 883 key = name + "%d"%subplot
882 884
883 885 xpos,ypos = self.setPowerprofileScreenPos()
884 886 xrange = [zmin,zmax]
885 887 yrange = [ymin,ymax]
886 888 powObj = BaseGraph(name,subplot,xpos,ypos,"dB","","Power Profile",szchar,xrange,yrange)
887 889 powObj.setLineStyle(2)
888 890 powObj.plotBox(powObj.xrange[0], powObj.xrange[1], powObj.yrange[0], powObj.yrange[1], "bcntg", "bc")
889 891 powObj.setLineStyle(1)
890 892 powObj.plotBox(powObj.xrange[0], powObj.xrange[1], powObj.yrange[0], powObj.yrange[1], "bc", "bc")
891 893 self.pwprofileObjDic[key] = powObj
892 894
893 895
894 896 def plot(self,subplot,x,y,z):
895 897 # RTI plot
896 898 name = "rti"
897 899 key = name + "%d"%subplot
898 900
899 901 data = numpy.reshape(z, (1,-1))
900 902 data = numpy.abs(data)
901 903 data = 10*numpy.log10(data)
902 904 newx = [x,x+1]
903 905
904 906 pcolorObj = self.pcolorObjDic[key]
905 907
906 908 if pcolorObj.xaxisIsTime:
907 909 xopt = "bcstd"
908 910 yopt = "bcst"
909 911 else:
910 912 xopt = "bcst"
911 913 yopt = "bcst"
912 914
913 915 pcolorObj.plotBox(pcolorObj.xrange[0], pcolorObj.xrange[1], pcolorObj.yrange[0], pcolorObj.yrange[1], xopt, yopt)
914 916
915 917 deltax = pcolorObj.deltax
916 918 deltay = None
917 919
918 920 if pcolorObj.xmin == None and pcolorObj.xmax == None:
919 921 pcolorObj.xmin = x
920 922 pcolorObj.xmax = x
921 923
922 924 if x >= pcolorObj.xmax:
923 925 xmin = x
924 926 xmax = x + deltax
925 927 x = [x]
926 928 pcolorObj.advPcolorPlot(data,
927 929 x,
928 930 y,
929 931 xmin=xmin,
930 932 xmax=xmax,
931 933 ymin=pcolorObj.yrange[0],
932 934 ymax=pcolorObj.yrange[1],
933 935 zmin=pcolorObj.zrange[0],
934 936 zmax=pcolorObj.zrange[1],
935 937 deltax=deltax,
936 938 deltay=deltay,
937 939 getGrid=pcolorObj.getGrid)
938 940
939 941 pcolorObj.xmin = xmin
940 942 pcolorObj.xmax = xmax
941 943
942 944
943 945 # Power Profile
944 946 if self.showPowerProfile:
945 947 data = numpy.reshape(data,(numpy.size(data)))
946 948 name = "pwprofile"
947 949 key = name + "%d"%subplot
948 950 powObj = self.pwprofileObjDic[key]
949 951
950 952 if powObj.setXYData() != None:
951 953 clearData(powObj)
952 954 powObj.setLineStyle(2)
953 955 powObj.plotBox(powObj.xrange[0], powObj.xrange[1], powObj.yrange[0], powObj.yrange[1], "bcntg", "bc")
954 956 powObj.setLineStyle(1)
955 957 else:
956 958 powObj.setXYData(data,y)
957 959
958 960 powObj.plotBox(powObj.xrange[0], powObj.xrange[1], powObj.yrange[0], powObj.yrange[1], "bc", "bc")
959 961 powObj.basicXYPlot(data,y)
960 962 powObj.setXYData(data,y)
961 963
962 964 def refresh(self):
963 965 plFlush() No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,213 +1,213
1 1 import numpy
2 2 from Model.Spectra import Spectra
3 3
4 4 def hildebrand_sekhon(Data, navg):
5 5 """
6 6 This method is for the objective determination of de noise level in Doppler spectra. This
7 7 implementation technique is based on the fact that the standard deviation of the spectral
8 8 densities is equal to the mean spectral density for white Gaussian noise
9 9
10 10 Inputs:
11 11 Data : heights
12 12 navg : numbers of averages
13 13
14 14 Return:
15 15 -1 : any error
16 16 anoise : noise's level
17 17 """
18 18
19 19 divisor = 8
20 20 ratio = 7 / divisor
21 21 data = Data.reshape(-1)
22 22 npts = data.size #numbers of points of the data
23 23
24 24 if npts < 32:
25 25 print "error in noise - requires at least 32 points"
26 26 return -1.0
27 27
28 28 # data sorted in ascending order
29 29 nmin = int(npts/divisor + ratio);
30 30 s = 0.0
31 31 s2 = 0.0
32 32 data2 = data[:npts]
33 33 data2.sort()
34 34
35 35 for i in range(nmin):
36 36 s += data2[i]
37 37 s2 += data2[i]**2;
38 38
39 39 icount = nmin
40 40 iflag = 0
41 41
42 42 for i in range(nmin, npts):
43 43 s += data2[i];
44 44 s2 += data2[i]**2
45 45 icount=icount+1;
46 46 p = s / float(icount);
47 47 p2 = p**2;
48 48 q = s2 / float(icount) - p2;
49 49 leftc = p2;
50 50 rightc = q * float(navg);
51 51
52 52 if leftc > rightc:
53 53 iflag = 1; #No weather signal
54 54 # Signal detect: R2 < 1 (R2 = leftc/rightc)
55 55 if(leftc < rightc):
56 56 if iflag:
57 57 break
58 58
59 59 anoise = 0.0;
60 60 for j in range(i):
61 61 anoise += data2[j];
62 62
63 63 anoise = anoise / float(i);
64 64
65 65 return anoise;
66 66
67 67 def sorting_bruce(Data, navg):
68 68 sortdata = numpy.sort(Data)
69 69 lenOfData = len(Data)
70 70 nums_min = lenOfData/10
71 71
72 72 if (lenOfData/10) > 0:
73 73 nums_min = lenOfData/10
74 74 else:
75 75 nums_min = 0
76 76
77 77 rtest = 1.0 + 1.0/navg
78 78
79 79 sum = 0.
80 80
81 81 sumq = 0.
82 82
83 83 j = 0
84 84
85 85 cont = 1
86 86
87 87 while((cont==1)and(j<lenOfData)):
88 88
89 89 sum += sortdata[j]
90 90
91 91 sumq += sortdata[j]**2
92 92
93 93 j += 1
94 94
95 95 if j > nums_min:
96 96 if ((sumq*j) <= (rtest*sum**2)):
97 97 lnoise = sum / j
98 98 else:
99 99 j = j - 1
100 100 sum = sum - sordata[j]
101 101 sumq = sumq - sordata[j]**2
102 102 cont = 0
103 103
104 104 if j == nums_min:
105 105 lnoise = sum /j
106 106
107 107 return lnoise
108 108
109 109 class Noise:
110 110 """
111 111 Clase que implementa los metodos necesarios para deternimar el nivel de ruido en un Spectro Doppler
112 112 """
113 113 data = None
114 114 noise = None
115 115 dim = None
116 116
117 117 def __init__(self, data=None):
118 118 """
119 119 Inicializador de la clase Noise para la la determinacion del nivel de ruido en un Spectro Doppler.
120 120
121 121 Inputs:
122 122 data: Numpy array de la forma nChan x nHeis x nProfiles
123 123
124 124 Affected:
125 125 self.noise
126 126
127 127 Return:
128 128 None
129 129 """
130 130
131 131 self.data = data
132 132 self.dim = None
133 133 self.nChannels = None
134 134 self.noise = None
135 135
136 136 def setNoise(self, data):
137 137 """
138 138 Inicializador de la clase Noise para la la determinacion del nivel de ruido en un Spectro Doppler.
139 139
140 140 Inputs:
141 141 data: Numpy array de la forma nChan x nHeis x nProfiles
142 142
143 143 Affected:
144 144 self.noise
145 145
146 146 Return:
147 147 None
148 148 """
149 149
150 150 if data == None:
151 return 0
151 raise ValueError, "The data value is not defined"
152 152
153 153 shape = data.shape
154 154 self.dim = len(shape)
155 155 if self.dim == 3:
156 156 nChan, nProfiles, nHeis = shape
157 157 elif self.dim == 2:
158 158 nChan, nHeis = shape
159 159 else:
160 160 raise ValueError, ""
161 161
162 162 self.nChannels = nChan
163 163 self.data = data.copy()
164 164 self.noise = numpy.zeros(nChan)
165 165
166 166 return 1
167 167
168 168
169 169 def byHildebrand(self, navg=1):
170 170 """
171 171 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
172 172
173 173 Return:
174 174 noiselevel
175 175 """
176 176
177 177 daux = None
178 178
179 179 for channel in range(self.nChannels):
180 180 daux = self.data[channel,:,:]
181 181 self.noise[channel] = hildebrand_sekhon(daux, navg)
182 182 return self.noise
183 183
184 184 def byWindow(self, heiIndexMin, heiIndexMax, freqIndexMin, freqIndexMax):
185 185 """
186 186 Determina el ruido del canal utilizando la ventana indicada con las coordenadas:
187 187 (heiIndexMIn, freqIndexMin) hasta (heiIndexMax, freqIndexMAx)
188 188
189 189 Inputs:
190 190 heiIndexMin: Limite inferior del eje de alturas
191 191 heiIndexMax: Limite superior del eje de alturas
192 192 freqIndexMin: Limite inferior del eje de frecuencia
193 193 freqIndexMax: Limite supoerior del eje de frecuencia
194 194 """
195 195
196 196 data = self.data[:, heiIndexMin:heiIndexMax, freqIndexMin:freqIndexMax]
197 197
198 198 for channel in range(self.nChannels):
199 199 daux = data[channel,:,:]
200 200 self.noise[channel] = numpy.average(daux)
201 201
202 202 return self.noise
203 203
204 204 def bySort(self,navg = 1):
205 205 daux = None
206 206
207 207 for channel in range(self.nChannels):
208 208 daux = self.data[channel,:,:]
209 209 self.noise[channel] = sorting_bruce(daux, navg)
210 210
211 211 return self.noise
212 212
213 213 No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,689 +1,693
1 1 '''
2 2 Created on Feb 7, 2012
3 3
4 4 @author $Author$
5 5 @version $Id$
6 6 '''
7 7 import os, sys
8 8 import numpy
9 9
10 10 path = os.path.split(os.getcwd())[0]
11 11 sys.path.append(path)
12 12
13 13 from Model.Spectra import Spectra
14 14 from IO.SpectraIO import SpectraWriter
15 15 from Graphics.SpectraPlot import Spectrum
16 16 from JRONoise import Noise
17 17
18 18 class SpectraProcessor:
19 19 '''
20 20 classdocs
21 21 '''
22 22
23 23 dataInObj = None
24 24
25 25 dataOutObj = None
26 26
27 27 noiseObj = None
28 28
29 29 integratorObjList = []
30 30
31 31 decoderObjList = []
32 32
33 33 writerObjList = []
34 34
35 35 plotterObjList = []
36 36
37 37 integratorObjIndex = None
38 38
39 39 decoderObjIndex = None
40 40
41 41 writerObjIndex = None
42 42
43 43 plotterObjIndex = None
44 44
45 45 buffer = None
46 46
47 47 profIndex = 0
48 48
49 49 nFFTPoints = None
50 50
51 51 nChannels = None
52 52
53 53 nHeights = None
54 54
55 55 nPairs = None
56 56
57 57 pairList = None
58 58
59 59
60 60 def __init__(self):
61 61 '''
62 62 Constructor
63 63 '''
64 64
65 65 self.integratorObjIndex = None
66 66 self.decoderObjIndex = None
67 67 self.writerObjIndex = None
68 68 self.plotterObjIndex = None
69 69
70 70 self.integratorObjList = []
71 71 self.decoderObjList = []
72 72 self.writerObjList = []
73 73 self.plotterObjList = []
74 74
75 75 self.noiseObj = Noise()
76 76 self.buffer = None
77 77 self.profIndex = 0
78 78
79 79 def setup(self, dataInObj=None, dataOutObj=None, nFFTPoints=None, pairList=None):
80 80
81 81 if dataInObj == None:
82 82 raise ValueError, ""
83 83
84 84 if nFFTPoints == None:
85 85 raise ValueError, ""
86 86
87 87 self.dataInObj = dataInObj
88 88
89 89 if dataOutObj == None:
90 90 dataOutObj = Spectra()
91 91
92 92 self.dataOutObj = dataOutObj
93 93 self.noiseObj = Noise()
94 94
95 95 ##########################################
96 96 self.nFFTPoints = nFFTPoints
97 97 self.nChannels = self.dataInObj.nChannels
98 98 self.nHeights = self.dataInObj.nHeights
99 99 self.pairList = pairList
100 100 if pairList != None:
101 101 self.nPairs = len(pairList)
102 102 else:
103 103 self.nPairs = 0
104 104
105 105 self.dataOutObj.heightList = self.dataInObj.heightList
106 106 self.dataOutObj.channelIndexList = self.dataInObj.channelIndexList
107 107 self.dataOutObj.m_BasicHeader = self.dataInObj.m_BasicHeader.copy()
108 108 self.dataOutObj.m_ProcessingHeader = self.dataInObj.m_ProcessingHeader.copy()
109 109 self.dataOutObj.m_RadarControllerHeader = self.dataInObj.m_RadarControllerHeader.copy()
110 110 self.dataOutObj.m_SystemHeader = self.dataInObj.m_SystemHeader.copy()
111 111
112 112 self.dataOutObj.dataType = self.dataInObj.dataType
113 113 self.dataOutObj.nPairs = self.nPairs
114 114 self.dataOutObj.nChannels = self.nChannels
115 115 self.dataOutObj.nProfiles = self.nFFTPoints
116 116 self.dataOutObj.nHeights = self.nHeights
117 117 self.dataOutObj.nFFTPoints = self.nFFTPoints
118 118 #self.dataOutObj.data = None
119 119
120 120 self.dataOutObj.m_SystemHeader.numChannels = self.nChannels
121 121 self.dataOutObj.m_SystemHeader.nProfiles = self.nFFTPoints
122 122
123 123 self.dataOutObj.m_ProcessingHeader.totalSpectra = self.nChannels + self.nPairs
124 124 self.dataOutObj.m_ProcessingHeader.profilesPerBlock = self.nFFTPoints
125 125 self.dataOutObj.m_ProcessingHeader.numHeights = self.nHeights
126 126 self.dataOutObj.m_ProcessingHeader.shif_fft = True
127 127
128 128 spectraComb = numpy.zeros( (self.nChannels+self.nPairs)*2,numpy.dtype('u1'))
129 129 k = 0
130 130 for i in range( 0,self.nChannels*2,2 ):
131 131 spectraComb[i] = k
132 132 spectraComb[i+1] = k
133 133 k += 1
134 134
135 135 k *= 2
136 136
137 137 if self.pairList != None:
138 138
139 139 for pair in self.pairList:
140 140 spectraComb[k] = pair[0]
141 141 spectraComb[k+1] = pair[1]
142 142 k += 2
143 143
144 144 self.dataOutObj.m_ProcessingHeader.spectraComb = spectraComb
145 145
146 146 return self.dataOutObj
147 147
148 148 def init(self):
149 149
150 150 self.integratorObjIndex = 0
151 151 self.decoderObjIndex = 0
152 152 self.writerObjIndex = 0
153 153 self.plotterObjIndex = 0
154 154
155 155 if self.dataInObj.type == "Voltage":
156 156
157 157 if self.buffer == None:
158 158 self.buffer = numpy.zeros((self.nChannels,
159 159 self.nFFTPoints,
160 160 self.nHeights),
161 161 dtype='complex')
162 162
163 163 self.buffer[:,self.profIndex,:] = self.dataInObj.data
164 164 self.profIndex += 1
165 165
166 166 if self.profIndex == self.nFFTPoints:
167 167 self.__getFft()
168 168 self.dataOutObj.flagNoData = False
169 169
170 170 self.buffer = None
171 171 self.profIndex = 0
172 172 return
173 173
174 174 self.dataOutObj.flagNoData = True
175 175
176 176 return
177 177
178 178 #Other kind of data
179 179 if self.dataInObj.type == "Spectra":
180 180 self.dataOutObj.copy(self.dataInObj)
181 181 self.dataOutObj.flagNoData = False
182 182 return
183 183
184 184 raise ValueError, "The datatype is not valid"
185 185
186 186 def __getFft(self):
187 187 """
188 188 Convierte valores de Voltaje a Spectra
189 189
190 190 Affected:
191 191 self.dataOutObj.data_spc
192 192 self.dataOutObj.data_cspc
193 193 self.dataOutObj.data_dc
194 194 self.dataOutObj.heightList
195 195 self.dataOutObj.m_BasicHeader
196 196 self.dataOutObj.m_ProcessingHeader
197 197 self.dataOutObj.m_RadarControllerHeader
198 198 self.dataOutObj.m_SystemHeader
199 199 self.profIndex
200 200 self.buffer
201 201 self.dataOutObj.flagNoData
202 202 self.dataOutObj.dataType
203 203 self.dataOutObj.nPairs
204 204 self.dataOutObj.nChannels
205 205 self.dataOutObj.nProfiles
206 206 self.dataOutObj.m_SystemHeader.numChannels
207 207 self.dataOutObj.m_ProcessingHeader.totalSpectra
208 208 self.dataOutObj.m_ProcessingHeader.profilesPerBlock
209 209 self.dataOutObj.m_ProcessingHeader.numHeights
210 210 self.dataOutObj.m_ProcessingHeader.spectraComb
211 211 self.dataOutObj.m_ProcessingHeader.shif_fft
212 212 """
213 213 if self.dataInObj.flagNoData:
214 214 return 0
215 215
216 216 fft_volt = numpy.fft.fft(self.buffer,axis=1)
217 217 dc = fft_volt[:,0,:]
218 218
219 219 #calculo de self-spectra
220 220 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
221 221 spc = numpy.abs(fft_volt * numpy.conjugate(fft_volt))
222 222
223 223 blocksize = 0
224 224 blocksize += dc.size
225 225 blocksize += spc.size
226 226
227 227 cspc = None
228 228 pairIndex = 0
229 229 if self.pairList != None:
230 230 #calculo de cross-spectra
231 231 cspc = numpy.zeros((self.nPairs, self.nFFTPoints, self.nHeights), dtype='complex')
232 232 for pair in self.pairList:
233 233 cspc[pairIndex,:,:] = numpy.abs(fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:]))
234 234 pairIndex += 1
235 235 blocksize += cspc.size
236 236
237 237 self.dataOutObj.data_spc = spc
238 238 self.dataOutObj.data_cspc = cspc
239 239 self.dataOutObj.data_dc = dc
240 240 self.dataOutObj.m_ProcessingHeader.blockSize = blocksize
241 241 self.dataOutObj.m_BasicHeader.utc = self.dataInObj.m_BasicHeader.utc
242
242
243 self.getNoise()
243 244
244 245 def addWriter(self,wrpath):
245 246 objWriter = SpectraWriter(self.dataOutObj)
246 247 objWriter.setup(wrpath)
247 248 self.writerObjList.append(objWriter)
248 249
249 250 def addPlotter(self,index=None):
250 251 if index==None:
251 252 index = self.plotterObjIndex
252 253
253 254 plotObj = Spectrum(self.dataOutObj, index)
254 255 self.plotterObjList.append(plotObj)
255 256
256 257 def addIntegrator(self,N,timeInterval):
257 258
258 259 objIncohInt = IncoherentIntegration(N,timeInterval)
259 260 self.integratorObjList.append(objIncohInt)
260 261
261 262 def writeData(self, wrpath):
262 263 if self.dataOutObj.flagNoData:
263 264 return 0
264 265
265 266 if len(self.writerObjList) <= self.writerObjIndex:
266 267 self.addWriter(wrpath)
267 268
268 269 self.writerObjList[self.writerObjIndex].putData()
269 270
270 271 self.writerObjIndex += 1
271 272
272 273 def plotData(self,
273 274 xmin=None,
274 275 xmax=None,
275 276 ymin=None,
276 277 ymax=None,
277 278 zmin=None,
278 279 zmax=None,
279 280 titleList=None,
280 281 xlabelList=None,
281 282 ylabelList=None,
282 283 winTitle='',
283 284 colormap="br_green",
284 285 showColorbar=False,
285 286 showPowerProfile=False,
286 287 XAxisAsTime=False,
287 288 save=False,
288 289 index=None):
289 290
290 291 if self.dataOutObj.flagNoData:
291 292 return 0
292 293
293 294 if len(self.plotterObjList) <= self.plotterObjIndex:
294 295 self.addPlotter(index)
295 296
296 297 self.plotterObjList[self.plotterObjIndex].plotData(xmin,
297 298 xmax,
298 299 ymin,
299 300 ymax,
300 301 zmin,
301 302 zmax,
302 303 titleList,
303 304 xlabelList,
304 305 ylabelList,
305 306 winTitle,
306 307 colormap,
307 308 showColorbar,
308 309 showPowerProfile,
309 310 XAxisAsTime,
310 311 save)
311 312
312 313 self.plotterObjIndex += 1
313 314
314 315 def integrator(self, N=None, timeInterval=None):
315 316
316 317 if self.dataOutObj.flagNoData:
317 318 return 0
318 319
319 320 if len(self.integratorObjList) <= self.integratorObjIndex:
320 321 self.addIntegrator(N,timeInterval)
321 322
322 323 myIncohIntObj = self.integratorObjList[self.integratorObjIndex]
323 324 myIncohIntObj.exe(data=self.dataOutObj.data_spc,timeOfData=self.dataOutObj.m_BasicHeader.utc)
324 325
325 326 if myIncohIntObj.isReady:
326 327 self.dataOutObj.data_spc = myIncohIntObj.data
327 328 self.dataOutObj.nAvg = myIncohIntObj.navg
328 329 self.dataOutObj.m_ProcessingHeader.incoherentInt *= myIncohIntObj.navg
329 330 #print "myIncohIntObj.navg: ",myIncohIntObj.navg
330 331 self.dataOutObj.flagNoData = False
331 332
332 self.getNoise(type="hildebrand",parm=myIncohIntObj.navg)
333 self.getNoise(type="hildebrand")
333 334 # self.getNoise(type="sort", parm=16)
334 335
335 336 else:
336 337 self.dataOutObj.flagNoData = True
337 338
338 339 self.integratorObjIndex += 1
339 340
340 341 """Calcular el ruido"""
341 342 # self.getNoise(type="hildebrand", parm=1)
342 343
343 344 def removeDC(self, type):
344 345
345 346 if self.dataOutObj.flagNoData:
346 347 return 0
347 348
348 349 def removeInterference(self):
349 350
350 351 if self.dataOutObj.flagNoData:
351 352 return 0
352 353
353 354 def removeSatellites(self):
354 355
355 356 if self.dataOutObj.flagNoData:
356 357 return 0
357 358
358 359 def getNoise(self, type="hildebrand", parm=None):
359 360
361 if parm == None:
362 parm =self.dataOutObj.m_ProcessingHeader.incoherentInt
363
360 364 self.noiseObj.setNoise(self.dataOutObj.data_spc)
361 365
362 366 if type == "hildebrand":
363 367 noise = self.noiseObj.byHildebrand(parm)
364 368
365 369 if type == "window":
366 370 noise = self.noiseObj.byWindow(parm)
367 371
368 372 if type == "sort":
369 373 noise = self.noiseObj.bySort(parm)
370 374
371 375 self.dataOutObj.noise = noise
372 376 # print 10*numpy.log10(noise)
373 377
374 378 def selectChannels(self, channelList, pairList=[]):
375 379
376 380 channelIndexList = []
377 381
378 382 for channel in channelList:
379 383 if channel in self.dataOutObj.channelList:
380 384 index = self.dataOutObj.channelList.index(channel)
381 385 channelIndexList.append(index)
382 386
383 387 pairIndexList = []
384 388
385 389 for pair in pairList:
386 390 if pair in self.dataOutObj.pairList:
387 391 index = self.dataOutObj.pairList.index(pair)
388 392 pairIndexList.append(index)
389 393
390 394 self.selectChannelsByIndex(channelIndexList, pairIndexList)
391 395
392 396 def selectChannelsByIndex(self, channelIndexList, pairIndexList=[]):
393 397 """
394 398 Selecciona un bloque de datos en base a canales y pares segun el
395 399 channelIndexList y el pairIndexList
396 400
397 401 Input:
398 402 channelIndexList : lista de indices de los canales a seleccionar por ej.
399 403
400 404 Si tenemos los canales
401 405
402 406 self.channelList = (2,3,5,7)
403 407
404 408 y deseamos escoger los canales (3,7)
405 409 entonces colocaremos el parametro
406 410
407 411 channelndexList = (1,3)
408 412
409 413 pairIndexList : tupla de indice depares que se desea selecionar por ej.
410 414
411 415 Si tenemos los pares :
412 416
413 417 ( (0,1), (0,2), (1,3), (2,5) )
414 418
415 419 y deseamos seleccionar los pares ((0,2), (2,5))
416 420 entonces colocaremos el parametro
417 421
418 422 pairIndexList = (1,3)
419 423
420 424 Affected:
421 425 self.dataOutObj.data_spc
422 426 self.dataOutObj.data_cspc
423 427 self.dataOutObj.data_dc
424 428 self.dataOutObj.nChannels
425 429 self.dataOutObj.nPairs
426 430 self.dataOutObj.m_ProcessingHeader.spectraComb
427 431 self.dataOutObj.m_SystemHeader.numChannels
428 432
429 433 self.dataOutObj.noise
430 434 Return:
431 435 None
432 436 """
433 437
434 438 if self.dataOutObj.flagNoData:
435 439 return 0
436 440
437 441 if pairIndexList == []:
438 442 pairIndexList = numpy.arange(len(self.dataOutObj.pairList))
439 443
440 444 nChannels = len(channelIndexList)
441 445 nPairs = len(pairIndexList)
442 446
443 447 blocksize = 0
444 448 #self spectra
445 449 spc = self.dataOutObj.data_spc[channelIndexList,:,:]
446 450 blocksize += spc.size
447 451
448 452 cspc = None
449 453 if pairIndexList != []:
450 454 cspc = self.dataOutObj.data_cspc[pairIndexList,:,:]
451 455 blocksize += cspc.size
452 456
453 457 #DC channel
454 458 dc = None
455 459 if self.dataOutObj.m_ProcessingHeader.flag_dc:
456 460 dc = self.dataOutObj.data_dc[channelIndexList,:]
457 461 blocksize += dc.size
458 462
459 463 #Almacenar las combinaciones de canales y cros espectros
460 464
461 465 spectraComb = numpy.zeros( (nChannels+nPairs)*2,numpy.dtype('u1'))
462 466 i = 0
463 467 for spcChannel in channelIndexList:
464 468 spectraComb[i] = spcChannel
465 469 spectraComb[i+1] = spcChannel
466 470 i += 2
467 471
468 472 if pairList != None:
469 473 for pair in pairList:
470 474 spectraComb[i] = pair[0]
471 475 spectraComb[i+1] = pair[1]
472 476 i += 2
473 477
474 478 #######
475 479
476 480 self.dataOutObj.data_spc = spc
477 481 self.dataOutObj.data_cspc = cspc
478 482 self.dataOutObj.data_dc = dc
479 483 self.dataOutObj.nChannels = nChannels
480 484 self.dataOutObj.nPairs = nPairs
481 485
482 486 self.dataOutObj.channelIndexList = channelIndexList
483 487
484 488 self.dataOutObj.m_ProcessingHeader.spectraComb = spectraComb
485 489 self.dataOutObj.m_ProcessingHeader.totalSpectra = nChannels + nPairs
486 490 self.dataOutObj.m_SystemHeader.numChannels = nChannels
487 491 self.dataOutObj.nChannels = nChannels
488 492 self.dataOutObj.m_ProcessingHeader.blockSize = blocksize
489 493
490 494 if cspc == None:
491 495 self.dataOutObj.m_ProcessingHeader.flag_dc = False
492 496 if dc == None:
493 497 self.dataOutObj.m_ProcessingHeader.flag_cpsc = False
494 498
495 499 def selectHeightsByValue(self, minHei, maxHei):
496 500 """
497 501 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
498 502 minHei <= height <= maxHei
499 503
500 504 Input:
501 505 minHei : valor minimo de altura a considerar
502 506 maxHei : valor maximo de altura a considerar
503 507
504 508 Affected:
505 509 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
506 510
507 511 Return:
508 512 None
509 513 """
510 514
511 515 if self.dataOutObj.flagNoData:
512 516 return 0
513 517
514 518 if (minHei < self.dataOutObj.heightList[0]) or (minHei > maxHei):
515 519 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
516 520
517 521 if (maxHei > self.dataOutObj.heightList[-1]):
518 522 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
519 523
520 524 minIndex = 0
521 525 maxIndex = 0
522 526 data = self.dataOutObj.heightList
523 527
524 528 for i,val in enumerate(data):
525 529 if val < minHei:
526 530 continue
527 531 else:
528 532 minIndex = i;
529 533 break
530 534
531 535 for i,val in enumerate(data):
532 536 if val <= maxHei:
533 537 maxIndex = i;
534 538 else:
535 539 break
536 540
537 541 self.selectHeightsByIndex(minIndex, maxIndex)
538 542
539 543 def selectHeightsByIndex(self, minIndex, maxIndex):
540 544 """
541 545 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
542 546 minIndex <= index <= maxIndex
543 547
544 548 Input:
545 549 minIndex : valor minimo de altura a considerar
546 550 maxIndex : valor maximo de altura a considerar
547 551
548 552 Affected:
549 553 self.dataOutObj.data_spc
550 554 self.dataOutObj.data_cspc
551 555 self.dataOutObj.data_dc
552 556 self.dataOutObj.heightList
553 557 self.dataOutObj.nHeights
554 558 self.dataOutObj.m_ProcessingHeader.numHeights
555 559 self.dataOutObj.m_ProcessingHeader.blockSize
556 560 self.dataOutObj.m_ProcessingHeader.firstHeight
557 561 self.dataOutObj.m_RadarControllerHeader.numHeights
558 562
559 563 Return:
560 564 None
561 565 """
562 566
563 567 if self.dataOutObj.flagNoData:
564 568 return 0
565 569
566 570 if (minIndex < 0) or (minIndex > maxIndex):
567 571 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
568 572
569 573 if (maxIndex >= self.dataOutObj.nHeights):
570 574 raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
571 575
572 576 nChannels = self.dataOutObj.nChannels
573 577 nPairs = self.dataOutObj.nPairs
574 578 nProfiles = self.dataOutObj.nProfiles
575 579 dataType = self.dataOutObj.dataType
576 580 nHeights = maxIndex - minIndex + 1
577 581 blockSize = 0
578 582
579 583 #self spectra
580 584 spc = self.dataOutObj.data_spc[:,:,minIndex:maxIndex+1]
581 585 blockSize += spc.size
582 586
583 587 #cross spectra
584 588 cspc = None
585 589 if self.dataOutObj.data_cspc != None:
586 590 cspc = self.dataOutObj.data_cspc[:,:,minIndex:maxIndex+1]
587 591 blockSize += cspc.size
588 592
589 593 #DC channel
590 594 dc = self.dataOutObj.data_dc[:,minIndex:maxIndex+1]
591 595 blockSize += dc.size
592 596
593 597 self.dataOutObj.data_spc = spc
594 598 if cspc != None:
595 599 self.dataOutObj.data_cspc = cspc
596 600 self.dataOutObj.data_dc = dc
597 601
598 602 firstHeight = self.dataOutObj.heightList[minIndex]
599 603
600 604 self.dataOutObj.nHeights = nHeights
601 605 self.dataOutObj.m_ProcessingHeader.blockSize = blockSize
602 606 self.dataOutObj.m_ProcessingHeader.numHeights = nHeights
603 607 self.dataOutObj.m_ProcessingHeader.firstHeight = firstHeight
604 608 self.dataOutObj.m_RadarControllerHeader.numHeights = nHeights
605 609
606 610 self.dataOutObj.heightList = self.dataOutObj.heightList[minIndex:maxIndex+1]
607 611
608 612
609 613 class IncoherentIntegration:
610 614
611 615 integ_counter = None
612 616 data = None
613 617 navg = None
614 618 buffer = None
615 619 nIncohInt = None
616 620
617 621 def __init__(self, N = None, timeInterval = None):
618 622 """
619 623 N
620 624 timeInterval - interval time [min], integer value
621 625 """
622 626
623 627 self.data = None
624 628 self.navg = None
625 629 self.buffer = None
626 630 self.timeOut = None
627 631 self.exitCondition = False
628 632 self.isReady = False
629 633 self.nIncohInt = N
630 634 self.integ_counter = 0
631 635 if timeInterval!=None:
632 636 self.timeIntervalInSeconds = timeInterval * 60. #if (type(timeInterval)!=integer) -> change this line
633 637
634 638 if ((timeInterval==None) and (N==None)):
635 639 print 'N = None ; timeInterval = None'
636 640 sys.exit(0)
637 641 elif timeInterval == None:
638 642 self.timeFlag = False
639 643 else:
640 644 self.timeFlag = True
641 645
642 646
643 647 def exe(self,data,timeOfData):
644 648 """
645 649 data
646 650
647 651 timeOfData [seconds]
648 652 """
649 653
650 654 if self.timeFlag:
651 655 if self.timeOut == None:
652 656 self.timeOut = timeOfData + self.timeIntervalInSeconds
653 657
654 658 if timeOfData < self.timeOut:
655 659 if self.buffer == None:
656 660 self.buffer = data
657 661 else:
658 662 self.buffer = self.buffer + data
659 663 self.integ_counter += 1
660 664 else:
661 665 self.exitCondition = True
662 666
663 667 else:
664 668 if self.integ_counter < self.nIncohInt:
665 669 if self.buffer == None:
666 670 self.buffer = data
667 671 else:
668 672 self.buffer = self.buffer + data
669 673
670 674 self.integ_counter += 1
671 675
672 676 if self.integ_counter == self.nIncohInt:
673 677 self.exitCondition = True
674 678
675 679 if self.exitCondition:
676 680 self.data = self.buffer
677 681 self.navg = self.integ_counter
678 682 self.isReady = True
679 683 self.buffer = None
680 684 self.timeOut = None
681 685 self.integ_counter = 0
682 686 self.exitCondition = False
683 687
684 688 if self.timeFlag:
685 689 self.buffer = data
686 690 self.timeOut = timeOfData + self.timeIntervalInSeconds
687 691 else:
688 692 self.isReady = False
689 693 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now