playground
/
gps-measurement-tools
mirror of https://github.com/google/gps-measurement-tools.git
1
0
Fork 0
Code Issues Releases Activity
Browse Source

GPS Tools Open Source, MATLAB tools to: 

read data from GnssLogger App,
  compute and visualize pseudoranges,
  compute weighted least squares position and velocity,
  view and analyze carrier phase (if it is present in the log file).
pull/1/merge
Frank van Diggelen 9 years ago
commit
9607dfad52
49 changed files with 20958 additions and 0 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 202
      LICENSE.txt
  2. 1
      README.md
  3. 27
      opensource/CONTRIBUTING.txt
  4. 37
      opensource/CheckGpsEphInputs.m
  5. 64
      opensource/ClosestGpsEph.m
  6. 44
      opensource/CompareVersions.m
  7. 72
      opensource/Contents.m
  8. 34
      opensource/DayOfYear.m
  9. 48
      opensource/FlightTimeCorrection.m
  10. 183
      opensource/GetNasaHourlyEphemeris.m
  11. 34
      opensource/GnssThresholds.m
  12. 148
      opensource/Gps2Utc.m
  13. 148
      opensource/GpsAdrResiduals.m
  14. 41
      opensource/GpsConstants.m
  15. 89
      opensource/GpsEph2Dtsv.m
  16. 66
      opensource/GpsEph2Pvt.m
  17. 152
      opensource/GpsEph2Xyz.m
  18. 124
      opensource/GpsWlsPvt.m
  19. 54
      opensource/JulianDay.m
  20. 59
      opensource/Kepler.m
  21. 202
      opensource/LICENSE.txt
  22. 94
      opensource/LeapSeconds.m
  23. 62
      opensource/Lla2Ned.m
  24. 54
      opensource/Lla2Xyz.m
  25. 101
      opensource/PlotAdr.m
  26. 94
      opensource/PlotAdrResids.m
  27. 76
      opensource/PlotCno.m
  28. 122
      opensource/PlotPseudorangeRates.m
  29. 99
      opensource/PlotPseudoranges.m
  30. 125
      opensource/PlotPvt.m
  31. 150
      opensource/PlotPvtStates.m
  32. 81
      opensource/ProcessAdr.m
  33. 244
      opensource/ProcessGnssMeas.m
  34. 84
      opensource/ProcessGnssMeasScript.m
  35. 94
      opensource/README.md
  36. 73
      opensource/README.txt
  37. 451
      opensource/ReadGnssLogger.m
  38. 257
      opensource/ReadRinexNav.m
  39. 60
      opensource/RotEcef2Ned.m
  40. 101
      opensource/SetDataFilter.m
  41. 111
      opensource/Utc2Gps.m
  42. 172
      opensource/WlsPvt.m
  43. 82
      opensource/Xyz2Lla.m
  44. 3352
      opensource/demoFiles/hour1820.16n
  45. 3360
      opensource/demoFiles/hour2350.16n
  46. 1380
      opensource/demoFiles/prs.csv
  47. 1384
      opensource/demoFiles/prsWithPrInMeters.csv
  48. 1607
      opensource/demoFiles/pseudoranges_log_2016_06_30_21_26_07.txt
  49. 5259
      opensource/demoFiles/pseudoranges_log_2016_08_22_14_45_50.txt

202
LICENSE.txt
Unescape View File

@ -0,0 +1,202 @@
Apache License
Version 2.0, January 2004
http://www.apache.org/licenses/
TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
1. Definitions.
"License" shall mean the terms and conditions for use, reproduction,
and distribution as defined by Sections 1 through 9 of this document.
"Licensor" shall mean the copyright owner or entity authorized by
the copyright owner that is granting the License.
"Legal Entity" shall mean the union of the acting entity and all
other entities that control, are controlled by, or are under common
control with that entity. For the purposes of this definition,
"control" means (i) the power, direct or indirect, to cause the
direction or management of such entity, whether by contract or
otherwise, or (ii) ownership of fifty percent (50%) or more of the
outstanding shares, or (iii) beneficial ownership of such entity.
"You" (or "Your") shall mean an individual or Legal Entity
exercising permissions granted by this License.
"Source" form shall mean the preferred form for making modifications,
including but not limited to software source code, documentation
source, and configuration files.
"Object" form shall mean any form resulting from mechanical
transformation or translation of a Source form, including but
not limited to compiled object code, generated documentation,
and conversions to other media types.
"Work" shall mean the work of authorship, whether in Source or
Object form, made available under the License, as indicated by a
copyright notice that is included in or attached to the work
(an example is provided in the Appendix below).
"Derivative Works" shall mean any work, whether in Source or Object
form, that is based on (or derived from) the Work and for which the
editorial revisions, annotations, elaborations, or other modifications
represent, as a whole, an original work of authorship. For the purposes
of this License, Derivative Works shall not include works that remain
separable from, or merely link (or bind by name) to the interfaces of,
the Work and Derivative Works thereof.
"Contribution" shall mean any work of authorship, including
the original version of the Work and any modifications or additions
to that Work or Derivative Works thereof, that is intentionally
submitted to Licensor for inclusion in the Work by the copyright owner
or by an individual or Legal Entity authorized to submit on behalf of
the copyright owner. For the purposes of this definition, "submitted"
means any form of electronic, verbal, or written communication sent
to the Licensor or its representatives, including but not limited to
communication on electronic mailing lists, source code control systems,
and issue tracking systems that are managed by, or on behalf of, the
Licensor for the purpose of discussing and improving the Work, but
excluding communication that is conspicuously marked or otherwise
designated in writing by the copyright owner as "Not a Contribution."
"Contributor" shall mean Licensor and any individual or Legal Entity
on behalf of whom a Contribution has been received by Licensor and
subsequently incorporated within the Work.
2. Grant of Copyright License. Subject to the terms and conditions of
this License, each Contributor hereby grants to You a perpetual,
worldwide, non-exclusive, no-charge, royalty-free, irrevocable
copyright license to reproduce, prepare Derivative Works of,
publicly display, publicly perform, sublicense, and distribute the
Work and such Derivative Works in Source or Object form.
3. Grant of Patent License. Subject to the terms and conditions of
this License, each Contributor hereby grants to You a perpetual,
worldwide, non-exclusive, no-charge, royalty-free, irrevocable
(except as stated in this section) patent license to make, have made,
use, offer to sell, sell, import, and otherwise transfer the Work,
where such license applies only to those patent claims licensable
by such Contributor that are necessarily infringed by their
Contribution(s) alone or by combination of their Contribution(s)
with the Work to which such Contribution(s) was submitted. If You
institute patent litigation against any entity (including a
cross-claim or counterclaim in a lawsuit) alleging that the Work
or a Contribution incorporated within the Work constitutes direct
or contributory patent infringement, then any patent licenses
granted to You under this License for that Work shall terminate
as of the date such litigation is filed.
4. Redistribution. You may reproduce and distribute copies of the
Work or Derivative Works thereof in any medium, with or without
modifications, and in Source or Object form, provided that You
meet the following conditions:
(a) You must give any other recipients of the Work or
Derivative Works a copy of this License; and
(b) You must cause any modified files to carry prominent notices
stating that You changed the files; and
(c) You must retain, in the Source form of any Derivative Works
that You distribute, all copyright, patent, trademark, and
attribution notices from the Source form of the Work,
excluding those notices that do not pertain to any part of
the Derivative Works; and
(d) If the Work includes a "NOTICE" text file as part of its
distribution, then any Derivative Works that You distribute must
include a readable copy of the attribution notices contained
within such NOTICE file, excluding those notices that do not
pertain to any part of the Derivative Works, in at least one
of the following places: within a NOTICE text file distributed
as part of the Derivative Works; within the Source form or
documentation, if provided along with the Derivative Works; or,
within a display generated by the Derivative Works, if and
wherever such third-party notices normally appear. The contents
of the NOTICE file are for informational purposes only and
do not modify the License. You may add Your own attribution
notices within Derivative Works that You distribute, alongside
or as an addendum to the NOTICE text from the Work, provided
that such additional attribution notices cannot be construed
as modifying the License.
You may add Your own copyright statement to Your modifications and
may provide additional or different license terms and conditions
for use, reproduction, or distribution of Your modifications, or
for any such Derivative Works as a whole, provided Your use,
reproduction, and distribution of the Work otherwise complies with
the conditions stated in this License.
5. Submission of Contributions. Unless You explicitly state otherwise,
any Contribution intentionally submitted for inclusion in the Work
by You to the Licensor shall be under the terms and conditions of
this License, without any additional terms or conditions.
Notwithstanding the above, nothing herein shall supersede or modify
the terms of any separate license agreement you may have executed
with Licensor regarding such Contributions.
6. Trademarks. This License does not grant permission to use the trade
names, trademarks, service marks, or product names of the Licensor,
except as required for reasonable and customary use in describing the
origin of the Work and reproducing the content of the NOTICE file.
7. Disclaimer of Warranty. Unless required by applicable law or
agreed to in writing, Licensor provides the Work (and each
Contributor provides its Contributions) on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
implied, including, without limitation, any warranties or conditions
of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
PARTICULAR PURPOSE. You are solely responsible for determining the
appropriateness of using or redistributing the Work and assume any
risks associated with Your exercise of permissions under this License.
8. Limitation of Liability. In no event and under no legal theory,
whether in tort (including negligence), contract, or otherwise,
unless required by applicable law (such as deliberate and grossly
negligent acts) or agreed to in writing, shall any Contributor be
liable to You for damages, including any direct, indirect, special,
incidental, or consequential damages of any character arising as a
result of this License or out of the use or inability to use the
Work (including but not limited to damages for loss of goodwill,
work stoppage, computer failure or malfunction, or any and all
other commercial damages or losses), even if such Contributor
has been advised of the possibility of such damages.
9. Accepting Warranty or Additional Liability. While redistributing
the Work or Derivative Works thereof, You may choose to offer,
and charge a fee for, acceptance of support, warranty, indemnity,
or other liability obligations and/or rights consistent with this
License. However, in accepting such obligations, You may act only
on Your own behalf and on Your sole responsibility, not on behalf
of any other Contributor, and only if You agree to indemnify,
defend, and hold each Contributor harmless for any liability
incurred by, or claims asserted against, such Contributor by reason
of your accepting any such warranty or additional liability.
END OF TERMS AND CONDITIONS
APPENDIX: How to apply the Apache License to your work.
To apply the Apache License to your work, attach the following
boilerplate notice, with the fields enclosed by brackets "[]"
replaced with your own identifying information. (Don't include
the brackets!) The text should be enclosed in the appropriate
comment syntax for the file format. We also recommend that a
file or class name and description of purpose be included on the
same "printed page" as the copyright notice for easier
identification within third-party archives.
Copyright [yyyy] [name of copyright owner]
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.

1
README.md
Unescape View File

@ -0,0 +1 @@
opensource/README.md

27
opensource/CONTRIBUTING.txt
Unescape View File

