Disambiguated dataFilter by defining two columns.

One for the field to be filtered, one for the MATLAB expression.
Fixed week rollover bug.
Other small changes.

opensource/CheckDataFilter.m

@@ -0,0 +1,86 @@
+function [bOk] = CheckDataFilter(dataFilter,header)
+% [bOk] = CheckDataFilter(dataFilter,[header])
+% Check that dataFilter is defined correctly
+%
+% Rule for setting dataFilter.
+%   dataFilter values must be defined in pairs:
+%   dataFilter{i,1} is a string with one of 'Raw' header values from the 
+%                    GnssLogger log file e.g. 'ConstellationType'
+%   dataFilter{i,2} is a string with a valid matlab expression, containing
+%                   the header value, e.g. 'ConstellationType==1'
+% Any comparison must be against a scalar value.
+% The heading type may be repeated, for example,
+% dataFilter{i,2} = 'ConstellationType==1 | ConstellationType==3';
+% but you may not have different heading types in a single dataFilter cell
+%
+% header, [optional input], is  an mx1 cell array of strings containing
+% 'Raw' header values from GnssLogger log file
+%
+% CheckDataFilter checks the consistency of dataFilter.
+% if header is provided, then dataFilter is checked against it.
+
+%Author: Frank van Diggelen
+%Open Source code for processing Android GNSS Measurements
+
+bOk = true;
+
+%% check the basic structure of dataFilter
+if isempty(dataFilter)
+    return
+end
+[N,M] = size(dataFilter);
+if M~=2 || ~iscell(dataFilter)
+    error('dataFilter must be an nx2 cell array\n')
+end
+for i=1:N
+    for j=1:2
+        if ~ischar(dataFilter{i,j})
+            error('dataFilter{%d,%d} is not a string\n',i,j);
+        end
+    end
+end
+
+%% Check that the value in dataFilter{i,1} occurs in dataFilter{i,2}
+for i=1:N
+    if ~any(strfind(dataFilter{i,2},dataFilter{i,1}))
+        error('dataFilter{%d,1} string, ''%s'' not found in dataFilter{%d,2}\n',...
+            i,dataFilter{i,1},i);
+    end
+end
+if nargin<2 || isempty(header)
+    return
+end
+
+%% check that dataFilter has a matching value in the header, 
+for i=1:N
+    iMatch = strcmp(dataFilter{i,1},header); %iMatch is logical array
+    if ~any(iMatch) %no match found
+        error('dataFilter value ''%s'' has no matches in log file header',...
+            dataFilter{i,1});
+    end
+end
+
+%TBD check for occurrence of two different header types in dataFilter{i,2}
+% this is  little tricky, because if you use strfind you will find
+% supersets of shorter header types: like FullBiasNanos and BiasNanos
+% so you have to take care of that.
+% Also, this is spoon feeding the user - if they follow the rules defined
+% above, they won't do this. So, we leave it for now and add later when we
+% have nothing else to do (ha ha).
+
+end% of function CheckDataFilter
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Copyright 2016 Google Inc.
+% 
+% Licensed under the Apache License, Version 2.0 (the "License");
+% you may not use this file except in compliance with the License.
+% You may obtain a copy of the License at
+% 
+%      http://www.apache.org/licenses/LICENSE-2.0
+% 
+% Unless required by applicable law or agreed to in writing, software
+% distributed under the License is distributed on an "AS IS" BASIS,
+% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+% See the License for the specific language governing permissions and
+% limitations under the License.

opensource/Gps2Utc.m

@@ -1,14 +1,18 @@
-function [utcTime] = Gps2Utc(gpsTime)
-% [utcTime] = Gps2Utc(gpsTime)
-%   Convert GPS time (week & seconds) to UTC
+function [utcTime] = Gps2Utc(gpsTime,fctSeconds)
+% [utcTime] = Gps2Utc(gpsTime,[fctSeconds])
+%   Convert GPS time (week & seconds), or Full Cycle Time (seconds) to UTC
 % 
 % Input: gpsTime, [mx2] matrix [gpsWeek, gpsSeconds], 
+%        fctSeconds, [optional] Full Cycle Time (seconds)
 %
 % Outputs: utcTime, [mx6] matrix = [year,month,day,hours,minutes,seconds]
 % 
+% If fctSeconds is provided, gpsTime is ignored
+%
 % Valid range of inputs: 
 %   gps times corresponding to 1980/6/1 <= time < 2100/1/1
 %   i.e. [0,0] <= gpsTime < [6260, 432000]
