1111
1212from gnss_lib_py .algorithms .residuals import solve_residuals
1313from gnss_lib_py .algorithms .snapshot import solve_wls
14+ from gnss_lib_py .utils .coordinates import geodetic_to_ecef
1415
1516def solve_fde (navdata , method = "residual" , remove_outliers = False ,
1617 max_faults = None , threshold = None ,verbose = False ,
@@ -49,6 +50,30 @@ def solve_fde(navdata, method="residual", remove_outliers=False,
4950
5051 """
5152
53+ # corr_pr_m_new = []
54+ # for _, _, navdata_subset in navdata.loop_time("gps_millis"):
55+ # # add ECEF coordinates
56+ # ecef_xyz = geodetic_to_ecef(navdata_subset[["lat_rx_gt_deg",
57+ # "lon_rx_gt_deg",
58+ # "alt_rx_gt_m"
59+ # ]])
60+ # navdata_subset["x_rx_gt_m"] = ecef_xyz[0,:]
61+ # navdata_subset["y_rx_gt_m"] = ecef_xyz[1,:]
62+ # navdata_subset["z_rx_gt_m"] = ecef_xyz[2,:]
63+ # # data.to_csv("example_smartloc_pre_wls.csv")
64+ #
65+ # # estimate receiver clock bias
66+ # wls_state_estimate = solve_wls(navdata_subset,
67+ # only_bias = True,
68+ # # delta_t_decimals=6,
69+ # receiver_state=navdata_subset,
70+ # )
71+ #
72+ # corr_pr_m_new += list(navdata_subset["corr_pr_m"] - wls_state_estimate["b_rx_wls_m"])
73+ #
74+ # navdata["corr_pr_m"] = corr_pr_m_new
75+
76+
5277 if method == "residual" :
5378 navdata = fde_residual_old (navdata , max_faults = max_faults ,
5479 threshold = threshold ,verbose = verbose ,
@@ -111,7 +136,7 @@ def fde_edm(navdata, max_faults=None, threshold=1.0, verbose=False,
111136 MIN_SATELLITES = 4
112137
113138 if threshold is None :
114- threshold = 1E6
139+ threshold = 0.6
115140
116141 # number of check indexes
117142 if max_faults is None :
@@ -127,17 +152,24 @@ def fde_edm(navdata, max_faults=None, threshold=1.0, verbose=False,
127152 time_start = time .time ()
128153
129154 nsl = len (navdata_subset )
130- if verbose :
131- print ("gt faults:" ,np .argwhere (navdata_subset ["fault_gt" ]== 1 )[:,0 ])
132155 sv_m = navdata_subset [["x_sv_m" ,"y_sv_m" ,"z_sv_m" ]]
133156 corr_pr_m = navdata_subset ["corr_pr_m" ]
134157
135158 # remove NaN indexes
136159 nan_idxs = sorted (list (set (np .arange (nsl )[np .isnan (sv_m ).any (axis = 0 )]).union ( \
137160 set (np .arange (nsl )[np .isnan (corr_pr_m )]).union ( \
138161 ))))[::- 1 ]
162+ navdata_subset .remove (cols = nan_idxs ,inplace = True )
163+ sv_m = navdata_subset [["x_sv_m" ,"y_sv_m" ,"z_sv_m" ]]
164+ corr_pr_m = navdata_subset ["corr_pr_m" ]
165+
166+ if verbose :
167+ print ("nan_idxs:" ,nan_idxs )
168+ print ("gt faults:" ,np .argwhere (navdata_subset ["NLOS (0 == no, 1 == yes, 2 == No Information)" ]== 1 )[:,0 ])
169+
139170 fault_idxs = []
140- orig_idxs = np .arange (nsl + 1 ) # add one for receiver index
171+ orig_idxs = np .arange (len (navdata_subset )+ 1 ) # add one for receiver index
172+ pre_nan_idxs = np .delete (np .arange (nsl + 1 ), nan_idxs )
141173
142174 edm = _edm_from_satellites_ranges (sv_m ,corr_pr_m )
143175
@@ -151,17 +183,19 @@ def fde_edm(navdata, max_faults=None, threshold=1.0, verbose=False,
151183 or (max_faults is not None and len (fault_idxs ) == max_faults ):
152184 break
153185
154-
155186 try :
156187 edm_detect = np .delete (np .delete (edm ,fault_idxs ,0 ),fault_idxs ,1 )
157188 detection_statistic_detect , svd_u , svd_s , svd_v = _edm_detection_statistic (edm_detect )
158189 detection_statistic = detection_statistic_detect [0 ]
159190 except Exception as exception :
160191 if verbose :
192+ print ("edm detect:" ,edm_detect )
161193 print (exception )
162- print (exception )
163194 break
164195
196+ if verbose :
197+ print ("detection statistic:" ,detection_statistic )
198+
165199 if detection_statistic < threshold :
166200 if verbose :
167201 print ("below threshold" )
@@ -290,9 +324,9 @@ def fde_edm(navdata, max_faults=None, threshold=1.0, verbose=False,
290324 detection_statistic = detection_statistic_exclude [min_idx ]
291325
292326 # important step! remove 1 since index 0 is the receiver index
293- fault_idxs = [i - 1 for i in fault_idxs ]
327+ fault_idxs = [pre_nan_idxs [ i - 1 ] for i in fault_idxs ]
294328
295- fault_edm_subset = np .array ([0 ] * len ( navdata_subset ) )
329+ fault_edm_subset = np .array ([0 ] * nsl )
296330 fault_edm_subset [fault_idxs ] = 1
297331 fault_edm_subset [nan_idxs ] = 2
298332 fault_edm += list (fault_edm_subset )
@@ -470,8 +504,6 @@ def fde_residual_old(navdata, max_faults, threshold,
470504 if verbose :
471505 print ("threshold fail:" )
472506 print ("r: " ,r )
473- for rri , rrr in enumerate (residuals ):
474- print (rri , rrr )
475507
476508 ri = list (ri )
477509
0 commit comments