北斗卫星导航信号串行捕获算法MATLAB仿真报告(附MATLAB程序)
北斗B1I信号的捕获算法
北斗B1I信号的捕获算法王丽黎;杨阳【摘要】北斗B1I信号的捕获是北斗2代接收机的核心模块,它是基于码相位和多普勒频移二维搜索的过程.对于捕获模块,通常采用并行码相位搜索捕获算法来实现对空中可见卫星的捕获.针对信号较弱情况下的卫星捕获,采用了非相干累加与并行码相位搜索捕获相结合的方法.测试结果表明,该捕获算法能够有效快速地实现弱信号的捕获.%BeiDou B1I signal acquisition is a hard core in a BeiDou receiver and is the two-dimension process of the search about the code delay and Doppler frequency. For acquisition, parallel search algorithm of code in frequency domain are usually adopted. For the acquisition of weak signal, a method combining non-correlation integration with parallel search algorithm of code in frequency domain was introduced. The test results show that the acquisition algorithm can effectively and quickly acquire the weak signal.【期刊名称】《计算机系统应用》【年(卷),期】2016(025)003【总页数】5页(P194-198)【关键词】北斗B1I信号;接收机;捕获;弱信号;非相干累加【作者】王丽黎;杨阳【作者单位】西安理工大学自动化与信息工程学院,西安 710048;西安理工大学自动化与信息工程学院,西安 710048【正文语种】中文随着我国成功将第16颗北斗导航卫星送入预定轨道, 北斗导航工程区域组网顺利完成[1]. 第2代北斗导航系统主要功能为定位、测速、单双向授时和短报文通信[2,3]. 因此, 针对北斗导航系统相应的接收技术的研究逐渐成为研究热点. 传统的接收机由射频前端、用于信号处理的ASIC以及高速运算的CPU核组成, 设计灵活性受到限制. 相比之下, 软件接收机只需对软件修改便可对接收机进行优化升级, 则更具有方便性和灵活性[4,5].在接收机内, 完成信号捕获是信号处理的第一步, 信号捕获重在估计两个重要参数: 一个是C/A码周期的开始, 另一个是输入信号的载波频率, 获得导航卫星信号的载波频率(或称为载波多普勒频移)和码相位这两个参数的粗略值, 然后跟踪过程则利用这些粗略值进一步得到频率和码相位的精确值并进而解算出导航电文. 因此对接收机性能好坏影响较大的是捕获跟踪算法的好坏[6,7].本文在详细分析了并行码相位搜索捕获算法原理的基础上, 为了能更好的实现对弱信号的捕获, 将并行码相位搜索捕获算法与非相干累加相结合, 不但提高捕获效率, 降低了噪声干扰, 而且有效的提高了弱信号的捕获.目前北斗2代播发B1、B2、B3三个频段的信号, 其信号复用方式为码分多址, 其中B1(1561.098 MHz)信号由I、Q两条支路的测距码和导航电文正交调制在载波上构成. 调制在B1频率上的信号可表示为式(1):其中上角标j表示卫星序号; AB1分别表示B1信号幅度; 分别表示B1信号测距码; 分别表示调制在B1测距码上的导航电文数据码; 表示B1信号载波频率; 表示B1信号载波初相.CB1I码的码速率为2.046 Mcps, 码长为2046, 码宽为488.7 ns(1/2.046 MHz). CB1I码发生的结构如图1所示, 其是由两个线性序列G1和G2模2和产生均衡Gold码后截短1个码片后生成. G1和G2序列分别由两个11级的线性移位寄存器生成, 其生成多项式如式(2)、式(3)所示.(2)(3)G1序列初始相位为************;G2序列初始相位为************.通过对产生G2序列的移位寄存器不同抽头的模2和可以实现G2序列相位的不同偏移, 与G1序列模2和后可生成不同卫星的测距码. 而之所以用这种码传输信号, 就是因为其良好的自相关和互相关特性.图1 CB1I码发生器示意图2 北斗2代B1信号CB1I码的捕获由于CB1I码除了自身完全对齐的情况外, 其余情况几乎是不相关的, 这种特性使得很容易找出两个完全对齐的相同的CB1I码. 捕获过程正是利用其这一特点.2.1 CB1I码的捕获接收机的信号捕获过程一般通过对卫星信号的载波频率和码相位进行扫描式搜索来完成. 捕获的目的就是为了对输入信号和一个测距码序列做相关运算. 并行码相位搜索捕获算法实际上是利用傅里叶变换这种数字信号处理技术来替代数字相关器的相关运算, 而我们需要证明一下两者的等价性.两个长度同为N的有限长序列和的离散傅里叶变换和计算如式(4)、式(5)所示.(4)(5)两个长度同为N的有限长序列和的循环互相关计算如式(6)所示.(6)下面的分析中均省略了中的缩放因子1/N, 的N点离散傅里叶变换计算如式(7)所示.(7)其中是的复共轭. 由式(7)可知, 两个序列与在时域内做相关运算, 相当于它们的离散傅里叶变换与在频域内做乘积运算. 于是倒过来, 乘积的离散傅里叶反变换正好是接收机需要进行检测的在各个码相位处的相关值. 一旦接收机通过傅里叶反变换计算得到相关值, 那么接下来的信号检测就同线性搜索捕获法一样, 即找出在所有搜索单元中自相关幅值的峰值, 并将该峰值与捕获门限值相比较. 若峰值超过捕获门限值, 则接收机捕获到了信号.图2 并行码相位搜索捕获算法原理图如图2所示为并行码相位搜索捕获算法的原理图, 考虑到导航数据位存在跳变的可能, 用含导航数据位的数据进行相关运算所获得的相关峰值将会有较大衰减从而造成漏捕, 在算法实现过程中, 总是采用两段连续数据进行同步相关运算, 在得到的两组相关结果中选择最大功率能量相关峰值较大者(认为该段数据中不包含导航数据位跳变)的相关结果作为捕获判断依据[8].2.2基于非相干累加的CB1I码捕获采用1 ms时长的数据进行上述捕获的时候, 由于噪声的作用可能导致误捕, 当信号较弱时, 甚至会出现漏捕. 而且在北斗MEO/IGSO卫星B1频点信号中, 因为NH 码调制的影响, 相干累加时间不能超过1 ms. 为了提高灵敏度, 只能通过提高非相干累加次数来捕获弱信号[9], 即将相干积分结果进行平方处理后再累加, 从而获得信号增益. 其非相关累加值[10]可表示为式(8):(8)其中: 为M ms的数据信息与测距码的非相干捕获相关值. 可以看出导航电文数据位翻转对积分结果的影响被平方运算有效的降低了. 因此非相干积分可以进行积分时间超过1 ms的积分.非相干累加法消除了导航电文数据位翻转造成的影响, 同时平方运算消除了相位误差造成的副作用, 而前面提到的并行码相位搜索捕获算法通过傅里叶变换实现循环相关, 将相位域捕获过程并行化, 使得搜索量减少到了只需搜索不同的载波频率, 提高运算效率. 将非相干累加法与并行码相位搜索捕获算法相结合, 在大幅降低捕获时间的前提下, 实现了对弱信号的捕获. 其原理图如图3所示.图3 非相干累加捕获原理图通过图3可知, 非相干累加捕获的方法是将并行码相位搜索捕获算法每毫秒的捕获结果按照预先设定好的累加时间进行累加, 其累加原理如式(8)所示, 由于噪声累加的结果增大不如信号累加的结果增加的快, 使得经过一段时间累加后, 可找出明显的相关峰值.3 仿真验证在Matlab环境下对本文研究的信号捕获方案进行仿真验证. 利用卫星信号模拟器对北斗B1频点中频信号进行仿真, 仿真信号的中频频率为2.098MHz, 采样频率为8.8MHz, 仿真产生60s的数据中频信号, 信号中共调制了1号、2号、3号、4号、7号、8号、10号、12号和13号等9颗卫星信号. 信号中加入的是高斯白噪声,信噪比为-35dB. 利用生成的信号, 就可以进行捕获的仿真, 捕获程序流程图如图4所示.图4 CB1I码捕获程序流程图图5是用图1所示的CB1I码发生器产生的对应10号卫星的本地伪码, 其是由两个11级移位寄存器进行模2和生成的. 不同的卫星编号对应不同的抽头, 不同卫星对应的CB1I码则通过查表的方式就可以实现. 横坐标表示采样点数, 截取了2046个码片的前100位, 纵坐标表示CB1I码的相位幅度.图5 10号卫星本地伪码部分截图CB1I码具有良好的自相关和互相关特性, 如图6和图7所示. 除了延迟为零外, 几乎没有自相关性. 只有当本地伪码与接收到的信号的伪码序列能够对齐时才可得到最大相关值. 根据这一特性可轻易找出何时两个码是严格对齐的, 本文采用并行码相位搜索捕获也正是基于此特性. 横坐标表示码片数, 纵坐标分别表示自相关值和互相关值.图6 测距码的自相关性图7 测距码的互相关性当检测门限选用最大峰值与次大峰值的比值(大于2.5)时, 首先使用两段连续的单位数据段, 本文以1 ms数据位为单位数据段, 对其做同步相关运算即对这两段数据进行并行码相位搜索捕获, 在得到的两组相关结果中选择较大的相关峰值作为捕获判断依据进行弱信号捕获的验证结果如图8所示. 其次使用5 ms和10 ms数据进行非相干积分的弱信号捕获验证结果如图9和图10所示.图8 并行码相位搜索捕获算法捕获情况图9 进行5 ms非相干累加的捕获情况图10 进行10 ms非相干累加的捕获情况由图可以看出, 在信噪比为-35dB的情况下, 图8中仅使用并行码相位搜索捕获算法得到的最大峰值与次大峰值的比值并不是很大; 而在图9和图10中采用将并行码相位搜索捕获算法与非相干累加相结合的捕获算法, 在累加时间增大时, 该比值结果明显增加. 证明当进行非相干积分所用数据长度从1ms增加到10ms时, 信号中所有可见卫星的最大峰值与次大峰值的比值增加的都很明显, 能够实现低信噪比信号的捕获, 提高接收机的灵敏度. 当然, 也可以根据实际的需要选择合适的相干累加时间, 达到设计目的.接着再从单颗卫星角度进行研究, 以10号卫星为例, 即PRN=10, 分别进行3ms、6ms和10ms数据的相干累加, 其捕获结果如图11、图12和图13所示.图11 10号卫星进行3 ms非相干累加结果图12 10号卫星进行6 ms非相干累加结果图13 10号卫星进行10 ms非相干累加结果由图11、图12和图13可以看出, 随着捕获所用数据长度的增加, 即非相干累加数据长度增加, 噪声得到一定的抑制, 可见卫星正确相位所对应的归一化相关值也更加明显, 该部分也达到了有效地捕获弱信号的预期效果.4 结语本文研究了北斗软件接收机捕获算法, 分析了并行码相位搜索捕获算法并将其与非相干积分相结合, 并通过仿真数据对算法进行了验证. 可见非相干积分与并行码相位搜索捕获算法相结合不但捕获效率高, 而且较好地抑制了噪声信号, 有效地实现了弱信号的捕获. 对于软件接收机相关模块的研究具有一定的意义, 能够使用户在接收机算法处理和软件更新等方面具有很大的灵活性.参考文献1 何敏,葛榜军.北斗卫星导航系统及应用.卫星应用,2012, (5):19–23.2 中国卫星导航系统管理办公室.北斗导航系统发展报告.国际太空,2012,(4):6–11.3 Meng WX, Liu E, Han Sh. Resaerch and development on satellite positioning and navigation in China. IEICE Trans. Commun, 2012(E95-B): 3385–3392.4 杨东凯,张飞舟,张波译.软件定义的GPS和伽利略接收机.北京:国防工业出版社,2009.5 杨俊,武奇生.GPS基本原理及其Matlab仿真.西安:西安电子科技大学出版社,2006.6 王冰.GPS信号捕获算法的研究.电子科技,2014,27(8): 154–156.7 谢刚.GPS原理与接收机设计.北京:电子工业出版社,2009.8 黄隽祎,李荣冰,王翌等.北斗B1 QPSK调制信号的高灵敏度捕获算法.航空计算技术,2012,42(5):38–42.9 史向男,巴晓辉,陈杰.北斗MEO/IGSO卫星B1频点信号捕获方法研究.国外电子测量技术,2013,32(4):19–21.10 陈军,潘高峰,李飞,余金峰,黄静华译. GPS软件接收机基础(第2版).北京:电子工业出版社,2007.Acquisition Algorithm of BeiDou B1I SignalWANG Li-Li, YANG Yang(Faculty of Automation and Information Engineering, Xi′an University of Technology, Xi’an 710048, China)Abstract：BeiDou B1I signal acquisition is a hard core in a BeiDou receiver and is the two-dimension process of the search about the code delay and Doppler frequency. For acquisition, parallel search algorithm of code in frequency domain are usually adopted. For the acquisition of weak signal, a method combining non-correlation integration with parallel search algorithm of code in frequency domain was introduced. The test results show that the acquisition algorithm can effectively and quickly acquire the weak signal.Key words：BeiDou B1I signal; receiver; acquisition; weak signal; non-coherent integration① 收稿时间：2015-07-06;收到修改稿时间:2015-09-06CB1I码的码速率为2.046 Mcps, 码长为2046, 码宽为488.7 ns(1/2.046 MHz). CB1I码发生的结构如图1所示, 其是由两个线性序列G1和G2模2和产生均衡Gold码后截短1个码片后生成. G1和G2序列分别由两个11级的线性移位寄存器生成, 其生成多项式如式(2)、式(3)所示.G1序列初始相位为************;G2序列初始相位为************.通过对产生G2序列的移位寄存器不同抽头的模2和可以实现G2序列相位的不同偏移, 与G1序列模2和后可生成不同卫星的测距码. 而之所以用这种码传输信号, 就是因为其良好的自相关和互相关特性.由于CB1I码除了自身完全对齐的情况外, 其余情况几乎是不相关的, 这种特性使得很容易找出两个完全对齐的相同的CB1I码. 捕获过程正是利用其这一特点.2.1 CB1I码的捕获接收机的信号捕获过程一般通过对卫星信号的载波频率和码相位进行扫描式搜索来完成. 捕获的目的就是为了对输入信号和一个测距码序列做相关运算. 并行码相位搜索捕获算法实际上是利用傅里叶变换这种数字信号处理技术来替代数字相关器的相关运算, 而我们需要证明一下两者的等价性.两个长度同为N的有限长序列和的离散傅里叶变换和计算如式(4)、式(5)所示. 两个长度同为N的有限长序列和的循环互相关计算如式(6)所示.下面的分析中均省略了中的缩放因子1/N, 的N点离散傅里叶变换计算如式(7)所示.其中是的复共轭. 由式(7)可知, 两个序列与在时域内做相关运算, 相当于它们的离散傅里叶变换与在频域内做乘积运算. 于是倒过来, 乘积的离散傅里叶反变换正好是接收机需要进行检测的在各个码相位处的相关值. 一旦接收机通过傅里叶反变换计算得到相关值, 那么接下来的信号检测就同线性搜索捕获法一样, 即找出在所有搜索单元中自相关幅值的峰值, 并将该峰值与捕获门限值相比较. 若峰值超过捕获门限值, 则接收机捕获到了信号.如图2所示为并行码相位搜索捕获算法的原理图, 考虑到导航数据位存在跳变的可能, 用含导航数据位的数据进行相关运算所获得的相关峰值将会有较大衰减从而造成漏捕, 在算法实现过程中, 总是采用两段连续数据进行同步相关运算, 在得到的两组相关结果中选择最大功率能量相关峰值较大者(认为该段数据中不包含导航数据位跳变)的相关结果作为捕获判断依据[8].2.2基于非相干累加的CB1I码捕获采用1 ms时长的数据进行上述捕获的时候, 由于噪声的作用可能导致误捕, 当信号较弱时, 甚至会出现漏捕. 而且在北斗MEO/IGSO卫星B1频点信号中, 因为NH码调制的影响, 相干累加时间不能超过1 ms. 为了提高灵敏度, 只能通过提高非相干累加次数来捕获弱信号[9], 即将相干积分结果进行平方处理后再累加, 从而获得信号增益. 其非相关累加值[10]可表示为式(8):其中: 为M ms的数据信息与测距码的非相干捕获相关值. 可以看出导航电文数据位翻转对积分结果的影响被平方运算有效的降低了. 因此非相干积分可以进行积分时间超过1 ms的积分.非相干累加法消除了导航电文数据位翻转造成的影响, 同时平方运算消除了相位误差造成的副作用, 而前面提到的并行码相位搜索捕获算法通过傅里叶变换实现循环相关, 将相位域捕获过程并行化, 使得搜索量减少到了只需搜索不同的载波频率, 提高运算效率. 将非相干累加法与并行码相位搜索捕获算法相结合, 在大幅降低捕获时间的前提下, 实现了对弱信号的捕获. 其原理图如图3所示.通过图3可知, 非相干累加捕获的方法是将并行码相位搜索捕获算法每毫秒的捕获结果按照预先设定好的累加时间进行累加, 其累加原理如式(8)所示, 由于噪声累加的结果增大不如信号累加的结果增加的快, 使得经过一段时间累加后, 可找出明显的相关峰值.在Matlab环境下对本文研究的信号捕获方案进行仿真验证. 利用卫星信号模拟器对北斗B1频点中频信号进行仿真, 仿真信号的中频频率为2.098MHz, 采样频率为8.8MHz, 仿真产生60s的数据中频信号, 信号中共调制了1号、2号、3号、4号、7号、8号、10号、12号和13号等9颗卫星信号. 信号中加入的是高斯白噪声,信噪比为-35dB. 利用生成的信号, 就可以进行捕获的仿真, 捕获程序流程图如图4所示.图5是用图1所示的CB1I码发生器产生的对应10号卫星的本地伪码, 其是由两个11级移位寄存器进行模2和生成的. 不同的卫星编号对应不同的抽头, 不同卫星对应的CB1I码则通过查表的方式就可以实现. 横坐标表示采样点数, 截取了2046个码片的前100位, 纵坐标表示CB1I码的相位幅度.CB1I码具有良好的自相关和互相关特性, 如图6和图7所示. 除了延迟为零外, 几乎没有自相关性. 只有当本地伪码与接收到的信号的伪码序列能够对齐时才可得到最大相关值. 根据这一特性可轻易找出何时两个码是严格对齐的, 本文采用并行码相位搜索捕获也正是基于此特性. 横坐标表示码片数, 纵坐标分别表示自相关值和互相关值.当检测门限选用最大峰值与次大峰值的比值(大于2.5)时, 首先使用两段连续的单位数据段, 本文以1 ms数据位为单位数据段, 对其做同步相关运算即对这两段数据进行并行码相位搜索捕获, 在得到的两组相关结果中选择较大的相关峰值作为捕获判断依据进行弱信号捕获的验证结果如图8所示. 其次使用5 ms和10 ms数据进行非相干积分的弱信号捕获验证结果如图9和图10所示.由图可以看出, 在信噪比为-35dB的情况下, 图8中仅使用并行码相位搜索捕获算法得到的最大峰值与次大峰值的比值并不是很大; 而在图9和图10中采用将并行码相位搜索捕获算法与非相干累加相结合的捕获算法, 在累加时间增大时, 该比值结果明显增加. 证明当进行非相干积分所用数据长度从1ms增加到10ms时, 信号中所有可见卫星的最大峰值与次大峰值的比值增加的都很明显, 能够实现低信噪比信号的捕获, 提高接收机的灵敏度. 当然, 也可以根据实际的需要选择合适的相干累加时间, 达到设计目的.接着再从单颗卫星角度进行研究, 以10号卫星为例, 即PRN=10, 分别进行3ms、6ms和10ms数据的相干累加, 其捕获结果如图11、图12和图13所示.由图11、图12和图13可以看出, 随着捕获所用数据长度的增加, 即非相干累加数据长度增加, 噪声得到一定的抑制, 可见卫星正确相位所对应的归一化相关值也更加明显, 该部分也达到了有效地捕获弱信号的预期效果.本文研究了北斗软件接收机捕获算法, 分析了并行码相位搜索捕获算法并将其与非相干积分相结合, 并通过仿真数据对算法进行了验证. 可见非相干积分与并行码相位搜索捕获算法相结合不但捕获效率高, 而且较好地抑制了噪声信号, 有效地实现了弱信号的捕获. 对于软件接收机相关模块的研究具有一定的意义, 能够使用户在接收机算法处理和软件更新等方面具有很大的灵活性.1 何敏,葛榜军.北斗卫星导航系统及应用.卫星应用,2012, (5):19–23.2 中国卫星导航系统管理办公室.北斗导航系统发展报告.国际太空,2012,(4):6–11.3 Meng WX, Liu E, Han Sh. Resaerch and development on satellite positioning and navigation in China. IEICE Trans. Commun, 2012(E95-B): 3385–3392.4 杨东凯,张飞舟,张波译.软件定义的GPS和伽利略接收机.北京:国防工业出版社,2009.5 杨俊,武奇生.GPS基本原理及其Matlab仿真.西安:西安电子科技大学出版社,2006.6 王冰.GPS信号捕获算法的研究.电子科技,2014,27(8): 154–156.7 谢刚.GPS原理与接收机设计.北京:电子工业出版社,2009.8 黄隽祎,李荣冰,王翌等.北斗B1 QPSK调制信号的高灵敏度捕获算法.航空计算技术,2012,42(5):38–42.9 史向男,巴晓辉,陈杰.北斗MEO/IGSO卫星B1频点信号捕获方法研究.国外电子测量技术,2013,32(4):19–21.10 陈军,潘高峰,李飞,余金峰,黄静华译. GPS软件接收机基础(第2版).北京:电子工业出版社,2007.。
基于某MATLAB地GPS信号仿真完整源代码123
配套毕业设计论文见百度文库请搜索《基于MATLAB的GPS信号仿真123》附录C 仿真程序代码1、数据码的产生function datacode=data(x)y=rand(1,x);for i=1:xif y(i)<0.5datacode(i)=0;elsedatacode(i)=1;endendy(1)=0;show2(1)=datacode(1);q=2;for i=1:length(datacode)for j=1:100y(q)=i-1+j*0.01;show2(q)=datacode(i);q=q+1;endendplot(y,show2);axis([0 length(datacode) -0.2 1.2]);1、C/A码的产生及扩频调制clc;c=input('请输入数据码的长度：c=');y=rand(1,c);for i=1:cif y(i)<0.5datacode(i)=0;elsedatacode(i)=1;endendx(1)=0;show(1)=datacode(1);p=2;for i=1:cfor j=1:100x(p)=i-1+j*0.01;show(p)=datacode(i);p=p+1;endendsubplot(4,1,1);plot(x,show);title('数据码');axis([0 c -0.2 1.2]);number=input('请输入卫星PRN号码:number=');cacode=CAgenerate(number);temp=cacode(1:100)x(1)=0;show(1)=temp(1);p=2;%下面的循环是为了将结果显示成方波形式for i=1:length(temp)for j=1:100x(p)=i-1+j*0.01;show(p)=temp(i);p=p+1;endend%画出仿真结果图subplot(4,1,2);plot(x,show);title('C/A码');axis([0 100 -0.2 1.2]);%截取CA码的前十个数据进行扩频，每个数据插入5个CA序列 cacode1=cacode(1:10);for i=1:cif datacode(i)==1datacodek((i-1)*50+1:i*50)=ones(1,50);elsedatacodek((i-1)*50+1:i*50)=zeros(1,50);endendfor i=1:cfor j=1:50addr=rem(((i-1)*50+j),10);if addr==0addr=10;endkuopindata((i-1)*50+j)=xor(datacodek((i-1)*50+j),cacode1(addr));endend%下面的循环是为了将结果显示成方波形式x(1)=0;show(1)=kuopindata(1);p=2;for i=1:length(kuopindata)for j=1:100x(p)=i-1+j*0.01;show(p)=kuopindata(i);p=p+1;endendsubplot(4,1,3);plot(x,show);title('扩频数据');axis([0 length(kuopindata) -0.2 1.2]);%每位数据通过正弦波来调制Sinwave=sin([0:2*pi/8:2*pi*7/8]);Sinwave=single(Sinwave);GPSsignal=zeros(1,1);Sinwave=[Sinwave Sinwave Sinwave Sinwave Sinwave];for i=1:length(kuopindata)GPSsignal=[GPSsignal kuopindata(i)*Sinwave];endGPSsignal=GPSsignal(2:length(GPSsignal));subplot(4,1,4);plot(GPSsignal(1:500));title('调制后数据');C/A码产生的子程序CAgenerate：function cacode=CAgenerate(number)if (number<1)|(number>37)disp('输入参数必须在1 ～ 37之间取值');returnendCACode=zeros(1,1023); %生成一个1*1023的零矩阵% 设置寄存器初相Reg1=[1,1,1,1,1,1,1,1,1,1];Reg2=[1,1,1,1,1,1,1,1,1,1];% 设置反馈点,1表示需要反馈gp1=[0,0,1,0,0,0,0,0,0,1];gp2=[0,1,1,0,0,1,0,1,1,1];% 抽头G2Table=[ 2,3,4,5,1,2,1,2,3,2,3,5,6,7,8,9,1,2,3,4,5,6,1,4,5,6,7,8,1,2,3 ,4,5,4,1,2,4;6,7,8,9,9,10,3,4,6,7,8,9,10,4,5,6,7,8,9,3,6,7,8,9,10,6,7,8,9,10,10,7,8,10;] % 生成一个周期的伪码序列for m=1:1023CACode(m)=mod(Reg1(10)+Reg2(G2Table(1,number))+Reg2(G2Table(2,number)),2);Reg1=[mod(Reg1*gp1',2),Reg1(1:9)];Reg2=[mod(Reg2*gp2',2),Reg2(1:9)];endcacode=CACode;2、C/A码的相关性分析clc;n=input('请输入卫星PRN号码:n=');cacode1=CAgenerate(n);%在G2序列中找出-1并转换为0,找出1并转换为1ind1=find(cacode1==1);ind2=find(cacode1==0);cacode1(ind1)=-ones(1,length(ind1)); cacode1(ind2)=ones(1,length(ind2)); N=1023;z=zeros(1,1023);for i=0:N-1for k=i+1:N-1z1(k)=cacode1(k)*cacode1(k-i); z(i+1)=z(i+1)+z1(k);endz(i+1)=z(i+1)/N;endsubplot(2,1,1);plot(z);title('自相关特性');axis([-50 1300 -0.5 1.2]);n=input('请输入卫星PRN号码:n='); cacode2=CAgenerate(n);ind1=find(cacode2==1);ind2=find(cacode2==0);cacode2(ind1)=-ones(1,length(ind1)); cacode2(ind2)=ones(1,length(ind2));N=1023;h=zeros(1,1023);for i=0:N-1for k=i+1:N-1h1(k)=cacode1(k)*cacode2(k-i); h(i+1)=h(i+1)+h1(k);endh(i+1)=h(i+1)/N;endsubplot(2,1,2);plot(h);title('互相关特性');axis([-50 1300 -0.5 1]);4、 P码的产生及扩频调制clc;c=input('请输入数据码的长度：c='); y=rand(1,c);for i=1:cif y(i)<0.5datacode(i)=0;elsedatacode(i)=1;endendx(1)=0;show(1)=datacode(1);p=2;for i=1:cfor j=1:100x(p)=i-1+j*0.01;show(p)=datacode(i);p=p+1;endendsubplot(4,1,1);plot(x,show);title('数据码');axis([0 c -0.2 1.2]);NumberPCode=input('enter the NumberPcode='); NumberShift=input('enter the NumberShift=');a=input('enter a=');pcode=Pcode(a,NumberPCode,NumberShift);x(1)=0;show(1)=pcode(1);p=2;for i=1:length(pcode)for j=1:100x(p)=i-1+j*0.01;show(p)=pcode(i);p=p+1;endendsubplot(4,1,2);plot(x,show);title('P码');axis([0 length(pcode) -0.2 1.2]);pcode=pcode(1:10);for i=1:cif datacode(i)==1datacodek((i-1)*50+1:i*50)=ones(1,50); elsedatacodek((i-1)*50+1:i*50)=zeros(1,50); endendfor i=1:cfor j=1:50addr=rem(((i-1)*50+j),10);if addr==0addr=10;endkuopindata((i-1)*50+j)=xor(datacodek((i-1)*50+j),pcode(addr)); endendx(1)=0;show(1)=kuopindata(1);p=2;%下面的循环是为了将结果显示成方波形式for i=1:length(kuopindata)for j=1:100x(p)=i-1+j*0.01;show(p)=kuopindata(i);p=p+1;endendsubplot(4,1,3);plot(x,show);title('扩频数据');axis([0 length(kuopindata) -0.2 1.2]);%每位数据通过正弦波来调制Sinwave=sin([0:2*pi/8:2*pi*7/8]);Sinwave=single(Sinwave);GPSsignal=zeros(1,1);Sinwave=[Sinwave Sinwave Sinwave Sinwave Sinwave];for i=1:length(kuopindata)GPSsignal=[GPSsignal kuopindata(i)*Sinwave];endGPSsignal=GPSsignal(2:length(GPSsignal));subplot(4,1,4);title('调制后数据');plot(GPSsignal(1:500));以下是P码产生的子程序Pcode:function pcode=Pcode(a,NumberPCode,NumberShift) % P码产生reg1a=[0 0 0 1 0 0 1 0 0 1 0 0];reg1b=[0 0 1 0 1 0 1 0 1 0 1 0];reg2a=[1 0 1 0 0 1 0 0 1 0 0 1];reg2b=[0 0 1 0 1 0 1 0 1 0 1 0];rx1a=0;rx1b=0;rx2a=0;rx2b=0;x1bWork=1;x2aWork=1;x2bWork=1;N=NumberShift;C1=4092*3750;C2=4093*3749;z1a=mod(N,4092);%取余数x1a=mod([(N-z1a)/4092],3750);y1a=(N-z1a-4092*x1a)/C1;if ((N-C1*y1a)>=C2)z1b=4092;x1bWork=0;x1b=3748;elsez1b=mod((N-C1*y1a),4093);x1bWork=1;x1b=(N-z1b-C1*y1a)/4093;endm=mod(N,(C1+37));y2a=(N-m)/(C1+37);if (m>=C1)dv=m-C1;elsedv=0;endz2a=mod((m-dv),4092);x2a=mod((((m-dv)-z2a)/4092),3750);z2b=mod((m-dv),4093);if (m>=C2)x2b=3748;elsex2b=(m-z2b)/4093;end%各移位寄存器的状态for i=1:z1aslave1a=mod(reg1a(6)+reg1a(8)+reg1a(11)+reg1a(12),2);reg1a(2:12)= reg1a(1:11);reg1a(1)=slave1a;endfor i=1:z1bslave1b=mod(reg1b(1)+reg1b(2)+reg1b(5)+reg1b(8)+reg1b(9)+reg1b(10)+reg1b(11 )+reg1b(12),2);reg1b(2:12)=reg1b(1:11);reg1b(1)=slave1b;endfor i=1:z2aslave2a=mod(reg2a(1)+reg2a(3)+reg2a(4)+reg2a(5)+reg2a(7)+reg2a(8)+reg2a(9)+ reg2a(10)+reg2a(11+reg2a(12)) ,2);reg2a(2:12)=reg2a(1:11);reg2a(1)=slave2a;endfor i=1:z2bslave2b=mod(reg2b(2)+reg2b(3)+reg2b(4)+reg2b(8)+reg2b(9)+reg2b(12) ,2);reg2b(2:12)=reg2b(1:11);reg2b(1)=slave2b;end%各控制变量的判断if z1a==4091rx1a=1;endif z1b==4092rx1b==1;endif z2a==4091rx2a=1;x2aWork=0;endif z2b==4092rx2b=1;x2bWork=0;end%开始产生P码p=zeros(NumberPCode,1);x1acou=0;x1bcou=0;x2acou=0;x2bcou=0;cou37=dv;x2(1:a)=1;for i=1:(NumberPCode+37)x1(i)=mod( reg1a(12)+reg1b(12),2);x2(i+a)=mod( reg2a(12)+reg2b(12),2);%寄存器x1b的移位函数if x1bWork==1if rx1b==0slave1b=mod(reg1b(1)+reg1b(2)+reg1b(5)+reg1b(8)+reg1b(9)+reg1b(10)+reg1b(11)+reg1b(12),2);reg1b(2:12)=reg1b(1:11);reg1b(1)=slave1b;else if rx1b==1reg1b=[0 0 1 0 1 0 1 0 1 0 1 0];rx1b=0;endendelse if x1bWork==0endendif reg1b==[0 1 0 1 0 1 0 1 0 1 0 0 ]rx1b=1;x1bcou=x1bcou+1;if x1bcou==3749x1bwork=0;x1bcou=0;endend%寄存器x1a的移位函数if rx1a==0slave1a=mod(reg1a(6)+reg1a(8)+reg1a(11)+reg1a(12),2);reg1a(2:12)=reg1a(1:11);reg1a(1)=slave1a;else if rx1a==1reg1a=[0 0 0 1 0 0 1 0 0 1 0 0];rx1a=0;endendif reg1a==[0 0 1 0 0 1 0 0 1 0 0 0]rx1a=1;x1acou=x1acou+1;if x1acou==3750x1bwork=1;x1acou=0;endend%寄存器x2b的移位函数if x2bWork==1if rx2b==0slave2b=mod(reg2b(2)+reg2b(3)+reg2b(4)+reg2b(8)+reg2b(9)+reg2b(12) ,2); reg2b(2:12)=reg2b(1:11);reg2b(1)=slave2b;else if rx2b==1x2bout=reg2b(11);reg2b=[0 0 1 0 1 0 1 0 1 0 1 0];rx2b=0;endendelse if x2bWork==0reg2b=[0 0 1 0 1 0 1 0 1 0 1 0];x2bWork=1;rx2b=0;endendif reg2b==[0 1 0 1 0 1 0 1 0 1 0 0]rx2b=1;x2bcou=x2bcou+1;if x2bcou==3749x2bWork=0;x2bcou=0;endend%寄存器x2a的移位函数if x2aWork==1if rx2a==0slave2a=mod(reg2a(1)+reg2a(3)+reg2a(4)+reg2a(5)+reg2a(7)+reg2a(8)+reg2a(9)+reg2a(10)+reg2a(11)+reg2a(12) ,2); reg2a(2:12)=reg2a(1:11);reg2a(1)=slave2a;else if rx2a==1reg2a=[1 0 1 0 0 1 0 0 1 0 0 1];rx2a=0;endendelse if x2aWork==0if rx2a==1cou37=cou37+1;if cou37==37rx2a=0;x2awork=1;cou37=0;endendendendif reg2a==[0 1 0 0 1 0 0 1 0 0 1 1] rx2a=1;x2acou=x2acou+1;if x2acou==3750x2awork=0;x2acou=0;endendendfor i=1:NumberPCodep(i)= mod( x1(i)+x2(i),2);endp=p';pcode=p';5、 P码的相关性分析clc;NumberPCode=input('enter the number Pcode='); NumberShift=input('enter the numbershift='); a=input('enter a=');pcode=Pcode(a,NumberPCode,NumberShift);ind1=find(pcode==1);ind2=find(pcode==0);pcode(ind1)=-ones(1,length(ind1));pcode(ind2)=ones(1,length(ind2));M=NumberPCode;z=zeros(1,M);for i=0:M-1for k=i+1:M-1z1(k)=pcode(k)*pcode(k-i);z(i+1)=z(i+1)+z1(k);endz(i+1)=z(i+1)/M;endsubplot(2,1,1);plot(z);title('自相关特性');axis([-50 M -0.5 1.2]);a=input('enter a=');NumberShift=input('enter the numbershift='); pcode2=Pcode(a,NumberPCode,NumberShift);h=zeros(1,M);for i=0:M-1for k=i+1:M-1h1(k)=pcode(k)*pcode2(k-i);h(i+1)=h(i+1)+h1(k);endh(i+1)=h(i+1)/M; endsubplot(2,1,2);plot(h);title('互相关特性'); axis([-50 M -0.5 1]);。
北斗卫星导航信号串行捕获算法MATLAB仿真报告(附MATLAB程序)
北斗卫星导航信号串行捕获算法MATLAB仿真报告一、原理卫星导航信号的串行捕获算法如图1所示。
图1 卫星导航信号的串行捕获算法接收机始终在本地不停地产生对应某特定卫星的本地伪码，并且接收机知道产生的伪码的相位，这个伪码按一定速率抽样后与接收的GPS中频信号相乘，然后再与同样知晓频率的本地产生的载波相乘。
GPS中频信号由接收机的射频前端将接收到的高频信号下边频得到。
实际产生对应相位相互正交的两个本地载波，分别称为同相载波和正交载波，信号与本地载波相乘后的信号分别成为，产生同相I支路信号和正交的Q 支路信号。
两支路信号分别经过一个码周期时间的积分后，平方相加。
分成两路是因为C/A码调制和P码支路正交的支路上，假设是I支路。
当然由于信号传输过程中引入了相位差，解调时的I支路不一定是调制时的I支路，Q支路也一样，二者不一定一一对应，因此为了确定是否检测到接收信号，需要同时对两支路信号进行研究。
相关后的积分是为了获取所有相关数据长度的值的相加结果，平方则是为了获得信号的功率。
最后将两个支路的功率相加，只有当本地伪码和本地载波的频率相位都与中频信号相同时，最后得到的功率才很大，否则结果近似为零。
根据这个结论考虑到噪声的干扰，在实际设计时应该设定一个判定门限，当两路信号功率和大于设定的门限时则判定为捕获成功，转入跟踪过程，否则继续扫描其它的频率或相位。
二、MATLAB仿真过程及结果仿真条件设置：抽样频率16MHz，中频5MHz，采样时间1ms，频率搜索步进1khz，相位搜索步进1chip，信号功率-200dBW，载噪比55dB（1）中频信号产生卫星导航信号采用数字nco的方式产生，如图2所示。
载波nco控制字为：carrier_nco_word=round(f_carrier*2^N/fs); 伪码nco控制字为：code_nco_word=round(f_code*2^N/fs);图 2其中载波rom存储的是正弦信号的2^12个采样点，伪码rom存储长度为2046的卫星伪码。
卫星导航定位算法与程序的设计——实验报告
2013 级测绘工程专业卫星导航定位算法与程序设计实验报告实验名称：卫星导航基本程序设计班级：学号：姓名：实验时间： 2016年6月28日~2016年6月30中国矿业大学目录实验一时空基准转换 (2)一、实验目的 (2)二、实验容 (2)三、实验过程 (2)四、实验感想 (6)实验二 RINEX文件读写 (7)一、实验目的 (7)二、实验容 (7)三、实验过程 (7)实验三卫星轨道计算 (12)一、实验目的 (12)二、实验容 (12)三、实验过程 (12)四、实验感想 (15)实验一时空基准转换一、实验目的1、加深对时空系统及其之间转换关系的理解2、掌握常用时空基准之间的转换模型与软件实现3、每人独立完成实验规定的容二、实验容本实验容包括：容一：编程实现GPS起点1980年1月6日0时对应的儒略日容二：编程实现2011年11月27日对应的GPS周数与一周的秒数容三：在WGS84椭球的条件下，编程实现当中央子午线为117度时，计算高斯坐标x = 3548910.811290287, y = 179854.6172135982 对应的经纬度坐标？容四：WGS84椭球下，表面x=-2408000; y=4698000;z= 3566000处的地平坐标系坐标为: e=704.8615;n=114.8683;u=751.9771的点对应的直角坐标为多少？三、实验过程1.针对第一、二部分容：1.1解决思路：先建立” TimeStruct.h”的头文件，将格里高利历、GPS 时间结构、儒略日时间结构共结构体的方式放在里面；在建立“TimeTr”的头文件，建立类“CTimeT r”，创建变量“GPS Time”、“Time”、”JulDay”，并且申明函数“TIME2JUL”、“TIME2GTIME”等，用这些函数分别实现所需要的转换。
1.2具体的实现函数：“TIME2JUL”函数：double CTimeTr::TIME2JUL()//TIME Time,JULIANDAY &JulDay{double m,y;double D;//h =Time.byHour+Time.byMinute/60.0+Time.dSecond/3600.00;if(Time.byMonth<=2){y=Time.wYear-1;m=Time.byMonth+12;}else{y=Time.wYear;m=Time.byMonth;}D=floor(365.25*(y+4716))+floor(30.6001*(m+1))+Time.byDay+Time.byHour/24.0-1537 .5;JulDay.lDay = int(D);JulDay.lSecond = D-int(JulDay.lDay);return 0;}“TIME2GTIME”：void CTimeTr::TIME2GTIME(){double JD;long m,y;int WN;double Wsecend;//UT=Time.byHour+Time.byMinute/60.0+Time.dSecond/3600.00;if(Time.byMonth<=2){y=Time.wYear-1;m=Time.byMonth+12;}else{y=Time.wYear;m=Time.byMonth;}JD=int(365.25*y)+int(30.6001*(m+1))+Time.byDay+Time.byHour/24.0+1720981.5;WN = floor((JD-2444244.5)/7.0);GpsTime.lWeek=WN;Wsecend=(JD-2444244.5-7*WN)*604800;GpsTime.lSecond=Wsecend;}1.3实验结果：2 针对第三部分容：2.1解决思路：运用实验指导书中提供的matlab高斯反算的代码，进行解算；将高斯反算的公式直接输成matlab代码，绕后在函数“function [B,L] = gauss_fansuan （x,y,L0）”中，将坐标x = 3548910.811290287，y = 179854.6172135982，L0 = 117，带入函数的坐边，即可得到所需要的经纬度。
北斗并行频率捕获算法研究与分析
北斗并行频率捕获算法研究与分析作者：陈伟强苏凯雄来源：《现代电子技术》2014年第19期摘要：采用FFT并行频率搜索的北斗二代卫星信号捕获算法，分析了并行频率捕获原理，利用Matlab对设计的捕获模型进行仿真，通过对比分析得到不同的FFT点数、子相关器个数及抽头设置等参数设置对捕获速度、频率搜索范围及灵敏度等性能的影响。
在实际设计中考虑到资源消耗和性能的平衡，给出改进后的捕获系统的FPGA实现方案。
改进传统的相关器结构，采用同一相关器分段积分方式模拟多个子相关器，减少硬件资源开销，提高了相关器使用效率。
关键词：北斗；捕获系统； FFT；并行频率搜索中图分类号： TN967⁃34 文献标识码： A 文章编号： 1004⁃373X（2014）19⁃0075⁃04 Study and analysis on parallel frequency acquisition algorithm of CompassCHEN Wei⁃qiang， SU Kai⁃xiong（Institute of Physics and Information Engineering， Fuzhou University， Fuzhou 350002，China）Abstract： The principle of parallel frequency acquisition is analyzed with Compass⁃2 satellite signal acquisition algorithm of FFT parallel frequency search. The designed acquisition model was simulated with Matlab. The influence of the different parameter settings of FFT spot number，sub⁃correlator number and tapping winning on acquisition speed， frequency search range and sensitivity performance were obtained by comparison and analysis. The balance of resource consumption and performance were considered in the actual design. The FPGA implementation scheme of the improved acquire system is given in this paper. Structure of the traditional correlator was improved， in which the multiple sub⁃correlators are simulated with the subsection integral mode of a same correlator to reduce the spending of hardware resource and improve the efficiency of the correlator.Keywords： Compass； acquisition system； FFT； parallel frequency searching0 引言由于卫星和接收机存在相对运动，接收机收到的卫星信号中会有一个多普勒频率偏移量。
数字中频GPS信号的MATLAB仿真
数字中频GPS信号的MATLAB仿真1杨勇，陈偲，王可东北京航空航天大学宇航学院，北京 (100083)E-mail：wangkd@摘要：文章以INS/GPS紧耦合为应用对象，在分析中频GPS信号结构的基础上，根据实际环境和载体运行状态，给出GPS信号延时、多普勒频移和钟差等参数，并应用中频信号解析表达式实现多颗卫星信号的合成。
最后，基于MATLAB语言进行了仿真计算，仿真结果表明信号符合实际情况，同时经过软件接收机的捕获、跟踪和解调计算，验证了信号的正确性。
关键词：GPS；高动态；紧耦合；中频；信号模拟器中图分类号：TP3911．引言随着固体弹道导弹射程的不断增加和打击精度的要求提高，纯惯性导航早已不能够满足要求。
全球定位系统（GPS）和惯性系统（INS）相结合是复合制导的重要发展方向之一，而对于GPS/INS组合导航来说，为了缩短研制周期，便于新信号开发及测试，软件信号模拟器和接收机的研究成为重要的研究方向之一。
GPS技术成长非常迅速，现在市场上的手持式GPS接收机已相当普遍，但是国内的自主知识产权的GPS技术产品的研发仍然比较薄弱，尤其是核心芯片的知识产权很少被国内所拥有。
国内的“北斗”、“GALILEO”导航定位都处在发展之中，信号模拟器的研究被越来越多的被重视。
信号模拟器具有成本低、可重复性好、数据完整等优点，不仅能用于组合导航技术研究，也能为新信号的验证研究提供支持，还可以为硬件接收机的接收性能测试提供有效的信号环境模拟。
GPS信号模拟器是软件无线电研究的一个方面，为处于设计阶段的GPS接收机提供仿真环境。
常见的GPS信号生成器产生的是射频信号，而目前接收机的设计重点侧重于基带信号处理，也就是本文中提到的数字中频GPS信号。
数字中频GPS信号模拟器目前主要是仿真载体运动、模拟时钟偏差、卫星星钟误差、电离层误差、对流层误差、多路经效应、天线的方向、弹体振动、以及云层、雷雨等实际环境对GPS信号的影响，并对接收机前端的下变频、滤波、采样和自动增益控制进行仿真，直到生成GPS接收机信号处理所需的数字信号。
基于MatLabSimulink的GPS系统仿真
(2) 我
δS = S − S 0 = δS d + δS w
引入的相位延迟为
δφ d =
(6)
们 在 仿 真 中 不 妨 采 用 一 9 位 LFSR
生成多项式为
其中 Tk 为绝对温度 P 为大气压(mbar) e0 为水汽风压(mbar) S 为实际传播路径 S0 为信号在真空中的传播路径 hs 为 90o 20o 15o 10o 5o 时 δS 的典型值分别为 2.51m 7.29 m 9.58 m 14.04 m 25.82 m 无线信道由于多径效应产生时延扩展 若收发信机处
k =1
N
ˆ(t) x
航电文经伪随机序列扩频并调制后的信号
则 L1 和 L2 载
N
波上的 GPS 信号可分别简单的建模为 S L1 (t ) = AP Pi (t )Wi (t )Di (t ) cos(w1t + ϕ1 ) + AC Ci (t )Di (t ) sin(w1t + ϕ1 )
S L 2 (t ) = B P Pi (t )W i (t ) D i (t ) cos( w 2 t + ϕ 2 )
• 1857 •
第 18 卷第 7 期 2006 年 7 月
系 统 仿 真 学 报 其中 hs 为 GPS 卫星相对观测站的高度角
11 12
Vol. 18 No. 7 July, 2006
组成 生成多项式分别为
X 1A ( x) = 1 + x + x + x + x
6 8
δS d = 1.552 × 10 − 5
[1]
ˆ (t ) − rq (t ) x ˆ (t ) 们采用 Jake 移动信道模型[7] y (t ) = ri (t ) x
自动实验一——典型环节的MATLAB仿真报告
自动实验一——典型环节的MATLAB仿真报告引言：典型环节的MATLAB仿真是一种常见的模拟实验方法，通过使用MATLAB软件进行建模和仿真，可以有效地研究和分析各种复杂的物理系统和控制系统。
本报告将介绍一个典型环节的MATLAB仿真实验，包括实验目的、实验原理、实验步骤、实验结果和讨论等内容。
一、实验目的本实验旨在通过MATLAB仿真实验，研究和分析一个典型环节的动态特性，深入了解其响应规律和控制方法，为实际系统的设计和优化提供理论支持。
二、实验原理典型环节是控制系统中的重要组成部分，一般包括惯性环节、惯性耦合和纯滞后等。
在本实验中，我们将重点研究一个惯性环节。
惯性环节是一种常见的动态系统，其特点是系统具有自身的动态惯性，对输入信号的响应具有一定的滞后效应，并且在输入信号发生变化时有一定的惯性。
三、实验步骤1.建立典型环节的数学模型。
根据实际情况，我们可以选择不同的数学模型描述典型环节的动态特性。
在本实验中，我们选择使用一阶惯性环节的传递函数模型进行仿真。
2.编写MATLAB程序进行仿真。
利用MATLAB软件的控制系统工具箱，我们可以方便地建立惯性环节的模型，并利用系统仿真和分析工具进行仿真实验和结果分析。
3.进行仿真实验。
选择合适的输入信号和参数设置，进行仿真实验，并记录仿真结果。
4.分析实验结果。
根据仿真结果，可以分析典型环节的动态响应特性，比较不同输入信号和控制方法对系统响应的影响。
四、实验结果和讨论通过以上步骤，我们成功地完成了典型环节的MATLAB仿真实验，并获得了仿真结果。
通过对仿真结果的分析，我们可以得到以下结论：1.惯性环节的响应规律。
惯性环节的响应具有一定的滞后效应，并且对输入信号的变化具有一定的惯性。
随着输入信号的变化速度增加，惯性环节的响应时间呈指数级减小。
2.稳态误差与控制增益的关系。
控制增益对稳态误差有重要影响，适当调整控制增益可以减小稳态误差。
3.不同输入信号的影响。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
北斗卫星导航信号串行捕获算法MATLAB仿真报告(附MATLAB 程序)北斗卫星导航信号串行捕获算法MATLAB仿真报告一、原理卫星导航信号的串行捕获算法如图1所示。
××∑∫( )2本地PRN发生器∫( )2本地载波发生器GPS中频信号×判决检波器≥VT？Yes 转跟踪NO 继续搜索图1 卫星导航信号的串行捕获算法接收机始终在本地不停地产生对应某特定卫星的本地伪码，并且接收机知道产生的伪码的相位，这个伪码按一定速率抽样后与接收的GPS中频信号相乘，然后再与同样知晓频率的本地产生的载波相乘。
GPS中频信号由接收机的射频前端将接收到的高频信号下边频得到。
实际产生对应相位相互正交的两个本地载波，分别称为同相载波和正交载波，信号与本地载波相乘后的信号分别成为，产生同相I支路信号和正交的Q 支路信号。
两支路信号分别经过一个码周期时间的积分后，平方相加。
分成两路是因为C/A码调制和P码支路正交的支路上，假设是I支路。
当然由于信号传输过程中引入了相位差，解调时的I支路不一定是调制时的I支路，Q支路也一样，二者不一定一一对应，因此为了确定是否检测到接收信号，需要同时对两支路信号进行研究。
相关后的积分是为了获取所有相关数据长度的值的相加结果，平方则是为了获得信号的功率。
最后将两个支路的功率相加，只有当本地伪码和本地载波的频率相位都与中频信号相同时，最后得到的功率才很大，否则结果近似为零。
根据这个结论考虑到噪声的干扰，在实际设计时应该设定一个判定门限，当两路信号功率和大于设定的门限时则判定为捕获成功，转入跟踪过程，否则继续扫描其它的频率或相位。
二、 MATLAB 仿真过程及结果仿真条件设置：抽样频率16MHz ，中频5MHz ，采样时间1ms ，频率搜索步进1khz ，相位搜索步进1chip ，信号功率-200dBW ，载噪比55dB（1） 中频信号产生卫星导航信号采用数字nco 的方式产生，如图2所示。
载波nco 控制字为：carrier_nco_word=round(f_carrier*2^N/fs); 伪码nco 控制字为：code_nco_word=round(f_code*2^N/fs);32位Adder12位载波rom模2046计数器伪码rom32位Adder Divide by 2^20溢出时输出脉冲carrier_nco_wordcode_nco_wordfsample××＋幅度加性噪声图 2其中载波rom 存储的是正弦信号的2^12个采样点，伪码rom 存储长度为2046的卫星伪码。
这样伪码采用2psk 的方式调制到射频，加性噪声很小是理想接收中频信号如图3所示。
-10图3 理想中频信号（2）噪声功率估计实际接收机接收到的导航信号淹没在噪声中，本程序对接收到的信号进行了噪声估计并进行了放大。
采用滑动平均估计法估计噪声功率，滑动平均估计法原理如图4所示。
图4 噪声功率滑动平均估计法原理迭代滤波器因子取0.8.功率估计结果是-191.48dBW 。
仿真中将接收中频信号放大到了signal_power_dB=-4.94dBW 。
这个功率与后面的判决门限有关系。
（2） 检测门限的确定常见的检测方法有幅度检波、平方检波和平方律检波。
幅度检波器的输出为在H0假设下，z(k) 服从瑞利分布，其概率密度函数为:幅度检波器NQ ∑=iN n niQ N 121迭代低通滤波器标度因子比较器+-σTV 比较结果XI psQ )()()(22k Q k I k Z +=2在H(1)假设下，z(k)服从莱斯分布，其概率密度函数为式中，为零阶修正的贝塞尔函数。
平方律检波输出为：在H0假设下，z(k)服从自由度为 2M 的 伽马分布，其概率密度函数为在H1假设下，z(k)服从自由度为 2M 的 卡方分布，其概率密度函数为当 M=1时，平方检波累积器就变成平方律检波器 ，可以计算出当归一化门限为Vt 时其虚警概率为：其中 采用恒虚警率检测，设虚警率为pfa ，本仿真取0.1 ，采用平方律检波，归一化判决门限为Vt=(-2*log2(pfa))^0.5，实际判决门限为VT=Vt*signal_power （3） 判决算法常见判决算法有单次判决、M/N 判决、（M/N+1）判决和Tong 判决，采用单次判决,虚警率为pfa=0.1 . 归一化检测门限为Vt=2.5776，判决门限为VT=0.8248； （5）仿真结果2210222000(/)exp[]() 02kk k k kk k Z Z D D Z f Z H I Z σσσ+=-⋅≥)(0⋅I ∑=+=Mi i Q i I k Z 122)]()([)(0)2exp()()2(1)/(21200≥-⋅Γ=-k kM k M k Z Z Z M H Z f σσ11/22112220001()(/)()exp()() 022M k k k k M k Z Z Z f Z H I Z λλσλσσ--+=⋅-⋅≥02201(/)exp[]221exp 2TTkfa k k k V V t Z P f Z H dZ dZ V σσ+∞+∞==-⎧⎫=-⎨⎬⎩⎭⎰⎰2t T V V σ=搜索21个多普勒频点和40和相位点，仿真设置接收中频为4.991MHz ，相位为2，结果如图5所示。
4.99555.0055.01x 1061020304050100150dopplercodephase图中最大处的相关结果是144，其他非峰值最大的约为1.约21dB。
部分相关值如下表：可以看出绝大部分数值都在门限之下，但也存在若干个在门限之上的数值，这些点可能造成虚警。
附：仿真主程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%´®ÐⶻñËã·¨·ÂÕæ%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%clear all;clc;% f_doppe=5000;fsample=21e6;f_m=5e6;t_sim=1e-3/1;df=1/t_sim;cnt_det=1; %Åоö´ÎÊýf_carrier=f_m-9*df; %½ÓÊÕÔز¨ÆµÂÊf_code=2.046e6; %αÂëËÙÂÊcode_phase_init=2; %Ê×´ÎËÑË÷ʱ½ÓÊÕÂëÏàλcode_phase_init_cmp=code_phase_init;f_local_init=f_m-10*df; %±¾µØÔز¨nco³õʼƵÂÊlocal_code_phase=3; %%%%±¾µØαÂëÏàλcnt_doppler=21;cnt_phase=40;dot_num=t_sim*fsample; %·ÂÕæÒ»¸öαÂëÖÜÆڵĵãÊýdBW_signal_pow=-200;dB_C_I=60;Am=10^(dBW_signal_pow/20)*2^0.5;pre_noise_power=dBW_signal_pow-dB_C_I+(10*log10(f_carrier));corr_result1=zeros(cnt_doppler,cnt_phase);for dect_num=1:cnt_det %Åоö´ÎÊýfor num_phase=1:cnt_phase %²éÕÒÏàλµãµãfor num_doppeler=1:cnt_doppler %µ¥´Î²éÕÒ¶àÆÕÀÕµãif num_doppeler==1f_local=f_local_init;elsef_local=f_local+df;endsignal_r=signal_gen(fsample,f_carrier,f_code,code_phase_init,dBW_sign al_pow,dB_C_I,dot_num); %%%%%ÕâÀïÔز¨Ö¸Ï±ßƵºóµÄÖÐƵ,ÊÕµ½ºó´¦Àí·¢Ë͵Ļ¬¶¯±êÖ¾£¬ËÑË÷ÏÂÒ»¸öÏàλ%%%%Ê×´ÎËÑË÷¹À¼ÆÐźŹ¦ÂÊif num_doppeler==1 & num_phase==1Ni=10;signal_r1=zeros(1,length(signal_r)+Ni);Qn1=zeros(1,length(signal_r)+Ni);Qn2=zeros(1,length(signal_r)+Ni);Qn3=zeros(1,length(signal_r)+Ni);pow_noise=zeros(1,10);for bbb=1:10for aaa=1:length(signal_r)/10signal_r1(aaa)= signal_r(aaa)*signal_r(aaa) ;a=0.8; %µü´úÂ˲¨Æ÷ϵÊý% Qn2(aaa+1)= signal_r1(aaa)*signal_r1(aaa);Qn1(aaa+Ni)=sum(signal_r1( (aaa):(aaa+Ni)))/Ni;%»¬¶¯Æ½¾ùQn(aaa)= Qn1(aaa+1)^0.5;Qn3(aaa+1)= a*Qn1(aaa)+(1-a)*Qn1(aaa+1);pow_noise1=Qn3(aaa+1);endpow_noise2(bbb)=pow_noise1;endpow_noise=sum(pow_noise2(1:10))/10;pow_noise_dB=10*log10(pow_noise);AD_min_volt=0.8;AD_R=1;AD_power=0.5*AD_min_volt*AD_min_volt/AD_R;AD_power_dB=10*log10(AD_power);Am1=AD_power_dB-pow_noise_dB;Am1=10^(Am1/20);endsignal_r=signal_r*Am1; %ÖÐƵ·Å´óAm_local=0.9;local_carrier=local_carrier_gen(fsample,dot_num,f_local,Am_local);%%%flocal°üº¬Á˶àÆÕÀÕlocal_code=local_code_gen(f_code,fsample,dot_num,local_code_phase)*Am _local;corr_result(num_phase,num_doppeler)=deal_local(signal_r,local_carrier ,local_code);% corr_result1(num_doppeler,cnt_phase)=corr_result;end% pulse_next_phase=1;code_phase_init=mod(code_phase_init+1,2046);% local_code_phase=local_code_phase+1;end%%%%²éÕÒµ¥´Î²¶»ñ×î´óÖµ¼«Î»ÖÃ[phase_max(dect_num)doppler_max(dect_num)]=find(corr_result==max(max(corr_result)));mod_max(dect_num)=max(max(corr_result));det_phase=local_code_phase-phase_max+1;det_doppler=(doppler_max-1)*df+f_m-10*df;disp(det_phase); disp(det_doppler);%%%%%%% %%%%%%% %%%%%%% %%%%%%% ÇóÅоöÃÅÏÞ ºãÐ龯ÂÊ%±ê¶ÈÒò×ÓÊÇÅоöÃÅÏÞÓëÔëÉù¹¦ÂʵıÈÖµfa=0.1;Vt=(-2*log2(fa))^0.5; %¹éÒ»»¯ÅоöÃÅÏÞ% Vt=0;%¾-¹ý·Å´óºó¹¦ÂÊΪAD_power% VT=( Vt*pow_noise );VT=Vt*AD_power ;if mod_max(dect_num)> VT %ƽ·½Âɼ첨flag_det=1;elseflag_det=0;end%Tong¼ì²âK=1;B=2;if flag_det==1K=K+1;elseK=K-1;endif K>=Bdisp('success!');elseif dect_num==cnt_detdisp('failed!');endendendcorr_result_dB=10*log10(corr_result/mod_max);cnt_doppler1=f_m-10*df:df:f_m+10*df;mesh(1:cnt_doppler,1:cnt_phase,corr_result);xlabel('doppler');ylabel( 'code_phase');11。