第5章GPS基本原理及其Matlab仿真课件

第5章GPS卫星信号的捕获5.1概述5.2GPS卫星信号的多普勒效应5.3GPS卫星信号捕获的考虑5.4GPS卫星信号的捕获方法5.5Matlab与Simulink仿真简介5.6GPS卫星信号捕获的例子5.7关于捕获的一些子程序第5章GPS卫星信号的捕获5.1概述15.1概

述

为了跟踪和解码GPS信号,首先要捕获到GPS信号。将捕获到的GPS信号的必要参数立刻传递给跟踪过程，再通过跟踪过程便可得到卫星的导航电文。GPS卫星处于高速运动中，因此，其频率会产生多普勒频移。载波频率与C/A码的多普勒频移将在下面详细讲述。为覆盖高速卫星预期中的所有多普勒频率范围，捕获方法覆盖的频率范围必须在±10kHz之内。一旦捕获到GPS信号，立刻去测量两个重要参数：C/A码的起始点和载波频率（因为多普勒频移而变化）。接收机接收到的一系列数据往往包含多个卫星信号，每个信号具有不同的C/A码的不同起始点和不同的多普勒频率。

5.1概述为了跟踪和解码GPS信号,首先要捕获2针对某个特定的卫星信号，捕获过程就是要找到C/A码的起始点，并利用找到的起始点展开C/A码频谱，一旦复现了C/A码的频谱，输出信号将变成连续波（ContinuousWave，CW），于是便得到其载波频率。也就是说，捕获过程就是要获得输入信号的C/A码的起始点和载波频率，然后传递给跟踪过程。捕获与跟踪过程所用到的数据都是从原始的卫星信号经过下变频器（即与中频混频）之后收集到的，其中频（IF）为21.25MHz,采样频率为5MHz，信号的中心频率为1.25MHz。这组数据来源于我们设计的卫星信号接收装置，但是用模拟仿真算法也可以产生类似的数据，以验证本章介绍的卫星信号的捕获算法。

针对某个特定的卫星信号，捕获过程就是要找到C/A码35.2GPS卫星信号的多普勒效应GPS卫星绕地球旋转一周的时间是11h58min2.05s，从卫星轨道的近似半径可以求得卫星的角速度dθ/dt和运动速度vs：

(5-1)5.2GPS卫星信号的多普勒效应GPS卫星绕地球旋转一4这里rs表卫星轨道的平均半径。我们知道，一个太阳日和一个恒星日之间相差3min55.91s，在这段时间里，卫星大约运行了914km(3874m/s×235.91s)，对应地球表面与卫星的最高点，相应的角度近似为0.045rad(914/20192)或2.6°。如果卫星接近地平线，相应的角度为0.035rad或者2°。因此我们可以看出，对于地球表面的固定一点，在每天的同一时间里，卫星位置大约改变2°～2.6°。这里rs表卫星轨道的平均半径。我们知道，一个太阳日和一个恒星5如图5-1所示，卫星在位置S处，用户在位置A处，相对用户的卫星角速度vd造成了多普勒频移，其中vd值为：(5-2)根据卫星轨道速度，取水平方向的最大值，

得到多普勒角速度的最大值为：

(5-3)

如图5-1所示，卫星在位置S处，用户在位置A处6图5-1卫星运动引起的多普勒频移