@ -0,0 +1,27 @@
Want to contribute? Great! First, read this page (including the small print at the end).
### Before you contribute
Before we can use your code, you must sign the
[Google Individual Contributor License Agreement]
(https://cla.developers.google.com/about/google-individual)
(CLA), which you can do online. The CLA is necessary mainly because you own the
copyright to your changes, even after your contribution becomes part of our
codebase, so we need your permission to use and distribute your code. We also
need to be sure of various other things—for instance that you'll tell us if you
know that your code infringes on other people's patents. You don't have to sign
the CLA until after you've submitted your code for review and a member has
approved it, but you must do it before we can put your code into our codebase.
Before you start working on a larger contribution, you should get in touch with
us first through the issue tracker with your idea so that we can help out and
possibly guide you. Coordinating up front makes it much easier to avoid
frustration later on.
### Code reviews
All submissions, including submissions by project members, require review. We
use Github pull requests for this purpose.
### The small print
Contributions made by corporations are covered by a different agreement than
the one above, the
[Software Grant and Corporate Contributor License Agreement]
(https://cla.developers.google.com/about/google-corporate).

37
opensource/CheckGpsEphInputs.m
Unescape View File

@ -0,0 +1,37 @@
function [bOk,gpsEph,gpsWeek,ttxSec] = CheckGpsEphInputs(gpsEph,gpsTime)
%[bOk,gpsEph,gpsWeek,ttxSec] = CheckGpsEphInputs(gpsEph,gpsTime)
%check the inputs for GpsEph2Pvt, GpsEph2Xyz, GpsEph2Dtsv
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
bOk=false;
if ~isstruct(gpsEph)
error('gpsEph input must be a structure, as defined by ReadRinexNav')
end
p=length(gpsEph);
%Check that gpsTime is a px2 vector
if any(size(gpsTime) ~= [p 2])
error('gpsTime must be px2 [gpsWeek, gpsSec], where p =length(gpsEph)')
end
gpsWeek = gpsTime(:,1);
ttxSec = gpsTime(:,2);
bOk = true;
end %end of function CheckGpsEphInputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

64
opensource/ClosestGpsEph.m
Unescape View File

@ -0,0 +1,64 @@
function [gpsEph,iSv] = ClosestGpsEph(allGpsEph,svIds,fctSeconds)
%[gpsEph,iSv] = ClosestGpsEph(allGpsEph,svIds,fctSeconds);
%find ephemeris in a GPS ephemeris structure allGpsEph for all svIds listed
%return gpsEph = unique ephemeris for svIds, with fctToe closest to fctSeconds,
% and |fctToe - fctSeconds| < fitInterval
%
%output: gpsEph, iSv
% gpsEph = valid ephemeris corresponding to svIds(iSv)
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
gpsEph = allGpsEph(1);%initialize gpsEph
numEph = 0;
%find all ephemeris corresponding to svIds(i)
for i=1:length(svIds)
iThisSv = [allGpsEph.PRN] == svIds(i);
if any(iThisSv)
ephThisSv = allGpsEph(iThisSv);
%find Toe within fit interval
%set fit interval
fitIntervalHours = [ephThisSv.Fit_interval];
%Rinex says "Zero if not known", so adjust for zeros
fitIntervalHours(fitIntervalHours == 0) = 4;
%full cycle time of ephemeris Toe
fctToe = [ephThisSv.GPS_Week]*GpsConstants.WEEKSEC + [ephThisSv.Toe];
%find freshest Toe
[ageSeconds, iMin] = min(abs(fctToe - fctSeconds));
if ageSeconds < (fitIntervalHours/2)*3600;
numEph = numEph+1;
gpsEph(numEph) = ephThisSv(iMin);
iSv(numEph) = i;%save index into svIds
else
fprintf('No valid ephemeris found for svId %d,', svIds(i))
ageHours = (fctToe(iMin)-fctSeconds)/3600;
fprintf(' closest Toe is %.1f hours away.\n',ageHours)
end
end
end
if numEph==0
gpsEph=[];iSv=[];%return empty matrices
end
end %of function ClosestGpsEph
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

44
opensource/CompareVersions.m
Unescape View File

@ -0,0 +1,44 @@
function s = CompareVersions(v,w)
%s = CompareVersions(v,w);
%compares version v to w
% s = 'before','equal','after'
% 'before <=> version v is before version w
%
% v and w must be vectors of the same length
%
% example: v = [1,0,0,0]; w = [1,4,0,0]; s = 'before'
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
if length(v)~=length(w)
error('The two inputs must be scalars or vectors of the same length')
end
d = v(:)-w(:);
%find first d that differs from zero
i = find(d~=0);
if isempty(i)
s = 'equal';
elseif d(i(1)) < 0 %the first field that differs is the most significant
s = 'before';
else
s = 'after';
end
end %end function CompareVersions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

72
opensource/Contents.m
Unescape View File

@ -0,0 +1,72 @@
% GNSS Tools Open Source, MATLAB tools to:
% read data from GnssLogger App,
% compute and visualize pseudoranges,
% compute weighted least squares position and velocity,
% view and analyze carrier phase (if it is present in the log file).
%
% ProcessGnssMeasScript - script to set file name and call all other functions
% Start with this file, and the run the demo log files provided with this code
%
% Coordinate transformations
% Lla2Ned - Difference of [lat,lon,alt], return NED coords in meters
% Lla2Xyz - lat,lon,alt transform to x,y,z (Earth Centered Earth Fixed)
% RotEcef2Ned - Rotation matrix to convert ECEF vector to NED, & vice versa
% Xyz2Lla - x,y,z (Earth Centered Earth Fixed) transform to lat,lon,alt
%
% Ephemeris and Orbit functions
% CheckGpsEphInputs - Check the inputs for all GpsEph2* functions
% ClosestGpsEph - find unique fresh ephemeris from a GPS ephemeris structure
% GpsEph2Dtsv - Satellite clock bias from GPS ephemeris
% GpsEph2Pvt - Satellite position, velocity and clock bias from GPS ephemeris
% GpsEph2Xyz - Satellite position from GPS ephemeris
% FlightTimeCorrection - Rotated coords from Earth rotation during flight time
% Kepler - Solve Kepler's equation for eccentric anomaly
% ReadRinexNav - Read ephemeris & iono data from an ASCII formatted RINEX2 file
%
% Navigation, Pseudorange and Accumulated Delta Range functions
% GpsAdrResiduals - Residuals from GPS Accumulated Delta Ranges (carrier)
% GpsWlsPvt - Position Velocity and Time from GPS measurements
% WlsPvt - Weighted least squares PVT solution from pr and prr
% ProcessAdr - Compute Delta PR minus ADR (carrier)
% ProcessGnssMeas - Process raw GnssLogger measurements and compute pseudoranges
%
% Plotting functions
% PlotAdrResids - Plot the Accumulated Delta Range (carrier) residuals
% PlotCno - Plot the carrier-to-noise-density ratio, C/No, from gnssMeas
% PlotPseudoranges - Plot the pseudoranges obtained from ProcessGnssMeas
% PlotPvt - Plot the results of GpsWlsPvt
% PlotPvtStates - Plot Position, Velocity and Time/clock states
%
% Time functions
% DayOfYear - Day number of the year
% Gps2Utc - Convert GPS time (week & seconds) to UTC
% JulianDay - Number of days since first GPS week
% LeapSeconds - Number of leap seconds since the first GPS week
% Utc2Gps - Convert UTC time to GPS time
%
% File reading
% GetNasaHourlyEphemeris - Read hourly ephemeris file
% ReadRinexNav - Read ephemeris & iono data from a RINEX2 Nav file
% ReadGnssLogger - Read the file created by Gnss Logger App in Android
%
% General functions and classes
% CompareVersions - Compare two version numbers
% GnssThresholds - GNSS validity thresholds we use in the code
% GpsConstants - GPS constants, from WGS84 and IS-GPS-200
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

34
opensource/DayOfYear.m
Unescape View File

@ -0,0 +1,34 @@
function dayNumber = DayOfYear(utcTime)
%dayNumber = DayOfYear(utcTime);
%
% Return the day number of the year
% utcTime: [1x6] [year,month,day,hours,minutes,seconds]
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
if any(size(utcTime)~=[1,6])
error('utcTime must be 1x6 for DayOfYear function')
end
jDay = JulianDay([utcTime(1:3),0,0,0]);%midnight morning of this day
jDayJan1 = JulianDay([utcTime(1),1,1,0,0,0]);%midnight morning of Jan 1st
dayNumber = jDay - jDayJan1 + 1;
end %end of function DayOfYear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

48
opensource/FlightTimeCorrection.m
Unescape View File

@ -0,0 +1,48 @@
function xERot = FlightTimeCorrection(xE, dTflightSeconds)
%xERot = FlightTimeCorrection(xE, dtflight);
% Compute rotated satellite ECEF coordinates caused by Earth
% rotation during signal flight time
%
% Inputs:
% xE - satellite ECEF position at time of transmission
% dtflight - signal flight time (seconds)
%
% Outputs:
% XeRot - rotated satelite position vector (ECEF at trx)
%
% Reference:
% IS GPS 200, 20.3.3.4.3.3.2 Earth-Centered, Inertial (ECI) Coordinate System
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%See also "A software receiver for GPS and Galileo", K. Borre at al.
%Rotation angle (radians):
theta = GpsConstants.WE * dTflightSeconds;
%Apply rotation from IS GPS 200-E, 20.3.3.4.3.3.2
%Note: IS GPS 200-E shows the rotation from ecef to eci
% so our rotation R3, is in the opposite direction:
R3 = [ cos(theta) sin(theta) 0;
-sin(theta) cos(theta) 0;
0 0 1];
xERot = (R3*xE(:))';
end %end of function FlightTimeCorrection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

183
opensource/GetNasaHourlyEphemeris.m
Unescape View File

@ -0,0 +1,183 @@
function [allGpsEph,allGloEph] = GetNasaHourlyEphemeris(utcTime,dirName)
%[allGpsEph,allGloEph] = GetNasaHourlyEphemeris(utcTime,dirName)
%Get hourly ephemeris files,
% If a GPS ephemeris file is in dirName, with valid ephemeris for at
% least 24 svs, then read it; else download from NASA's archive of
% Space Geodesy Data
%
%Output allGpsEph structure of ephemerides in format defined in ReadRinexNav.m
% TBD: allGloEph, and any other supported by the NSA ftp site
%
%examples of how to call GetNasaHourlyEphemeris:
%1) with /demoFiles and utcTime where ephemeris has already been downloaded:
% replace '~/...' with actual path:
% dirName = '~/Documents/MATLAB/gpstools/opensource/demoFiles';
% utcTime = [2016,8,22,22,50,00];
% [allGpsEph,allGloEph] = GetNasaHourlyEphemeris(utcTime,dirName)
%
%2) with utcTime for which ephemeris has not been downloaded,
% this exercises the automatic ftp download.
% Replace '~' with actual path:
% dirName = '~/Documents/MATLAB/gpstools/opensource/demoFiles';
% utcTime = [2016,5,ceil(rand(1,2).*[30,24]),0,0],%a random Hour and Day in May
% [allGpsEph,allGloEph] = GetNasaHourlyEphemeris(utcTime,dirName)
%
% More details:
% The Nasa ephemeris Unix-compressed.
% GetNasaHourlyEphemeris will automatically uncompress it, if you have the right
% uncompress function on you computer. If you need to install an unzip utility,
% see http://cddis.nasa.gov/ and http://www.gpzip.org
% Search for 'uncompress' in the GetNasaHourlyEphemeris function to find and
% edit the name of the unzip utility.
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
allGpsEph=[];allGloEph=[];
if nargin<2
dirName = [];
end
[bOk,dirName] = CheckInputs(utcTime,dirName);
if ~bOk, return, end
%Description of file names, see:
%http://cddis.gsfc.nasa.gov/Data_and_Derived_Products/GNSS/...
% broadcast_ephemeris_data.html#GPShourly
yearNumber4Digit = utcTime(1);
yearNumber2Digit = rem(utcTime(1),100);
dayNumber = DayOfYear(utcTime);
hourlyZFile = sprintf('hour%03d0.%02dn.Z',dayNumber,yearNumber2Digit);
ephFilename = hourlyZFile(1:end-2);
fullEphFilename = [dirName,ephFilename]; %full name (with directory specified)
%check if ephemeris file already exists (e.g. you downloaded it 'by hand')
%and if there are fresh ephemeris for lotsa sats within 2 hours of fctSeconds
bGotGpsEph = false;
if exist(fullEphFilename,'file')==2
%% file exists locally, so read it
fprintf('Reading GPS ephemeris from ''%s'' file in local directory\n',...
ephFilename);
fprintf('%s\n',dirName);
allGpsEph = ReadRinexNav(fullEphFilename);
[~,fctSeconds] = Utc2Gps(utcTime);
ephAge = [allGpsEph.GPS_Week]*GpsConstants.WEEKSEC + [allGpsEph.Toe] - ...
fctSeconds;
%get index into fresh and healthy ephemeris (health bit set => unhealthy)
iFreshAndHealthy = abs(ephAge)<GpsConstants.EPHVALIDSECONDS & ...
~[allGpsEph.health];
%TBD look at allEph.Fit_interval, and deal with values > 0
goodEphSvs = unique([allGpsEph(iFreshAndHealthy).PRN]);
if length(goodEphSvs)>=GnssThresholds.MINNUMGPSEPH
bGotGpsEph = true;
end
end
if ~bGotGpsEph
%% get ephemeris from Nasa site
ftpSiteName = 'cddis.gsfc.nasa.gov';
hourlyDir=sprintf('/gnss/data/hourly/%4d/%03d/',yearNumber4Digit,dayNumber);
fprintf('\nGetting GPS ephemeris ''%s'' from NASA ftp site ...',hourlyZFile)
fprintf('\nftp://%s%s\n',ftpSiteName,hourlyDir);
%check that the dirName has been properly specified
if strfind(dirName,'~')
fprintf('\nYou set: dirName = ''%s''\n',dirName)
fprintf('To download ephemeris from ftp,')
fprintf(' specify ''dirName'' with full path, and no tilde ''~''\n')
fprintf('Change ''dirName'' or download ephemeris file ''%s'' by hand\n',...
hourlyZFile);
fprintf('To unzip the *.Z file, see http://www.gzip.org/\n')
return
end
try
nasaFtp=ftp(ftpSiteName); %connect to ftp server and create ftp object
cd(nasaFtp,hourlyDir);
zF = mget(nasaFtp,hourlyZFile,dirName); %get compressed (.Z) file
catch
%failed automatic ftp, ask user to do this by hand
fprintf('\nAutomatic ftp download failed.\n')
fprintf('Please go to this ftp, ftp://cddis.gsfc.nasa.gov/ \n')
fprintf('and get this file:%s%s \n',hourlyDir,hourlyZFile);
fprintf('Extract contents to the directory with your pseudorange log file:\n')
fprintf('%s\n',dirName)
fprintf('To unzip the *.Z file, see http://www.gzip.org/\n')
fprintf('Then run this function again, it will read from that directory\n')
return
end
%Now we have the zipped file from nasa. Unzip it:
unzipCommand='uncompress';%edit if your platform uses something different
if any(strfind(computer,'PCWIN'))
unzipCommand='gunzip';
end
s = sprintf('%s "%s"',unzipCommand,zF{1}); %string for system command
[sysFlag,result] = system(s);
if any(sysFlag)
fprintf('\nError in GetNasaHourlyEphemeris.m\n')
fprintf('%s',result);
fprintf('System command ''%s'' failed\n',unzipCommand)
fprintf('Unzip contents of %s by hand. See http://www.gzip.org/\n',zF{1})
fprintf('Then run this function again, it will read the uncompressed file\n')
return
end
bOk = bOk && ~sysFlag;
if bOk
fprintf('\nSuccessfully downloaded ephemeris file ''%s''\n',ephFilename)
if ~isempty(dirName)
fprintf('to: %s\n',dirName);
end
else
end
end
allGpsEph = ReadRinexNav(fullEphFilename);
end %end of function GetNasaHourlyEphemeris
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bOk,dirName] = CheckInputs(utcTime,dirName)
%% check we have the utcTime and right kind of dirName
if any(size(utcTime)~=[1,6])
error('utcTime must be a (1x6) vector')
end
if isempty(dirName), return, end
bOk = true;
if ~exist(dirName,'dir')
bOk = false;
fprintf('Error: directory ''%s'' not found\n',dirName);
if any(strfind(computer,'PCWIN')) && any(strfind(dirName,'/'))
fprintf('Looks like you''re using a PC, be sure your directory is specified with back-slashes ''\\''\n');
end
else
% check if there is a slash at the end of dirName
%decide what type of slash to use:
slashChar = '/'; %default
if any(strfind(dirName,'\'))
slashChar = '\';
end
if dirName(end)~=slashChar
dirName = [dirName,slashChar]; %add slash
end
end
end %end of function CheckInputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

34
opensource/GnssThresholds.m
Unescape View File

@ -0,0 +1,34 @@
classdef GnssThresholds
% GNSS validity thresholds we use in the code
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
properties (Constant) %listed alphabetically
MAXDELPOSFORNAVM = 20; %maximum position can change on one iteration of
% nav solution without los vector changing by more than 1 microradian
MAXPRRUNCMPS = 10; %max pseudorange rate (Doppler) uncertainty.
%bigger values may just be the search bin size, thus not valid for nav.
MAXTOWUNCNS = 500; %maximum Tow uncertainty, 500 ns. Satellite range
%can change by about half a millimeter in this time
MINNUMGPSEPH = 24; %minimum number of GPS ephemeris considered OK
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

148
opensource/Gps2Utc.m
Unescape View File

@ -0,0 +1,148 @@
function [utcTime] = Gps2Utc(gpsTime)
% [utcTime] = Gps2Utc(gpsTime)
% Convert GPS time (week & seconds) to UTC
%
% Input: gpsTime, [mx2] matrix [gpsWeek, gpsSeconds],
%
% Outputs: utcTime, [mx6] matrix = [year,month,day,hours,minutes,seconds]
%
% Valid range of inputs:
% gps times corresponding to 1980/6/1 <= time < 2100/1/1
% i.e. [0,0] <= gpsTime < [6260, 432000]
%
% See also: Utc2Gps
% Other functions needed: LeapSeconds.m
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%Algorithm for handling leap seconds:
% When a leap second happens, utc time stands still for one second, so
% gps seconds get further ahead, so we subtract leapSecs to move from gps time
%
% 1) convert gpsTime to time = [yyyy,mm,dd,hh,mm,ss] (with no leap seconds)
% 2) look up leap seconds for time: ls = LeapSeconds(time);
% This is (usually) the correct leap second value. Unless:
% If (utcTime=time) and (utcTime=time+ls) straddle a leap second
% then we need to add 1 to ls
% So, after step 2) ...
% 3) convert gpsTime-ls to timeMLs
% 4) look up leap seconds: ls1 = LeapSeconds(timeMLs);
% 5) if ls1~=ls, convert (gpsTime-ls1) to UTC Time
%% Check inputs
if size(gpsTime,2)~=2
error('gpsTime must have two columns')
end
fctSeconds = gpsTime*[GpsConstants.WEEKSEC; 1];
%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)
error('gpsTime must be in this range: [0,0] <= gpsTime < [6260, 432000]');
end
%% Finished checks
%from now on work with fct
%% Apply algorithm for handling leaps seconds
% 1) convert gpsTime to time = [yyyy,mm,dd,hh,mm,ss] (with no leap seconds)
time = Fct2Ymdhms(fctSeconds);
% 2) look up leap seconds for time: ls = LeapSeconds(time);
ls = LeapSeconds(time);
% 3) convert gpsTime-ls to timeMLs
timeMLs = Fct2Ymdhms(fctSeconds-ls);
% 4) look up leap seconds: ls1 = LeapSeconds(timeMLs);
ls1 = LeapSeconds(timeMLs);
% 5) if ls1~=ls, convert (gpsTime-ls1) to UTC Time
if all(ls1==ls)
utcTime = timeMLs;
else
utcTime = Fct2Ymdhms(fctSeconds-ls1);
end
% NOTE:
% Gps2Utc.m doesn't produce 23:59:60, at a leap second.
% Instead, as the leap second occurs, the Gps2Utc.m sequence of
% UTC hh:mm:ss is 23:59:59, 00:00:00, 00:00:00
% and we keep it like that for code simplicity.
% Here are a sequence of UTC and GPS times around a leap second:
% formalUtcTimes = [1981 12 31 23 59 59; 1981 12 31 23 59 60; 1982 1 1 0 0 0];
% gpsTimes = [103 431999; 103 432000; 103 432001];
% >> Gps2Utc(gpsTimes)
% ans =
% 1981 12 31 23 59 59
% 1982 1 1 0 0 0
% 1982 1 1 0 0 0
% If you want to change this you could check LeapSeconds.m to see if you're
% exactly on the addition of a leap second, and then change the UTC format
% to include the '60' seconds
end %end of function Gps2Utc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function time = Fct2Ymdhms(fctSeconds)
%Utility function for Gps2Utc
%Convert GPS full cycle time to [yyyy,mm,dd,hh,mm,ss.s] format
HOURSEC = 3600; MINSEC = 60;
monthDays=[31,28,31,30,31,30,31,31,30,31,30,31];%days each month (not leap year)
m=length(fctSeconds);
days = floor(fctSeconds / GpsConstants.DAYSEC) + 6;%days since 1980/1/1
years=zeros(m,1)+1980;
%decrement days by a year at a time, until we have calculated the year:
leap=ones(m,1); %1980 was a leap year
while (any(days > (leap+365)))
I = find(days > (leap+365) );
days(I) = days(I) - (leap(I) + 365);
years(I)=years(I)+1;
leap(I) = (rem(years(I),4) == 0); % leap = 1 on a leap year, 0 otherwise
% This works from 1901 till 2099, 2100 isn't a leap year (2000 is).
% Calculate the year, ie time(1)
end,
time=zeros(m,6);%initialize matrix
time(:,1)=years;
%decrement days by a month at a time, until we have calculated the month
% Calculate the month, ie time(:,2)
% Loop through m:
for i=1:m
month = 1;
if (rem(years(i),4) == 0) %This works from 1901 till 2099
monthDays(2)=29; %Make February have 29 days in a leap year
else
monthDays(2)=28;
end
while (days(i) > monthDays(month))
days(i) = days(i)-monthDays(month);
month = month+1;
end
time(i,2)=month;
end
time(:,3) = days;
sinceMidnightSeconds = rem(fctSeconds, GpsConstants.DAYSEC);
time(:,4) = fix(sinceMidnightSeconds/HOURSEC);
lastHourSeconds = rem(sinceMidnightSeconds, HOURSEC);
time(:,5) = fix(lastHourSeconds/MINSEC);
time(:,6) = rem(lastHourSeconds,MINSEC);
end %end of function Fct2Ymdhms
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

148
opensource/GpsAdrResiduals.m
Unescape View File

@ -0,0 +1,148 @@
function [adrResid]= GpsAdrResiduals(gnssMeas,allGpsEph,llaDegDegM)
% [adrResid]= GpsAdrResiduals(gnssMeas,allGpsEph,llaDegDegM)
% Compute residuals from GPS Accumulated Delta Ranges
%
% Inputs:
% gnssMeas.FctSeconds = Nx1 vector. Rx time tag of measurements.
% .ClkDCount = Nx1 vector. Hw clock discontinuity count
% .Svid = 1xM vector of all svIds found in gnssRaw.
% ...
% .tRxSeconds = NxM time of reception, seconds of gps week
% .tTxSeconds = NxM time of tranmission, seconds of gps week
% .AdrM = NxM accumulated delta range (= -k*carrier phase)
% ...
%
% allGpsEph, structure with all ephemeris
% llaDegDegM [1x3] true position
%
% Output:
% adrResid.FctSeconds = Nx1 time vector, same as gnssMeas.FctSeconds
% .Svid0 = reference satellite for single differences
% .Svid = 1xM vector of all svid
% .ResidM = [NxM] adr residuals
%
%Algorithm: compute single difference from sv to reference satellite svid0, then
% diff from reference time: tk - t0 (where t0 is the first common epoch for
% sv & svid0), then subtract expected values
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
adrResid = [];
if ~any(any(isfinite(gnssMeas.AdrM) & gnssMeas.AdrM~=0))
%Nothing in AdrM but NaNs and zeros
return
end
if nargin<3 || isempty(llaDegDegM)
fprintf('\nGpsAdrWlsPvt needs the true position: llaDegDegM')
return
end
xyz0M = Lla2Xyz(llaDegDegM);
M = length(gnssMeas.Svid);
N = length(gnssMeas.FctSeconds);
weekNum = floor(gnssMeas.FctSeconds/GpsConstants.WEEKSEC);
adrResid.FctSeconds = gnssMeas.FctSeconds;
adrResid.Svid0 = [];
adrResid.Svid = gnssMeas.Svid;
adrResid.ResidM = zeros(N,M)+NaN;
%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)
%choose Svid0 as the satellite that has most valid adr
numValidAdr = zeros(1,M);
for j=1:M
numValidAdr(j) = length(find(bitand(gnssMeas.AdrState(:,j),2^0)));
end
[~,j0] = max(numValidAdr);
adrResid.Svid0 = gnssMeas.Svid(j0);
svid = gnssMeas.Svid;
%% Compute expected pseudoranges
prHatM = zeros(N,M)+NaN; %to store expected pseudoranges
%"pseudo" here refers to the clock error in the satellite, not the receiver
%compute expected pr at each epoch
for i=1:N
for j=1:M
ttxSeconds = gnssMeas.tTxSeconds(i,j);
if isnan(ttxSeconds)
continue %skip to next
end
[gpsEph,iSv]= ClosestGpsEph(allGpsEph,svid(j),gnssMeas.FctSeconds(i));
if isempty(iSv)
continue; %skip to next
end
%compute pr for this sv
dtsv = GpsEph2Dtsv(gpsEph,ttxSeconds);
ttxSeconds = ttxSeconds - dtsv;%subtract dtsv from sv time to get true gps time
%calculate satellite position at ttx:
[svXyzTtxM,dtsv]=GpsEph2Xyz(gpsEph,[weekNum(i),ttxSeconds]);
%in ECEF coordinates at trx:
dtflightSeconds = norm(xyz0M - svXyzTtxM)/GpsConstants.LIGHTSPEED;
svXyzTrxM = FlightTimeCorrection(svXyzTtxM, dtflightSeconds);
prHatM(i,j) = norm(xyz0M - svXyzTrxM) - GpsConstants.LIGHTSPEED*dtsv;
% Use of dtsv: dtsv>0 <=> pr too small
end
end
%% Compute single difference, then delta from t0, then residuals
iT0 = zeros(1,M); %to store common start index for svid and svid0
for i=1:N %loop over time
if ~bitand(gnssMeas.AdrState(i,j0),2^0)
continue; %skip to next epoch
end
%now we have valid Adr for reference satellite, svid(j0) at this epoch
for j=1:M %loop over svid
if j==j0
continue; %skip over reference satellite
end
if ~bitand(gnssMeas.AdrState(i,j),2^0) %no valid adr for this svid
iT0(j) = 0;%reset common start index for this svid
continue %skip to next svid
elseif ~iT0(j) %we haven't yet set a common start,
if all(isfinite(prHatM(i,[j,j0]))) %we have valid prHat at this epoch
iT0(j) = i;%set common start index
end
end
%Now we have valid Adr for svid(j) and svid(j0)
i0 = iT0(j); %common start index for svid(j) and svid(j0)
if i>i0 && all(isfinite(prHatM(i,[j,j0])))%valid prHat at this epoch
delAdrM = (gnssMeas.AdrM(i,j) - gnssMeas.AdrM(i,j0)) - ...
(gnssMeas.AdrM(i0,j) - gnssMeas.AdrM(i0,j0));
delPrHatM = (prHatM(i,j) - prHatM(i,j0)) - ...
(prHatM(i0,j) - prHatM(i0,j0));
adrResid.ResidM(i,j) = delAdrM - delPrHatM;
%by all the above 'isfinite' checks, residM should not be NaN
assert(isfinite(adrResid.ResidM(i,j)),...
'Residual should be finite, check the above logic');
end
end
end
end %end of function GpsAdrResiduals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

41
opensource/GpsConstants.m
Unescape View File

@ -0,0 +1,41 @@
classdef GpsConstants
% GPS constants, defined by WGS84 and IS-GPS-200, and derived from those
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
properties (Constant) %listed alphabetically
EARTHECCEN2 = 6.69437999014e-3; %WGS 84 (Earth eccentricity)^2 (m^2)
EARTHMEANRADIUS = 6371009;%Mean R of ellipsoid(m) IU Gedosey& Geophysics
EARTHSEMIMAJOR = 6378137; %WGS 84 Earth semi-major axis (m)
EPHVALIDSECONDS = 7200; %+- 2 hours ephemeris validity
DAYSEC = 86400; %number of seconds in a day
FREL = -4.442807633e-10; %Clock relativity parameter, (s/m^1/2)
GPSEPOCHJD = 2444244.5; %GPS Epoch in Julian Days
HORIZDEG = 5; %angle above horizon at which GPS models break down
LIGHTSPEED = 2.99792458e8; %WGS-84 Speed of light in a vacuum (m/s)
%mean time of flight btwn closest GPS sat (~66 ms) & furthest (~84 ms):
MEANTFLIGHTSECONDS = 75e-3;
mu = 3.986005e14; %WGS-84 Universal gravitational parameter (m^3/sec^2)
WE = 7.2921151467e-5; % WGS 84 value of earth's rotation rate (rad/s)
WEEKSEC = 604800; %number of seconds in a week
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

89
opensource/GpsEph2Dtsv.m
Unescape View File

