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