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

1、第5章 GPS卫星信号的捕获 5.1 概述 5.2 GPS卫星信号的多普勒效应5.3 GPS卫星信号捕获的考虑5.4 GPS卫星信号的捕获方法 5.5 Matlab与Simulink仿真简介 5.6 GPS卫星信号捕获的例子 5.7 关于捕获的一些子程序 5.1 概 述 为了跟踪和解码GPS信号, 首先要捕获到GPS信号。将捕获到的GPS信号的必要参数立刻传递给跟踪过程，再通过跟踪过程便可得到卫星的导航电文。GPS卫星处于高速运动中，因此，其频率会产生多普勒频移。载波频率与C/A码的多普勒频移将在下面详细讲述。为覆盖高速卫星预期中的所有多普勒频率范围，捕获方法覆盖的频率范围必须在10 kHz之

2、内。一旦捕获到GPS信号，立刻去测量两个重要参数： C/A码的起始点和载波频率（因为多普勒频移而变化）。接收机接收到的一系列数据往往包含多个卫星信号，每个信号具有不同的C/A码的不同起始点和不同的多普勒频率。 针对某个特定的卫星信号， 捕获过程就是要找到C/A码的起始点，并利用找到的起始点展开C/A码频谱，一旦复现了C/A码的频谱， 输出信号将变成连续波（Continuous Wave， CW），于是便得到其载波频率。也就是说，捕获过程就是要获得输入信号的C/A码的起始点和载波频率，然后传递给跟踪过程。 捕获与跟踪过程所用到的数据都是从原始的卫星信号经过下变频器（即与中频混频）之后收集到的，其

3、中频（IF）为21.25 MHz, 采样频率为5 MHz， 信号的中心频率为1.25 MHz。这组数据来源于我们设计的卫星信号接收装置，但是用模拟仿真算法也可以产生类似的数据，以验证本章介绍的卫星信号的捕获算法。 5.2 GPS卫星信号的多普勒效应GPS卫星绕地球旋转一周的时间是11 h 58 min 2.05 s， 从卫星轨道的近似半径可以求得卫星的角速度d/dt和运动速度vs： (5-1) 这里rs表卫星轨道的平均半径。我们知道，一个太阳日和一个恒星日之间相差3 min55.91 s，在这段时间里，卫星大约运行了914 km(3874 m/s235.91 s)，对应地球表面与卫星的最高点，

4、相应的角度近似为0.045 rad(914/20 192)或2.6。如果卫星接近地平线，相应的角度为0.035 rad或者2。 因此我们可以看出，对于地球表面的固定一点，在每天的同一时间里， 卫星位置大约改变22.6。 如图5-1所示，卫星在位置S处，用户在位置A处，相对用户的卫星角速度vd造成了多普勒频移， 其中vd值为： (5-2) 根据卫星轨道速度，取水平方向的最大值， 得到多普勒角速度的最大值为： (5-3) 图5-1 卫星运动引起的多普勒频移 这个值与高速航空卫星的速度相等。由地面设施引起的多普勒频移通常非常小，即使它直接相对卫星运动引起了最大多普勒效应。对由C/A码调制过的频率L1