@ -0,0 +1,89 @@
function [dtsvS] = GpsEph2Dtsv(gpsEph,tS)
%[dtsvS] = GpsEph2Dtsv(gpsEph,tS)
%
% Calculate satellite clock bias
% Inputs:
% gpsEph = vector of GPS ephemeris structures defined in ReadRinexNav
% tS GPS time of week (secs) at time of transmission, calculate this as:
% tS = trx - PR/c,
% where trx = time of reception (receiver clock)
% PR = pseudorange
%
% tS may be a vector, gpsEph is a structure or vector of structures
% If gpsEph is a vector then t must have the same number of rows as gpsEph
% tS(i) is interpreted as being associated with gpsEph(i,:)
% If gpsEph is a single structure then tS may be any length
%
% Outputs:
% dtsvS = sat clock bias (seconds). GPS time = satellite time - dtsvS
% length(dtsvS) = length(tS)
%
% See IS-GPS-200 and RTCM Paper 99-92 SC104-88 for details of data
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%% Check size of inputs
if min(size(tS))>1
error('tS must be a vector or a scalar, not a matrix')
end
tS=tS(:)';%make t a row vector, to match [gpsEph.*] vectors below
pt=length(tS);
[p]=length(gpsEph);
if (p>1 && pt~=p),
error('If gpsEph is a vector tS must be a vector with #rows = length(gpsEph),\n')
end
%%
%%Extract the necessary variables from gpsEph
TGD = [gpsEph.TGD];
Toc = [gpsEph.Toc];
af2 = [gpsEph.af2];
af1 = [gpsEph.af1];
af0 = [gpsEph.af0];
Delta_n = [gpsEph.Delta_n];
M0 = [gpsEph.M0];
e = [gpsEph.e];
Asqrt = [gpsEph.Asqrt];
Toe = [gpsEph.Toe];
%%
%%Calculate dependent variables ------------------------------------------------
tk = tS - Toe; %time since time of applicability
I = find(tk > 302400.0);
if any(I), tk(I) = tk(I)-GpsConstants.WEEKSEC; end,
I = find(tk < -302400.0);
if (I), tk(I) = tk(I)+GpsConstants.WEEKSEC; end,
A = Asqrt.^2; %semi-major axis of orbit
n0=sqrt(GpsConstants.mu./(A.^3)); %Computed mean motion (rad/sec)
n=n0+Delta_n; %Corrected Mean Motion
Mk=M0+n.*tk; %Mean Anomaly
Ek=Kepler(Mk,e); %Solve Kepler's equation for eccentric anomaly
%%
%%Calculate satellite clock bias (See ICD-GPS-200 20.3.3.3.3.1) ----------------
dt = tS - Toc;
I = find(dt > 302400.0);
if any(I), dt(I) = dt(I)-GpsConstants.WEEKSEC; end,
I = find(dt < -302400.0);
if (I), dt(I) = dt(I)+GpsConstants.WEEKSEC; end,
dtsvS = af0 + af1.*dt + af2.*(dt.^2) + ...
GpsConstants.FREL.*e.*Asqrt.*sin(Ek) -TGD;
end % end of function GpsEph2Dtsv
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

66
opensource/GpsEph2Pvt.m
Unescape View File

@ -0,0 +1,66 @@
function [xM,dtsvS,vMps,dtsvSDot] = GpsEph2Pvt(gpsEph,gpsTime)
%[xM,dtsvS,vMps,dtsvSDot] = GpsEph2Pvt(gpsEph,gpsTime)
%
% Calculate sv coordinates, in ECEF frame, sv clock bias, and sv velocity
% Inputs:
% gpsEph: vector of ephemeris structures, as defined in ReadRinexEph.m
%
% gpsTime = [gpsWeek, ttxSec]: GPS time at time of transmission, ttx
% ttx = trx - PR/c - dtsvS, where trx is time of reception (receiver clock),
% dtsvS is the satellite clock error (seconds), can be computed in advance
% using eph2dtsv.m or iterating this function: gps time = sat time - dtsvS
% gpsWeek, ttxSec must be vectors of length(gpsEph),
%
% outputs:
% xM = [i,j,k] matrix of coordinates of satellites (ecef meters)
% dtsvS = vector of satellite clock error (seconds)
% vMps = [i,j,k] matrix of satellite velocity (ecef m/s)
% dtsvSDot = vector of satellite clock error rate (seconds/second)
% The row dimension of xM, dtsvS, vMps, dtsvSDot = length(gpsEph)
%
% xM and vMps are the satellite positions and velocities
% at time ttxSec, in terms of ecef coords at the same time
% Use FlightTimeCorrection.m to get xM & vMps in ecef coord at time of reception
%
% functions called: GpsEph2Xyz.m
%
% See IS-GPS-200 for details of data
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
vMps = []; dtsvSDot = [];
[xM,dtsvS]=GpsEph2Xyz(gpsEph,gpsTime);
if isempty(xM)
return
end
%compute velocity from delta position & dtsvS at (t+0.5) - (t-0.5)
%This is better than differentiating, because both the orbital and relativity
%terms have nonlinearities that are not easily differentiable
t1 = [gpsTime(:,1), gpsTime(:,2)+0.5]; %time + 0.5 seconds
[xPlus,dtsvPlus] = GpsEph2Xyz(gpsEph,t1);
t1 = [gpsTime(:,1), gpsTime(:,2)-0.5]; %time - 0.5 seconds
[xMinus,dtsvMinus]= GpsEph2Xyz(gpsEph,t1);
vMps = xPlus - xMinus;
dtsvSDot = dtsvPlus - dtsvMinus;
end %end of function GpsEph2Pvt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

152
opensource/GpsEph2Xyz.m
Unescape View File

@ -0,0 +1,152 @@
function [xyzM,dtsvS] = GpsEph2Xyz(gpsEph,gpsTime)
%[xyzM,dtsvS] = GpsEph2Xyz(gpsEph,gpsTime)
% Calculate sv coordinates, in ECEF frame, at ttx = gpsTime
%
% Inputs:
% gpsEph: vector of GPS ephemeris structures, as defined in ReadRinexNav.m
%
% gpsTime = [gpsWeek, ttxSec]: GPS time at time of transmission, ttx
% ttx = trx - PR/c - dtsvS, where trx is time of reception (receiver clock),
% dtsvS is the satellite clock error (seconds), can be computed in advance
% using Eph2Dtsv.m or iterating this function: gps time = sat time - dtsvS
% gpsWeek, ttxSec must be vectors of length(gpsEph),
%
% outputs:
% xyzM = [i,j,k] matrix of coordinates of satellites (ecef meters)
% dtsvS = vector of satellite clock error (seconds)
% The row dimension of xyzM, dtsvS = length(gpsEph)
%
% xyzM = sat positions at time ttxSec, in terms of ecef coords at the same time
% Use FlightTimeCorrection.m to get xyzM & v in ecef coords at time of reception
%
% See IS-GPS-200 for details of data
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
xyzM=[]; dtsvS=[];
[bOk,gpsEph,gpsWeek,ttxSec] = CheckGpsEphInputs(gpsEph,gpsTime);
if ~bOk
return
end
p=length(gpsEph);
%Now we are done checking and manipulating the inputs
% the time vectors: gpsWeek, ttxSec are the same length as gpsEph
%-------------------------------------------------------------------------------
%set fitIntervalSeconds
fitIntervalHours = [gpsEph.Fit_interval]';
%Rinex says "Zero if not known", so adjust for zeros
fitIntervalHours(fitIntervalHours == 0) = 2;
fitIntervalSeconds = fitIntervalHours*3600;
%Extract variables from gpsEph, into column vectors
%orbit variable names follow RINEX 2.1 Table A4
%clock variables af0, af1, af2 follow IS GPS 200
TGD = [gpsEph.TGD]';
Toc = [gpsEph.Toc]';
af2 = [gpsEph.af2]';
af1 = [gpsEph.af1]';
af0 = [gpsEph.af0]';
Crs = [gpsEph.Crs]';
Delta_n = [gpsEph.Delta_n]';
M0 = [gpsEph.M0]';
Cuc = [gpsEph.Cuc]';
e = [gpsEph.e]';
Cus = [gpsEph.Cus]';
Asqrt = [gpsEph.Asqrt]';
Toe = [gpsEph.Toe]';
Cic = [gpsEph.Cic]';
OMEGA = [gpsEph.OMEGA]';
Cis = [gpsEph.Cis]';
i0 = [gpsEph.i0]';
Crc = [gpsEph.Crc]';
omega = [gpsEph.omega]';
OMEGA_DOT=[gpsEph.OMEGA_DOT]';
IDOT = [gpsEph.IDOT]';
ephGpsWeek = [gpsEph.GPS_Week]';
%Calculate dependent variables -------------------------------------------------
%Time since time of applicability accounting for weeks and therefore rollovers
%subtract weeks first, to avoid precision errors:
tk =(gpsWeek-ephGpsWeek)*GpsConstants.WEEKSEC + (ttxSec-Toe);
I = find(abs(tk)>fitIntervalSeconds);
if ~isempty(I)
numTimes = length(I);
fprintf(sprintf('WARNING in GpsEph2Xyz.m, %d times outside fit interval.',...
numTimes));
end
A = Asqrt.^2; %semi-major axis of orbit
n0=sqrt(GpsConstants.mu./(A.^3)); %Computed mean motion (rad/sec)
n=n0+Delta_n; %Corrected Mean Motion
h = sqrt(A.*(1-e.^2).*GpsConstants.mu);
Mk=M0+n.*tk; %Mean Anomaly
Ek=Kepler(Mk,e); %Solve Kepler's equation for eccentric anomaly
%Calculate satellite clock bias (See ICD-GPS-200 20.3.3.3.3.1)
%subtract weeks first, to avoid precision errors:
dt =(gpsWeek-ephGpsWeek)*GpsConstants.WEEKSEC + (ttxSec-Toc);
%Calculate satellite clock bias
sin_Ek=sin(Ek);
cos_Ek=cos(Ek);
dtsvS = af0 + af1.*dt + af2.*(dt.^2) + ...
GpsConstants.FREL.*e.*Asqrt.*sin_Ek -TGD;
%true anomaly:
vk=atan2(sqrt(1-e.^2).*sin_Ek./(1-e.*cos_Ek),(cos_Ek-e)./(1-e.*cos_Ek));
Phik=vk + omega; %Argument of latitude
sin_2Phik=sin(2*Phik);
cos_2Phik=cos(2*Phik);
% The next three terms are the second harmonic perturbations
duk = Cus.*sin_2Phik + Cuc.*cos_2Phik; %Argument of latitude correction
drk = Crc.*cos_2Phik + Crs.*sin_2Phik; %Radius Correction
dik = Cic.*cos_2Phik + Cis.*sin_2Phik; %Correction to Inclination
uk = Phik + duk; %Corrected argument of latitude
rk = A.*((1-e.^2)./(1+e.*cos(vk))) + drk; %Corrected radius
ik = i0 + IDOT.*tk + dik; %Corrected inclination
sin_uk=sin(uk);
cos_uk=cos(uk);
xkp = rk.*cos_uk; % Position in orbital plane
ykp = rk.*sin_uk; % Position in orbital plane
% Wk = corrected longitude of ascending node,
Wk = OMEGA + (OMEGA_DOT - GpsConstants.WE).*tk - GpsConstants.WE*[gpsEph.Toe]';
%for dtflight, see FlightTimeCorrection.m
sin_Wk = sin(Wk);
cos_Wk = cos(Wk);
xyzM=zeros(p,3);
% The matrix xyzM contains the ECEF coordinates of position
sin_ik=sin(ik);
cos_ik=cos(ik);
xyzM(:,1)=xkp.*cos_Wk-ykp.*cos_ik.*sin_Wk;
xyzM(:,2)=xkp.*sin_Wk+ykp.*cos_ik.*cos_Wk;
xyzM(:,3)=ykp.*sin_ik;
end %end of function GpsEph2Xyz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

124
opensource/GpsWlsPvt.m
Unescape View File