+%        0 <= fctSeconds < 3786480000
 %
 % See also: Utc2Gps
 
@@ -32,10 +36,13 @@ function [utcTime] = Gps2Utc(gpsTime)
 % 5) if ls1~=ls, convert (gpsTime-ls1) to UTC Time
 
 %% Check inputs
-if size(gpsTime,2)~=2
+if nargin<2 && size(gpsTime,2)~=2
     error('gpsTime must have two columns')
 end
-fctSeconds = gpsTime*[GpsConstants.WEEKSEC; 1];
+if nargin<2
+    fctSeconds = gpsTime*[GpsConstants.WEEKSEC; 1];
+end
+
 %fct at 2100/1/1 00:00:00, not counting leap seconds:
 fct2100 =  [6260, 432000]*[GpsConstants.WEEKSEC; 1];
 if any(fctSeconds<0) || any (fctSeconds >= fct2100)

opensource/GpsAdrResiduals.m

@@ -34,7 +34,7 @@ if ~any(any(isfinite(gnssMeas.AdrM) & gnssMeas.AdrM~=0))
     return
 end
 if nargin<3 || isempty(llaDegDegM)
-    fprintf('\nGpsAdrWlsPvt needs the true position: llaDegDegM')
+    fprintf('GpsAdrResiduals needs the true position: llaDegDegM\n')
     return
 end
 xyz0M = Lla2Xyz(llaDegDegM);

opensource/GpsWlsPvt.m

@@ -33,6 +33,16 @@ end
 xo =zeros(8,1);%initial state: [center of the Earth, bc=0, velocities = 0]'
 
 weekNum     = floor(gnssMeas.FctSeconds/GpsConstants.WEEKSEC);
+%TBD check for week rollover here (it is checked in ProcessGnssMeas, but
+%this function should stand alone, so we should check again, and adjust 
+%tRxSeconds by +- a week if necessary)
+%btw, Q. why not just carry around fct and not worry about the hassle of
+%weeknumber, and the associated week rollover problems?
+% A. because you cannot get better than 1000ns (1 microsecond) precsision
+% when you put fct into a double. And that would cause errors of ~800m/s * 1us
+% (satellite range rate * time error) ~ 1mm in the range residual computation
+% So what? well, if you start processing with carrier phase, these errors
+% could accumulate.
 
 N = length(gnssMeas.FctSeconds);
 gpsPvt.FctSeconds      = gnssMeas.FctSeconds;

opensource/Lla2Ned.m

@@ -34,7 +34,7 @@ refXyz  = (xyz1M+xyz2M)/2;
 northM  = zeros(m1,1);
 eastM   = zeros(m1,1);
 for i=1:m1
-    Ce2n = RotEcef2Ned(llaDegDegM(i,1),llaDegDegM(1,2));
+    Ce2n = RotEcef2Ned(llaDegDegM(i,1),llaDegDegM(i,2));
     v = Ce2n*(xyz1M(i,:)-xyz2M(i,:))';
     northM(i)=v(1);
     eastM(i)=v(2);

opensource/PlotAdr.m

@@ -1,6 +1,6 @@
 function [colorsOut]= PlotAdr(gnssMeas,prFileName,colors)
 % [colors] = PlotAdr(gnssMeas,[prFileName],[colors])
-% plot the Accumulated Delta Ranges obtained from ProcessAdr
+% plot Valid Accumulated Delta Ranges obtained from ProcessAdr
 %
 %gnssMeas.FctSeconds = Nx1 vector. Rx time tag of measurements.
 %        .ClkDCount  = Nx1 vector. Hw clock discontinuity count
@@ -38,8 +38,19 @@ timeSeconds =gnssMeas.FctSeconds-gnssMeas.FctSeconds(1);%elapsed time in seconds
 for j=1:M %loop over Svid
     %% plot AdrM
     h123(1) = subplot(5,1,1:2); grid on, hold on,
-    AdrM = gnssMeas.AdrM(:,j);%local variable for convenience
-    iFi = find(isfinite(AdrM));
+    AdrM = gnssMeas.AdrM(:,j);%local variables for convenience
+    AdrState = gnssMeas.AdrState(:,j);
+
+    %From gps.h:
+    %/* However, it is expected that the data is only accurate when:
+    % *  'accumulated delta range state' == GPS_ADR_STATE_VALID.
+    %*/
+    % #define GPS_ADR_STATE_UNKNOWN                       0
+    % #define GPS_ADR_STATE_VALID                     (1<<0)
+    % #define GPS_ADR_STATE_RESET                     (1<<1)
+    % #define GPS_ADR_STATE_CYCLE_SLIP                (1<<2)
+    iValid = bitand(AdrState,2^0);
+    iFi = find(isfinite(AdrM) & iValid);
     if ~isempty(iFi)
         ti = timeSeconds(iFi(end));
         h=plot(timeSeconds,AdrM); set(h,'Marker','.','MarkerSize',4)
