##// END OF EJS Templates
schain
RSS
Update Clase Oblique_Gauss_Fit in jroproc_parameters.py
Alexander Valdez -
r1686:b55c513e2e62
parent child
Show More
Side by Side Unified
Context file:
r1686:b55c513e2e62 Loading diff...:
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1075 lines changed) Show them Hide them
@@ -23,7 +23,6 import warnings
23 from numpy import NaN
23 from numpy import NaN
24 from scipy.optimize.optimize import OptimizeWarning
24 from scipy.optimize.optimize import OptimizeWarning
25 warnings.filterwarnings('ignore')
25 warnings.filterwarnings('ignore')
26 import pdb
27
26
28
27
29 SPEED_OF_LIGHT = 299792458
28 SPEED_OF_LIGHT = 299792458
@@ -212,7 +211,7 class RemoveWideGC(Operation):
212 self.i = 0
211 self.i = 0
213 self.ich = 0
212 self.ich = 0
214 self.ir = 0
213 self.ir = 0
215
214
216 def run(self, dataOut, ClutterWidth=2.5):
215 def run(self, dataOut, ClutterWidth=2.5):
217
216
218 self.spc = dataOut.data_pre[0].copy()
217 self.spc = dataOut.data_pre[0].copy()
@@ -238,11 +237,11 class RemoveWideGC(Operation):
238 junk = numpy.append(numpy.insert(numpy.squeeze(self.spc[ich,gc_values,ir]),0,HSn),HSn)
237 junk = numpy.append(numpy.insert(numpy.squeeze(self.spc[ich,gc_values,ir]),0,HSn),HSn)
239 j1index = numpy.squeeze(numpy.where(numpy.diff(junk)>0))
238 j1index = numpy.squeeze(numpy.where(numpy.diff(junk)>0))
240 j2index = numpy.squeeze(numpy.where(numpy.diff(junk)<0))
239 j2index = numpy.squeeze(numpy.where(numpy.diff(junk)<0))
241 if ((numpy.size(j1index)<=1) | (numpy.size(j2index)<=1)) :
240 if ((numpy.size(j1index)<=1) | (numpy.size(j2index)<=1)) :
242 continue
241 continue
243 junk3 = numpy.squeeze(numpy.diff(j1index))
242 junk3 = numpy.squeeze(numpy.diff(j1index))
244 junk4 = numpy.squeeze(numpy.diff(j2index))
243 junk4 = numpy.squeeze(numpy.diff(j2index))
245
244
246 valleyindex = j2index[numpy.where(junk4>1)]
245 valleyindex = j2index[numpy.where(junk4>1)]
247 peakindex = j1index[numpy.where(junk3>1)]
246 peakindex = j1index[numpy.where(junk3>1)]
248
247
@@ -252,7 +251,7 class RemoveWideGC(Operation):
252 if numpy.size(isvalid) >1 :
251 if numpy.size(isvalid) >1 :
253 vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])
252 vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])
254 isvalid = isvalid[vindex]
253 isvalid = isvalid[vindex]
255
254
256 # clutter peak
255 # clutter peak
257 gcpeak = peakindex[isvalid]
256 gcpeak = peakindex[isvalid]
258 vl = numpy.where(valleyindex < gcpeak)
257 vl = numpy.where(valleyindex < gcpeak)
@@ -274,7 +273,7 class RemoveWideGC(Operation):
274 return dataOut
273 return dataOut
275
274
276 class SpectralFilters(Operation):
275 class SpectralFilters(Operation):
277 ''' This class allows to replace the novalid values with noise for each channel
276 ''' This class allows to replace the novalid values with noise for each channel
278 This applies to CLAIRE RADAR
277 This applies to CLAIRE RADAR
279
278
280 PositiveLimit : RightLimit of novalid data
279 PositiveLimit : RightLimit of novalid data
@@ -663,7 +662,9 class GaussianFit(Operation):
663 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
662 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
664
663
665 class Oblique_Gauss_Fit(Operation):
664 class Oblique_Gauss_Fit(Operation):
666
665 '''
666 Written by R. Flores
667 '''
667 def __init__(self):
668 def __init__(self):
668 Operation.__init__(self)
669 Operation.__init__(self)
669
670
@@ -714,7 +715,6 class Oblique_Gauss_Fit(Operation):
714
715
715 return gaussian(x, popt[0], popt[1], popt[2], popt[3]), popt[0], popt[1], popt[2], popt[3]
716 return gaussian(x, popt[0], popt[1], popt[2], popt[3]), popt[0], popt[1], popt[2], popt[3]
716
717
717
718 def Gauss_fit_2(self,spc,x,nGauss):
718 def Gauss_fit_2(self,spc,x,nGauss):
719
719
720
720
@@ -750,13 +750,7 class Oblique_Gauss_Fit(Operation):
750 print("ERROR")
750 print("ERROR")
751
751
752 d = numpy.mean(y[-100:])
752 d = numpy.mean(y[-100:])
753
754 # define a least squares function to optimize
755 popt,pcov = curve_fit(gaussian,x,y,p0=[a,b,c,d])
753 popt,pcov = curve_fit(gaussian,x,y,p0=[a,b,c,d])
756 #popt,fopt,niter,funcalls = fmin(minfunc,[a,b,c,d])
757
758
759 #return gaussian(x, popt[0], popt[1], popt[2], popt[3]), popt[0], popt[1], popt[2], popt[3]
760 return gaussian(x, popt[0], popt[1], popt[2], popt[3]),popt[0], popt[1], popt[2], popt[3]
754 return gaussian(x, popt[0], popt[1], popt[2], popt[3]),popt[0], popt[1], popt[2], popt[3]
761
755
762 def Double_Gauss_fit(self,spc,x,A1,B1,C1,A2,B2,C2,D):
756 def Double_Gauss_fit(self,spc,x,A1,B1,C1,A2,B2,C2,D):
@@ -807,14 +801,675 class Oblique_Gauss_Fit(Operation):
807 d = D
801 d = D
808
802
809 # fit
803 # fit
810
811 popt,pcov = curve_fit(double_gaussian,x,y,p0=[a1,b1,c1,a2,b2,c2,d])
804 popt,pcov = curve_fit(double_gaussian,x,y,p0=[a1,b1,c1,a2,b2,c2,d])
805 error = numpy.sqrt(numpy.diag(pcov))
806
807 return popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6], error[0], error[1], error[2], error[3], error[4], error[5], error[6]
808
809 def windowing_double(self,spc,x,A1,B1,C1,A2,B2,C2,D):
810 from scipy.optimize import curve_fit,fmin
811
812 def R_gaussian(x, a, b, c):
813 N = int(numpy.shape(x)[0])
814 val = a * numpy.exp(-((x)*c*2*2*numpy.pi)**2 / (2))* numpy.exp(1.j*b*x*4*numpy.pi)
815 return val
816
817 def T(x,N):
818 T = 1-abs(x)/N
819 return T
820
821 def R_T_spc_fun(x, a1, b1, c1, a2, b2, c2, d):
822
823 N = int(numpy.shape(x)[0])
824
825 x_max = x[-1]
826
827 x_pos = x[1600:]
828 x_neg = x[:1600]
829
830 R_T_neg_1 = R_gaussian(x, a1, b1, c1)[:1600]*T(x_neg,-x[0])
831 R_T_pos_1 = R_gaussian(x, a1, b1, c1)[1600:]*T(x_pos,x[-1])
832 R_T_sum_1 = R_T_pos_1 + R_T_neg_1
833 R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real
834 R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1)
835 max_val_1 = numpy.max(R_T_spc_1)
836 R_T_spc_1 = R_T_spc_1*a1/max_val_1
837
838 R_T_neg_2 = R_gaussian(x, a2, b2, c2)[:1600]*T(x_neg,-x[0])
839 R_T_pos_2 = R_gaussian(x, a2, b2, c2)[1600:]*T(x_pos,x[-1])
840 R_T_sum_2 = R_T_pos_2 + R_T_neg_2
841 R_T_spc_2 = numpy.fft.fft(R_T_sum_2).real
842 R_T_spc_2 = numpy.fft.fftshift(R_T_spc_2)
843 max_val_2 = numpy.max(R_T_spc_2)
844 R_T_spc_2 = R_T_spc_2*a2/max_val_2
845
846 R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N))
847 R_T_d_neg = R_T_d[:1600]*T(x_neg,-x[0])
848 R_T_d_pos = R_T_d[1600:]*T(x_pos,x[-1])
849 R_T_d_sum = R_T_d_pos + R_T_d_neg
850 R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real
851 R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3)
852
853 R_T_final = R_T_spc_1 + R_T_spc_2 + R_T_spc_3
854
855 return R_T_final
856
857 y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))
858
859 from scipy.stats import norm
860 mean,std=norm.fit(spc)
861
862 # estimate starting values from the data
863 a1 = A1
864 b1 = B1
865 c1 = C1#numpy.std(spc)
866
867 a2 = A2#y.max()
868 b2 = B2#x[numpy.argmax(y)]
869 c2 = C2#numpy.std(spc)
870 d = D
871
872 ippSeconds = 250*20*1.e-6/3
873
874 x_t = ippSeconds * (numpy.arange(1600) -1600 / 2.)
875
876 x_t = numpy.linspace(x_t[0],x_t[-1],3200)
877
878 x_freq = numpy.fft.fftfreq(1600,d=ippSeconds)
879 x_freq = numpy.fft.fftshift(x_freq)
880
881 # define a least squares function to optimize
882 def minfunc(params):
883 return sum((y-R_T_spc_fun(x_t,params[0],params[1],params[2],params[3],params[4],params[5],params[6]))**2/1)#y**2)
884
885 # fit
886 popt_full = fmin(minfunc,[a1,b1,c1,a2,b2,c2,d],full_output=True)
887 popt = popt_full[0]
888
889 return popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]
890
891 def Double_Gauss_fit_weight(self,spc,x,A1,B1,C1,A2,B2,C2,D):
892 from scipy.optimize import curve_fit,fmin
893
894 def double_gaussian(x, a1, b1, c1, a2, b2, c2, d):
895 val = a1 * numpy.exp(-(x - b1)**2 / (2*c1**2)) + a2 * numpy.exp(-(x - b2)**2 / (2*c2**2)) + d
896 return val
897
898 y = spc
899
900 from scipy.stats import norm
901 mean,std=norm.fit(spc)
902
903 # estimate starting values from the data
904 a1 = A1
905 b1 = B1
906 c1 = C1#numpy.std(spc)
907
908 a2 = A2#y.max()
909 b2 = B2#x[numpy.argmax(y)]
910 c2 = C2#numpy.std(spc)
911 d = D
912
913 y_clean = signal.medfilt(y)
914 # define a least squares function to optimize
915 def minfunc(params):
916 return sum((y-double_gaussian(x,params[0],params[1],params[2],params[3],params[4],params[5],params[6]))**2/(y_clean**2/1))
917
918 # fit
919 popt_full = fmin(minfunc,[a1,b1,c1,a2,b2,c2,d], disp =False, full_output=True)
920 #print("nIter", popt_full[2])
921 popt = popt_full[0]
922 #popt,pcov = curve_fit(double_gaussian,x,y,p0=[a1,b1,c1,a2,b2,c2,d])
923
924 #return double_gaussian(x, popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]), popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]
925 return popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]
926
927 def DH_mode(self,spectra,VelRange):
928
929 from scipy.optimize import curve_fit
930
931 def double_gauss(x, a1,b1,c1, a2,b2,c2, d):
932 val = a1 * numpy.exp(-(x - b1)**2 / (2*c1**2)) + a2 * numpy.exp(-(x - b2)**2 / (2*c2**2)) + d
933 return val
934
935 spec = (spectra.copy()).flatten()
936 amp=spec.max()
937 params=numpy.array([amp,-400,30,amp/4,-200,150,1.0e7])
938 #try:
939 popt,pcov=curve_fit(double_gauss, VelRange, spec, p0=params,bounds=([0,-460,0,0,-400,120,0],[numpy.inf,-340,50,numpy.inf,0,250,numpy.inf]))
812
940
813 error = numpy.sqrt(numpy.diag(pcov))
941 error = numpy.sqrt(numpy.diag(pcov))
942 #doppler_2=popt[4]
943 #err_2 = numpy.sqrt(pcov[4][4])
944
945 #except:
946 #pass
947 #doppler_2=numpy.NAN
948 #err_2 = numpy.NAN
949
950 #return doppler_2, err_2
814
951
815 return popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6], error[0], error[1], error[2], error[3], error[4], error[5], error[6]
952 return popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6], error[0], error[1], error[2], error[3], error[4], error[5], error[6]
816
953
817 def run(self, dataOut):
954 def Tri_Marco(self,spc,freq,a1,b1,c1,a2,b2,c2,d):
955
956 from scipy.optimize import least_squares
957
958 freq_max = numpy.max(numpy.abs(freq))
959 spc_max = numpy.max(spc)
960
961 def tri_gaussian(x, a1, b1, c1, a2, b2, c2, a3, b3, c3, d):
962 z1 = (x-b1)/c1
963 z2 = (x-b2)/c2
964 z3 = (x-b3)/c3
965 val = a1 * numpy.exp(-z1**2/2) + a2 * numpy.exp(-z2**2/2) + a3 * numpy.exp(-z3**2/2) + d
966 return val
967
968 from scipy.signal import medfilt
969 Nincoh = 20
970 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
971 c1 = abs(c1)
972 c2 = abs(c2)
973
974 # define a least squares function to optimize
975 def lsq_func(params):
976 return (spc-tri_gaussian(freq,params[0],params[1],params[2],params[3],params[4],params[5],params[6],params[7],params[8],params[9]))/spcm
977
978 # fit
979 bounds=([0,-numpy.inf,0,0,-numpy.inf,0,0,0,0,0],[numpy.inf,-100,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf,600,numpy.inf,numpy.inf])
980
981 params_scale = [spc_max,freq_max,freq_max,spc_max,freq_max,freq_max,spc_max,freq_max,freq_max,spc_max]
982 #print(a1,b1,c1,a2,b2,c2,d)
983 popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,a2/4,-b1,c1,d],x_scale=params_scale,bounds=bounds)
984
985 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]
986 A2f = popt.x[3]; B2f = popt.x[4]; C2f = popt.x[5]
987 A3f = popt.x[6]; B3f = popt.x[7]; C3f = popt.x[8]
988 Df = popt.x[9]
989
990 return A1f, B1f, C1f, A2f, B2f, C2f, Df
991
992 def Tri_Marco(self,spc,freq,a1,b1,c1,a2,b2,c2,d):
993
994 from scipy.optimize import least_squares
995
996 freq_max = numpy.max(numpy.abs(freq))
997 spc_max = numpy.max(spc)
998
999 def duo_gaussian(x, a1, b1, c1, a2, b2, c2, d):
1000 z1 = (x-b1)/c1
1001 z2 = (x-b2)/c2
1002 #z3 = (x-b3)/c3
1003 val = a1 * numpy.exp(-z1**2/2) + a2 * numpy.exp(-z2**2/2) + d
1004 return val
1005
1006 from scipy.signal import medfilt
1007 Nincoh = 20
1008 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1009 c1 = abs(c1)
1010 c2 = abs(c2)
1011
1012 # define a least squares function to optimize
1013 def lsq_func(params):
1014 return (spc-tri_gaussian(freq,params[0],params[1],params[2],params[3],params[4],params[5],params[6]))/spcm
1015
1016 # fit
1017 bounds=([0,-numpy.inf,0,0,-numpy.inf,0,0],[numpy.inf,-100,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf])
1018
1019 params_scale = [spc_max,freq_max,freq_max,spc_max,freq_max,freq_max,spc_max]
1020 popt = least_squares(lsq_func,[a1,b1,c1,a2,b2,c2,d],x_scale=params_scale,bounds=bounds)
1021
1022 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]
1023 A2f = popt.x[3]; B2f = popt.x[4]; C2f = popt.x[5]
1024 #A3f = popt.x[6]; B3f = popt.x[7]; C3f = popt.x[8]
1025 Df = popt.x[9]
1026
1027 return A1f, B1f, C1f, A2f, B2f, C2f, Df
1028
1029 def double_gaussian_skew(self,x, a1, b1, c1, a2, b2, c2, k2, d):
1030 z1 = (x-b1)/c1
1031 z2 = (x-b2)/c2
1032 h2 = 1-k2*z2
1033 h2[h2<0] = 0
1034 y2 = -1/k2*numpy.log(h2)
1035 val = a1 * numpy.exp(-z1**2/2) + a2 * numpy.exp(-y2**2/2)/(1-k2*z2) + d
1036 return val
1037
1038 def gaussian(self, x, a, b, c, d):
1039 z = (x-b)/c
1040 val = a * numpy.exp(-z**2/2) + d
1041 return val
1042
1043 def double_gaussian(self, x, a1, b1, c1, a2, b2, c2, d):
1044 z1 = (x-b1)/c1
1045 z2 = (x-b2)/c2
1046 val = a1 * numpy.exp(-z1**2/2) + a2 * numpy.exp(-z2**2/2) + d
1047 return val
1048
1049 def double_gaussian_double_skew(self,x, a1, b1, c1, k1, a2, b2, c2, k2, d):
1050
1051 z1 = (x-b1)/c1
1052 h1 = 1-k1*z1
1053 h1[h1<0] = 0
1054 y1 = -1/k1*numpy.log(h1)
1055
1056 z2 = (x-b2)/c2
1057 h2 = 1-k2*z2
1058 h2[h2<0] = 0
1059 y2 = -1/k2*numpy.log(h2)
1060
1061 val = a1 * numpy.exp(-y1**2/2)/(1-k1*z1) + a2 * numpy.exp(-y2**2/2)/(1-k2*z2) + d
1062 return val
1063
1064 def gaussian_skew(self,x, a2, b2, c2, k2, d):
1065 z2 = (x-b2)/c2
1066 h2 = 1-k2*z2
1067 h2[h2<0] = 0
1068 y2 = -1/k2*numpy.log(h2)
1069 val = a2 * numpy.exp(-y2**2/2)/(1-k2*z2) + d
1070 return val
1071
1072 def triple_gaussian_skew(self,x, a1, b1, c1, a2, b2, c2, k2, a3, b3, c3, k3, d):
1073 z1 = (x-b1)/c1
1074 z2 = (x-b2)/c2
1075 z3 = (x-b3)/c3
1076 h2 = 1-k2*z2
1077 h2[h2<0] = 0
1078 y2 = -1/k2*numpy.log(h2)
1079 h3 = 1-k3*z3
1080 h3[h3<0] = 0
1081 y3 = -1/k3*numpy.log(h3)
1082 val = a1 * numpy.exp(-z1**2/2) + a2 * numpy.exp(-y2**2/2)/(1-k2*z2) + a3 * numpy.exp(-y3**2/2)/(1-k3*z3) + d
1083 return val
1084
1085 def Double_Gauss_Skew_fit_weight_bound_no_inputs(self,spc,freq):
1086
1087 from scipy.optimize import least_squares
1088
1089 freq_max = numpy.max(numpy.abs(freq))
1090 spc_max = numpy.max(spc)
1091
1092 from scipy.signal import medfilt
1093 Nincoh = 20
1094 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1095
1096 # define a least squares function to optimize
1097 def lsq_func(params):
1098 return (spc-self.double_gaussian_skew(freq,params[0],params[1],params[2],params[3],params[4],params[5],params[6],params[7]))/spcm
1099
1100 # fit
1101 bounds=([0,-numpy.inf,0,0,-400,0,0,0],[numpy.inf,-340,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1102
1103 params_scale = [spc_max,freq_max,freq_max,spc_max,freq_max,freq_max,1,spc_max]
1104 x0_value = numpy.array([spc_max,-400,30,spc_max/4,-200,150,1,1.0e7])
1105 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1106 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]
1107 A2f = popt.x[3]; B2f = popt.x[4]; C2f = popt.x[5]; K2f = popt.x[6]
1108 Df = popt.x[7]
1109
1110 aux = self.gaussian_skew(freq, A2f, B2f, C2f, K2f, Df)
1111 doppler = freq[numpy.argmax(aux)]
1112
1113 return A1f, B1f, C1f, A2f, B2f, C2f, K2f, Df, doppler
1114
1115 def Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(self,spc,freq,Nincoh,hei):
1116
1117 from scipy.optimize import least_squares
1118
1119 freq_max = numpy.max(numpy.abs(freq))
1120 spc_max = numpy.max(spc)
1121
1122 #from scipy.signal import medfilt
1123 #Nincoh = 20
1124 #Nincoh = 80
1125 Nincoh = Nincoh
1126 #spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1127 spcm = spc/numpy.sqrt(Nincoh)
1128
1129 # define a least squares function to optimize
1130 def lsq_func(params):
1131 return (spc-self.double_gaussian_double_skew(freq,params[0],params[1],params[2],params[3],params[4],params[5],params[6],params[7],params[8]))/spcm
1132
1133 # fit
1134 bounds=([0,-numpy.inf,0,-5,0,-400,0,0,0],[numpy.inf,-200,numpy.inf,5,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1135
1136 params_scale = [spc_max,freq_max,freq_max,1,spc_max,freq_max,freq_max,1,spc_max]
1137
1138 dop1_x0 = freq[numpy.argmax(spc)]
1139 if dop1_x0 < 0:
1140 dop2_x0 = dop1_x0 + 100
1141 if dop1_x0 > 0:
1142 dop2_x0 = dop1_x0 - 100
1143
1144 x0_value = numpy.array([spc_max,dop1_x0,30,-.1,spc_max/4, dop2_x0,150,1,1.0e7])
1145 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1146 J = popt.jac
1147
1148 try:
1149 cov = numpy.linalg.inv(J.T.dot(J))
1150 error = numpy.sqrt(numpy.diagonal(cov))
1151 except:
1152 error = numpy.ones((9))*numpy.NAN
1153
1154 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]; K1f = popt.x[3]
1155 A2f = popt.x[4]; B2f = popt.x[5]; C2f = popt.x[6]; K2f = popt.x[7]
1156 Df = popt.x[8]
1157 aux1 = self.gaussian_skew(freq, A1f, B1f, C1f, K1f, Df)
1158 doppler1 = freq[numpy.argmax(aux1)]
1159
1160 aux2 = self.gaussian_skew(freq, A2f, B2f, C2f, K2f, Df)
1161 doppler2 = freq[numpy.argmax(aux2)]
1162 #print("error",error)
1163 #exit(1)
1164
1165
1166 return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler1, doppler2, error
1167
1168 def Double_Gauss_fit_weight_bound_no_inputs(self,spc,freq,Nincoh):
1169
1170 from scipy.optimize import least_squares
1171
1172 freq_max = numpy.max(numpy.abs(freq))
1173 spc_max = numpy.max(spc)
1174
1175 from scipy.signal import medfilt
1176 Nincoh = 20
1177 Nincoh = 80
1178 Nincoh = Nincoh
1179 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1180
1181 # define a least squares function to optimize
1182 def lsq_func(params):
1183 return (spc-self.double_gaussian(freq,params[0],params[1],params[2],params[3],params[4],params[5],params[6]))/spcm
1184
1185 # fit
1186 # bounds=([0,-460,0,0,-400,120,0],[numpy.inf,-340,50,numpy.inf,0,250,numpy.inf])
1187 # bounds=([0,-numpy.inf,0,0,-numpy.inf,0,-numpy.inf,0],[numpy.inf,-200,numpy.inf,numpy.inf,0,numpy.inf,0,numpy.inf])
1188 #print(a1,b1,c1,a2,b2,c2,k2,d)
1189
1190 dop1_x0 = freq[numpy.argmax(spcm)]
1191
1192 bounds=([0,-numpy.inf,0,0,dop1_x0-50,0,0],[numpy.inf,-300,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf])
1193 params_scale = [spc_max,freq_max,freq_max,spc_max,freq_max,freq_max,spc_max]
1194 x0_value = numpy.array([spc_max,-400.5,30,spc_max/4,dop1_x0,150,1.0e7])
1195 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1196 J = popt.jac
1197
1198 try:
1199 cov = numpy.linalg.inv(J.T.dot(J))
1200 error = numpy.sqrt(numpy.diagonal(cov))
1201 except:
1202 error = numpy.ones((7))*numpy.NAN
1203
1204 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]
1205 A2f = popt.x[3]; B2f = popt.x[4]; C2f = popt.x[5]
1206 Df = popt.x[6]
1207 return A1f, B1f, C1f, A2f, B2f, C2f, Df, error
1208
1209 def Double_Gauss_Double_Skew_fit_weight_bound_with_inputs(self, spc, freq, a1, b1, c1, a2, b2, c2, k2, d):
1210
1211 from scipy.optimize import least_squares
1212
1213 freq_max = numpy.max(numpy.abs(freq))
1214 spc_max = numpy.max(spc)
1215
1216 from scipy.signal import medfilt
1217 Nincoh = dataOut.nIncohInt
1218 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1219
1220 # define a least squares function to optimize
1221 def lsq_func(params):
1222 return (spc-self.double_gaussian_double_skew(freq,params[0],params[1],params[2],params[3],params[4],params[5],params[6],params[7],params[8]))/spcm
1223
1224
1225 bounds=([0,-numpy.inf,0,-numpy.inf,0,-400,0,0,0],[numpy.inf,-340,numpy.inf,0,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf])
1226
1227 params_scale = [spc_max,freq_max,freq_max,1,spc_max,freq_max,freq_max,1,spc_max]
1228
1229 x0_value = numpy.array([a1,b1,c1,-.1,a2,b2,c2,k2,d])
1230
1231 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1232
1233 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]; K1f = popt.x[3]
1234 A2f = popt.x[4]; B2f = popt.x[5]; C2f = popt.x[6]; K2f = popt.x[7]
1235 Df = popt.x[8]
1236
1237 aux = self.gaussian_skew(freq, A2f, B2f, C2f, K2f, Df)
1238 doppler = x[numpy.argmax(aux)]
1239
1240 return A1f, B1f, C1f, K1f, A2f, B2f, C2f, K2f, Df, doppler
1241
1242 def Triple_Gauss_Skew_fit_weight_bound_no_inputs(self,spc,freq):
1243
1244 from scipy.optimize import least_squares
1245
1246 freq_max = numpy.max(numpy.abs(freq))
1247 spc_max = numpy.max(spc)
1248
1249 from scipy.signal import medfilt
1250 Nincoh = 20
1251 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1252
1253 # define a least squares function to optimize
1254 def lsq_func(params):
1255 return (spc-self.triple_gaussian_skew(freq,params[0],params[1],params[2],params[3],params[4],params[5],params[6],params[7],params[8],params[9],params[10],params[11]))/spcm
1256
1257 # fit
1258 bounds=([0,-numpy.inf,0,0,-400,0,0,0,0,0,0,0],[numpy.inf,-340,numpy.inf,numpy.inf,0,numpy.inf,numpy.inf,numpy.inf,numpy.inf,numpy.inf,numpy.inf,numpy.inf])
1259
1260 params_scale = [spc_max,freq_max,freq_max,spc_max,freq_max,freq_max,1,spc_max,freq_max,freq_max,1,spc_max]
1261 x0_value = numpy.array([spc_max,-400,30,spc_max/4,-200,150,1,spc_max/4,400,150,1,1.0e7])
1262 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1263
1264 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]
1265 A2f = popt.x[3]; B2f = popt.x[4]; C2f = popt.x[5]; K2f = popt.x[6]
1266 A3f = popt.x[7]; B3f = popt.x[8]; C3f = popt.x[9]; K3f = popt.x[10]
1267 Df = popt.x[11]
1268
1269 aux = self.gaussian_skew(freq, A2f, B2f, C2f, K2f, Df)
1270 doppler = freq[numpy.argmax(aux)]
1271
1272 return A1f, B1f, C1f, A2f, B2f, C2f, K2f, A3f, B3f, C3f, K3f, Df, doppler
1273
1274 def CEEJ_Skew_fit_weight_bound_no_inputs(self,spc,freq,Nincoh):
1275
1276 from scipy.optimize import least_squares
1277
1278 freq_max = numpy.max(numpy.abs(freq))
1279 spc_max = numpy.max(spc)
1280
1281 from scipy.signal import medfilt
1282 Nincoh = 20
1283 Nincoh = 80
1284 Nincoh = Nincoh
1285 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1286
1287 # define a least squares function to optimize
1288 def lsq_func(params):
1289 return (spc-self.gaussian_skew(freq,params[0],params[1],params[2],params[3],params[4]))#/spcm
1290
1291
1292 bounds=([0,0,0,-numpy.inf,0],[numpy.inf,numpy.inf,numpy.inf,0,numpy.inf])
1293
1294 params_scale = [spc_max,freq_max,freq_max,1,spc_max]
1295
1296 x0_value = numpy.array([spc_max,freq[numpy.argmax(spc)],30,-.1,numpy.mean(spc[:50])])
1297
1298 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1299
1300 J = popt.jac
1301
1302 try:
1303 error = numpy.ones((9))*numpy.NAN
1304 cov = numpy.linalg.inv(J.T.dot(J))
1305 error[:4] = numpy.sqrt(numpy.diagonal(cov))[:4]
1306 error[-1] = numpy.sqrt(numpy.diagonal(cov))[-1]
1307 except:
1308 error = numpy.ones((9))*numpy.NAN
1309
1310 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]; K1f = popt.x[3]
1311 Df = popt.x[4]
1312
1313 aux1 = self.gaussian_skew(freq, A1f, B1f, C1f, K1f, Df)
1314 doppler1 = freq[numpy.argmax(aux1)]
1315 #print("CEEJ ERROR:",error)
1316
1317 return A1f, B1f, C1f, K1f, numpy.NAN, numpy.NAN, numpy.NAN, numpy.NAN, Df, doppler1, numpy.NAN, error
1318
1319 def CEEJ_fit_weight_bound_no_inputs(self,spc,freq,Nincoh):
1320
1321 from scipy.optimize import least_squares
1322
1323 freq_max = numpy.max(numpy.abs(freq))
1324 spc_max = numpy.max(spc)
1325
1326 from scipy.signal import medfilt
1327 Nincoh = 20
1328 Nincoh = 80
1329 Nincoh = Nincoh
1330 spcm = medfilt(spc,11)/numpy.sqrt(Nincoh)
1331
1332 # define a least squares function to optimize
1333 def lsq_func(params):
1334 return (spc-self.gaussian(freq,params[0],params[1],params[2],params[3]))#/spcm
1335
1336
1337 bounds=([0,0,0,0],[numpy.inf,numpy.inf,numpy.inf,numpy.inf])
1338
1339 params_scale = [spc_max,freq_max,freq_max,spc_max]
1340
1341 x0_value = numpy.array([spc_max,freq[numpy.argmax(spcm)],30,numpy.mean(spc[:50])])
1342
1343 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1344
1345 J = popt.jac
1346
1347 try:
1348 error = numpy.ones((4))*numpy.NAN
1349 cov = numpy.linalg.inv(J.T.dot(J))
1350 error = numpy.sqrt(numpy.diagonal(cov))
1351 except:
1352 error = numpy.ones((4))*numpy.NAN
1353
1354 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]
1355 Df = popt.x[3]
1356
1357 return A1f, B1f, C1f, Df, error
1358
1359 def Simple_fit_bound(self,spc,freq,Nincoh):
1360
1361 freq_max = numpy.max(numpy.abs(freq))
1362 spc_max = numpy.max(spc)
1363
1364 Nincoh = Nincoh
1365
1366 def lsq_func(params):
1367 return (spc-self.gaussian(freq,params[0],params[1],params[2],params[3]))
1368
1369 bounds=([0,-50,0,0],[numpy.inf,+50,numpy.inf,numpy.inf])
1370
1371 params_scale = [spc_max,freq_max,freq_max,spc_max]
1372
1373 x0_value = numpy.array([spc_max,-20.5,5,1.0e7])
1374
1375 popt = least_squares(lsq_func,x0=x0_value,x_scale=params_scale,bounds=bounds,verbose=0)
1376
1377 J = popt.jac
1378
1379 try:
1380 cov = numpy.linalg.inv(J.T.dot(J))
1381 error = numpy.sqrt(numpy.diagonal(cov))
1382 except:
1383 error = numpy.ones((4))*numpy.NAN
1384
1385 A1f = popt.x[0]; B1f = popt.x[1]; C1f = popt.x[2]
1386 Df = popt.x[3]
1387
1388 return A1f, B1f, C1f, Df, error
1389
1390 def clean_outliers(self,param):
1391
1392 threshold = 700
1393
1394 param = numpy.where(param < -threshold, numpy.nan, param)
1395 param = numpy.where(param > +threshold, numpy.nan, param)
1396
1397 return param
1398
1399 def windowing_single(self,spc,x,A,B,C,D,nFFTPoints):
1400 from scipy.optimize import curve_fit,fmin
1401
1402 def R_gaussian(x, a, b, c):
1403 N = int(numpy.shape(x)[0])
1404 val = a * numpy.exp(-((x)*c*2*2*numpy.pi)**2 / (2))* numpy.exp(1.j*b*x*4*numpy.pi)
1405 return val
1406
1407 def T(x,N):
1408 T = 1-abs(x)/N
1409 return T
1410
1411 def R_T_spc_fun(x, a, b, c, d, nFFTPoints):
1412
1413 N = int(numpy.shape(x)[0])
1414
1415 x_max = x[-1]
1416
1417 x_pos = x[int(nFFTPoints/2):]
1418 x_neg = x[:int(nFFTPoints/2)]
1419
1420 R_T_neg_1 = R_gaussian(x, a, b, c)[:int(nFFTPoints/2)]*T(x_neg,-x[0])
1421 R_T_pos_1 = R_gaussian(x, a, b, c)[int(nFFTPoints/2):]*T(x_pos,x[-1])
1422 R_T_sum_1 = R_T_pos_1 + R_T_neg_1
1423 R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real
1424 R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1)
1425 max_val_1 = numpy.max(R_T_spc_1)
1426 R_T_spc_1 = R_T_spc_1*a/max_val_1
1427
1428 R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N))
1429 R_T_d_neg = R_T_d[:int(nFFTPoints/2)]*T(x_neg,-x[0])
1430 R_T_d_pos = R_T_d[int(nFFTPoints/2):]*T(x_pos,x[-1])
1431 R_T_d_sum = R_T_d_pos + R_T_d_neg
1432 R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real
1433 R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3)
1434
1435 R_T_final = R_T_spc_1 + R_T_spc_3
1436
1437 return R_T_final
1438
1439 y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))
1440
1441 from scipy.stats import norm
1442 mean,std=norm.fit(spc)
1443
1444 # estimate starting values from the data
1445 a = A
1446 b = B
1447 c = C#numpy.std(spc)
1448 d = D
1449 '''
1450 ippSeconds = 250*20*1.e-6/3
1451
1452 x_t = ippSeconds * (numpy.arange(1600) -1600 / 2.)
1453
1454 x_t = numpy.linspace(x_t[0],x_t[-1],3200)
1455
1456 x_freq = numpy.fft.fftfreq(1600,d=ippSeconds)
1457 x_freq = numpy.fft.fftshift(x_freq)
1458 '''
1459 # define a least squares function to optimize
1460 def minfunc(params):
1461 return sum((y-R_T_spc_fun(x,params[0],params[1],params[2],params[3],params[4],params[5],params[6]))**2/1)#y**2)
1462
1463 # fit
1464
1465 popt_full = fmin(minfunc,[a,b,c,d],full_output=True)
1466 #print("nIter", popt_full[2])
1467 popt = popt_full[0]
1468
1469 #return R_T_spc_fun(x_t,popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]), popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]
1470 return popt[0], popt[1], popt[2], popt[3]
1471
1472 def run(self, dataOut, mode = 0, Hmin1 = None, Hmax1 = None, Hmin2 = None, Hmax2 = None, Dop = 'Shift'):
818
1473
819 pwcode = 1
1474 pwcode = 1
820
1475
@@ -828,39 +1483,381 class Oblique_Gauss_Fit(Operation):
828 dataOut.power = numpy.average(z, axis=1)
1483 dataOut.power = numpy.average(z, axis=1)
829 dataOut.powerdB = 10 * numpy.log10(dataOut.power)
1484 dataOut.powerdB = 10 * numpy.log10(dataOut.power)
830
1485
831
832 x = dataOut.getVelRange(0)
1486 x = dataOut.getVelRange(0)
833
1487
834 dataOut.Oblique_params = numpy.ones((1,7,dataOut.nHeights))*numpy.NAN
1488 dataOut.Oblique_params = numpy.ones((1,7,dataOut.nHeights))*numpy.NAN
835 dataOut.Oblique_param_errors = numpy.ones((1,7,dataOut.nHeights))*numpy.NAN
1489 dataOut.Oblique_param_errors = numpy.ones((1,7,dataOut.nHeights))*numpy.NAN
1490 dataOut.dplr_2_u = numpy.ones((1,1,dataOut.nHeights))*numpy.NAN
1491
1492 if mode == 6:
1493 dataOut.Oblique_params = numpy.ones((1,9,dataOut.nHeights))*numpy.NAN
1494 elif mode == 7:
1495 dataOut.Oblique_params = numpy.ones((1,13,dataOut.nHeights))*numpy.NAN
1496 elif mode == 8:
1497 dataOut.Oblique_params = numpy.ones((1,10,dataOut.nHeights))*numpy.NAN
1498 elif mode == 9:
1499 dataOut.Oblique_params = numpy.ones((1,11,dataOut.nHeights))*numpy.NAN
1500 dataOut.Oblique_param_errors = numpy.ones((1,9,dataOut.nHeights))*numpy.NAN
1501 elif mode == 11:
1502 dataOut.Oblique_params = numpy.ones((1,7,dataOut.nHeights))*numpy.NAN
1503 dataOut.Oblique_param_errors = numpy.ones((1,7,dataOut.nHeights))*numpy.NAN
1504 elif mode == 10: #150 km
1505 dataOut.Oblique_params = numpy.ones((1,4,dataOut.nHeights))*numpy.NAN
1506 dataOut.Oblique_param_errors = numpy.ones((1,4,dataOut.nHeights))*numpy.NAN
1507 dataOut.snr_log10 = numpy.ones((1,dataOut.nHeights))*numpy.NAN
836
1508
837 dataOut.VelRange = x
1509 dataOut.VelRange = x
838
1510
839
1511
840 l1=range(22,36)
841 l2=range(58,99)
842
1512
843 for hei in itertools.chain(l1, l2):
1513 #l1=range(22,36) #+62
1514 #l1=range(32,36)
1515 #l2=range(58,99) #+62
844
1516
845 try:
1517 #if Hmin1 == None or Hmax1 == None or Hmin2 == None or Hmax2 == None:
846 spc = dataOut.data_spc[0,:,hei]
847
1518
848 spc_fit, A1, B1, C1, D1 = self.Gauss_fit_2(spc,x,'first')
1519 minHei1 = 105.
1520 maxHei1 = 122.5
1521 maxHei1 = 130.5
849
1522
850 spc_diff = spc - spc_fit
1523 if mode == 10: #150 km
851 spc_diff[spc_diff < 0] = 0
1524 minHei1 = 100
1525 maxHei1 = 100
852
1526
853 spc_fit_diff, A2, B2, C2, D2 = self.Gauss_fit_2(spc_diff,x,'second')
1527 inda1 = numpy.where(dataOut.heightList >= minHei1)
1528 indb1 = numpy.where(dataOut.heightList <= maxHei1)
854
1529
855 D = (D1+D2)
1530 minIndex1 = inda1[0][0]
1531 maxIndex1 = indb1[0][-1]
1532
1533 minHei2 = 150.
1534 maxHei2 = 201.25
1535 maxHei2 = 225.3
1536
1537 if mode == 10: #150 km
1538 minHei2 = 110
1539 maxHei2 = 165
1540
1541 inda2 = numpy.where(dataOut.heightList >= minHei2)
1542 indb2 = numpy.where(dataOut.heightList <= maxHei2)
1543
1544 minIndex2 = inda2[0][0]
1545 maxIndex2 = indb2[0][-1]
1546
1547 l1=range(minIndex1,maxIndex1)
1548 l2=range(minIndex2,maxIndex2)
1549
1550 if mode == 4:
1551 '''
1552 for ind in range(dataOut.nHeights):
1553 if(dataOut.heightList[ind]>=168 and dataOut.heightList[ind]<188):
1554 try:
1555 dataOut.Oblique_params[0,0,ind],dataOut.Oblique_params[0,1,ind],dataOut.Oblique_params[0,2,ind],dataOut.Oblique_params[0,3,ind],dataOut.Oblique_params[0,4,ind],dataOut.Oblique_params[0,5,ind],dataOut.Oblique_params[0,6,ind],dataOut.Oblique_param_errors[0,0,ind],dataOut.Oblique_param_errors[0,1,ind],dataOut.Oblique_param_errors[0,2,ind],dataOut.Oblique_param_errors[0,3,ind],dataOut.Oblique_param_errors[0,4,ind],dataOut.Oblique_param_errors[0,5,ind],dataOut.Oblique_param_errors[0,6,ind] = self.DH_mode(dataOut.data_spc[0,:,ind],dataOut.VelRange)
1556 except:
1557 pass
1558 '''
1559 for ind in itertools.chain(l1, l2):
1560
1561 try:
1562 dataOut.Oblique_params[0,0,ind],dataOut.Oblique_params[0,1,ind],dataOut.Oblique_params[0,2,ind],dataOut.Oblique_params[0,3,ind],dataOut.Oblique_params[0,4,ind],dataOut.Oblique_params[0,5,ind],dataOut.Oblique_params[0,6,ind],dataOut.Oblique_param_errors[0,0,ind],dataOut.Oblique_param_errors[0,1,ind],dataOut.Oblique_param_errors[0,2,ind],dataOut.Oblique_param_errors[0,3,ind],dataOut.Oblique_param_errors[0,4,ind],dataOut.Oblique_param_errors[0,5,ind],dataOut.Oblique_param_errors[0,6,ind] = self.DH_mode(dataOut.data_spc[0,:,ind],dataOut.VelRange)
1563 dataOut.dplr_2_u[0,0,ind] = dataOut.Oblique_params[0,4,ind]/numpy.sin(numpy.arccos(102/dataOut.heightList[ind]))
1564 except:
1565 pass
1566
1567 else:
1568 for hei in itertools.chain(l1, l2):
1569 if numpy.isnan(dataOut.snl[0,hei]) or dataOut.snl[0,hei]<.0:
1570
1571 continue #Avoids the analysis when there is only noise
1572
1573 try:
1574 spc = dataOut.data_spc[0,:,hei]
1575
1576 if mode == 6: #Skew Weighted Bounded
1577 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_params[0,7,hei],dataOut.Oblique_params[0,8,hei] = self.Double_Gauss_Skew_fit_weight_bound_no_inputs(spc,x)
1578 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,8,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1579
1580 elif mode == 7: #Triple Skew Weighted Bounded
1581 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_params[0,7,hei],dataOut.Oblique_params[0,8,hei],dataOut.Oblique_params[0,9,hei],dataOut.Oblique_params[0,10,hei],dataOut.Oblique_params[0,11,hei],dataOut.Oblique_params[0,12,hei] = self.Triple_Gauss_Skew_fit_weight_bound_no_inputs(spc,x)
1582 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,12,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1583
1584 elif mode == 8: #Double Skewed Weighted Bounded with inputs
1585 a1, b1, c1, a2, b2, c2, k2, d, dopp = self.Double_Gauss_Skew_fit_weight_bound_no_inputs(spc,x)
1586 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_params[0,7,hei],dataOut.Oblique_params[0,8,hei],dataOut.Oblique_params[0,9,hei] = self.Double_Gauss_Skew_fit_weight_bound_no_inputs(spc,x, a1, b1, c1, a2, b2, c2, k2, d)
1587 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,9,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1588
1589 elif mode == 9: #Double Skewed Weighted Bounded no inputs
1590 #if numpy.max(spc) <= 0:
1591 from scipy.signal import medfilt
1592 spcm = medfilt(spc,11)
1593 if x[numpy.argmax(spcm)] <= 0:
1594 #print("EEJ", dataOut.heightList[hei], hei)
1595 #if hei != 70:
1596 #continue
1597 #else:
1598 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_params[0,7,hei],dataOut.Oblique_params[0,8,hei],dataOut.Oblique_params[0,9,hei],dataOut.Oblique_params[0,10,hei],dataOut.Oblique_param_errors[0,:,hei] = self.Double_Gauss_Double_Skew_fit_weight_bound_no_inputs(spcm,x,dataOut.nIncohInt,dataOut.heightList[hei])
1599 #if dataOut.Oblique_params[0,-2,hei] < -500 or dataOut.Oblique_params[0,-2,hei] > 500 or dataOut.Oblique_params[0,-1,hei] < -500 or dataOut.Oblique_params[0,-1,hei] > 500:
1600 # dataOut.Oblique_params[0,:,hei] *= numpy.NAN
1601 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1602
1603 else:
1604 #print("CEEJ")
1605 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_params[0,7,hei],dataOut.Oblique_params[0,8,hei],dataOut.Oblique_params[0,9,hei],dataOut.Oblique_params[0,10,hei],dataOut.Oblique_param_errors[0,:,hei] = self.CEEJ_Skew_fit_weight_bound_no_inputs(spcm,x,dataOut.nIncohInt)
1606 #if dataOut.Oblique_params[0,-2,hei] < -500 or dataOut.Oblique_params[0,-2,hei] > 500 or dataOut.Oblique_params[0,-1,hei] < -500 or dataOut.Oblique_params[0,-1,hei] > 500:
1607 # dataOut.Oblique_params[0,:,hei] *= numpy.NAN
1608 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,10,hei]/numpy.sin(numpy.arccos(100./dataOut.heightList[hei]))
1609 elif mode == 11: #Double Weighted Bounded no inputs
1610 #if numpy.max(spc) <= 0:
1611 from scipy.signal import medfilt
1612 spcm = medfilt(spc,11)
1613
1614 if x[numpy.argmax(spcm)] <= 0:
1615 #print("EEJ")
1616 #print("EEJ",dataOut.heightList[hei])
1617 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_param_errors[0,:,hei] = self.Double_Gauss_fit_weight_bound_no_inputs(spc,x,dataOut.nIncohInt)
1618 #if dataOut.Oblique_params[0,-2,hei] < -500 or dataOut.Oblique_params[0,-2,hei] > 500 or dataOut.Oblique_params[0,-1,hei] < -500 or dataOut.Oblique_params[0,-1,hei] > 500:
1619 # dataOut.Oblique_params[0,:,hei] *= numpy.NAN
1620 else:
1621 #print("CEEJ",dataOut.heightList[hei])
1622 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_param_errors[0,:,hei] = self.CEEJ_fit_weight_bound_no_inputs(spc,x,dataOut.nIncohInt)
1623
1624 elif mode == 10: #150km
1625 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_param_errors[0,:,hei] = self.Simple_fit_bound(spc,x,dataOut.nIncohInt)
1626 snr = (dataOut.power[0,hei]*factor - dataOut.Oblique_params[0,3,hei])/dataOut.Oblique_params[0,3,hei]
1627 dataOut.snr_log10[0,hei] = numpy.log10(snr)
1628
1629 else:
1630 spc_fit, A1, B1, C1, D1 = self.Gauss_fit_2(spc,x,'first')
1631
1632 spc_diff = spc - spc_fit
1633 spc_diff[spc_diff < 0] = 0
1634
1635 spc_fit_diff, A2, B2, C2, D2 = self.Gauss_fit_2(spc_diff,x,'second')
1636
1637 D = (D1+D2)
1638
1639 if mode == 0: #Double Fit
1640 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_param_errors[0,0,hei],dataOut.Oblique_param_errors[0,1,hei],dataOut.Oblique_param_errors[0,2,hei],dataOut.Oblique_param_errors[0,3,hei],dataOut.Oblique_param_errors[0,4,hei],dataOut.Oblique_param_errors[0,5,hei],dataOut.Oblique_param_errors[0,6,hei] = self.Double_Gauss_fit_2(spc,x,A1,B1,C1,A2,B2,C2,D)
1641 #spc_double_fit,dataOut.Oblique_params = self.Double_Gauss_fit(spc,x,A1,B1,C1,A2,B2,C2,D)
1642
1643 elif mode == 1: #Double Fit Windowed
1644 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei] = self.windowing_double(spc,dataOut.getFreqRange(0),A1,B1,C1,A2,B2,C2,D)
1645
1646 elif mode == 2: #Double Fit Weight
1647 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei] = self.Double_Gauss_fit_weight(spc,x,A1,B1,C1,A2,B2,C2,D)
1648
1649 elif mode == 3: #Simple Fit
1650 dataOut.Oblique_params[0,0,hei] = A1
1651 dataOut.Oblique_params[0,1,hei] = B1
1652 dataOut.Oblique_params[0,2,hei] = C1
1653 dataOut.Oblique_params[0,3,hei] = A2
1654 dataOut.Oblique_params[0,4,hei] = B2
1655 dataOut.Oblique_params[0,5,hei] = C2
1656 dataOut.Oblique_params[0,6,hei] = D
1657
1658 elif mode == 5: #Triple Fit Weight
1659 if hei in l1:
1660 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei] = self.duo_Marco(spc,x,A1,B1,C1,A2,B2,C2,D)
1661 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,4,hei]/numpy.sin(numpy.arccos(102/dataOut.heightList[hei]))
1662 #print(dataOut.Oblique_params[0,0,hei])
1663 #print(dataOut.dplr_2_u[0,0,hei])
1664 else:
1665 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei] = self.Double_Gauss_fit_weight(spc,x,A1,B1,C1,A2,B2,C2,D)
1666 dataOut.dplr_2_u[0,0,hei] = dataOut.Oblique_params[0,4,hei]/numpy.sin(numpy.arccos(102/dataOut.heightList[hei]))
1667
1668
1669 except:
1670 ###dataOut.Oblique_params[0,:,hei] = dataOut.Oblique_params[0,:,hei]*numpy.NAN
1671 pass
1672
1673 #exit(1)
1674 dataOut.paramInterval = dataOut.nProfiles*dataOut.nCohInt*dataOut.ippSeconds
1675 dataOut.lat=-11.95
1676 dataOut.lon=-76.87
1677 '''
1678 dataOut.Oblique_params = numpy.where(dataOut.Oblique_params<-700, numpy.nan, dop_t1)
1679 dataOut.Oblique_params = numpy.where(dataOut.Oblique_params<+700, numpy.nan, dop_t1)
1680 Aquí debo exceptuar las amplitudes
1681 '''
1682 if mode == 9: #Double Skew Gaussian
1683 #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,-2,:] #Pos[Max_value]
1684 #dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,1,:] #Shift
1685 dataOut.Spec_W_T1 = dataOut.Oblique_params[:,2,:]
1686 #dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,-1,:] #Pos[Max_value]
1687 #dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,5,:] #Shift
1688 dataOut.Spec_W_T2 = dataOut.Oblique_params[:,6,:]
1689 if Dop == 'Shift':
1690 dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,1,:] #Shift
1691 dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,5,:] #Shift
1692 elif Dop == 'Max':
1693 dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,-2,:] #Pos[Max_value]
1694 dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,-1,:] #Pos[Max_value]
1695
1696 dataOut.Err_Dop_EEJ_T1 = dataOut.Oblique_param_errors[:,1,:] #En realidad este es el error?
1697 dataOut.Err_Spec_W_T1 = dataOut.Oblique_param_errors[:,2,:]
1698 dataOut.Err_Dop_EEJ_T2 = dataOut.Oblique_param_errors[:,5,:] #En realidad este es el error?
1699 dataOut.Err_Spec_W_T2 = dataOut.Oblique_param_errors[:,6,:]
1700
1701 elif mode == 11: #Double Gaussian
1702 dataOut.Dop_EEJ_T1 = dataOut.Oblique_params[:,1,:]
1703 dataOut.Spec_W_T1 = dataOut.Oblique_params[:,2,:]
1704 dataOut.Dop_EEJ_T2 = dataOut.Oblique_params[:,4,:]
1705 dataOut.Spec_W_T2 = dataOut.Oblique_params[:,5,:]
1706
1707 dataOut.Err_Dop_EEJ_T1 = dataOut.Oblique_param_errors[:,1,:]
1708 dataOut.Err_Spec_W_T1 = dataOut.Oblique_param_errors[:,2,:]
1709 dataOut.Err_Dop_EEJ_T2 = dataOut.Oblique_param_errors[:,4,:]
1710 dataOut.Err_Spec_W_T2 = dataOut.Oblique_param_errors[:,5,:]
1711
1712 #print("Before: ", dataOut.Dop_EEJ_T2)
1713 dataOut.Spec_W_T1 = self.clean_outliers(dataOut.Spec_W_T1)
1714 dataOut.Spec_W_T2 = self.clean_outliers(dataOut.Spec_W_T2)
1715 dataOut.Dop_EEJ_T1 = self.clean_outliers(dataOut.Dop_EEJ_T1)
1716 dataOut.Dop_EEJ_T2 = self.clean_outliers(dataOut.Dop_EEJ_T2)
1717 #print("After: ", dataOut.Dop_EEJ_T2)
1718 dataOut.Err_Spec_W_T1 = self.clean_outliers(dataOut.Err_Spec_W_T1)
1719 dataOut.Err_Spec_W_T2 = self.clean_outliers(dataOut.Err_Spec_W_T2)
1720 dataOut.Err_Dop_EEJ_T1 = self.clean_outliers(dataOut.Err_Dop_EEJ_T1)
1721 dataOut.Err_Dop_EEJ_T2 = self.clean_outliers(dataOut.Err_Dop_EEJ_T2)
1722 #print("Before data_snr: ", dataOut.data_snr)
1723 #dataOut.data_snr = numpy.where(numpy.isnan(dataOut.Dop_EEJ_T1), numpy.nan, dataOut.data_snr)
1724 dataOut.snl = numpy.where(numpy.isnan(dataOut.Dop_EEJ_T1), numpy.nan, dataOut.snl)
1725
1726 #print("After data_snr: ", dataOut.data_snr)
1727 dataOut.mode = mode
1728 dataOut.flagNoData = numpy.all(numpy.isnan(dataOut.Dop_EEJ_T1)) #Si todos los valores son NaN no se prosigue
1729 ###dataOut.flagNoData = False #Descomentar solo para ploteo sino mantener comentado (para guardado)
1730
1731 return dataOut
1732
1733 class Gaussian_Windowed(Operation):
1734 '''
1735 Written by R. Flores
1736 '''
1737 def __init__(self):
1738 Operation.__init__(self)
1739
1740 def windowing_single(self,spc,x,A,B,C,D,nFFTPoints):
1741 from scipy.optimize import curve_fit,fmin
1742
1743 def gaussian(x, a, b, c, d):
1744 val = a * numpy.exp(-(x - b)**2 / (2*c**2)) + d
1745 return val
1746
1747 def R_gaussian(x, a, b, c):
1748 N = int(numpy.shape(x)[0])
1749 val = a * numpy.exp(-((x)*c*2*2*numpy.pi)**2 / (2))* numpy.exp(1.j*b*x*4*numpy.pi)
1750 return val
1751
1752 def T(x,N):
1753 T = 1-abs(x)/N
1754 return T
1755
1756 def R_T_spc_fun(x, a, b, c, d, nFFTPoints):
1757
1758 N = int(numpy.shape(x)[0])
1759
1760 x_max = x[-1]
1761
1762 x_pos = x[nFFTPoints:]
1763 x_neg = x[:nFFTPoints]
1764 #print([int(nFFTPoints/2))
1765 #print("x: ", x)
1766 #print("x_neg: ", x_neg)
1767 #print("x_pos: ", x_pos)
1768
1769
1770 R_T_neg_1 = R_gaussian(x, a, b, c)[:nFFTPoints]*T(x_neg,-x[0])
1771 R_T_pos_1 = R_gaussian(x, a, b, c)[nFFTPoints:]*T(x_pos,x[-1])
1772 #print(T(x_pos,x[-1]),x_pos,x[-1])
1773 #print(R_T_neg_1.shape,R_T_pos_1.shape)
1774 R_T_sum_1 = R_T_pos_1 + R_T_neg_1
1775 R_T_spc_1 = numpy.fft.fft(R_T_sum_1).real
1776 R_T_spc_1 = numpy.fft.fftshift(R_T_spc_1)
1777 max_val_1 = numpy.max(R_T_spc_1)
1778 R_T_spc_1 = R_T_spc_1*a/max_val_1
1779
1780 R_T_d = d*numpy.fft.fftshift(signal.unit_impulse(N))
1781 R_T_d_neg = R_T_d[:nFFTPoints]*T(x_neg,-x[0])
1782 R_T_d_pos = R_T_d[nFFTPoints:]*T(x_pos,x[-1])
1783 R_T_d_sum = R_T_d_pos + R_T_d_neg
1784 R_T_spc_3 = numpy.fft.fft(R_T_d_sum).real
1785 R_T_spc_3 = numpy.fft.fftshift(R_T_spc_3)
1786
1787 R_T_final = R_T_spc_1 + R_T_spc_3
1788
1789 return R_T_final
1790
1791 y = spc#gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))
1792
1793 from scipy.stats import norm
1794 mean,std=norm.fit(spc)
1795
1796 # estimate starting values from the data
1797 a = A
1798 b = B
1799 c = C#numpy.std(spc)
1800 d = D
1801 #'''
1802 #ippSeconds = 250*20*1.e-6/3
1803
1804 #x_t = ippSeconds * (numpy.arange(nFFTPoints) - nFFTPoints / 2.)
1805
1806 #x_t = numpy.linspace(x_t[0],x_t[-1],3200)
1807 #print("x_t: ", x_t)
1808 #print("nFFTPoints: ", nFFTPoints)
1809 x_vel = numpy.linspace(x[0],x[-1],int(2*nFFTPoints))
1810 #print("x_vel: ", x_vel)
1811 #x_freq = numpy.fft.fftfreq(1600,d=ippSeconds)
1812 #x_freq = numpy.fft.fftshift(x_freq)
1813 #'''
1814 # define a least squares function to optimize
1815 def minfunc(params):
1816 #print("y.shape: ", numpy.shape(y))
1817 return sum((y-R_T_spc_fun(x_vel,params[0],params[1],params[2],params[3],nFFTPoints))**2/1)#y**2)
1818
1819 # fit
1820
1821 popt_full = fmin(minfunc,[a,b,c,d], disp=False)
1822 #print("nIter", popt_full[2])
1823 popt = popt_full#[0]
1824
1825 fun = gaussian(x, popt[0], popt[1], popt[2], popt[3])
1826
1827 #return R_T_spc_fun(x_t,popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]), popt[0], popt[1], popt[2], popt[3], popt[4], popt[5], popt[6]
1828 return fun, popt[0], popt[1], popt[2], popt[3]
1829
1830 def run(self, dataOut):
1831
1832 from scipy.signal import medfilt
1833 import matplotlib.pyplot as plt
1834 dataOut.moments = numpy.ones((dataOut.nChannels,4,dataOut.nHeights))*numpy.NAN
1835 dataOut.VelRange = dataOut.getVelRange(0)
1836 for nChannel in range(dataOut.nChannels):
1837 for hei in range(dataOut.heightList.shape[0]):
1838 #print("ipp: ", dataOut.ippSeconds)
1839 spc = numpy.copy(dataOut.data_spc[nChannel,:,hei])
1840
1841 #print(VelRange)
1842 #print(dataOut.getFreqRange(64))
1843 spcm = medfilt(spc,11)
1844 spc_max = numpy.max(spcm)
1845 dop1_x0 = dataOut.VelRange[numpy.argmax(spcm)]
1846 D = numpy.min(spcm)
1847
1848 fun, A, B, C, D = self.windowing_single(spc,dataOut.VelRange,spc_max,dop1_x0,abs(dop1_x0),D,dataOut.nFFTPoints)
1849 dataOut.moments[nChannel,0,hei] = A
1850 dataOut.moments[nChannel,1,hei] = B
1851 dataOut.moments[nChannel,2,hei] = C
1852 dataOut.moments[nChannel,3,hei] = D
1853 '''
1854 plt.figure()
1855 plt.plot(VelRange,spc,marker='*',linestyle='')
1856 plt.plot(VelRange,fun)
1857 plt.title(dataOut.heightList[hei])
1858 plt.show()
1859 '''
856
1860
857 dataOut.Oblique_params[0,0,hei],dataOut.Oblique_params[0,1,hei],dataOut.Oblique_params[0,2,hei],dataOut.Oblique_params[0,3,hei],dataOut.Oblique_params[0,4,hei],dataOut.Oblique_params[0,5,hei],dataOut.Oblique_params[0,6,hei],dataOut.Oblique_param_errors[0,0,hei],dataOut.Oblique_param_errors[0,1,hei],dataOut.Oblique_param_errors[0,2,hei],dataOut.Oblique_param_errors[0,3,hei],dataOut.Oblique_param_errors[0,4,hei],dataOut.Oblique_param_errors[0,5,hei],dataOut.Oblique_param_errors[0,6,hei] = self.Double_Gauss_fit_2(spc,x,A1,B1,C1,A2,B2,C2,D)
858 #spc_double_fit,dataOut.Oblique_params = self.Double_Gauss_fit(spc,x,A1,B1,C1,A2,B2,C2,D)
859
860 except:
861 ###dataOut.Oblique_params[0,:,hei] = dataOut.Oblique_params[0,:,hei]*numpy.NAN
862 pass
863
864 return dataOut
1861 return dataOut
865
1862
866 class PrecipitationProc(Operation):
1863 class PrecipitationProc(Operation):
@@ -1223,7 +2220,7 class FullSpectralAnalysis(Operation):
1223 # spwd limit - updated by D. Scipión 30.03.2021
2220 # spwd limit - updated by D. Scipión 30.03.2021
1224 widthlimit = 10
2221 widthlimit = 10
1225 '''************************* SPC is normalized ********************************'''
2222 '''************************* SPC is normalized ********************************'''
1226 spc_norm = spc.copy()
2223 spc_norm = spc.copy()
1227 # For each channel
2224 # For each channel
1228 for i in range(nChan):
2225 for i in range(nChan):
1229 spc_sub = spc_norm[i,:] - noise[i] # only the signal power
2226 spc_sub = spc_norm[i,:] - noise[i] # only the signal power
@@ -1242,9 +2239,9 class FullSpectralAnalysis(Operation):
1242 >= 0, as it is the modulus squared of the signals (complex * it's conjugate)
2239 >= 0, as it is the modulus squared of the signals (complex * it's conjugate)
1243 """
2240 """
1244 # initial conditions
2241 # initial conditions
1245 popt = [1e-10,0,1e-10]
2242 popt = [1e-10,0,1e-10]
1246 # Spectra average
2243 # Spectra average
1247 SPCMean = numpy.average(SPC_Samples,0)
2244 SPCMean = numpy.average(SPC_Samples,0)
1248 # Moments in frequency
2245 # Moments in frequency
1249 SPCMoments = self.Moments(SPCMean[xvalid], xSamples_zoom)
2246 SPCMoments = self.Moments(SPCMean[xvalid], xSamples_zoom)
1250
2247
@@ -1621,6 +2618,8 class SpectralMoments(Operation):
1621 if (nicoh is None): nicoh = 1
2618 if (nicoh is None): nicoh = 1
1622 if (graph is None): graph = 0
2619 if (graph is None): graph = 0
1623 if (smooth is None): smooth = 0
2620 if (smooth is None): smooth = 0
2621 elif (self.smooth < 3): smooth = 0
2622
1624 if (type1 is None): type1 = 0
2623 if (type1 is None): type1 = 0
1625 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
2624 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
1626 if (snrth is None): snrth = -20.0
2625 if (snrth is None): snrth = -20.0
General Comments 0
You need to be logged in to leave comments. Login now