5、(f=1575.42 MHz），最大的频移为： (5-4) 这里c是光速。 因而， 对于一固定观测器来说，最大的多普勒频移是5 kHz。 如果地面GPS接收设备高速移动，就要考虑进多普勒效应了。要单独产生一个5 kHz的频率范围，接收设备需以2078 m/h的速度相对卫星运动，这个速度将包含最高速度的航空器。 因而，在设计GPS接收机时，如果接收机用在低速设备，则认为载波频率的多普勒频移范围在5 kHz；相反地，如果用在高速设备，就要合理假定其多普勒频移的范围在10 kHz。 这些值对于确定捕获过程的搜索频率范围是至关重要的。因此，就有了引言中我们提到的为覆盖高速卫星预期中的所有多普勒频率范围

6、，捕获方法覆盖的频率范围必须在10 kHz之内之说。 由于C/A码频率非常低，所以C/A码上的多普勒频移也是非常小的。C/A码的频率是1.023 MHz，比载波频率低1540(1575.42/1.023)倍。 C/A码的多普勒频率为： (5-5) 如果接收机高速移动，这个值就得取双倍，为6.4 Hz，得到的这个值对跟踪过程（BASS方法）来说是至关重要的，这点我们将在下一章讨论。在BASS跟踪过程中，输入信号与本地产生数据必须要紧密对齐，而C/A码的多普勒频移可能造成输入码与本地生成码的不对齐。 在数字化卫星信号中，如果数据用5 MHz采样(称为采样频率)， 则每个采样之间相隔200 ns（称

7、为采样时间）。 在跟踪过程中，我们期望本地生成信号与输入信号未对齐的长度在半个采样时间或近似100 ns之内，若两个信号之间相差超过这个长度，跟踪将失锁，即失去跟踪灵敏性。C/A码的基波时间为977.5 ns或1/1.023106 s，它移动一个周期需花费156.3 ms（1/6.4）。对这句话的理解是：多普勒频移使C/A码的频率由1.023106变为JP21.023106+6.4，即每秒内多变化了6.4个周期，则多变化一个周期的时间为(1/6.4) ms。 因此，它移动100ns的数据长度，近似花费16 ms(100156.3/977.5)。在高速导航器中，需每16 ms选择一批数据，以保证

8、输入信号与本地生成码更好地匹配。 但是由于输入信号中有噪声，用1 ms的数据可能不能非常精确地匹配， 因而可以拓宽输入信号的调整时间，每20 ms选择一批输入信号。 对一个低速导航器来说，这个时间可拓宽到40 ms。 从以上讨论我们知道，输入信号的调整取决于采样频率。 采样频率高将缩短调整时间，因为采样时间变短，在半个周期内匹配输入码与本地生成码的时间变短。如果输入信号的强弱和跟踪灵敏度不成问题，则输入信号的调整时间可以拓宽。然而，输入信号与本地生成信号相离必须在半个基波或488.75 ns(977.5/2)之内，这个时间看做是输入信号与本地生成信号相离允许的最大值。C/A码的多普勒频率是6.

9、4 Hz， 调整时间可以拓宽到78.15 ms（1/26.4）。 5.3 GPS卫星信号捕获的考虑 5.3.1 捕获时的最大电文长度 C/A码长1 ms，那么至少要用1ms的电文来捕获，甚至只用1 ms的电文来捕获时，都可能发生导航电文相位偏移。如果电文中有数据偏移，电文的下一个1ms将不会含有偏移。因此，为了保证捕获电文中不含有数据偏移，需要用两组连续的电文来捕获， 这个电文最大长度是10 ms。如果使用两组连续的10 ms电文来捕获， 就保证了在某一组电文中不含相位偏移。 限制电文长度的第二个因素是C/A码的多普勒效应。如果理想相关峰值是1，当C/A码超前或滞后1/2码时，相关的峰值将降为

10、0.5，对应地，幅值下降6 dB。如上节所讨论的，假定未对准的C/A码允许在半个基波（0.489 s）内, 基波频率是1.023 MHz，C/A码上预期的最大多普勒频移是6.4 Hz， 那么两个相差6.4 Hz的频率交换半个基波需要78 ms（1/26.4）， 这个数据远大于10 ms。 5.3.2 捕获中的频率步长 捕获时的另一个考虑因素是捕获中剥离载波所需的频率。 在前文曾提过，需被覆盖的多普勒频率范围是10 kHz， 决定覆盖这20 kHz范围的步进频率非常重要。步进频率的大小与捕获中的电文长度紧密相关。如果输入信号与本地混频信号相距1个周期，它们两者没有相关性；如果两者小于1个周期，

11、则它们有部分相关性。随意选取两个信号允许的最大频率相离是0.5个周期， 如果电文记录长度为1 ms，一个1 kHz的信号将在1 ms内变换一个周期，为了保持1 ms内的最大频率相离在0.5个周期，步进频率须为1 kHz。在这种情况下，在输入信号与卷积信号之间最大限度的频率相离是500 Hz/s或0.5 Hz/ms， 且输入信号刚好在两个频率点之间。如果电文记录长度为10 ms，100 Hz的搜索步进频率就可以满足要求。看待这个问题的一个简单方法是频率相离与电文长度成反向关系。 由以上的讨论可以得出如下结论：当捕获输入电文为1 ms长时，步进频率是1kHz；当电文为10 ms时，步进频率是100

12、 Hz。从这个简单的讨论可以明显看出，捕获中的操作执行次数与总的数据点并不是成线性比例关系，当电文长度从1 ms上升到10 ms时，电文长度上升10倍，频率点数目也上升了10倍， 捕获所需的操作次数上升了不止10倍。因此，实际捕获时， 如果强调捕获的速度的话，电文长度需保持在最小值。 执行次数的增加，取决于实际捕获所用的方法，下面将讨论具体的捕获方法。 5.4 GPS卫星信号的捕获方法5.4.1 传统捕获方法 捕获的基本思想就是展开输入信号，找到载波频率。如果相位正确的C/A码与输入信号相乘，输入信号将如图5-2中所示变成一个连续信号。图5-2中最上端是输入信号，这个信号是C/A码编码的射频信

13、号相位，注意到射频信号与C/A码是任意选择的，它们并不代表卫星传输过来的真正的信号。中间的图表示C/A码，其值在1之间。 最下端的图是连续信号，代表输入信号与C/A码的乘积，其对应的频谱不再是扩频， 而是一个连续信号，这个过程有时称做输入信号的C/A码剥离。 图5-2 C/A编码的输入信号与本地C/A码相乘 一旦输入变成连续信号，通过傅立叶变换可以知晓其频率。 如果输入电文长是1 ms，则傅立叶变换的频率分辨率为1 kHz。 可以设置一个阈值，即一个门限来判断一个频率分量是否足够大。高于这个门限的最高频率分量就是我们要找的频率。 如果用5 MHz使输入信号数字化，1 ms的电文将含有5000个

14、数据点， 一个5000点的FFT(快速傅立叶变换)将产生5000个频率分量。 然而， 5000个频率分量中只有前2500个频率分量包含有用信息， 后2500个是前2500个频率分量的复共轭。频率分辨率是1 kHz， 因此， FFT覆盖的总频率范围是2.5 MHz，是采样频率的一半。 但我们感兴趣的频率范围只有20 kHz，而非2.5 MHz，因此， 为了节省计算时间， 可以用DFT(离散傅立叶)方法， 只计算相隔1 kHz的21个频率分量。 由于输入信号中C/A码的起始点还未知，因此，必须先找到这个点。为了找到这个点，我们产生本地的一个C/A码，将其数字化为5000个点，且与输入信号点对点相乘

15、，对此乘积的结果进行FFT或DFT，可以找到其频率。为了搜索1 ms的电文， 输入电文与本地产生数字化的C/A码必须彼此相对移动5000次。 如果使用FFT，需操作5000次，每次操作包含5000个点对点的乘积和一个5000点的FFT。由于只有2500个频率分量提供信息， 其余2500个分量提供的是多余的信息，因此，输出的5000点数据每个含有2500个频率分量，这样，在频域共有1.25107(50002500)个输出。 在这1.25107个输出中最高幅值如果大于预先设定的门限值，则该分量就是期望值。在如此多的数据中寻找最高频率分量也是非常费时的。由于在期望的20 kHz的FFT输出中，只有2

16、1个频率分量是我们感兴趣的， 总输出可以降低到105 000(500021)个。通过这种方法，可以用时间分辨率为200 ns(1/5 MHz)和1 kHz的频率分辨率找出C/A码的起始点。 如果使用10 ms的电文，因为信号只需要卷积1 ms，因此需要5000次操作。每次操作包含50 000个点对点乘积和50 000次FFT，总共有1.25108(500025 000)个输出结果。如果在预期的20 kHz范围内只考虑201个频率分量，需从1 005 000(5000201)个结果中挑选，因此可以看出，从1 ms到10 ms操作的增长是非常显著的。 寻找C/A码起始点的时间分辨率仍是200 ns

17、， 但频率分辨率提高到了100 Hz。 传统捕获方法是由硬件实现的，硬件主要用以执行以上所讲的处理。假定输入信号用5 MHz数字化，一个可行方法是产生5000点的数字C/A码，并与输入信号点对点相乘，每200 ns执行一次5000点相乘，每200 ns对一个5000点相乘的结果进行FFT频率分析，图5-3示意了该过程。如果C/A码与输入电文匹配， FFT结果中将含有一个强分量，如前所述，这个方法将产生 1.25107(50002500)个输出结果，然而，只有在10 kHz频率范围内的输出结果中挑选，这个约束使得排序过程变简单了。 图5-3 用C/A码和频率分析进行捕获 完成这个操作的另一个方法

18、是通过离散傅立叶变换(DFT), 将本地产生码变换成含有C/A码和射频信号的电文。射频是复数，可用ejwt表示。本地码信号由射频与C/A码乘积获得，因此， 本地码也是复数。假定L1频率(1575.42 MHz)被转换到21.25 MHz，并用5 MHz数字化，输出频率将会是1.25 MHz; 同时假定捕获过程以1 kHz的步长搜索125010 kHz 的频率范围，且总共有21个频率分量，本地码lsi可用下式表示： (5-6) 其中，s表示卫星的编号；i=1, 2, 3, , 21；Cs表示编号为s的卫星的C/A码；fi=1250-10,1250-9,1250-8,1250+10 kHz。 本地

19、信号须用5 MHz数字化，产生5000个数据点， 这21组数据代表相距1 kHz的21个频率。这些数据与输入信号有关，如果本地产生信号包含正确的C/A码和正确的频率分量，当遇到正确的C/A码相位时，输出结果将会特别高。 图5-4 通过本地产生C/A码和射频码捕获 5.4.2 循环相关捕获方法 1. 循环卷积与循环相关 在讲述循环相关捕获方法之前， 首先简单介绍执行相关操作的主要数学运算。 如果一个输入信号经过一个线性时不变系统，可由卷积运算得到时域结果，也可由傅立叶变换得到频域结果。如果系统脉冲响应是h(t)，输入信号x(t)可通过卷积产生如下的输出y(t)： (5-7) 通过傅立叶变换，y(

20、t)的频域响应为： (5-8) 令变量t-=u，那么 (5-9) 为得到时域输出，需对Y(f)进行反傅立叶变换， 结果可写成： (5-10) 此处，*表示卷积，F-1表示傅立叶反变换。 可以发现一个简单的关系，时域的卷积等效于频域的傅立叶变换，或者说频域卷积运算等效于时域傅立叶变换， 这两种关系可以写成： (5-11) 通常称其为傅立叶变换中卷积的二元性。 这个思想还可用于离散系统。与连续时域表达式不同，响应y(n)可以表达成： (5-12) 此处，x(m)表示输入信号，h(n-m)是离散时域系统响应。注意到，该式子中，h(n-m)中的时移是循环的，因为离散操作是周期性的。对上式进行离散傅立叶

21、变换(DFT)， 结果为： (5-13) 然而，捕获时用到的并不是卷积运算，而是相关运算， 它与卷积运算不同。x(n)与h(n)之间的相关表示为： (5-14) 上式与式(5-12)的惟一不同点就在于h(n+m)中的m前面的符号不同。h(n)是另一信号， 而非线性系统的脉冲响应。 z(n)的DFT变换为： (5-15) 这里X-1(k)表示离散傅立叶反变换。上式还可写成： (5-16) 如果x(n)是实数，x(n)*=x(n)，其中*表示复共轭。 有了这个关系， Z(k)的幅值为： (5-17) 利用上式的这个关系，可得到输入信号与本地产生信号的相关性。如前所述，这个式子提供了循环相关性。 2

22、. 循环相关捕获的方法 上面讨论了循环相关的运算法则，利用这些运算法则就可以执行捕获了。然而实际捕获中， 输入电文不是连续到达接收机的，因此循环相关操作适合一组或一批电文。输入电文经过5 MHz的ADC(模拟/数字转换器)采样存储在存储器中，只有1 ms的输入电文用来寻找C/A码的起始点，其搜索频率分辨率即步进频率是1 MHz。 捕获输入电文， 需经过以下几个步骤： (1) 对1ms的输入电文x(n)进行快速傅立叶变换FFT， 将输入转换到频域，值为X(k)，这里， n=k=0, 1, 2, , 4999。 (2) 取X(k)的复共轭，值为X(k)*。 (3) 利用式(5-6)，产生21个本地

23、码lsi(n)(i=1, 2, , 21)，本地码包含了卫星s的C/A码和一个复射频信号的乘积， 然后用5 MHz采样本地码，本地码的频率fi相距1MHz。 (4) 对lsi(n)进行FFT，变换到频域中，值为Lsi(k)。 (5) 将X(k)*与Lsi(k)点对点相乘，结果为Rsi(k)。 (6) 对Rsi(k)进行反傅立叶变换，变换到时域中的值为rsi(n)，得到其绝对值|rsi(n)|，总共有105 000(500021)个|rsi(n)|。 (7) 在输入电文200ns的时间分辨率和载波频率为1 kHz分辨率的条件下，|rsi(n)|最大值中的第n位和第i个载波频率芯给出了C/A码的初

24、始点。 以上的操作可由图5-5说明。 图5-5所示结果是时域中的， 且只给出21个本地码的其中一个。可认为输入电文与本地电文位于两个圆柱体表面，为了去匹配输入电文， 本地码要旋转5000次。换句话说，一个圆柱体相对于另一个圆柱体旋转5000次。在每一步，5000个输入电文与5000个本地电文点对点相乘，相乘结果加到一起。 包含本地码与输入码所有可能的乘积需5000步，乘积中最高幅值将被记录下。最高幅值如果同时大于预置门限值的话， 就是我们的期望值。 图5-5 循环相关捕获示意图 3. 循环相关捕获的改进方法 循环相关捕获的改进方法与上面介绍的捕获方法的惟一区别就是傅立叶变换FFT的长度减小了一

25、半。上一节的循环相关方法的第(3)步产生本地码lsi(n)，由于lsi(n)是复量，故其频谱是不对称的， 如图5-6所示。 图5-6 本地产生信号频谱图 从图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)