@ -0,0 +1,124 @@
function gpsPvt = GpsWlsPvt(gnssMeas,allGpsEph,bRaw)
%gpsPvt = GpsWlsPvt(gnssMeas,allGpsEph,bRaw)
%compute PVT from gnssMeas
% Input: gnssMeas, structure of pseudoranges, etc. from ProcessGnssMeas
% allGpsEph, structure with all ephemeris
% [bRaw], default true, true => use raw pr, false => use smoothed
%
% Output:
% gpsPvt.FctSeconds Nx1 time vector, same as gnssMeas.FctSeconds
% .allLlaDegDegM Nx3 matrix, (i,:) = [lat (deg), lon (deg), alt (m)]
% .sigmaLlaM Nx3 standard deviation of [lat,lon,alt] (m)
% .allBcMeters Nx1 common bias computed with llaDegDegM
% .allVelMps Nx3 (i,:) = velocity in NED coords
% .sigmaVelMps Nx3 standard deviation of velocity (m/s)
% .allBcDotMps Nx1 common freq bias computed with velocity
% .numSvs Nx1 number of satellites used in corresponding llaDegDegM
% .hdop Nx1 hdop of corresponding fix
%
%Algorithm: Weighted Least Squares
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
if nargin<3
bRaw = true;
else
%check that smoothed pr fields exists in input
if any(~isfield(gnssMeas,{'PrSmM','PrSmSigmaM'}))
error('If bRaw is false, gnssMeas must have fields gnssMeas.PrSmM and gnssMeas.PrSmSigmaM')
end
end
xo =zeros(8,1);%initial state: [center of the Earth, bc=0, velocities = 0]'
weekNum = floor(gnssMeas.FctSeconds/GpsConstants.WEEKSEC);
N = length(gnssMeas.FctSeconds);
gpsPvt.FctSeconds = gnssMeas.FctSeconds;
gpsPvt.allLlaDegDegM = zeros(N,3)+NaN;
gpsPvt.sigmaLLaM = zeros(N,3)+NaN;
gpsPvt.allBcMeters = zeros(N,1)+NaN;
gpsPvt.allVelMps = zeros(N,3)+NaN;
gpsPvt.sigmaVelMps = zeros(N,3)+NaN;
gpsPvt.allBcDotMps = zeros(N,1)+NaN;
gpsPvt.numSvs = zeros(N,1);
gpsPvt.hdop = zeros(N,1)+inf;
for i=1:N
iValid = find(isfinite(gnssMeas.PrM(i,:))); %index into valid svid
svid = gnssMeas.Svid(iValid)';
[gpsEph,iSv] = ClosestGpsEph(allGpsEph,svid,gnssMeas.FctSeconds(i));
svid = svid(iSv); %svid for which we have ephemeris
numSvs = length(svid); %number of satellites this epoch
gpsPvt.numSvs(i) = numSvs;
if numSvs<4
continue;%skip to next epoch
end
%% WLS PVT -----------------------------------------------------------------
%for those svIds with valid ephemeris, pack prs matrix for WlsNav
prM = gnssMeas.PrM(i,iValid(iSv))';
prSigmaM= gnssMeas.PrSigmaM(i,iValid(iSv))';
prrMps = gnssMeas.PrrMps(i,iValid(iSv))';
prrSigmaMps = gnssMeas.PrrSigmaMps(i,iValid(iSv))';
tRx = [ones(numSvs,1)*weekNum(i),gnssMeas.tRxSeconds(i,iValid(iSv))'];
prs = [tRx, svid, prM, prSigmaM, prrMps, prrSigmaMps];
xo(5:7) = zeros(3,1); %initialize speed to zero
[xHat,~,~,H,Wpr,Wrr] = WlsPvt(prs,gpsEph,xo);%compute WLS solution
xo = xo + xHat;
%extract position states
llaDegDegM = Xyz2Lla(xo(1:3)');
gpsPvt.allLlaDegDegM(i,:) = llaDegDegM;
gpsPvt.allBcMeters(i) = xo(4);
%extract velocity states
RE2N = RotEcef2Ned(llaDegDegM(1),llaDegDegM(2));
%NOTE: in real-time code compute RE2N once until position changes
vNed = RE2N*xo(5:7); %velocity in NED
gpsPvt.allVelMps(i,:) = vNed;
gpsPvt.allBcDotMps(i) = xo(8);
%compute HDOP
H = [H(:,1:3)*RE2N', ones(numSvs,1)]; %observation matrix in NED
P = inv(H'*H);%unweighted covariance
gpsPvt.hdop(i) = sqrt(P(1,1)+P(2,2));
%compute variance of llaDegDegM
%inside LsPvt the weights are used like this:
% z = Hx, premultiply by W: Wz = WHx, and solve for x:
% x = pinv(Wpr*H)*Wpr*zPr;
% the point of the weights is to make sigma(Wz) = 1
% therefore, the variances of x come from diag(inv(H'Wpr'WprH))
P = inv(H'*(Wpr'*Wpr)*H); %weighted covariance
gpsPvt.sigmaLLaM(i,:) = sqrt(diag(P(1:3,1:3)));
%similarly, compute variance of velocity
P = inv(H'*(Wrr'*Wrr)*H); %weighted covariance
gpsPvt.sigmaVelMps(i,:) = sqrt(diag(P(1:3,1:3)));
%%end WLS PVT --------------------------------------------------------------
end
end %end of function GpsWlsPvt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

54
opensource/JulianDay.m
Unescape View File

@ -0,0 +1,54 @@
function jDays = JulianDay(utcTime)
% jDays = JulianDay(utcTime);
%
% input: utcTime [mx6] matrix [year,month,day,hours,minutes,seconds]
%
% output: totalDays in Julian Days [mx1] vector (real number of days)
%
% Valid input range: 1900 < year < 2100
%Algorithm from Meeus, (1991) Astronomical Algorithms,
%see http://www.geoastro.de/elevaz/basics/meeus.htm for online summary
% valid range 1900/3/1 to 2100/2/28
% but we limit inputs to 1901 through 2099, because it's simpler
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
% check inputs
if size(utcTime,2)~=6
error('utcTime must have 6 columns')
end
y = utcTime(:,1);
m = utcTime(:,2);
d = utcTime(:,3);
h = utcTime(:,4) + utcTime(:,5)/60 + utcTime(:,5)/3600;
%check that date is in valid range
if ( any(y<1901) || any (y>2099) )
error('utcTime(:,1) not in allowed range: 1900 < year < 2100')
end
i2 = m<=2; %index into months <=2
m(i2) = m(i2)+12;
y(i2) = y(i2)-1;
jDays = floor(365.25*y) + floor(30.6001*(m+1)) - 15 + 1720996.5 + d + h/24;
end %end of function JulianDay
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

59
opensource/Kepler.m
Unescape View File

@ -0,0 +1,59 @@
function ek = Kepler(mk,e)
% ek = Kepler(mk,e)
% Kepler - Solves Kepler's equation for ek through iteration.
%
% Inputs: mk: mean anomaly (rad)
% e: eccentricity
% Outputs: ek: eccentric anomaly
%
% NOTE: mk and e may be a vectors of the same dimensions, or one may be a scalar
% the output is a vector of the same dimensions as the input vector
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%Check inputs size
if min(size(mk))>1
error('mk must be a vector or a scalar, not a matrix')
end
if min(size(e))>1
error('e must be a vector or a scalar, not a matrix')
end
if length(mk)>1 && length(e)>1 && any(size(mk)~=size(e))
%both are vectors, they must have the same dimensions
error('If mk and e are both vectors they must have the same dimensions')
end
err = 1;
ek=mk;
iterCount = 0;
maxIterCount = 20;
while any(abs(err) > 1e-8) && iterCount < maxIterCount
err = ek - mk - (e.*sin(ek));
ek = ek - err;
iterCount = iterCount + 1;
if iterCount == maxIterCount
fprintf('Failed convergence on Kepler''s equation.\n')
end
end
end %end of function Kepler
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

202
opensource/LICENSE.txt
Unescape View File

@ -0,0 +1,202 @@
Apache License
Version 2.0, January 2004
http://www.apache.org/licenses/
TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
1. Definitions.
"License" shall mean the terms and conditions for use, reproduction,
and distribution as defined by Sections 1 through 9 of this document.
"Licensor" shall mean the copyright owner or entity authorized by
the copyright owner that is granting the License.
"Legal Entity" shall mean the union of the acting entity and all
other entities that control, are controlled by, or are under common
control with that entity. For the purposes of this definition,
"control" means (i) the power, direct or indirect, to cause the
direction or management of such entity, whether by contract or
otherwise, or (ii) ownership of fifty percent (50%) or more of the
outstanding shares, or (iii) beneficial ownership of such entity.
"You" (or "Your") shall mean an individual or Legal Entity
exercising permissions granted by this License.
"Source" form shall mean the preferred form for making modifications,
including but not limited to software source code, documentation
source, and configuration files.
"Object" form shall mean any form resulting from mechanical
transformation or translation of a Source form, including but
not limited to compiled object code, generated documentation,
and conversions to other media types.
"Work" shall mean the work of authorship, whether in Source or
Object form, made available under the License, as indicated by a
copyright notice that is included in or attached to the work
(an example is provided in the Appendix below).
"Derivative Works" shall mean any work, whether in Source or Object
form, that is based on (or derived from) the Work and for which the
editorial revisions, annotations, elaborations, or other modifications
represent, as a whole, an original work of authorship. For the purposes
of this License, Derivative Works shall not include works that remain
separable from, or merely link (or bind by name) to the interfaces of,
the Work and Derivative Works thereof.
"Contribution" shall mean any work of authorship, including
the original version of the Work and any modifications or additions
to that Work or Derivative Works thereof, that is intentionally
submitted to Licensor for inclusion in the Work by the copyright owner
or by an individual or Legal Entity authorized to submit on behalf of
the copyright owner. For the purposes of this definition, "submitted"
means any form of electronic, verbal, or written communication sent
to the Licensor or its representatives, including but not limited to
communication on electronic mailing lists, source code control systems,
and issue tracking systems that are managed by, or on behalf of, the
Licensor for the purpose of discussing and improving the Work, but
excluding communication that is conspicuously marked or otherwise
designated in writing by the copyright owner as "Not a Contribution."
"Contributor" shall mean Licensor and any individual or Legal Entity
on behalf of whom a Contribution has been received by Licensor and
subsequently incorporated within the Work.
2. Grant of Copyright License. Subject to the terms and conditions of
this License, each Contributor hereby grants to You a perpetual,
worldwide, non-exclusive, no-charge, royalty-free, irrevocable
copyright license to reproduce, prepare Derivative Works of,
publicly display, publicly perform, sublicense, and distribute the
Work and such Derivative Works in Source or Object form.
3. Grant of Patent License. Subject to the terms and conditions of
this License, each Contributor hereby grants to You a perpetual,
worldwide, non-exclusive, no-charge, royalty-free, irrevocable
(except as stated in this section) patent license to make, have made,
use, offer to sell, sell, import, and otherwise transfer the Work,
where such license applies only to those patent claims licensable
by such Contributor that are necessarily infringed by their
Contribution(s) alone or by combination of their Contribution(s)
with the Work to which such Contribution(s) was submitted. If You
institute patent litigation against any entity (including a
cross-claim or counterclaim in a lawsuit) alleging that the Work
or a Contribution incorporated within the Work constitutes direct
or contributory patent infringement, then any patent licenses
granted to You under this License for that Work shall terminate
as of the date such litigation is filed.
4. Redistribution. You may reproduce and distribute copies of the
Work or Derivative Works thereof in any medium, with or without
modifications, and in Source or Object form, provided that You
meet the following conditions:
(a) You must give any other recipients of the Work or
Derivative Works a copy of this License; and
(b) You must cause any modified files to carry prominent notices
stating that You changed the files; and
(c) You must retain, in the Source form of any Derivative Works
that You distribute, all copyright, patent, trademark, and
attribution notices from the Source form of the Work,
excluding those notices that do not pertain to any part of
the Derivative Works; and
(d) If the Work includes a "NOTICE" text file as part of its
distribution, then any Derivative Works that You distribute must
include a readable copy of the attribution notices contained
within such NOTICE file, excluding those notices that do not
pertain to any part of the Derivative Works, in at least one
of the following places: within a NOTICE text file distributed
as part of the Derivative Works; within the Source form or
documentation, if provided along with the Derivative Works; or,
within a display generated by the Derivative Works, if and
wherever such third-party notices normally appear. The contents
of the NOTICE file are for informational purposes only and
do not modify the License. You may add Your own attribution
notices within Derivative Works that You distribute, alongside
or as an addendum to the NOTICE text from the Work, provided
that such additional attribution notices cannot be construed
as modifying the License.
You may add Your own copyright statement to Your modifications and
may provide additional or different license terms and conditions
for use, reproduction, or distribution of Your modifications, or
for any such Derivative Works as a whole, provided Your use,
reproduction, and distribution of the Work otherwise complies with
the conditions stated in this License.
5. Submission of Contributions. Unless You explicitly state otherwise,
any Contribution intentionally submitted for inclusion in the Work
by You to the Licensor shall be under the terms and conditions of
this License, without any additional terms or conditions.
Notwithstanding the above, nothing herein shall supersede or modify
the terms of any separate license agreement you may have executed
with Licensor regarding such Contributions.
6. Trademarks. This License does not grant permission to use the trade
names, trademarks, service marks, or product names of the Licensor,
except as required for reasonable and customary use in describing the
origin of the Work and reproducing the content of the NOTICE file.
7. Disclaimer of Warranty. Unless required by applicable law or
agreed to in writing, Licensor provides the Work (and each
Contributor provides its Contributions) on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
implied, including, without limitation, any warranties or conditions
of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
PARTICULAR PURPOSE. You are solely responsible for determining the
appropriateness of using or redistributing the Work and assume any
risks associated with Your exercise of permissions under this License.
8. Limitation of Liability. In no event and under no legal theory,
whether in tort (including negligence), contract, or otherwise,
unless required by applicable law (such as deliberate and grossly
negligent acts) or agreed to in writing, shall any Contributor be
liable to You for damages, including any direct, indirect, special,
incidental, or consequential damages of any character arising as a
result of this License or out of the use or inability to use the
Work (including but not limited to damages for loss of goodwill,
work stoppage, computer failure or malfunction, or any and all
other commercial damages or losses), even if such Contributor
has been advised of the possibility of such damages.
9. Accepting Warranty or Additional Liability. While redistributing
the Work or Derivative Works thereof, You may choose to offer,
and charge a fee for, acceptance of support, warranty, indemnity,
or other liability obligations and/or rights consistent with this
License. However, in accepting such obligations, You may act only
on Your own behalf and on Your sole responsibility, not on behalf
of any other Contributor, and only if You agree to indemnify,
defend, and hold each Contributor harmless for any liability
incurred by, or claims asserted against, such Contributor by reason
of your accepting any such warranty or additional liability.
END OF TERMS AND CONDITIONS
APPENDIX: How to apply the Apache License to your work.
To apply the Apache License to your work, attach the following
boilerplate notice, with the fields enclosed by brackets "[]"
replaced with your own identifying information. (Don't include
the brackets!) The text should be enclosed in the appropriate
comment syntax for the file format. We also recommend that a
file or class name and description of purpose be included on the
same "printed page" as the copyright notice for easier
identification within third-party archives.
Copyright [yyyy] [name of copyright owner]
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.

94
opensource/LeapSeconds.m
Unescape View File

@ -0,0 +1,94 @@
function [leapSecs] = LeapSeconds(utcTime)
% leapSecs = LeapSeconds(utcTime)
% find the number of leap seconds since the GPS Epoch
%
% utcTime: [mx6] matrix
% utcTime(i,:) = [year,month,day,hours,minutes,seconds]
% year must be specified using four digits, e.g. 1994
% year valid range: 1980 <= year <= 2099
%
% Output: leapSecs,
% leapSecs(i) = number of leap seconds between the GPS Epoch and utcTime(i,:)
%
% The function looks up the number of leap seconds from a UTC time table
%
% LATEST LEAP SECOND IN THE TABLE = 31 Dec 2016.
% On 1 Jan 2017: GPS-UTC=18s
% See IERS Bulletin C, https://hpiers.obspm.fr/iers/bul/bulc/bulletinc.dat
% and http://tycho.usno.navy.mil/leapsec.html
%
% Aren't Leap Seconds a pain? Yes, and very costly. Ban the Leap Second:
% Leap seconds occur on average once every two years.
% What would happen if we had no leap seconds:
% 1) Thousand of engineers would NOT spend several weeks each two years fixing
% or planning for leap second bugs.
% 2) About seven thousand years from now, solar time would be 1 hour behind UTC
% and we would need a 1 hour adjustment - similar to daylight savings time.
% 3) GMT (which is solar time) will lose its significance.
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
[m,n] = size(utcTime);
if n~=6
error('utcTime input must have 6 columns');
end
% UTC table contains UTC times (in the form of [year,month,day,hours,mins,secs])
% At each of these times a leap second had just occurred
utcTable = [1982 1 1 0 0 0;
1982 7 1 0 0 0;
1983 7 1 0 0 0;
1985 7 1 0 0 0;
1988 1 1 0 0 0;
1990 1 1 0 0 0;
1991 1 1 0 0 0;
1992 7 1 0 0 0;
1993 7 1 0 0 0;
1994 7 1 0 0 0;
1996 1 1 0 0 0;
1997 7 1 0 0 0;
1999 1 1 0 0 0;
2006 1 1 0 0 0;
2009 1 1 0 0 0;
2012 7 1 0 0 0;
2015 7 1 0 0 0;
2017 1 1 0 0 0
];
%when a new leap second is announced in IERS Bulletin C
%update the table with the UTC time right after the new leap second
tableJDays = JulianDay(utcTable) - GpsConstants.GPSEPOCHJD; %days since GPS Epoch
%tableSeconds = tableJDays*GpsConstants.DAYSEC + utcTable(:,4:6)*[3600;60;1];
%NOTE: JulianDay returns a realed value number, corresponding to days and
% fractions thereof, so we multiply it by DAYSEC to get the full time in seconds
tableSeconds = tableJDays*GpsConstants.DAYSEC;
jDays = JulianDay(utcTime)- GpsConstants.GPSEPOCHJD; %days since GPS Epoch
%timeSeconds = jDays*GpsConstants.DAYSEC + utcTime(:,4:6)*[3600;60;1];
timeSeconds = jDays*GpsConstants.DAYSEC;
% tableSeconds and timeSeconds now contain number of seconds since the GPS epoch
leapSecs=zeros(m,1);
for i=1:m
%add up the number of leap seconds that have occured by timeSeconds(i)
leapSecs(i) = sum(tableSeconds<=timeSeconds(i));
end
end %end of function LeapSeconds
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

62
opensource/Lla2Ned.m
Unescape View File

@ -0,0 +1,62 @@
function [nedM] = Lla2Ned(lla1DegDegM,lla2DegDegM)
% [nedM] = Lla2Ned(lla1DegDegM,lla2DegDegM)
% function to difference latitude, longitude and altitude
% and provide an answer in NED coordinates in meters.
%
% Inputs: lla1DegDegM: mx3 matrix, [latitude(deg),longitude(deg),altitude(m)]
% lla2DegDegM: mx3 or 1x3, [latitude(deg),longitude(deg),altitude(m)]
% Output: nedM = lla1DegDegM - lla2DegDegM in NED coords (meters),
%
% Useful rules of thumb for quick conversions:
% 1e-5 (5th decimal place) of a degree of latitude approx= 1.1 meters
% 1e-5 of a degree of longitude approx= cos(latitude) * 1.1 meters
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
if nargin<2, error('Two inputs arguments needed'),end
[m1,n1]=size(lla1DegDegM);
[m2,n2]=size(lla2DegDegM);
if m2==1
lla2DegDegM=ones(m1,1)*lla2DegDegM;
else
if m1~=m2,
error('Second input must have one row or same number of rows as first input')
end
end
if n1~=3 || n2~=3, error('Both input matrices must have 3 columns'),end
[xyz1M] = Lla2Xyz(lla1DegDegM);
[xyz2M] = Lla2Xyz(lla2DegDegM);
refXyz = (xyz1M+xyz2M)/2;
[llaDegDegM] = Xyz2Lla(refXyz);
northM = zeros(m1,1);
eastM = zeros(m1,1);
for i=1:m1
Ce2n = RotEcef2Ned(llaDegDegM(i,1),llaDegDegM(1,2));
v = Ce2n*(xyz1M(i,:)-xyz2M(i,:))';
northM(i)=v(1);
eastM(i)=v(2);
end
downM = -lla1DegDegM(:,3)+lla2DegDegM(:,3);
nedM = [northM,eastM,downM];
end %end of function Lla2Ned
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

54
opensource/Lla2Xyz.m
Unescape View File

@ -0,0 +1,54 @@
function [xyzM] = Lla2Xyz(llaDegDegM)
% [xyzM] = Lla2Xyz(llaDegDegM)
% Transform latitude, longitude, altitude to ECEF coordinates.
%
% input: llaDegDegM = [mx3] matrix = [latDeg,lonDeg,altM]
% latitude, longitude are in degrees, altitude in meters
% output: xyzM = [mx3] matrix of ECEF coordinates in meters
%
% See also Xyz2Lla
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
% check inputs
if size(llaDegDegM,2)~=3
error('Input llaDegDegM must have three columns');
end
latDeg = llaDegDegM(:,1); lonDeg = llaDegDegM(:,2); altM = llaDegDegM(:,3);
%No rotation of longitude, by definition of ECEF
% Compute sines and cosines.
D2R = pi/180;
clat = cos(latDeg*D2R);
clon = cos(lonDeg*D2R);
slat = sin(latDeg*D2R);
slon = sin(lonDeg*D2R);
% Compute position vector in ECEF coordinates.
r0 = GpsConstants.EARTHSEMIMAJOR * ...
(sqrt(1.0 - GpsConstants.EARTHECCEN2 .* slat .* slat)).^(-1);
xM = (altM + r0) .* clat .* clon; % x coordinate
yM = (altM + r0) .* clat .* slon; % y coordinate
zM = (altM + r0 .* (1.0 - GpsConstants.EARTHECCEN2)).* slat;% z coordinate
[xyzM] = [xM,yM,zM];
end %end of function Lla2Xyz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

101
opensource/PlotAdr.m
Unescape View File

@ -0,0 +1,101 @@
function [colorsOut]= PlotAdr(gnssMeas,prFileName,colors)
% [colors] = PlotAdr(gnssMeas,[prFileName],[colors])
% plot the Accumulated Delta Ranges obtained from ProcessAdr
%
%gnssMeas.FctSeconds = Nx1 vector. Rx time tag of measurements.
% .ClkDCount = Nx1 vector. Hw clock discontinuity count
% .Svid = 1xM vector of all svIds found in gnssRaw.
% ...
% .AdrM = NxM accumulated delta range (= -k*carrier phase)
% .DelPrMinusAdrM = NxM DelPrM - AdrM
%
% Optional inputs: prFileName = string with file name
% colors, Mx3 color matrix
% if colors is not Mx3, it will be ignored
%
%Output: colors, color matrix, so we match colors each time we plot the same sv
if ~any(any(isfinite(gnssMeas.AdrM) & gnssMeas.AdrM~=0))
%Nothing in AdrM but NaNs and zeros
fprintf(' No ADR to plot\n'), return
end
if nargin<2
prFileName = '';
end
M = length(gnssMeas.Svid);
N = length(gnssMeas.FctSeconds);
if nargin<3 || any(size(colors)~=[M,3])
colors = zeros(M,3); %initialize color matrix for storing colors
bGotColors = false;
else
bGotColors = true;
end
gray = [.5 .5 .5];
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));
if ~isempty(iFi)
ti = timeSeconds(iFi(end));
h=plot(timeSeconds,AdrM); set(h,'Marker','.','MarkerSize',4)
if bGotColors
set(h,'Color',colors(j,:));
else
colors(j,:) = get(h,'Color');
end
text(ti,AdrM(iFi(end)),int2str(gnssMeas.Svid(j)),'Color',colors(j,:));
h123(2) = subplot(5,1,3:4); grid on,
h=plot(timeSeconds,gnssMeas.DelPrMinusAdrM(:,j)); hold on
set(h,'Marker','.','MarkerSize',4)
set(h,'Color',colors(j,:));
end
end
subplot(5,1,1:2); ax=axis;
title('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)')
bClockDis = [0;diff(gnssMeas.ClkDCount)~=0];%binary, 1 <=> clock discontinuity
%plot Clock discontinuity
h123(3) = subplot(5,1,5);
iCont = ~bClockDis;%index into where clock is continuous
plot(timeSeconds(iCont),bClockDis(iCont),'.b');%blue dots for continuous
hold on
plot(timeSeconds(~iCont),bClockDis(~iCont),'.r');%red dots for discontinuous
set(gca,'XLim',ax(1:2),'YLim',[-0.5 1.5]); grid on
set(gca,'YTick',[0 1],'YTickLabel',{'continuous ','discontinuous'})
set(gca,'YTickLabelRotation',45)
title('HW Clock Discontinuity')
xs = sprintf('time (seconds)\n%s',prFileName);
xlabel(xs,'Interpreter','none')
linkaxes(h123,'x');
if nargout>1
colorsOut = colors;
end
end %end of function PlotPseudorangeRates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

94
opensource/PlotAdrResids.m
Unescape View File

@ -0,0 +1,94 @@
function PlotAdrResids(adrResid,gnssMeas,prFileName,colors)
% PlotAdrResids(adrResid,gnssMeas,[prFileName],[colors])
% plot the Accumulated Delta Range residuals for the 5 sats with most valid adr
%
% adrResid.FctSeconds = Nx1 time vector, same as gnssMeas.FctSeconds
% .Svid0 = reference satellite for single differences
% .Svid = 1xM vector of all svid
% .ResidM = [NxM] adr residuals (adr = -k*carrier)
%
% gnssMeas, the measurements (from ProcessGnssMeas) used to get adrResids
%
% Optional inputs: prFileName = string with file name
% colors, Mx3 color matrix
% if colors is not Mx3, it will be ignored
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
K = 5; %number of satellites to plot
if ~any(any(isfinite(adrResid.ResidM))) %Nothing but NaNs
fprintf(' No adr residuals to plot\n'), return
end
if nargin<2
prFileName = '';
end
M = length(adrResid.Svid);
N = length(adrResid.FctSeconds);
if nargin<3 || any(size(colors)~=[M,3])
bGotColors = false;
else
bGotColors = true;
end
timeSeconds =adrResid.FctSeconds-adrResid.FctSeconds(1);%elapsed time in seconds
% find the K sats with most data
numValid = zeros(1,M);
for j=1:M
numValid(j) = sum(isfinite(adrResid.ResidM(:,j)));
end
[~,jSorted] = sort(numValid,'descend');
hK=zeros(1,K); %initialize plot handle
for k=1:K
hK(k) = subplot(K,1,k);
jSv = jSorted(k); %index into correct adrResid.Svid, and columns of .ResidM
svid = adrResid.Svid(jSv);
h = plot(timeSeconds,adrResid.ResidM(:,jSv)); grid on, hold on
j = find(gnssMeas.Svid == svid); %index into gnssMeas columns
if bGotColors
set(h,'Color',colors(j,:));
end
%get cycle slip flags
%From gps.h:
% #define GPS_ADR_STATE_CYCLE_SLIP (1<<2)
iCs = find(bitand(gnssMeas.AdrState(:,j),2^2));
numCs = length(iCs);
if numCs
h=plot(timeSeconds(iCs),zeros(numCs,1),'xk');
set(h,'MarkerSize',5)
end
ts=sprintf('Svids %d - %d',svid,adrResid.Svid0);
title(ts)
ylabel('(meters)');
end
xs = sprintf('time (seconds)\n%s',prFileName);
xlabel(xs,'Interpreter','none')
linkaxes(hK,'x');
subplot(K,1,1)
svid = adrResid.Svid(jSorted(1));
ts=sprintf('ADR single difference residuals. No iono or tropo correction. Svids %d - %d',...
svid,adrResid.Svid0);
title(ts)
ax=axis;
ht=text(ax(1),ax(3),' ''x'' = declared cycle slip');
set(ht,'VerticalAlignment','bottom','FontSize',8);
end %end of function PlotAdrResids
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

76
opensource/PlotCno.m
Unescape View File