@@ -57,7 +68,7 @@ for j=1:M %loop over Svid
     end
 end
 subplot(5,1,1:2); ax=axis;
-title('Accumulated Delta Range (= -k*carrier phase) vs time'), ylabel('(meters)')
+title('Valid Accumulated Delta Range (= -k*carrier phase) vs time'), ylabel('(meters)')
 subplot(5,1,3:4); set(gca,'XLim',ax(1:2));
 title('DelPrM - AdrM'),ylabel('(meters)')
 
@@ -83,7 +94,7 @@ if nargout>1
     colorsOut = colors;
 end
 
-end %end of function PlotPseudorangeRates
+end %end of function PlotAdr
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % Copyright 2016 Google Inc.

opensource/PlotAdrResids.m

@@ -17,7 +17,7 @@ function PlotAdrResids(adrResid,gnssMeas,prFileName,colors)
 %Open Source code for processing Android GNSS Measurements
 
 K = 5; %number of satellites to plot 
-if ~any(any(isfinite(adrResid.ResidM))) %Nothing but NaNs
+if isempty(adrResid) || ~any(any(isfinite(adrResid.ResidM))) %Nothing but NaNs
     fprintf(' No adr residuals to plot\n'), return
 end
 if nargin<2

opensource/PlotPvt.m

@@ -34,10 +34,13 @@ if bGotLlaTrue
     llaRef = llaTrueDegDegM;
 else
     llaRef = llaMed;
+    %print median lla so user can use it as reference position if they want
+    fprintf('Median llaDegDegM = [%.7f %.7f %.2f]\n',llaMed)
 end
 
 %% plot ne errors vs llaTrueDegDegM --------------------------------------------
-nedM = Lla2Ned(gpsPvt.allLlaDegDegM,llaRef);
+nedM = Lla2Ned(gpsPvt.allLlaDegDegM,llaRef);%keep the NaNs in for the plot
+%so we see a break in the lines where there was no position
 h123=subplot(4,1,1:2);
 h1 = plot(nedM(:,2),nedM(:,1));
 set(h1,'LineStyle','-','LineWidth',0.1,'Color',ltgray)
@@ -71,7 +74,7 @@ title(titleString);
 axis equal, grid on
 ylabel('North (m)'),xlabel('East (m)')
 % compute error distribution and plot circle
-distM = sqrt(sum(nedM(:,1:2).^2,2));
+distM = sqrt(sum(nedM(iFi,1:2).^2,2));%use only finite values here
 medM = median(distM);
 %plot a circle using 'rectangle' (yes really :)
 hr=rectangle('Position',[-1 -1 2 2]*medM,'Curvature',[1 1]);

opensource/PlotPvtStates.m

@@ -51,7 +51,7 @@ set(h(iMin),'VerticalAlignment','top');
 
 %plot common bias, in microseconds and meters ---------------------------------
 h1234(2)=subplot(4,1,2);
-iFi = isfinite(gnssPvt.allBcMeters);%index into finite results
+iFi = find(isfinite(gnssPvt.allBcMeters));%index into finite results
 if any(iFi)
     plot(tSeconds,gnssPvt.allBcMeters - gnssPvt.allBcMeters(iFi(1)))
     grid on

opensource/ProcessGnssMeas.m

@@ -81,8 +81,9 @@ State = gnssRaw.State(1);
 assert(bitand(State,2^0) &  bitand(State,2^3),...
   'gnssRaw.State(1) must have bits 0 and 3 true before calling ProcessGnssMeas')
 
-%tRxNanos is now since the beginning of the week
-assert(all(tRxNanos <= WEEKNANOS),'tRxNanos should be <= WEEKNANOS')
+%tRxNanos now since beginning of the week, unless we had a week rollover
+%assert(all(tRxNanos <= WEEKNANOS),'tRxNanos should be <= WEEKNANOS')
+%TBD check week rollover code, and add assert tRxNanos <= WEEKNANOS after
 assert(all(tRxNanos >= 0),'tRxNanos should be >= 0')
 
 %subtract the fractional offsets TimeOffsetNanos and BiasNanos:
@@ -90,8 +91,7 @@ tRxSeconds  = (double(tRxNanos)-gnssRaw.TimeOffsetNanos-gnssRaw.BiasNanos)*1e-9;
 tTxSeconds  = double(gnssRaw.ReceivedSvTimeNanos)*1e-9;
 
 %check for week rollover in tRxSeconds
-prSeconds  = tRxSeconds - tTxSeconds;
-prSeconds  = CheckGpsWeekRollover(prSeconds);
+[prSeconds,tRxSeconds]  = CheckGpsWeekRollover(tRxSeconds,tTxSeconds);
 %we are ready to compute pseudorange in meters:
 PrM         = prSeconds*GpsConstants.LIGHTSPEED;
 
@@ -203,22 +203,27 @@ gnssMeas.DelPrM = delPrM;
 end %of function GetDelPr
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-function prSeconds  = CheckGpsWeekRollover(prSeconds)
+function [prSeconds,tRxSeconds]  = CheckGpsWeekRollover(tRxSeconds,tTxSeconds)
 %utility function for ProcessGnssMeas
 
+prSeconds  = tRxSeconds - tTxSeconds;
+
 iRollover = prSeconds > GpsConstants.WEEKSEC/2;
 if any(iRollover)
-    fprintf('\nWARNING: week rollover detected in time tags. Adjusting ...')
+    fprintf('\nWARNING: week rollover detected in time tags. Adjusting ...\n')
     prS = prSeconds(iRollover);
-    prS = prS - round(prS/GpsConstants.WEEKSEC)*GpsConstants.WEEKSEC;
+    delS = round(prS/GpsConstants.WEEKSEC)*GpsConstants.WEEKSEC;
+    prS = prS - delS;
     %prS are in the range [-WEEKSEC/2 : WEEKSEC/2];
     %check that common bias is not huge (like, bigger than 10s)
     maxBiasSeconds = 10; 
     if any(prS>maxBiasSeconds)
-        error('Failed to correct week rollover')
+        error('Failed to correct week rollover\n')
     else
         prSeconds(iRollover) = prS; %put back into prSeconds vector
-        fprintf('\nCorrected week rollover')
+        %Now adjust tRxSeconds by the same amount:
+        tRxSeconds(iRollover) = tRxSeconds(iRollover) - delS;
+        fprintf('Corrected week rollover\n')
     end
 end
 %TBD Unit test this

opensource/ProcessGnssMeasScript.m

@@ -10,6 +10,7 @@ prFileName = 'pseudoranges_log_2016_06_30_21_26_07.txt'; %with duty cycling, no
 % 2) change 'dirName = ...' to match the local directory you are using:
 dirName = '~/Documents/MATLAB/gpstools/opensource/demoFiles';
 % 3) run ProcessGnssMeasScript.m script file 
+param.llaTrueDegDegM = [];
 
 %Author: Frank van Diggelen
 %Open Source code for processing Android GNSS Measurements
@@ -18,24 +19,23 @@ dirName = '~/Documents/MATLAB/gpstools/opensource/demoFiles';
 %To add your own data:
 % save data from GnssLogger App, and edit dirName and prFileName appropriately
 %dirName = 'put the full path for your directory here';
-%prFileName = 'put the pseudoranges log file name here';
+%prFileName = 'put the pseuoranges log file name here';
 
 %% parameters
-param.llaTrueDegDegM = [];
+%param.llaTrueDegDegM = [];
 %enter true WGS84 lla, if you know it:
 param.llaTrueDegDegM = [37.422578, -122.081678, -28];%Charleston Park Test Site
 
 %% Set the data filter and Read log file
 dataFilter = SetDataFilter;
 [gnssRaw,gnssAnalysis] = ReadGnssLogger(dirName,prFileName,dataFilter,param);
+if isempty(gnssRaw), return, end
 
 %% Get online ephemeris from Nasa ftp, first compute UTC Time from gnssRaw:
 fctSeconds = 1e-3*double(gnssRaw.allRxMillis(end));
-gpsTime = [0,0];
-gpsTime(1) = floor(fctSeconds./GpsConstants.WEEKSEC);
-gpsTime(2) = fctSeconds - gpsTime(1)*GpsConstants.WEEKSEC;
-utcTime = Gps2Utc(gpsTime);
+utcTime = Gps2Utc([],fctSeconds);
 allGpsEph = GetNasaHourlyEphemeris(utcTime,dirName);