26、， 每个lsi(n)有5000个数据点。 (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) |，总共有52 500(250021)个| rsi(n) |。 (7) |rsi(n) |的最大值如果大于预置门限值的话就是我们的期望值。第i个频率就是分辨率1kHz的载波频率，第n位就是分辨率400ns的C/A码的初

27、始点。 (8) 由于这种方法寻找C/A码起始点的时间分辨率是400 ns，可以通过比较第n位的幅值与(n-1)和(n+1)位幅值， 将分辨率提高到200ns。 该方法中，从第(5)步到第(7)步，只有2500个数据点的操作，而非5000个。第(7)步挑选最大幅值的过程由于只用到一半的输出而简化了。因此这种方法节省了操作时间。仿真结果表明，这种方法的信噪比略低于常规循环相关捕获方法1.1dB，这可能是由于丢弃频域内后半部分信号引起的。 5.4.3 延迟与累积捕获方法 基于延迟与累积的捕获方法是为了去除输入信号中的频率信息。去除了频率信息，只需用本地C/A产生码去寻找输入信号C/A码的起始点，一旦

28、找到其起始点，就可由FFT或DFT得到其频率。从理论角度看，这是个非常吸引人的方法，然而，运用到实际中处理GPS信号，还需进一步研究。这种方法讨论如下。 首先，我们假定输入信号s(t)是复量，因此表示为： (5-18) 这里，Cs(t)表示卫星s的C/A码。 输入信号的延迟表达式为： (5-19) 表示延迟时间，s(t)与延迟表达式的复共轭的乘积为： (5-20) 这里的 (5-21) 戈尔德码(Gold Code)与其延迟的乘积可以看做是一个“新码”。这个“新码”也属于戈尔德码家族。仿真结果表明，利用其自相关和交叉相关可找到“新码”的起始点。“新码”的起始点与C/A码的起始点相同。式（5-2