@ -0,0 +1,76 @@
function [colorsOut] = PlotCno(gnssMeas,prFileName,colors)
%[colors] = PlotCno(gnssMeas,[prFileName],[colors])
% Plot the C/No from gnssMeas
%
%gnssMeas.FctSeconds = Nx1 vector. Rx time tag of measurements.
% .Svid = 1xM vector of all svIds found in gpsRaw.
% ...
% .Cn0DbHz = NxM
%
% Optional inputs: prFileName = string with file name
% colors, Mx3 color matrix
%
%Output: colors, color matrix, so we match colors each time we plot the same sv
M = length(gnssMeas.Svid);
N = length(gnssMeas.FctSeconds);
if nargin<2
prFileName = '';
end
if nargin<3 || any(size(colors)~=[M,3])
colors = zeros(M,3); %initialize color matrix for storing colors
bGotColors = false;
else
bGotColors = true;
end
timeSeconds = gnssMeas.FctSeconds-gnssMeas.FctSeconds(1);%elapsed time in seconds
%Plot C/No
for i=1:M
Cn0iDbHz = gnssMeas.Cn0DbHz(:,i);
iF = find(isfinite(Cn0iDbHz));
if ~isempty(iF)
ti = timeSeconds(iF(end));
h = plot(timeSeconds,Cn0iDbHz);
hold on
if bGotColors
set(h,'Color',colors(i,:));
else
colors(i,:) = get(h,'Color');
end
ts = int2str(gnssMeas.Svid(i));
if isfinite(gnssMeas.AzDeg(i))
ts = sprintf('%s, %03.0f^o, %02.0f^o',ts,...
gnssMeas.AzDeg(i),gnssMeas.ElDeg(i));
end
text(ti,Cn0iDbHz(iF(end)),ts,'Color',colors(i,:));
end
end
title('C/No in dB.Hz'),ylabel('(dB.Hz)')
xs = sprintf('time (seconds)\n%s',prFileName);
xlabel(xs,'Interpreter','none')
grid on
if nargout
colorsOut = colors;
end
end %end of function PlotCno
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

122
opensource/PlotPseudorangeRates.m
Unescape View File

@ -0,0 +1,122 @@
function [colorsOut]= PlotPseudorangeRates(gnssMeas,prFileName,colors)
% [colors] = PlotPseudorangeRates(gnssMeas,[prFileName],[colors])
% plot the Pseudorange Rates obtained from ProcessGnssMeas
%
%gnssMeas.FctSeconds = Nx1 vector. Rx time tag of measurements.
% .ClkDCount = Nx1 vector. Hw clock discontinuity count
% .Svid = 1xM vector of all svIds found in gnssRaw.
% ...
% .DelPrM = NxM change in pr while clock continuous
% .PrrMps = NxM pseudorange rate
% ...
%
% Optional inputs: prFileName = string with file name
% colors, Mx3 color matrix
% if colors is not Mx3, it will be ignored
%
%Output: colors, color matrix, so we match colors each time we plot the same sv
if nargin<2
prFileName = '';
end
M = length(gnssMeas.Svid);
N = length(gnssMeas.FctSeconds);
if nargin<3 || any(size(colors)~=[M,3])
colors = zeros(M,3); %initialize color matrix for storing colors
bGotColors = false;
else
bGotColors = true;
end
gray = [.5 .5 .5];
timeSeconds =gnssMeas.FctSeconds-gnssMeas.FctSeconds(1);%elapsed time in seconds
%plot slope of pr and prr
delPr = gnssMeas.DelPrM;
ts = ('diff(raw pr)/diff(time) and reported prr');
h12(1) = subplot(5,1,1:4); hold on
deltaMeanM = zeros(1,M)+NaN; %store mean diff(pr) - mean prr
for i=1:M
%plot prr
y = gnssMeas.PrrMps(:,i);
iFi = find(isfinite(y));
if ~isempty(iFi)
h = plot(timeSeconds,y); set(h,'Color',gray);
ti = timeSeconds(iFi(1));
ht=text(ti,y(iFi(1)),int2str(gnssMeas.Svid(i)),'Color',gray);
set(ht,'HorizontalAlignment','right')
meanPrrM = mean(y(iFi));%store for analysing delta prr dpr
%plot delta pr
y = delPr(:,i);
iFi = find(isfinite(y));
if any(iFi)
ti = timeSeconds(iFi(end));
y = diff(y)./diff(timeSeconds);%slope of pr (m/s)
h = plot(timeSeconds(2:end),y); set(h,'Marker','.','MarkerSize',4)
if bGotColors
set(h,'Color',colors(i,:));
else
colors(i,:) = get(h,'Color');
end
text(ti,y(iFi(end)-1),int2str(gnssMeas.Svid(i)),'Color',colors(i,:));
meanDprM = mean(y(isfinite(y)));%store for analysing delta prr dpr
deltaMeanM(i) = meanPrrM - meanDprM;
end
end
end
title(ts),ylabel('(m/s)')
ax = axis; %remember axis
yLimMps = [ax(3)-200, ax(4)];%make an extra 200 m/s on axis, to add text
set(gca,'YLim',yLimMps);
ts = ' For Svids [';
ds = sprintf('%.0f, ',gnssMeas.Svid);
ht=text(ax(1),ax(3),[ts,ds(1:end-2),'],']);
set(ht,'VerticalAlignment','top','Color',gray)
ts = ' mean(prr) - mean (diff(pr)/diff(time)) = [';
ds = sprintf('%.2f, ',deltaMeanM);
ht=text(ax(1),ax(3)-100,[ts,ds(1:end-2),'] (m/s)']);
set(ht,'VerticalAlignment','top','Color',gray)
bClockDis = [0;diff(gnssMeas.ClkDCount)~=0];%binary, 1 <=> clock discontinuity
%plot Clock discontinuity
h12(2) = subplot(5,1,5);
iCont = ~bClockDis;%index into where clock is continuous
plot(timeSeconds(iCont),bClockDis(iCont),'.b');%blue dots for continuous
hold on
plot(timeSeconds(~iCont),bClockDis(~iCont),'.r');%red dots for discontinuous
set(gca,'XLim',ax(1:2),'YLim',[-0.5 1.5]); grid on
set(gca,'YTick',[0 1],'YTickLabel',{'continuous ','discontinuous'})
set(gca,'YTickLabelRotation',45)
title('HW Clock Discontinuity')
xs = sprintf('time (seconds)\n%s',prFileName);
xlabel(xs,'Interpreter','none')
linkaxes(h12,'x')
if nargout
colorsOut = colors;
end
end %end of function PlotPseudorangeRates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

99
opensource/PlotPseudoranges.m
Unescape View File

@ -0,0 +1,99 @@
function [colorsOut] = PlotPseudoranges(gnssMeas,prFileName,colors)
%[colors] = PlotPseudoranges(gnssMeas,[prFileName],[colors])
% plot the Pseudoranges obtained from ProcessGnssMeas
%
%gnssMeas.FctSeconds = Nx1 vector. Rx time tag of measurements.
% .ClkDCount = Nx1 vector. Hw clock discontinuity count
% .Svid = 1xM vector of all svIds found in gnssRaw.
% ...
% .PrM = NxM matrix, row i corresponds to FctSeconds(i)
% ...
% .Cn0DbHz = NxM
%
% Optional inputs: prFileName = string with file name
% colors, Mx3 color matrix
% if colors is not Mx3, it will be ignored
%
%Output: colors, color matrix, so we match colors each time we plot the same sv
if nargin<2
prFileName = '';
end
M = length(gnssMeas.Svid);
N = length(gnssMeas.FctSeconds);
if nargin<3 || any(size(colors)~=[M,3])
colors = zeros(M,3); %initialize color matrix for storing colors
bGotColors = false;
else
bGotColors = true;
end
timeSeconds =gnssMeas.FctSeconds-gnssMeas.FctSeconds(1);%elapsed time in seconds
for i=1:M
%plot pseudoranges in meters
h123(1) = subplot(5,1,1:2); hold on,
priM = gnssMeas.PrM(:,i);
iF = find(isfinite(priM));
if ~isempty(iF)
ti = timeSeconds(iF(end));
h=plot(timeSeconds,priM); set(h,'Marker','.','MarkerSize',4)
if bGotColors
set(h,'Color',colors(i,:));
else
colors(i,:) = get(h,'Color');
end
text(ti,priM(iF(end)),int2str(gnssMeas.Svid(i)),'Color',colors(i,:));
%plot change in pseudoranges since first epoch
h123(2) = subplot(5,1,3:4); hold on
y = priM-priM(iF(1));
h=plot(timeSeconds,y);set(h,'Marker','.','MarkerSize',4)
set(h,'Color',colors(i,:));
text(ti,y(iF(end)),int2str(gnssMeas.Svid(i)),'Color',colors(i,:));
end
end
subplot(5,1,1:2); ax=axis;
title('Pseudoranges vs time'), ylabel('(meters)')
subplot(5,1,3:4); set(gca,'XLim',ax(1:2));
title('Pseudoranges change from initial value'),ylabel('(meters)')
bClockDis = [0;diff(gnssMeas.ClkDCount)~=0];%binary, 1 <=> clock discontinuity
%plot Clock discontinuity
h123(3) = subplot(5,1,5);
iCont = ~bClockDis;%index into where clock is continuous
plot(timeSeconds(iCont),bClockDis(iCont),'.b');%blue dots for continuous
hold on
plot(timeSeconds(~iCont),bClockDis(~iCont),'.r');%red dots for discontinuous
set(gca,'XLim',ax(1:2),'YLim',[-0.5 1.5]); grid on
set(gca,'YTick',[0 1],'YTickLabel',{'continuous ','discontinuous'})
set(gca,'YTickLabelRotation',45)
title('HW Clock Discontinuity')
xs = sprintf('time (seconds)\n%s',prFileName);
xlabel(xs,'Interpreter','none')
linkaxes(h123,'x');
if nargout
colorsOut = colors;
end
end %end of function PlotPseudoranges
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

125
opensource/PlotPvt.m
Unescape View File

