1   /* Copyright 2002-2025 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.orekit.estimation.measurements.filtering;
19  
20  import java.io.BufferedReader;
21  import java.io.File;
22  import java.io.IOException;
23  import java.nio.charset.StandardCharsets;
24  import java.nio.file.Files;
25  import java.nio.file.Path;
26  import java.nio.file.Paths;
27  import java.util.ArrayList;
28  import java.util.List;
29  
30  import org.hipparchus.stat.descriptive.DescriptiveStatistics;
31  import org.junit.jupiter.api.Assertions;
32  import org.junit.jupiter.api.BeforeEach;
33  import org.junit.jupiter.api.Test;
34  import org.orekit.Utils;
35  import org.orekit.data.DataFilter;
36  import org.orekit.data.DataSource;
37  import org.orekit.files.rinex.HatanakaCompressFilter;
38  import org.orekit.files.rinex.observation.ObservationData;
39  import org.orekit.files.rinex.observation.ObservationDataSet;
40  import org.orekit.files.rinex.observation.RinexObservationParser;
41  import org.orekit.gnss.MeasurementType;
42  import org.orekit.gnss.ObservationType;
43  import org.orekit.gnss.PredefinedObservationType;
44  import org.orekit.gnss.SatInSystem;
45  import org.orekit.gnss.SatelliteSystem;
46  
47  public class PseudoRangeFilteringTest {
48  
49      private String baseName;
50      @BeforeEach
51      public void setUp() {
52          baseName = "src/test/resources/gnss/filtering/";
53          Utils.setDataRoot("regular-data:gnss/filtering");
54      }
55  
56      @Test
57      public void testHatchDoppler() throws IOException {
58  
59          // Definition of the ObservationTypes to study
60          ObservationType dopplerType = PredefinedObservationType.D1C;
61          ObservationType rangeType = PredefinedObservationType.C1C;
62  
63          // Definition of the Satellite to study
64          SatelliteSystem system = SatelliteSystem.GPS;
65          int prnNumber = 1;
66  
67  
68          String fileName = "AGGO00ARG_S_20190250000_15M_01S_MO.crx";
69          File file  = new File(baseName + fileName);
70          DataFilter filter = new HatanakaCompressFilter();
71          DataSource nd = new DataSource(file);
72          nd = filter.filter(nd);
73          final RinexObservationParser parser = new RinexObservationParser();
74          List<ObservationDataSet> listObsDataSet = parser.parse(nd).getObservationDataSets();
75          ObservationDataSet lastObsDataSet = listObsDataSet.get(listObsDataSet.size() - 1);
76  
77          // Test reset and null condition on doppler
78          ObservationData obsDataRange = new ObservationData(rangeType, 10, 0, 7);
79          ObservationData obsDataDopplerNull = new ObservationData(dopplerType, Double.NaN, 0, 7);
80          List<ObservationData> listObsData = new ArrayList<>();
81          listObsData.add(obsDataDopplerNull);
82          listObsData.add(obsDataRange);
83          ObservationDataSet obsDataSetNullDoppler = new ObservationDataSet(new SatInSystem(system, prnNumber),
84                                                                            lastObsDataSet.getDate(), 0, 0.0, listObsData);
85  
86          ObservationData obsDataRangeNull = new ObservationData(rangeType, 10, 0, 0);
87          ObservationData obsDataDoppler = new ObservationData(dopplerType, Double.NaN, 0, 0);
88          List<ObservationData> listObsData2 = new ArrayList<>();
89          listObsData2.add(obsDataDoppler);
90          listObsData2.add(obsDataRangeNull);
91          ObservationDataSet obsDataSetNullRange= new ObservationDataSet(new SatInSystem(system, prnNumber),
92                                                                         lastObsDataSet.getDate(), 0, 0.0, listObsData2);
93  
94          List<ObservationDataSet> copiedListObsDataSet = new ArrayList<>(listObsDataSet);
95          copiedListObsDataSet.add(obsDataSetNullRange);
96          copiedListObsDataSet.add(obsDataSetNullDoppler);
97  
98          SingleFrequencySmoother prs = new SingleFrequencySmoother(MeasurementType.DOPPLER, 100.0, 1, 50.0);
99          prs.filterDataSet(copiedListObsDataSet, system, prnNumber, PredefinedObservationType.D1C);
100 
101         List<SmoothedObservationDataSet> listObsDataSetUpdate = prs.getFilteredDataMap().get(rangeType);
102 
103         double lastUpdatedValue = listObsDataSetUpdate.get(listObsDataSetUpdate.size() - 1).getSmoothedData().getValue();
104         Assertions.assertEquals(2.0650729099E7, lastUpdatedValue, 1E-6);
105 
106     }
107 
108     @Test
109     public void testHatchCarrierPhaseValues() throws IOException {
110         /* The data used for the validation was generated following the gLAB tutorial available at :
111          * https://gssc.esa.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_II.pdf
112          * The truncated gLAB processed file is given for tests.
113          * The following commands were used for generation :
114          * gLAB_linux -input:cfg meas.cfg -input:obs UPC33510.08O_trunc | gawk '{if ($6==3) print $0}' > upc3.meas
115             gawk '{print $4,$11-$14-3.09*($14-$16)}' upc3.meas > upc3.C1
116             gawk 'BEGIN{Ts=60}{if (NR>Ts){n=Ts}else{n=NR}; C1s=$11/n+(n-1)/n*(C1s+($14-L1p));L1p=$14; print $4,C1s-$14-3.09*($14-$16)}' upc3.meas > upc3.C1S_60
117             gawk 'BEGIN{Ts=60}{if (NR>Ts){n=Ts}else{n=NR};L1f=$14+3.09*($14-$16);C1s=$11/n+(n-1)/n*(C1s+(L1f-L1fp));L1fp=L1f;print $4,C1s-L1f}' upc3.meas > upc3.C1DFreeS_60
118             graph.py -f upc3.C1 -s- -l "C1: PRN G03" -f upc3.C1S_60 -s.- -l "C1 smoothed" -f upc3.C1DFreeS_60 -s- -l "C1 DFree smoothed" --xn 35000 --xx 40000 --yn 17 --yx 25
119          *
120          * The commands were slightly modified, but not for the core formula.
121          * The behaviour is graphically consistent with that of gLAB, still the difference
122          * in the mean result and RMS error might be due to the prealignment step realized
123          * in gLAB but not in Orekit.
124          * Therefore the data will be filtered and checked against a constant value.
125          *
126          * This test uses observations that produces low divergence on the single frequency
127          * Hatch filter.
128          */
129 
130         ObservationType rangeType = PredefinedObservationType.C1;
131         ObservationType phaseTypeF1 = PredefinedObservationType.L1;
132         ObservationType phaseTypeF2 = PredefinedObservationType.L2;
133 
134         SatelliteSystem system = SatelliteSystem.GPS;
135         int prnNumber = 3;
136 
137         String fileName = "UPC33510.08O_trunc";
138         String fileName_gLAB_SF = "upc3.C1S_60";
139         String fileName_gLAB_DF = "upc3.C1DFreeS_60";
140         File file  = new File(baseName+fileName);
141 
142         DataSource nd = new DataSource(file);
143         final RinexObservationParser parser = new RinexObservationParser();
144 
145         // Test SatelliteSystem / SNR
146         List<ObservationDataSet> listObsDataSet = parser.parse(nd).getObservationDataSets();
147         ObservationDataSet lastObsDataSet = listObsDataSet.get(listObsDataSet.size() - 1);
148 
149         ObservationData obsDataRange = new ObservationData(rangeType, 0, 0, 0);
150         ObservationData obsDataRangeSNR = new ObservationData(rangeType, 0, 0, 1);
151         ObservationData obsDataF1 = new ObservationData(phaseTypeF1, 0, 0, 0);
152         ObservationData obsDataF2 = new ObservationData(phaseTypeF2, 0, 0, 0);
153 
154         List<ObservationData> listObsDataSatSystem = new ArrayList<>();
155         listObsDataSatSystem.add(obsDataF1);
156         listObsDataSatSystem.add(obsDataF2);
157         listObsDataSatSystem.add(obsDataRange);
158         List<ObservationData> listObsDataSNR = new ArrayList<>();
159         listObsDataSNR.add(obsDataF1);
160         listObsDataSNR.add(obsDataF2);
161         listObsDataSNR.add(obsDataRangeSNR);
162         ObservationDataSet obsDataSetRangeGLONASS = new ObservationDataSet(new SatInSystem(SatelliteSystem.GLONASS, prnNumber),
163                                                                            lastObsDataSet.getDate(), 0, 0.0, listObsDataSatSystem);
164 
165         ObservationDataSet obsDataSetRangeSNR = new ObservationDataSet(new SatInSystem(system, prnNumber),
166                                                                        lastObsDataSet.getDate(), 0, 0.0, listObsDataSNR);
167         //
168         List<ObservationDataSet> copiedListObsDataSet = new ArrayList<>(listObsDataSet);
169         copiedListObsDataSet.add(obsDataSetRangeGLONASS);
170         copiedListObsDataSet.add(obsDataSetRangeSNR);
171 
172         //
173         DualFrequencySmoother prs = new DualFrequencySmoother(100.0, 60);
174         prs.filterDataSet(copiedListObsDataSet, system, prnNumber, phaseTypeF1, phaseTypeF2);
175         SingleFrequencySmoother prsSF = new SingleFrequencySmoother(MeasurementType.CARRIER_PHASE, 100.0, 60, 50.0);
176         prsSF.filterDataSet(copiedListObsDataSet, system, prnNumber, phaseTypeF1);
177 
178         DualFrequencyHatchFilter filter = prs.getMapFilters().get(rangeType);
179         SingleFrequencyHatchFilter filterSF = prsSF.getMapFilters().get(rangeType);
180 
181         ArrayList<Double> filteredSF = filterSF.getSmoothedCodeHistory();
182         ArrayList<Double> filteredDF = filter.getSmoothedCodeHistory();
183         ArrayList<Double> gLAB_SF = readFile(baseName + fileName_gLAB_SF);
184         ArrayList<Double> gLAB_DF = readFile(baseName + fileName_gLAB_DF);
185 
186         ArrayList<Double> phaseArrayF1 = filter.getFirstFrequencyPhaseHistory();
187         ArrayList<Double> phaseArrayF2 = filter.getSecondFrequencyPhaseHistory();
188 
189         DescriptiveStatistics dSF = new DescriptiveStatistics();
190         DescriptiveStatistics dDF = new DescriptiveStatistics();
191         for (int i = 0; i < 5000; i++) {
192             double diffSF = gLAB_SF.get(i)-(filteredSF.get(i)-2.4931100000e7 - phaseArrayF1.get(i) - 3.09*(phaseArrayF1.get(i)-phaseArrayF2.get(i)));
193             double diffDF = gLAB_DF.get(i)-(filteredDF.get(i)-2.4931100000e7 - phaseArrayF1.get(i) - 3.09*(phaseArrayF1.get(i)-phaseArrayF2.get(i)));
194             dSF.addValue(diffSF);
195             dDF.addValue(diffDF);
196         }
197 
198         double rmsSF = dSF.getQuadraticMean();
199         double rmsDF = dDF.getQuadraticMean();
200 
201         // Non-Regression test : The value is above one due to a constant bias between the 2.
202         // The reason of the bias is to be explored, but might be related the pre alignement process
203         // performed by gLAB.
204         Assertions.assertTrue(rmsSF < 1.0063);
205         Assertions.assertTrue(rmsDF < 1.0060);
206     }
207 
208     @Test
209     public void testHatchCarrierPhase2() {
210         ObservationType rangeType = PredefinedObservationType.C1;
211         ObservationType phaseTypeF1 = PredefinedObservationType.L1;
212         ObservationType phaseTypeF2 = PredefinedObservationType.L2;
213 
214         SatelliteSystem system = SatelliteSystem.GPS;
215         int prnNumber = 7;
216 
217         String fileName = "irkm0440.16o";
218 
219         File file  = new File(baseName+fileName);
220 
221         DataSource nd = new DataSource(file);
222         final RinexObservationParser parser = new RinexObservationParser();
223 
224         DualFrequencySmoother prs = new DualFrequencySmoother(100.0, 60);
225         prs.filterDataSet(parser.parse(nd).getObservationDataSets(), system, prnNumber, phaseTypeF1, phaseTypeF2);
226 
227         DualFrequencyHatchFilter filterDF = prs.getMapFilters().get(rangeType);
228         ArrayList<Double> codeDFArray = filterDF.getCodeHistory();
229         ArrayList<Double> smoothedDFArray = filterDF.getSmoothedCodeHistory();
230 
231         double lastValueSmoothed = smoothedDFArray.get(smoothedDFArray.size()-1);
232         double lastValueCode = codeDFArray.get(codeDFArray.size()-1);
233 
234         // Regression test
235         Assertions.assertEquals(2.4715822416833777E7, lastValueSmoothed, 1e-6);
236         Assertions.assertEquals(2.4715823158E7, lastValueCode, 1e-4);
237 
238         ///// Test CarrierHatchFilterSingleFrequency
239 
240         SingleFrequencySmoother prsSF = new SingleFrequencySmoother(MeasurementType.CARRIER_PHASE, 100.0, 60, 50.0);
241         prsSF.filterDataSet(parser.parse(nd).getObservationDataSets(), system, prnNumber, phaseTypeF1);
242 
243         SingleFrequencyHatchFilter filterSF = prsSF.getMapFilters().get(rangeType);
244 
245         ArrayList<Double> codeSFArray = filterSF.getCodeHistory();
246         ArrayList<Double> smoothedSFArray = filterSF.getSmoothedCodeHistory();
247 
248         lastValueSmoothed = smoothedSFArray.get(smoothedSFArray.size()-1);
249         lastValueCode = codeSFArray.get(codeSFArray.size()-1);
250 
251         // Regression test
252         Assertions.assertEquals(2.4715820677129257E7, lastValueSmoothed, 1e-6);
253         Assertions.assertEquals(2.4715823158E7, lastValueCode, 1e-4);
254 
255         // Threshold test
256         Assertions.assertEquals(100.0, filterDF.getThreshold(), 1e-12);
257 
258         // Test getFilteredDataMap
259         List<SmoothedObservationDataSet> listObsDataSetUpdateDF = prs.getFilteredDataMap().get(rangeType);
260         List<SmoothedObservationDataSet> listObsDataSetUpdateSF = prsSF.getFilteredDataMap().get(rangeType);
261 
262         double lastUpdatedValueDF = listObsDataSetUpdateDF.get(listObsDataSetUpdateDF.size() - 1).getSmoothedData().getValue();
263         Assertions.assertEquals(2.4715822416833777E7, lastUpdatedValueDF, 1E-6);
264 
265         double lastUpdatedValueSF = listObsDataSetUpdateSF.get(listObsDataSetUpdateSF.size() - 1).getSmoothedData().getValue();
266         Assertions.assertEquals(2.4715820677129257E7, lastUpdatedValueSF, 1E-6);
267 
268     }
269 
270     public ArrayList<Double> readFile(final String fileName) throws IOException {
271 
272         final ArrayList<Double> valueArray = new ArrayList<>();
273         int cpt = 0;
274         Path pathToFile = Paths.get(fileName);
275 
276         try (BufferedReader br = Files.newBufferedReader(pathToFile, StandardCharsets.UTF_8)) {
277             String line = br.readLine();
278             while (line != null) {
279                 String[] splitLine = line.split("\\s+");
280                 valueArray.add(Double.parseDouble(splitLine[1]));
281                 cpt = cpt+1;
282                 line = br.readLine();
283             }
284         }
285         return valueArray;
286     }
287 
288 }