##// END OF EJS Templates
schain
RSS
A SpectraProc se agrega:...
Daniel Valdez -
r445:ff1e0d352ade
parent child
Show More
Show file before | Show file after | Hide comments
@@ -78,7 +78,7 def sorting_bruce(data, navg):
78 78
79 79 rtest = 1.0 + 1.0/navg
80 80
81 sum = 0.
81 sump = 0.
82 82
83 83 sumq = 0.
84 84
@@ -88,23 +88,23 def sorting_bruce(data, navg):
88 88
89 89 while((cont==1)and(j<lenOfData)):
90 90
91 sum += sortdata[j]
91 sump += sortdata[j]
92 92
93 93 sumq += sortdata[j]**2
94 94
95 95 j += 1
96 96
97 97 if j > nums_min:
98 if ((sumq*j) <= (rtest*sum**2)):
99 lnoise = sum / j
98 if ((sumq*j) <= (rtest*sump**2)):
99 lnoise = sump / j
100 100 else:
101 101 j = j - 1
102 sum = sum - sordata[j]
103 sumq = sumq - sordata[j]**2
102 sump = sump - sortdata[j]
103 sumq = sumq - sortdata[j]**2
104 104 cont = 0
105 105
106 106 if j == nums_min:
107 lnoise = sum /j
107 lnoise = sump /j
108 108
109 109 return lnoise
110 110
@@ -370,6 +370,8 class Spectra(JROData):
370 370
371 371 nCohInt = None #se requiere para determinar el valor de timeInterval
372 372
373 ippFactor = None
374
373 375 def __init__(self):
374 376 '''
375 377 Constructor
@@ -479,14 +481,14 class Spectra(JROData):
479 481
480 482 def getFreqRange(self, extrapoints=0):
481 483
482 deltafreq = self.getFmax() / self.nFFTPoints
484 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
483 485 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
484 486
485 487 return freqrange
486 488
487 489 def getVelRange(self, extrapoints=0):
488 490
489 deltav = self.getVmax() / self.nFFTPoints
491 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
490 492 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
491 493
492 494 return velrange
Show file before | Show file after | Hide comments
@@ -2766,69 +2766,6 class SpectraHeisWriter(Operation):
2766 2766 self.putData()
2767 2767
2768 2768
2769 class FITS:
2770 name=None
2771 format=None
2772 array =None
2773 data =None
2774 thdulist=None
2775 prihdr=None
2776 hdu=None
2777
2778 def __init__(self):
2779
2780 pass
2781
2782 def setColF(self,name,format,array):
2783 self.name=name
2784 self.format=format
2785 self.array=array
2786 a1=numpy.array([self.array],dtype=numpy.float32)
2787 self.col1 = pyfits.Column(name=self.name, format=self.format, array=a1)
2788 return self.col1
2789
2790 # def setColP(self,name,format,data):
2791 # self.name=name
2792 # self.format=format
2793 # self.data=data
2794 # a2=numpy.array([self.data],dtype=numpy.float32)
2795 # self.col2 = pyfits.Column(name=self.name, format=self.format, array=a2)
2796 # return self.col2
2797
2798
2799 def writeData(self,name,format,data):
2800 self.name=name
2801 self.format=format
2802 self.data=data
2803 a2=numpy.array([self.data],dtype=numpy.float32)
2804 self.col2 = pyfits.Column(name=self.name, format=self.format, array=a2)
2805 return self.col2
2806
2807 def cFImage(self,idblock,year,month,day,hour,minute,second):
2808 self.hdu= pyfits.PrimaryHDU(idblock)
2809 self.hdu.header.set("Year",year)
2810 self.hdu.header.set("Month",month)
2811 self.hdu.header.set("Day",day)
2812 self.hdu.header.set("Hour",hour)
2813 self.hdu.header.set("Minute",minute)
2814 self.hdu.header.set("Second",second)
2815 return self.hdu
2816
2817
2818 def Ctable(self,colList):
2819 self.cols=pyfits.ColDefs(colList)
2820 self.tbhdu = pyfits.new_table(self.cols)
2821 return self.tbhdu
2822
2823
2824 def CFile(self,hdu,tbhdu):
2825 self.thdulist=pyfits.HDUList([hdu,tbhdu])
2826
2827 def wFile(self,filename):
2828 if os.path.isfile(filename):
2829 os.remove(filename)
2830 self.thdulist.writeto(filename)
2831
2832 2769
2833 2770 class ParameterConf:
2834 2771 ELEMENTNAME = 'Parameter'
Show file before | Show file after | Hide comments
@@ -781,7 +781,7 class SpectraProc(ProcessingUnit):
781 781 self.dataOut.blockSize = blocksize
782 782 self.dataOut.flagShiftFFT = False
783 783
784 def init(self, nProfiles=None, nFFTPoints=None, pairsList=None):
784 def init(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None):
785 785
786 786 self.dataOut.flagNoData = True
787 787
@@ -799,6 +799,10 class SpectraProc(ProcessingUnit):
799 799 else:
800 800 nPairs = len(pairsList)
801 801
802 if ippFactor == None:
803 ippFactor = 1
804 self.dataOut.ippFactor = ippFactor
805
802 806 self.dataOut.nFFTPoints = nFFTPoints
803 807 self.dataOut.pairsList = pairsList
804 808 self.dataOut.nPairs = nPairs
@@ -828,7 +832,7 class SpectraProc(ProcessingUnit):
828 832
829 833 return
830 834
831 raise ValuError, "The type object %s is not valid"%(self.dataIn.type)
835 raise ValueError, "The type object %s is not valid"%(self.dataIn.type)
832 836
833 837 def selectChannels(self, channelList):
834 838
@@ -962,35 +966,228 class SpectraProc(ProcessingUnit):
962 966
963 967 return 1
964 968
965 def removeDC(self, mode = 1):
969 def removeDC(self, mode = 2):
970 jspectra = self.dataOut.data_spc
971 jcspectra = self.dataOut.data_cspc
966 972
967 dc_index = 0
968 freq_index = numpy.array([-2,-1,1,2])
969 data_spc = self.dataOut.data_spc
970 data_cspc = self.dataOut.data_cspc
971 data_dc = self.dataOut.data_dc
973
974 num_chan = jspectra.shape[0]
975 num_hei = jspectra.shape[2]
976
977 if jcspectra != None:
978 jcspectraExist = True
979 num_pairs = jcspectra.shape[0]
980 else: jcspectraExist = False
981
982 freq_dc = jspectra.shape[1]/2
983 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
984
985 if ind_vel[0]<0:
986 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
987
988 if mode == 1:
989 jspectra[:freq_dc,:] = (jspectra[:ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
990
991 if jcspectraExist:
992 jcspectra[:freq_dc,:] = (jcspectra[:ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
972 993
973 if self.dataOut.flagShiftFFT:
974 dc_index += self.dataOut.nFFTPoints/2
975 freq_index += self.dataOut.nFFTPoints/2
994 if mode == 2:
976 995
977 if mode == 1:
978 data_spc[dc_index] = (data_spc[:,freq_index[1],:] + data_spc[:,freq_index[2],:])/2
979 if data_cspc != None:
980 data_cspc[dc_index] = (data_cspc[:,freq_index[1],:] + data_cspc[:,freq_index[2],:])/2
981 return 1
996 vel = numpy.array([-2,-1,1,2])
997 xx = numpy.zeros([4,4])
998
999 for fil in range(4):
1000 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1001
1002 xx_inv = numpy.linalg.inv(xx)
1003 xx_aux = xx_inv[0,:]
1004
1005 for ich in range(num_chan):
1006 yy = jspectra[ich,ind_vel,:]
1007 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1008
1009 junkid = jspectra[ich,freq_dc,:]<=0
1010 cjunkid = sum(junkid)
1011
1012 if cjunkid.any():
1013 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1014
1015 if jcspectraExist:
1016 for ip in range(num_pairs):
1017 yy = jcspectra[ip,ind_vel,:]
1018 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
982 1019
983 if mode == 2:
984 pass
985 1020
986 if mode == 3:
987 pass
988
989 raise ValueError, "mode parameter has to be 1, 2 or 3"
1021 self.dataOut.data_spc = jspectra
1022 self.dataOut.data_cspc = jcspectra
1023
1024 return 1
990 1025
991 def removeInterference(self):
1026 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
1027
1028 jspectra = self.dataOut.data_spc
1029 jcspectra = self.dataOut.data_cspc
1030 jnoise = self.dataOut.getNoise()
1031 num_incoh = self.dataOut.nIncohInt
1032
1033 num_channel = jspectra.shape[0]
1034 num_prof = jspectra.shape[1]
1035 num_hei = jspectra.shape[2]
1036
1037 #hei_interf
1038 if hei_interf == None:
1039 count_hei = num_hei/2 #Como es entero no importa
1040 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
1041 hei_interf = numpy.asarray(hei_interf)[0]
1042 #nhei_interf
1043 if (nhei_interf == None):
1044 nhei_interf = 5
1045 if (nhei_interf < 1):
1046 nhei_interf = 1
1047 if (nhei_interf > count_hei):
1048 nhei_interf = count_hei
1049 if (offhei_interf == None):
1050 offhei_interf = 0
1051
1052 ind_hei = range(num_hei)
1053 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1054 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1055 mask_prof = numpy.asarray(range(num_prof))
1056 num_mask_prof = mask_prof.size
1057 comp_mask_prof = [0, num_prof/2]
1058
1059
1060 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1061 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1062 jnoise = numpy.nan
1063 noise_exist = jnoise[0] < numpy.Inf
1064
1065 #Subrutina de Remocion de la Interferencia
1066 for ich in range(num_channel):
1067 #Se ordena los espectros segun su potencia (menor a mayor)
1068 power = jspectra[ich,mask_prof,:]
1069 power = power[:,hei_interf]
1070 power = power.sum(axis = 0)
1071 psort = power.ravel().argsort()
1072
1073 #Se estima la interferencia promedio en los Espectros de Potencia empleando
1074 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
1075
1076 if noise_exist:
1077 # tmp_noise = jnoise[ich] / num_prof
1078 tmp_noise = jnoise[ich]
1079 junkspc_interf = junkspc_interf - tmp_noise
1080 #junkspc_interf[:,comp_mask_prof] = 0
1081
1082 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
1083 jspc_interf = jspc_interf.transpose()
1084 #Calculando el espectro de interferencia promedio
1085 noiseid = numpy.where(jspc_interf <= tmp_noise/ math.sqrt(num_incoh))
1086 noiseid = noiseid[0]
1087 cnoiseid = noiseid.size
1088 interfid = numpy.where(jspc_interf > tmp_noise/ math.sqrt(num_incoh))
1089 interfid = interfid[0]
1090 cinterfid = interfid.size
1091
1092 if (cnoiseid > 0): jspc_interf[noiseid] = 0
1093
1094 #Expandiendo los perfiles a limpiar
1095 if (cinterfid > 0):
1096 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
1097 new_interfid = numpy.asarray(new_interfid)
1098 new_interfid = {x for x in new_interfid}
1099 new_interfid = numpy.array(list(new_interfid))
1100 new_cinterfid = new_interfid.size
1101 else: new_cinterfid = 0
1102
1103 for ip in range(new_cinterfid):
1104 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
1105 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
1106
1107
1108 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
1109
1110 #Removiendo la interferencia del punto de mayor interferencia
1111 ListAux = jspc_interf[mask_prof].tolist()
1112 maxid = ListAux.index(max(ListAux))
1113
1114
1115 if cinterfid > 0:
1116 for ip in range(cinterfid*(interf == 2) - 1):
1117 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/math.sqrt(num_incoh))).nonzero()
1118 cind = len(ind)
1119
1120 if (cind > 0):
1121 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/math.sqrt(num_incoh))
1122
1123 ind = numpy.array([-2,-1,1,2])
1124 xx = numpy.zeros([4,4])
1125
1126 for id1 in range(4):
1127 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
1128
1129 xx_inv = numpy.linalg.inv(xx)
1130 xx = xx_inv[:,0]
1131 ind = (ind + maxid + num_mask_prof)%num_mask_prof
1132 yy = jspectra[ich,mask_prof[ind],:]
1133 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
1134
1135
1136 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/math.sqrt(num_incoh))).nonzero()
1137 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/math.sqrt(num_incoh))
1138
1139 #Remocion de Interferencia en el Cross Spectra
1140 if jcspectra == None: return jspectra, jcspectra
1141 num_pairs = jcspectra.size/(num_prof*num_hei)
1142 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
992 1143
993 pass
1144 for ip in range(num_pairs):
1145
1146 #-------------------------------------------
1147
1148 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
1149 cspower = cspower[:,hei_interf]
1150 cspower = cspower.sum(axis = 0)
1151
1152 cspsort = cspower.ravel().argsort()
1153 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
1154 junkcspc_interf = junkcspc_interf.transpose()
1155 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
1156
1157 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1158
1159 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
1160 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
1161 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
1162
1163 for iprof in range(num_prof):
1164 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
1165 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
1166
1167 #Removiendo la Interferencia
1168 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
1169
1170 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1171 maxid = ListAux.index(max(ListAux))
1172
1173 ind = numpy.array([-2,-1,1,2])
1174 xx = numpy.zeros([4,4])
1175
1176 for id1 in range(4):
1177 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
1178
1179 xx_inv = numpy.linalg.inv(xx)
1180 xx = xx_inv[:,0]
1181
1182 ind = (ind + maxid + num_mask_prof)%num_mask_prof
1183 yy = jcspectra[ip,mask_prof[ind],:]
1184 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
1185
1186 #Guardar Resultados
1187 self.dataOut.data_spc = jspectra
1188 self.dataOut.data_cspc = jcspectra
1189
1190 return 1
994 1191
995 1192 def setRadarFrequency(self, frequency=None):
996 1193 if frequency != None:
@@ -1470,7 +1667,7 class SpectraHeisProc(ProcessingUnit):
1470 1667
1471 1668 return
1472 1669
1473 raise ValuError, "The type object %s is not valid"%(self.dataIn.type)
1670 raise ValueError, "The type object %s is not valid"%(self.dataIn.type)
1474 1671
1475 1672
1476 1673 def selectChannels(self, channelList):
General Comments 0
You need to be logged in to leave comments. Login now