+if isempty(allGpsEph), return, end
 
 %% process raw measurements, compute pseudoranges:
 [gnssMeas] = ProcessGnssMeas(gnssRaw);
@@ -75,7 +75,7 @@ end
 % you may not use this file except in compliance with the License.
 % You may obtain a copy of the License at
 % 
-%     http://www.apache.org/licenses/LICENSE-2.0
+%     http://www.apache.org/licenses/LICENSE-2.0
 % 
 % Unless required by applicable law or agreed to in writing, software
 % distributed under the License is distributed on an "AS IS" BASIS,

opensource/ReadGnssLogger.m

@@ -1,4 +1,4 @@
-function [gnssRaw,gnssAnalysis] = ReadGnssLogger(dirName,fileName,dataFilterIn,gnssAnalysis)
+function [gnssRaw,gnssAnalysis] = ReadGnssLogger(dirName,fileName,dataFilter,gnssAnalysis)
 %% [gnssRaw,gnssAnalysis]=ReadGnssLogger(dirName,fileName,dataFilter,gnssAnalysis);
 % Read the log file created by Gnss Logger App in Android
 % Compatible with Android release N
@@ -7,8 +7,12 @@ function [gnssRaw,gnssAnalysis] = ReadGnssLogger(dirName,fileName,dataFilterIn,g
 %  dirName = string with directory of fileName, 
 %                e.g. '~/Documents/MATLAB/Pseudoranges/2016-03-28'
 %  fileName = string with filename
-%  dataFilter = cell array, dataFilter.{i}=string with a valid matlab expression
-%              e.g. dataFilter{1} = 'ConstellationType==1'
+%  dataFilter = nx2 cell array of pairs of strings, 
+%      dataFilter{i,1} is a string with one of 'Raw' header values from the 
+%                      GnssLogger log file e.g. 'ConstellationType'
+%      dataFilter{i,2} is a string with a valid matlab expression, containing
+%                      the header value, e.g. 'ConstellationType==1'
+%  See SetDataFilter.m for full rules and examples of dataFilter.
 %
 % Output: 
 %  gnssRaw, all GnssClock and GnssMeasurement fields from log file, including:
@@ -44,7 +48,7 @@ function [gnssRaw,gnssAnalysis] = ReadGnssLogger(dirName,fileName,dataFilterIn,g
 gnssAnalysis.GnssClockErrors = 'GnssClock Errors.';
 gnssAnalysis.GnssMeasurementErrors = 'GnssMeasurement Errors.';
 gnssAnalysis.ApiPassFail = '';
-if nargin<3, dataFilterIn = []; end
+if nargin<3, dataFilter = []; end
 
 %% check we have the right kind of fileName
 extension = fileName(end-3:end);
@@ -57,7 +61,8 @@ rawCsvFile = MakeCsv(dirName,fileName);
 [header,C] = ReadRawCsv(rawCsvFile);
 
 %% apply dataFilter 
-[dataFilter] = CheckDataFilter(dataFilterIn,header);
+[bOk] = CheckDataFilter(dataFilter,header);
+if ~bOk, return, end
 C = FilterData(C,dataFilter,header);
 
 %% pack data into gnssRaw structure
@@ -67,8 +72,6 @@ C = FilterData(C,dataFilter,header);
 [gnssRaw,gnssAnalysis] = CheckGnssClock(gnssRaw,gnssAnalysis);
 gnssAnalysis = ReportMissingFields(gnssAnalysis,missing);
 
-%TBD on any early return, return gnssAnalysis.ApiPassFail  = 'explanation'
-% so that reporting tool reports why
 
 end %end of function ReadGnssLogger
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -84,7 +87,7 @@ function csvFileName = MakeCsv(dirName,fileName)
 if dirName(end)~='/'
     dirName = [dirName,'/']; %add /
 end
-csvFileName = [dirName,'prs.csv'];
+csvFileName = [dirName,'raw.csv'];
 if strcmp(fileName(end-3:end),'.csv')
     return %input file is a csv file, nothing more to do here
 end
@@ -146,6 +149,7 @@ end
 %     ' | sed -e ''s/true/1/'' -e ''s/false/0/'' -e ''s/# //'' ',...
 %     ' -e ''s/Raw,//'' ',... %replace "Raw," with nothing
 %     '-e ''s/(//g'' -e ''s/)//g'' > ',csvFileName]);