29、0）有意思的地方在于它的频率自变性，即频率变化不受其他变量支配，因为f和t都是常量， 因而ej2f项是常量，因此，只需搜索“新码”的起始点。尽管这种方法非常吸引人，但它要求输入信号必须是复量，由于收集到的输入电文都是实数， 因而必须将它们转化成复数。 为避免将实数信号转化成复数信号的额外计算， 我们可以采用如下方法。 输入信号为 这里，Cs(t)表示卫星s的C/A码。输入信号的延迟表达式为： (5-23) s(t)与延迟s(t-)的乘积是： (5-24) Cn(t)是式(5-21)的定义。式(5-24)含有两项：一个直流项和一个高频项。 通常，可以滤去高频。为使式(5-24)可用， |cos(

30、2f)|必须接近单位1。理论上，很难做到这点，因为频率f未知，然而，由于频率在125010 kHz之内，选择一个延迟时间来满足这个条件是完全可能的。例如，可选择21250103=，如此=0.410-6s=400 ns。输入电文用5MHz数字化，采样时间是200 ns(1/5 MHz)。 在这个条件下，|cos(2f)|=|cos()|=1。如果频率偏离10 kHz，其相关值|cos(2f)|=|cos(2126010340010-6)|=0.9997， 非常接近单位1。 因此这个方法可以应用到实际电文中，惟一的约束是不可以像式(5-20)中那样任意选择延迟时间，但可以选择延迟不是很长的0.4s