@ -0,0 +1,125 @@
function PlotPvt(gpsPvt,prFileName,llaTrueDegDegM,titleString)
%PlotGpsPvt(gpsPvt,prFileName,[llaTrueDegDegM],[titleString])
%Plot the results of GpsLsPvt:
%
% gpsPvt.FctSeconds Nx1 time vector, same as gpsMeas.FctSeconds
% .allLlaDegDegM Nx3 matrix, (i,:) = [lat (deg), lon (deg), alt (m)]
% .sigmaLlaM Nx3 standard deviation of [lat,lon,alt] (m)
% .allBcMeters Nx1 common bias computed with lla
% .allVelMps Nx3 (i,:) = velocity in NED coords
% .sigmaVelMps Nx3 standard deviation of velocity (m/s)
% .allBcDotMps Nx1 common freq bias computed with velocity
% .numSvs Nx1 number of satellites used in corresponding lla
% .hdop Nx1 hdop of corresponding fix
%
% Optional inputs: [llaTrueDegDegM] = reference position, [ts] = title srtring
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
gray = [.5 .5 .5];
ltgray = [.8 .8 .8];
%set reference lla for plots
iFi = isfinite(gpsPvt.allLlaDegDegM(:,1));%index into finite results
if ~any(iFi)
return
end
llaMed = median(gpsPvt.allLlaDegDegM(iFi,:));%use median position as reference
if nargin < 3, llaTrueDegDegM = []; end
if nargin < 4. titleString = 'PVT solution'; end
bGotLlaTrue = ~isempty(llaTrueDegDegM) && any(llaTrueDegDegM);
%not empty and not all zeros
if bGotLlaTrue
llaRef = llaTrueDegDegM;
else
llaRef = llaMed;
end
%% plot ne errors vs llaTrueDegDegM --------------------------------------------
nedM = Lla2Ned(gpsPvt.allLlaDegDegM,llaRef);
h123=subplot(4,1,1:2);
h1 = plot(nedM(:,2),nedM(:,1));
set(h1,'LineStyle','-','LineWidth',0.1,'Color',ltgray)
hold on, plot(nedM(:,2),nedM(:,1),'cx');
lls = sprintf(' [%.6f^o, %.6f^o]',llaMed(1:2));
nedMedM = Lla2Ned(llaMed,llaRef);
h=plot(nedMedM(2),nedMedM(1),'+k','MarkerSize',18);
ts = [' median ',lls];
ht1 = text(nedMedM(2),nedMedM(1),ts,'color','k');
lls = sprintf(' [%.6f^o, %.6f^o]',llaRef(1:2));
if bGotLlaTrue
h=plot(0,0,'+r','MarkerSize',18);
ts = [' true pos ',lls];
ht2 = text(0,0,ts,'color','r');
%adjust VerticalAlignment to avoid overwriting previous text
if nedMedM(1)<0
set(ht1,'VerticalAlignment','top');%moves the 'median' label down
set(ht2,'VerticalAlignment','bottom');
else
set(ht1,'VerticalAlignment','bottom');
set(ht2,'VerticalAlignment','top');
end
%print median position error
ts = sprintf('|median - true pos| = %.1f m',norm(nedMedM(1:2)));
xm = min(0,nedMedM(2)); %align label with median x value, or zero
ym = min(nedM(:,1)); %align with smallest y
ht3 = text(xm,ym,ts);
set(ht3,'Color','k','VerticalAlignment','bottom');
end
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));
medM = median(distM);
%plot a circle using 'rectangle' (yes really :)
hr=rectangle('Position',[-1 -1 2 2]*medM,'Curvature',[1 1]);
set(hr,'EdgeColor',gray)
ts = sprintf('50%% distribution = %.1f m',medM);
text(0,medM,ts,'VerticalAlignment','bottom','Color',gray)
%% end of plot ne errors vs llaTrueDegDegM -------------------------------------
%time variable for plots
tSeconds = gpsPvt.FctSeconds-gpsPvt.FctSeconds(1);
%% plot speed
h123(2)=subplot(4,1,3);
N = length(gpsPvt.FctSeconds);
iGood = isfinite(gpsPvt.allVelMps(:,1));
speedMps = zeros(1,N)+NaN;
speedMps(iGood) = sqrt(sum(gpsPvt.allVelMps(iGood,1:2)'.^2)); %horizontal speed
plot(tSeconds,speedMps);grid on
ylabel('Horiz. speed (m/s)')
%% plot hdop & # sats
h123(3)=subplot(4,1,4);
[hyy,h1]=plotyy(tSeconds,gpsPvt.hdop,tSeconds,gpsPvt.numSvs,'plot','stairs');
grid on
set(h1,'LineWidth',1)
ylabel(hyy(1),'HDOP'); ylabel(hyy(2),'# sats');
xs = sprintf('time(seconds)\n%s',prFileName);
xlabel(xs,'Interpreter','none')
set(hyy,'Box','off')
linkaxes(h123(2:3),'x');
linkaxes(hyy,'x')
end %end of function PlotPvt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

150
opensource/PlotPvtStates.m
Unescape View File

@ -0,0 +1,150 @@
function PlotPvtStates(gnssPvt,prFileName)
%PlotPvtStates(gnssPvt,prFileName);
%Plot the Position, Velocity and Time/clock states in gnssPvt
%
% gnssPvt.FctSeconds Nx1 time vector, same as gpsMeas.FctSeconds
% .allLlaDegDegM Nx3 matrix, (i,:) = [lat (deg), lon (deg), alt (m)]
% .sigmaLlaM Nx3 standard deviation of [lat,lon,alt] (m)
% .allBcMeters Nx1 common bias computed with lla
% .allVelMps Nx3 (i,:) = velocity in NED coords
% .sigmaVelMps Nx3 standard deviation of velocity (m/s)
% .allBcDotMps Nx1 common freq bias computed with velocity
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%find median values
iFi = isfinite(gnssPvt.allLlaDegDegM(:,1));%index into finite results
if ~any(iFi)
return
end
llaMed = median(gnssPvt.allLlaDegDegM(iFi,:));%use medians as reference
%time variable for plots
tSeconds = gnssPvt.FctSeconds-gnssPvt.FctSeconds(1);
h1234 = zeros(1,4); %handles for subplots
%plot ned errors vs medians ----------------------------------------------------
ned = Lla2Ned(gnssPvt.allLlaDegDegM,llaMed);
h1234(1)=subplot(4,1,1);
plot(tSeconds,ned(:,1),'r');hold on
plot(tSeconds,ned(:,2),'g');
plot(tSeconds,ned(:,3),'b');
title('WLS: Position states offset from medians [Lat,Lon,Alt]');
grid on, ylabel('(meters)'),
%label Latitude, Longitude and Altitude
iFi = isfinite(ned(:,1));%index into finite results
h=zeros(1,3); %handles for Lat, Lon, Alt
h(1)=text(tSeconds(end),ned(iFi(end),1),'Lat');
set(h(1),'Color','r')
h(2)=text(tSeconds(end),ned(iFi(end),2),'Lon');
set(h(2),'Color','g')
h(3)=text(tSeconds(end),ned(iFi(end),3),'Alt');
set(h(3),'Color','b')
%shift the highest a little higher, so it doesnt overwrite the others
[~,iMax] = max(ned(iFi(end),:));
set(h(iMax),'VerticalAlignment','bottom');
%shift the lowest a little lower
[~,iMin] = min(ned(iFi(end),:));
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
if any(iFi)
plot(tSeconds,gnssPvt.allBcMeters - gnssPvt.allBcMeters(iFi(1)))
grid on
bc0Microseconds = gnssPvt.allBcMeters(iFi(1))/GpsConstants.LIGHTSPEED*1e6;
bc0Text = sprintf('%.1f',bc0Microseconds);
title(['Common bias ''clock'' offset from initial value of ',...
bc0Text,' {\mu}s'])
ylabel(('meters'))
end
%add microseconds label on right
ax=axis;
axMeters = ax(3:4);
set(gca,'YLim',axMeters); %fix this, else rescaling fig breaks axis proportion
axMicroseconds = axMeters/GpsConstants.LIGHTSPEED*1e6;
set(gca,'Box','off'); %# Turn off the box surrounding the whole axes
axesPosition = get(gca,'Position'); %# Get the current axes position
hNewAxes = axes('Position',axesPosition,... %# Place a new axes on top...
'Color','none',... %# ... with no background color
'YLim',axMicroseconds,... %# ... and a different scale
'YAxisLocation','right',... %# ... located on the right
'XTick',[],... %# ... with no x tick marks
'Box','off'); %# ... and no surrounding box
ylabel(hNewAxes,'(microseconds)'); %# Add label to right
%you must link the axes, else proportion will change when you scale figure:
linkaxes([hNewAxes, h1234(2)],'x');
%plot three components of speed ------------------------------------------------
h1234(3)=subplot(4,1,3);
vel = gnssPvt.allVelMps;
plot(tSeconds,vel(:,1),'r');hold on
plot(tSeconds,vel(:,2),'g');
plot(tSeconds,vel(:,3),'b');
title('Velocity states [North,East,Down]');
grid on, ylabel('(m/s)'),
%label North, East, Down
iFi = isfinite(vel(:,1));%index into finite results
h=zeros(1,3); %handles for Lat, Lon, Alt
h(1)=text(tSeconds(end),vel(iFi(end),1),'North');
set(h(1),'Color','r')
h(2)=text(tSeconds(end),vel(iFi(end),2),'East');
set(h(2),'Color','g')
h(3)=text(tSeconds(end),vel(iFi(end),3),'Down');
set(h(3),'Color','b')
%shift the highest a little higher, so it doesnt overwrite the others
[~,iMax] = max(vel(iFi(end),:));
set(h(iMax),'VerticalAlignment','bottom');
%shift the lowest a little lower
[~,iMin] = min(vel(iFi(end),:));
set(h(iMin),'VerticalAlignment','top');
%plot common frequency offset -------------------------------------------------
h1234(4)=subplot(4,1,4);
plot(tSeconds,gnssPvt.allBcDotMps)
grid on
title('Common frequency offset'); ylabel(('m/s'))
%add microseconds label on right
ax=axis;
axMps = ax(3:4);
set(gca,'YLim',axMps);%fix this, else rescaling fig breaks axis proportion
PPMPERMPS = 1/GpsConstants.LIGHTSPEED*1E6; %true for any frequency
axPpm = axMps*PPMPERMPS;
set(gca,'Box','off'); %# Turn off the box surrounding the whole axes
axesPosition = get(gca,'Position'); %# Get the current axes position
hNewAxes = axes('Position',axesPosition,... %# Place a new axes on top...
'Color','none',... %# ... with no background color
'YLim',axPpm,... %# ... and a different scale
'YAxisLocation','right',... %# ... located on the right
'XTick',[],... %# ... with no x tick marks
'Box','off'); %# ... and no surrounding box
ylabel(hNewAxes,'(ppm)'); %# Add label to right
linkaxes([hNewAxes, h1234(4)],'x')
xs = sprintf('\ntime(seconds)\n%s',prFileName);
xlabel(xs,'Interpreter','none')
linkaxes(h1234,'x'); %link all x axes
end %end of function PlotPvtStates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

81
opensource/ProcessAdr.m
Unescape View File

@ -0,0 +1,81 @@
function [gnssMeas]= ProcessAdr(gnssMeas)
% [gnssMeas]= ProcessAdr(gnssMeas)
% process the Accumulated Delta Ranges obtained from ProcessGnssMeas
%
%gnssMeas.FctSeconds = Nx1 vector. Rx time tag of measurements.
% .ClkDCount = Nx1 vector. Hw clock discontinuity count
% .Svid = 1xM vector of all svIds found in gnssRaw.
% ...
% .PrM = NxM pseudoranges, row i corresponds to FctSeconds(i)
% .DelPrM = NxM change in pr while clock continuous
% .AdrM = NxM accumulated delta range (= -k*carrier phase)
% ...
%
% output:
% gnssMeas.DelPrMinusAdrM = NxM DelPrM - AdrM, re-initialized to zero at each
% discontinuity or reset of DelPrM or AdrM
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
if ~any(any(isfinite(gnssMeas.AdrM) & gnssMeas.AdrM~=0))
%Nothing in AdrM but NaNs and zeros
fprintf(' No ADR recorded\n'), return
end
M = length(gnssMeas.Svid);
N = length(gnssMeas.FctSeconds);
DelPrMinusAdrM = zeros(N,M)+NaN;
for j=1:M %loop over Svid
AdrM = gnssMeas.AdrM(:,j); %make local variables for readability
DelPrM = gnssMeas.DelPrM(:,j);
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)
%keep valid values of AdrM only
iValid = bitand(AdrState,2^0);
iReset = bitand(AdrState,2^1);
AdrM(~iValid) = NaN;
%% work out DelPrM - AdrM since last discontinuity, plot DelPrM-AdrM
DelPrM0 = NaN; %to store initial offset from AdrM
for i=1:N %loop over time
if isfinite(AdrM(i)) && (AdrM(i)~=0) && isfinite(DelPrM(i)) && ...
~iReset(i)
%reinitialize after NaNs or AdrM zero or AdrState reset
if isnan(DelPrM0)
DelPrM0 = DelPrM(i) - AdrM(i);
end
else %reset at NaNs or AdrM zero
DelPrM0 = NaN;
end
DelPrMinusAdrM(i,j) = DelPrM(i) - DelPrM0 - AdrM(i);
end
end
gnssMeas.DelPrMinusAdrM = DelPrMinusAdrM;
end %end of function ProcessAdr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

244
opensource/ProcessGnssMeas.m
Unescape View File

@ -0,0 +1,244 @@
function gnssMeas = ProcessGnssMeas(gnssRaw)
% gnssMeas = ProcessGnssMeas(gnssRaw)
% Process raw measurements read from ReadGnssLogger
% Using technique explained in "Raw GNSS Measurements from Android" tutorial
%
% Input: gnssRaw, output from ReadGnssLogger
% Output: gnssMeas structure formatted conveniently for batch processing:
%gnssMeas.FctSeconds = Nx1 Full cycle time tag of M batched measurements.
% .ClkDCount = Nx1 Hw clock discontinuity count
% .HwDscDelS = Nx1 Hw clock change during each discontiuity (seconds)
% .Svid = 1xM all svIds found in gnssRaw.
% .AzDeg = 1xM azimuth in degrees at last valid epoch
% .ElDeg = 1xM elevation, ditto
% .tRxSeconds = NxM time of reception, seconds of gps week
% .tTxSeconds = NxM time of transmission, seconds of gps week
% .PrM = NxM pseudoranges, row i corresponds to FctSeconds(i)
% .PrSigmaM = NxM pseudorange error estimate (1-sigma)
% .DelPrM = NxM change in pr while clock continuous
% .PrrMps = NxM pseudorange rate
% .PrrSigmaMps= NxM
% .AdrM = NxM accumulated delta range (= -k * carrier phase)
% .AdrSigmaM = NxM
% .AdrState = NxM
% .Cn0DbHz = NxM
%
% all gnssMeas values are doubles
%
% Az and El are returned NaN. Compute Az,El using ephemeris and lla,
% or read from NMEA or GnssStatus
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
% Filter valid values first, so that rollover checks, etc, are on valid data
gnssRaw = FilterValid(gnssRaw);
%anything within 1ms is considered same epoch:
allRxMilliseconds = double(gnssRaw.allRxMillis);
gnssMeas.FctSeconds = (unique(allRxMilliseconds))*1e-3;
N = length(gnssMeas.FctSeconds);
gnssMeas.ClkDCount = zeros(N,1);
gnssMeas.HwDscDelS = zeros(N,1);
gnssMeas.Svid = unique(gnssRaw.Svid)'; %all the sv ids found in gnssRaw
M = length(gnssMeas.Svid);
gnssMeas.AzDeg = zeros(1,M)+NaN;
gnssMeas.ElDeg = zeros(1,M)+NaN;
gnssMeas.tRxSeconds = zeros(N,M)+NaN; %time of reception, seconds of gps week
gnssMeas.tTxSeconds = zeros(N,M)+NaN; %time of transmission, seconds of gps week
gnssMeas.PrM = zeros(N,M)+NaN;
gnssMeas.PrSigmaM = zeros(N,M)+NaN;
gnssMeas.DelPrM = zeros(N,M)+NaN;
gnssMeas.PrrMps = zeros(N,M)+NaN;
gnssMeas.PrrSigmaMps= zeros(N,M)+NaN;
gnssMeas.AdrM = zeros(N,M)+NaN;
gnssMeas.AdrSigmaM = zeros(N,M)+NaN;
gnssMeas.AdrState = zeros(N,M);
gnssMeas.Cn0DbHz = zeros(N,M)+NaN;
%GPS Week number:
weekNumber = floor(-double(gnssRaw.FullBiasNanos)*1e-9/GpsConstants.WEEKSEC);
%check for fields that are commonly all zero and may be missing from gnssRaw
if ~isfield(gnssRaw,'BiasNanos')
gnssRaw.BiasNanos = 0;
end
if ~isfield(gnssRaw,'TimeOffsetNanos')
gnssRaw.TimeOffsetNanos = 0;
end
%compute time of measurement relative to start of week
%subtract big longs (i.e. time from 1980) before casting time of week as double
WEEKNANOS = int64(GpsConstants.WEEKSEC*1e9);
weekNumberNanos = int64(weekNumber)*int64(GpsConstants.WEEKSEC*1e9);
%compute tRxNanos using gnssRaw.FullBiasNanos(1), so that
% tRxNanos includes rx clock drift since the first epoch:
tRxNanos = gnssRaw.TimeNanos -gnssRaw.FullBiasNanos(1) - weekNumberNanos;
%Assert if Tow state ~=1, because then gnssRaw.FullBiasNanos(1) might be wrong
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')
assert(all(tRxNanos >= 0),'tRxNanos should be >= 0')
%subtract the fractional offsets TimeOffsetNanos and BiasNanos:
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);
%we are ready to compute pseudorange in meters:
PrM = prSeconds*GpsConstants.LIGHTSPEED;
PrSigmaM = double(gnssRaw.ReceivedSvTimeUncertaintyNanos)*1e-9*...
GpsConstants.LIGHTSPEED;
PrrMps = gnssRaw.PseudorangeRateMetersPerSecond;
PrrSigmaMps = gnssRaw.PseudorangeRateUncertaintyMetersPerSecond;
AdrM = gnssRaw.AccumulatedDeltaRangeMeters;
AdrSigmaM = gnssRaw.AccumulatedDeltaRangeUncertaintyMeters;
AdrState = gnssRaw.AccumulatedDeltaRangeState;
Cn0DbHz = gnssRaw.Cn0DbHz;
%Now pack these vectors into the NxM matrices
for i=1:N %i is index into gnssMeas.FctSeconds and matrix rows
%get index of measurements within 1ms of this time tag
J = find(abs(gnssMeas.FctSeconds(i)*1e3 - allRxMilliseconds)<1);
for j=1:length(J) %J(j) is index into gnssRaw.*
k = find(gnssMeas.Svid==gnssRaw.Svid(J(j)));
%k is the index into gnssMeas.Svid and matrix columns
gnssMeas.tRxSeconds(i,k) = tRxSeconds(J(j));
gnssMeas.tTxSeconds(i,k) = tTxSeconds(J(j));
gnssMeas.PrM(i,k) = PrM(J(j));
gnssMeas.PrSigmaM(i,k) = PrSigmaM(J(j));
gnssMeas.PrrMps(i,k) = PrrMps(J(j));
gnssMeas.PrrSigmaMps(i,k)= PrrSigmaMps(J(j));
gnssMeas.AdrM(i,k) = AdrM(J(j));
gnssMeas.AdrSigmaM(i,k) = AdrSigmaM(J(j));
gnssMeas.AdrState(i,k) = AdrState(J(j));
gnssMeas.Cn0DbHz(i,k) = Cn0DbHz(J(j));
end
%save the hw clock discontinuity count for this epoch:
gnssMeas.ClkDCount(i) = gnssRaw.HardwareClockDiscontinuityCount(J(1));
if gnssRaw.HardwareClockDiscontinuityCount(J(1)) ~= ...
gnssRaw.HardwareClockDiscontinuityCount(J(end))
error('HardwareClockDiscontinuityCount changed within the same epoch');
end
end
gnssMeas = GetDelPr(gnssMeas);
end %of function ProcessGnssMeas
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gnssRaw = FilterValid(gnssRaw)
%utility function for ProcessGnssMeas,
%remove fields corresponding to measurements that are invalid
%NOTE: this makes it simpler to process data. But it removes data,
% so if you want to investigate *why* fields are invalid, then do so
% before calling this function
%check ReceivedSvTimeUncertaintyNanos, PseudorangeRateUncertaintyMetersPerSecond
%for now keep only Svid with towUnc<0.5 microseconds and prrUnc < 10 mps
iTowUnc = gnssRaw.ReceivedSvTimeUncertaintyNanos > GnssThresholds.MAXTOWUNCNS;
iPrrUnc = gnssRaw.PseudorangeRateUncertaintyMetersPerSecond > ...
GnssThresholds.MAXPRRUNCMPS;
iBad = iTowUnc | iPrrUnc;
if any(iBad)
numBad = sum(iBad);
%assert if we're about to remove everything:
assert(numBad<length(iBad),'Removing all measurements in gnssRaw')
names = fieldnames(gnssRaw);
for i=1:length(names)
ts = sprintf('gnssRaw.%s(iBad) = [];',names{i});
eval(ts); %remove fields for invalid meas
end
%explain to user what happened:
fprintf('\nRemoved %d bad meas inside ProcessGnssMeas>FilterValid because:\n',...
sum(iBad))
if any(iTowUnc)
fprintf('towUnc > %.0f ns\n',GnssThresholds.MAXTOWUNCNS)
end
if any(iPrrUnc)
fprintf('prrUnc > %.0f m/s\n',GnssThresholds.MAXPRRUNCMPS)
end
end
end %end of function FilterValid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gnssMeas = GetDelPr(gnssMeas)
% utility function for ProcessGnssMeas, compute DelPr
% gnssMeas.DelPrM = NxM, change in pr while clock is continuous
N = length(gnssMeas.FctSeconds);
M = length(gnssMeas.Svid);
bClockDis = [0;diff(gnssMeas.ClkDCount)~=0];%binary, 1 <=> clock discontinuity
%initialize first epoch to zero (by definition), rest to NaN
delPrM = zeros(N,M); delPrM(2:end,:) = NaN;
for j=1:M
i0=1; %i0 = index from where we compute DelPr
for i=2:N
if bClockDis(i) || isnan(gnssMeas.PrM(i0,j))
i0 = i; %reset to i if clock discont or a break in tracking
end
if bClockDis(i)
delPrM(i,j) = NaN;
else
delPrM(i,j) = gnssMeas.PrM(i,j) - gnssMeas.PrM(i0,j);
end
end
end
gnssMeas.DelPrM = delPrM;
end %of function GetDelPr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function prSeconds = CheckGpsWeekRollover(prSeconds)
%utility function for ProcessGnssMeas
iRollover = prSeconds > GpsConstants.WEEKSEC/2;
if any(iRollover)
fprintf('\nWARNING: week rollover detected in time tags. Adjusting ...')
prS = prSeconds(iRollover);
prS = prS - round(prS/GpsConstants.WEEKSEC)*GpsConstants.WEEKSEC;
%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')
else
prSeconds(iRollover) = prS; %put back into prSeconds vector
fprintf('\nCorrected week rollover')
end
end
%TBD Unit test this
end %end of function CheckGpsWeekRollover
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

84
opensource/ProcessGnssMeasScript.m
Unescape View File

@ -0,0 +1,84 @@
%ProcessGnssMeasScript.m, script to read GnssLogger output, compute and plot:
% pseudoranges, C/No, and weighted least squares PVT solution
%
% you can run the data in pseudoranges log files provided for you:
prFileName = 'pseudoranges_log_2016_06_30_21_26_07.txt'; %with duty cycling, no carrier phase
% prFileName = 'pseudoranges_log_2016_08_22_14_45_50.txt'; %no duty cycling, with carrier phase
% as follows
% 1) copy everything from GitHub google/gps-measurement-tools/ to
% a local directory on your machine
% 2) change 'dirName = ...' to match the local directory you are using:
dirName = '~/Documents/MATLAB/gpstools/opensource/demoFiles';
% 3) run ProcessGnssMeasScript.m script file
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%% data
%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';
%% parameters
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);
%% 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);
allGpsEph = GetNasaHourlyEphemeris(utcTime,dirName);
%% process raw measurements, compute pseudoranges:
[gnssMeas] = ProcessGnssMeas(gnssRaw);
%% plot pseudoranges and pseudorange rates
h1 = figure;
[colors] = PlotPseudoranges(gnssMeas,prFileName);
h2 = figure;
PlotPseudorangeRates(gnssMeas,prFileName,colors);
h3 = figure;
PlotCno(gnssMeas,prFileName,colors);
%% compute WLS position and velocity
gpsPvt = GpsWlsPvt(gnssMeas,allGpsEph);
%% plot Pvt results
h4 = figure;
ts = 'Raw Pseudoranges, Weighted Least Squares solution';
PlotPvt(gpsPvt,prFileName,param.llaTrueDegDegM,ts); drawnow;
h5 = figure;
PlotPvtStates(gpsPvt,prFileName);
%% Plot Accumulated Delta Range
if any(any(isfinite(gnssMeas.AdrM) & gnssMeas.AdrM~=0))
[gnssMeas]= ProcessAdr(gnssMeas);
h6 = figure;
PlotAdr(gnssMeas,prFileName,colors);
[adrResid]= GpsAdrResiduals(gnssMeas,allGpsEph,param.llaTrueDegDegM);drawnow
h7 = figure;
PlotAdrResids(adrResid,gnssMeas,prFileName,colors);
end
%% end of ProcessGnssMeasScript
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

94
opensource/README.md
Unescape View File