+%
 % On versions from v1.4.0.0 N:
 % grep on "Raw," replace alpha characters amongst the numbers,
 % remove parentheses in the header,
@@ -172,7 +176,6 @@ if isempty(line) %line should be -1 at eof
     error('\nError occurred while reading file %s\n',fileName)
 end
 
-
 end %end of function MakeCsv
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -227,73 +230,42 @@ for i=1:M
         formatSpec = sprintf('%s %%f',formatSpec);
     end
 end
+%for empty fields, enter 'NaN' into csv
 C = textscan(fid,formatSpec,'Delimiter',',','EmptyValue',NaN);
 fclose(fid);
 
 end% of function ReadRawCsv
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-
-function [dataFilter] = CheckDataFilter(dataFilterIn,header)
-%% check that dataFilter has matching values in the header, 
-% extract this header value and add it as a second column for datafilter
-dataFilter = dataFilterIn(:); %make dataFilter a column array
-if isempty(dataFilterIn)
-    return
-end
-N = length(dataFilter);
-%check that each value of dataFilter is valid: i.e. it contains one of the
-%header types from the GnssLogger file
-L = length(header);
-for i=1:N
-    foundInHeader = zeros(1,L);
-    for j=1:L
-        foundInHeader(j) = any(strfind(dataFilter{i,1},header{j}));
-    end
-    if ~any(foundInHeader) %no match found (too coold)
-        error('dataFilter value ''%s'' has no matches in log file header',...
-            dataFilter{i,1});
-    elseif sum(foundInHeader)>1 % too many matches found (tooo hot)
-        error('dataFilter value ''%s'' has >1 match in log file header',...
-            dataFilter{i,1});
-    else %one match found (juust right)
-        k = find(foundInHeader);%index into where we found the matching header
-        dataFilter{i,2} = header{k};%save in second column
-    end
-end
-%dataFilter now has two columns: second one contains the matching header type
-end% of function CheckDataFilter
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 function C = FilterData(C,dataFilter,header)
 %% filter C based on contents of dataFilter
 
 iS = ones(size(C{1})); %initialize index into rows of C
 for i=1:size(dataFilter,1)
-    j=find(strcmp(header,dataFilter{i,2}));%j = index into header
+    j=find(strcmp(header,dataFilter{i,1}));%j = index into header
     %we should always be a value of j, because checkDataFilter checks for this:
     assert(any(j),'dataFilter{i} = %s not found in header\n',dataFilter{i,1})
     
-    %now we must evaluate the expression in dataFilter{i}, for example:
+    %now we must evaluate the expression in dataFilter{i,2}, for example:
     % 'BiasUncertaintyNanos < 1e7'
     %assign the relevant cell of C to a variable with same name as the header
     ts = sprintf('%s = C{%d};',header{j},j);
     eval(ts);
-    %create an index vector from the expression in dataFilter{i}
-    ts = sprintf('iSi = %s;',dataFilter{i,1});
+    %create an index vector from the expression in dataFilter{i,2}
+    ts = sprintf('iSi = %s;',dataFilter{i,2});
     eval(ts);
     
     %AND the iS index values on each iteration of i
     iS = iS & iSi;
 end
-%Check if filter removes all values,
+% Check if filter removes all values
 if ~any(iS) %if all zeros
-    fprintf('\nAll measurements removed. Specify dataFilter less strictly, ')
-    dataFilter(:,1)
+    fprintf('\nAll measurements removed. Specify dataFilter less strictly than this:, ')
+    dataFilter(:,2)
     return
 end
 
-%Now remove all values of C indexed by iS
+% Keep only those values of C indexed by iS
 for i=1:length(C)
     C{i} = C{i}(iS);
 end
@@ -355,6 +327,9 @@ for j = 1:length(header)
         missing.MeasurementFields{end+1} = header{j};
     end
 end
+%So, if a field is not reported, it will be all NaNs from makeCsv, and the above
+%code will not load it into gnssRaw. So when we call 'CheckGnssClock' it can
+%check for missing fields in gnssRaw.
 
 %TBD look for all zeros that can not legitimately be all zero, 
 %e.g. AccumulatedDeltaRangeMeters, and report these as missing data
@@ -362,42 +337,55 @@ end
 end %end of function PackGnssRaw
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-function [gnssRaw,gnssAnalysis] = CheckGnssClock(gnssRaw,gnssAnalysis)
+function [gnssRaw,gnssAnalysis,bOk] = CheckGnssClock(gnssRaw,gnssAnalysis)
 %% check clock values in gnssRaw