31、的倍数，例如=1.6 s，当频率偏移10kHz时，|cos(2f)|=0.995。可以看出，如果对中心频率1250kHz用一个长延迟的话，|cos(2f)|会下降很快。 如果延迟过长，|cos(2f)|可能不再接近1。 这个方法的问题是当两个带噪声的信号相乘时，噪底会上升，搜索1ms的电文将捕获不到某个卫星。 图 5-7 相位偏移对延迟和累积捕获方法的影响(a) 无相位偏移； (b) 相位偏移两位 5.4.4 长记录电文的相干处理寻找弱信号的方法通常是增加捕获的电文长度，这样做的好处是提高了信噪比， 对此的一个简单解释是，长2 ms的电文的FFT产生500 Hz的频率分辨率，而1 ms的电文的

32、FFT产生1 kHz的频率分辨率。由于信号在频谱展开后带宽窄，因而信号强度并不会因为相对窄的频率分辨率而降低，降低频率分辨率带宽反而将噪声降低了一半，因此，信噪比上升了3 dB。 长记录电文相干处理的思想是利用很少的数据点执行快速傅立叶变换(FFT)。利用10ms的电文(或50 000个数据点)来说明这个思想。电文的中心频率是1.25 MHz，输入信号将被转换成一个基带信号和一个2.5 MHz的高频带。如果将高频带用滤波器滤去，将只需处理基带信号。假定此时已经滤去高频信号， 基带信号是带C/A码的输入的下变频信号，将这个信号与C/A码点对点相乘。如果到达了C/A码的正确相位，输出将变成连续信号

33、，由多普勒效应引起的最大频率范围在10 kHz。既然这个信号的带宽是20 kHz，我们就可以用其2.5倍的频率来采样， 即用50 kHz采样该信号。10 ms的电文在此采样频率下只产生500个数据点。然而，我们是用5 MHz采样的， 产生了50 000个数据点。我们可以每100个数据点取一次平均值作为一个新的数据点，这个取均值的过程等价于一个低通滤波器，因此， 它去除了乘法运算后的1.25 MHz连续信号中的高频分量，同样， 也去除了接收机收集到的信号中含有的噪声信号。 这个方法可以用如下10ms的电文的例子来说明： (1) 将10ms长的输入信号与本地生成1.25 MHz的复合连续信号相乘，

34、并用5 MHz采样。因为输出最大频率是10 kHz，我们可以将此时的输出称做低频输出。频率近似为2.5 MHz的高频分量将被滤波除去，因而可以忽略它。这时，得到了50 000个数据点。 (2) 将这些输出数据与10个C/A码点对点相乘，总共得到50 000个数据点。 (3) 将每100个相邻数据点取其均值，变成一个新的数据点，这个过程滤去了频率近似为2.5 MHz的高频分量。 (4) 对新产生的500个数据点进行FFT变换，找到频率内的最高输出，这个操作产生的频率输出中，只有250个输出有用。 (5) 将本地生长码相对低频输出移动一个数据点，然后重复步骤(3)和(4)。由于每毫秒产生一次C/A

35、码，我们需要操作5000次，而非50 000次。 (6) 在频域中，总共得到1.250106(2505000)个输出， 大于预置门限值的最高幅值就是期望值，从这个值可以得到C/A码的起始点和多普勒频移。 由此得到的频率分辨率是100 Hz。 尽管上述方法简单易懂，但循环相关的捕获方法也可以用很少的操作来达到同样的目的。 5.4.5 精频估计的基本概念 从1 ms电文中得到的频率分辨率大约是1 kHz，对跟踪环来说，这个值太粗糙了。适合跟踪过程的频率必须在几十Hz之内。通常，跟踪环的频宽是几Hz，用DFT(或FFT)计算精细频率（Fine Frequency）是不可取的。因为为了得到10 Hz的

