##// 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, (1027 lines changed) Show them Hide them
@@ -23,7 +23,6 import warnings
23 23 from numpy import NaN
24 24 from scipy.optimize.optimize import OptimizeWarning
25 25 warnings.filterwarnings('ignore')
26 import pdb
27 26
28 27
29 28 SPEED_OF_LIGHT = 299792458
@@ -663,7 +662,9 class GaussianFit(Operation):
663 662 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
664 663
665 664 class Oblique_Gauss_Fit(Operation):
666
665 '''
666 Written by R. Flores
667 '''
667 668 def __init__(self):
668 669 Operation.__init__(self)
669 670
@@ -714,7 +715,6 class Oblique_Gauss_Fit(Operation):
714 715
715 716 return gaussian(x, popt[0], popt[1], popt[2], popt[3]), popt[0], popt[1], popt[2], popt[3]
716 717
717
718 718 def Gauss_fit_2(self,spc,x,nGauss):
719 719
720 720
@@ -750,13 +750,7 class Oblique_Gauss_Fit(Operation):
750 750 print("ERROR")
751 751
752 752 d = numpy.mean(y[-100:])
753
754 # define a least squares function to optimize
755 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 754 return gaussian(x, popt[0], popt[1], popt[2], popt[3]),popt[0], popt[1], popt[2], popt[3]
761 755
762 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 801 d = D
808 802
809 803 # fit
810
811 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 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 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 1474 pwcode = 1
820 1475
@@ -828,23 +1483,150 class Oblique_Gauss_Fit(Operation):
828 1483 dataOut.power = numpy.average(z, axis=1)
829 1484 dataOut.powerdB = 10 * numpy.log10(dataOut.power)
830 1485
831
832 1486 x = dataOut.getVelRange(0)
833 1487
834 1488 dataOut.Oblique_params = numpy.ones((1,7,dataOut.nHeights))*numpy.NAN
835 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 1509 dataOut.VelRange = x
838 1510
839 1511
840 l1=range(22,36)
841 l2=range(58,99)
842 1512
1513 #l1=range(22,36) #+62
1514 #l1=range(32,36)
1515 #l2=range(58,99) #+62
1516
1517 #if Hmin1 == None or Hmax1 == None or Hmin2 == None or Hmax2 == None:
1518
1519 minHei1 = 105.
1520 maxHei1 = 122.5
1521 maxHei1 = 130.5
1522
1523 if mode == 10: #150 km
1524 minHei1 = 100
1525 maxHei1 = 100
1526
1527 inda1 = numpy.where(dataOut.heightList >= minHei1)
1528 indb1 = numpy.where(dataOut.heightList <= maxHei1)
1529
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:
843 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
844 1572
845 1573 try:
846 1574 spc = dataOut.data_spc[0,:,hei]
847 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:
848 1630 spc_fit, A1, B1, C1, D1 = self.Gauss_fit_2(spc,x,'first')
849 1631
850 1632 spc_diff = spc - spc_fit
@@ -854,13 +1636,228 class Oblique_Gauss_Fit(Operation):
854 1636
855 1637 D = (D1+D2)
856 1638
1639 if mode == 0: #Double Fit
857 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)
858 1641 #spc_double_fit,dataOut.Oblique_params = self.Double_Gauss_fit(spc,x,A1,B1,C1,A2,B2,C2,D)
859 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
860 1669 except:
861 1670 ###dataOut.Oblique_params[0,:,hei] = dataOut.Oblique_params[0,:,hei]*numpy.NAN
862 1671 pass
863 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 '''
1860
864 1861 return dataOut
865 1862
866 1863 class PrecipitationProc(Operation):
@@ -1621,6 +2618,8 class SpectralMoments(Operation):
1621 2618 if (nicoh is None): nicoh = 1
1622 2619 if (graph is None): graph = 0
1623 2620 if (smooth is None): smooth = 0
2621 elif (self.smooth < 3): smooth = 0
2622
1624 2623 if (type1 is None): type1 = 0
1625 2624 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
1626 2625 if (snrth is None): snrth = -20.0
General Comments 0
You need to be logged in to leave comments. Login now