-
+bOk = true;
+sFail = ''; %initialize string to record failure messafes
 N = length(gnssRaw.ReceivedSvTimeNanos);
 %Insist on the presence of TimeNanos (time from hw clock)
 if ~isfield(gnssRaw,'TimeNanos')
-    error('TimeNanos data missing from GnssLogger file\n');
+    s = ' TimeNanos  missing from GnssLogger File.';
+    fprintf('WARNING: %s\n',s);
+    sFail = [sFail,s];
+    bOk = false;
 end
 if ~isfield(gnssRaw,'FullBiasNanos')
-    error('FullBiasNanos is missing or zero, we need it to get gnssRaw week\n')
-    %TBD change to fatal warning, so a report is still generated, with warning
+    s = 'FullBiasNanos missing from GnssLogger file.';
+    fprintf('WARNING: %s, we need it to get the week number\n',s);
+    sFail = [sFail,s];
+    bOk = false;
 end
 if ~isfield(gnssRaw,'BiasNanos')
     gnssRaw.BiasNanos = zeros(N,1);
 end
 if ~isfield(gnssRaw,'HardwareClockDiscontinuityCount')
+    %possibly non fatal error, we assume there is no hardware clock discontinuity
+    %so we set to zero and move on, but we print a warning
     gnssRaw.HardwareClockDiscontinuityCount = zeros(N,1);
-    fprintf('WARNING: Added HardwareClockDiscontinuityCount because it is missing from GNSS Logger file\n');
+    fprintf('WARNING: Added HardwareClockDiscontinuityCount=0 because it is missing from GNSS Logger file\n');
 end
-%auto-detect sign of FullBiasNanos, if it is positive, give warning and change
-if ~all(gnssRaw.FullBiasNanos<=0)
+
+%check FullBiasNanos, it should be negative values
+bChangeSign = any(gnssRaw.FullBiasNanos<0) & any(gnssRaw.FullBiasNanos>0);
+assert(~bChangeSign,...
+    'FullBiasNanos changes sign within log file, this should never happen');
+%Now we know FullBiasNanos doesnt change sign,auto-detect sign of FullBiasNanos, 
+%if it is positive, give warning and change
+if any(gnssRaw.FullBiasNanos>0)
     gnssRaw.FullBiasNanos = -1*gnssRaw.FullBiasNanos;
     fprintf('WARNING: FullBiasNanos wrong sign. Should be negative. Auto changing inside ReadGpsLogger\n');
     gnssAnalysis.GnssClockErrors = [gnssAnalysis.GnssClockErrors,...
         sprintf(' FullBiasNanos wrong sign.')];
 end