36、分辨率，需100 ms的电文记录，如果每毫秒有5000个数据点，100 ms就有500 000个数据点，其FFT操作将是非常费时的。另外， 100 ms电文的相位偏移的概率也相对高。 得到适合跟踪环的精细频率分辨率(Fine Frequency Resolution)（以后简称精频）的方法是通过相位关系。 一旦输入信号中剥离了C/A码，输入就将变成连续波形。如果在m时刻，1ms电文中最高频率分量是Xm(k)(k表示输入信号的频率分量), 则可由DFT输出得到输入信号的初始相位m(k)： (5-25) Im和Re分别表示虚部和实部。让我们假定在m时刻之后很短时间的n时刻，1ms电文的DFT分量X

37、n(k)也是最强分量，因为输入分量在很短时间内不会迅速变化。n时刻输入信号的初始相位角和频率分量k为： (5-26) 这两个相位角可用来计算精频： (5-27) 这个方程式给出了一个比从DFT得到的要精确得多的频率分辨率。为了保持其值的惟一性，n(k)-m(k)的相位差必须小于2。如果相位差是最大值2，带宽就是1/(n-m)。这里的(n-m)是两组连续电文之间的延时。 5.4.6 消除精频测量中的不确定性尽管求解精频的基本方法都是基于式(5-27)，但是具体实现它，则有几种略微相异的方法。如果每毫秒取DFT的第k个分量，频率分辨率是1 kHz，精频测量的不确定宽度也是1 kHz。 在图5-8(

38、a)中，给出了5个频率分量，它们之间相隔1 kHz。 如果输入信号如图5-8(b)所示落在了两个频率分量之间，由于系统中噪声的影响， 将不能确定其相位。 图5-8 频域内的不确定区域（a） 频率相离; （b） 输入信号 为了去除这种不确定性，有两个解决方法。第一个方法是加速DFT操作过程。如果每0.5 ms执行一次DFT操作，不确定带宽是2 kHz。在频率分辨率为1 kHz和不确定带宽为2 kHz的条件下，确定一个精细频率就不存在不确定的问题了。 但是这种方法使DFT操作增加了一倍。解决这个问题的第二个方法是， 如果输入是连续信号的话，采用一种幅值比较方案，而非加快DFT操作速度。 如图5-8

39、(b)所示，输入信号落在两个频率槽之间。假定X(k)的幅值略高于X(k-1)，将X(k)代入式(5-26)和式(5-27)，得到精细分辨率。输入频率与X(k)的值相差接近500 Hz。 正确的结果是输入频率比X(k)低500Hz，由于噪声的存在， 我们可能错误地认定是输入频率比X(k)高了500 Hz。然而， 输入频率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)大于

40、X(k+1)，输入频率就是低于X(k)的;否则，就认为输入频率高于X(k)。如此，精频的准确性由相位决定，但是频率差值的符号由最高频率分量邻近的两个频率分量的幅值决定。 然而，问题比这要略为复杂一些，因为两组连续数据间可能由于导航电文引起了180的相位偏移。如果发生了这种情况，输入信号就不能被看做连续信号对待了。这种可能性使不确定宽度延迟了1ms，且限制在250Hz。频率相距250Hz对应的相位角是/2。如果频率相距+250 Hz，相位角是+/2，然而导航电文引起的相位偏移将使相位角变成+/2(+/2-)，相应地引起-250 kHz的频率变换。如果在寻找精频时不考虑相位偏移，这个结果将变成50

41、0 kHz。 为了避免这个问题，频率不确定的最大值必须小于250Hz。 如果频率最大相差是200 Hz(这个值是实验选取的)，则相应的相位角的差是2/5，如图5-9所示。如果存在相位偏移，相位差了3/5(|(2/5)|)，这个值要比2/5大得多。 这样来看，相位差可以用来确定精频而不产生错误的频率偏移。 如果相位差大于2/5，就从结果中减去以保证频率差小于2/5。为了保持频率在200 kHz之内，X(k)中的k值之间的最大距离将是400 kHz。如果输入是在两个邻近的k值当中， 则输入信号跟两边的k值都相距200 kHz。图 5-9 相位差小于2/5将不会造成频率误差 图 5-10 从两个相位

42、角得到的角度差 讨论的最后一点是将X(k)的实部和虚部转化成相位角。 通常， 相位角的范围在之间。式(5-26)和式(5-27)中的两个相位角都将用这种方式得到，两个相位角的角度差可以是02间的任何值，如图5-10所示。由于200 kHz允许的最大相位角差是2/5，这两个相位角的角度差必须等于或小于2/5。如果结果大于2/5，可以从结果中加上或减去2，其绝对值就肯定小于2/5了。如果将噪声考虑进去，2/5门限值可以适当放宽，如用2.3/5，意味着其差必须等于或小于2.3/5这个值。如果这样处理后的值仍然大于这个门限值，就意味着两个连续毫秒数据间有相位偏移了，必须从结果中减去。当然，得到的角度值

