##// END OF EJS Templates
schain
RSS
Miguel Valdez -
r225:89c13353a3f3
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,511 +1,512
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 os, sys
8 8 import copy
9 9 import numpy
10 10 import datetime
11 11
12 12 from jroheaderIO import SystemHeader, RadarControllerHeader
13 13
14 14 def hildebrand_sekhon(data, navg):
15 15 """
16 16 This method is for the objective determination of de noise level in Doppler spectra. This
17 17 implementation technique is based on the fact that the standard deviation of the spectral
18 18 densities is equal to the mean spectral density for white Gaussian noise
19 19
20 20 Inputs:
21 21 Data : heights
22 22 navg : numbers of averages
23 23
24 24 Return:
25 25 -1 : any error
26 26 anoise : noise's level
27 27 """
28 28
29 29 dataflat = data.copy().reshape(-1)
30 30 dataflat.sort()
31 31 npts = dataflat.size #numbers of points of the data
32 npts_noise = 0.2*npts
32 33
33 34 if npts < 32:
34 35 print "error in noise - requires at least 32 points"
35 36 return -1.0
36 37
37 38 dataflat2 = numpy.power(dataflat,2)
38 39
39 40 cs = numpy.cumsum(dataflat)
40 41 cs2 = numpy.cumsum(dataflat2)
41 42
42 43 # data sorted in ascending order
43 44 nmin = int((npts + 7.)/8)
44 45
45 46 for i in range(nmin, npts):
46 47 s = cs[i]
47 48 s2 = cs2[i]
48 49 p = s / float(i);
49 50 p2 = p**2;
50 51 q = s2 / float(i) - p2;
51 52 leftc = p2;
52 53 rightc = q * float(navg);
53 54 R2 = leftc/rightc
54 55
55 56 # Signal detect: R2 < 1 (R2 = leftc/rightc)
56 57 if R2 < 1:
57 58 npts_noise = i
58 59 break
59 60
60 61
61 62 anoise = numpy.average(dataflat[0:npts_noise])
62 63
63 64 return anoise;
64 65
65 66 def sorting_bruce(data, navg):
66 67
67 68 data = data.copy()
68 69
69 70 sortdata = numpy.sort(data)
70 71 lenOfData = len(data)
71 72 nums_min = lenOfData/10
72 73
73 74 if (lenOfData/10) > 0:
74 75 nums_min = lenOfData/10
75 76 else:
76 77 nums_min = 0
77 78
78 79 rtest = 1.0 + 1.0/navg
79 80
80 81 sum = 0.
81 82
82 83 sumq = 0.
83 84
84 85 j = 0
85 86
86 87 cont = 1
87 88
88 89 while((cont==1)and(j<lenOfData)):
89 90
90 91 sum += sortdata[j]
91 92
92 93 sumq += sortdata[j]**2
93 94
94 95 j += 1
95 96
96 97 if j > nums_min:
97 98 if ((sumq*j) <= (rtest*sum**2)):
98 99 lnoise = sum / j
99 100 else:
100 101 j = j - 1
101 102 sum = sum - sordata[j]
102 103 sumq = sumq - sordata[j]**2
103 104 cont = 0
104 105
105 106 if j == nums_min:
106 107 lnoise = sum /j
107 108
108 109 return lnoise
109 110
110 111 class JROData:
111 112
112 113 # m_BasicHeader = BasicHeader()
113 114 # m_ProcessingHeader = ProcessingHeader()
114 115
115 116 systemHeaderObj = SystemHeader()
116 117
117 118 radarControllerHeaderObj = RadarControllerHeader()
118 119
119 120 # data = None
120 121
121 122 type = None
122 123
123 124 dtype = None
124 125
125 126 # nChannels = None
126 127
127 128 # nHeights = None
128 129
129 130 nProfiles = None
130 131
131 132 heightList = None
132 133
133 134 channelList = None
134 135
135 136 flagNoData = True
136 137
137 138 flagTimeBlock = False
138 139
139 140 utctime = None
140 141
141 142 blocksize = None
142 143
143 144 nCode = None
144 145
145 146 nBaud = None
146 147
147 148 code = None
148 149
149 150 flagDecodeData = True #asumo q la data esta decodificada
150 151
151 152 flagDeflipData = True #asumo q la data esta sin flip
152 153
153 154 flagShiftFFT = False
154 155
155 156 ippSeconds = None
156 157
157 158 timeInterval = None
158 159
159 160 nCohInt = None
160 161
161 162 noise = None
162 163
163 164 #Speed of ligth
164 165 C = 3e8
165 166
166 167 frequency = 49.92e6
167 168
168 169 def __init__(self):
169 170
170 171 raise ValueError, "This class has not been implemented"
171 172
172 173 def copy(self, inputObj=None):
173 174
174 175 if inputObj == None:
175 176 return copy.deepcopy(self)
176 177
177 178 for key in inputObj.__dict__.keys():
178 179 self.__dict__[key] = inputObj.__dict__[key]
179 180
180 181 def deepcopy(self):
181 182
182 183 return copy.deepcopy(self)
183 184
184 185 def isEmpty(self):
185 186
186 187 return self.flagNoData
187 188
188 189 def getNoise(self):
189 190
190 191 raise ValueError, "Not implemented"
191 192
192 193 def getNChannels(self):
193 194
194 195 return len(self.channelList)
195 196
196 197 def getChannelIndexList(self):
197 198
198 199 return range(self.nChannels)
199 200
200 201 def getNHeights(self):
201 202
202 203 return len(self.heightList)
203 204
204 205 def getHeiRange(self, extrapoints=0):
205 206
206 207 heis = self.heightList
207 208 # deltah = self.heightList[1] - self.heightList[0]
208 209 #
209 210 # heis.append(self.heightList[-1])
210 211
211 212 return heis
212 213
213 214 def getDatatime(self):
214 215
215 216 datatime = datetime.datetime.utcfromtimestamp(self.utctime)
216 217 return datatime
217 218
218 219 def getTimeRange(self):
219 220
220 221 datatime = []
221 222
222 223 datatime.append(self.utctime)
223 224 datatime.append(self.utctime + self.timeInterval)
224 225
225 226 datatime = numpy.array(datatime)
226 227
227 228 return datatime
228 229
229 230 def getFmax(self):
230 231
231 232 PRF = 1./(self.ippSeconds * self.nCohInt)
232 233
233 234 fmax = PRF/2.
234 235
235 236 return fmax
236 237
237 238 def getVmax(self):
238 239
239 240 _lambda = self.C/self.frequency
240 241
241 242 vmax = self.getFmax() * _lambda / 2.
242 243
243 244 return vmax
244 245
245 246 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
246 247 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
247 248 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
248 249 noise = property(getNoise, "I'm the 'nHeights' property.")
249 250 datatime = property(getDatatime, "I'm the 'datatime' property")
250 251
251 252 class Voltage(JROData):
252 253
253 254 #data es un numpy array de 2 dmensiones (canales, alturas)
254 255 data = None
255 256
256 257 def __init__(self):
257 258 '''
258 259 Constructor
259 260 '''
260 261
261 262 self.radarControllerHeaderObj = RadarControllerHeader()
262 263
263 264 self.systemHeaderObj = SystemHeader()
264 265
265 266 self.type = "Voltage"
266 267
267 268 self.data = None
268 269
269 270 self.dtype = None
270 271
271 272 # self.nChannels = 0
272 273
273 274 # self.nHeights = 0
274 275
275 276 self.nProfiles = None
276 277
277 278 self.heightList = None
278 279
279 280 self.channelList = None
280 281
281 282 # self.channelIndexList = None
282 283
283 284 self.flagNoData = True
284 285
285 286 self.flagTimeBlock = False
286 287
287 288 self.utctime = None
288 289
289 290 self.nCohInt = None
290 291
291 292 self.blocksize = None
292 293
293 294 def getNoisebyHildebrand(self):
294 295 """
295 296 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
296 297
297 298 Return:
298 299 noiselevel
299 300 """
300 301
301 302 for channel in range(self.nChannels):
302 303 daux = self.data_spc[channel,:,:]
303 304 self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt)
304 305
305 306 return self.noise
306 307
307 308 def getNoise(self, type = 1):
308 309
309 310 self.noise = numpy.zeros(self.nChannels)
310 311
311 312 if type == 1:
312 313 noise = self.getNoisebyHildebrand()
313 314
314 315 return 10*numpy.log10(noise)
315 316
316 317 class Spectra(JROData):
317 318
318 319 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
319 320 data_spc = None
320 321
321 322 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
322 323 data_cspc = None
323 324
324 325 #data es un numpy array de 2 dmensiones (canales, alturas)
325 326 data_dc = None
326 327
327 328 nFFTPoints = None
328 329
329 330 nPairs = None
330 331
331 332 pairsList = None
332 333
333 334 nIncohInt = None
334 335
335 336 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
336 337
337 338 nCohInt = None #se requiere para determinar el valor de timeInterval
338 339
339 340 def __init__(self):
340 341 '''
341 342 Constructor
342 343 '''
343 344
344 345 self.radarControllerHeaderObj = RadarControllerHeader()
345 346
346 347 self.systemHeaderObj = SystemHeader()
347 348
348 349 self.type = "Spectra"
349 350
350 351 # self.data = None
351 352
352 353 self.dtype = None
353 354
354 355 # self.nChannels = 0
355 356
356 357 # self.nHeights = 0
357 358
358 359 self.nProfiles = None
359 360
360 361 self.heightList = None
361 362
362 363 self.channelList = None
363 364
364 365 # self.channelIndexList = None
365 366
366 367 self.flagNoData = True
367 368
368 369 self.flagTimeBlock = False
369 370
370 371 self.utctime = None
371 372
372 373 self.nCohInt = None
373 374
374 375 self.nIncohInt = None
375 376
376 377 self.blocksize = None
377 378
378 379 self.nFFTPoints = None
379 380
380 381 self.wavelength = None
381 382
382 383 def getNoisebyHildebrand(self):
383 384 """
384 385 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
385 386
386 387 Return:
387 388 noiselevel
388 389 """
389 390
390 391 for channel in range(self.nChannels):
391 392 daux = self.data_spc[channel,:,:]
392 393 self.noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
393 394
394 395 return self.noise
395 396
396 397 def getNoisebyWindow(self, heiIndexMin=0, heiIndexMax=-1, freqIndexMin=0, freqIndexMax=-1):
397 398 """
398 399 Determina el ruido del canal utilizando la ventana indicada con las coordenadas:
399 400 (heiIndexMIn, freqIndexMin) hasta (heiIndexMax, freqIndexMAx)
400 401
401 402 Inputs:
402 403 heiIndexMin: Limite inferior del eje de alturas
403 404 heiIndexMax: Limite superior del eje de alturas
404 405 freqIndexMin: Limite inferior del eje de frecuencia
405 406 freqIndexMax: Limite supoerior del eje de frecuencia
406 407 """
407 408
408 409 data = self.data_spc[:, heiIndexMin:heiIndexMax, freqIndexMin:freqIndexMax]
409 410
410 411 for channel in range(self.nChannels):
411 412 daux = data[channel,:,:]
412 413 self.noise[channel] = numpy.average(daux)
413 414
414 415 return self.noise
415 416
416 417 def getNoisebySort(self):
417 418
418 419 for channel in range(self.nChannels):
419 420 daux = self.data_spc[channel,:,:]
420 421 self.noise[channel] = sorting_bruce(daux, self.nIncohInt)
421 422
422 423 return self.noise
423 424
424 425 def getNoise(self, type = 1):
425 426
426 427 self.noise = numpy.zeros(self.nChannels)
427 428
428 429 if type == 1:
429 430 noise = self.getNoisebyHildebrand()
430 431
431 432 if type == 2:
432 433 noise = self.getNoisebySort()
433 434
434 435 if type == 3:
435 436 noise = self.getNoisebyWindow()
436 437
437 438 return 10*numpy.log10(noise)
438 439
439 440
440 441 def getFreqRange(self, extrapoints=0):
441 442
442 443 delfreq = 2 * self.getFmax() / self.nFFTPoints
443 444 freqrange = deltafreqs*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
444 445
445 446 return freqrange
446 447
447 448 def getVelRange(self, extrapoints=0):
448 449
449 450 deltav = 2 * self.getVmax() / self.nFFTPoints
450 451 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
451 452
452 453 return velrange
453 454
454 455 def getNPairs(self):
455 456
456 457 return len(self.pairsList)
457 458
458 459 def getPairsIndexList(self):
459 460
460 461 return range(self.nPairs)
461 462
462 463 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
463 464 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
464 465
465 466 class SpectraHeis(JROData):
466 467
467 468 data_spc = None
468 469
469 470 data_cspc = None
470 471
471 472 data_dc = None
472 473
473 474 nFFTPoints = None
474 475
475 476 nPairs = None
476 477
477 478 pairsList = None
478 479
479 480 nIncohInt = None
480 481
481 482 def __init__(self):
482 483
483 484 self.radarControllerHeaderObj = RadarControllerHeader()
484 485
485 486 self.systemHeaderObj = SystemHeader()
486 487
487 488 self.type = "SpectraHeis"
488 489
489 490 self.dtype = None
490 491
491 492 # self.nChannels = 0
492 493
493 494 # self.nHeights = 0
494 495
495 496 self.nProfiles = None
496 497
497 498 self.heightList = None
498 499
499 500 self.channelList = None
500 501
501 502 # self.channelIndexList = None
502 503
503 504 self.flagNoData = True
504 505
505 506 self.flagTimeBlock = False
506 507
507 508 self.nPairs = 0
508 509
509 510 self.utctime = None
510 511
511 512 self.blocksize = None
General Comments 0
You need to be logged in to leave comments. Login now