| @@ -706,7 +706,7 class ProcUnitConf(): | |||
|
|
706 | 706 | #print "\tRunning the '%s' operation with %s" %(opConfObj.name, opConfObj.id) |
|
|
707 | 707 | sts = self.procUnitObj.call(opType = opConfObj.type, |
|
|
708 | 708 | opName = opConfObj.name, |
|
|
709 |
opId = opConfObj.id |
|
|
|
709 | opId = opConfObj.id | |
|
|
710 | 710 | ) |
|
|
711 | 711 | |
|
|
712 | 712 | # total_time = time.time() - ini |
| @@ -82,7 +82,7 class FitGauPlot(Figure): | |||
|
|
82 | 82 | save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1, |
|
|
83 | 83 | server=None, folder=None, username=None, password=None, |
|
|
84 | 84 | ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False, |
|
|
85 |
xaxis="frequency", colormap='jet', normFactor=None , GauSelector = |
|
|
|
85 | xaxis="frequency", colormap='jet', normFactor=None , GauSelector = 0): | |
|
|
86 | 86 | |
|
|
87 | 87 | """ |
|
|
88 | 88 | |
| @@ -125,7 +125,7 class FitGauPlot(Figure): | |||
|
|
125 | 125 | x = dataOut.spc_range[1] |
|
|
126 | 126 | xlabel = "Time (ms)" |
|
|
127 | 127 | |
|
|
128 | else: | |
|
|
128 | else: | |
|
|
129 | 129 | x = dataOut.spc_range[2] |
|
|
130 | 130 | xlabel = "Velocity (m/s)" |
|
|
131 | 131 | |
| @@ -667,13 +667,14 class WindProfilerPlot(Figure): | |||
|
|
667 | 667 | #If there is a SNR function defined |
|
|
668 | 668 | if dataOut.data_SNR is not None: |
|
|
669 | 669 | nplots += 1 |
|
|
670 | SNR = dataOut.data_SNR | |
|
|
671 |
SNRavg = |
|
|
|
670 | SNR = dataOut.data_SNR[0] | |
|
|
671 | SNRavg = SNR#numpy.average(SNR, axis=0) | |
|
|
672 | 672 | |
|
|
673 | 673 | SNRdB = 10*numpy.log10(SNR) |
|
|
674 | 674 | SNRavgdB = 10*numpy.log10(SNRavg) |
|
|
675 | 675 | |
|
|
676 |
if SNRthresh == None: |
|
|
|
676 | if SNRthresh == None: | |
|
|
677 | SNRthresh = -5.0 | |
|
|
677 | 678 | ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0] |
|
|
678 | 679 | |
|
|
679 | 680 | for i in range(nplotsw): |
| @@ -27,8 +27,8 class SpectraPlot(Figure): | |||
|
|
27 | 27 | self.isConfig = False |
|
|
28 | 28 | self.__nsubplots = 1 |
|
|
29 | 29 | |
|
|
30 |
self.WIDTH = |
|
|
|
31 |
self.HEIGHT = |
|
|
|
30 | self.WIDTH = 300 | |
|
|
31 | self.HEIGHT = 300 | |
|
|
32 | 32 | self.WIDTHPROF = 120 |
|
|
33 | 33 | self.HEIGHTPROF = 0 |
|
|
34 | 34 | self.counter_imagwr = 0 |
| @@ -128,6 +128,10 class SpectraPlot(Figure): | |||
|
|
128 | 128 | factor = normFactor |
|
|
129 | 129 | if xaxis == "frequency": |
|
|
130 | 130 | x = dataOut.getFreqRange(1)/1000. |
|
|
131 | print '#######################################################' | |
|
|
132 | print 'xlen', len(x) | |
|
|
133 | print x | |
|
|
134 | print '#######################################################' | |
|
|
131 | 135 | xlabel = "Frequency (kHz)" |
|
|
132 | 136 | |
|
|
133 | 137 | elif xaxis == "time": |
| @@ -439,24 +443,45 class CrossSpectraPlot(Figure): | |||
|
|
439 | 443 | xlabel=xlabel, ylabel=ylabel, title=title, |
|
|
440 | 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 | 447 | coherence = numpy.abs(coherenceComplex) |
|
|
444 | 448 | # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi |
|
|
445 | 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 | 455 | # fig = plt.figure(10+self.indice) |
|
|
450 |
# plt.plot( x[0: |
|
|
|
456 | # #plt.plot( x[0:256],coherence[:,25] ) | |
|
|
457 | # cohAv = numpy.average(coherence,1) | |
|
|
458 | # | |
|
|
459 | # plt.plot( x[0:256],cohAv ) | |
|
|
451 | 460 | # #plt.axis([-12, 12, 15, 50]) |
|
|
452 | 461 | # plt.title("%s" %( '%s %s, Channel %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S") , i))) |
|
|
453 | 462 | # plt.ylabel('Desfase [grados]') |
|
|
454 | 463 | # plt.xlabel('Velocidad [m/s]') |
|
|
455 | 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 | 479 | # plt.show() |
|
|
458 | 480 | # self.indice=self.indice+1 |
|
|
459 | 481 | |
|
|
482 | ||
|
|
483 | ||
|
|
484 | ||
|
|
460 | 485 | title = "Coherence Ch%d * Ch%d" %(pair[0], pair[1]) |
|
|
461 | 486 | axes0 = self.axesList[i*self.__nsubplots+2] |
|
|
462 | 487 | axes0.pcolor(x, y, coherence, |
| @@ -500,7 +525,7 class RTIPlot(Figure): | |||
|
|
500 | 525 | self.__nsubplots = 1 |
|
|
501 | 526 | |
|
|
502 | 527 | self.WIDTH = 800 |
|
|
503 |
self.HEIGHT = |
|
|
|
528 | self.HEIGHT = 250 | |
|
|
504 | 529 | self.WIDTHPROF = 120 |
|
|
505 | 530 | self.HEIGHTPROF = 0 |
|
|
506 | 531 | self.counter_imagwr = 0 |
| @@ -110,6 +110,7 class BLTRParamReader(JRODataReader, ProcessingUnit): | |||
|
|
110 | 110 | status_value=0, |
|
|
111 | 111 | **kwargs): |
|
|
112 | 112 | |
|
|
113 | self.verbose = kwargs.get('verbose', True) | |
|
|
113 | 114 | self.path = path |
|
|
114 | 115 | self.startTime = startTime |
|
|
115 | 116 | self.endTime = endTime |
| @@ -214,17 +215,19 class BLTRParamReader(JRODataReader, ProcessingUnit): | |||
|
|
214 | 215 | self.readBlock() |
|
|
215 | 216 | |
|
|
216 | 217 | if (self.datatime.time() < self.startTime) or (self.datatime.time() > self.endTime): |
|
|
217 | print "[Reading] Record No. %d/%d -> %s [Skipping]" %( | |
|
|
218 | self.counter_records, | |
|
|
219 | self.nrecords, | |
|
|
220 |
self. |
|
|
|
218 | if self.verbose: | |
|
|
219 | print "[Reading] Record No. %d/%d -> %s [Skipping]" %( | |
|
|
220 | self.counter_records, | |
|
|
221 | self.nrecords, | |
|
|
222 | self.datatime.ctime()) | |
|
|
221 | 223 | continue |
|
|
222 | 224 | break |
|
|
223 | 225 | |
|
|
224 | print "[Reading] Record No. %d/%d -> %s" %( | |
|
|
225 | self.counter_records, | |
|
|
226 | self.nrecords, | |
|
|
227 | self.datatime.ctime()) | |
|
|
226 | if self.verbose: | |
|
|
227 | print "[Reading] Record No. %d/%d -> %s" %( | |
|
|
228 | self.counter_records, | |
|
|
229 | self.nrecords, | |
|
|
230 | self.datatime.ctime()) | |
|
|
228 | 231 | |
|
|
229 | 232 | return 1 |
|
|
230 | 233 | |
| @@ -1784,7 +1784,7 class JRODataWriter(JRODataIO): | |||
|
|
1784 | 1784 | |
|
|
1785 | 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 | 1789 | if not(self.isConfig): |
|
|
1790 | 1790 | |
| @@ -300,6 +300,7 class SpectraReader(JRODataReader, ProcessingUnit): | |||
|
|
300 | 300 | self.data_cspc = cspc['real'] + cspc['imag']*1j |
|
|
301 | 301 | else: |
|
|
302 | 302 | self.data_cspc = None |
|
|
303 | print 'SALE NONE ***********************************************************' | |
|
|
303 | 304 | |
|
|
304 | 305 | |
|
|
305 | 306 | if self.processingHeaderObj.flag_dc: |
| @@ -434,7 +435,7 class SpectraWriter(JRODataWriter, Operation): | |||
|
|
434 | 435 | |
|
|
435 | 436 | # dataOut = None |
|
|
436 | 437 | |
|
|
437 | def __init__(self): | |
|
|
438 | def __init__(self, **kwargs): | |
|
|
438 | 439 | """ |
|
|
439 | 440 | Inicializador de la clase SpectraWriter para la escritura de datos de espectros. |
|
|
440 | 441 | |
| @@ -449,7 +450,7 class SpectraWriter(JRODataWriter, Operation): | |||
|
|
449 | 450 | Return: None |
|
|
450 | 451 | """ |
|
|
451 | 452 | |
|
|
452 | Operation.__init__(self) | |
|
|
453 | Operation.__init__(self, **kwargs) | |
|
|
453 | 454 | |
|
|
454 | 455 | self.isConfig = False |
|
|
455 | 456 | |
| @@ -518,7 +519,7 class SpectraWriter(JRODataWriter, Operation): | |||
|
|
518 | 519 | |
|
|
519 | 520 | |
|
|
520 | 521 | def writeBlock(self): |
|
|
521 | """ | |
|
|
522 | """processingHeaderObj | |
|
|
522 | 523 | Escribe el buffer en el file designado |
|
|
523 | 524 | |
|
|
524 | 525 | |
| @@ -542,8 +543,10 class SpectraWriter(JRODataWriter, Operation): | |||
|
|
542 | 543 | data.tofile(self.fp) |
|
|
543 | 544 | |
|
|
544 | 545 | if self.data_cspc is not None: |
|
|
545 | data = numpy.zeros( self.shape_cspc_Buffer, self.dtype ) | |
|
|
546 | ||
|
|
546 | 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 | 550 | if not( self.processingHeaderObj.shif_fft ): |
|
|
548 | 551 | cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones |
|
|
549 | 552 | data['real'] = cspc.real |
| @@ -553,8 +556,9 class SpectraWriter(JRODataWriter, Operation): | |||
|
|
553 | 556 | |
|
|
554 | 557 | |
|
|
555 | 558 | if self.data_dc is not None: |
|
|
556 | data = numpy.zeros( self.shape_dc_Buffer, self.dtype ) | |
|
|
559 | ||
|
|
557 | 560 | dc = self.data_dc |
|
|
561 | data = numpy.zeros( numpy.shape(dc), self.dtype ) | |
|
|
558 | 562 | data['real'] = dc.real |
|
|
559 | 563 | data['imag'] = dc.imag |
|
|
560 | 564 | data = data.reshape((-1)) |
| 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 | 52 | break |
|
|
53 | 53 | return func.__get__(obj, cls) |
|
|
54 | 54 | |
|
|
55 | ||
|
|
56 | ||
|
|
57 | ||
|
|
58 | ||
|
|
59 | ||
|
|
60 | ||
|
|
61 | ||
|
|
62 | 55 | class ParametersProc(ProcessingUnit): |
|
|
63 | 56 | |
|
|
64 | 57 | nSeconds = None |
| @@ -125,11 +118,11 class ParametersProc(ProcessingUnit): | |||
|
|
125 | 118 | |
|
|
126 | 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 | 122 | print 'self.dataIn.data_spc', self.dataIn.data_spc.shape |
|
|
130 | 123 | self.dataOut.abscissaList = self.dataIn.getVelRange(1) |
|
|
131 | 124 | self.dataOut.spc_noise = self.dataIn.getNoise() |
|
|
132 |
self.dataOut.spc_range = (self.dataIn.getFreqRange(1) |
|
|
|
125 | self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1) ) | |
|
|
133 | 126 | |
|
|
134 | 127 | self.dataOut.normFactor = self.dataIn.normFactor |
|
|
135 | 128 | #self.dataOut.outputInterval = self.dataIn.outputInterval |
| @@ -142,9 +135,9 class ParametersProc(ProcessingUnit): | |||
|
|
142 | 135 | |
|
|
143 | 136 | print 'datain chandist ',self.dataOut.ChanDist |
|
|
144 | 137 | |
|
|
145 | if hasattr(self.dataIn, 'VelRange'): #Velocities range | |
|
|
146 | self.dataOut.VelRange = self.dataIn.VelRange | |
|
|
147 | else: self.dataOut.VelRange = None | |
|
|
138 | #if hasattr(self.dataIn, 'VelRange'): #Velocities range | |
|
|
139 | # self.dataOut.VelRange = self.dataIn.VelRange | |
|
|
140 | #else: self.dataOut.VelRange = None | |
|
|
148 | 141 | |
|
|
149 | 142 | if hasattr(self.dataIn, 'RadarConst'): #Radar Constant |
|
|
150 | 143 | self.dataOut.RadarConst = self.dataIn.RadarConst |
| @@ -193,6 +186,7 def target(tups): | |||
|
|
193 | 186 | #print 'TARGETTT', obj, args |
|
|
194 | 187 | return obj.FitGau(args) |
|
|
195 | 188 | |
|
|
189 | ||
|
|
196 | 190 | class GaussianFit(Operation): |
|
|
197 | 191 | |
|
|
198 | 192 | ''' |
| @@ -212,7 +206,7 class GaussianFit(Operation): | |||
|
|
212 | 206 | self.i=0 |
|
|
213 | 207 | |
|
|
214 | 208 | |
|
|
215 |
def run(self, dataOut, num_intg=7, pnoise=1., |
|
|
|
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 | 210 | """This routine will find a couple of generalized Gaussians to a power spectrum |
|
|
217 | 211 | input: spc |
|
|
218 | 212 | output: |
| @@ -239,12 +233,9 class GaussianFit(Operation): | |||
|
|
239 | 233 | #self.Num_Bin = len(self.spc) |
|
|
240 | 234 | self.Num_Bin = self.spc.shape[1] |
|
|
241 | 235 | self.Num_Chn = self.spc.shape[0] |
|
|
242 | ||
|
|
243 | 236 | Vrange = dataOut.abscissaList |
|
|
244 | 237 | |
|
|
245 | #print 'self.spc2', numpy.asarray(self.spc).shape | |
|
|
246 | ||
|
|
247 | GauSPC = numpy.empty([2,self.Num_Bin,self.Num_Hei]) | |
|
|
238 | GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei]) | |
|
|
248 | 239 | SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei]) |
|
|
249 | 240 | SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei]) |
|
|
250 | 241 | SPC_ch1[:] = numpy.NaN |
| @@ -256,264 +247,14 class GaussianFit(Operation): | |||
|
|
256 | 247 | noise_ = dataOut.spc_noise[0].copy() |
|
|
257 | 248 | |
|
|
258 | 249 | |
|
|
259 | ||
|
|
260 | 250 | pool = Pool(processes=self.Num_Chn) |
|
|
261 | 251 | args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)] |
|
|
262 | 252 | objs = [self for __ in range(self.Num_Chn)] |
|
|
263 | 253 | attrs = zip(objs, args) |
|
|
264 | 254 | gauSPC = pool.map(target, attrs) |
|
|
265 | 255 | dataOut.GauSPC = numpy.asarray(gauSPC) |
|
|
266 | # ret = [] | |
|
|
267 | # for n in range(self.Num_Chn): | |
|
|
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 | |
|
|
256 | ||
|
|
257 | ||
|
|
517 | 258 | |
|
|
518 | 259 | print '========================================================' |
|
|
519 | 260 | print 'total_time: ', time.time()-start_time |
| @@ -560,20 +301,22 class GaussianFit(Operation): | |||
|
|
560 | 301 | # normalizing spc and noise |
|
|
561 | 302 | # This part differs from gg1 |
|
|
562 | 303 | spc_norm_max = max(spc) |
|
|
563 | spc = spc / spc_norm_max | |
|
|
564 | pnoise = pnoise / spc_norm_max | |
|
|
304 | #spc = spc / spc_norm_max | |
|
|
305 | pnoise = pnoise #/ spc_norm_max | |
|
|
565 | 306 | ############################################# |
|
|
566 |
|
|
|
|
307 | ||
|
|
567 | 308 | fatspectra=1.0 |
|
|
568 | 309 | |
|
|
569 | wnoise = noise_ / spc_norm_max | |
|
|
310 | wnoise = noise_ #/ spc_norm_max | |
|
|
570 | 311 | #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used |
|
|
571 | 312 | #if wnoise>1.1*pnoise: # to be tested later |
|
|
572 | 313 | # wnoise=pnoise |
|
|
573 |
noisebl=wnoise*0.9; |
|
|
|
314 | noisebl=wnoise*0.9; | |
|
|
315 | noisebh=wnoise*1.1 | |
|
|
574 | 316 | spc=spc-wnoise |
|
|
575 | 317 | # print 'wnoise', noise_[0], spc_norm_max, wnoise |
|
|
576 | 318 | minx=numpy.argmin(spc) |
|
|
319 | #spcs=spc.copy() | |
|
|
577 | 320 | spcs=numpy.roll(spc,-minx) |
|
|
578 | 321 | cum=numpy.cumsum(spcs) |
|
|
579 | 322 | tot_noise=wnoise * self.Num_Bin #64; |
| @@ -597,7 +340,8 class GaussianFit(Operation): | |||
|
|
597 | 340 | #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4: |
|
|
598 | 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 | 345 | cumlo=cummax*epsi; |
|
|
602 | 346 | cumhi=cummax*(1-epsi) |
|
|
603 | 347 | powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0]) |
| @@ -619,7 +363,7 class GaussianFit(Operation): | |||
|
|
619 | 363 | x=numpy.arange( self.Num_Bin ) |
|
|
620 | 364 | y_data=spc+wnoise |
|
|
621 | 365 | |
|
|
622 |
|
|
|
|
366 | ''' single Gaussian ''' | |
|
|
623 | 367 | shift0=numpy.mod(midpeak+minx, self.Num_Bin ) |
|
|
624 | 368 | width0=powerwidth/4.#Initialization entire power of spectrum divided by 4 |
|
|
625 | 369 | power0=2. |
| @@ -629,16 +373,7 class GaussianFit(Operation): | |||
|
|
629 | 373 | lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True) |
|
|
630 | 374 | |
|
|
631 | 375 | chiSq1=lsq1[1]; |
|
|
632 | jack1= self.y_jacobian1(x,lsq1[0]) | |
|
|
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 | ||
|
|
376 | ||
|
|
642 | 377 | |
|
|
643 | 378 | if fatspectra<1.0 and powerwidth<4: |
|
|
644 | 379 | choice=0 |
| @@ -654,7 +389,7 class GaussianFit(Operation): | |||
|
|
654 | 389 | #return (numpy.array([shift0,width0,Amplitude0,p0]), |
|
|
655 | 390 | # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice) |
|
|
656 | 391 | |
|
|
657 |
|
|
|
|
392 | ''' two gaussians ''' | |
|
|
658 | 393 | #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64) |
|
|
659 | 394 | shift0=numpy.mod(firstpeak+minx, self.Num_Bin ); |
|
|
660 | 395 | shift1=numpy.mod(secondpeak+minx, self.Num_Bin ) |
| @@ -669,24 +404,16 class GaussianFit(Operation): | |||
|
|
669 | 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 | 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) | |
|
|
673 | ||
|
|
674 | ||
|
|
675 |
chiSq2=lsq2[1]; |
|
|
|
407 | lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True ) | |
|
|
408 | ||
|
|
409 | ||
|
|
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 | 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 | ||
|
|
689 |
if snrdB>- |
|
|
|
415 | ||
|
|
416 | if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error | |
|
|
690 | 417 | if oneG: |
|
|
691 | 418 | choice=0 |
|
|
692 | 419 | else: |
| @@ -696,7 +423,7 class GaussianFit(Operation): | |||
|
|
696 | 423 | s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; |
|
|
697 | 424 | s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2; |
|
|
698 | 425 | gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling |
|
|
699 | ||
|
|
426 | ||
|
|
700 | 427 | if gp1>gp2: |
|
|
701 | 428 | if a1>0.7*a2: |
|
|
702 | 429 | choice=1 |
| @@ -716,13 +443,15 class GaussianFit(Operation): | |||
|
|
716 | 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]) | |
|
|
720 |
|
|
|
|
446 | shift0=lsq2[0][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 = |
|
|
|
451 | max_vel = 1.0 | |
|
|
723 | 452 | |
|
|
724 | 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 | 455 | shift0=lsq2[0][0] |
|
|
727 | 456 | width0=lsq2[0][1] |
|
|
728 | 457 | Amplitude0=lsq2[0][2] |
| @@ -745,12 +474,12 class GaussianFit(Operation): | |||
|
|
745 | 474 | p0=lsq2[0][7] |
|
|
746 | 475 | noise=lsq2[0][8] |
|
|
747 | 476 | |
|
|
748 |
if Amplitude0<0. |
|
|
|
749 |
shift0,width0,Amplitude0,p0 = |
|
|
|
750 |
if Amplitude1<0. |
|
|
|
751 |
shift1,width1,Amplitude1,p1 = |
|
|
|
752 |
|
|
|
|
753 |
|
|
|
|
477 | if Amplitude0<0.05: # in case the peak is noise | |
|
|
478 | shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN] | |
|
|
479 | if Amplitude1<0.05: | |
|
|
480 | shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN] | |
|
|
481 | ||
|
|
482 | ||
|
|
754 | 483 | # if choice==0: # pick the single gaussian fit |
|
|
755 | 484 | # Amplitude0=lsq1[0][2] |
|
|
756 | 485 | # shift0=lsq1[0][0] |
| @@ -802,63 +531,8 class GaussianFit(Operation): | |||
|
|
802 | 531 | # print 'p1', p1 |
|
|
803 | 532 | # print 'noise', noise |
|
|
804 | 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 | |
|
|
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 | |
|
|
535 | return GauSPC | |
|
|
862 | 536 | |
|
|
863 | 537 | def y_model1(self,x,state): |
|
|
864 | 538 | shift0,width0,amplitude0,power0,noise=state |
| @@ -890,6 +564,7 class GaussianFit(Operation): | |||
|
|
890 | 564 | def misfit2(self,state,y_data,x,num_intg): |
|
|
891 | 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 | 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= |
|
|
|
912 |
tauW= |
|
|
|
586 | def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118, | |
|
|
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 | 590 | Velrange = dataOut.abscissaList |
|
|
920 | 591 | |
|
|
921 | 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 | 598 | Ze = self.dBZeMODE2(dataOut) |
|
|
924 | 599 | |
|
|
925 | 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 | 611 | self.Pt = Pt |
|
|
928 | 612 | self.Gt = Gt |
|
|
929 | 613 | self.Gr = Gr |
| @@ -933,48 +617,63 class PrecipitationProc(Operation): | |||
|
|
933 | 617 | self.ThetaT = ThetaT |
|
|
934 | 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 | 628 | SPCmean = numpy.mean(self.spc,0) |
|
|
938 | 629 | ETA = numpy.zeros(self.Num_Hei) |
|
|
939 | 630 | Pr = numpy.sum(SPCmean,0) |
|
|
940 | 631 | |
|
|
941 | #for R in range(self.Num_Hei): | |
|
|
942 | # ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA) | |
|
|
943 |
|
|
|
|
632 | VelMeteoro = numpy.mean(SPCmean,axis=0) | |
|
|
633 | print '==================== Vel SHAPE',VelMeteoro | |
|
|
634 | ||
|
|
944 | 635 | D_range = numpy.zeros(self.Num_Hei) |
|
|
636 | SIGMA = numpy.zeros(self.Num_Hei) | |
|
|
637 | N_dist = numpy.zeros(self.Num_Hei) | |
|
|
945 | 638 | EqSec = numpy.zeros(self.Num_Hei) |
|
|
946 | 639 | del_V = numpy.zeros(self.Num_Hei) |
|
|
947 | 640 | |
|
|
641 | ||
|
|
948 | 642 | for R in range(self.Num_Hei): |
|
|
949 | 643 | ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA) |
|
|
950 | 644 | |
|
|
951 | 645 | h = R + Altitude #Range from ground to radar pulse altitude |
|
|
952 | 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 - (Vel |
|
|
|
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 | 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 | 652 | N_dist[R] = ETA[R] / SIGMA[R] |
|
|
653 | ||
|
|
654 | ||
|
|
958 | 655 | |
|
|
959 | 656 | Ze = (ETA * Lambda**4) / (numpy.pi * Km) |
|
|
960 | 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 | 664 | dBRR = 10*numpy.log10(RR) |
|
|
966 | 665 | |
|
|
967 | 666 | dBZe = 10*numpy.log10(Ze) |
|
|
968 | 667 | dataOut.data_output = Ze |
|
|
969 |
dataOut.data_param = numpy.ones([ |
|
|
|
970 | dataOut.channelList = [0,1] | |
|
|
668 | dataOut.data_param = numpy.ones([3,self.Num_Hei]) | |
|
|
669 | dataOut.channelList = [0,1,2] | |
|
|
971 | 670 | print 'channelList', dataOut.channelList |
|
|
972 | 671 | dataOut.data_param[0]=dBZe |
|
|
973 | 672 | dataOut.data_param[1]=dBRR |
|
|
974 | 673 | print 'RR SHAPE', dBRR.shape |
|
|
975 | 674 | print 'Ze SHAPE', dBZe.shape |
|
|
976 | 675 | print 'dataOut.data_param SHAPE', dataOut.data_param.shape |
|
|
977 | ||
|
|
676 | ||
|
|
978 | 677 | |
|
|
979 | 678 | def dBZeMODE2(self, dataOut): # Processing for MIRA35C |
|
|
980 | 679 | |
| @@ -1001,26 +700,27 class PrecipitationProc(Operation): | |||
|
|
1001 | 700 | |
|
|
1002 | 701 | return Ze |
|
|
1003 | 702 | |
|
|
1004 | def GetRadarConstant(self): | |
|
|
1005 | ||
|
|
1006 | """ | |
|
|
1007 | Constants: | |
|
|
1008 | ||
|
|
1009 | Pt: Transmission Power dB 5kW | |
|
|
1010 | Gt: Transmission Gain dB 24.7 dB | |
|
|
1011 | Gr: Reception Gain dB 18.5 dB | |
|
|
1012 | Lambda: Wavelenght m 0.6741 m | |
|
|
1013 | aL: Attenuation loses dB | |
|
|
1014 | tauW: Width of transmission pulse s | |
|
|
1015 | ThetaT: Transmission antenna bean angle rad 0.1656317 rad | |
|
|
1016 | ThetaR: Reception antenna beam angle rad 0.36774087 rad | |
|
|
1017 | ||
|
|
1018 | """ | |
|
|
1019 | Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) | |
|
|
1020 | Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR) | |
|
|
1021 | RadarConstant = Numerator / Denominator | |
|
|
1022 | ||
|
|
1023 | return RadarConstant | |
|
|
703 | # def GetRadarConstant(self): | |
|
|
704 | # | |
|
|
705 | # """ | |
|
|
706 | # Constants: | |
|
|
707 | # | |
|
|
708 | # Pt: Transmission Power dB 5kW 5000 | |
|
|
709 | # Gt: Transmission Gain dB 24.7 dB 295.1209 | |
|
|
710 | # Gr: Reception Gain dB 18.5 dB 70.7945 | |
|
|
711 | # Lambda: Wavelenght m 0.6741 m 0.6741 | |
|
|
712 | # aL: Attenuation loses dB 4dB 2.5118 | |
|
|
713 | # tauW: Width of transmission pulse s 4us 4e-6 | |
|
|
714 | # ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317 | |
|
|
715 | # ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087 | |
|
|
716 | # | |
|
|
717 | # """ | |
|
|
718 | # | |
|
|
719 | # Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) | |
|
|
720 | # Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR) | |
|
|
721 | # RadarConstant = Numerator / Denominator | |
|
|
722 | # | |
|
|
723 | # return RadarConstant | |
|
|
1024 | 724 | |
|
|
1025 | 725 | |
|
|
1026 | 726 | |
| @@ -1045,8 +745,10 class FullSpectralAnalysis(Operation): | |||
|
|
1045 | 745 | """ |
|
|
1046 | 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() | |
|
|
1049 | cspc = dataOut.data_pre[1].copy() | |
|
|
748 | self.indice=int(numpy.random.rand()*1000) | |
|
|
749 | ||
|
|
750 | spc = dataOut.data_pre[0] | |
|
|
751 | cspc = dataOut.data_pre[1] | |
|
|
1050 | 752 | |
|
|
1051 | 753 | nChannel = spc.shape[0] |
|
|
1052 | 754 | nProfiles = spc.shape[1] |
| @@ -1060,10 +762,10 class FullSpectralAnalysis(Operation): | |||
|
|
1060 | 762 | |
|
|
1061 | 763 | #print 'ChanDist', ChanDist |
|
|
1062 | 764 | |
|
|
1063 | if dataOut.VelRange is not None: | |
|
|
1064 |
|
|
|
|
1065 | else: | |
|
|
1066 |
|
|
|
|
765 | #if dataOut.VelRange is not None: | |
|
|
766 | AbbsisaRange= dataOut.spc_range#dataOut.VelRange | |
|
|
767 | #else: | |
|
|
768 | # AbbsisaRange= dataOut.spc_range#dataOut.abscissaList | |
|
|
1067 | 769 | |
|
|
1068 | 770 | ySamples=numpy.ones([nChannel,nProfiles]) |
|
|
1069 | 771 | phase=numpy.ones([nChannel,nProfiles]) |
| @@ -1080,14 +782,16 class FullSpectralAnalysis(Operation): | |||
|
|
1080 | 782 | #print 'noise',noise |
|
|
1081 | 783 | #SNRdB = 10*numpy.log10(dataOut.data_SNR) |
|
|
1082 | 784 | |
|
|
1083 |
FirstMoment = |
|
|
|
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 | 788 | #SNRdBMean = [] |
|
|
1085 | 789 | |
|
|
1086 | 790 | |
|
|
1087 | 791 | #for j in range(nHeights): |
|
|
1088 | 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 | 793 | # SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]])) |
|
|
1090 |
|
|
|
|
794 | ||
|
|
1091 | 795 | data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN |
|
|
1092 | 796 | |
|
|
1093 | 797 | velocityX=[] |
| @@ -1097,21 +801,21 class FullSpectralAnalysis(Operation): | |||
|
|
1097 | 801 | |
|
|
1098 | 802 | dbSNR = 10*numpy.log10(dataSNR) |
|
|
1099 | 803 | dbSNR = numpy.average(dbSNR,0) |
|
|
804 | ||
|
|
1100 | 805 | for Height in range(nHeights): |
|
|
1101 | 806 | |
|
|
1102 |
[Vzon,Vmer,Vver, GaussCenter, PhaseSlope]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, |
|
|
|
1103 | ||
|
|
807 | [Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC]= self.WindEstimation(dataOut.data_param, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR[Height], SNRlimit) | |
|
|
1104 | 808 | PhaseLine = numpy.append(PhaseLine, PhaseSlope) |
|
|
1105 | 809 | |
|
|
1106 | 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 | 813 | else: |
|
|
1110 | 814 | #print 'Vzon',Vzon |
|
|
1111 | 815 | velocityX=numpy.append(velocityX, numpy.NaN) |
|
|
1112 | 816 | |
|
|
1113 | 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 | 820 | else: |
|
|
1117 | 821 | #print 'Vmer',Vmer |
| @@ -1129,24 +833,27 class FullSpectralAnalysis(Operation): | |||
|
|
1129 | 833 | |
|
|
1130 | 834 | |
|
|
1131 | 835 | |
|
|
1132 | data_output[0]=numpy.array(velocityX) | |
|
|
1133 | data_output[1]=numpy.array(velocityY) | |
|
|
1134 | data_output[2]=-velocityV#FirstMoment | |
|
|
836 | data_output[0] = numpy.array(velocityX) #self.moving_average(numpy.array(velocityX) , N=1) | |
|
|
837 | data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1) | |
|
|
838 | data_output[2] = -velocityV#FirstMoment | |
|
|
1135 | 839 | |
|
|
1136 | 840 | print ' ' |
|
|
1137 | 841 | #print 'FirstMoment' |
|
|
1138 | 842 | #print FirstMoment |
|
|
843 | print 'velocityX',numpy.shape(data_output[0]) | |
|
|
1139 | 844 | print 'velocityX',data_output[0] |
|
|
1140 | 845 | print ' ' |
|
|
846 | print 'velocityY',numpy.shape(data_output[1]) | |
|
|
1141 | 847 | print 'velocityY',data_output[1] |
|
|
848 | print 'velocityV',data_output[2] | |
|
|
1142 | 849 | print 'PhaseLine',PhaseLine |
|
|
1143 | 850 | #print numpy.array(velocityY) |
|
|
1144 | print ' ' | |
|
|
1145 | 851 | #print 'SNR' |
|
|
1146 | 852 | #print 10*numpy.log10(dataOut.data_SNR) |
|
|
1147 | 853 | #print numpy.shape(10*numpy.log10(dataOut.data_SNR)) |
|
|
1148 | 854 | print ' ' |
|
|
1149 | 855 | |
|
|
856 | xFrec=AbbsisaRange[0][0:spc.shape[1]] | |
|
|
1150 | 857 | |
|
|
1151 | 858 | dataOut.data_output=data_output |
|
|
1152 | 859 | |
| @@ -1156,15 +863,26 class FullSpectralAnalysis(Operation): | |||
|
|
1156 | 863 | def moving_average(self,x, N=2): |
|
|
1157 | 864 | return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):] |
|
|
1158 | 865 | |
|
|
1159 |
def gaus(self,xSamples, |
|
|
|
1160 |
return |
|
|
|
866 | def gaus(self,xSamples,Amp,Mu,Sigma): | |
|
|
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 | 887 | ySamples=numpy.ones([spc.shape[0],spc.shape[1]]) |
|
|
1170 | 888 | phase=numpy.ones([spc.shape[0],spc.shape[1]]) |
| @@ -1172,7 +890,14 class FullSpectralAnalysis(Operation): | |||
|
|
1172 | 890 | coherence=numpy.ones([spc.shape[0],spc.shape[1]]) |
|
|
1173 | 891 | PhaseSlope=numpy.zeros(spc.shape[0]) |
|
|
1174 | 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 | 902 | '''Getting Eij and Nij''' |
|
|
1178 | 903 | |
| @@ -1191,150 +916,243 class FullSpectralAnalysis(Operation): | |||
|
|
1191 | 916 | for i in range(spc.shape[0]): |
|
|
1192 | 917 | |
|
|
1193 | 918 | '''****** Line of Data SPC ******''' |
|
|
1194 | zline=z[i,:,Height] | |
|
|
919 | zline=z[i,:,Height].copy() - noise[i] # Se resta ruido | |
|
|
1195 | 920 | |
|
|
1196 | 921 | '''****** SPC is normalized ******''' |
|
|
1197 | FactNorm= (zline.copy()-noise[i]) / numpy.sum(zline.copy()) | |
|
|
1198 | FactNorm= FactNorm/numpy.sum(FactNorm) | |
|
|
1199 | ||
|
|
1200 | SmoothSPC=self.moving_average(FactNorm,N=3) | |
|
|
922 | SmoothSPC =self.moving_average(zline.copy(),N=1) # Se suaviza el ruido | |
|
|
923 | FactNorm = SmoothSPC/numpy.nansum(SmoothSPC) # SPC Normalizado y suavizado | |
|
|
1201 | 924 | |
|
|
1202 | xSamples = ar(range(len(SmoothSPC))) | |
|
|
1203 | ySamples[i] = SmoothSPC | |
|
|
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) | |
|
|
925 | xSamples = xFrec # Se toma el rango de frecuncias | |
|
|
926 | ySamples[i] = FactNorm # Se toman los valores de SPC normalizado | |
|
|
1216 | 927 | |
|
|
1217 | 928 | for i in range(spc.shape[0]): |
|
|
1218 | 929 | |
|
|
1219 | 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 | 936 | chan_index0 = pairsList[i][0] |
|
|
1224 | 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 | 942 | CSPCSamples[i] = CSPCNorm |
|
|
943 | ||
|
|
1230 | 944 | coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor) |
|
|
1231 | 945 | |
|
|
1232 |
coherence[i]= self.moving_average(coherence[i],N= |
|
|
|
946 | #coherence[i]= self.moving_average(coherence[i],N=1) | |
|
|
1233 | 947 | |
|
|
1234 | 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] | |
|
|
1237 | #print 'CSPCFactor', CSPCFactor#, CSPCFactor[0:20] | |
|
|
1238 | #print numpy.sum(ySamples[chan_index0]), numpy.sum(ySamples[chan_index1]), -noise[i] | |
|
|
1239 | #print 'CSPCNorm', CSPCNorm.shape, CSPCNorm[0:20] | |
|
|
1240 | #print 'CSPCNorm suma', numpy.sum(CSPCNorm) | |
|
|
1241 | #print 'CSPCSamples', CSPCSamples.shape, CSPCSamples[0,0:20] | |
|
|
950 | CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPCSamples[0]), xSamples), | |
|
|
951 | self.Moments(numpy.abs(CSPCSamples[1]), xSamples), | |
|
|
952 | self.Moments(numpy.abs(CSPCSamples[2]), xSamples)]) | |
|
|
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=[] | |
|
|
1246 | yMean2=[] | |
|
|
972 | mask01 = ~numpy.isnan(CSPCMask01) | |
|
|
973 | mask02 = ~numpy.isnan(CSPCMask02) | |
|
|
974 | mask12 = ~numpy.isnan(CSPCMask12) | |
|
|
1247 | 975 | |
|
|
1248 | for j in range(len(ySamples[1])): | |
|
|
1249 | yMean=numpy.append(yMean,numpy.mean([ySamples[0,j],ySamples[1,j],ySamples[2,j]])) | |
|
|
976 | #mask = ~numpy.isnan(CSPCMask01) | |
|
|
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): | |
|
|
1265 | if dbSNR > SNRlimit and abs(meanGauss/sigma**2) > 0.0001: | |
|
|
1266 | try: | |
|
|
1267 | popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=[1,meanGauss,sigma]) | |
|
|
985 | ||
|
|
986 | '''***Fit Gauss CSPC01***''' | |
|
|
987 | if dbSNR > SNRlimit: | |
|
|
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 | 1018 | except :#RuntimeError: |
|
|
1276 | 1019 | FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) |
|
|
1277 | 1020 | |
|
|
1278 | ||
|
|
1021 | ||
|
|
1279 | 1022 | else: |
|
|
1280 | 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 | 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] | |
|
|
1295 | if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0): | |
|
|
1296 | Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001 | |
|
|
1041 | if xSamples[PointFij] > xSamples[PointGauCenter]: | |
|
|
1042 | Fij = xSamples[PointFij] - xSamples[PointGauCenter] | |
|
|
1043 | ||
|
|
1297 | 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 ' ' | |
|
|
1312 | #print 'GCpos',GCpos, ( len(xFrec)- len(xFrec)*0.1) | |
|
|
1313 | #print 'Rangpos',Rangpos | |
|
|
1314 | print 'RANGE: ', Range | |
|
|
1315 | FrecRange=xFrec[Range[0]:Range[1]] | |
|
|
1100 | #GaussCenter = popt[1] #Primer momento 01 | |
|
|
1101 | GauWidth = popt[2] *3/2 #Ancho de banda de Gau01 | |
|
|
1102 | Range = numpy.empty(2) | |
|
|
1103 | Range[0] = GaussCenter - GauWidth | |
|
|
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 | 1122 | '''****** Getting SCPC Slope ******''' |
|
|
1318 | 1123 | |
|
|
1319 | 1124 | for i in range(spc.shape[0]): |
|
|
1320 | 1125 | |
|
|
1321 |
if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0. |
|
|
|
1126 | if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.6: | |
|
|
1322 | 1127 | PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3) |
|
|
1323 | 1128 | |
|
|
1324 | print 'FrecRange', len(FrecRange) , FrecRange | |
|
|
1325 | print 'PhaseRange', len(PhaseRange), PhaseRange | |
|
|
1326 | print ' ' | |
|
|
1129 | #print 'Ancho espectral Frecuencias', FrecRange[-1]-FrecRange[0], 'Hz' | |
|
|
1130 | #print 'Ancho espectral Velocidades', VelRange[-1]-VelRange[0], 'm/s' | |
|
|
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 | 1140 | if len(FrecRange) == len(PhaseRange): |
|
|
1328 | slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange) | |
|
|
1329 | PhaseSlope[i]=slope | |
|
|
1330 |
|
|
|
|
1141 | try: | |
|
|
1142 | slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange[mask], PhaseRange[mask]) | |
|
|
1143 | PhaseSlope[i]=slope | |
|
|
1144 | PhaseInter[i]=intercept | |
|
|
1145 | except: | |
|
|
1146 | PhaseSlope[i]=0 | |
|
|
1147 | PhaseInter[i]=0 | |
|
|
1331 | 1148 | else: |
|
|
1332 | 1149 | PhaseSlope[i]=0 |
|
|
1333 | 1150 | PhaseInter[i]=0 |
|
|
1334 | 1151 | else: |
|
|
1335 | 1152 | PhaseSlope[i]=0 |
|
|
1336 | 1153 | PhaseInter[i]=0 |
|
|
1337 | ||
|
|
1154 | ||
|
|
1155 | ||
|
|
1338 | 1156 | '''Getting constant C''' |
|
|
1339 | 1157 | cC=(Fij*numpy.pi)**2 |
|
|
1340 | 1158 | |
| @@ -1346,9 +1164,9 class FullSpectralAnalysis(Operation): | |||
|
|
1346 | 1164 | (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults) |
|
|
1347 | 1165 | |
|
|
1348 | 1166 | '''****** Getting constants A, B and H ******''' |
|
|
1349 | W01=numpy.amax(coherence[0]) | |
|
|
1350 | W02=numpy.amax(coherence[1]) | |
|
|
1351 | W12=numpy.amax(coherence[2]) | |
|
|
1167 | W01=numpy.nanmax( FitGauss01 ) #numpy.abs(CSPCSamples[0])) | |
|
|
1168 | W02=numpy.nanmax( FitGauss02 ) #numpy.abs(CSPCSamples[1])) | |
|
|
1169 | W12=numpy.nanmax( FitGauss12 ) #numpy.abs(CSPCSamples[2])) | |
|
|
1352 | 1170 | |
|
|
1353 | 1171 | WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC)) |
|
|
1354 | 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 | 1182 | VxVyResults=numpy.array([-cF,-cG]) |
|
|
1365 | 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 | 1190 | Vzon = Vy |
|
|
1368 | 1191 | Vmer = Vx |
|
|
1369 | 1192 | Vmag=numpy.sqrt(Vzon**2+Vmer**2) |
|
|
1370 | 1193 | Vang=numpy.arctan2(Vmer,Vzon) |
|
|
1371 |
Vver= |
|
|
|
1372 | print 'Height',Height | |
|
|
1373 | print 'vzon y vmer', Vzon, Vmer | |
|
|
1374 | return Vzon, Vmer, Vver, GaussCenter, PhaseSlope | |
|
|
1194 | Vver=SPCmoments[1] | |
|
|
1195 | FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12]) | |
|
|
1196 | ||
|
|
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 | 1236 | class SpectralMoments(Operation): |
|
|
1377 | 1237 | |
|
|
1378 | 1238 | ''' |
| @@ -1433,7 +1293,7 class SpectralMoments(Operation): | |||
|
|
1433 | 1293 | if (aliasing == None): aliasing = 0 |
|
|
1434 | 1294 | if (oldfd == None): oldfd = 0 |
|
|
1435 | 1295 | if (wwauto == None): wwauto = 0 |
|
|
1436 |
|
|
|
|
1296 | ||
|
|
1437 | 1297 | if (n0 < 1.e-20): n0 = 1.e-20 |
|
|
1438 | 1298 | |
|
|
1439 | 1299 | freq = oldfreq |
| @@ -1441,6 +1301,8 class SpectralMoments(Operation): | |||
|
|
1441 | 1301 | vec_fd = numpy.zeros(oldspec.shape[1]) |
|
|
1442 | 1302 | vec_w = numpy.zeros(oldspec.shape[1]) |
|
|
1443 | 1303 | vec_snr = numpy.zeros(oldspec.shape[1]) |
|
|
1304 | ||
|
|
1305 | oldspec = numpy.ma.masked_invalid(oldspec) | |
|
|
1444 | 1306 | |
|
|
1445 | 1307 | for ind in range(oldspec.shape[1]): |
|
|
1446 | 1308 | |
| @@ -1475,11 +1337,11 class SpectralMoments(Operation): | |||
|
|
1475 | 1337 | if (ss1 > m): ss1 = m |
|
|
1476 | 1338 | |
|
|
1477 | 1339 | valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1 |
|
|
1478 | power = ((spec2[valid] - n0)*fwindow[valid]).sum() | |
|
|
1479 | fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power | |
|
|
1340 | power = ( (spec2[valid] - n0) * fwindow[valid] ).sum() | |
|
|
1341 | fd = ( (spec2[valid]- n0) * freq[valid] * fwindow[valid] ).sum() / power | |
|
|
1480 | 1342 | w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power) |
|
|
1481 | 1343 | snr = (spec2.mean()-n0)/n0 |
|
|
1482 |
|
|
|
|
1344 | ||
|
|
1483 | 1345 | if (snr < 1.e-20) : |
|
|
1484 | 1346 | snr = 1.e-20 |
|
|
1485 | 1347 | |
| @@ -1487,7 +1349,7 class SpectralMoments(Operation): | |||
|
|
1487 | 1349 | vec_fd[ind] = fd |
|
|
1488 | 1350 | vec_w[ind] = w |
|
|
1489 | 1351 | vec_snr[ind] = snr |
|
|
1490 |
|
|
|
|
1352 | ||
|
|
1491 | 1353 | moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w)) |
|
|
1492 | 1354 | return moments |
|
|
1493 | 1355 | |
| @@ -1744,7 +1606,7 class SpectralFitting(Operation): | |||
|
|
1744 | 1606 | |
|
|
1745 | 1607 | fm = self.dataOut.library.modelFunction(p, constants) |
|
|
1746 | 1608 | fmp=numpy.dot(LT,fm) |
|
|
1747 |
|
|
|
|
1609 | ||
|
|
1748 | 1610 | return dp-fmp |
|
|
1749 | 1611 | |
|
|
1750 | 1612 | def __getSNR(self, z, noise): |
| @@ -1778,8 +1640,8 class WindProfiler(Operation): | |||
|
|
1778 | 1640 | |
|
|
1779 | 1641 | n = None |
|
|
1780 | 1642 | |
|
|
1781 | def __init__(self): | |
|
|
1782 | Operation.__init__(self) | |
|
|
1643 | def __init__(self, **kwargs): | |
|
|
1644 | Operation.__init__(self, **kwargs) | |
|
|
1783 | 1645 | |
|
|
1784 | 1646 | def __calculateCosDir(self, elev, azim): |
|
|
1785 | 1647 | zen = (90 - elev)*numpy.pi/180 |
| @@ -2293,7 +2155,7 class WindProfiler(Operation): | |||
|
|
2293 | 2155 | |
|
|
2294 | 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 | 2160 | param = dataOut.data_param |
|
|
2299 | 2161 | if dataOut.abscissaList != None: |
| @@ -4,6 +4,8 from jroproc_base import ProcessingUnit, Operation | |||
|
|
4 | 4 | from schainpy.model.data.jrodata import Spectra |
|
|
5 | 5 | from schainpy.model.data.jrodata import hildebrand_sekhon |
|
|
6 | 6 | |
|
|
7 | import matplotlib.pyplot as plt | |
|
|
8 | ||
|
|
7 | 9 | class SpectraProc(ProcessingUnit): |
|
|
8 | 10 | |
|
|
9 | 11 | def __init__(self, **kwargs): |
| @@ -275,7 +277,46 class SpectraProc(ProcessingUnit): | |||
|
|
275 | 277 | self.__selectPairsByChannel(self.dataOut.channelList) |
|
|
276 | 278 | |
|
|
277 | 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 | 320 | def selectHeights(self, minHei, maxHei): |
|
|
280 | 321 | """ |
|
|
281 | 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 | 333 | 1 si el metodo se ejecuto con exito caso contrario devuelve 0 |
|
|
293 | 334 | """ |
|
|
294 | 335 | |
|
|
336 | ||
|
|
295 | 337 | if (minHei > maxHei): |
|
|
296 | 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 | 361 | maxIndex = len(heights) |
|
|
320 | 362 | |
|
|
321 | 363 | self.selectHeightsByIndex(minIndex, maxIndex) |
|
|
364 | ||
|
|
322 | 365 | |
|
|
323 | 366 | return 1 |
|
|
324 | 367 | |
| @@ -361,6 +404,41 class SpectraProc(ProcessingUnit): | |||
|
|
361 | 404 | |
|
|
362 | 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 | 443 | def selectHeightsByIndex(self, minIndex, maxIndex): |
|
|
366 | 444 | """ |
| @@ -405,7 +483,8 class SpectraProc(ProcessingUnit): | |||
|
|
405 | 483 | self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1] |
|
|
406 | 484 | |
|
|
407 | 485 | return 1 |
|
|
408 | ||
|
|
486 | ||
|
|
487 | ||
|
|
409 | 488 | def removeDC(self, mode = 2): |
|
|
410 | 489 | jspectra = self.dataOut.data_spc |
|
|
411 | 490 | jcspectra = self.dataOut.data_cspc |
| @@ -463,6 +542,86 class SpectraProc(ProcessingUnit): | |||
|
|
463 | 542 | |
|
|
464 | 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 | 625 | def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None): |
|
|
467 | 626 | |
|
|
468 | 627 | jspectra = self.dataOut.data_spc |
General Comments 0
You need to be logged in to leave comments.
Login now