43、也还要经过加减2，以保证最后得到的值小于门限值。 由以上讨论可知，在某一确定卫星中需要经过下面几个步骤来找出它的C/A码起始点和载波频率： (1) 对1 ms的输入执行循环相关操作，某个确定的C/A码起始点可以从这些循环相关里找到，载波频率可以以1 kHz的分辨率得到。 (2) 找到最高频率分量X(k)，在同一个毫秒数据内的两个分量，一个比X(k)中的k值低400 kHz，一个比k值高400 kHz， 对这两个分量执行DFT操作。三个输出X(k-1), X(k), X(k+1)中的最高输出将被指定为新的X(k)，并用这个新的X(k)作为DFT分量来求解精频。 (3) 从C/A码起始点处开始选择

44、连续几个毫秒的数据，我们随意选择为5ms，将这些数据与5组连续C/A码相乘，结果肯定是一个5 ms长的连续信号。但是在任意1 ms数据中，都有可能含有一个相位偏移。 (4) 从所有输入数据中找出Xn(k)，其中n=1, 2, 3, 4, 5， 然后由式(5-25)找出相位角。 角度差值定位为： (5-28) (5) 角度差的绝对值必须小于门限值(2.3/5)，如果不能实现，就要从上加上或减去2。如果结果仍然大于门限值，就要调整相位偏移，在上加上或减去，得到的值还要与门限值2.3/5进行比较，来决定是否要再次加减2。经过这些调整后，最终的角度值就是期望值。 (6) 式(5-27)用来计算精频。由

45、于有5ms的数据，将得到4个精频值。为了提高精确度，我们将这4个精频的平均值作为要求解的值。 5.5 Matlab与Simulink仿真简介 5.5.1 仿真工具简介 Matlab是由美国Mathsworks公司开发的大型软件。在Matlab软件中，包括了两大部分：数学计算和工程仿真。其数学计算部分提供了强大的矩阵处理和绘图功能。在工程仿真方面，Matlab提供的软件几乎遍布各个工程领域，并且不断加以完善，现已成为国际上最为流行的科学计算与工程计算软件之一。5.5.2 Simulink的使用本小节将利用BPSK传输作为例子对Simulink进行简单介绍。但是，我们只局限在Simulink的原理

46、及在这个仿真辅导中所要求的特点上，至于更为详尽的内容，建议读者参考Simulink手册。这里给出的Windows界面的图保持了原样。 Simulink是Matlab软件的扩展，它是实现动态建模和仿真的一个软件包。它与Matlab 语言的主要区别在于，它与用户接口是基于Windows的模型化图形输入的，从而使得用户可以把更多的精力投入到系统模型的构建而非语言的编程上。 图 5-11 Simulink Library Browser窗口 图 5-12 BPSK传输的Simulink模型 为了创建一个新的Simulink模型，可以通过FileNewModel菜单打开一个新的模型。某个库模块可以从库窗

47、口简单地拖放到Simulink模型窗口中。利用Uniform Random Number Generator模块以及随后的Sign模块，可以做图5-12中的二进制信号源。Uniform Random Number模块在指定的区间上产生均匀分布的随机数，在本例中该随机数在-1, +1内产生。跟在其后的Sign模块的输出是： 当输入大于0时输出为-1。为了产生输出比特为等概率的+1和-1二进制信号源，可以加入Lookup Table模块，它的输入输出映射关系为：1到1，-1到-1，而0到+1。为了将模块Uniform Random Number加到模块中，可以在Simulink Block Lib

48、rary Browser中双击模块Sources库(见图513)。将模块Uniform Random Number拖入Simulink模型并双击它。 图 5-13 Sources 库 打开的Block Parameters窗口如图5-14所示，在这个窗口中即可输入如该图所示的有关参数。现将Math Operations库中的Sign模块拖入模型。利用鼠标将模块Uniform Random Number和Sign连接起来：在模块Uniform Random Number的输出上按下鼠标左键并拖到Sign模块的输入端，然后释放鼠标键。Lookup Table模块从Lookup Tables库中得到

49、，在模型中双击Lookup Table模块并输入如图5-15所示的模块参数。 图 5-14 模块参数： Uniform Random Number 图 5-15 模块参数： Lookup Table 在开始仿真之前，为了观察源信号，必须调整Configuration Parameters。在菜单栏中使用命令SimulationConfiguration Parameters。即可打开如图5-16所示的Configuration Parameters窗口。 就目前而言， 仅涉及参数Start time和End time， 稍后为了求得准确的结果， 可能需要调整用于数值积分的Max step si

50、ze 和Min step size。在菜单栏中使用命令SimulationStart启动仿真。类似地， 可以使用SimulationStop命令结束或停止仿真。为了观察信号， 双击Scope即可打开如图5-17所示的窗口。 图 5-16 仿真参数窗口 图 5-17 源信号(Scope) 在窗口中单击鼠标右键会打开一个弹出菜单，通过该菜单可以打开Axes properties 窗口，如图5-18所示，这是为了设定y轴的范围。 图 5-18 Axes properties窗口（Scope） 在图5-17中，Scope窗口的时间范围可以通过菜单栏中的Parameters框来改变（见图5-19）。 图