@ -0,0 +1,94 @@
# GPS Measurement Tools
The GNSS Measurement Tools code is provided for you to:
* read data from GnssLogger App,
* compute and visualize pseudoranges,
* compute weighted least squares position and velocity,
* view and analyze carrier phase (if it is present in the log file).
# Origin
This code is maintained on GitHub at the following link:
https://github.com/google/gps-measurement-tools
# Initial setup:
1. Extract the contents of the zip file to a directory, for example:
~/gpstools/*
and include the directory '~/gpstools/opensource' in your matlab path:
addpath('~/gpstools/opensource');
(Note: the tilde '~' is a place holder, don't actually use it, fill in
the actual complete path)
2. Edit ProcessGnssMeasScript.m to add the demoFiles directory, as follows:
dirName = '~/gpstools/opensource/demoFiles'
(again, replace tilde '~' with actual complete path)
3. Run ProcessGnssMeasScript.m, it will run with pre-recorded log files.
## To process a log file you collected from GnssLogger:
1. save the log file in a directory
2. edit ProcessGpsMeasScript.m, specify the file name and directory path
3. run ProcessGpsMeasScript.m
The code includes a function (GetNasaHourlyEphemeris.m) to read ephemeris
files from the NASA's archive of Space Geodesy Data, ftp://cddis.gsfc.nasa.gov
It will automatically go to the ftp when you have a new log file.
On some systems you need to use passive mode FTP; if this is required, see
The Mathworks site for how to do it.
Or (simpler): get the appropriate ephemeris file 'by hand' from the Nasa ftp
site (GetNasaHourlyEphemeris.m will tell you the correct url and filename),
copy the file to the directory where your log file is,
and GetNasaHourlyEphemeris.m will read it from there.
## For a summary of the open source GNSS Measurements Tools
See ~/gpstools/opensource/Contents.m or type 'help opensource' in matlab
command window.
# Platform specific notes:
For Windows: use '\' (backslash), instead of '/' for directories.
For Mac: when installing MATLAB.
System Preferences --> Security & Privacy -->
Allow Apps to be downloaded from: Mac App Store and identified developers
Uncompress/Unzip utility called from GetNasaHourlyEphemeris.m:
The ephemeris on the Nasa ftp is Unix-compressed. GetNasaHourlyEphemeris will
automatically uncompress it, if you have the right uncompress function on your
computer. If you need to install an unzip utility, see http://www.gpzip.org
Then search for 'uncompress' in the GetNasaHourlyEphemeris function to find and
edit the name of the unzip utility:
unzipCommand='uncompress';%edit if your platform uses something different
If you uncompress the file 'by hand' and rerun GetNasaHourlyEphemeris.m, it will
read the uncompressed file.
# Copyright Notice
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.

73
opensource/README.txt
Unescape View File

@ -0,0 +1,73 @@
The GNSS Measurement Tools code is provided for you to:
read data from GnssLogger App,
compute and visualize pseudoranges,
compute weighted least squares position and velocity,
view and analyze carrier phase (if it is present in the log file).
Initial setup:
1) Extract the contents of the zip file to a directory, for example:
~/gpstools/*
and include the directory '~/gpstools/opensource' in your matlab path:
addpath('~/gpstools/opensource');
(Note: the tilde '~' is a place holder, don't actually use it,
fill in the actual complete path)
2) Edit ProcessGnssMeasScript.m to add the demoFiles directory, as follows:
dirName = '~/gpstools/opensource/demoFiles'
(again, replace tilde '~' with actual complete path)
3) Run ProcessGnssMeasScript.m, it will run with pre-recorded log files.
To process a log file you collected from GnssLogger:
1) save the log file in a directory
2) edit ProcessGpsMeasScript.m, specify the file name and directory path
3) run ProcessGpsMeasScript.m
The code includes a function (GetNasaHourlyEphemeris.m) to read ephemeris
files from the NASA's archive of Space Geodesy Data, ftp://cddis.gsfc.nasa.gov
It will automatically go to the ftp when you have a new log file.
On some systems you need to use passive mode FTP; if this is required, see
The Mathworks site for how to do it.
Or (simpler): get the appropriate ephemeris file 'by hand' from the Nasa ftp
site (GetNasaHourlyEphemeris.m will tell you the correct url and filename),
copy the file to the directory where your log file is,
and GetNasaHourlyEphemeris.m will read it from there.
For a summary of the open source GNSS Measurements Tools,
see ~/gpstools/opensource/Contents.m
or type 'help opensource' in matlab command window
--------------------------------------------------------------------------------
Platform specific notes:
For Windows: use '\' (backslash), instead of '/' for directories.
For Mac: when installing MATLAB.
System Preferences --> Security & Privacy -->
Allow Apps to be downloaded from: Mac App Store and identified developers
Uncompress/Unzip utility called from GetNasaHourlyEphemeris.m:
The ephemeris on the Nasa ftp is Unix-compressed. GetNasaHourlyEphemeris will
automatically uncompress it, if you have the right uncompress function on your
computer. If you need to install an unzip utility, see http://www.gpzip.org
Then search for 'uncompress' in the GetNasaHourlyEphemeris function to find and
edit the name of the unzip utility:
unzipCommand='uncompress';%edit if your platform uses something different
If you uncompress the file 'by hand' and rerun GetNasaHourlyEphemeris.m, it will
read the uncompressed file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

451
opensource/ReadGnssLogger.m
Unescape View File

@ -0,0 +1,451 @@
function [gnssRaw,gnssAnalysis] = ReadGnssLogger(dirName,fileName,dataFilterIn,gnssAnalysis)
%% [gnssRaw,gnssAnalysis]=ReadGnssLogger(dirName,fileName,dataFilter,gnssAnalysis);
% Read the log file created by Gnss Logger App in Android
% Compatible with Android release N
%
% Input:
% 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'
%
% Output:
% gnssRaw, all GnssClock and GnssMeasurement fields from log file, including:
% .TimeNanos (int64)
% .FullBiasNanos (int64)
% ...
% .Svid
% .ReceivedSvTimeNanos (int64)
% .PseudorangeRateMetersPerSecond
% ...
% and data fields created by this function:
% .allRxMillis (int64), full cycle time of measurement (milliseconds)
% accurate to one millisecond, it is convenient for matching up time
% tags. For computing accurate location, etc, you must still use
% TimeNanos and gnssMeas.tRxSeconds
%
% gnssAnalysis, structure containing analysis, including list of missing fields
%
% see also: SetDataFilter, ProcessGnssMeas
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%factored into a few main sub-functions:
% MakeCsv()
% ReadRawCsv()
% FilterData()
% PackGnssRaw()
% CheckGnssClock()
% ReportMissingFields()
%% Initialize outputs and inputs
gnssAnalysis.GnssClockErrors = 'GnssClock Errors.';
gnssAnalysis.GnssMeasurementErrors = 'GnssMeasurement Errors.';
gnssAnalysis.ApiPassFail = '';
if nargin<3, dataFilterIn = []; end
%% check we have the right kind of fileName
extension = fileName(end-3:end);
if ~any(strcmp(extension,{'.txt','.csv'}))
error('Expecting file name of the form "*.txt", or "*.csv"');
end
%% read log file into a numeric matrix 'S', and a cell array 'header'
rawCsvFile = MakeCsv(dirName,fileName);
[header,C] = ReadRawCsv(rawCsvFile);
%% apply dataFilter
[dataFilter] = CheckDataFilter(dataFilterIn,header);
C = FilterData(C,dataFilter,header);
%% pack data into gnssRaw structure
[gnssRaw,missing] = PackGnssRaw(C,header);
%% check clock and measurements
[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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function csvFileName = MakeCsv(dirName,fileName)
%% make csv file, if necessary.
%And return extended csv file name (i.e. with full path in the name)
%TBD, maybe, do this entirely with Matlab read/write functions, make independent
%from grep and sed
%make extended file name
if dirName(end)~='/'
dirName = [dirName,'/']; %add /
end
csvFileName = [dirName,'prs.csv'];
if strcmp(fileName(end-3:end),'.csv')
return %input file is a csv file, nothing more to do here
end
extendedFileName = [dirName,fileName];
fprintf('\nReading file %s\n',extendedFileName)
%% read version
txtfileID = fopen(extendedFileName,'r');
if txtfileID<0
error('file ''%s'' not found',extendedFileName);
end
line='';
while isempty(strfind(lower(line),'version'))
line = fgetl(txtfileID);
if ~ischar(line) %eof or error occurred
if isempty(line)
error('\nError occurred while reading file %s\n',fileName)
end
break
end
end
if line==-1
fprintf('\nCould not find "Version" in input file %s\n',fileName)
return
end
%look for the beginning of the version number, e.g. 1.4.0.0
iDigits = regexp(line,'\d'); %index into the first number found in line
v = sscanf(line(iDigits(1):end),'%d.%d.%d.%d',4);
if length(v)<4
v(end+1:4,1)=0; %make v into a length 4 column vector
end
%Now extract the platform
k = strfind(line,'Platform:');
if any(k)
sPlatform = line(k+9:end);
else
sPlatform = '';%set empty if 'Platform:' not found
end
if isempty(strfind(sPlatform,'N'))
%add || strfind(platform,'O') and so on for future platforms
fprintf('\nThis version of ReadGnssLogger supports Android N\n')
error('Found "%s" in log file, expected "Platform: N"',line)
end
v1 = [1;4;0;0];
sCompare = CompareVersions(v,v1);
%Note, we need to check both the logger version (e.g. v1.0.0.0) and the
%Platform version "e.g. Platform: N" for any logic based on version
if strcmp(sCompare,'before')
fprintf('\nThis version of ReadGnssLogger supports v1.4.0.0 onwards\n')
error('Found "%s" in log file',line)
end
%% write csv file with header and numbers
%We could use grep and sed to make a csv file
%fclose(txtfileID);
% system(['grep -e ''Raw,'' ',extendedFileName,...
% ' | 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,
% note use of /g for "global" so sed acts on every occurrence in each line
% csv file "prs.csv" now contains a header row followed by numerical data
%
%But we'll do the same thing using Matlab, so people don't need grep/sed:
csvfileID = fopen(csvFileName,'w');
while ischar(line)
line = fgetl(txtfileID);
if isempty(strfind(line,'Raw,'))
continue %skip to next line
end
%Now 'line' contains the raw measurements header or data
line = strrep(line,'Raw,','');
line = strrep(line,'#',''); line = strrep(line,' ','');%remove '#' and spaces
%from versions v1.4.0.0 N we actually dont need to look for '(',')','true'
%or 'false' anymore. So we are done with replacing. That was easy.
fprintf(csvfileID,'%s\n',line);
end
fclose(txtfileID);
fclose(csvfileID);
if isempty(line) %line should be -1 at eof
error('\nError occurred while reading file %s\n',fileName)
end
end %end of function MakeCsv
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [header,C] = ReadRawCsv(rawCsvFile)
%% read data from csv file into a numerical matrix 'S' and cell array 'header'
S = csvread(rawCsvFile,1,0);%read numerical data from second row onwards
%Note csvread fills ,, with zero, so we will need a lower level read function to
%tell the difference between empty fields and valid zeros
%T = readtable(csvFileName,'FileType','text'); %use this to debug
%read header row:
fid = fopen(rawCsvFile);
if fid<0
error('file ''%s'' not found',rawCsvFile);
end
headerString = fgetl(fid);
if isempty(strfind(headerString,'TimeNanos'))
error('\n"TimeNanos" string not found in file %s\n',fileName)
end
header=textscan(headerString,'%s','Delimiter',',');
header = header{1}; %this makes header a numFieldsx1 cell array
numFields = size(header,1);
%check that numFields == size(S,2)
[~,M] = size(S); %M = number of columns
assert(numFields==M,...
'# of header names is different from # of columns of numerical data')
%read lines using formatSpec so we get TimeNanos and FullBiasNanos as
%int64, everything else as doubles, and empty values as NaN
formatSpec='';
for i=1:M
%lotsa || here, because we are comparing a vector, 'header'
%to a specific string each time. Not sure how to do this another way
%and still be able to easily read and debug. Better safe than too clever.
%longs
if i == find(strcmp(header,'TimeNanos')) || ...
i == find(strcmp(header,'FullBiasNanos')) || ...
i == find(strcmp(header,'ReceivedSvTimeNanos')) || ...
i == find(strcmp(header,'ReceivedSvTimeUncertaintyNanos')) || ...
i == find(strcmp(header,'CarrierCycles'))
formatSpec = sprintf('%s %%d64',formatSpec);
elseif 0
%ints
% TBD maybe %d32 for ints: AccumulatedDeltaRangeState, ...
% ConstellationType, MultipathIndicator, State, Svid
formatSpec = sprintf('%s %%d32',formatSpec);
else
%everything else doubles
formatSpec = sprintf('%s %%f',formatSpec);
end
end
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
%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:
% '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});
eval(ts);
%AND the iS index values on each iteration of i
iS = iS & iSi;
end
%Check if filter removes all values,
if ~any(iS) %if all zeros
fprintf('\nAll measurements removed. Specify dataFilter less strictly, ')
dataFilter(:,1)
return
end
%Now remove all values of C indexed by iS
for i=1:length(C)
C{i} = C{i}(iS);
end
end %end of function FilterDataS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [gnssRaw,missing] = PackGnssRaw(C,header)
%% pack data into gnssRaw, and report missing fields
assert(length(C)==length(header),...
'length(C) ~= length(header). This should have been checked before here')
gnssRaw = [];
%report clock fields present/missing, based on:
gnssClockFields = {...
'TimeNanos'
'TimeUncertaintyNanos'
'LeapSecond'
'FullBiasNanos'
'BiasUncertaintyNanos'
'DriftNanosPerSecond'
'DriftUncertaintyNanosPerSecond'
'HardwareClockDiscontinuityCount'
'BiasNanos'
};
missing.ClockFields = {};
%report measurements fields present/missing, based on:
gnssMeasurementFields = {...
'Cn0DbHz'
'ConstellationType'
'MultipathIndicator'
'PseudorangeRateMetersPerSecond'
'PseudorangeRateUncertaintyMetersPerSecond'
'ReceivedSvTimeNanos'
'ReceivedSvTimeUncertaintyNanos'
'State'
'Svid'
'AccumulatedDeltaRangeMeters'
'AccumulatedDeltaRangeUncertaintyMeters'
};
%leave these out for now, 'cause we dont care (for now), or they're deprecated,
% or they could legitimately be left out (because they are not computed in
% a particular GNSS implementation)
% SnrInDb, TimeOffsetNanos, CarrierFrequencyHz, CarrierCycles, CarrierPhase,
% CarrierPhaseUncertainty
missing.MeasurementFields = {};
%pack data into vector variables, if the fields are not NaNs
for j = 1:length(header)
if any(isfinite(C{j})) %not all NaNs
%TBD what if there are some NaNs, but not all. i.e. some missing
%data in the log file - TBD deal with this
eval(['gnssRaw.',header{j}, '=C{j};']);
elseif any(strcmp(header{j},gnssClockFields))
missing.ClockFields{end+1} = header{j};
elseif any(strcmp(header{j},gnssMeasurementFields))
missing.MeasurementFields{end+1} = header{j};
end
end
%TBD look for all zeros that can not legitimately be all zero,
%e.g. AccumulatedDeltaRangeMeters, and report these as missing data
end %end of function PackGnssRaw
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [gnssRaw,gnssAnalysis] = CheckGnssClock(gnssRaw,gnssAnalysis)
%% check clock values in gnssRaw
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');
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
end
if ~isfield(gnssRaw,'BiasNanos')
gnssRaw.BiasNanos = zeros(N,1);
end
if ~isfield(gnssRaw,'HardwareClockDiscontinuityCount')
gnssRaw.HardwareClockDiscontinuityCount = zeros(N,1);
fprintf('WARNING: Added HardwareClockDiscontinuityCount 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)
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)
end %end of function CheckGnssClock
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gnssAnalysis = ReportMissingFields(gnssAnalysis,missing)
%% report missing clock and measurement fields in gnssAnalysis
%report missing clock fields
if ~isempty(missing.ClockFields)
gnssAnalysis.GnssClockErrors = sprintf(...
'%s Missing Fields:',gnssAnalysis.GnssClockErrors);
for i=1:length(missing.ClockFields)
gnssAnalysis.GnssClockErrors = sprintf(...
'%s %s,',gnssAnalysis.GnssClockErrors,missing.ClockFields{i});
end
gnssAnalysis.GnssClockErrors(end) = '.';%replace final comma with period
end
%report missing measurement fields
if ~isempty(missing.MeasurementFields)
gnssAnalysis.GnssMeasurementErrors = sprintf(...
'%s Missing Fields:',gnssAnalysis.GnssMeasurementErrors);
for i=1:length(missing.MeasurementFields)
gnssAnalysis.GnssMeasurementErrors = sprintf(...
'%s %s,',gnssAnalysis.GnssMeasurementErrors,...
missing.MeasurementFields{i});
end
gnssAnalysis.GnssMeasurementErrors(end) = '.';%replace last comma with period
end
%assign pass/fail
if isempty(missing.ClockFields) && isempty(missing.MeasurementFields)
gnssAnalysis.ApiPassFail = 'PASS';
else
gnssAnalysis.ApiPassFail = 'FAIL BECAUSE OF MISSING FIELDS';
end
end %end of function ReportMissingFields
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

257
opensource/ReadRinexNav.m
Unescape View File

@ -0,0 +1,257 @@
function [gpsEph,iono] = ReadRinexNav(fileName)
% [gpsEph,iono] = ReadRinexNav(fileName)
%
% Read GPS ephemeris and iono data from an ASCII formatted RINEX 2.10 Nav file.
% Input:
% fileName - string containing name of RINEX formatted navigation data file
% Output:
% gpsEph: vector of ephemeris data, each element is an ephemeris structure
% structure order and orbit variable names follow RINEX 2.1 Table A4
% clock variable names af0, af1, af2 follow IS GPS 200
% gpsEph(i).PRN % SV PRN number
% gpsEph(i).Toc % Time of clock (seconds)
% gpsEph(i).af0 % SV clock bias (seconds)
% gpsEph(i).af1 % SV clock drift (sec/sec)
% gpsEph(i).af2 % SV clock drift rate (sec/sec2)
% gpsEph(i).IODE % Issue of data, ephemeris
% gpsEph(i).Crs % Sine harmonic correction to orbit radius (meters)
% gpsEph(i).Delta_n % Mean motion difference from computed value (radians/sec)
% gpsEph(i).M0 % Mean anomaly at reference time (radians)
% gpsEph(i).Cuc % Cosine harmonic correction to argument of lat (radians)
% gpsEph(i).e % Eccentricity (dimensionless)
% gpsEph(i).Cus % Sine harmonic correction to argument of latitude (radians)
% gpsEph(i).Asqrt % Square root of semi-major axis (meters^1/2)
% gpsEph(i).Toe % Reference time of ephemeris (seconds)
% gpsEph(i).Cic % Sine harmonic correction to angle of inclination (radians)
% gpsEph(i).OMEGA % Longitude of ascending node at weekly epoch (radians)
% gpsEph(i).Cis % Sine harmonic correction to angle of inclination (radians)
% gpsEph(i).i0 % Inclination angle at reference time (radians)
% gpsEph(i).Crc % Cosine harmonic correction to the orbit radius (meters)
% gpsEph(i).omega % Argument of perigee (radians)
% gpsEph(i).OMEGA_DOT% Rate of right ascension (radians/sec)
% gpsEph(i).IDOT % Rate of inclination angle (radians/sec)
% gpsEph(i).codeL2 % codes on L2 channel
% gpsEph(i).GPS_Week % GPS week (to go with Toe), (NOT Mod 1024)
% gpsEph(i).L2Pdata % L2 P data flag
% gpsEph(i).accuracy % SV user range accuracy (meters)
% gpsEph(i).health % Satellite health
% gpsEph(i).TGD % Group delay (seconds)
% gpsEph(i).IODC % Issue of Data, Clock
% gpsEph(i).ttx % Transmission time of message (seconds)
% gpsEph(i).Fit_interval %fit interval (hours), zero if not known
%
% iono: ionospheric parameter structure
% iono.alpha = [alpha0, alpha1, alpha2, alpha3]
% iono.beta = [ beta0, beta1, beta2, beta3]
% if iono data is not present in the Rinex file, iono is returned empty.
fidEph = fopen(fileName);
[numEph,numHdrLines] = countEph(fidEph);
%Now read from the begining again, looking for iono parameters
frewind(fidEph);
iono = readIono(fidEph,numHdrLines);
%initialize ephemeris structure array:
gpsEph = InitializeGpsEph;
gpsEph = repmat(gpsEph,1,numEph);
%now read each ephemeris into gpsEph(j)
%RINEX defines the format in terms of numbers of characters, so that's how we
%read it, e.g. "gpsEph(j).PRN = str2num(line(1:2));" and so on
for j = 1:numEph
line = fgetl(fidEph);
gpsEph(j).PRN = str2num(line(1:2));
%NOTE: we use str2num, not str2double, since str2num handles 'D' for exponent
%% get Toc (Rinex gives this as UTC time yy,mm,dd,hh,mm,ss)
year = str2num(line(3:6));
%convert year to a 4-digit year, this code is good to the year 2080.
%From 2080 RINEX 2.1 is ambiguous and shouldnt be used, because is has a
%2-digit year, and 100 years will have passed since the GPS Epoch.
if year < 80,
year = 2000+year;
else
year = 1900+year;
end
month = str2num(line(7:9));
day = str2num(line(10:12));
hour = str2num(line(13:15));
minute = str2num(line(16:18));
second = str2num(line(19:22));
%convert Toc to gpsTime
gpsTime = Utc2Gps([year,month,day,hour,minute,second]);
gpsEph(j).Toc = gpsTime(2);
%% get all other parameters
gpsEph(j).af0 = str2num(line(23:41));
gpsEph(j).af1 = str2num(line(42:60));
gpsEph(j).af2 = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).IODE = str2num(line(4:22));
gpsEph(j).Crs = str2num(line(23:41));
gpsEph(j).Delta_n = str2num(line(42:60));
gpsEph(j).M0 = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).Cuc = str2num(line(4:22));
gpsEph(j).e = str2num(line(23:41));
gpsEph(j).Cus = str2num(line(42:60));
gpsEph(j).Asqrt = str2num(line(61:79));
line=fgetl(fidEph);
gpsEph(j).Toe = str2num(line(4:22));
gpsEph(j).Cic = str2num(line(23:41));
gpsEph(j).OMEGA = str2num(line(42:60));
gpsEph(j).Cis = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).i0 = str2num(line(4:22));
gpsEph(j).Crc = str2num(line(23:41));
gpsEph(j).omega = str2num(line(42:60));
gpsEph(j).OMEGA_DOT = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).IDOT = str2num(line(4:22));
gpsEph(j).codeL2 = str2num(line(23:41));
gpsEph(j).GPS_Week = str2num(line(42:60));
gpsEph(j).L2Pdata = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).accuracy = str2num(line(4:22));
gpsEph(j).health = str2num(line(23:41));
gpsEph(j).TGD = str2num(line(42:60));
gpsEph(j).IODC = str2num(line(61:79));
line = fgetl(fidEph);
gpsEph(j).ttx = str2num(line(4:22));
gpsEph(j).Fit_interval = str2num(line(23:41));
end
fclose(fidEph);
end %end of function ReadRinexNav
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [numEph,numHdrLines] = countEph(fidEph,fileName)
%utility function for ReadRinexNav
%Read past the header, and then read to the end, counting ephemerides:
numHdrLines = 0;
bFoundHeader = false;
while ~bFoundHeader %Read past the header
numHdrLines = numHdrLines+1;
line = fgetl(fidEph);
if ~ischar(line), break, end
k = strfind(line,'END OF HEADER');
if ~isempty(k),
bFoundHeader = true;
break
end
end
if ~bFoundHeader
error('Error reading file: %s\nExpected RINEX header not found',fileName);
end
%Now read to the end of the file
numEph = -1;
while 1
numEph = numEph+1;
line = fgetl(fidEph);
if line == -1,
break;
elseif length(line)~=79
%use this opportunity to check line is the right length
%because in the rest of ReadRinexNav we depend on line being this length
error('Incorrect line length encountered in RINEX file');
end
end;
%check that we read the expected number of lines:
if rem(numEph,8)~=0
error('Number of nav lines in %s should be divisible by 8',fileName);
end
numEph = numEph/8; %8 lines per ephemeris
end %end of function countEph
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function iono = readIono(fidEph,numHdrLines)
%utility function to read thru the header lines, and find iono parameters
iono = []; %return empty if iono not found
bIonoAlpha=false; bIonoBeta=false;
for i = 1:numHdrLines,
line = fgetl(fidEph);
% Look for iono parameters, and read them in
if ~isempty(strfind(line,'ION ALPHA')) %line contains iono alpha parameters
ii = strfind(line,'ION ALPHA');
iono.alpha=str2num(line(1:ii-1));
%If we have 4 parameters then we have the complete iono alpha
bIonoAlpha = (length(iono.alpha)==4);
end
if ~isempty(strfind(line,'ION BETA'))%line contains the iono beta parameters
ii = strfind(line,'ION BETA');
iono.beta=str2num(line(1:ii-1));
%If we have 4 parameters then we have the complete iono beta
bIonoBeta = (length(iono.beta)==4);
end
end
if ~(bIonoAlpha && bIonoBeta)
%if we didn't get both alpha and beta iono correctly, then return empty iono
iono=[];
end
end %end of function readIono
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gpsEph = InitializeGpsEph
%utility function to initialize the ephemeris structure
gpsEph.PRN = 0;
gpsEph.Toc = 0;
gpsEph.af0 = 0;
gpsEph.af1 = 0;
gpsEph.af2 = 0;
gpsEph.IODE = 0;
gpsEph.Crs = 0;
gpsEph.Delta_n = 0;
gpsEph.M0 = 0;
gpsEph.Cuc = 0;
gpsEph.e = 0;
gpsEph.Cus = 0;
gpsEph.Asqrt = 0;
gpsEph.Toe = 0;
gpsEph.Cic = 0;
gpsEph.OMEGA = 0;
gpsEph.Cis = 0;
gpsEph.i0 = 0;
gpsEph.Crc = 0;
gpsEph.omega = 0;
gpsEph.OMEGA_DOT = 0;
gpsEph.IDOT = 0;
gpsEph.codeL2 = 0;
gpsEph.GPS_Week = 0;
gpsEph.L2Pdata = 0;
gpsEph.accuracy = 0;
gpsEph.health = 0;
gpsEph.TGD = 0;
gpsEph.IODC = 0;
gpsEph.ttx = 0;
gpsEph.Fit_interval= 0;
end %end of function InitializeEph
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

60
opensource/RotEcef2Ned.m
Unescape View File

@ -0,0 +1,60 @@
function Re2n = RotEcef2Ned(latDeg, lonDeg)
% Re2n = RotEcef2Ned(latDeg, lonDeg)
% Rotation matrix to convert an ECEF vector to
% North, East, Down coordinates, and vice-versa
%
% inputs: latDeg, lonDeg (degrees)
% output: Re2n, 3x3 unitary rotation matrix =
% [-sin(lat)*cos(lon), -sin(lat)*sin(lon), cos(lat);
% -sin(lon), cos(lon), 0 ;
% -cos(lat)*cos(lon), -cos(lat)*sin(lon),-sin(lat)]
%
% Example: vNed = Re2n*vEcef,
% Re2n'*vNed = vEcef
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
%CHECK INPUTS
if any(size(latDeg)~=[1,1]) || any(size(lonDeg)~=[1,1])
error('Inputs latDeg, lonDeg must be scalars')
end
D2R = pi/180; %degrees to radians scale factor
latRad=D2R*latDeg(:); lonRad=D2R*lonDeg(:);
clat = cos(latRad);
slat = sin(latRad);
clon = cos(lonRad);
slon = sin(lonRad);
Re2n = zeros(3,3);
Re2n(1,1) = -slat.*clon;
Re2n(1,2) = -slat.*slon;
Re2n(1,3) = clat;
Re2n(2,1) = -slon;
Re2n(2,2) = clon;
Re2n(2,3) = 0;
Re2n(3,1) = -clat.*clon;
Re2n(3,2) = -clat.*slon;
Re2n(3,3) = -slat;
end %end of function RotEcef2Ned
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

101
opensource/SetDataFilter.m
Unescape View File

@ -0,0 +1,101 @@
function 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
%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';
%you can create other filters in the same way ...
%for example, suppose you want to remove Svid 23:
% dataFilter{end+1} = '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';
% 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';
%limit to GPS only:
dataFilter{end+1} = 'ConstellationType==1';
%ConstellationType values are defined in Android HAL Documentation, gps.h,
% typedef uint8_t GnssConstellationType;
% #define GNSS_CONSTELLATION_UNKNOWN 0
% #define GNSS_CONSTELLATION_GPS 1
% #define GNSS_CONSTELLATION_SBAS 2
% #define GNSS_CONSTELLATION_GLONASS 3
% #define GNSS_CONSTELLATION_QZSS 4
% #define GNSS_CONSTELLATION_BEIDOU 5
% #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,
%so long as you dont use different heading types in one dataFilter{} entry
%bitwise data filters
%some fields are defined bitwise, including: State, AccumulatedDeltaRangeState
%
% GnssMeasurementState values are defined in Android HAL Documentation, gps.h,
% typedef uint32_t GnssMeasurementState;
% #define GNSS_MEASUREMENT_STATE_UNKNOWN 0
% #define GNSS_MEASUREMENT_STATE_CODE_LOCK (1<<0)
% #define GNSS_MEASUREMENT_STATE_BIT_SYNC (1<<1)
% #define GNSS_MEASUREMENT_STATE_SUBFRAME_SYNC (1<<2)
% #define GNSS_MEASUREMENT_STATE_TOW_DECODED (1<<3)
% #define GNSS_MEASUREMENT_STATE_MSEC_AMBIGUOUS (1<<4)
% #define GNSS_MEASUREMENT_STATE_SYMBOL_SYNC (1<<5)
% #define GNSS_MEASUREMENT_STATE_GLO_STRING_SYNC (1<<6)
% #define GNSS_MEASUREMENT_STATE_GLO_TOD_DECODED (1<<7)
% #define GNSS_MEASUREMENT_STATE_BDS_D2_BIT_SYNC (1<<8)
% #define GNSS_MEASUREMENT_STATE_BDS_D2_SUBFRAME_SYNC (1<<9)
% #define GNSS_MEASUREMENT_STATE_GAL_E1BC_CODE_LOCK (1<<10)
% #define GNSS_MEASUREMENT_STATE_GAL_E1C_2ND_CODE_LOCK (1<<11)
% #define GNSS_MEASUREMENT_STATE_GAL_E1B_PAGE_SYNC (1<<12)
% #define GNSS_MEASUREMENT_STATE_SBAS_SYNC (1<<13)
%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)';
% GNSS_MEASUREMENT_STATE_CODE_LOCK & GNSS_MEASUREMENT_STATE_TOW_DECODED
% GnssAccumulatedDeltaRangeState values are defined gps.h,
% typedef uint16_t GnssAccumulatedDeltaRangeState;
% #define GNSS_ADR_STATE_UNKNOWN 0
% #define GNSS_ADR_STATE_VALID (1<<0)
% #define GNSS_ADR_STATE_RESET (1<<1)
% #define GNSS_ADR_STATE_CYCLE_SLIP (1<<2)
%
%Example of how to use dataFilter for ADR State bit fields:
% get AccumulatedDeltaRangeState values that are valid
% dataFilter{end+1} = 'bitand(AccumulatedDeltaRangeState,1)';
end %end of function SetDataFilter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

111
opensource/Utc2Gps.m
Unescape View File

@ -0,0 +1,111 @@
function [gpsTime,fctSeconds] = Utc2Gps(utcTime)
% [gpsTime,fctSeconds] = Utc2Gps(utcTime)
% Convert the UTC date and time to GPS week & seconds
%
% Inputs:
% utcTime: [mx6] matrix
% utcTime(i,:) = [year,month,day,hours,minutes,seconds]
% year must be specified using four digits, e.g. 1994
% year valid range: 1980 <= year <= 2099
%
% Outputs:
% gpsTime: [mx2] matrix [gpsWeek, gpsSeconds],
% gpsWeek = number of weeks since GPS epoch
% gpsSeconds = number of seconds into gpsWeek,
% fctSeconds: full cycle time = seconds since GPS epoch (1980/01/06 00:00 UTC)
% Other functions needed: JulianDay.m, LeapSeconds.m
%initialize outputs
gpsTime=[];
fctSeconds=[];
[bOk]=CheckUtcTimeInputs(utcTime);
if ~bOk
return
end
HOURSEC = 3600; MINSEC = 60;
daysSinceEpoch = floor(JulianDay(utcTime) - GpsConstants.GPSEPOCHJD);
gpsWeek = fix(daysSinceEpoch/7);
dayofweek = rem(daysSinceEpoch,7);
% calculate the number of seconds since Sunday at midnight:
gpsSeconds = dayofweek*GpsConstants.DAYSEC + utcTime(:,4)*HOURSEC + ...
utcTime(:,5)*MINSEC + utcTime(:,6);
gpsWeek = gpsWeek + fix(gpsSeconds/GpsConstants.WEEKSEC);
gpsSeconds = rem(gpsSeconds,GpsConstants.WEEKSEC);
% now add leapseconds
leapSecs = LeapSeconds(utcTime);
fctSeconds = gpsWeek(:)*GpsConstants.WEEKSEC + gpsSeconds(:) + leapSecs(:);
% when a leap second happens, utc time stands still for one second,
% so gps seconds get further ahead, so we add leapsecs in going to gps time
gpsWeek = fix(fctSeconds/GpsConstants.WEEKSEC);
iZ = gpsWeek==0;
gpsSeconds(iZ) = fctSeconds(iZ); %set gpsSeconds directly, because rem(x,0)=NaN
gpsSeconds(~iZ) = rem(fctSeconds(~iZ),gpsWeek(~iZ)*GpsConstants.WEEKSEC);
gpsTime=[gpsWeek,gpsSeconds];
assert(all(fctSeconds==gpsWeek*GpsConstants.WEEKSEC+gpsSeconds),...
'Error in computing gpsWeek, gpsSeconds');
end %end of function Utc2Gps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bOk] = CheckUtcTimeInputs(utcTime)
%utility function for Utc2Gps
%check inputs
if size(utcTime,2)~=6
error('utcTime must have 6 columns')
end
%check that year, month & day are integers
x = utcTime(:,1:3);
if any(any( (x-fix(x)) ~= 0 ))
error('year,month & day must be integers')
end
%check that year is in valid range
if ( any(utcTime(:,1)<1980) || any(utcTime(:,1)>2099) )
error('year must have values in the range: [1980:2099]')
end
%check validity of month, day and time
if (any(utcTime(:,2))<1 || any(utcTime(:,2))>12 )
error('The month in utcTime must be a number in the set [1:12]')
end
if (any(utcTime(:,3)<1) || any(utcTime(:,3)>31))
error('The day in utcTime must be a number in the set [1:31]')
end
if (any(utcTime(:,4)<0) || any(utcTime(:,4)>=24))
error('The hour in utcTime must be in the range [0,24)')
end
if (any(utcTime(:,5)<0) || any(utcTime(:,5)>=60))
error('The minutes in utcTime must be in the range [0,60)')
end
if (any(utcTime(:,6)<0) || any(utcTime(:,6)>60))
%Note: seconds can equal 60 exactly, on the second a leap second is added
error('The seconds in utcTime must be in the range [0,60]')
end
bOk = true;
end %end of function CheckUtcTimeInputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

172
opensource/WlsPvt.m
Unescape View File

@ -0,0 +1,172 @@
function [xHat,z,svPos,H,Wpr,Wrr] = WlsPvt(prs,gpsEph,xo)
% [xHat,z,svPos,H,Wpr,Wrr] = WlsPvt(prs,gpsEph,xo)
% calculate a weighted least squares PVT solution, xHat
% given pseudoranges, pr rates, and initial state
%
% Inputs:
% prs: matrix of raw pseudoranges, and pr rates, each row of the form:
% [trxWeek,trxSeconds,sv,prMeters,prSigmaMeters,prrMps,prrSigmaMps]
% trxWeek, trxSeconds: Rx time of measurement
% where trxSeconds = seconds in the current GPS week
% sv: satellite id number
% prMeters, prSigmaMeters: pseudorange and standard deviation (meters)
% prrMps, prrSigmaMps: pseudorange rate and standard deviation (m/s)
% gpsEph: matching vector of GPS ephemeris struct, defined in ReadRinexNav
% xo: initial (previous) state, [x,y,z,bc,xDot,yDot,xDot,bcDot]'
% in ECEF coordinates(meters and m/s)
%
% Outputs: xHat: estimate of state update
% z = [zPr; zPrr] a-posteriori residuals (measured-calculated)
% svPos: matrix of calculated sv positions and sv clock error:
% [sv prn, x,y,z (ecef m), dtsv (s),xDot,yDot,zDot, dtsvDot]
% H: H observation matrix corresponding to svs in svPos(:,1)
% Wpr,Wrr Weights used in WlsPvt = 1/diag(sigma measurements)
% use these matrices to compute variances of xHat
%
% The PVT solution = xo + xHat, in ECEF coordinates
% For unweighted solution, set all sigmas = 1
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
jWk=1; jSec=2; jSv=3; jPr=4; jPrSig=5; jPrr=6; jPrrSig=7;%index of columns
[bOk,numVal] = checkInputs(prs, gpsEph, xo);
if ~bOk
error('inputs not right size, or not properly aligned with each other')
end
xHat=[]; z=[]; H=[]; svPos=[];
xyz0 = xo(1:3);
bc = xo(4);
if numVal<4
return
end
ttxWeek = prs(:,jWk); %week of tx. Note - we could get a rollover, when ttx_sv
%goes negative, and it is handled in GpsEph2Pvt, where we work with fct
ttxSeconds = prs(:,jSec) - prs(:,jPr)/GpsConstants.LIGHTSPEED; %ttx by sv clock
% this is accurate satellite time of tx, because we use actual pseudo-ranges
% here, not corrected ranges
% write the equation for pseudorange to see the rx clock error exactly cancel
% to get precise GPS time: we subtract the satellite clock error from sv time,
% as done next:
dtsv = GpsEph2Dtsv(gpsEph,ttxSeconds);
dtsv = dtsv(:); %make into a column for compatibility with other time vectors
ttx = ttxSeconds - dtsv; %subtract dtsv from sv time to get true gps time
%calculate satellite position at ttx
[svXyzTtx,dtsv,svXyzDot,dtsvDot]=GpsEph2Pvt(gpsEph,[ttxWeek,ttx]);
svXyzTrx = svXyzTtx; %initialize svXyz at time of reception
%Compute weights ---------------------------------------------------
Wpr = diag(1./prs(:,jPrSig));
Wrr = diag(1./prs(:,jPrrSig));
%iterate on this next part tilL change in pos & line of sight vectors converge
xHat=zeros(4,1);
dx=xHat+inf;
whileCount=0; maxWhileCount=100;
%we expect the while loop to converge in < 10 iterations, even with initial
%position on other side of the Earth (see Stanford course AA272C "Intro to GPS")
while norm(dx) > GnssThresholds.MAXDELPOSFORNAVM
whileCount=whileCount+1;
assert(whileCount < maxWhileCount,...
'while loop did not converge after %d iterations',whileCount);
for i=1:length(gpsEph)
% calculate tflight from, bc and dtsv
dtflight = (prs(i,jPr)-bc)/GpsConstants.LIGHTSPEED + dtsv(i);
% Use of bc: bc>0 <=> pr too big <=> tflight too big.
% i.e. trx = trxu - bc/GpsConstants.LIGHTSPEED
% Use of dtsv: dtsv>0 <=> pr too small <=> tflight too small.
% i.e ttx = ttxsv - dtsv
svXyzTrx(i,:) = FlightTimeCorrection(svXyzTtx(i,:), dtflight);
end
%calculate line of sight vectors and ranges from satellite to xo
v = xyz0(:)*ones(1,numVal,1) - svXyzTrx';%v(:,i) = vector from sv(i) to xyz0
range = sqrt( sum(v.^2) );
v = v./(ones(3,1)*range); % line of sight unit vectors from sv to xo
svPos=[prs(:,3),svXyzTrx,dtsv(:)];
%calculate the a-priori range residual
prHat = range(:) + bc -GpsConstants.LIGHTSPEED*dtsv;
% Use of bc: bc>0 <=> pr too big <=> rangehat too big
% Use of dtsv: dtsv>0 <=> pr too small
zPr = prs(:,jPr)-prHat;
H = [v', ones(numVal,1)]; % H matrix = [unit vector,1]
%z = Hx, premultiply by W: Wz = WHx, and solve for x:
dx = pinv(Wpr*H)*Wpr*zPr;
% update xo, xhat and bc
xHat=xHat+dx;
xyz0=xyz0(:)+dx(1:3);
bc=bc+dx(4);
%Now calculate the a-posteriori range residual
zPr = zPr-H*dx;
end
% Compute velocities ---------------------------------------------------------
rrMps = zeros(numVal,1);
for i=1:numVal
%range rate = [satellite velocity] dot product [los from xo to sv]
rrMps(i) = -svXyzDot(i,:)*v(:,i);
end
prrHat = rrMps + xo(8) - GpsConstants.LIGHTSPEED*dtsvDot;
zPrr = prs(:,jPrr)-prrHat;
%z = Hx, premultiply by W: Wz = WHx, and solve for x:
vHat = pinv(Wrr*H)*Wrr*zPrr;
xHat = [xHat;vHat];
z = [zPr;zPrr];
end %end of function WlsPvt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bOk,numVal] = checkInputs(prs, gpsEph, xo)
%utility function for WlsPvt
jWk=1; jSec=2; jSv=3; jPr=4; jPrSig=5; jPrr=6; jPrrSig=7;%index of columns
bOk=false;
%check inputs
numVal=size(prs,1);
if (max(prs(:,jSec))-min(prs(:,jSec)))> eps
return
elseif length(gpsEph)~=numVal
return
elseif any(prs(:,jSv) ~= [gpsEph.PRN]')
return
elseif any(size(xo) ~= [8,1])
return
elseif size(prs,2)~=7
return
else
bOk = true;
end
%We insist that gpsEph and prs are aligned first.
%ClosestGpsEph.m does this, and passes back indices for prs - this is the way to
%do it right, so we don't have nested searches for svId
end %end of function checkInputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

82
opensource/Xyz2Lla.m
Unescape View File

@ -0,0 +1,82 @@
function [llaDegDegM] = Xyz2Lla(xyzM)
% [llaDegDegM] = Xyz2Lla(xyzM)
%
% Transform Earth-Centered-Earth-Fixed x,y,z, coordinates to lat,lon,alt.
% Input: xyzM = [mx3] matrix of ECEF coordinates in meters
% Output: llaDegDegM = [mx3] matrix = [latDeg,lonDeg,altM]
% latitude, longitude are returned in degrees and altitude in meters
%
% See also Lla2Xyz
%Author: Frank van Diggelen
%Open Source code for processing Android GNSS Measurements
% check inputs
if size(xyzM,2)~=3
error('Input xyzM must have three columns');
end
% algorithm: Hoffman-Wellenhof, Lichtenegger & Collins "GPS Theory & Practice"
R2D = 180/pi;
%if x and y ecef positions are both zero then lla is undefined
iZero = ( xyzM(:,1)==0 & xyzM(:,2)==0);
xyzM(iZero,:) = NaN; %set to NaN, so lla will also be NaN
xM = xyzM(:,1); yM = xyzM(:,2); zM = xyzM(:,3);
%following algorithm from Hoffman-Wellenhof, et al. "GPS Theory & Practice":
a = GpsConstants.EARTHSEMIMAJOR;
a2 = a^2;
b2 = a2*(1-GpsConstants.EARTHECCEN2);
b = sqrt(b2);
e2 = GpsConstants.EARTHECCEN2;
ep2 = (a2-b2)/b2;
p=sqrt(xM.^2 + yM.^2);
% two sides and hypotenuse of right angle triangle with one angle = theta:
s1 = zM*a;
s2 = p*b;
h = sqrt(s1.^2 + s2.^2);
sin_theta = s1./h;
cos_theta = s2./h;
%theta = atan(s1./s2);
% two sides and hypotenuse of right angle triangle with one angle = lat:
s1 = zM+ep2*b*(sin_theta.^3);
s2 = p-a*e2.*(cos_theta.^3);
h = sqrt(s1.^2 + s2.^2);
tan_lat = s1./s2;
sin_lat = s1./h;
cos_lat = s2./h;
latDeg = atan(tan_lat);
latDeg = latDeg*R2D;
N = a2*((a2*(cos_lat.^2) + b2*(sin_lat.^2)).^(-0.5));
altM = p./cos_lat - N;
% rotate longitude to where it would be for a fixed point in ECI
%lonDeg = atan2(yM,xM) - GpsConstants.WE*deltaTime;
lonDeg = atan2(yM,xM); %since deltaTime = 0 for ECEF
lonDeg = rem(lonDeg,2*pi)*R2D;
if (lonDeg>180)
lonDeg = lonDeg-360;
end
llaDegDegM = [latDeg, lonDeg,altM];
end %end of function Xyz2Lla
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

3352
opensource/demoFiles/hour1820.16n
View File

File diff suppressed because it is too large Load Diff

3360
opensource/demoFiles/hour2350.16n
View File

File diff suppressed because it is too large Load Diff

1380
opensource/demoFiles/prs.csv
View File

File diff suppressed because it is too large Load Diff

1384
opensource/demoFiles/prsWithPrInMeters.csv
View File

File diff suppressed because it is too large Load Diff

1607
opensource/demoFiles/pseudoranges_log_2016_06_30_21_26_07.txt
View File

File diff suppressed because it is too large Load Diff

5259
opensource/demoFiles/pseudoranges_log_2016_08_22_14_45_50.txt
View File

File diff suppressed because it is too large Load Diff
Write Preview
Loading…
Cancel
Save