-%now all FullBiasNanos should be negative - if there are any are > 0 it means
-%something is very wrong with the log file, because FullBiasNanos has changed
-%sign from a large negative to large positive number, and we must assert
-assert(all(gnssRaw.FullBiasNanos<=0),...
-    'FullBiasNanos changes sign within log file, this should never happen')
 
 %compute full cycle time of measurement, in milliseonds
 gnssRaw.allRxMillis = int64((gnssRaw.TimeNanos - gnssRaw.FullBiasNanos)*1e-6);
 %allRxMillis is now accurate to one millisecond (because it's an integer)
 
+if ~bOk
+    gnssAnalysis.ApiPassFail = ['FAIL ',sFail]; 
+end
+
 end %end of function CheckGnssClock
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -428,10 +416,12 @@ if ~isempty(missing.MeasurementFields)
 end
 
 %assign pass/fail
-if isempty(missing.ClockFields) && isempty(missing.MeasurementFields)
-    gnssAnalysis.ApiPassFail = 'PASS';
-else
-    gnssAnalysis.ApiPassFail = 'FAIL BECAUSE OF MISSING FIELDS';
+if ~any(strfind(gnssAnalysis.ApiPassFail,'FAIL')) %didn't already fail
+    if isempty(missing.ClockFields) && isempty(missing.MeasurementFields)
+        gnssAnalysis.ApiPassFail = 'PASS';
+    else
+        gnssAnalysis.ApiPassFail = 'FAIL BECAUSE OF MISSING FIELDS';
+    end
 end
 end %end of function ReportMissingFields
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

opensource/SetDataFilter.m

@@ -1,36 +1,38 @@
 function dataFilter = SetDataFilter
-% dataFilter = SetDataFilter
+% dataFilter = SetDataFilter;
 % Function to set data filter for use with ReadGnssLogger
 % This function has no inputs. Edit it directly to change the data filter
 %
-% Rule for setting dataFilter:
-% dataFilter{i} is a cell with a string with a valid matlab expression, 
-% containing one of the 'Raw' headings in the GnssLogger log file
-% and where any comparison is against a scalar value
-%
-% the heading type may be repeated, for example,
-% dataFilter{i} = 'ConstellationType==1 | ConstellationType==3';
-% but you may not have more than one heading type in a single dataFilter{i} cell
+% Rule for setting dataFilter: see CheckDataFilter.m
 
 %Author: Frank van Diggelen
 %Open Source code for processing Android GNSS Measurements
 
 %filter for fine time measurements only  <=> uncertainty < 10 ms = 1e7 ns
-dataFilter{1} = 'BiasUncertaintyNanos < 1e7'; 
+dataFilter{1,1} = 'BiasUncertaintyNanos'; 
+dataFilter{1,2} = 'BiasUncertaintyNanos < 1e7'; 
+
+%filter out FullBiasNanos == 0
+dataFilter{end+1,1} = 'FullBiasNanos'; 
+dataFilter{end,2}   = 'FullBiasNanos ~= 0'; 
 
 %you can create other filters in the same way ...
 %for example, suppose you want to remove Svid 23:
-% dataFilter{end+1} = 'Svid~=23';
+% dataFilter{end+1,1} = 'Svid'; 
+% dataFilter{end,2}   = 'Svid ~= 23';
+% 
 %or suppose you want to keep only Svid 2,5,10, & 17
-% dataFilter{end+1} = 'Svid==2 | Svid==5 | Svid==10 | Svid==17';
+% dataFilter{end,2} = 'Svid==2 | Svid==5 | Svid==10 | Svid==17';
 % NOTE: you *cannot* use 'any(Svid)==[2,5,10,17]' because Svid refers to a 
 % vector variable and you must compare it to a scalar.
 
 %keep only Svid 2
-% dataFilter{end+1} = 'Svid==2';
+% dataFilter{end+1,1} = 'Svid'; 
+% dataFilter{end,2}   = 'Svid==2';
 
 %limit to GPS only:
-dataFilter{end+1} = 'ConstellationType==1'; 
+dataFilter{end+1,1} = 'ConstellationType'; 
+dataFilter{end,2}   = 'ConstellationType==1'; 
 %ConstellationType values are defined in Android HAL Documentation, gps.h, 
 %   typedef uint8_t                         GnssConstellationType;
 %   #define GNSS_CONSTELLATION_UNKNOWN      0
@@ -42,8 +44,8 @@ dataFilter{end+1} = 'ConstellationType==1';
 %   #define GNSS_CONSTELLATION_GALILEO      6
 
 %Example of how to select satellites from GPS and GLONASS:
-% dataFilter{end+1} = 'ConstellationType==1 | ConstellationType==3'; 
-%You may use the heading value (e.g. 'ConstellationType') more than once,
+% dataFilter{end+1} = '(ConstellationType)==1 | (ConstellationType)==3'; 
+%You may use the heading value e.g. '(ConstellationType)' more than once,
 %so long as you dont use different heading types in one dataFilter{} entry
 
 %bitwise data filters
@@ -69,7 +71,8 @@ dataFilter{end+1} = 'ConstellationType==1';
 
 %Example of how to use dataFilter for GnssMeasurementState 'State' bit fields: 
 %filter on GPS measurements with Code lock & TOW decoded:
-dataFilter{end+1} = 'bitand(State,2^0) &  bitand(State,2^3)';
+dataFilter{end+1,1} = 'State';
+dataFilter{end,2}   = 'bitand(State,2^0) & bitand(State,2^3)';
 % GNSS_MEASUREMENT_STATE_CODE_LOCK & GNSS_MEASUREMENT_STATE_TOW_DECODED
 
 % GnssAccumulatedDeltaRangeState values are defined gps.h, 
@@ -81,7 +84,8 @@ dataFilter{end+1} = 'bitand(State,2^0) &  bitand(State,2^3)';
 %
 %Example of how to use dataFilter for ADR State bit fields: 
 % get AccumulatedDeltaRangeState values that are valid
-% dataFilter{end+1} = 'bitand(AccumulatedDeltaRangeState,1)'; 
+% dataFilter{end+1,1} = 'AccumulatedDeltaRangeState'; 
+% dataFilter{end,2}   = 'bitand(AccumulatedDeltaRangeState,1)'; 
 end %end of function SetDataFilter
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%