##// END OF EJS Templates
schain-jc
RSS
23/11/2017
ebocanegra -
r1123:72ed20327727
parent child
Show More
Side by Side Unified
Context file:
r1123:72ed20327727 Loading diff...:
Show file before | Show file after | Hide comments
@@ -706,7 +706,7 class ProcUnitConf():
706 #print "\tRunning the '%s' operation with %s" %(opConfObj.name, opConfObj.id)
706 #print "\tRunning the '%s' operation with %s" %(opConfObj.name, opConfObj.id)
707 sts = self.procUnitObj.call(opType = opConfObj.type,
707 sts = self.procUnitObj.call(opType = opConfObj.type,
708 opName = opConfObj.name,
708 opName = opConfObj.name,
709 opId = opConfObj.id,
709 opId = opConfObj.id
710 )
710 )
711
711
712 # total_time = time.time() - ini
712 # total_time = time.time() - ini
Show file before | Show file after | Hide comments
@@ -82,7 +82,7 class FitGauPlot(Figure):
82 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
82 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
83 server=None, folder=None, username=None, password=None,
83 server=None, folder=None, username=None, password=None,
84 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
84 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
85 xaxis="frequency", colormap='jet', normFactor=None , GauSelector = 1):
85 xaxis="frequency", colormap='jet', normFactor=None , GauSelector = 0):
86
86
87 """
87 """
88
88
@@ -125,7 +125,7 class FitGauPlot(Figure):
125 x = dataOut.spc_range[1]
125 x = dataOut.spc_range[1]
126 xlabel = "Time (ms)"
126 xlabel = "Time (ms)"
127
127
128 else:
128 else:
129 x = dataOut.spc_range[2]
129 x = dataOut.spc_range[2]
130 xlabel = "Velocity (m/s)"
130 xlabel = "Velocity (m/s)"
131
131
@@ -667,13 +667,14 class WindProfilerPlot(Figure):
667 #If there is a SNR function defined
667 #If there is a SNR function defined
668 if dataOut.data_SNR is not None:
668 if dataOut.data_SNR is not None:
669 nplots += 1
669 nplots += 1
670 SNR = dataOut.data_SNR
670 SNR = dataOut.data_SNR[0]
671 SNRavg = numpy.average(SNR, axis=0)
671 SNRavg = SNR#numpy.average(SNR, axis=0)
672
672
673 SNRdB = 10*numpy.log10(SNR)
673 SNRdB = 10*numpy.log10(SNR)
674 SNRavgdB = 10*numpy.log10(SNRavg)
674 SNRavgdB = 10*numpy.log10(SNRavg)
675
675
676 if SNRthresh == None: SNRthresh = -5.0
676 if SNRthresh == None:
677 SNRthresh = -5.0
677 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
678 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
678
679
679 for i in range(nplotsw):
680 for i in range(nplotsw):
Show file before | Show file after | Hide comments
@@ -27,8 +27,8 class SpectraPlot(Figure):
27 self.isConfig = False
27 self.isConfig = False
28 self.__nsubplots = 1
28 self.__nsubplots = 1
29
29
30 self.WIDTH = 250
30 self.WIDTH = 300
31 self.HEIGHT = 250
31 self.HEIGHT = 300
32 self.WIDTHPROF = 120
32 self.WIDTHPROF = 120
33 self.HEIGHTPROF = 0
33 self.HEIGHTPROF = 0
34 self.counter_imagwr = 0
34 self.counter_imagwr = 0
@@ -128,6 +128,10 class SpectraPlot(Figure):
128 factor = normFactor
128 factor = normFactor
129 if xaxis == "frequency":
129 if xaxis == "frequency":
130 x = dataOut.getFreqRange(1)/1000.
130 x = dataOut.getFreqRange(1)/1000.
131 print '#######################################################'
132 print 'xlen', len(x)
133 print x
134 print '#######################################################'
131 xlabel = "Frequency (kHz)"
135 xlabel = "Frequency (kHz)"
132
136
133 elif xaxis == "time":
137 elif xaxis == "time":
@@ -439,24 +443,45 class CrossSpectraPlot(Figure):
439 xlabel=xlabel, ylabel=ylabel, title=title,
443 xlabel=xlabel, ylabel=ylabel, title=title,
440 ticksize=9, colormap=power_cmap, cblabel='')
444 ticksize=9, colormap=power_cmap, cblabel='')
441
445
442 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:])
446 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:] / numpy.sqrt( dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:] )
443 coherence = numpy.abs(coherenceComplex)
447 coherence = numpy.abs(coherenceComplex)
444 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
448 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
445 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
449 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
446
450
447
451
448 # print 'FASE', numpy.shape(phase), y[10]
452
453
454 # #print 'FASE', numpy.shape(phase), y[25]
449 # fig = plt.figure(10+self.indice)
455 # fig = plt.figure(10+self.indice)
450 # plt.plot( x[0:128],phase[:,10] )
456 # #plt.plot( x[0:256],coherence[:,25] )
457 # cohAv = numpy.average(coherence,1)
458 #
459 # plt.plot( x[0:256],cohAv )
451 # #plt.axis([-12, 12, 15, 50])
460 # #plt.axis([-12, 12, 15, 50])
452 # plt.title("%s" %( '%s %s, Channel %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S") , i)))
461 # plt.title("%s" %( '%s %s, Channel %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S") , i)))
453 # plt.ylabel('Desfase [grados]')
462 # plt.ylabel('Desfase [grados]')
454 # plt.xlabel('Velocidad [m/s]')
463 # plt.xlabel('Velocidad [m/s]')
455 # fig.savefig('/home/erick/Documents/Pics/to{}.png'.format(self.indice))
464 # fig.savefig('/home/erick/Documents/Pics/to{}.png'.format(self.indice))
456 #
465 #
466 # plt.show()
467 # self.indice=self.indice+1
468
469
470 # print 'FASE', numpy.shape(phase), y[25]
471 # fig = plt.figure(10+self.indice)
472 # plt.plot( x[0:256],phase[:,25] )
473 # #plt.axis([-12, 12, 15, 50])
474 # plt.title("%s" %( '%s %s, Channel %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S") , i)))
475 # plt.ylabel('Desfase [grados]')
476 # plt.xlabel('Velocidad [m/s]')
477 # fig.savefig('/home/erick/Documents/Pics/to{}.png'.format(self.indice))
478 #
457 # plt.show()
479 # plt.show()
458 # self.indice=self.indice+1
480 # self.indice=self.indice+1
459
481
482
483
484
460 title = "Coherence Ch%d * Ch%d" %(pair[0], pair[1])
485 title = "Coherence Ch%d * Ch%d" %(pair[0], pair[1])
461 axes0 = self.axesList[i*self.__nsubplots+2]
486 axes0 = self.axesList[i*self.__nsubplots+2]
462 axes0.pcolor(x, y, coherence,
487 axes0.pcolor(x, y, coherence,
@@ -500,7 +525,7 class RTIPlot(Figure):
500 self.__nsubplots = 1
525 self.__nsubplots = 1
501
526
502 self.WIDTH = 800
527 self.WIDTH = 800
503 self.HEIGHT = 180
528 self.HEIGHT = 250
504 self.WIDTHPROF = 120
529 self.WIDTHPROF = 120
505 self.HEIGHTPROF = 0
530 self.HEIGHTPROF = 0
506 self.counter_imagwr = 0
531 self.counter_imagwr = 0
Show file before | Show file after | Hide comments
@@ -110,6 +110,7 class BLTRParamReader(JRODataReader, ProcessingUnit):
110 status_value=0,
110 status_value=0,
111 **kwargs):
111 **kwargs):
112
112
113 self.verbose = kwargs.get('verbose', True)
113 self.path = path
114 self.path = path
114 self.startTime = startTime
115 self.startTime = startTime
115 self.endTime = endTime
116 self.endTime = endTime
@@ -214,17 +215,19 class BLTRParamReader(JRODataReader, ProcessingUnit):
214 self.readBlock()
215 self.readBlock()
215
216
216 if (self.datatime.time() < self.startTime) or (self.datatime.time() > self.endTime):
217 if (self.datatime.time() < self.startTime) or (self.datatime.time() > self.endTime):
217 print "[Reading] Record No. %d/%d -> %s [Skipping]" %(
218 if self.verbose:
218 self.counter_records,
219 print "[Reading] Record No. %d/%d -> %s [Skipping]" %(
219 self.nrecords,
220 self.counter_records,
220 self.datatime.ctime())
221 self.nrecords,
222 self.datatime.ctime())
221 continue
223 continue
222 break
224 break
223
225
224 print "[Reading] Record No. %d/%d -> %s" %(
226 if self.verbose:
225 self.counter_records,
227 print "[Reading] Record No. %d/%d -> %s" %(
226 self.nrecords,
228 self.counter_records,
227 self.datatime.ctime())
229 self.nrecords,
230 self.datatime.ctime())
228
231
229 return 1
232 return 1
230
233
Show file before | Show file after | Hide comments
@@ -1784,7 +1784,7 class JRODataWriter(JRODataIO):
1784
1784
1785 return 1
1785 return 1
1786
1786
1787 def run(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1787 def run(self, dataOut, path, blocksPerFile=100, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1788
1788
1789 if not(self.isConfig):
1789 if not(self.isConfig):
1790
1790
Show file before | Show file after | Hide comments
@@ -300,6 +300,7 class SpectraReader(JRODataReader, ProcessingUnit):
300 self.data_cspc = cspc['real'] + cspc['imag']*1j
300 self.data_cspc = cspc['real'] + cspc['imag']*1j
301 else:
301 else:
302 self.data_cspc = None
302 self.data_cspc = None
303 print 'SALE NONE ***********************************************************'
303
304
304
305
305 if self.processingHeaderObj.flag_dc:
306 if self.processingHeaderObj.flag_dc:
@@ -434,7 +435,7 class SpectraWriter(JRODataWriter, Operation):
434
435
435 # dataOut = None
436 # dataOut = None
436
437
437 def __init__(self):
438 def __init__(self, **kwargs):
438 """
439 """
439 Inicializador de la clase SpectraWriter para la escritura de datos de espectros.
440 Inicializador de la clase SpectraWriter para la escritura de datos de espectros.
440
441
@@ -449,7 +450,7 class SpectraWriter(JRODataWriter, Operation):
449 Return: None
450 Return: None
450 """
451 """
451
452
452 Operation.__init__(self)
453 Operation.__init__(self, **kwargs)
453
454
454 self.isConfig = False
455 self.isConfig = False
455
456
@@ -518,7 +519,7 class SpectraWriter(JRODataWriter, Operation):
518
519
519
520
520 def writeBlock(self):
521 def writeBlock(self):
521 """
522 """processingHeaderObj
522 Escribe el buffer en el file designado
523 Escribe el buffer en el file designado
523
524
524
525
@@ -542,8 +543,10 class SpectraWriter(JRODataWriter, Operation):
542 data.tofile(self.fp)
543 data.tofile(self.fp)
543
544
544 if self.data_cspc is not None:
545 if self.data_cspc is not None:
545 data = numpy.zeros( self.shape_cspc_Buffer, self.dtype )
546
546 cspc = numpy.transpose( self.data_cspc, (0,2,1) )
547 cspc = numpy.transpose( self.data_cspc, (0,2,1) )
548 data = numpy.zeros( numpy.shape(cspc), self.dtype )
549 print 'data.shape', self.shape_cspc_Buffer
547 if not( self.processingHeaderObj.shif_fft ):
550 if not( self.processingHeaderObj.shif_fft ):
548 cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
551 cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
549 data['real'] = cspc.real
552 data['real'] = cspc.real
@@ -553,8 +556,9 class SpectraWriter(JRODataWriter, Operation):
553
556
554
557
555 if self.data_dc is not None:
558 if self.data_dc is not None:
556 data = numpy.zeros( self.shape_dc_Buffer, self.dtype )
559
557 dc = self.data_dc
560 dc = self.data_dc
561 data = numpy.zeros( numpy.shape(dc), self.dtype )
558 data['real'] = dc.real
562 data['real'] = dc.real
559 data['imag'] = dc.imag
563 data['imag'] = dc.imag
560 data = data.reshape((-1))
564 data = data.reshape((-1))
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (934 lines changed) Show them Hide them
@@ -52,13 +52,6 def _unpickle_method(func_name, obj, cls):
52 break
52 break
53 return func.__get__(obj, cls)
53 return func.__get__(obj, cls)
54
54
55
56
57
58
59
60
61
62 class ParametersProc(ProcessingUnit):
55 class ParametersProc(ProcessingUnit):
63
56
64 nSeconds = None
57 nSeconds = None
@@ -125,11 +118,11 class ParametersProc(ProcessingUnit):
125
118
126 if self.dataIn.type == "Spectra":
119 if self.dataIn.type == "Spectra":
127
120
128 self.dataOut.data_pre = (self.dataIn.data_spc,self.dataIn.data_cspc)
121 self.dataOut.data_pre = (self.dataIn.data_spc , self.dataIn.data_cspc)
129 print 'self.dataIn.data_spc', self.dataIn.data_spc.shape
122 print 'self.dataIn.data_spc', self.dataIn.data_spc.shape
130 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
123 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
131 self.dataOut.spc_noise = self.dataIn.getNoise()
124 self.dataOut.spc_noise = self.dataIn.getNoise()
132 self.dataOut.spc_range = (self.dataIn.getFreqRange(1)/1000. , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1) )
125 self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1) )
133
126
134 self.dataOut.normFactor = self.dataIn.normFactor
127 self.dataOut.normFactor = self.dataIn.normFactor
135 #self.dataOut.outputInterval = self.dataIn.outputInterval
128 #self.dataOut.outputInterval = self.dataIn.outputInterval
@@ -142,9 +135,9 class ParametersProc(ProcessingUnit):
142
135
143 print 'datain chandist ',self.dataOut.ChanDist
136 print 'datain chandist ',self.dataOut.ChanDist
144
137
145 if hasattr(self.dataIn, 'VelRange'): #Velocities range
138 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
146 self.dataOut.VelRange = self.dataIn.VelRange
139 # self.dataOut.VelRange = self.dataIn.VelRange
147 else: self.dataOut.VelRange = None
140 #else: self.dataOut.VelRange = None
148
141
149 if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
142 if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
150 self.dataOut.RadarConst = self.dataIn.RadarConst
143 self.dataOut.RadarConst = self.dataIn.RadarConst
@@ -193,6 +186,7 def target(tups):
193 #print 'TARGETTT', obj, args
186 #print 'TARGETTT', obj, args
194 return obj.FitGau(args)
187 return obj.FitGau(args)
195
188
189
196 class GaussianFit(Operation):
190 class GaussianFit(Operation):
197
191
198 '''
192 '''
@@ -212,7 +206,7 class GaussianFit(Operation):
212 self.i=0
206 self.i=0
213
207
214
208
215 def run(self, dataOut, num_intg=7, pnoise=1., vel_arr=None, SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
209 def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
216 """This routine will find a couple of generalized Gaussians to a power spectrum
210 """This routine will find a couple of generalized Gaussians to a power spectrum
217 input: spc
211 input: spc
218 output:
212 output:
@@ -239,12 +233,9 class GaussianFit(Operation):
239 #self.Num_Bin = len(self.spc)
233 #self.Num_Bin = len(self.spc)
240 self.Num_Bin = self.spc.shape[1]
234 self.Num_Bin = self.spc.shape[1]
241 self.Num_Chn = self.spc.shape[0]
235 self.Num_Chn = self.spc.shape[0]
242
243 Vrange = dataOut.abscissaList
236 Vrange = dataOut.abscissaList
244
237
245 #print 'self.spc2', numpy.asarray(self.spc).shape
238 GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei])
246
247 GauSPC = numpy.empty([2,self.Num_Bin,self.Num_Hei])
248 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
239 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
249 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
240 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
250 SPC_ch1[:] = numpy.NaN
241 SPC_ch1[:] = numpy.NaN
@@ -256,264 +247,14 class GaussianFit(Operation):
256 noise_ = dataOut.spc_noise[0].copy()
247 noise_ = dataOut.spc_noise[0].copy()
257
248
258
249
259
260 pool = Pool(processes=self.Num_Chn)
250 pool = Pool(processes=self.Num_Chn)
261 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
251 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
262 objs = [self for __ in range(self.Num_Chn)]
252 objs = [self for __ in range(self.Num_Chn)]
263 attrs = zip(objs, args)
253 attrs = zip(objs, args)
264 gauSPC = pool.map(target, attrs)
254 gauSPC = pool.map(target, attrs)
265 dataOut.GauSPC = numpy.asarray(gauSPC)
255 dataOut.GauSPC = numpy.asarray(gauSPC)
266 # ret = []
256
267 # for n in range(self.Num_Chn):
257
268 # self.FitGau(args[n])
269 # dataOut.GauSPC = ret
270
271
272
273 # for ch in range(self.Num_Chn):
274 #
275 # for ht in range(self.Num_Hei):
276 # #print (numpy.asarray(self.spc).shape)
277 # spc = numpy.asarray(self.spc)[ch,:,ht]
278 #
279 # #############################################
280 # # normalizing spc and noise
281 # # This part differs from gg1
282 # spc_norm_max = max(spc)
283 # spc = spc / spc_norm_max
284 # pnoise = pnoise / spc_norm_max
285 # #############################################
286 #
287 # if abs(vel_arr[0])<15.0: # this switch is for spectra collected with different length IPP's
288 # fatspectra=1.0
289 # else:
290 # fatspectra=0.5
291 #
292 # wnoise = noise_ / spc_norm_max
293 # #print 'wnoise', noise_, dataOut.spc_noise[0], wnoise
294 # #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
295 # #if wnoise>1.1*pnoise: # to be tested later
296 # # wnoise=pnoise
297 # noisebl=wnoise*0.9; noisebh=wnoise*1.1
298 # spc=spc-wnoise
299 #
300 # minx=numpy.argmin(spc)
301 # spcs=numpy.roll(spc,-minx)
302 # cum=numpy.cumsum(spcs)
303 # tot_noise=wnoise * self.Num_Bin #64;
304 # #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
305 # #snr=tot_signal/tot_noise
306 # #snr=cum[-1]/tot_noise
307 #
308 # #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
309 #
310 # snr = sum(spcs)/tot_noise
311 # snrdB=10.*numpy.log10(snr)
312 #
313 # #if snrdB < -9 :
314 # # snrdB = numpy.NaN
315 # # continue
316 #
317 # #print 'snr',snrdB # , sum(spcs) , tot_noise
318 #
319 #
320 # #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
321 # # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
322 #
323 # cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
324 # cumlo=cummax*epsi;
325 # cumhi=cummax*(1-epsi)
326 # powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
327 #
328 # #if len(powerindex)==1:
329 # ##return [numpy.mod(powerindex[0]+minx,64),None,None,None,],[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
330 # #return [numpy.mod(powerindex[0]+minx, self.Num_Bin ),None,None,None,],[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
331 # #elif len(powerindex)<4*fatspectra:
332 # #return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
333 #
334 # if len(powerindex) < 1:# case for powerindex 0
335 # continue
336 # powerlo=powerindex[0]
337 # powerhi=powerindex[-1]
338 # powerwidth=powerhi-powerlo
339 #
340 # firstpeak=powerlo+powerwidth/10.# first gaussian energy location
341 # secondpeak=powerhi-powerwidth/10.#second gaussian energy location
342 # midpeak=(firstpeak+secondpeak)/2.
343 # firstamp=spcs[int(firstpeak)]
344 # secondamp=spcs[int(secondpeak)]
345 # midamp=spcs[int(midpeak)]
346 # #x=numpy.spc.shape[1]
347 #
348 # #x=numpy.arange(64)
349 # x=numpy.arange( self.Num_Bin )
350 # y_data=spc+wnoise
351 #
352 # # single gaussian
353 # #shift0=numpy.mod(midpeak+minx,64)
354 # shift0=numpy.mod(midpeak+minx, self.Num_Bin )
355 # width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
356 # power0=2.
357 # amplitude0=midamp
358 # state0=[shift0,width0,amplitude0,power0,wnoise]
359 # #bnds=((0,63),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
360 # bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
361 # #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(0.1,0.5))
362 # # bnds = range of fft, power width, amplitude, power, noise
363 # lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
364 #
365 # chiSq1=lsq1[1];
366 # jack1= self.y_jacobian1(x,lsq1[0])
367 #
368 #
369 # try:
370 # sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
371 # except:
372 # std1=32.; sigmas1=numpy.ones(5)
373 # else:
374 # std1=sigmas1[0]
375 #
376 #
377 # if fatspectra<1.0 and powerwidth<4:
378 # choice=0
379 # Amplitude0=lsq1[0][2]
380 # shift0=lsq1[0][0]
381 # width0=lsq1[0][1]
382 # p0=lsq1[0][3]
383 # Amplitude1=0.
384 # shift1=0.
385 # width1=0.
386 # p1=0.
387 # noise=lsq1[0][4]
388 # #return (numpy.array([shift0,width0,Amplitude0,p0]),
389 # # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
390 #
391 # # two gaussians
392 # #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
393 # shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
394 # shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
395 # width0=powerwidth/6.;
396 # width1=width0
397 # power0=2.;
398 # power1=power0
399 # amplitude0=firstamp;
400 # amplitude1=secondamp
401 # state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
402 # #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
403 # bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
404 # #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
405 #
406 # lsq2=fmin_l_bfgs_b(self.misfit2,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
407 #
408 #
409 # chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
410 #
411 #
412 # try:
413 # sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
414 # except:
415 # std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
416 # else:
417 # std2a=sigmas2[0]; std2b=sigmas2[4]
418 #
419 #
420 #
421 # oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
422 #
423 # if snrdB>-9: # when SNR is strong pick the peak with least shift (LOS velocity) error
424 # if oneG:
425 # choice=0
426 # else:
427 # w1=lsq2[0][1]; w2=lsq2[0][5]
428 # a1=lsq2[0][2]; a2=lsq2[0][6]
429 # p1=lsq2[0][3]; p2=lsq2[0][7]
430 # s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
431 # gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
432 #
433 # if gp1>gp2:
434 # if a1>0.7*a2:
435 # choice=1
436 # else:
437 # choice=2
438 # elif gp2>gp1:
439 # if a2>0.7*a1:
440 # choice=2
441 # else:
442 # choice=1
443 # else:
444 # choice=numpy.argmax([a1,a2])+1
445 # #else:
446 # #choice=argmin([std2a,std2b])+1
447 #
448 # else: # with low SNR go to the most energetic peak
449 # choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
450 #
451 # #print 'choice',choice
452 #
453 # if choice==0: # pick the single gaussian fit
454 # Amplitude0=lsq1[0][2]
455 # shift0=lsq1[0][0]
456 # width0=lsq1[0][1]
457 # p0=lsq1[0][3]
458 # Amplitude1=0.
459 # shift1=0.
460 # width1=0.
461 # p1=0.
462 # noise=lsq1[0][4]
463 # elif choice==1: # take the first one of the 2 gaussians fitted
464 # Amplitude0 = lsq2[0][2]
465 # shift0 = lsq2[0][0]
466 # width0 = lsq2[0][1]
467 # p0 = lsq2[0][3]
468 # Amplitude1 = lsq2[0][6] # This is 0 in gg1
469 # shift1 = lsq2[0][4] # This is 0 in gg1
470 # width1 = lsq2[0][5] # This is 0 in gg1
471 # p1 = lsq2[0][7] # This is 0 in gg1
472 # noise = lsq2[0][8]
473 # else: # the second one
474 # Amplitude0 = lsq2[0][6]
475 # shift0 = lsq2[0][4]
476 # width0 = lsq2[0][5]
477 # p0 = lsq2[0][7]
478 # Amplitude1 = lsq2[0][2] # This is 0 in gg1
479 # shift1 = lsq2[0][0] # This is 0 in gg1
480 # width1 = lsq2[0][1] # This is 0 in gg1
481 # p1 = lsq2[0][3] # This is 0 in gg1
482 # noise = lsq2[0][8]
483 #
484 # #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0)
485 # SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
486 # SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
487 # #print 'SPC_ch1.shape',SPC_ch1.shape
488 # #print 'SPC_ch2.shape',SPC_ch2.shape
489 # #dataOut.data_param = SPC_ch1
490 # GauSPC[0] = SPC_ch1
491 # GauSPC[1] = SPC_ch2
492
493 # #plt.gcf().clear()
494 # plt.figure(50+self.i)
495 # self.i=self.i+1
496 # #plt.subplot(121)
497 # plt.plot(self.spc,'k')#,label='spc(66)')
498 # plt.plot(SPC_ch1[ch,ht],'b')#,label='gg1')
499 # #plt.plot(SPC_ch2,'r')#,label='gg2')
500 # #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
501 # #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
502 # #plt.plot(xFrec,FitGauss,'yo:',label='fit')
503 # plt.legend()
504 # plt.title('DATOS A ALTURA DE 7500 METROS')
505 # plt.show()
506 # print 'shift0', shift0
507 # print 'Amplitude0', Amplitude0
508 # print 'width0', width0
509 # print 'p0', p0
510 # print '========================'
511 # print 'shift1', shift1
512 # print 'Amplitude1', Amplitude1
513 # print 'width1', width1
514 # print 'p1', p1
515 # print 'noise', noise
516 # print 's_noise', wnoise
517
258
518 print '========================================================'
259 print '========================================================'
519 print 'total_time: ', time.time()-start_time
260 print 'total_time: ', time.time()-start_time
@@ -560,20 +301,22 class GaussianFit(Operation):
560 # normalizing spc and noise
301 # normalizing spc and noise
561 # This part differs from gg1
302 # This part differs from gg1
562 spc_norm_max = max(spc)
303 spc_norm_max = max(spc)
563 spc = spc / spc_norm_max
304 #spc = spc / spc_norm_max
564 pnoise = pnoise / spc_norm_max
305 pnoise = pnoise #/ spc_norm_max
565 #############################################
306 #############################################
566
307
567 fatspectra=1.0
308 fatspectra=1.0
568
309
569 wnoise = noise_ / spc_norm_max
310 wnoise = noise_ #/ spc_norm_max
570 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
311 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
571 #if wnoise>1.1*pnoise: # to be tested later
312 #if wnoise>1.1*pnoise: # to be tested later
572 # wnoise=pnoise
313 # wnoise=pnoise
573 noisebl=wnoise*0.9; noisebh=wnoise*1.1
314 noisebl=wnoise*0.9;
315 noisebh=wnoise*1.1
574 spc=spc-wnoise
316 spc=spc-wnoise
575 # print 'wnoise', noise_[0], spc_norm_max, wnoise
317 # print 'wnoise', noise_[0], spc_norm_max, wnoise
576 minx=numpy.argmin(spc)
318 minx=numpy.argmin(spc)
319 #spcs=spc.copy()
577 spcs=numpy.roll(spc,-minx)
320 spcs=numpy.roll(spc,-minx)
578 cum=numpy.cumsum(spcs)
321 cum=numpy.cumsum(spcs)
579 tot_noise=wnoise * self.Num_Bin #64;
322 tot_noise=wnoise * self.Num_Bin #64;
@@ -597,7 +340,8 class GaussianFit(Operation):
597 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
340 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
598 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
341 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
599
342
600 cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
343 cummax=max(cum);
344 epsi=0.08*fatspectra # cumsum to narrow down the energy region
601 cumlo=cummax*epsi;
345 cumlo=cummax*epsi;
602 cumhi=cummax*(1-epsi)
346 cumhi=cummax*(1-epsi)
603 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
347 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
@@ -619,7 +363,7 class GaussianFit(Operation):
619 x=numpy.arange( self.Num_Bin )
363 x=numpy.arange( self.Num_Bin )
620 y_data=spc+wnoise
364 y_data=spc+wnoise
621
365
622 # single gaussian
366 ''' single Gaussian '''
623 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
367 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
624 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
368 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
625 power0=2.
369 power0=2.
@@ -629,16 +373,7 class GaussianFit(Operation):
629 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
373 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
630
374
631 chiSq1=lsq1[1];
375 chiSq1=lsq1[1];
632 jack1= self.y_jacobian1(x,lsq1[0])
376
633
634
635 try:
636 sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
637 except:
638 std1=32.; sigmas1=numpy.ones(5)
639 else:
640 std1=sigmas1[0]
641
642
377
643 if fatspectra<1.0 and powerwidth<4:
378 if fatspectra<1.0 and powerwidth<4:
644 choice=0
379 choice=0
@@ -654,7 +389,7 class GaussianFit(Operation):
654 #return (numpy.array([shift0,width0,Amplitude0,p0]),
389 #return (numpy.array([shift0,width0,Amplitude0,p0]),
655 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
390 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
656
391
657 # two gaussians
392 ''' two gaussians '''
658 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
393 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
659 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
394 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
660 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
395 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
@@ -669,24 +404,16 class GaussianFit(Operation):
669 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
404 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
670 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
405 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
671
406
672 lsq2=fmin_l_bfgs_b(self.misfit2,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
407 lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
673
408
674
409
675 chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
410 chiSq2=lsq2[1];
411
676
412
677
413
678 try:
679 sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
680 except:
681 std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
682 else:
683 std2a=sigmas2[0]; std2b=sigmas2[4]
684
685
686
687 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
414 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
688
415
689 if snrdB>-6: # when SNR is strong pick the peak with least shift (LOS velocity) error
416 if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
690 if oneG:
417 if oneG:
691 choice=0
418 choice=0
692 else:
419 else:
@@ -696,7 +423,7 class GaussianFit(Operation):
696 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
423 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
697 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
424 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
698 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
425 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
699
426
700 if gp1>gp2:
427 if gp1>gp2:
701 if a1>0.7*a2:
428 if a1>0.7*a2:
702 choice=1
429 choice=1
@@ -716,13 +443,15 class GaussianFit(Operation):
716 choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
443 choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
717
444
718
445
719 shift0=lsq2[0][0]; vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
446 shift0=lsq2[0][0];
720 shift1=lsq2[0][4]; vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
447 vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
448 shift1=lsq2[0][4];
449 vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
721
450
722 max_vel = 20
451 max_vel = 1.0
723
452
724 #first peak will be 0, second peak will be 1
453 #first peak will be 0, second peak will be 1
725 if vel0 > 0 and vel0 < max_vel : #first peak is in the correct range
454 if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range
726 shift0=lsq2[0][0]
455 shift0=lsq2[0][0]
727 width0=lsq2[0][1]
456 width0=lsq2[0][1]
728 Amplitude0=lsq2[0][2]
457 Amplitude0=lsq2[0][2]
@@ -745,12 +474,12 class GaussianFit(Operation):
745 p0=lsq2[0][7]
474 p0=lsq2[0][7]
746 noise=lsq2[0][8]
475 noise=lsq2[0][8]
747
476
748 if Amplitude0<0.1: # in case the peak is noise
477 if Amplitude0<0.05: # in case the peak is noise
749 shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
478 shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN]
750 if Amplitude1<0.1:
479 if Amplitude1<0.05:
751 shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
480 shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN]
752
481
753
482
754 # if choice==0: # pick the single gaussian fit
483 # if choice==0: # pick the single gaussian fit
755 # Amplitude0=lsq1[0][2]
484 # Amplitude0=lsq1[0][2]
756 # shift0=lsq1[0][0]
485 # shift0=lsq1[0][0]
@@ -802,63 +531,8 class GaussianFit(Operation):
802 # print 'p1', p1
531 # print 'p1', p1
803 # print 'noise', noise
532 # print 'noise', noise
804 # print 's_noise', wnoise
533 # print 's_noise', wnoise
805
806 return GauSPC
807
808
809 def y_jacobian1(self,x,state): # This function is for further analysis of generalized Gaussians, it is not too importan for the signal discrimination.
810 y_model=self.y_model1(x,state)
811 s0,w0,a0,p0,n=state
812 e0=((x-s0)/w0)**2;
813
814 e0u=((x-s0-self.Num_Bin)/w0)**2;
815
816 e0d=((x-s0+self.Num_Bin)/w0)**2
817 m0=numpy.exp(-0.5*e0**(p0/2.));
818 m0u=numpy.exp(-0.5*e0u**(p0/2.));
819 m0d=numpy.exp(-0.5*e0d**(p0/2.))
820 JA=m0+m0u+m0d
821 JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d)
822
823 JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)
824
825 JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2
826 jack1=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,1./y_model])
827 return jack1.T
828
829 def y_jacobian2(self,x,state):
830 y_model=self.y_model2(x,state)
831 s0,w0,a0,p0,s1,w1,a1,p1,n=state
832 e0=((x-s0)/w0)**2;
833
834 e0u=((x-s0- self.Num_Bin )/w0)**2;
835
534
836 e0d=((x-s0+ self.Num_Bin )/w0)**2
535 return GauSPC
837 e1=((x-s1)/w1)**2;
838
839 e1u=((x-s1- self.Num_Bin )/w1)**2;
840
841 e1d=((x-s1+ self.Num_Bin )/w1)**2
842 m0=numpy.exp(-0.5*e0**(p0/2.));
843 m0u=numpy.exp(-0.5*e0u**(p0/2.));
844 m0d=numpy.exp(-0.5*e0d**(p0/2.))
845 m1=numpy.exp(-0.5*e1**(p1/2.));
846 m1u=numpy.exp(-0.5*e1u**(p1/2.));
847 m1d=numpy.exp(-0.5*e1d**(p1/2.))
848 JA=m0+m0u+m0d
849 JA1=m1+m1u+m1d
850 JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d)
851 JP1=(-1/4.)*a1*m1*e1**(p1/2.)*numpy.log(e1)+(-1/4.)*a1*m1u*e1u**(p1/2.)*numpy.log(e1u)+(-1/4.)*a1*m1d*e1d**(p1/2.)*numpy.log(e1d)
852
853 JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)
854
855 JS1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)
856
857 JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2
858
859 JW1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)**2+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)**2+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)**2
860 jack2=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,JS1/y_model,JW1/y_model,JA1/y_model,JP1/y_model,1./y_model])
861 return jack2.T
862
536
863 def y_model1(self,x,state):
537 def y_model1(self,x,state):
864 shift0,width0,amplitude0,power0,noise=state
538 shift0,width0,amplitude0,power0,noise=state
@@ -890,6 +564,7 class GaussianFit(Operation):
890 def misfit2(self,state,y_data,x,num_intg):
564 def misfit2(self,state,y_data,x,num_intg):
891 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
565 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
892
566
567
893
568
894 class PrecipitationProc(Operation):
569 class PrecipitationProc(Operation):
895
570
@@ -908,22 +583,31 class PrecipitationProc(Operation):
908 '''
583 '''
909
584
910
585
911 def run(self, dataOut, radar=None, Pt=None, Gt=None, Gr=None, Lambda=None, aL=None,
586 def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
912 tauW=None, ThetaT=None, ThetaR=None, Km = 0.93, Altitude=None):
587 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350):
913
588
914 self.spc = dataOut.data_pre[0].copy()
915 self.Num_Hei = self.spc.shape[2]
916 self.Num_Bin = self.spc.shape[1]
917 self.Num_Chn = self.spc.shape[0]
918
589
919 Velrange = dataOut.abscissaList
590 Velrange = dataOut.abscissaList
920
591
921 if radar == "MIRA35C" :
592 if radar == "MIRA35C" :
922
593
594 self.spc = dataOut.data_pre[0].copy()
595 self.Num_Hei = self.spc.shape[2]
596 self.Num_Bin = self.spc.shape[1]
597 self.Num_Chn = self.spc.shape[0]
923 Ze = self.dBZeMODE2(dataOut)
598 Ze = self.dBZeMODE2(dataOut)
924
599
925 else:
600 else:
926
601
602 self.spc = dataOut.GauSPC[1] #dataOut.data_pre[0].copy()
603 self.Num_Hei = self.spc.shape[2]
604 self.Num_Bin = self.spc.shape[1]
605 self.Num_Chn = self.spc.shape[0]
606 print '==================== SPC SHAPE',numpy.shape(self.spc)
607
608
609 ''' Se obtiene la constante del RADAR '''
610
927 self.Pt = Pt
611 self.Pt = Pt
928 self.Gt = Gt
612 self.Gt = Gt
929 self.Gr = Gr
613 self.Gr = Gr
@@ -933,48 +617,63 class PrecipitationProc(Operation):
933 self.ThetaT = ThetaT
617 self.ThetaT = ThetaT
934 self.ThetaR = ThetaR
618 self.ThetaR = ThetaR
935
619
936 RadarConstant = GetRadarConstant()
620 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
621 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
622 RadarConstant = Numerator / Denominator
623 print '***'
624 print '*** RadarConstant' , RadarConstant, '****'
625 print '***'
626 ''' ============================= '''
627
937 SPCmean = numpy.mean(self.spc,0)
628 SPCmean = numpy.mean(self.spc,0)
938 ETA = numpy.zeros(self.Num_Hei)
629 ETA = numpy.zeros(self.Num_Hei)
939 Pr = numpy.sum(SPCmean,0)
630 Pr = numpy.sum(SPCmean,0)
940
631
941 #for R in range(self.Num_Hei):
632 VelMeteoro = numpy.mean(SPCmean,axis=0)
942 # ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
633 print '==================== Vel SHAPE',VelMeteoro
943
634
944 D_range = numpy.zeros(self.Num_Hei)
635 D_range = numpy.zeros(self.Num_Hei)
636 SIGMA = numpy.zeros(self.Num_Hei)
637 N_dist = numpy.zeros(self.Num_Hei)
945 EqSec = numpy.zeros(self.Num_Hei)
638 EqSec = numpy.zeros(self.Num_Hei)
946 del_V = numpy.zeros(self.Num_Hei)
639 del_V = numpy.zeros(self.Num_Hei)
947
640
641
948 for R in range(self.Num_Hei):
642 for R in range(self.Num_Hei):
949 ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
643 ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
950
644
951 h = R + Altitude #Range from ground to radar pulse altitude
645 h = R + Altitude #Range from ground to radar pulse altitude
952 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
646 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
953
647
954 D_range[R] = numpy.log( (9.65 - (Velrange[R]/del_V[R])) / 10.3 ) / -0.6 #Range of Diameter of drops related to velocity
648 D_range[R] = numpy.log( (9.65 - (VelMeteoro[R]/del_V[R])) / 10.3 ) / -0.6 #Range of Diameter of drops related to velocity
955 SIGMA[R] = numpy.pi**5 / Lambda**4 * Km * D_range[R]**6 #Equivalent Section of drops (sigma)
649 SIGMA[R] = numpy.pi**5 / Lambda**4 * Km * D_range[R]**6 #Equivalent Section of drops (sigma)
956
650 #print '******* D_range ********', self.Num_Hei
651 #print D_range
957 N_dist[R] = ETA[R] / SIGMA[R]
652 N_dist[R] = ETA[R] / SIGMA[R]
653
654
958
655
959 Ze = (ETA * Lambda**4) / (numpy.pi * Km)
656 Ze = (ETA * Lambda**4) / (numpy.pi * Km)
960 Z = numpy.sum( N_dist * D_range**6 )
657 Z = numpy.sum( N_dist * D_range**6 )
961 RR = 6*10**-4*numpy.pi * numpy.sum( D_range**3 * N_dist * Velrange ) #Rainfall rate
658 #print 'Velrange',len(Velrange), 'D_range', len(D_range)
659 RR = 6*10**-4*numpy.pi * numpy.sum( D_range[R]**3 * N_dist[R] * VelMeteoro[R] ) #Rainfall rate
660
962
661
963
662
964 RR = (Ze/200)**(1/1.6)
663 RR2 = (Ze/200)**(1/1.6)
965 dBRR = 10*numpy.log10(RR)
664 dBRR = 10*numpy.log10(RR)
966
665
967 dBZe = 10*numpy.log10(Ze)
666 dBZe = 10*numpy.log10(Ze)
968 dataOut.data_output = Ze
667 dataOut.data_output = Ze
969 dataOut.data_param = numpy.ones([2,self.Num_Hei])
668 dataOut.data_param = numpy.ones([3,self.Num_Hei])
970 dataOut.channelList = [0,1]
669 dataOut.channelList = [0,1,2]
971 print 'channelList', dataOut.channelList
670 print 'channelList', dataOut.channelList
972 dataOut.data_param[0]=dBZe
671 dataOut.data_param[0]=dBZe
973 dataOut.data_param[1]=dBRR
672 dataOut.data_param[1]=dBRR
974 print 'RR SHAPE', dBRR.shape
673 print 'RR SHAPE', dBRR.shape
975 print 'Ze SHAPE', dBZe.shape
674 print 'Ze SHAPE', dBZe.shape
976 print 'dataOut.data_param SHAPE', dataOut.data_param.shape
675 print 'dataOut.data_param SHAPE', dataOut.data_param.shape
977
676
978
677
979 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
678 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
980
679
@@ -1001,26 +700,27 class PrecipitationProc(Operation):
1001
700
1002 return Ze
701 return Ze
1003
702
1004 def GetRadarConstant(self):
703 # def GetRadarConstant(self):
1005
704 #
1006 """
705 # """
1007 Constants:
706 # Constants:
1008
707 #
1009 Pt: Transmission Power dB 5kW
708 # Pt: Transmission Power dB 5kW 5000
1010 Gt: Transmission Gain dB 24.7 dB
709 # Gt: Transmission Gain dB 24.7 dB 295.1209
1011 Gr: Reception Gain dB 18.5 dB
710 # Gr: Reception Gain dB 18.5 dB 70.7945
1012 Lambda: Wavelenght m 0.6741 m
711 # Lambda: Wavelenght m 0.6741 m 0.6741
1013 aL: Attenuation loses dB
712 # aL: Attenuation loses dB 4dB 2.5118
1014 tauW: Width of transmission pulse s
713 # tauW: Width of transmission pulse s 4us 4e-6
1015 ThetaT: Transmission antenna bean angle rad 0.1656317 rad
714 # ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317
1016 ThetaR: Reception antenna beam angle rad 0.36774087 rad
715 # ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087
1017
716 #
1018 """
717 # """
1019 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
718 #
1020 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
719 # Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
1021 RadarConstant = Numerator / Denominator
720 # Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
1022
721 # RadarConstant = Numerator / Denominator
1023 return RadarConstant
722 #
723 # return RadarConstant
1024
724
1025
725
1026
726
@@ -1045,8 +745,10 class FullSpectralAnalysis(Operation):
1045 """
745 """
1046 def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None, SNRlimit=7):
746 def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None, SNRlimit=7):
1047
747
1048 spc = dataOut.data_pre[0].copy()
748 self.indice=int(numpy.random.rand()*1000)
1049 cspc = dataOut.data_pre[1].copy()
749
750 spc = dataOut.data_pre[0]
751 cspc = dataOut.data_pre[1]
1050
752
1051 nChannel = spc.shape[0]
753 nChannel = spc.shape[0]
1052 nProfiles = spc.shape[1]
754 nProfiles = spc.shape[1]
@@ -1060,10 +762,10 class FullSpectralAnalysis(Operation):
1060
762
1061 #print 'ChanDist', ChanDist
763 #print 'ChanDist', ChanDist
1062
764
1063 if dataOut.VelRange is not None:
765 #if dataOut.VelRange is not None:
1064 VelRange= dataOut.VelRange
766 AbbsisaRange= dataOut.spc_range#dataOut.VelRange
1065 else:
767 #else:
1066 VelRange= dataOut.abscissaList
768 # AbbsisaRange= dataOut.spc_range#dataOut.abscissaList
1067
769
1068 ySamples=numpy.ones([nChannel,nProfiles])
770 ySamples=numpy.ones([nChannel,nProfiles])
1069 phase=numpy.ones([nChannel,nProfiles])
771 phase=numpy.ones([nChannel,nProfiles])
@@ -1080,14 +782,16 class FullSpectralAnalysis(Operation):
1080 #print 'noise',noise
782 #print 'noise',noise
1081 #SNRdB = 10*numpy.log10(dataOut.data_SNR)
783 #SNRdB = 10*numpy.log10(dataOut.data_SNR)
1082
784
1083 FirstMoment = numpy.average(dataOut.data_param[:,1,:],0)
785 FirstMoment = dataOut.data_param[0,1,:]#numpy.average(dataOut.data_param[:,1,:],0)
786 SecondMoment = numpy.average(dataOut.data_param[:,2,:],0)
787
1084 #SNRdBMean = []
788 #SNRdBMean = []
1085
789
1086
790
1087 #for j in range(nHeights):
791 #for j in range(nHeights):
1088 # FirstMoment = numpy.append(FirstMoment,numpy.mean([dataOut.data_param[0,1,j],dataOut.data_param[1,1,j],dataOut.data_param[2,1,j]]))
792 # FirstMoment = numpy.append(FirstMoment,numpy.mean([dataOut.data_param[0,1,j],dataOut.data_param[1,1,j],dataOut.data_param[2,1,j]]))
1089 # SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]]))
793 # SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]]))
1090
794
1091 data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN
795 data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN
1092
796
1093 velocityX=[]
797 velocityX=[]
@@ -1097,21 +801,21 class FullSpectralAnalysis(Operation):
1097
801
1098 dbSNR = 10*numpy.log10(dataSNR)
802 dbSNR = 10*numpy.log10(dataSNR)
1099 dbSNR = numpy.average(dbSNR,0)
803 dbSNR = numpy.average(dbSNR,0)
804
1100 for Height in range(nHeights):
805 for Height in range(nHeights):
1101
806
1102 [Vzon,Vmer,Vver, GaussCenter, PhaseSlope]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR[Height], SNRlimit)
807 [Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC]= self.WindEstimation(dataOut.data_param, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR[Height], SNRlimit)
1103
1104 PhaseLine = numpy.append(PhaseLine, PhaseSlope)
808 PhaseLine = numpy.append(PhaseLine, PhaseSlope)
1105
809
1106 if abs(Vzon)<100. and abs(Vzon)> 0.:
810 if abs(Vzon)<100. and abs(Vzon)> 0.:
1107 velocityX=numpy.append(velocityX, Vzon)#Vmag
811 velocityX=numpy.append(velocityX, -Vzon)#Vmag
1108
812
1109 else:
813 else:
1110 #print 'Vzon',Vzon
814 #print 'Vzon',Vzon
1111 velocityX=numpy.append(velocityX, numpy.NaN)
815 velocityX=numpy.append(velocityX, numpy.NaN)
1112
816
1113 if abs(Vmer)<100. and abs(Vmer) > 0.:
817 if abs(Vmer)<100. and abs(Vmer) > 0.:
1114 velocityY=numpy.append(velocityY, Vmer)#Vang
818 velocityY=numpy.append(velocityY, -Vmer)#Vang
1115
819
1116 else:
820 else:
1117 #print 'Vmer',Vmer
821 #print 'Vmer',Vmer
@@ -1129,24 +833,27 class FullSpectralAnalysis(Operation):
1129
833
1130
834
1131
835
1132 data_output[0]=numpy.array(velocityX)
836 data_output[0] = numpy.array(velocityX) #self.moving_average(numpy.array(velocityX) , N=1)
1133 data_output[1]=numpy.array(velocityY)
837 data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1)
1134 data_output[2]=-velocityV#FirstMoment
838 data_output[2] = -velocityV#FirstMoment
1135
839
1136 print ' '
840 print ' '
1137 #print 'FirstMoment'
841 #print 'FirstMoment'
1138 #print FirstMoment
842 #print FirstMoment
843 print 'velocityX',numpy.shape(data_output[0])
1139 print 'velocityX',data_output[0]
844 print 'velocityX',data_output[0]
1140 print ' '
845 print ' '
846 print 'velocityY',numpy.shape(data_output[1])
1141 print 'velocityY',data_output[1]
847 print 'velocityY',data_output[1]
848 print 'velocityV',data_output[2]
1142 print 'PhaseLine',PhaseLine
849 print 'PhaseLine',PhaseLine
1143 #print numpy.array(velocityY)
850 #print numpy.array(velocityY)
1144 print ' '
1145 #print 'SNR'
851 #print 'SNR'
1146 #print 10*numpy.log10(dataOut.data_SNR)
852 #print 10*numpy.log10(dataOut.data_SNR)
1147 #print numpy.shape(10*numpy.log10(dataOut.data_SNR))
853 #print numpy.shape(10*numpy.log10(dataOut.data_SNR))
1148 print ' '
854 print ' '
1149
855
856 xFrec=AbbsisaRange[0][0:spc.shape[1]]
1150
857
1151 dataOut.data_output=data_output
858 dataOut.data_output=data_output
1152
859
@@ -1156,15 +863,26 class FullSpectralAnalysis(Operation):
1156 def moving_average(self,x, N=2):
863 def moving_average(self,x, N=2):
1157 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
864 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
1158
865
1159 def gaus(self,xSamples,a,x0,sigma):
866 def gaus(self,xSamples,Amp,Mu,Sigma):
1160 return a*numpy.exp(-(xSamples-x0)**2/(2*sigma**2))
867 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
1161
868
1162 def Find(self,x,value):
1163 for index in range(len(x)):
1164 if x[index]==value:
1165 return index
1166
869
1167 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR, SNRlimit):
870
871 def Moments(self, ySamples, xSamples):
872 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
873 yNorm = ySamples / Pot
874
875 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
876 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
877 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
878
879 return numpy.array([Pot, Vr, Desv])
880
881 def WindEstimation(self, data_param, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
882
883 print ' '
884 print '######################## Height',Height, (1000 + 75*Height), '##############################'
885 print ' '
1168
886
1169 ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
887 ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
1170 phase=numpy.ones([spc.shape[0],spc.shape[1]])
888 phase=numpy.ones([spc.shape[0],spc.shape[1]])
@@ -1172,7 +890,14 class FullSpectralAnalysis(Operation):
1172 coherence=numpy.ones([spc.shape[0],spc.shape[1]])
890 coherence=numpy.ones([spc.shape[0],spc.shape[1]])
1173 PhaseSlope=numpy.zeros(spc.shape[0])
891 PhaseSlope=numpy.zeros(spc.shape[0])
1174 PhaseInter=numpy.ones(spc.shape[0])
892 PhaseInter=numpy.ones(spc.shape[0])
1175 xFrec=VelRange
893 xFrec=AbbsisaRange[0][0:spc.shape[1]]
894 xVel =AbbsisaRange[2][0:spc.shape[1]]
895 Vv=numpy.empty(spc.shape[2])*0
896 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) #spc[0]-noise[0]#
897
898 SPCmoments = self.Moments(SPCav[:,Height], xVel )
899 CSPCmoments = []
900 cspcNoise = numpy.empty(3)
1176
901
1177 '''Getting Eij and Nij'''
902 '''Getting Eij and Nij'''
1178
903
@@ -1191,150 +916,243 class FullSpectralAnalysis(Operation):
1191 for i in range(spc.shape[0]):
916 for i in range(spc.shape[0]):
1192
917
1193 '''****** Line of Data SPC ******'''
918 '''****** Line of Data SPC ******'''
1194 zline=z[i,:,Height]
919 zline=z[i,:,Height].copy() - noise[i] # Se resta ruido
1195
920
1196 '''****** SPC is normalized ******'''
921 '''****** SPC is normalized ******'''
1197 FactNorm= (zline.copy()-noise[i]) / numpy.sum(zline.copy())
922 SmoothSPC =self.moving_average(zline.copy(),N=1) # Se suaviza el ruido
1198 FactNorm= FactNorm/numpy.sum(FactNorm)
923 FactNorm = SmoothSPC/numpy.nansum(SmoothSPC) # SPC Normalizado y suavizado
1199
1200 SmoothSPC=self.moving_average(FactNorm,N=3)
1201
924
1202 xSamples = ar(range(len(SmoothSPC)))
925 xSamples = xFrec # Se toma el rango de frecuncias
1203 ySamples[i] = SmoothSPC
926 ySamples[i] = FactNorm # Se toman los valores de SPC normalizado
1204
1205 #dbSNR=10*numpy.log10(dataSNR)
1206 print ' '
1207 print ' '
1208 print ' '
1209
1210 #print 'dataSNR', dbSNR.shape, dbSNR[0,40:120]
1211 #print 'SmoothSPC', SmoothSPC.shape, SmoothSPC[0:20]
1212 #print 'noise',noise
1213 #print 'zline',zline.shape, zline[0:20]
1214 #print 'FactNorm',FactNorm.shape, FactNorm[0:20]
1215 #print 'FactNorm suma', numpy.sum(FactNorm)
1216
927
1217 for i in range(spc.shape[0]):
928 for i in range(spc.shape[0]):
1218
929
1219 '''****** Line of Data CSPC ******'''
930 '''****** Line of Data CSPC ******'''
1220 cspcLine=cspc[i,:,Height].copy()
931 cspcLine = ( cspc[i,:,Height].copy())# - noise[i] ) # no! Se resta el ruido
932 SmoothCSPC =self.moving_average(cspcLine,N=1) # Se suaviza el ruido
933 cspcNorm = SmoothCSPC/numpy.nansum(SmoothCSPC) # CSPC normalizado y suavizado
1221
934
1222 '''****** CSPC is normalized ******'''
935 '''****** CSPC is normalized with respect to Briggs and Vincent ******'''
1223 chan_index0 = pairsList[i][0]
936 chan_index0 = pairsList[i][0]
1224 chan_index1 = pairsList[i][1]
937 chan_index1 = pairsList[i][1]
1225 CSPCFactor= abs(numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])) #
1226
938
1227 CSPCNorm = (cspcLine.copy() -noise[i]) / numpy.sqrt(CSPCFactor)
939 CSPCFactor= numpy.abs(numpy.nansum(ySamples[chan_index0]))**2 * numpy.abs(numpy.nansum(ySamples[chan_index1]))**2
940 CSPCNorm = cspcNorm / numpy.sqrt(CSPCFactor)
1228
941
1229 CSPCSamples[i] = CSPCNorm
942 CSPCSamples[i] = CSPCNorm
943
1230 coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
944 coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
1231
945
1232 coherence[i]= self.moving_average(coherence[i],N=2)
946 #coherence[i]= self.moving_average(coherence[i],N=1)
1233
947
1234 phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
948 phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
1235
949
1236 #print 'cspcLine', cspcLine.shape, cspcLine[0:20]
950 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPCSamples[0]), xSamples),
1237 #print 'CSPCFactor', CSPCFactor#, CSPCFactor[0:20]
951 self.Moments(numpy.abs(CSPCSamples[1]), xSamples),
1238 #print numpy.sum(ySamples[chan_index0]), numpy.sum(ySamples[chan_index1]), -noise[i]
952 self.Moments(numpy.abs(CSPCSamples[2]), xSamples)])
1239 #print 'CSPCNorm', CSPCNorm.shape, CSPCNorm[0:20]
1240 #print 'CSPCNorm suma', numpy.sum(CSPCNorm)
1241 #print 'CSPCSamples', CSPCSamples.shape, CSPCSamples[0,0:20]
1242
953
1243 '''****** Getting fij width ******'''
954 #print '##### SUMA de SPC #####', len(ySamples)
955 #print numpy.sum(ySamples[0])
956 #print '##### SUMA de CSPC #####', len(coherence)
957 #print numpy.sum(numpy.abs(CSPCNorm))
958 #print numpy.sum(coherence[0])
959 print 'len',len(xSamples)
960 print 'CSPCmoments', numpy.shape(CSPCmoments)
961 print CSPCmoments
962 print '#######################'
963
964 popt=[1e-10,1e-10,1e-10]
965 popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
966 FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
967
968 CSPCMask01 = numpy.abs(CSPCSamples[0])
969 CSPCMask02 = numpy.abs(CSPCSamples[1])
970 CSPCMask12 = numpy.abs(CSPCSamples[2])
1244
971
1245 yMean=[]
972 mask01 = ~numpy.isnan(CSPCMask01)
1246 yMean2=[]
973 mask02 = ~numpy.isnan(CSPCMask02)
974 mask12 = ~numpy.isnan(CSPCMask12)
1247
975
1248 for j in range(len(ySamples[1])):
976 #mask = ~numpy.isnan(CSPCMask01)
1249 yMean=numpy.append(yMean,numpy.mean([ySamples[0,j],ySamples[1,j],ySamples[2,j]]))
977 CSPCMask01 = CSPCMask01[mask01]
978 CSPCMask02 = CSPCMask02[mask02]
979 CSPCMask12 = CSPCMask12[mask12]
980 #CSPCMask01 = numpy.ma.masked_invalid(CSPCMask01)
1250
981
1251 '''******* Getting fitting Gaussian ******'''
1252 meanGauss=sum(xSamples*yMean) / len(xSamples)
1253 sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples)
1254
982
1255 #print '****************************'
1256 #print 'len(xSamples): ',len(xSamples)
1257 #print 'yMean: ', yMean.shape, yMean[0:20]
1258 #print 'ySamples', ySamples.shape, ySamples[0,0:20]
1259 #print 'xSamples: ',xSamples.shape, xSamples[0:20]
1260
983
1261 #print 'meanGauss',meanGauss
1262 #print 'sigma',sigma
1263
984
1264 #if (abs(meanGauss/sigma**2) > 0.0001) : #0.000000001):
985
1265 if dbSNR > SNRlimit and abs(meanGauss/sigma**2) > 0.0001:
986 '''***Fit Gauss CSPC01***'''
1266 try:
987 if dbSNR > SNRlimit:
1267 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
988 try:
989 popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
990 popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
991 popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2])
992 FitGauss01 = self.gaus(xSamples,*popt01)
993 FitGauss02 = self.gaus(xSamples,*popt02)
994 FitGauss12 = self.gaus(xSamples,*popt12)
995 except:
996 FitGauss01=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[0]))
997 FitGauss02=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[1]))
998 FitGauss12=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[2]))
999
1000
1001 CSPCopt = numpy.vstack([popt01,popt02,popt12])
1002
1003 '''****** Getting fij width ******'''
1004
1005 yMean = numpy.average(ySamples, axis=0) # ySamples[0]
1006
1007 '''******* Getting fitting Gaussian *******'''
1008 meanGauss = sum(xSamples*yMean) / len(xSamples) # Mu, velocidad radial (frecuencia)
1009 sigma2 = sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) # Varianza, Ancho espectral (frecuencia)
1010
1011 yMoments = self.Moments(yMean, xSamples)
1012
1013 if dbSNR > SNRlimit: # and abs(meanGauss/sigma2) > 0.00001:
1014 try:
1015 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments)
1016 FitGauss=self.gaus(xSamples,*popt)
1268
1017
1269 if numpy.amax(popt)>numpy.amax(yMean)*0.3:
1270 FitGauss=self.gaus(xSamples,*popt)
1271
1272 else:
1273 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1274 print 'Verificador: Dentro', Height
1275 except :#RuntimeError:
1018 except :#RuntimeError:
1276 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1019 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1277
1020
1278
1021
1279 else:
1022 else:
1280 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1023 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1281
1024
1282 Maximun=numpy.amax(yMean)
1283 eMinus1=Maximun*numpy.exp(-1)#*0.8
1284
1025
1285 HWpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1)))
1286 HalfWidth= xFrec[HWpos]
1287 GCpos=self.Find(FitGauss, numpy.amax(FitGauss))
1288 Vpos=self.Find(FactNorm, numpy.amax(FactNorm))
1289
1290 #Vpos=FirstMoment[]
1291
1026
1292 '''****** Getting Fij ******'''
1027 '''****** Getting Fij ******'''
1028 Fijcspc = CSPCopt[:,2]/2*3
1029
1030 #GauWidth = (popt[2]/2)*3
1031 GaussCenter = popt[1] #xFrec[GCpos]
1032 #Punto en Eje X de la Gaussiana donde se encuentra el centro
1033 ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
1034 PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0]
1035
1036 #Punto e^-1 hubicado en la Gaussiana
1037 PeMinus1 = numpy.max(FitGauss)* numpy.exp(-1)
1038 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss"
1039 PointFij = numpy.where(FitGauss==FijClosest)[0][0]
1293
1040
1294 GaussCenter=xFrec[GCpos]
1041 if xSamples[PointFij] > xSamples[PointGauCenter]:
1295 if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
1042 Fij = xSamples[PointFij] - xSamples[PointGauCenter]
1296 Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
1043
1297 else:
1044 else:
1298 Fij=abs(GaussCenter-HalfWidth)+0.0000001
1045 Fij = xSamples[PointGauCenter] - xSamples[PointFij]
1046
1047 print 'CSPCopt'
1048 print CSPCopt
1049 print 'popt'
1050 print popt
1051 print '#######################################'
1052 #print 'dataOut.data_param', numpy.shape(data_param)
1053 #print 'dataOut.data_param0', data_param[0,0,Height]
1054 #print 'dataOut.data_param1', data_param[0,1,Height]
1055 #print 'dataOut.data_param2', data_param[0,2,Height]
1056
1057
1058 print 'yMoments', yMoments
1059 print 'Moments', SPCmoments
1060 print 'Fij2 Moment', Fij
1061 #print 'Fij', Fij, 'popt[2]/2',popt[2]/2
1062 print 'Fijcspc',Fijcspc
1063 print '#######################################'
1064
1065
1066 # Range = numpy.empty([3,2])
1067 # GaussCenter = numpy.empty(3)
1068 # FrecRange, VelRange = [[],[],[]],[[],[],[]]
1069 # FrecRange01, FrecRange02, FrecRange12 = numpy.empty(3), numpy.empty(3), numpy.empty(3)
1070 # VelRange01, VelRange02, VelRange12 = numpy.empty(3), numpy.empty(3), numpy.empty(3)
1071 #
1072 # GauWidth = numpy.empty(3)
1073 # ClosRangeMin, ClosRangeMax = numpy.empty(3), numpy.empty(3)
1074 # PointRangeMin, PointRangeMax = numpy.empty(3), numpy.empty(3)
1075 # '''****** Taking frequency ranges from CSPCs ******'''
1076 # for i in range(spc.shape[0]):
1077 #
1078 # GaussCenter[i] = CSPCopt[i][1] #Primer momento 01
1079 # GauWidth[i] = CSPCopt[i][2]*2/2 #Ancho de banda de Gau01
1080 #
1081 # Range[i][0] = GaussCenter[i] - GauWidth[i]
1082 # Range[i][1] = GaussCenter[i] + GauWidth[i]
1083 # #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max)
1084 # ClosRangeMin[i] = xSamples[numpy.abs(xSamples-Range[i][0]).argmin()]
1085 # ClosRangeMax[i] = xSamples[numpy.abs(xSamples-Range[i][1]).argmin()]
1086 #
1087 # PointRangeMin[i] = numpy.where(xSamples==ClosRangeMin[i])[0][0]
1088 # PointRangeMax[i] = numpy.where(xSamples==ClosRangeMax[i])[0][0]
1089 #
1090 # Range[i]=numpy.array([ PointRangeMin[i], PointRangeMax[i] ])
1091 #
1092 # FrecRange[i] = xFrec[ Range[i][0] : Range[i][1] ]
1093 # VelRange[i] = xVel[ Range[i][0] : Range[i][1] ]
1094
1299
1095
1300 '''****** Getting Frecuency range of significant data ******'''
1301
1096
1302 Rangpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10)))
1097 '''****** Taking frequency ranges from SPCs ******'''
1303
1098
1304 if Rangpos<GCpos:
1305 Range=numpy.array([Rangpos,2*GCpos-Rangpos])
1306 elif Rangpos< ( len(xFrec)- len(xFrec)*0.1):
1307 Range=numpy.array([2*GCpos-Rangpos,Rangpos])
1308 else:
1309 Range = numpy.array([0,0])
1310
1099
1311 #print ' '
1100 #GaussCenter = popt[1] #Primer momento 01
1312 #print 'GCpos',GCpos, ( len(xFrec)- len(xFrec)*0.1)
1101 GauWidth = popt[2] *3/2 #Ancho de banda de Gau01
1313 #print 'Rangpos',Rangpos
1102 Range = numpy.empty(2)
1314 print 'RANGE: ', Range
1103 Range[0] = GaussCenter - GauWidth
1315 FrecRange=xFrec[Range[0]:Range[1]]
1104 Range[1] = GaussCenter + GauWidth
1105 #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max)
1106 ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()]
1107 ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()]
1108
1109 PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
1110 PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
1111
1112 Range=numpy.array([ PointRangeMin, PointRangeMax ])
1113
1114 FrecRange = xFrec[ Range[0] : Range[1] ]
1115 VelRange = xVel[ Range[0] : Range[1] ]
1116
1117
1118 #print 'RANGE: ', Range
1119 #print 'FrecRange', numpy.shape(FrecRange)#,FrecRange
1120 #print 'len: ', len(FrecRange)
1316
1121
1317 '''****** Getting SCPC Slope ******'''
1122 '''****** Getting SCPC Slope ******'''
1318
1123
1319 for i in range(spc.shape[0]):
1124 for i in range(spc.shape[0]):
1320
1125
1321 if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.5:
1126 if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.6:
1322 PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
1127 PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
1323
1128
1324 print 'FrecRange', len(FrecRange) , FrecRange
1129 #print 'Ancho espectral Frecuencias', FrecRange[-1]-FrecRange[0], 'Hz'
1325 print 'PhaseRange', len(PhaseRange), PhaseRange
1130 #print 'Ancho espectral Velocidades', VelRange[-1]-VelRange[0], 'm/s'
1326 print ' '
1131 #print 'FrecRange', len(FrecRange) , FrecRange
1132 #print 'VelRange', len(VelRange) , VelRange
1133 #print 'PhaseRange', numpy.shape(PhaseRange), PhaseRange
1134 #print ' '
1135
1136 '''***********************VelRange******************'''
1137
1138 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
1139
1327 if len(FrecRange) == len(PhaseRange):
1140 if len(FrecRange) == len(PhaseRange):
1328 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
1141 try:
1329 PhaseSlope[i]=slope
1142 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange[mask], PhaseRange[mask])
1330 PhaseInter[i]=intercept
1143 PhaseSlope[i]=slope
1144 PhaseInter[i]=intercept
1145 except:
1146 PhaseSlope[i]=0
1147 PhaseInter[i]=0
1331 else:
1148 else:
1332 PhaseSlope[i]=0
1149 PhaseSlope[i]=0
1333 PhaseInter[i]=0
1150 PhaseInter[i]=0
1334 else:
1151 else:
1335 PhaseSlope[i]=0
1152 PhaseSlope[i]=0
1336 PhaseInter[i]=0
1153 PhaseInter[i]=0
1337
1154
1155
1338 '''Getting constant C'''
1156 '''Getting constant C'''
1339 cC=(Fij*numpy.pi)**2
1157 cC=(Fij*numpy.pi)**2
1340
1158
@@ -1346,9 +1164,9 class FullSpectralAnalysis(Operation):
1346 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1164 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1347
1165
1348 '''****** Getting constants A, B and H ******'''
1166 '''****** Getting constants A, B and H ******'''
1349 W01=numpy.amax(coherence[0])
1167 W01=numpy.nanmax( FitGauss01 ) #numpy.abs(CSPCSamples[0]))
1350 W02=numpy.amax(coherence[1])
1168 W02=numpy.nanmax( FitGauss02 ) #numpy.abs(CSPCSamples[1]))
1351 W12=numpy.amax(coherence[2])
1169 W12=numpy.nanmax( FitGauss12 ) #numpy.abs(CSPCSamples[2]))
1352
1170
1353 WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1171 WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1354 WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1172 WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
@@ -1363,16 +1181,58 class FullSpectralAnalysis(Operation):
1363
1181
1364 VxVyResults=numpy.array([-cF,-cG])
1182 VxVyResults=numpy.array([-cF,-cG])
1365 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1183 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1366
1184 #print 'MijResults, cC, PhaseSlope', MijResults, cC, PhaseSlope
1185 #print 'W01,02,12', W01, W02, W12
1186 #print 'WijResult0,1,2',WijResult0, WijResult1, WijResult2, 'Results', WijResults
1187 #print 'cA,cB,cH, cF, cG', cA, cB, cH, cF, cG
1188 #print 'VxVy', VxVyResults
1189 #print '###########################****************************************'
1367 Vzon = Vy
1190 Vzon = Vy
1368 Vmer = Vx
1191 Vmer = Vx
1369 Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1192 Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1370 Vang=numpy.arctan2(Vmer,Vzon)
1193 Vang=numpy.arctan2(Vmer,Vzon)
1371 Vver=xFrec[Vpos]
1194 Vver=SPCmoments[1]
1372 print 'Height',Height
1195 FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12])
1373 print 'vzon y vmer', Vzon, Vmer
1196
1374 return Vzon, Vmer, Vver, GaussCenter, PhaseSlope
1197
1198 ''' Ploteo por altura '''
1199 if Height == 28:
1200 for i in range(3):
1201 #print 'FASE', numpy.shape(phase), y[25]
1202 #print numpy.shape(coherence)
1203 fig = plt.figure(10+self.indice)
1204 #plt.plot( x[0:256],coherence[:,25] )
1205 #cohAv = numpy.average(coherence[i],1)
1206 Pendiente = FrecRange * PhaseSlope[i]
1207 plt.plot( FrecRange, Pendiente)
1208 plt.plot( xFrec,phase[i])
1209
1210 CSPCmean = numpy.mean(numpy.abs(CSPCSamples),0)
1211 #plt.plot(xFrec, FitGauss01)
1212 #plt.plot(xFrec, CSPCmean)
1213 #plt.plot(xFrec, numpy.abs(CSPCSamples[0]))
1214 #plt.plot(xFrec, FitGauss)
1215 #plt.plot(xFrec, yMean)
1216 #plt.plot(xFrec, numpy.abs(coherence[0]))
1217
1218 #plt.axis([-12, 12, 15, 50])
1219 #plt.title("%s" %( '%s %s, Channel %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S") , i)))
1220 plt.ylabel('Desfase [rad]')
1221 #plt.ylabel('CSPC normalizado')
1222 plt.xlabel('Frec range [Hz]')
1223
1224 #fig.savefig('/home/erick/Documents/Pics/to{}.png'.format(self.indice))
1225
1226 plt.show()
1227 self.indice=self.indice+1
1375
1228
1229
1230
1231
1232
1233 print 'vzon y vmer', Vzon, Vmer
1234 return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC
1235
1376 class SpectralMoments(Operation):
1236 class SpectralMoments(Operation):
1377
1237
1378 '''
1238 '''
@@ -1433,7 +1293,7 class SpectralMoments(Operation):
1433 if (aliasing == None): aliasing = 0
1293 if (aliasing == None): aliasing = 0
1434 if (oldfd == None): oldfd = 0
1294 if (oldfd == None): oldfd = 0
1435 if (wwauto == None): wwauto = 0
1295 if (wwauto == None): wwauto = 0
1436
1296
1437 if (n0 < 1.e-20): n0 = 1.e-20
1297 if (n0 < 1.e-20): n0 = 1.e-20
1438
1298
1439 freq = oldfreq
1299 freq = oldfreq
@@ -1441,6 +1301,8 class SpectralMoments(Operation):
1441 vec_fd = numpy.zeros(oldspec.shape[1])
1301 vec_fd = numpy.zeros(oldspec.shape[1])
1442 vec_w = numpy.zeros(oldspec.shape[1])
1302 vec_w = numpy.zeros(oldspec.shape[1])
1443 vec_snr = numpy.zeros(oldspec.shape[1])
1303 vec_snr = numpy.zeros(oldspec.shape[1])
1304
1305 oldspec = numpy.ma.masked_invalid(oldspec)
1444
1306
1445 for ind in range(oldspec.shape[1]):
1307 for ind in range(oldspec.shape[1]):
1446
1308
@@ -1475,11 +1337,11 class SpectralMoments(Operation):
1475 if (ss1 > m): ss1 = m
1337 if (ss1 > m): ss1 = m
1476
1338
1477 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
1339 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
1478 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
1340 power = ( (spec2[valid] - n0) * fwindow[valid] ).sum()
1479 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
1341 fd = ( (spec2[valid]- n0) * freq[valid] * fwindow[valid] ).sum() / power
1480 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
1342 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
1481 snr = (spec2.mean()-n0)/n0
1343 snr = (spec2.mean()-n0)/n0
1482
1344
1483 if (snr < 1.e-20) :
1345 if (snr < 1.e-20) :
1484 snr = 1.e-20
1346 snr = 1.e-20
1485
1347
@@ -1487,7 +1349,7 class SpectralMoments(Operation):
1487 vec_fd[ind] = fd
1349 vec_fd[ind] = fd
1488 vec_w[ind] = w
1350 vec_w[ind] = w
1489 vec_snr[ind] = snr
1351 vec_snr[ind] = snr
1490
1352
1491 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
1353 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
1492 return moments
1354 return moments
1493
1355
@@ -1744,7 +1606,7 class SpectralFitting(Operation):
1744
1606
1745 fm = self.dataOut.library.modelFunction(p, constants)
1607 fm = self.dataOut.library.modelFunction(p, constants)
1746 fmp=numpy.dot(LT,fm)
1608 fmp=numpy.dot(LT,fm)
1747
1609
1748 return dp-fmp
1610 return dp-fmp
1749
1611
1750 def __getSNR(self, z, noise):
1612 def __getSNR(self, z, noise):
@@ -1778,8 +1640,8 class WindProfiler(Operation):
1778
1640
1779 n = None
1641 n = None
1780
1642
1781 def __init__(self):
1643 def __init__(self, **kwargs):
1782 Operation.__init__(self)
1644 Operation.__init__(self, **kwargs)
1783
1645
1784 def __calculateCosDir(self, elev, azim):
1646 def __calculateCosDir(self, elev, azim):
1785 zen = (90 - elev)*numpy.pi/180
1647 zen = (90 - elev)*numpy.pi/180
@@ -2293,7 +2155,7 class WindProfiler(Operation):
2293
2155
2294 return data_output
2156 return data_output
2295
2157
2296 def run(self, dataOut, technique, **kwargs):
2158 def run(self, dataOut, technique, positionY, positionX, azimuth, **kwargs):
2297
2159
2298 param = dataOut.data_param
2160 param = dataOut.data_param
2299 if dataOut.abscissaList != None:
2161 if dataOut.abscissaList != None:
Show file before | Show file after | Hide comments
@@ -4,6 +4,8 from jroproc_base import ProcessingUnit, Operation
4 from schainpy.model.data.jrodata import Spectra
4 from schainpy.model.data.jrodata import Spectra
5 from schainpy.model.data.jrodata import hildebrand_sekhon
5 from schainpy.model.data.jrodata import hildebrand_sekhon
6
6
7 import matplotlib.pyplot as plt
8
7 class SpectraProc(ProcessingUnit):
9 class SpectraProc(ProcessingUnit):
8
10
9 def __init__(self, **kwargs):
11 def __init__(self, **kwargs):
@@ -275,7 +277,46 class SpectraProc(ProcessingUnit):
275 self.__selectPairsByChannel(self.dataOut.channelList)
277 self.__selectPairsByChannel(self.dataOut.channelList)
276
278
277 return 1
279 return 1
280
281
282 def selectFFTs(self, minFFT, maxFFT ):
283 """
284 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
285 minFFT<= FFT <= maxFFT
286 """
287
288 if (minFFT > maxFFT):
289 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT)
290
291 if (minFFT < self.dataOut.getFreqRange()[0]):
292 minFFT = self.dataOut.getFreqRange()[0]
293
294 if (maxFFT > self.dataOut.getFreqRange()[-1]):
295 maxFFT = self.dataOut.getFreqRange()[-1]
296
297 minIndex = 0
298 maxIndex = 0
299 FFTs = self.dataOut.getFreqRange()
300
301 inda = numpy.where(FFTs >= minFFT)
302 indb = numpy.where(FFTs <= maxFFT)
278
303
304 try:
305 minIndex = inda[0][0]
306 except:
307 minIndex = 0
308
309 try:
310 maxIndex = indb[0][-1]
311 except:
312 maxIndex = len(FFTs)
313
314 self.selectFFTsByIndex(minIndex, maxIndex)
315
316 return 1
317
318
319
279 def selectHeights(self, minHei, maxHei):
320 def selectHeights(self, minHei, maxHei):
280 """
321 """
281 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
322 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
@@ -292,6 +333,7 class SpectraProc(ProcessingUnit):
292 1 si el metodo se ejecuto con exito caso contrario devuelve 0
333 1 si el metodo se ejecuto con exito caso contrario devuelve 0
293 """
334 """
294
335
336
295 if (minHei > maxHei):
337 if (minHei > maxHei):
296 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
338 raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
297
339
@@ -319,6 +361,7 class SpectraProc(ProcessingUnit):
319 maxIndex = len(heights)
361 maxIndex = len(heights)
320
362
321 self.selectHeightsByIndex(minIndex, maxIndex)
363 self.selectHeightsByIndex(minIndex, maxIndex)
364
322
365
323 return 1
366 return 1
324
367
@@ -361,6 +404,41 class SpectraProc(ProcessingUnit):
361
404
362 return 1
405 return 1
363
406
407 def selectFFTsByIndex(self, minIndex, maxIndex):
408 """
409
410 """
411
412 if (minIndex < 0) or (minIndex > maxIndex):
413 raise ValueError, "Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex)
414
415 if (maxIndex >= self.dataOut.nProfiles):
416 maxIndex = self.dataOut.nProfiles-1
417
418 #Spectra
419 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
420
421 data_cspc = None
422 if self.dataOut.data_cspc is not None:
423 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
424
425 data_dc = None
426 if self.dataOut.data_dc is not None:
427 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
428
429 self.dataOut.data_spc = data_spc
430 self.dataOut.data_cspc = data_cspc
431 self.dataOut.data_dc = data_dc
432
433 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
434 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
435 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
436
437 #self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
438
439 return 1
440
441
364
442
365 def selectHeightsByIndex(self, minIndex, maxIndex):
443 def selectHeightsByIndex(self, minIndex, maxIndex):
366 """
444 """
@@ -405,7 +483,8 class SpectraProc(ProcessingUnit):
405 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
483 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
406
484
407 return 1
485 return 1
408
486
487
409 def removeDC(self, mode = 2):
488 def removeDC(self, mode = 2):
410 jspectra = self.dataOut.data_spc
489 jspectra = self.dataOut.data_spc
411 jcspectra = self.dataOut.data_cspc
490 jcspectra = self.dataOut.data_cspc
@@ -463,6 +542,86 class SpectraProc(ProcessingUnit):
463
542
464 return 1
543 return 1
465
544
545 def removeInterference2(self):
546
547 cspc = self.dataOut.data_cspc
548 spc = self.dataOut.data_spc
549 print numpy.shape(spc)
550 Heights = numpy.arange(cspc.shape[2])
551 realCspc = numpy.abs(cspc)
552
553 for i in range(cspc.shape[0]):
554 LinePower= numpy.sum(realCspc[i], axis=0)
555 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
556 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
557 #print numpy.shape(realCspc)
558 #print '',SelectedHeights, '', numpy.shape(realCspc[i,:,SelectedHeights])
559 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
560 print SelectedHeights
561 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
562 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
563
564
565 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
566 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
567 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
568 cspc[i,InterferenceRange,:] = numpy.NaN
569
570 print '########################################################################################'
571 print 'Len interference sum',len(InterferenceSum)
572 print 'InterferenceThresholdMin', InterferenceThresholdMin, 'InterferenceThresholdMax', InterferenceThresholdMax
573 print 'InterferenceRange',InterferenceRange
574 print '########################################################################################'
575
576 ''' Ploteo '''
577
578 #for i in range(3):
579 #print 'FASE', numpy.shape(phase), y[25]
580 #print numpy.shape(coherence)
581 #fig = plt.figure(10+ int(numpy.random.rand()*100))
582 #plt.plot( x[0:256],coherence[:,25] )
583 #cohAv = numpy.average(coherence[i],1)
584 #Pendiente = FrecRange * PhaseSlope[i]
585 #plt.plot( InterferenceSum)
586 #plt.plot( numpy.sort(InterferenceSum))
587 #plt.plot( LinePower )
588 #plt.plot( xFrec,phase[i])
589
590 #CSPCmean = numpy.mean(numpy.abs(CSPCSamples),0)
591 #plt.plot(xFrec, FitGauss01)
592 #plt.plot(xFrec, CSPCmean)
593 #plt.plot(xFrec, numpy.abs(CSPCSamples[0]))
594 #plt.plot(xFrec, FitGauss)
595 #plt.plot(xFrec, yMean)
596 #plt.plot(xFrec, numpy.abs(coherence[0]))
597
598 #plt.axis([-12, 12, 15, 50])
599 #plt.title("%s" %( '%s %s, Channel %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S") , i)))
600
601
602 #fig.savefig('/home/erick/Documents/Pics/nom{}.png'.format(int(numpy.random.rand()*100)))
603
604 #plt.show()
605 #self.indice=self.indice+1
606 #raise
607
608
609 self.dataOut.data_cspc = cspc
610
611 # for i in range(spc.shape[0]):
612 # LinePower= numpy.sum(spc[i], axis=0)
613 # Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
614 # SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
615 # #print numpy.shape(realCspc)
616 # #print '',SelectedHeights, '', numpy.shape(realCspc[i,:,SelectedHeights])
617 # InterferenceSum = numpy.sum( spc[i,:,SelectedHeights], axis=0 )
618 # InterferenceThreshold = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
619 # InterferenceRange = numpy.where( InterferenceSum > InterferenceThreshold )
620 # if len(InterferenceRange)<int(spc.shape[1]*0.03):
621 # spc[i,InterferenceRange,:] = numpy.NaN
622
623 #self.dataOut.data_spc = spc
624
466 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
625 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
467
626
468 jspectra = self.dataOut.data_spc
627 jspectra = self.dataOut.data_spc
General Comments 0
You need to be logged in to leave comments. Login now