51、 5-19 Scope的参数窗口 图5-12中的其他模块也能用相同的方式包括在内。利用UserDefined Functions 模块可以产生载波cos(4*pi*u)，该模块将某个指定的表达式（如cos(u)）加到输入变量u上。将模块Fcn拖入仿真模型并双击它，会立即打开Block Parameters窗口， 如图5-20所示， 在这里可以输入表达式cos(4*pi*u)。 图 5-20 模块参数： Fcn 时间变量u可以从Sources库中的Clock模块得到。Math Operations库中的Product模块用来产生已调信号。由于考虑的是无噪声传输仿真，因此就将发送端和接收端直接相连

52、。 接收端的匹配滤波器是利用Continuous库中的Integrator模块按积分清除来实现的。复位信号由Sources库中的Pulse Generator模块产生，其参数如图5-21所示。最后的仿真结果如图5-22所示。 图 5-21 模块参数： Pulse Generator 图 5-22 宿信号（Scope 8） 5.6 GPS卫星信号捕获的例子 程序中所用到的输入电文在时域内的情形如图5-23所示。 输入信号看起来像噪声，其频率描绘可通过FFT求得，如图5-24所示。正如预料的，FFT变换后带宽是2.5 MHz，其频谱形状与射频链的滤波器形状相似，通过循环相关后如图5-25所示， 其

53、中给出了6号卫星的C/A码起始点，其起始点在2884。图5-26给出了相距1 kHz的21个频率分量，最高值发生在k=7处。从图5-25和图5-26我们可以轻易地看出C/A码的起始点和频率。由于输入电文是实际收集到的数据，精频的精确度就难以测定， 因为其多普勒频率未知。精频还取决于在下变频中用到的本地振荡器的频率精确性和采样频率的精确性。 图5-23 5000个输入数据点 图5-24 输入信号的快速傅立叶变换 图 5-25 6号卫星的C/A码起始点 图5-26 卫星6的展开的各个频率分量 感知计算精频准确性的一个方法是利用数据的不同部分。 从输入电文的不同部分计算出六个精频，用到的电文部分是：

54、 125 000， 500130 001， 10 00135 001，15 00140 001，20 00145 001，25 00150 001。这些电文都是5 ms长，且起始点偏离1 ms，在相邻两组间，有4/5的数据相同。 因此，得到的五个精频是非常接近的，这六个精频之间相差是-2.4，9.0，-8.2，5.4和2.3 Hz。数据都是从固定的接收机上收集到的，因此，每毫秒刷新一次频率的变化是非常小的。 从而，这些频差可以看做是捕获方法中的误差。当信号长度改变时，精频的差也跟着变化，对一个弱信号来说，如果使用相同的捕获方法， 其精频可相差几十赫兹。 对弱信号（24号卫星）执行捕获，如图5-

55、27和图5-28所示。从图中我们很难判断C/A 码的起始点和频率是否为正确值，需要通过另外的相关运算才能找出它们，它的正确值与我们在图5-27和图5-28中看到的不同。 图5-27 24号卫星的C/A码起始点 图5-28 24号卫星的展开信号的频率分量 5.7 关于捕获的一些子程序 5.7.1 随机编码过程仿真打开Matlab，新建一个MFile，输入以下的文本，保存为GenerateCode.m，并调试运行。可以单步或多步运行M文件。%该程序仿真产生三组伪随机数， 模拟随机编码过程close all;clear;clc;t=0:5:1000; %时间变量为1 s， 步长5 ms%x1=ran

56、d(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(inde

57、x1)=-ones(1,length(index1); y1=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;for i=2: length(t) if (x1(i)=x1(i-1) y1(k)=

58、x1(i); t1(k)=t(i); k=k+1; else y1(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; end if (x2(i)=x2(i-1) y2(m)=x2(i); t2(m)=t(i); m=m+1; else y2(m)=x2(i-1); y2(m+1)=x2(i); t2(m)=t(i); t2(m+1)=t(i)+0.01; m=m+2; end if (x3(i)=x3(i-1) y3(n)=x3(i); t3(n)=t(i); n=n+

59、1; else y3(n)=x3(i-1); y3(n+1)=x3(i); t3(n)=t(i); t3(n+1)=t(i)+0.01; n=n+2; endend %画图说明产生的三组伪随机码%subplot(3,1,1);plot(t1,y1,y);grid on;axis(-1,1001,-1.5,1.5);subplot(3,1,2);plot(t2,y2,y);grid on;axis(-1,1001,-1.5,1.5);subplot(3,1,3);plot(t3,y3,y);grid on;axis(-1,1001,-1.5,1.5); whitebg(black);%将三组伪随

60、机值存放到新矩阵并传送到Workspace中， 方便以后调用数据%simin1=t1;y1;simin2=t2;y2;simin3=t3;y3; 输入上面的文本后，运行得到如图5-29所示的伪随机码。 其结果还可以通过Workspace查看，之所以要传送到Workspace中， 是因为我们在后面的仿真中还要用到这些伪随机码。 图 5.29 模拟三组伪随机码的仿真示意图 5.7.2 获取导航信息的仿真 在这个仿真中，我们要得到卫星的导航信息，其中忽略了频率偏移，且卫星信号的数据是基于前面章节中得到的卫星数据，存放在GPSsignal.mat文件中。在运行程序之前， 我们必须将GPSsignal.