##// 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 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 17 data = data.copy()
69 18
70 sortdata = numpy.sort(data)
71 lenOfData = len(data)
19 sortdata = numpy.sort(data,axis=None)
20 lenOfData = len(sortdata)
72 21 nums_min = lenOfData/10
73 22
74 if (lenOfData/10) > 0:
23 if (lenOfData/10) > 2:
75 24 nums_min = lenOfData/10
76 25 else:
77 nums_min = 0
78
79 rtest = 1.0 + 1.0/navg
80
26 nums_min = 2
27
81 28 sump = 0.
82 29
83 30 sumq = 0.
@@ -95,17 +42,15 def sorting_bruce(data, navg):
95 42 j += 1
96 43
97 44 if j > nums_min:
98 if ((sumq*j) <= (rtest*sump**2)):
99 lnoise = sump / j
100 else:
45 rtest = float(j)/(j-1) + 1.0/navg
46 if ((sumq*j) > (rtest*sump**2)):
101 47 j = j - 1
102 48 sump = sump - sortdata[j]
103 49 sumq = sumq - sortdata[j]**2
104 50 cont = 0
105
106 if j == nums_min:
107 lnoise = sump /j
108
51
52 lnoise = sump /j
53 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
109 54 return lnoise
110 55
111 56 class JROData:
@@ -422,6 +367,8 class Spectra(JROData):
422 367 self.flagShiftFFT = False
423 368
424 369 self.ippFactor = 1
370
371 self.noise = None
425 372
426 373 def getNoisebyHildebrand(self):
427 374 """
@@ -436,48 +383,12 class Spectra(JROData):
436 383 self.noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
437 384
438 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 387 def getNoise(self, type = 1):
469 388 if self.noise == None:
470 389 self.noise = numpy.zeros(self.nChannels)
471
472 if type == 1:
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
390 self.noise = self.getNoisebyHildebrand()
391
481 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 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 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 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 1236 minIndex = 0
1207 1237 maxIndex = 0
1208 1238 heights = self.dataOut.heightList
@@ -1226,8 +1256,22 class SpectraProc(ProcessingUnit):
1226 1256 if (maxIndex >= self.dataOut.nHeights):
1227 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 1275 noise = numpy.zeros(self.dataOut.nChannels)
1232 1276
1233 1277 for channel in range(self.dataOut.nChannels):
General Comments 0
You need to be logged in to leave comments. Login now