图5-1卫星运动引起的多普勒频移7这个值与高速航空卫星的速度相等。由地面设施引起的多普勒频移通常非常小，即使它直接相对卫星运动引起了最大多普勒效应。对由C/A码调制过的频率L1(f=1575.42MHz），最大的频移为：

(5-4)

这里c是光速。

因而，

对于一固定观测器来说，最大的多普勒频移是±5kHz。

这个值与高速航空卫星的速度相等。由地面设施引起的多普勒频移通8如果地面GPS接收设备高速移动，就要考虑进多普勒效应了。要单独产生一个±5kHz的频率范围，接收设备需以2078m/h的速度相对卫星运动，这个速度将包含最高速度的航空器。因而，在设计GPS接收机时，如果接收机用在低速设备，则认为载波频率的多普勒频移范围在±5kHz；相反地，如果用在高速设备，就要合理假定其多普勒频移的范围在±10kHz。这些值对于确定捕获过程的搜索频率范围是至关重要的。因此，就有了引言中我们提到的为覆盖高速卫星预期中的所有多普勒频率范围，捕获方法覆盖的频率范围必须在±10kHz之内之说。

如果地面GPS接收设备高速移动，就要考虑进多普勒效应9由于C/A码频率非常低，所以C/A码上的多普勒频移也是非常小的。C/A码的频率是1.023MHz，比载波频率低1540(1575.42/1.023)倍。C/A码的多普勒频率为：

(5-5)

如果接收机高速移动，这个值就得取双倍，为6.4Hz，得到的这个值对跟踪过程（BASS方法）来说是至关重要的，这点我们将在下一章讨论。在BASS跟踪过程中，输入信号与本地产生数据必须要紧密对齐，而C/A码的多普勒频移可能造成输入码与本地生成码的不对齐。

由于C/A码频率非常低，所以C/A码上的多普勒频移也10在数字化卫星信号中，如果数据用5MHz采样(称为采样频率)，则每个采样之间相隔200ns（称为采样时间）。在跟踪过程中，我们期望本地生成信号与输入信号未对齐的长度在半个采样时间或近似100ns之内，若两个信号之间相差超过这个长度，跟踪将失锁，即失去跟踪灵敏性。C/A码的基波时间为977.5ns或1/1.023×106s，它移动一个周期需花费156.3ms（1/6.4）。对这句话的理解是：多普勒频移使C/A码的频率由1.023×106变为[JP2]1.023×106+6.4，即每秒内多变化了6.4个周期，则多变化一个周期的时间为(1/6.4)ms。在数字化卫星信号中，如果数据用5MHz采样11因此，它移动100ns的数据长度，近似花费16ms(100×156.3/977.5)。在高速导航器中，需每16ms选择一批数据，以保证输入信号与本地生成码更好地匹配。但是由于输入信号中有噪声，用1ms的数据可能不能非常精确地匹配，因而可以拓宽输入信号的调整时间，每20ms选择一批输入信号。

对一个低速导航器来说，这个时间可拓宽到40ms。

因此，它移动100ns的数据长度，近似花费1612从以上讨论我们知道，输入信号的调整取决于采样频率。采样频率高将缩短调整时间，因为采样时间变短，在半个周期内匹配输入码与本地生成码的时间变短。如果输入信号的强弱和跟踪灵敏度不成问题，则输入信号的调整时间可以拓宽。然而，输入信号与本地生成信号相离必须在半个基波或488.75ns(977.5/2)之内，这个时间看做是输入信号与本地生成信号相离允许的最大值。C/A码的多普勒频率是6.4Hz，调整时间可以拓宽到78.15ms（1/2×6.4）。

从以上讨论我们知道，输入信号的调整取决于采样频率。135.3GPS卫星信号捕获的考虑

5.3.1捕获时的最大电文长度

C/A码长1ms，那么至少要用1ms的电文来捕获，甚至只用1ms的电文来捕获时，都可能发生导航电文相位偏移。如果电文中有数据偏移，电文的下一个1ms将不会含有偏移。因此，为了保证捕获电文中不含有数据偏移，需要用两组连续的电文来捕获，这个电文最大长度是10ms。如果使用两组连续的10ms电文来捕获，

就保证了在某一组电文中不含相位偏移。

5.3GPS卫星信号捕获的考虑5.3.1捕获时的最大电14限制电文长度的第二个因素是C/A码的多普勒效应。如果理想相关峰值是1，当C/A码超前或滞后1/2码时，相关的峰值将降为0.5，对应地，幅值下降6dB。如上节所讨论的，假定未对准的C/A码允许在半个基波（0.489μs）内,基波频率是1.023MHz，C/A码上预期的最大多普勒频移是6.4Hz，那么两个相差6.4Hz的频率交换半个基波需要78ms（1/2×6.4），这个数据远大于10ms。限制电文长度的第二个因素是C/A码的多普勒效155.3.2捕获中的频率步长捕获时的另一个考虑因素是捕获中剥离载波所需的频率。在前文曾提过，需被覆盖的多普勒频率范围是±10kHz，决定覆盖这20kHz范围的步进频率非常重要。步进频率的大小与捕获中的电文长度紧密相关。如果输入信号与本地混频信号相距1个周期，它们两者没有相关性；如果两者小于1个周期，则它们有部分相关性。随意选取两个信号允许的最大频率相离是0.5个周期，如果电文记录长度为1ms，一个1kHz的信号将在1ms内变换一个周期，为了保持1ms内的最大频率相离在0.5个周期，步进频率须为1kHz。在这种情况下，在输入信号与卷积信号之间最大限度的频率相离是500Hz/s或0.5Hz/ms，且输入信号刚好在两个频率点之间。如果电文记录长度为10ms，100Hz的搜索步进频率就可以满足要求。看待这个问题的一个简单方法是频率相离与电文长度成反向关系。

5.3.2捕获中的频率步长16由以上的讨论可以得出如下结论：当捕获输入电文为1ms长时，步进频率是1kHz；当电文为10ms时，步进频率是100Hz。从这个简单的讨论可以明显看出，捕获中的操作执行次数与总的数据点并不是成线性比例关系，当电文长度从1ms上升到10ms时，电文长度上升10倍，频率点数目也上升了10倍，捕获所需的操作次数上升了不止10倍。因此，实际捕获时，如果强调捕获的速度的话，电文长度需保持在最小值。执行次数的增加，取决于实际捕获所用的方法，下面将讨论具体的捕获方法。

由以上的讨论可以得出如下结论：当捕获输入电文为1m175.4GPS卫星信号的捕获方法5.4.1传统捕获方法捕获的基本思想就是展开输入信号，找到载波频率。如果相位正确的C/A码与输入信号相乘，输入信号将如图5-2中所示变成一个连续信号。图5-2中最上端是输入信号，这个信号是C/A码编码的射频信号相位，注意到射频信号与C/A码是任意选择的，它们并不代表卫星传输过来的真正的信号。中间的图表示C/A码，其值在±1之间。最下端的图是连续信号，代表输入信号与C/A码的乘积，其对应的频谱不再是扩频，而是一个连续信号，这个过程有时称做输入信号的C/A码剥离。

5.4GPS卫星信号的捕获方法5.4.1传统捕获方法18图5-2C/A编码的输入信号与本地C/A码相乘

图5-2C/A编码的输入信号与本地C/A码相乘19一旦输入变成连续信号，通过傅立叶变换可以知晓其频率。如果输入电文长是1ms，则傅立叶变换的频率分辨率为1kHz。可以设置一个阈值，即一个门限来判断一个频率分量是否足够大。高于这个门限的最高频率分量就是我们要找的频率。如果用5MHz使输入信号数字化，1ms的电文将含有5000个数据点，一个5000点的FFT(快速傅立叶变换)将产生5000个频率分量。然而，5000个频率分量中只有前2500个频率分量包含有用信息，后2500个是前2500个频率分量的复共轭。频率分辨率是1kHz，因此，FFT覆盖的总频率范围是2.5MHz，是采样频率的一半。但我们感兴趣的频率范围只有20kHz，而非2.5MHz，因此，为了节省计算时间，

可以用DFT(离散傅立叶)方法，

只计算相隔1kHz的21个频率分量。

一旦输入变成连续信号，通过傅立叶变换可以知晓其频率。20由于输入信号中C/A码的起始点还未知，因此，必须先找到这个点。为了找到这个点，我们产生本地的一个C/A码，将其数字化为5000个点，且与输入信号点对点相乘，对此乘积的结果进行FFT或DFT，可以找到其频率。为了搜索1ms的电文，输入电文与本地产生数字化的C/A码必须彼此相对移动5000次。如果使用FFT，需操作5000次，每次操作包含5000个点对点的乘积和一个5000点的FFT。由于只有2500个频率分量提供信息，其余2500个分量提供的是多余的信息，因此，输出的5000点数据每个含有2500个频率分量，这样，在频域共有1.25×107(5000×2500)个输出。在这1.25×107个输出中最高幅值如果大于预先设定的门限值，则该分量就是期望值。在如此多的数据中寻找最高频率分量也是非常费时的。由于在期望的20kHz的FFT输出中，只有21个频率分量是我们感兴趣的，总输出可以降低到105000(5000×21)个。通过这种方法，可以用时间分辨率为200ns(1/5MHz)和1kHz的频率分辨率找出C/A码的起始点。

由于输入信号中C/A码的起始点还未知，因此，必须先找21如果使用10ms的电文，因为信号只需要卷积1ms，因此需要5000次操作。每次操作包含50000个点对点乘积和50000次FFT，总共有1.25×108(5000×25000)个输出结果。如果在预期的20kHz范围内只考虑201个频率分量，需从1005000(5000×201)个结果中挑选，因此可以看出，从1ms到10ms操作的增长是非常显著的。

寻找C/A码起始点的时间分辨率仍是200ns，

但频率分辨率提高到了100Hz。

如果使用10ms的电文，因为信号只需要卷积1ms22传统捕获方法是由硬件实现的，硬件主要用以执行以上所讲的处理。假定输入信号用5MHz数字化，一个可行方法是产生5000点的数字C/A码，并与输入信号点对点相乘，每200ns执行一次5000点相乘，每200ns对一个5000点相乘的结果进行FFT频率分析，图5-3示意了该过程。如果C/A码与输入电文匹配，FFT结果中将含有一个强分量，如前所述，这个方法将产生1.25×107(5000×2500)个输出结果，然而，只有在±10kHz频率范围内的输出结果中挑选，这个约束使得排序过程变简单了。

传统捕获方法是由硬件实现的，硬件主要用以执行以上所讲23图5-3用C/A码和频率分析进行捕获

图5-3用C/A码和频率分析进行捕获24完成这个操作的另一个方法是通过离散傅立叶变换(DFT),将本地产生码变换成含有C/A码和射频信号的电文。射频是复数，可用ejwt表示。本地码信号由射频与C/A码乘积获得，因此，本地码也是复数。假定L1频率(1575.42MHz)被转换到21.25MHz，并用5MHz数字化，输出频率将会是1.25MHz;同时假定捕获过程以1kHz的步长搜索1250±10kHz的频率范围，且总共有21个频率分量，本地码lsi可用下式表示：

(5-6)完成这个操作的另一个方法是通过离散傅立叶变换(DFT25其中，s表示卫星的编号；i=1,2,3,…,21；Cs表示编号为s的卫星的C/A码；fi=1250-10,1250-9,1250-8,…,1250+10kHz。本地信号须用5MHz数字化，产生5000个数据点，这21组数据代表相距1kHz的21个频率。这些数据与输入信号有关，如果本地产生信号包含正确的C/A码和正确的频率分量，当遇到正确的C/A码相位时，输出结果将会特别高。其中，s表示卫星的编号；i=1,2,3,…,21；C26图5-4通过本地产生C/A码和射频码捕获

图5-4通过本地产生C/A码和射频码捕获275.4.2循环相关捕获方法1.循环卷积与循环相关在讲述循环相关捕获方法之前，首先简单介绍执行相关操作的主要数学运算。如果一个输入信号经过一个线性时不变系统，可由卷积运算得到时域结果，也可由傅立叶变换得到频域结果。如果系统脉冲响应是h(t)，输入信号x(t)可通过卷积产生如下的输出y(t)：(5-7)5.4.2循环相关捕获方法(5-7)28通过傅立叶变换，y(t)的频域响应为：

(5-8)令变量t-τ=u，那么

(5-9)通过傅立叶变换，y(t)的频域响应为：(5-8)令变量t29为得到时域输出，需对Y(f)进行反傅立叶变换，

结果可写成：

(5-10)

此处，*表示卷积，F-1表示傅立叶反变换。可以发现一个简单的关系，时域的卷积等效于频域的傅立叶变换，或者说频域卷积运算等效于时域傅立叶变换，

这两种关系可以写成：

(5-11)

通常称其为傅立叶变换中卷积的二元性。

为得到时域输出，需对Y(f)进行反傅立叶变换，结果可写成：30这个思想还可用于离散系统。与连续时域表达式不同，响应y(n)可以表达成：

(5-12)此处，x(m)表示输入信号，h(n-m)是离散时域系统响应。注意到，该式子中，h(n-m)中的时移是循环的，因为离散操作是周期性的。对上式进行离散傅立叶变换(DFT)，结果为：(5-13)这个思想还可用于离散系统。与连续时域表达式不同，响应31然而，捕获时用到的并不是卷积运算，而是相关运算，它与卷积运算不同。x(n)与h(n)之间的相关表示为：

(5-14)上式与式(5-12)的惟一不同点就在于h(n+m)中的m前面的符号不同。h(n)是另一信号，而非线性系统的脉冲响应。z(n)的DFT变换为：(5-15)

然而，捕获时用到的并不是卷积运算，而是相关运算，它32这里X-1(k)表示离散傅立叶反变换。上式还可写成：

(5-16)

如果x(n)是实数，x(n)*=x(n)，其中*表示复共轭。有了这个关系，Z(k)的幅值为：(5-17)利用上式的这个关系，可得到输入信号与本地产生信号的相关性。如前所述，这个式子提供了循环相关性。

这里X-1(k)表示离散傅立叶反变换。上式还可写成：(5-33

2.循环相关捕获的方法上面讨论了循环相关的运算法则，利用这些运算法则就可以执行捕获了。然而实际捕获中，输入电文不是连续到达接收机的，因此循环相关操作适合一组或一批电文。输入电文经过5MHz的ADC(模拟/数字转换器)采样存储在存储器中，只有1ms的输入电文用来寻找C/A码的起始点，其搜索频率分辨率即步进频率是1MHz。

2.循环相关捕获的方法34捕获输入电文，需经过以下几个步骤：(1)对1ms的输入电文x(n)进行快速傅立叶变换FFT，将输入转换到频域，值为X(k)，这里，n=k=0,1,2,…,4999。(2)取X(k)的复共轭，值为X(k)*。(3)利用式(5-6)，产生21个本地码lsi(n)(i=1,2,…,21)，本地码包含了卫星s的C/A码和一个复射频信号的乘积，然后用5MHz采样本地码，本地码的频率fi相距1MHz。(4)对lsi(n)进行FFT，变换到频域中，值为Lsi(k)。(5)将X(k)*与Lsi(k)点对点相乘，结果为Rsi(k)。捕获输入电文，需经过以下几个步骤：35(6)对Rsi(k)进行反傅立叶变换，变换到时域中的值为rsi(n)，得到其绝对值|rsi(n)|，总共有105000(5000×21)个|rsi(n)|。(7)在输入电文200ns的时间分辨率和载波频率为1kHz分辨率的条件下，|rsi(n)|最大值中的第n位和第i个载波频率芯给出了C/A码的初始点。(6)对Rsi(k)进行反傅立叶变换，变换到时域中的值36以上的操作可由图5-5说明。图5-5所示结果是时域中的，且只给出21个本地码的其中一个。可认为输入电文与本地电文位于两个圆柱体表面，为了去匹配输入电文，本地码要旋转5000次。换句话说，一个圆柱体相对于另一个圆柱体旋转5000次。在每一步，5000个输入电文与5000个本地电文点对点相乘，相乘结果加到一起。包含本地码与输入码所有可能的乘积需5000步，乘积中最高幅值将被记录下。最高幅值如果同时大于预置门限值的话，就是我们的期望值。

以上的操作可由图5-5说明。图5-5所示结果是时域37图5-5循环相关捕获示意图

图5-5循环相关捕获示意图38

3.循环相关捕获的改进方法循环相关捕获的改进方法与上面介绍的捕获方法的惟一区别就是傅立叶变换FFT的长度减小了一半。上一节的循环相关方法的第(3)步产生本地码lsi(n)，由于lsi(n)是复量，故其频谱是不对称的，如图5-6所示。3.循环相关捕获的改进方法39图5-6本地产生信号频谱图

图5-6本地产生信号频谱图40从图5-6中显而易见，其信息包含在前半频谱图中，后半部分只含有非常少量的信息，从而，循环相关的捕获方法可以做如下修改：(1)对1ms的输入电文x(n)进行快速傅立叶变换FFT，将输入转换到频域，值为X(k)，这里，n=k=0,1,2,…,4999。(2)取X(k)的前2500(k=0,1,2,…,4999)个值的复共轭，值为X(k)*。(3)利用式(5-6)，产生21个本地码lsi(n)(i=1,2,…,21)，每个lsi(n)有5000个数据点。从图5-6中显而易见，其信息包含在前半频谱图中，后半部分41(4)对lsi(n)进行FFT，变换到频域中，值为Lsi(k)。(5)Lsi(k)的后半部只含有少量信息，因此，取前半部分的Lsi(k)，与X(k)*点对点相乘，结果为Rsi(k)(k=0,1,2,…,4999)。(6)对Rsi(k)进行反傅立叶变换，变换到时域中的值为rsi(n)，并求出其绝对值|rsi(n)|，总共有52500(2500×21)个|rsi(n)|。(4)对lsi(n)进行FFT，变换到频域42(7)|rsi(n)|的最大值如果大于预置门限值的话就是我们的期望值。第i个频率就是分辨率1kHz的载波频率，第n位就是分辨率400ns的C/A码的初始点。(8)由于这种方法寻找C/A码起始点的时间分辨率是400ns，可以通过比较第n位的幅值与(n-1)和(n+1)位幅值，将分辨率提高到200ns。该方法中，从第(5)步到第(7)步，只有2500个数据点的操作，而非5000个。第(7)步挑选最大幅值的过程由于只用到一半的输出而简化了。因此这种方法节省了操作时间。仿真结果表明，这种方法的信噪比略低于常规循环相关捕获方法1.1dB，这可能是由于丢弃频域内后半部分信号引起的。(7)|rsi(n)|的最大值如果大于预435.4.3延迟与累积捕获方法基于延迟与累积的捕获方法是为了去除输入信号中的频率信息。去除了频率信息，只需用本地C/A产生码去寻找输入信号C/A码的起始点，一旦找到其起始点，就可由FFT或DFT得到其频率。从理论角度看，这是个非常吸引人的方法，然而，运用到实际中处理GPS信号，还需进一步研究。这种方法讨论如下。

首先，我们假定输入信号s(t)是复量，因此表示为：

(5-18)

5.4.3延迟与累积捕获方法(5-18)44这里，Cs(t)表示卫星s的C/A码。

输入信号的延迟表达式为：

(5-19)

τ表示延迟时间，s(t)与延迟表达式的复共轭的乘积为：

(5-20)这里的

(5-21)这里，Cs(t)表示卫星s的C/A码。输入信号的延迟表达式45戈尔德码(GoldCode)与其延迟的乘积可以看做是一个“新码”。这个“新码”也属于戈尔德码家族。仿真结果表明，利用其自相关和交叉相关可找到“新码”的起始点。“新码”的起始点与C/A码的起始点相同。式（5-20）有意思的地方在于它的频率自变性，即频率变化不受其他变量支配，因为f和t都是常量，因而ej2πfτ项是常量，因此，只需搜索“新码”的起始点。尽管这种方法非常吸引人，但它要求输入信号必须是复量，由于收集到的输入电文都是实数，因而必须将它们转化成复数。戈尔德码(GoldCode)与其延迟46为避免将实数信号转化成复数信号的额外计算，

我们可以采用如下方法。

输入信号为

这里，Cs(t)表示卫星s的C/A码。输入信号的延迟表达式为：

(5-23)s(t)与延迟s(t-τ)的乘积是：(5-24)为避免将实数信号转化成复数信号的额外计算，我们可以47

Cn(t)是式(5-21)的定义。式(5-24)含有两项：一个直流项和一个高频项。通常，可以滤去高频。为使式(5-24)可用，|cos(2πfτ)|必须接近单位1。理论上，很难做到这点，因为频率f未知，然而，由于频率在1250±10kHz之内，选择一个延迟时间来满足这个条件是完全可能的。例如，可选择2×π×1250×103τ=π，如此τ=0.4×10-6s=400ns。输入电文用5MHz数字化，采样时间是200ns(1/5MHz)。在这个条件下，|cos(2πfτ)|=|cos(π)|=1。如果频率偏离10kHz，其相关值|cos(2πfτ)|=|cos(2π×1260×103×400×10-6)|=0.9997，非常接近单位1。Cn(t)是式(5-21)的定义。式(5-24)含有两项48因此这个方法可以应用到实际电文中，惟一的约束是不可以像式(5-20)中那样任意选择延迟时间，但可以选择延迟不是很长的0.4μs的倍数，例如τ=1.6μs，当频率偏移±10kHz时，|cos(2πfτ)|=0.995。可以看出，如果对中心频率1250kHz用一个长延迟的话，|cos(2πfτ)|会下降很快。如果延迟过长，|cos(2πfτ)|可能不再接近1。这个方法的问题是当两个带噪声的信号相乘时，噪底会上升，搜索1ms的电文将捕获不到某个卫星。因此这个方法可以应用到实际电文中，惟一的约束是不可以像式49图5-7相位偏移对延迟和累积捕获方法的影响(a)无相位偏移；(b)相位偏移两位

图5-7相位偏移对延迟和累积捕获方法的影响505.4.4长记录电文的相干处理]寻找弱信号的方法通常是增加捕获的电文长度，这样做的好处是提高了信噪比，对此的一个简单解释是，长2ms的电文的FFT产生500Hz的频率分辨率，而1ms的电文的FFT产生1kHz的频率分辨率。由于信号在频谱展开后带宽窄，因而信号强度并不会因为相对窄的频率分辨率而降低，降低频率分辨率带宽反而将噪声降低了一半，因此，信噪比上升了3dB。5.4.4长记录电文的相干处理]51长记录电文相干处理的思想是利用很少的数据点执行快速傅立叶变换(FFT)。利用10ms的电文(或50000个数据点)来说明这个思想。电文的中心频率是1.25MHz，输入信号将被转换成一个基带信号和一个2.5MHz的高频带。如果将高频带用滤波器滤去，将只需处理基带信号。假定此时已经滤去高频信号，基带信号是带C/A码的输入的下变频信号，将这个信号与C/A码点对点相乘。如果到达了C/A码的正确相位，输出将变成连续信号，由多普勒效应引起的最大频率范围在±10kHz。既然这个信号的带宽是20kHz，我们就可以用其2.5倍的频率来采样，即用50kHz采样该信号。10ms的电文在此采样频率下只产生500个数据点。然而，我们是用5MHz采样的，产生了50000个数据点。我们可以每100个数据点取一次平均值作为一个新的数据点，这个取均值的过程等价于一个低通滤波器，因此，它去除了乘法运算后的1.25MHz连续信号中的高频分量，同样，

也去除了接收机收集到的信号中含有的噪声信号。

长记录电文相干处理的思想是利用很少的数据点执行快速傅立叶52这个方法可以用如下10ms的电文的例子来说明：(1)将10ms长的输入信号与本地生成1.25MHz的复合连续信号相乘，并用5MHz采样。因为输出最大频率是10kHz，我们可以将此时的输出称做低频输出。频率近似为2.5MHz的高频分量将被滤波除去，因而可以忽略它。这时，得到了50000个数据点。(2)将这些输出数据与10个C/A码点对点相乘，总共得到50000个数据点。(3)将每100个相邻数据点取其均值，变成一个新的数据点，这个过程滤去了频率近似为2.5MHz的高频分量。

这个方法可以用如下10ms的电文的例子来说明：53(4)对新产生的500个数据点进行FFT变换，找到频率内的最高输出，这个操作产生的频率输出中，只有250个输出有用。(5)将本地生长码相对低频输出移动一个数据点，然后重复步骤(3)和(4)。由于每毫秒产生一次C/A码，我们需要操作5000次，而非50000次。(6)在频域中，总共得到1.250×106(250×5000)个输出，大于预置门限值的最高幅值就是期望值，从这个值可以得到C/A码的起始点和多普勒频移。由此得到的频率分辨率是100Hz。尽管上述方法简单易懂，但循环相关的捕获方法也可以用很少的操作来达到同样的目的。

(4)对新产生的500个数据点进行FFT变换，找到频率545.4.5精频估计的基本概念从1ms电文中得到的频率分辨率大约是1kHz，对跟踪环来说，这个值太粗糙了。适合跟踪过程的频率必须在几十Hz之内。通常，跟踪环的频宽是几Hz，用DFT(或FFT)计算精细频率（FineFrequency）是不可取的。因为为了得到10Hz的分辨率，需100ms的电文记录，如果每毫秒有5000个数据点，100ms就有500000个数据点，其FFT操作将是非常费时的。另外，100ms电文的相位偏移的概率也相对高。

5.4.5精频估计的基本概念55得到适合跟踪环的精细频率分辨率(FineFrequencyResolution)（以后简称精频）的方法是通过相位关系。一旦输入信号中剥离了C/A码，输入就将变成连续波形。如果在m时刻，1ms电文中最高频率分量是Xm(k)(k表示输入信号的频率分量),则可由DFT输出得到输入信号的初始相位θm(k)：

(5-25)得到适合跟踪环的精细频率分辨率(FineFrequen56Im和Re分别表示虚部和实部。让我们假定在m时刻之后很短时间的n时刻，1ms电文的DFT分量Xn(k)也是最强分量，因为输入分量在很短时间内不会迅速变化。n时刻输入信号的初始相位角和频率分量k为：(5-26)这两个相位角可用来计算精频：

(5-27)Im和Re分别表示虚部和实部。让我们假定在m时刻之后很短57这个方程式给出了一个比从DFT得到的要精确得多的频率分辨率。为了保持其值的惟一性，θn(k)-θm(k)的相位差必须小于2π。如果相位差是最大值2π，带宽就是1/(n-m)。这里的(n-m)是两组连续电文之间的延时。这个方程式给出了一个比从DFT得到的要精确得多的频率分辨585.4.6消除精频测量中的不确定性尽管求解精频的基本方法都是基于式(5-27)，但是具体实现它，则有几种略微相异的方法。如果每毫秒取DFT的第k个分量，频率分辨率是1kHz，精频测量的不确定宽度也是1kHz。在图5-8(a)中，给出了5个频率分量，它们之间相隔1kHz。如果输入信号如图5-8(b)所示落在了两个频率分量之间，由于系统中噪声的影响，

将不能确定其相位。

5.4.6消除精频测量中的不确定性59图5-8频域内的不确定区域（a）

频率相离;（b）

输入信号

图5-8频域内的不确定区域60为了去除这种不确定性，有两个解决方法。第一个方法是加速DFT操作过程。如果每0.5ms执行一次DFT操作，不确定带宽是2kHz。在频率分辨率为1kHz和不确定带宽为2kHz的条件下，确定一个精细频率就不存在不确定的问题了。但是这种方法使DFT操作增加了一倍。解决这个问题的第二个方法是，如果输入是连续信号的话，采用一种幅值比较方案，而非加快DFT操作速度。如图5-8(b)所示，输入信号落在两个频率槽之间。假定X(k)的幅值略高于X(k-1)，将X(k)代入式(5-26)和式(5-27)，得到精细分辨率。输入频率与X(k)的值相差接近500Hz。为了去除这种不确定性，有两个解决方法。第一个方法是加速D61正确的结果是输入频率比X(k)低500Hz，由于噪声的存在，我们可能错误地认定是输入频率比X(k)高了500Hz。然而，输入频率X(k)和X(k-1)是非常接近的，且它们比X(k-1)要强得多。因此，如果输入最高频率槽是X(k)，且相位计算在不确定区域范围内，这个区域接近于X(k)与X(k-1)的中心或者X(k)与X(k+1)的中心，这时就比较X(k-1)与X(k+1)两者之间的相位。如果X(k-1)大于X(k+1)，输入频率就是低于X(k)的;否则，就认为输入频率高于X(k)。如此，精频的准确性由相位决定，但是频率差值的符号由最高频率分量邻近的两个频率分量的幅值决定。正确的结果是输入频率比X(k)低500Hz，由于噪声的存62然而，问题比这要略为复杂一些，因为两组连续数据间可能由于导航电文引起了180°的相位偏移。如果发生了这种情况，输入信号就不能被看做连续信号对待了。这种可能性使不确定宽度延迟了1ms，且限制在250Hz。频率相距±250Hz对应的相位角是±π/2。如果频率相距+250Hz，相位角是+π/2，然而导航电文引起的相位偏移π将使相位角变成+π/2(+π/2-π)，相应地引起-250kHz的频率变换。如果在寻找精频时不考虑相位偏移，这个结果将变成500kHz。然而，问题比这要略为复杂一些，因为两组连续数据间可能由于63为了避免这个问题，频率不确定的最大值必须小于250Hz。如果频率最大相差是±200Hz(这个值是实验选取的)，则相应的相位角的差是±2π/5，如图5-9所示。如果存在相位偏移π，相位差了3π/5(|±(2π/5)π|)，这个值要比2π/5大得多。这样来看，相位差可以用来确定精频而不产生错误的频率偏移。如果相位差大于2π/5，就从结果中减去π以保证频率差小于2π/5。为了保持频率在200kHz之内，X(k)中的k值之间的最大距离将是400kHz。如果输入是在两个邻近的k值当中，则输入信号跟两边的k值都相距200kHz。为了避免这个问题，频率不确定的最大值必须小于250Hz。64图5-9相位差小于2π/5将不会造成频率误差

图5-9相位差小于2π/5将不会造成频率误差65图5-10从两个相位角得到的角度差

图5-10从两个相位角得到的角度差66讨论的最后一点是将X(k)的实部和虚部转化成相位角。通常，相位角的范围在±π之间。式(5-26)和式(5-27)中的两个相位角都将用这种方式得到，两个相位角的角度差可以是0～2π间的任何值，如图5-10所示。由于200kHz允许的最大相位角差是2π/5，这两个相位角的角度差必须等于或小于2π/5。如果结果大于2π/5，可以从结果中加上或减去2π，其绝对值就肯定小于2π/5了。如果将噪声考虑进去，2π/5门限值可以适当放宽，如用2.3π/5，意味着其差必须等于或小于2.3π/5这个值。如果这样处理后的值仍然大于这个门限值，就意味着两个连续毫秒数据间有相位偏移了，必须从结果中减去π。当然，得到的角度值也还要经过加减2π，以保证最后得到的值小于门限值。讨论的最后一点是将X(k)的实部和虚部转化成相位角。通67由以上讨论可知，在某一确定卫星中需要经过下面几个步骤来找出它的C/A码起始点和载波频率：(1)对1ms的输入执行循环相关操作，某个确定的C/A码起始点可以从这些循环相关里找到，载波频率可以以1kHz的分辨率得到。(2)找到最高频率分量X(k)，在同一个毫秒数据内的两个分量，一个比X(k)中的k值低400kHz，一个比k值高400kHz，对这两个分量执行DFT操作。三个输出［X(k-1),X(k),X(k+1)］中的最高输出将被指定为新的X(k)，并用这个新的X(k)作为DFT分量来求解精频。由以上讨论可知，在某一确定卫星中需要经过下面几个步骤来找68(3)从C/A码起始点处开始选择连续几个毫秒的数据，我们随意选择为5ms，将这些数据与5组连续C/A码相乘，结果肯定是一个5ms长的连续信号。但是在任意1ms数据中，都有可能含有一个π相位偏移。(4)从所有输入数据中找出Xn(k)，其中n=1,2,3,4,5，然后由式(5-25)找出相位角。角度差值定位为：(5-28)(3)从C/A码起始点处开始选择连续几个毫秒的数据，69(5)角度差的绝对值必须小于门限值(2.3π/5)，如果不能实现，就要从Δθ上加上或减去2π。如果结果仍然大于门限值，就要调整π相位偏移，在Δθ上加上或减去π，得到的值还要与门限值2.3π/5进行比较，来决定是否要再次加减2π。经过这些调整后，最终的角度值就是期望值。(6)式(5-27)用来计算精频。由于有5ms的数据，将得到4个精频值。为了提高精确度，我们将这4个精频的平均值作为要求解的值。

(5)角度差的绝对值必须小于门限值(2.3π/5)，如705.5Matlab与Simulink仿真简介

5.5.1仿真工具简介Matlab是由美国Mathsworks公司开发的大型软件。在Matlab软件中，包括了两大部分：数学计算和工程仿真。其数学计算部分提供了强大的矩阵处理和绘图功能。在工程仿真方面，Matlab提供的软件几乎遍布各个工程领域，并且不断加以完善，现已成为国际上最为流行的科学计算与工程计算软件之一。5.5Matlab与Simulink仿真简介5.5.1715.5.2Simulink的使用本小节将利用BPSK传输作为例子对Simulink进行简单介绍。但是，我们只局限在Simulink的原理及在这个仿真辅导中所要求的特点上，至于更为详尽的内容，建议读者参考Simulink手册。这里给出的Windows界面的图保持了原样。Simulink是Matlab软件的扩展，它是实现动态建模和仿真的一个软件包。它与Matlab语言的主要区别在于，它与用户接口是基于Windows的模型化图形输入的，从而使得用户可以把更多的精力投入到系统模型的构建而非语言的编程上。

5.5.2Simulink的使用72图5-11SimulinkLibraryBrowser窗口

图5-11SimulinkLibraryBrowse73图5-12BPSK传输的Simulink模型

图5-12BPSK传输的Simulink模型74为了创建一个新的Simulink模型，可以通过File→New→Model菜单打开一个新的模型。某个库模块可以从库窗口简单地拖放到Simulink模型窗口中。利用UniformRandomNumberGenerator模块以及随后的Sign模块，可以做图5-12中的二进制信号源。UniformRandomNumber模块在指定的区间上产生均匀分布的随机数，在本例中该随机数在［-1,+1］内产生。跟在其后的Sign模块的输出是：当输入大于0时输出为-1。为了产生输出比特为等概率的+1和-1二进制信号源，可以加入LookupTable模块，它的输入—输出映射关系为：＋1到＋1，-1到-1，而0到+1。为了将模块UniformRandomNumber加到模块中，可以在SimulinkBlockLibraryBrowser中双击模块Sources库(见图5－13)。将模块UniformRandomNumber拖入Simulink模型并双击它。为了创建一个新的Simulink模型，可以通过File→75图5-13Sources库

图5-13Sources库76打开的BlockParameters窗口如图5-14所示，在这个窗口中即可输入如该图所示的有关参数。现将MathOperations库中的Sign模块拖入模型。利用鼠标将模块UniformRandomNumber和Sign连接起来：在模块UniformRandomNumber的输出上按下鼠标左键并拖到Sign模块的输入端，然后释放鼠标键。LookupTable模块从LookupTables库中得到，在模型中双击LookupTable模块并输入如图5-15所示的模块参数。

打开的BlockParameters窗口如图5-14所77图5-14模块参数：UniformRandomNumber图5-14模块参数：UniformRandomNu78图5-15模块参数：LookupTable图5-15模块参数：LookupTable79在开始仿真之前，为了观察源信号，必须调整ConfigurationParameters。在菜单栏中使用命令Simulation→ConfigurationParameters。即可打开如图5-16所示的ConfigurationParameters窗口。就目前而言，仅涉及参数Starttime和Endtime，稍后为了求得准确的结果，可能需要调整用于数值积分的Maxstepsize和Minstepsize。在菜单栏中使用命令Simulation→Start启动仿真。类似地，可以使用Simulation→Stop命令结束或停止仿真。为了观察信号，

双击Scope即可打开如图5-17所示的窗口。

在开始仿真之前，为了观察源信号，必须调整Configur80图5-16仿真参数窗口

图5-16仿真参数窗口81图5-17源信号(Scope)图5-17源信号(Scope)82在窗口中单击鼠标右键会打开一个弹出菜单，通过该菜单可以打开Axesproperties窗口，如图5-18所示，这是为了设定y轴的范围。

图5-18Axesproperties窗口（Scope）

在窗口中单击鼠标右键会打开一个弹出菜单，通过该菜单可83在图5-17中，Scope窗口的时间范围可以通过菜单栏中的Parameters框来改变（见图5-19）。

图5-19Scope的参数窗口

在图5-17中，Scope窗口的时间范围可以通过菜单栏中84图5-12中的其他模块也能用相同的方式包括在内。利用UserDefinedFunctions模块可以产生载波cos(4*pi*u)，该模块将某个指定的表达式（如cos(u)）加到输入变量u上。将模块Fcn拖入仿真模型并双击它，会立即打开BlockParameters窗口，如图5-20所示，在这里可以输入表达式cos(4*pi*u)。图5-12中的其他模块也能用相同的方式包括在内。利用Us85图5-20模块参数：Fcn图5-20模块参数：Fcn86时间变量u可以从Sources库中的Clock模块得到。MathOperations库中的Product模块用来产生已调信号。由于考虑的是无噪声传输仿真，因此就将发送端和接收端直接相连。接收端的匹配滤波器是利用Continuous库中的Integrator模块按积分清除来实现的。复位信号由Sources库中的PulseGenerator模块产生，其参数如图5-21所示。最后的仿真结果如图5-22所示。

时间变量u可以从Sources库中的Clock模块得到。87图5-21模块参数：PulseGenerator图5-21模块参数：PulseGenerator88图5-22宿信号（Scope8）

图5-22宿信号（Scope8）895.6GPS卫星信号捕获的例子

程序中所用到的输入电文在时域内的情形如图5-23所示。输入信号看起来像噪声，其频率描绘可通过FFT求得，如图5-24所示。正如预料的，FFT变换后带宽是2.5MHz，其频谱形状与射频链的滤波器形状相似，通过循环相关后如图5-25所示，其中给出了6号卫星的C/A码起始点，其起始点在2884。图5-26给出了相距1kHz的21个频率分量，最高值发生在k=7处。从图5-25和图5-26我们可以轻易地看出C/A码的起始点和频率。由于输入电文是实际收集到的数据，精频的精确度就难以测定，因为其多普勒频率未知。精频还取决于在下变频中用到的本地振荡器的频率精确性和采样频率的精确性。

5.6GPS卫星信号捕获的例子程序中所用到的输入电文90图5-235000个输入数据点

图5-235000个输入数据点91图5-24输入信号的快速傅立叶变换

图5-24输入信号的快速傅立叶变换92图5-256号卫星的C/A码起始点

图5-256号卫星的C/A码起始点93图5-26卫星6的展开的各个频率分量

图5-26卫星6的展开的各个频率分量94感知计算精频准确性的一个方法是利用数据的不同部分。从输入电文的不同部分计算出六个精频，用到的电文部分是：1～25000，5001～30001，10001～35001，15001～40001，20001～45001，25001～50001。这些电文都是5ms长，且起始点偏离1ms，在相邻两组间，有4/5的数据相同。因此，得到的五个精频是非常接近的，这六个精频之间相差是-2.4，9.0，-8.2，5.4和2.3Hz。数据都是从固定的接收机上收集到的，因此，每毫秒刷新一次频率的变化是非常小的。从而，这些频差可以看做是捕获方法中的误差。当信号长度改变时，精频的差也跟着变化，对一个弱信号来说，如果使用相同的捕获方法，其精频可相差几十赫兹。感知计算精频准确性的一个方法是利用数据的不同部分。从输95对弱信号（24号卫星）执行捕获，如图5-27和图5-28所示。从图中我们很难判断C/A码的起始点和频率是否为正确值，需要通过另外的相关运算才能找出它们，它的正确值与我们在图5-27和图5-28中看到的不同。

对弱信号（24号卫星）执行捕获，如图5-27和图5-2896图5-2724号卫星的C/A码起始点

图5-2724号卫星的C/A码起始点97图5-2824号卫星的展开信号的频率分量

图5-2824号卫星的展开信号的频率分量985.7关于捕获的一些子程序

5.7.1随机编码过程仿真打开Matlab，新建一个MFile，输入以下的文本，保存为GenerateCode.m，并调试运行。可以单步或多步运行M文件。%该程序仿真产生三组伪随机数，模拟随机编码过程closeall;clear;clc;t=0:5:1000; %时间变量为1s，步长5ms%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5.7关于捕获的一些子程序5.7.1随机编码过程仿真99x1=rand(1,length(t));%生成一组1*length(t)维的伪随机数x1=round(x1); %将x1数组中的元素转换到与其最近的整数index1=find(x1==0);%找出x1数组中为0的元素x1(index1)=-ones(1,length(index1));%返回一个index1长度的一维数组x2=rand(1,length(t));x2=round(x2);index1=find(x2==0);x2(index1)=-ones(1,length(index1));x3=rand(1,length(t));x3=round(x3);index1=find(x3==0);x3(index1)=-ones(1,length(index1));x1=rand(1,length(t));%生成一组100y1=zeros(1,1);y2=zeros(1,1);y3=zeros(1,1);%生成三个1*1的零矩阵，存放伪随机值t1=zeros(1,1);t2=zeros(1,1);t3=zeros(1,1);%生成三个1*1的零矩阵，存放时间值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%将生成的伪随机值存放到对应的矩阵中%%%%%%y1(1)=x1(1);y2(1)=x2(1);y3(1)=x3(1);t1(1)=t(1);t2(1)=t(1);t3(1)=t(1);k=2;m=2;n=2;fori=2:length(t)if(x1(i)==x1(i-1))y1(k)=x1(i);t1(k)=t(i);k=k+1;y1=zeros(1,1);y2=zeros(1,1);101elsey1(k)=x1(i-1);y1(k+1)=x1(i);t1(k)=t(i);t1(k+1)=t(i)+0.01;%如果该时刻阶跃变化，则右移0.01来记录k=k+2;endif(x2(i)==x2(i-1))y2(m)=x2(i);t2(m)=t(i);m=m+1;elsey2(m)=x2(i-1);y2(m+1)=x2(i);t2(m)=t(i);t2(m+1)=t(i)+0.01;else102m=m+2;endif(x3(i)==x3(i-1))y3(n)=x3(i);t3(n)=t(i);n=n+1;elsey3(n)=x3(i-1);y3(n+1)=x3(i);t3(n)=t(i);t3(n+1)=t(i)+0.01;n=n+2;endendm=m+2;103%%%%%%画图说明产生的三组伪随机码%%%%%subplot(3,1,1);plot(t1,y1,′y′);gridon;axis(［-1,1001,-1.5,1.5］);subplot(3,1,2);plot(t2,y2,′y′);gridon;axis(［-1,1001,-1.5,1.5］);subplot(3,1,3);plot(t3,y3,′y′);gridon;axis(［-1,1001,-1.5,1.5］);%%%%%%画图说明产生的三组伪随机码%%%%%104whitebg(′black′);%将三组伪随机值存放到新矩阵并传送到Workspace中，方便以后调用数据%simin1=［t1;y1］′;simin2=［t2;y2］′;simin3=［t3;y3］′;whitebg(′black′);105输入上面的文本后，运行得到如图5-29所示的伪随机码。其结果还可以通过Workspace查看，之所以要传送到Workspace中，

是因为我们在后面的仿真中还要用到这些伪随机码。

输入上面的文本后，运行得到如图5-29所示的伪随机码。106图5.29模拟三组伪随机码的仿真示意图

图5.29模拟三组伪随机码的仿真示意图1075.7.2获取导航信息的仿真在这个仿真中，我们要得到卫星的导航信息，其中忽略了频率偏移，且卫星信号的数据是基于前面章节中得到的卫星数据，存放在GPSsignal.mat文件中。在运行程序之前，我们必须将GPSsignal.mat中的数据传送到Workspace中。

5.7.2获取导航信息的仿真108clear;%清除所有的变量clc;%关闭所有指令窗closeall;%关闭所有未隐藏的窗口loadGPSsignal.mat%从GPSsignal.mat中获得所有变量到Workspace中 %n=length(GPSsignal);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%SvNum=12;%调用fGenerateCAcode3.m文件，获得卫星编号为12的C/A码Temp=fGenerateCAcode3(SvNum);index1=find(Temp==0);%找出C/A码中的低电平，形成列向量存放在index1中Temp(index1)=-ones(1,length(index1));%返回一个与index1相同长度的一维数组

clear;%清除所有的变量109SinWave=sin(［0:2*pi/8:2*pi*7/8］);%产生一个步进2*pi/8的正弦波SinWave=single(SinWave);%将SinWave矩阵转换成单精度矩阵GpsMatch=zeros(1,1);%生成一个1*1的零矩阵SinWave=［SinWaveSinWaveSinWaveSinWaveSinWave］;GPSsignal.mat %生成一个5ms的连续信号fori=1:length(Temp)GpsMatch=［GpsMatchTemp(1,i)*SinWave］;endGpsMatch=GpsMatch(2:length(GpsMatch));%生成本地匹配信号并去掉第一个元素n=length(GpsMatch);m=50000;fori=1:mRes(i)=GpsMatch*GPSsignal(1,i:i+n-1)′;%将本地生成匹配信号GpsMatch与GPSsignal循环卷积

SinWave=sin(［0:2*pi/8:2*pi*7/8110endplot(1:m,Res);%图示卷积的幅值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Res2=abs(Res);fori=1:100%该循环目的是找出Res中的最大值［CI］=max(Res2);Res2(1,I)=0;Index(1,i)=I;end;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%下面是要找到GPSsignal中的导航电文%

end111%方法是检验GPSsignal(已被C/Acode调制过)的每个最大值处是否具有相位偏移%w=length(GpsMatch);m=1;fori=Index(1,1):w:(length(GPSsignal)-w+1)NavigationBit(m)=(GPSsignal(i:i+w-1)*GpsMatch′)/w;m=m+1;endNavigationCode=zeros(1,1);NavigationCode(1,1)=NavigationBit(1,1);m=2;Count=0;fori=2:length(NavigationBit)if(NavigationBit(1,i)～=NavigationBit(1,i-1))NavigationCode(1,m)=NavigationBit(1,i);%方法是检验GPSsignal(已被C/Acode调制过112 m=m+1;Count=0;elseCount=Count+1;if(Count>=5)NavigationCode(1,m)=NavigationBit(1,i);m=m+1;Count=0;endendend m=m+1;113%%%%%我们已得到导航码NavigationCode,下面将它转化成数字%%%%NavigationCode=NavigationCode./abs(NavigationCode);index1=find(NavigationCode<0);NavigationCode(index1)=zeros(1,length(index1));Table1=［0000;0001;0010;0011;0100;0101;0110;0111;1000;1001;1010;%%%%%我们已得到导航码NavigationCode,下面1141011;1100;1101;1110;1111;］;Result=″;fori=1:4:length(NavigationCode)TT=NavigationCode(i:i+3);MatchTable=Table1(:,1:4)-［TT;TT;TT;TT;TT;TT;TT;TT; TT;TT;TT;TT;TT;TT;TT;TT］;MatchTable=sum(abs(MatchTable)′);TempChar=′0′;BestMatch=find(MatchTable==min(MatchTable));1011;115%用查表的方式将导航码NavigationCode转换成十六进制数switch(BestMatch)case1TempChar=′0′;case2TempChar=′1′;case3TempChar=′2′;case4TempChar=′3′;case5TempChar=′4′;