##// END OF EJS Templates
schain
RSS
En jrodata.py se eliminan los metodos sorting_bruce, getNoisebySort, getNoisebyWindow....
Daniel Valdez -
r450:3d2d18302dc5
parent child
Show More
Show file before | Show file after | Hide comments
@@ -11,73 +11,20 import datetime
11
11
12 from jroheaderIO import SystemHeader, RadarControllerHeader
12 from jroheaderIO import SystemHeader, RadarControllerHeader
13
13
14 def hildebrand_sekhon(data, navg):
15 """
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
18 densities is equal to the mean spectral density for white Gaussian noise
19
20 Inputs:
21 Data : heights
22 navg : numbers of averages
23
24 Return:
25 -1 : any error
26 anoise : noise's level
27 """
28
29 dataflat = data.copy().reshape(-1)
30 dataflat.sort()
31 npts = dataflat.size #numbers of points of the data
32 npts_noise = 0.2*npts
33
34 if npts < 32:
35 print "error in noise - requires at least 32 points"
36 return -1.0
37
38 dataflat2 = numpy.power(dataflat,2)
39
40 cs = numpy.cumsum(dataflat)
41 cs2 = numpy.cumsum(dataflat2)
42
43 # data sorted in ascending order
44 nmin = int((npts + 7.)/8)
45
46 for i in range(nmin, npts):
47 s = cs[i]
48 s2 = cs2[i]
49 p = s / float(i);
50 p2 = p**2;
51 q = s2 / float(i) - p2;
52 leftc = p2;
53 rightc = q * float(navg);
54 R2 = leftc/rightc
55
56 # Signal detect: R2 < 1 (R2 = leftc/rightc)
57 if R2 < 1:
58 npts_noise = i
59 break
60
61
62 anoise = numpy.average(dataflat[0:npts_noise])
63
64 return anoise;
65
14
66 def sorting_bruce(data, navg):
15 def hildebrand_sekhon(data, navg):
67
16
68 data = data.copy()
17 data = data.copy()
69
18
70 sortdata = numpy.sort(data)
19 sortdata = numpy.sort(data,axis=None)
71 lenOfData = len(data)
20 lenOfData = len(sortdata)
72 nums_min = lenOfData/10
21 nums_min = lenOfData/10
73
22
74 if (lenOfData/10) > 0:
23 if (lenOfData/10) > 2:
75 nums_min = lenOfData/10
24 nums_min = lenOfData/10
76 else:
25 else:
77 nums_min = 0
26 nums_min = 2
78
27
79 rtest = 1.0 + 1.0/navg
80
81 sump = 0.
28 sump = 0.
82
29
83 sumq = 0.
30 sumq = 0.
@@ -95,17 +42,15 def sorting_bruce(data, navg):
95 j += 1
42 j += 1
96
43
97 if j > nums_min:
44 if j > nums_min:
98 if ((sumq*j) <= (rtest*sump**2)):
45 rtest = float(j)/(j-1) + 1.0/navg
99 lnoise = sump / j
46 if ((sumq*j) > (rtest*sump**2)):
100 else:
101 j = j - 1
47 j = j - 1
102 sump = sump - sortdata[j]
48 sump = sump - sortdata[j]
103 sumq = sumq - sortdata[j]**2
49 sumq = sumq - sortdata[j]**2
104 cont = 0
50 cont = 0
105
51
106 if j == nums_min:
52 lnoise = sump /j
107 lnoise = sump /j
53 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
108
109 return lnoise
54 return lnoise
110
55
111 class JROData:
56 class JROData:
@@ -422,6 +367,8 class Spectra(JROData):
422 self.flagShiftFFT = False
367 self.flagShiftFFT = False
423
368
424 self.ippFactor = 1
369 self.ippFactor = 1
370
371 self.noise = None
425
372
426 def getNoisebyHildebrand(self):
373 def getNoisebyHildebrand(self):
427 """
374 """
@@ -436,48 +383,12 class Spectra(JROData):
436 self.noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
383 self.noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
437
384
438 return self.noise
385 return self.noise
439
440 def getNoisebyWindow(self, heiIndexMin=0, heiIndexMax=-1, freqIndexMin=0, freqIndexMax=-1):
441 """
442 Determina el ruido del canal utilizando la ventana indicada con las coordenadas:
443 (heiIndexMIn, freqIndexMin) hasta (heiIndexMax, freqIndexMAx)
444
445 Inputs:
446 heiIndexMin: Limite inferior del eje de alturas
447 heiIndexMax: Limite superior del eje de alturas
448 freqIndexMin: Limite inferior del eje de frecuencia
449 freqIndexMax: Limite supoerior del eje de frecuencia
450 """
451
452 data = self.data_spc[:, heiIndexMin:heiIndexMax, freqIndexMin:freqIndexMax]
453
386
454 for channel in range(self.nChannels):
455 daux = data[channel,:,:]
456 self.noise[channel] = numpy.average(daux)
457
458 return self.noise
459
460 def getNoisebySort(self):
461
462 for channel in range(self.nChannels):
463 daux = self.data_spc[channel,:,:]
464 self.noise[channel] = sorting_bruce(daux, self.nIncohInt)
465
466 return self.noise
467
468 def getNoise(self, type = 1):
387 def getNoise(self, type = 1):
469 if self.noise == None:
388 if self.noise == None:
470 self.noise = numpy.zeros(self.nChannels)
389 self.noise = numpy.zeros(self.nChannels)
471
390 self.noise = self.getNoisebyHildebrand()
472 if type == 1:
391
473 self.noise = self.getNoisebyHildebrand()
474
475 if type == 2:
476 self.noise = self.getNoisebySort()
477
478 if type == 3:
479 self.noise = self.getNoisebyWindow()
480
481 return self.noise
392 return self.noise
482
393
483
394
Show file before | Show file after | Hide comments
@@ -1195,14 +1195,44 class SpectraProc(ProcessingUnit):
1195
1195
1196 return 1
1196 return 1
1197
1197
1198 def getNoise(self, minHei, maxHei):
1198 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
1199 #validacion de rango
1200 if minHei == None:
1201 minHei = self.dataOut.heightList[0]
1199
1202
1203 if maxHei == None:
1204 maxHei = self.dataOut.heightList[-1]
1205
1200 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
1206 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
1201 raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
1207 print 'minHei: %.2f is out of the heights range'%(minHei)
1208 print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
1209 minHei = self.dataOut.heightList[0]
1202
1210
1203 if (maxHei > self.dataOut.heightList[-1]):
1211 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
1212 print 'maxHei: %.2f is out of the heights range'%(maxHei)
1213 print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
1204 maxHei = self.dataOut.heightList[-1]
1214 maxHei = self.dataOut.heightList[-1]
1205
1215
1216 # validacion de velocidades
1217 velrange = self.dataOut.getVelRange(1)
1218
1219 if minVel == None:
1220 minVel = velrange[0]
1221
1222 if maxVel == None:
1223 maxVel = velrange[-1]
1224
1225 if (minVel < velrange[0]) or (minVel > maxVel):
1226 print 'minVel: %.2f is out of the velocity range'%(minVel)
1227 print 'minVel is setting to %.2f'%(velrange[0])
1228 minVel = velrange[0]
1229
1230 if (maxVel > velrange[-1]) or (maxVel < minVel):
1231 print 'maxVel: %.2f is out of the velocity range'%(maxVel)
1232 print 'maxVel is setting to %.2f'%(velrange[-1])
1233 maxVel = velrange[-1]
1234
1235 # seleccion de indices para rango
1206 minIndex = 0
1236 minIndex = 0
1207 maxIndex = 0
1237 maxIndex = 0
1208 heights = self.dataOut.heightList
1238 heights = self.dataOut.heightList
@@ -1226,8 +1256,22 class SpectraProc(ProcessingUnit):
1226 if (maxIndex >= self.dataOut.nHeights):
1256 if (maxIndex >= self.dataOut.nHeights):
1227 maxIndex = self.dataOut.nHeights-1
1257 maxIndex = self.dataOut.nHeights-1
1228
1258
1229 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
1259 # seleccion de indices para velocidades
1260 indminvel = numpy.where(velrange >= minVel)
1261 indmaxvel = numpy.where(velrange <= maxVel)
1262 try:
1263 minIndexVel = indminvel[0][0]
1264 except:
1265 minIndexVel = 0
1266
1267 try:
1268 maxIndexVel = indmaxvel[0][-1]
1269 except:
1270 maxIndexVel = len(velrange)
1230
1271
1272 #seleccion del espectro
1273 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
1274 #estimacion de ruido
1231 noise = numpy.zeros(self.dataOut.nChannels)
1275 noise = numpy.zeros(self.dataOut.nChannels)
1232
1276
1233 for channel in range(self.dataOut.nChannels):
1277 for channel in range(self.dataOut.nChannels):
General Comments 0
You need to be logged in to leave comments. Login now