继续来看伪距测量函数postNavigation.m函数和plotNavigation.m，其中伪距测量函数又包含了位置计算函数和伪距计算函数（leastSquarePos.m和calculatePseudoranges.m）
在先前运行postProcessing.m时，每当运行到postNavigation就会报错，那么今天就来寻找一下报错的原因吧~
目录
calculatePseudoranges()函数_伪距计算函数
leastSquarePos.m函数_使用最小二乘法计算接收机位置
postNavigation.m函数
该函数功能：生成接收机的导航解决方案(伪距，位置)。最后，它将WGS84系统的坐标转换为UTM、地心或任何其他坐标系统。
function [navSolutions, eph] = postNavigation(trackResults, settings)输入输出的参数：
Inputs:
trackResults - 跟踪函数(结构数组)的结果
settings - receiver settings.
Outputs:
navSolutions - 包含测量的伪距，接收机时钟误差，多个坐标系中的接收机坐标(至少是ECEF和UTM)。
eph - 接收到所有SV的星历表(结构数组)。
纵观该函数，发现计算位置这个函数调用了很多函数，其中包括但不限于：
findPreambles函数（功能：寻找遥测字和交接字）
ephemeris函数（功能：解码星历表和第一子帧的周内时）
calculatePseudoranges函数（功能：计算伪距）
satpos函数（功能：找到卫星位置并完成时钟校正）
leastSquarePos函数（功能：计算接收机位置）
cart2geo函数（功能：转换为大地坐标）
findUtmZone函数
cart2utm函数（功能：转换为UTM坐标系）
为了方便理解，本文先系统的说明postNavigation.m函数，然后对涉及到的各个函数依次说明，在有涉及相关函数时，会附带链接。
postNavigation.m函数从确定位转变和帧头位置开始，得到每位的值， 并解码获得星历，包括子帧1、子帧2和子帧3。子帧4和子帧5里的信息也可能包括。
接下来该函数调用伪距测量函数并计算位置坐标。伪距和位置的计算覆盖接收机sttings结构变量中描述的指定时间段。
postNavigation函数读取以下变量:
(1) navSolPeriod :计算伪距和位置的频率。
(2) elevationMask: 卫星屏蔽角,设置用于位置计算的卫星的最小仰角。低仰角卫星信号其大气误差大。
(3) UTMzone: UTM时区用来转换坐标( ECEF到UTM) ,这是一个与接收机位置有关的整数。
(4) truePosition: 如果接收机天线的精确位置已知，则可指定天线的(东、北、天)坐标,此坐标从软件接收机计算的结果中提取,并画出其图。
也可输人(E .N、U)做为近似坐标或零点。
postNavigation.m函数主要包含以下几个部分：
① 判断是否符合进行位置计算的条件
② 寻找帧头位置的开始点
③ 解码星历表
④ 检测卫星的数量是否一直在3颗以上并初始化
⑤ 计算伪距
⑥ 找到卫星位置并完成时钟校正
⑦ 计算接收机位置
⑧ 坐标变换
以下分别对这8个部分一一介绍，设涉及到的函数，为不影响postNavigation.m函数的连续性，会有跳转链接，点击跳转即可参考。
① 判断是否符合进行位置计算的条件
必须至少有三个子帧(编号1、2和3)才能找到卫星坐标。然后也可以找到接收器的位置。因为跟踪开始于任意点，因此该函数需要所有5个子帧的数据。
一个子帧长度是6秒，因此我们需要至少30秒长的记录(5 * 6 = 30秒= 30000ms)。当跟踪在一个子帧的中间开始时，还需要为这种情况添加额外的秒数。
% 判断是否符合执行条件
if (settings.msToProcess < 36000) || (sum([trackResults.status] ~= '-') < 4)
% Show the error message and exit
disp('Record is to short or too few satellites tracked. Exiting!');
navSolutions = [];
eph = [];
return
end② 寻找帧头位置的开始点
GPS导航数据解码中的首要问题是确定子帧的起始位置，该语句调用了findPreambles函数，findPreambles函数的功能是：在每个信道的比特流中查找遥测字与交接字。
通过检查时间间隔(6秒)和对子帧中的前两个字进行奇偶校验来验证是否是帧头的位置。同时函数返回处于跟踪状态的通道列表，以及具有有效导航数据的导航帧头的位置。
具体的解读参考下文函数。
[subFrameStart, activeChnList] = findPreambles(trackResults, settings);③ 解码星历表
每个正确的同步头都标记了导航数据子帧的起始位置，每个子帧都包含300bit的字。理论上，每1ms可以得到1bit，但是，所处理的是带有噪声的微弱信号，因此每20ms计算一个平均值，并二值化为1和一1。每一导航比特的持续时间即为20ms。导航数据的比特率为50b/s。跟踪模块输出值的采样频率为1000sps，相当于1ms输出一一个值。在导航数据解码之前，跟踪模块输出的信号必须从1000sps 转换到50b/s。也就是说，每20个相邻值必须被替换为1个值。
该模块中涉及到ephemeris.m函数，其作用是解码星历表，详细解释见下文。
%% Decode ephemerides 解码星历表======================================
for channelNr = activeChnList
%=== 将tracking函数的输出转换为导航数据 =======================
%--- 从tracking函数的输出中，复制5个子帧长的记录 ---------------
navBitsSamples = trackResults(channelNr).I_P(subFrameStart(channelNr) - 20 : ...
subFrameStart(channelNr) + (1500 * 20) -1)';
%--- 每20ms作为一个值 ------------------------
navBitsSamples = reshape(navBitsSamples, ...
20, (size(navBitsSamples, 1) / 20));
%--- 将每20ms的平均值作为一个替换值 -------------
navBits = sum(navBitsSamples);
%--- Now threshold and make 1 and 0 -----------------------------------
% 表达式(navBits > 0)返回一个数组，如果满足条件则元素设置为1，如果不满足则设置为0。
navBits = (navBits > 0);
%--- Convert from decimal to binary -----------------------------------
% 函数的星历表需要二进制形式的输入。在Matlab中，它是一个只包含“0”和“1”字符的字符串数组。
navBitsBin = dec2bin(navBits);
%=== 解码星历表和第一子帧的转换字 ================
[eph(trackResults(channelNr).PRN), TOW] = ...
ephemeris(navBitsBin(2:1501)', navBitsBin(1));
%--- 排除数据不充足的导航数据 -----
if (isempty(eph(trackResults(channelNr).PRN).IODC) || ...
isempty(eph(trackResults(channelNr).PRN).IODE_sf2) || ...
isempty(eph(trackResults(channelNr).PRN).IODE_sf3))
%--- Exclude channel from the list (from further processing) ------
activeChnList = setdiff(activeChnList, channelNr);
%setdiff 两个数组的差集
end
end④ 检测卫星的数量是否一直在3颗以上并初始化
首先检测卫星的数量
%% Check if the number of satellites is still above 3 =====================
% 检测卫星的数量是否一直在3颗以上
if (isempty(activeChnList) || (size(activeChnList, 2) < 4))
% Show error message and exit
disp('Too few satellites with ephemeris data for postion calculations. Exiting!');
navSolutions = [];
eph = [];
return
end初始化卫星的仰角阵列，设置为无穷大，并保存有效频道的列表。
%% Initialization 初始化 =============================================
% 将卫星仰角阵列设置为无穷大，以包括第一次计算接收机位置时的所有卫星。
% 没有参考点来找到仰角，因为没有在这一点上的接收机位置估计。
satElev = inf(1, settings.numberOfChannels);
% 保存有效频道的列表。该列表包含被跟踪的卫星，并具有所需的星历数据。
% 在接下来的步骤中，列表包含的卫星将取决于每颗卫星的仰角，它将随着时间的变化而变化。
readyChnList = activeChnList;
transmitTime = TOW;
%% Initialization of current measurement ==================================
% 初始化当前
for currMeasNr = 1:fix((settings.msToProcess - max(subFrameStart)) / ...
settings.navSolPeriod)
% fix四舍五入 settings.navSolPeriod计算伪距和位置的周期
% Exclude satellites, that are belove elevation mask
activeChnList = intersect(find(satElev >= settings.elevationMask), ...
readyChnList);
% Save list of satellites used for position calculation
navSolutions.channel.PRN(activeChnList, currMeasNr) = ...
[trackResults(activeChnList).PRN];
% 这两条线有助于skyPlot函数。被排除的卫星将不会“跳”到天空图中的位置(0,0)。
navSolutions.channel.el(:, currMeasNr) = ...
NaN(settings.numberOfChannels, 1);
navSolutions.channel.az(:, currMeasNr) = ...
NaN(settings.numberOfChannels, 1);
⑤ 计算伪距
计算伪距需要调用calculatePseudoranges()函数，运行过程的错误就发生在calculatePseudoranges()函数中，具体错误分析详见calculatePseudoranges()函数。
%% Find pseudoranges ======================================================
% 计算伪距
navSolutions.channel.rawP(:, currMeasNr) = calculatePseudoranges(...
trackResults, ...
subFrameStart + settings.navSolPeriod * (currMeasNr-1), ...
activeChnList, settings);
% subFrameStart：某一通道有效导航数据的起始位置⑥ 找到卫星位置并完成时钟校正
找到卫星位置并完成时钟校正，调用satpos函数
%% Find satellites positions and clocks corrections =======================
% 找到卫星位置并完成时钟校正，调用satpos函数
[satPositions, satClkCorr] = satpos(transmitTime, ...
[trackResults(activeChnList).PRN], ...
eph, settings);
⑦ 计算接收机位置
计算过程主要调用了leastSquarePos（）函数，其作用是：计算最小二乘解。
利用伪距计算位置的最常用算法是最小二乘法，当观测量多于未知量时常使用该方法。程序描述了如何利用最小二乘法，从四颗或者更多卫星的伪距求解接收机的位置。具体运算步骤见《软件定义的..》P110
%% Find receiver position =================================================
% 计算接收机位置
% 只有当3颗以上卫星的信号可用时，才能找到接收器的三维位置
if size(activeChnList, 2) > 3
%=== Calculate receiver position ==================================
[xyzdt, ...
navSolutions.channel.el(activeChnList, currMeasNr), ...
navSolutions.channel.az(activeChnList, currMeasNr), ...
navSolutions.DOP(:, currMeasNr)] = ...
leastSquarePos(satPositions, ...
navSolutions.channel.rawP(activeChnList, currMeasNr)' + satClkCorr * settings.c, ...
settings);
%--- Save results -------------------------------------------------
navSolutions.X(currMeasNr) = xyzdt(1);
navSolutions.Y(currMeasNr) = xyzdt(2);
navSolutions.Z(currMeasNr) = xyzdt(3);
navSolutions.dt(currMeasNr) = xyzdt(4);
% Update the satellites elevations vector
satElev = navSolutions.channel.el(:, currMeasNr);
%=== Correct pseudorange measurements for clocks errors ===========
navSolutions.channel.correctedP(activeChnList, currMeasNr) = ...
navSolutions.channel.rawP(activeChnList, currMeasNr) + ...
satClkCorr' * settings.c + navSolutions.dt(currMeasNr);⑧ 坐标变换
通过调用cart2geo和cart2utm函数，可以将ECEF地心地固坐标系转换为大地坐标系和UTM坐标系。此外将结果定位结果放至navSolutions中。程序如下：
%% Coordinate conversion ==================================================
%=== Convert to geodetic coordinates ==============================
[navSolutions.latitude(currMeasNr), ...
navSolutions.longitude(currMeasNr), ...
navSolutions.height(currMeasNr)] = cart2geo(...
navSolutions.X(currMeasNr), ...
navSolutions.Y(currMeasNr), ...
navSolutions.Z(currMeasNr), ...
5);
%=== Convert to UTM coordinate system =============================
navSolutions.utmZone = findUtmZone(navSolutions.latitude(currMeasNr), ...
navSolutions.longitude(currMeasNr));
[navSolutions.E(currMeasNr), ...
navSolutions.N(currMeasNr), ...
navSolutions.U(currMeasNr)] = cart2utm(xyzdt(1), xyzdt(2), ...
xyzdt(3), ...
navSolutions.utmZone);
else % if size(activeChnList, 2) > 3
%--- There are not enough satellites to find 3D position ----------
disp([' Measurement No. ', num2str(currMeasNr), ...
': Not enough information for position solution.']);
%--- Set the missing solutions to NaN. These results will be
%excluded automatically in all plots. For DOP it is easier to use
%zeros. NaN values might need to be excluded from results in some
%of further processing to obtain correct results.
navSolutions.X(currMeasNr) = NaN;
navSolutions.Y(currMeasNr) = NaN;
navSolutions.Z(currMeasNr) = NaN;
navSolutions.dt(currMeasNr) = NaN;
navSolutions.DOP(:, currMeasNr) = zeros(5, 1);
navSolutions.latitude(currMeasNr) = NaN;
navSolutions.longitude(currMeasNr) = NaN;
navSolutions.height(currMeasNr) = NaN;
navSolutions.E(currMeasNr) = NaN;
navSolutions.N(currMeasNr) = NaN;
navSolutions.U(currMeasNr) = NaN;
navSolutions.channel.az(activeChnList, currMeasNr) = ...
NaN(1, length(activeChnList));
navSolutions.channel.el(activeChnList, currMeasNr) = ...
NaN(1, length(activeChnList));
% TODO: Know issue. Satellite positions are not updated if the
% satellites are excluded do to elevation mask. Therefore rasing
% satellites will be not included even if they will be above
% elevation mask at some point. This would be a good place to
% update positions of the excluded satellites.
end % if size(activeChnList, 2) > 3
%=== Update the transmit time ("measurement time") ====================
transmitTime = transmitTime + settings.navSolPeriod / 1000;
end %for currMeasNr...至此，介绍了GPS的定位算法，同时，已基本学习完所有的软件接收机函数，最后简要介绍一下plotNavigation.m函数。
plotNavigation.m函数
下图为经过定位后绘制出的图形，其主要描述了静态定位的精度，左下角的图形表示各个卫星的仰角。
Inputs:
navSolutions - 数据来源于PostNavigation函数，它包含测量的伪距和接收机坐标。
settings - 接收机设置。真正的接收机坐标包含在这个结构中。
findPreambles.m寻找帧头（同步头）函数
GPS导航数据解码中的首要问题是确定子帧的起始位置，该位置由8bit长的同步头标记，格式为10001011。由于科斯塔斯环具有跟踪信号180°相位翻转的能力，该同步头可能以反码01110100的形式出现。事实上，这两组数据可能出现在接收数据的任意位置，所以必须进行额外检查对同步头进行确认。确认过程检查重复出现同步头的时间是否为6s,对应于两个连续子帧播发的时间间隔。
同步头的搜索通过相关完成，相关函数的第一个输 人是接收到的导航数据比特序列,用1和-1表示。第二个输人是8bit的同步头，同样用1和-1表示。用-1和1代替0和1时,同步头出现时相关性函数的输出为8,同步头的反码出现时相关函数输出为-8。
区分哪一个最大相关值是子帧起始点的方法，包括对两个连续的最大相关值间隔的确定，只有最大相关值间的延迟为6s,且奇偶校验正确，才可以确定为子帧的起始位置。
找到准确的同步头后，每个子帧的数据便可以提取出来。如果相关性显示同步头为反码，则全部导航序列必须反转。
此函数对应postNavigation.m函数中模块②中的函数调用，其功能：找到每个信道的遥测字与交接字。
function [firstSubFrame, activeChnList] = findPreambles(trackResults, settings)Inputs:
trackResults - 追踪函数的输出
settings - Receiver settings.Outputs:
firstSubframe - 该数组包含每个通道中第一个帧头的位置，该位置是自跟踪开始以来的毫秒数。如果在通道中没有检测到有效的帧头，则将相应的值设置为0。
activeChnList - 包含有效帧头信息的通道列表
首先该函数初始化了数组并给出遥测字的矩阵[1 -1 -1 -1 1 -1 1 1]：
% 该函数可以延迟到跟踪结果中的一个较晚的点，以避免跟踪环路瞬变引起的噪声
searchStartOffset = 0;
%--- Initialize the firstSubFrame array -----------------------------------
% 初始化firstSubFrame数组
firstSubFrame = zeros(1, settings.numberOfChannels);
%--- Generate the preamble pattern ----------------------------------------
% 生成preamble模式，每一帧的第一个字为遥测字，遥测字的首八个比特固定位10001011的同步码
preamble_bits = [1 -1 -1 -1 1 -1 1 1];
% "Upsample" the preamble - make 20 vales per one bit. The preamble must be
% found with precision of a sample.
preamble_ms = kron(preamble_bits, ones(1, 20));
% kron函数：如果 A 是 m×n 矩阵，而 B 是 p×q 矩阵，则 kron(A,B) 是通过获取 A 元素与
% 矩阵 B 元素之间的所有可能积而形成的一个 m*p×n*q 矩阵。
%--- Make a list of channels excluding not tracking channels --------------
% 列出只包含已追踪到的卫星通道
activeChnList = find([trackResults.status] ~= '-');然后开始寻找与遥测字矩阵完全一致或完全相反的导航数据：
%=== For all tracking channels ...
for channelNr = activeChnList
%% 追踪函数的输出，与遥测字求相关 ================================
% trackResults包含导航位。相关记录直接从searchStartOffset这里开始，以避免跟踪环路瞬变。
bits = trackResults(channelNr).I_P(1 + searchStartOffset : end);
% 理论上1ms得到1bit。但是，所处理的是带有噪声的微弱信号，
% 因此每20ms计算一个平均值，并二值化为1或-1
bits(bits > 0) = 1;
bits(bits <= 0) = -1;
% 求相关
tlmXcorrResult = xcorr(bits, preamble_ms);
之后分析并确定同步头的位置
%% Find all starting points off all preamble like patterns ================
% 找到所有与遥测字相类似的起始点
clear index
clear index2
xcorrLength = (length(tlmXcorrResult) + 1) /2;
%--- Find at what index/ms the preambles start ------------------------
% 找到遥测字开始的地方/毫秒
index = find(...
abs(tlmXcorrResult(xcorrLength : xcorrLength * 2 - 1)) > 153)' + ...
searchStartOffset;对每个可能的起点进行分析，并使用奇偶校验的方法（navPartyChk函数判定）确定。
%% 分析检测到的遥测字疑似起始点 ================================
for i = 1:size(index) % For each occurrence 对每个疑似点进行排查
% 找出这个起点和其他疑似起点间的时间间隔。如果间隔是6000毫秒(一个子帧)，
% 则通过验证两个GPS子帧的奇偶性来进一步验证
index2 = index - index(i);
if (~isempty(find(index2 == 6000)))
% 通过检查子帧中前两个单词的奇偶性来验证起始点的真实性。假定数据跳变的边界是已知的。因此，将超过20ms的位值组合起来，
% 以提高接收机在噪声信号下的性能。共读取62位:需要前一帧的2位进行奇偶校验;
% 前两个30位字(遥测字TLM和交接字HOW)用60位。索引指向遥测字字的开始。
bits = trackResults(channelNr).I_P(index(i)-40 : ...
index(i) + 20 * 60 -1)';
%--- 将20ms的值平均后变为1个值 ------------------------
bits = reshape(bits, 20, (size(bits, 1) / 20));
bits = sum(bits);
% Now threshold and make it -1 and +1
bits(bits > 0) = 1;
bits(bits <= 0) = -1;
%--- 进行奇偶校验 ----------------
if (navPartyChk(bits(1:32)) ~= 0) && ...
(navPartyChk(bits(31:62)) ~= 0)
% 判断奇偶校验是否正常，记录遥测字的起始位置，并跳过此通道的遥测字
% 检查的其余部分，处理下一个通道。
firstSubFrame(channelNr) = index(i);
break;
end % if parity is OK ...
end % if (~isempty(find(index2 == 6000)))
end % for i = 1:size(index)
% 如果没有检测到有效的起始点，则从有效通道列表中排除该内容
if firstSubFrame(channelNr) == 0
% 进一步处理排除通道。由于不包含任何有效的起始位，因此无法执行
activeChnList = setdiff(activeChnList, channelNr);
disp(['Could not find valid preambles in channel ', ...
num2str(channelNr),'!']);
end
end % for channelNr = activeChnList在进行寻找帧头（同步头）时，又牵涉到了一个函数——navPartyChk函数（奇偶校验函数），下面就对该函数简要分析。
navPartyChk.m奇偶校验函数
每一字中的奇偶检验码(汉明编码)可以帮助用户接收机检查经解调得到的字中是否包含错误比特，并且它还有一定的比特纠错功能。一个字中的奇偶检验码是通过对该字的前24上一字的最后两个比特按照以下公式计算产生的:
其中，d1,d2,...,d24是24个原始数据比特，D1,D2,...,D24是卫星世纪发社的24个数据比特，D25,D26,...,D30是卫星实际发射的6个奇偶校验码，D29-和D30-是卫星实际发射的上一字中的最后两位奇偶校验码。
同时也可以表达成如下的矩阵运算：
矩阵H中的元素值1与0分别代表相应的比特值D29-，D30-和di(i=1,2, .，,24)是否参加异或相加运算，其中1代表参加，0代表不参加。观察矩阵H的元素分布，我们从中可以发现如下一个编码规律:矩阵H中的第2至第5行基本上分别是由它们相应的前一-行向右旋转平移一个元素得到的。另外，上述的矩阵相乘运算可非常有效地利用整数位操作加以完成。
ephemeris.m函数_获取星历表
函数从给定的位流中解码遥测字和TOW。位流(数组)必须包含1500位。数组中的第一个元素必须是子帧的第一个比特，其中调用了checkPhase.m函数，随后进行介绍。
function [eph, TOW] = ephemeris(bits, D30Star)每一颗卫星都是一帧接着一帧地发送导航电文，而在发送每帧电文时，卫星又以一子帧一子帧的形式发送。每帧电文长1500比特包含5个子帧。（参考文献：《GPS原理与接收机设计》P35-37）
第一数据块
第一数据块包括第1子帧中的第3字至第10字，它又常称为时钟数据块。由一颗卫星播发的时钟数据块提供该卫星的时钟校正参数和健康状态等如下内容。
(1)星期数(WN)：星期数来自Z计数器的高10位，指代当前的GPS星期。每当周内时计数在星期六午夜零时从最大值跳回零的同时，星期数的值加1。
(2)用户测距精度(URA)：用户测距精度是对所有由GPS地面监控部分和空间星座部分引起的测距误差大小的一个统计值，它是通过导航电文中的一个由4比特表示的用户测距精度因子N而提供给非特许用户的。用户测距精度因子N的值在0至15之间，每一个值对应于一个用户测距精度URA,URA值越大，则表示从该卫星信号中得到的GPS距离测量值的精度越低。当N等于15时，URA的估算值缺省，此时用户要自己承担使用该卫星的风险。
(3)卫星健康状况：若共计6比特的卫星健康状况的最高位是0，则表示导航电文全部正确。著它的最高位是1,则表示导航电文出错，而低5位又具体指出信号各部分的出错情况。
(4)时钟校正参数：a_f0, a_f1 和a_f2是卫星时钟校正模型方程中的三个系数。另外，参数1称为第一数据块的参考时间，它在时钟校正模型中被用做时间参考点。
(5)群波延时校正值：群波延时校正值TeD只适合于单频(L1或L2)接收机，而双频接收机则无需此项校正。单频接收机之所以有此项校正，是因为时钟校正参数a_f0是针对双频测量值而言的。
(6)时钟数据期号( IODC):用10比特表示的IODC是时钟数据块的“期刊号”，-个IODC值对应一套时钟校正参数。
第二数据块
卫星信号的第2子帧和第3子帧数据块起组成第二数据块，它提供该卫星自身的星历参数。
星历的原意是一张用来精确描述卫星在各个时刻的空间位置和运行速度的大表格。为了减少需要播发的数据量，GPS用开普勒方程来描述卫星的运行轨道，并通过最小二乘法逼近来求解方程中的各个系数。表2.4列出了由GPS卫星播发的共计16个开普勒系数的一套星历参数。
一套星历参数的有效期一般是以参考时间toe为中心的4小时之内，而超过此有效时段的星历经常被认为是过期且无效的。因为由过期星历参数计算得到的卫星轨道值一般会存在一个较大的误差，所以它们通常不能用于GPS的正常定位计算中。为了判断一套星历参数是否有效，正如上一小节所指出的那样,我们通常需要同时查看它的参考时间和星期数。
第三数据块
第三数据块由第4子帧和第5子帧的数据块组成。每颗卫星播发的第三数据块主要提供所有(自身和其他)卫星的历书参数、电离层延时校正参数、GPS 时间与UTC之间的关系以及卫尾健康状况等数据信息。与前两个数据块不同，第三数据块的内容并不是接收机在实现定位前所急需获得的。
下面需要对接收到的比特流信号进行解码，主要分为两大部分，分别为检验是否符合解码标准、进行分割解码两部分。
第一部分：需要检查数据数否充足、确定星历表是否为strings类型
第二部分：需要将数据流进行分组、纠正数据流的极性、识别子帧序号，解码相应数据。
由于程序过长且易懂，故将程序保存至文末网盘中。
checkPhase.m函数
该函数的主要功能是检验数据流中的相位是否反相，如果反相，那么将相位翻转。检查所提供的30位字的奇偶校验。前一个字的最后一个奇偶校验位用于计算。GPS标准定位服务信号规范提供了有关程序的说明。程序如下：
function word = checkPhase(word, D30Star)
if D30Star == '1'
% Data bits must be inverted
word(1:24) = invert(word(1:24));
endcalculatePseudoranges()函数_伪距计算函数
伪距计算的流程图如下所示，该函数计算所有追踪卫星的伪距。为跟踪卫星的通道均不会计算伪距，未找到信号帧头或者有奇偶校验错误无的通道也不会计算伪距，此函数的输入是桐乡即时支路相关器的输出，同时该函数也需要精确的码相位信息，以上两者均来自结构变量trackResult。
function [pseudoranges] = calculatePseudoranges(trackResults, ...
msOfTheSignal, ...
channelList, settings)该函数的输入输出要素如下：
Inputs:
trackResults - tracking函数的输出结果
msOfTheSignal - trackResults结构体中的伪距测量起始点(毫秒)
channelList - 经过处理后的通道列表
settings - receiver settingsOutputs:
pseudoranges - relative pseudoranges to the satellites. 与卫星的相对伪距。
在运行该函数时，程序发生报错：
错误提示：
需要大括号或点索引表达式中的一个输出，但结果有 8 个。
错误发生在该函数的65行，即：
for channelNr = channelList
%--- Compute the travel times -----------------------------------------
travelTime(channelNr) = ...
trackResults(channelNr).absoluteSample(msOfTheSignal(channelNr)) / samplesPerCode;
end该程序段中，for循环条件给出的channelNr = channelList格式不正确，应该将其改为channelNr = 1:length(channelList)，这样就可以顺利地将值付给channelNr。
同时，附上详细报错解释：软件定义的GPS和伽利略接收机这本书提供的matlab代码运行出错 - 手持设备 - EETOP 创芯网论坛 (原名：电子顶级开发网)
该程序的主要思想：
在识别出同步头后，就确定了所有可用卫星子帧的起点。信号从卫星到地球的传输时间为65ms~83ms，此值可用于确定初始伪距。距离地球最近的卫星到达的子帧也最早。因此通道N的卫星信号传输时间为65ms,剩余通道的传输时间可以根据通道N计算出来。
程序代码注释如下：
%--- Set initial travel time to infinity 将初始化时间设置为无穷大-----------
% 因为需要在代码中将选择最短的伪距。因此，来自非跟踪信道的伪距必须是最长的，所以选择无穷大。
travelTime = inf(1, settings.numberOfChannels);
% Find number of samples per spreading code 找出每个码片的样本数量
samplesPerCode = round(settings.samplingFreq / ...
(settings.codeFreqBasis / settings.codeLength));
%--- For all channels in the list ...
for channelNr = channelList
%--- Compute the travel times -----------------------------------------
% 计算采样传输时间
travelTime(channelNr) = ...
trackResults(channelNr).absoluteSample(msOfTheSignal(channelNr)) / samplesPerCode;
end
%--- Truncate the travelTime and compute pseudoranges ---------------------
minimum = floor(min(travelTime));
travelTime = travelTime - minimum + settings.startOffset;
%--- Convert travel time to a distance ------------------------------------
% The speed of light must be converted from meters per second to meters
% per millisecond. 计算出每个卫星距离接收机的伪距
pseudoranges = travelTime * (settings.c / 1000);satpos.m函数_计算卫星位置
其函数功能是：给定星历表相位，计算X,Y,Z卫星在传送时的坐标，为PRNLIST列表中的每个卫星计算坐标。
根据前文的epChemeris（）函数，可以得到eph的星历表，进而求解出卫星的时钟校正和位置。
Inputs:
transmitTime - transmission time
prnList - list of PRN-s to be processed
eph - ephemeridies of satellites
settings - receiver settingsOutputs:
satPositions - 卫星位置(in ECEF system [X; Y; Z;])
satClkCorr - 卫星的时钟校正
leastSquarePos.m函数_使用最小二乘法计算接收机位置
leastSquarePos.m函数通过计算伪距计算接收机的位置。如果一颗卫星的完整星历无法得到，那么，该卫星的伪距也无法使用。在postNavigation函数的最后，ECEF坐标被转化为UTM何打的坐标系统，结果保存在navSolutions结构变量中。最小二乘位置结算流程如下：
function [pos, el, az, dop] = leastSquarePos(satpos, obs, settings)Inputs:
satpos - 卫星的位置
obs - Observations - 每颗卫星的伪距测量值: (e.g. [20000000 21000000 .... .... .... .... ....])
settings - receiver settingsOutputs:
pos - receiver position and receiver clock error
(in ECEF system: [X, Y, Z, dt])
el - Satellites elevation angles (degrees)
az - Satellites azimuth angles (degrees)
dop - Dilutions Of Precision ([GDOP PDOP HDOP VDOP TDOP])
扫一扫
425
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
