GPS系统仿真
第一章 系统仿真与GPSS
第一章系统仿真与GPSS系统仿真是利用系统模型的试验方法来研究现有的或计划的或设计的系统的有效性、合理性或经济性的过程。
仿真技术涉及到管理工程、系统工程、信息技术、控制理论、统计学、计算方法和计算机科学与技术等众多学科和领域的知识与应用，并随着这些学科领域的发展而不断发展和完善。
1.1基本概念1.系统与元素系统：系统是一个整体，由相互联系和相互依存的事物组成，它的范围由研究者根据所要解决问题的复杂性而定。
元素：也称为实体，即组成系统的基本事物。
在定义一个系统时，首先要确定系统的边界。
边界确定了系统的范围，边界以外对系统的作用称为系统的输入，系统对边界以外的环境的作用称为系统的输出。
尽管系统千差万别，但人们总结出描述系统的三个要素，即实体、属性和活动。
实体确定了系统的构成，属性描述每一实体的特征；活动表示了实体的行为以及它们之间的联系或作用。
例如：一个理发馆可以定义为一个管理系统。
构成理发馆系统的元素为理发师和顾客。
理发师与顾客相互联系而又相互依存，组成一个统一的整体。
理发师与顾客存在着服务与被服务的关系，同时理发师和顾客的存在不可或缺。
现实中的系统千差万别，常见的分类方式有以下3种：1.静态（static）和动态（dynamic）：静态系统与时间没有关系，但在动态系统中，时间却扮演着不可或缺的角色。
在绪论中介绍的蒲丰投针问题就属于静态系统的仿真，而第2章和第3章描述的理发馆系统则属于动态系统。
GPSS/JA V A就是为动态系统建立计算机模型的，因此，本书的主要研究动态模型。
2.连续（continuous）与离散（discrete）：在连续系统中，系统状态随时间连续变化，例如水库蓄水、放水以及出现降水和蒸发时水位的变化等。
我国发射的“嫦娥1号”探月飞船在发射过程中和绕月飞行过程中的状态也是时刻变化的，它也属于连续系统。
在离散模型中，系统的状态仅在离散的时间点上发生变化，例如在制造系统中，零件会在特定的时间到达和离开，机器会在特定的时间出现故障及被维修，工人也会在特定的时间开始休息。
GPS信号仿真器校准方法
校准 规范 。我们在 参 照这些相关 校准规 范及 生产 其 电平会被放大 ，增加 １ Ｂ ０ 。既然 内部噪声主要 ｄ 厂家 的校准手册 的基 础上 ，摸索 行之有效 的校 准 由中放 的第 一级产 生 ，因而输 人衰减器 不影 响内
方 法 ，对 信 号仿 真 器 的主要 项 目：频 率 、电平 、 部 噪声 电平 ，但 输入衰 减器影 响混频器 的信号 电 调制 等参数进行校准 。 平 ，并降低了信 噪 比。
对 于 Ｇ Ｓ仿 真器通 常选 取 １５ ５４ Ｐ ７ ．２ＭＨｚ 和 １２． ｚ ７０ ２ ６ ＭＨ 两个 常用 的频率校准点 。
个 别仪器 仍 旧会 碰到信 号抖动 的情况 ，可 以 使用求平均值的功能来读得相对稳定的信号 。
２ ． 电平准确度校准 ２
通 过 以 上 方 法 ，可 以准 确 、 快 速 地 在 指 定 频 与测 试 合 成 频 率 类 似 ，为 了模 拟 地 面 收 到 率点上测得电平值 。 的 Ｇ Ｓ 信号 ，Ｇ Ｓ 真器 的最 大输 出电平一般 ２３ 数字调制参数校准 Ｐ 小 Ｐ 仿 －
工 业 出版 社 ，２ ０ ． ０ ６
常 ；反 之 ，Ｅ Ｍ 正 常 ，Ｒ ｏ 多 正 常。 因此 ，对 Ｖ ｈ大 ＱＳ Ｐ Ｋ调制而 言 ，Ｒ ｏ只是 反映调制 质量 的一个 侧 ｈ 面 ，仅 以 Ｒ ｏ 判 断 调制 是 否 正 常 是不 全 面 的。 ｈ来
［ 吴幼璋 ，赵海 宁 ，于汇东 ，等． ２ ］ 数字调制质量参数 的校准 和量值
溯源 北京： 计量学报 ，２ ｏ ，２ ３ ０ ５ ６（ ）：２ １２ ４ ７— ７．
如果 Ｒ ｏ ｈ 不与 Ｅ Ｍ及 星座 图相结合 ，往往 不能 比 Ｖ 较客 观地 、完 整地反 映 总 的调制 质 量 。在 Ｃ ＭＡ Ｄ 网络 的测量 中 ，普遍存在 Ｒ ｏ ｈ 合格而 Ｅ Ｍ不合格 Ｖ
第二章系统仿真与GPSS幻灯片1
5）离散系统与连续系统
当所有能够引起系统件时 则该系统就称为离散系统； 若所有这些事件都是连续的 ，则该系统就称为离散系统 则该系统称为连续系统；若一个系统发生的事件。部分为离 散的而部分为连续的，则称其为混合系统。大多数管理 系统可以视为动态离散随机的排队系统。 为动态离散随机的排队系统 系统
2.2 离散随机排队系统的模拟方法
蒙特卡洛法能够解决的二类静态或动态问题：
问题－ 问题－：某人每轮向靶子射10箭，已知其击中靶心 的概率为25％，问一轮中射中7箭的概率为多少？ 问题二：有一银行营业点打算添置一台自动存取款机 问题二 （12小时服务），顾客按一定的间隔时间到来，排 队接受服务，先来者先用，后来者后用，顾客不愿 在队列中等待太久，否则会离去。管理人员想了解 等待时间超过3分钟的顾客的比例为多少，若该比 例太大，则考虑再增设一台机器。
2.1.6 系统模型与仿真
什么是仿真： 什么是仿真：
仿真是对系统模型进行实验的过程，是一种基于模型 仿真是对系统模型进行实验的过程 的活动。 仿真的基本过程是：确立研究的目的 确立研究的目的后，先建立系统 仿真的基本过程是 确立研究的目的 建立系统 的模型，然后使用某种计算工具求解模型 求解模型，最后对模 的模型 求解模型 对模 拟的结果进行分析，以获得答案。 拟的结果进行分析 仿真技术即是指以相似原理、系统理论与方法、信息 指以相似原理、 仿真技术 指以相似原理 系统理论与方法、 技术以及仿真应用领域的有关专业知识为基础， 技术以及仿真应用领域的有关专业知识为基础，以计 算机系统或与应用有关的物理效应设备及仿真器为工 具，利用模型对系统（己有的或设想的）进行研究的 利用模型对系统（己有的或设想的） 一门多学科的综合性的技术。 一门多学科的综合性的技术
gpss仿真例子
GPSS基本模块分类(1) 与活动实体有关的模块A.产生活动实体的模块GENERATE A,B,C,D,E,F,G(分别为：到达间隔时间均值，到达间隔时间方差，第一个动态实体产生的时间，应产生动态实体的总数，优先级)注：必选A或DB.活动实体延时模块ADVANCE A, B（延迟时间均值必选项，延迟时间方差可选项）C.活动实体结束模块TERMINATE A（离开系统的活动实体个数）D.活动实体的转向模块TRANSFER A,B,C,D 注意: A,B,C,D 之间不能有空格A的转向方式：①，逗号（即缺省）为无条件转向B域给出的地址②小数表示去C域的百分比,1-A 为去 B 址的百分比, 通路是随机的。
③BOTH 先去B 域地址, 不成就去C 域地址，都不成就保留在原模块中。
B为域地址。
(1) 活动实体参数的赋值语句ASSIGN A,B,CA为需要赋值或改变参数的参数号或参数名，A域中的+，-号表示从原有参数中加上或减去B。
B为需要赋给参数A或从参数A中加上或减去的数值。
(1) 比较测试模块TEST O A,B,CO 为辅助码,必选项.必须是E(等于),G(大于),GE(大于等于), L(小于),LE(小于等于),或NE(不等于)。
A，B为要比较的内容，可以是名称,数字,字符串,SNA或SNA的参数。
C是比较结果为假时，活动实体要进入的模块号。
(2) 逻辑开关设置语句LOGIC O A （O表示辅助码，A为逻辑开关号）逻辑开关三种形式：LOGIC R A ;将第A号逻辑开关置1LOGIC S A ;将第A号逻辑开关置0LOGIC I A ;将第A号逻辑开关反转(3) 初始化语句(逻辑开关,矩阵实体) INITIAL A,B（B缺省为1）注：若只是初始化逻辑开关，则必须以LS开始模式，即INITIAL LS$one；初始化逻辑开关one为1。
当然也可初始化系统内其他保存值。
多入口多设施（涉及到存储器）某有2个入口,2台加油器的汽车加油站, 每个入口处加油的汽车到达间隔时间都为均匀分布，均值为100秒，方差为10秒。
《GPS基本原理及其Matlab仿真》课件第3章
图3-6 岁差旋转
首先作顺时针Givens转动：
cos sin 0
RZ () sin cos 0
0
0
1
（3-1）
其意义是以Z0轴为旋转轴，顺时针转动ξ角，使X0轴旋转并到达
通过（观测）平天极的子午面上。
第二步是作逆时针Givens转动：
cos 0 sin RY () 0 1 0
5）
黄道是指地球绕太阳公转时的轨道平面和天球表面相交的 大圆，即当地球绕太阳公转时，地球上的观测者所看到的太阳
在天球面上作视运动的轨迹。黄道平面和天球赤道面的夹角ε 称为黄赤交角，ε≈23.5°。
6）
黄极是指过天球中心且垂直于黄道平面的直线和天球表面
的交点。黄极也有黄北极（KN）和黄南极（KS）的区分。
大地高H——地面点沿椭球面法线到椭球面的距离。
因此，地面任意一点P的位置， 在地球坐标系中可表示为 地心空间直角坐标（X， Y， Z）或地心大地坐标（B，L， H）。
图3-9 地心大地坐标系
这两种坐标系的换算关系为：
X (N H ) cos B • cos L Y (N H ) cos B • sin L Z [N (1 e2 ) H ]sin B
Z
L arcsin Y X2 Y2
ce2 tan B 1 e2 tan 2
B
H X2 Y2 N cos B
（3-16）
式中，c=a2/b为极点处的子午线曲率半径；e′2=(a2-b2)/b2 为椭球第二偏心率。式（3-16）中，大地纬度B需迭代计算，
但其收敛速度很快，迭代4次后，大地纬度B的精度可达0.000
232621.448 46.815T 0.00059T 2 0.001813T 3
GPS软件接收机的仿真与实现汇总
毕业论文GPS软件接收机的仿真与实现学院：地质工程与测绘学院专业：测绘工程摘要随着GPS的升级和新的卫星导航系统的发展，相比较传统GPS接收机，GPS软件接收机具有的成本低、灵活性高等优点越来越突出。
它使用软件方法和少量硬件即可实现信号接收处理，可以直接由运行在微处理器上的Matlab程序完成信号处理，因此具有良好的灵活性、可移植性及可扩展性。
因此，研究GPS软件接收机的仿真平台具有重要意义。
本文重点对GPS软件接收机的捕获和跟踪部分进行了研究，并在Matlab 中进行了定位解算。
本文在掌握GPS软件接收机原理的基础上，实现了对信号的仿真、捕获、跟踪及定位。
捕获部分为了提高GPS软件接收机的定位速度和定位精度，选用了在Matlab 环境下执行时间短、性能高的并行码相位搜索捕获算法。
跟踪部分将码跟踪环和载波跟踪环组合在一起，降低了跟踪环路的复杂度。
载波跟踪环路则选用了对1800相位转换不敏感的Costas环，以保证载波跟踪环路对信号的正确跟踪。
最后在Matlab环境下，编写了捕获、跟踪和数据处理等程序，用软件方式实现了对用户的定位，并对定位结果进行了分析概括，验证了所有算法的可行性，讨论了不足之处，为后续软件接收机的相关研究工作奠定了良好的基础。
关键词：GPS，软件接收机，仿真，捕获，跟踪，同步目录摘要 (II)ABSTRACT ......................................................................................... 错误！未定义书签。
第一章绪论 . (1)1.1GPS发展概况与组成 (1)1.2GPS的组成 (1)1.2.1 GPS空间卫星星座部分 (2)1.2.2 地面控制部分 (2)1.2.3 用户设备部分 (3)1.3GPS接收机的发展概况 (3)1.4软件接收机的特点及国内外研究现状 (4)1.4.1 软件接收机的结构特点 (4)1.4.2 软件接收机的发展现状 (6)1.5课题研究的意义 (6)1.6论文研究的主要内容 (7)第二章GPS信号的产生和结构 (8)2.1GPS信号的产生 (8)2.2GPS信号结构 (9)2.2.1 载波信号 (9)2.2.2 C/A码和P码 (9)2.2.3 导航电文 (12)2.2.4 GPS卫星信号的调制 (13)2.3本章小结 (13)第三章GPS信号的捕获 (14)3.1GPS信号捕获原理 (14)3.2GPS软件接收机捕获算法 (14)3.2.1 串行搜索捕获算法 (14)3.2.2 并行频域搜索捕获算法 (15)3.2.3 并行码相位搜索捕获算法 (16)3.3本章小结 (18)第四章GPS信号的跟踪 (19)4.1解调过程 (19)4.2锁相环原理 (20)4.3载波跟踪 (21)4.4码跟踪 (22)4.5本章小结 (23)第五章GPS软件接收机的MATLAB实现 (25)5.1并行码相位搜索捕获算法的MATLAB实现及捕获结果 (25)5.2GPS信号跟踪的MATLAB实现及跟踪结果 (27)5.3软件接收机的定位结果 (33)5.4本章小结 (34)总结与展望 (36)致谢 ................................................................................................. 错误！未定义书签。
基于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
GPS卫星运动及定位matlab仿真_毕业设计论文
毕业设计论文GPS卫星运动及定位matlab仿真摘要全球定位系统是具有全球性、全能性、全天候优势的导航定位、定时和测速系统，现在在全球很多领域获得了应用。
GPS卫星的定位是一个比较复杂的系统，其包含参数众多，如时间系统、空间坐标系统等。
此次设计是针对卫星运动定位的matlab仿真实现，因要求不高，所以对卫星运动做了理想化处理，摄动力对卫星的影响忽略不计（所以为无摄运动），采用开普勒定律及最小二乘法计算其轨道参数，对其运动规律进行简略分析，并使用matlab编程仿真实现了卫星的运功轨道平面、运动动态、可见卫星的分布及利用可见卫星计算出用户位置。
通过此次设计，对于GPS卫星有了初步的认识，对于静态单点定位、伪距等相关概念有一定了解。
关键字：GPS卫星无摄运动伪距matlab仿真The movement and location of GPS satellite onMA TLABAbstract：Global positioning system is a global, versatility, all-weather advantage of navigation and positioning, timing and speed system, now there has many application in many fields.GPS satellite positioning is a complex system, which includes many parameters, such as time and space coordinates system.This design is based on the matlab simulation of satellite motion and location, because demand is not high, so to do the idealized satellite movement, and ignore the disturbed motion ( so call it non-disturbed motion ).Using the Kepler and least-square method for calculating the parameters of orbital motion, for the characteristics of motion to make a simple analysis, and use the matlab simulation to program achieve the orbital plane of satellite, the dynamic motion, the distribution of visible satellites and using visible satellites to calculate the users‟ home.Through the design have primary understanding for the GPS satellite, and understanding the static single-point, pseudorange and so on.Key words：GPS satellite non-disturbed motion pseudorange matlab simulation目录第一章前言 (1)1.1课题背景 (1)1.2本课题研究的意义和方法 (2)1.3GPS前景 (2)第二章 GPS测量原理 (4)2.1伪距测量的原理 (4)2.1.1 计算卫星位置 (5)2.1.2 用户位置的计算 (5)2.1.3 最小二乘法介绍 (5)2.2载波相位测量原理 (6)第三章 GPS的坐标、时间系统 (10)3.1坐标系统 (10)3.1.1 天球坐标系 (10)3.1.2 地球坐标系 (12)3.2时间系统 (13)3.2.1 世界时系统 (14)3.2.2 原子时系统 (15)3.2.3动力学时系统 (16)3.2.4协调世界时 (16)3.2.5 GPS时间系统 (16)第四章卫星运动基本定律及其求解 (18)4.1开普勒第一定律 (18)4.2开普勒第二定律 (19)4.3开普勒第三定律 (20)4.4卫星的无摄运动参数 (20)4.5真近点角的概念及其求解 (21)4.6卫星瞬时位置的求解 (22)第五章 GPS的MATLAB仿真 (25)5.1卫星可见性的估算 (25)5.2GPS卫星运动的MATLAB仿真 (26)结论 (38)致谢 (40)参考文献 (1)附录 (2)第一章前言1.1 课题背景GPS系统的前身为美军研制的一种子午仪卫星定位系统(Transit），1958年研制，64年正式投入使用。
GPSS管理仿真实验五 公交车站运行模拟
一、实验目的练习高级GPSS语言主要模块的使用，其中包括LOGIC、GATE、TEST、TABLE、TABULATE、QTABLE、LINK、UNLINK、SPLIT、ASSEMNLE、MATCH、GATHER等的应用。
二、实验环境Excel2007三、实验内容及步骤任务一以例8.3为基础，模拟一公共汽车站的等车过程。
（一）用文本编辑器编写GPSS程序打开文本编辑器，输入程序代码。
如图所示（二）使用GPSS软件安装GPSS软件后，可以通过“开始”按钮，选择“程序”菜单。
然后，打开GPSS软件程序，启动后的窗口如图所示。
点击“file”按钮，新建一个model。
将在文本编辑器中的代码复制在GPSS软件中的代码区域。
如图所示运行模拟程序，选择Command->Create Simulation创建模拟，选择Command->START在对话框中取代1，输入300选择ok当300个实体进入TERMINATE模块时，模拟就会结束。
这代表300个人通过了十字转门。
产生的标准输出结果如图所示。
从标准输入表中，可以看到当第300个观众通过十字转门时，时间过去了2134.023秒。
模拟程序多次执行会由于随机数的不同产生微小的差异。
可以看到在整个模拟过程中排队等待的人最多有3个人。
（三）GPSS World程序的基本组成从结构上而言，程序由四个部分组成，它们是模拟开始语句，定义语句，模拟模块和模拟结束控制语句。
1.模拟开始控制语句的两种实现方式（1）在程序的最前面加入Simulation（2）由Command->Create Simulation命令执行。
2.模拟结束控制语句的两种实现方式（1）由Command->START命令执行（2）在程序的结尾加入Start n，通过n来控制模拟的长度3.定义语句定义语句相当于说明语句，置于模拟程序之前。
4.模拟模块部分模拟模块是程序的主体，反应了模拟的过程。
基于MATLAB的GPS信号的仿真
摘要扩频通信是近几年来迅速发展起来的一种通信技术。
在早期研究这种技术的主要目的是为提高军事通信的保密和抗干扰性能，因此这种技术的开发和应用一直是处于保密状态。
扩频技术在军事应用上的最成功范例可以以美国和俄国的全球定位系统（GPS和GLONASS）为代表；在民用上GPS和GLONASS也都得到了广泛的应用，这些系统的基础就是扩频技术。
全球定位系统（GPS）用于对全球的民用及军用飞机、舰船、人员、车辆等提供实时导航定位服务。
GPS系统采用典型的CDMA体制，这种扩频调制信号具有低截获概率特性。
该系统主要利用直接序列扩频调制技术，采用的伪码有C/A码和P(Y)码两种。
本文讲述的是直接序列扩频通信技术在全球定位系统（GPS）中的应用。
主要介绍扩频通信中的伪码仿真，简要论述M序列和伪随机噪声码（P码和C/A码）及其产生，并使用MATLAB7.0仿真M序列、P码和C/A码的编码过程和仿真结果，介绍直扩频技术伪码的相关知识，重点介绍P码。
关键字：全球定位系统；直接扩频通信；伪码仿真AbstractSpread spectrum communication is a communications technology developed rapidly in recent years. In early studies the main purpose of this technology is to improve the military communications confidential and anti-jamming performance, therefore the development and application of this technology is always in secret state. Spread spectrum technology in the most successful military application examples are the United States and Russia could the global positioning system (GPS and GLONASS) for representative; In civil GPS and GLONASS also have been widely used，which foundation of system is the spread spectrum technology.Global positioning system (GPS) is used to provide real-time navigation and positioning services for global civil and military aircraft, ships, personnel, vehicles and so on. GPS system adopts the typical CDMA system, which kind of spread spectrummodulation signals have low intercept probability characteristic. This system mainly used the direct sequence spread spectrum modulation technology, using the PRN code including C/A code, P codes and Y codes.This article tells the direct sequence spread spectrum communication technology applied in global positioning system (GPS) .The article mainly introduces the pn code spread spectrum communication simulation, briefly discussing M sequence and pseudo random noise code (P yards and C/A yards) and its produce and use MATLAB7.0 simulate M series, P yards and C/A yards of encoding process and the simulation results, introducing pn code straight spread-spectrum technology knowledge, especially P yards.Key: GPS; DS-SS;Pn code simulation目录引言 (4)1GPS理论及其特性 (5)1.1GPS系统概述 (5)1.2GPS信号构成 (7)1.2.1M序列 (10)1.2.2C/A码 (13)1.2.3P码 (15)1.3小结 (20)2MATLAB软件 (21)2.1MATLAB软件简介 (21)2.2MATLAB应用概述 (22)2.2.1 MATLAB功能介绍 (22)2.2.2 MATLAB使用方法 (24)3 GPS卫星导航信号算法及其MATLAB仿真 (31)3.1 C/A码仿真代码及其仿真结果 (31)3.2 P码的仿真代码及其仿真结果 (34)3.3 结果分析及其相关性分析 (38)结论 (39)致谢 (40)参考文献 (41)附录A 英文原文 .................................................................................................. 错误!未定义书签。
INSGPS组合导航系统仿真分析
W E , N , U , a E , a N , 0, 0, 0, 0, 0, 0, 0, 0
T
E , N , U 陀 螺 沿 东 ， 北 ， 天 方 向 的 随 机 漂 移 aE , aN 加 速 度 计 沿 东 ， 北 方 向 的 随 机 漂 移
.
vE tan L ie sin L Rn h 0 vN Rm h
vE ie co s L Rn h vN Rm h 0
vN Rm h E E vN N N ie sin L L R h U n U vE vE 2 tan L ie co s L sec L L R h Rn h n
p
p
陀 螺 漂 移 引 起 的 平 台 等 效 误 差 转 动
上式有关项可写成：
t t t t p t t 于是平台误差角方程的矢量表达式为： ＝ － ie et ie et 综合以上平台误差角方程的矩阵表达式可写成：
it ie et
系统量测方程
INS的位置量测信息可表达为地理系下的真值与相应误差之和：
LI I hI Lt L t ht h
GPS接收机给出的位置量测信息，可表示为地理系下的真值与相应误差之差： 式中：
在程序开始时，必须要有初值才能进行卡尔曼滤波。卡尔曼滤波要求初值条件为：
X 0 m X m X 0 0 P0 C X P 0 C X 0 0
GPS卫星信号仿真和捕获算法研究
摘要GPS（Global Positioning System）是基于卫星的定位系统，被认定是当前定位导航设备中最重要的发展，在军事和民用中有着广泛的应用。
传统GPS接收机基于专用硬件，具有较高的灵敏度和信号处理速度，但是灵活性差，定型后便无法改动。
GPS软件接收机是利用通用处理器，对实时采集的GPS信号进行处理和计算。
其优点为具有高度的灵活性：可以方便地研究和实现新算法；也可以容易地对软件进行修改来适应未来的GPS信号，节约硬件成本。
因具有以上优点，GPS软件接收机的研发成为了当今世界的热门研究课题。
而研究GPS软件接收机的关键在于其接收设备的基带处理单元，包括捕获和跟踪环路两个部分。
本文的重点就是在Matlab平台上对捕获进行研究和实现。
文中首先介绍了GPS系统的组成：GPS系统由GPS卫星、地面监控系统、GPS 接收机三部分组成。
之后阐述了直接扩频通信以及GPS信号的结构与特性，接着在本文的主体部分，着重研究了GPS软件接收机的信号捕获技术。
对串行时域捕获算法，并行频域捕获算法，并行码相位捕获算法进行了比较，并选择并行码相位捕获算法实现软件接收而串行实现硬件接收。
由于文中所研究的接收机是静态接收机，因此捕获模块的捕获范围是载波多普勒频移5KHz，但由于C/A码多普勒频移非常小，故本仿真是忽略多普勒频移。
然后在Matlab软件平台上实现对GPS 信号捕获的仿真。
关键字：Simulink，GPS信号捕获，Matlab ，C/A码ABSTRACTGPS (Global optimisation techniques based on satellite Positioning System) is the Positioning System, was considered the current orientation navigation equipment, the most important development in military and civilian has extensive applications. Traditional GPS receiver special hardware, based on high sensitivity and signal processing speed, but low agility, finalize the design and can change. GPS software receiver is using gm to collect real-time processor, the GPS signal processing and calculation. Its advantage for a high degree of flexibility: can easily research and to realize the new algorithms; Also can easily modify of software to adapt to future GPS signal, managing hardware cost.Because of its advantages, GPS software receiver above the world development become the hot research topic. And the key research GPS software receiver lies in its receiving equipment, including the baseband processing unit capture and tracking loop two parts. This paper focuses on the Matlab is studied and realization of capture.This paper firstly introduces the GPS system composition: GPS system consists of GPS satellite, ground monitoring system, GPS receiver of three parts. Then expounds the spread spectrum communication and GPS signal directly, then structure and characteristics of the main body part in this article, this paper studies the GPS software receiver signal capturing technologies. Time domain capture of serial algorithm, and HangPin domain capture algorithm, parallel yards compares the phase capture algorithm, and select parallel code phase capture algorithm software to accept and serial realize hardware accept. Because the study is static receiver, so receiver module capture range is captured the doppler frequency modulated 5KHz, but because C/A code doppler frequency is very small, so the simulation is to ignore the doppler frequency. Then in Matlab software platform to achieve captured GPS signals simulation.KEY WORDS：Simulink, GPS signal acquisition, Matlab, C/A目录摘要 (I)第一章绪论 (1)1.1引言 (1)1.11 GPS的发展概况 (1)1.2课题研究的意义和目的 (2)1.3本文的内容和组织结构 (2)第二章GPS系统原理概述 (4)2.1GPS系统组成 (4)2.2GPS信号的组成 (5)2.2.1 载波信号 (5)2.2.2 伪随机码信号 (6)2.2.3 导航电文 (7)2.3GPS扩频通信系统基本原理 (8)2.4GPS信号的MATLAB仿真及实现 (9)2.4.1数据码产生器 (9)2.4.2 C/A码产生器 (10)2.4.3 C/A码互相关性和自相关性验证 (11)第三章GPS信号的捕获 (13)3.1多普勒效应对C/A码相位和载波频率的影响 (13)3.2GPS信号的捕获 (14)3.2.1 串行捕获算法 (15)3.2.2 并行频率捕获算法 (16)3.2.3 并行码相位捕获算法 (17)并行码相位捕获算法，是在频域里执行循环卷积。
基于MATLAB的GPS信号的仿真研究
基于MATLAB的GPS信号的仿真研究一、本文概述随着全球定位系统（GPS）技术的广泛应用，其在导航、定位、授时等领域的重要性日益凸显。
为了更好地理解GPS信号的特性，提高GPS接收机的设计水平和性能，对GPS信号进行仿真研究显得尤为重要。
本文旨在探讨基于MATLAB的GPS信号仿真方法，分析GPS信号的特点，以及如何利用MATLAB这一强大的数值计算环境和图形化编程工具，对GPS信号进行高效、精确的仿真。
文章首先介绍了GPS系统的发展历程、基本原理和信号特性，为后续的信号仿真提供了理论基础。
随后，详细阐述了GPS信号仿真的一般流程，包括信号生成、传播模型、噪声添加等关键环节。
在此基础上，重点介绍了如何利用MATLAB编写GPS信号仿真程序，包括信号生成、传播模型建立、噪声模拟等方面的具体实现方法。
文章还通过实际案例，展示了基于MATLAB的GPS信号仿真在接收机设计、性能评估等方面的应用。
通过仿真实验，可以深入了解GPS信号在不同环境下的传播特性，为接收机算法优化和性能提升提供有力支持。
本文的研究不仅有助于加深对GPS信号特性和仿真方法的理解，也为GPS接收机的研究和开发提供了一种有效的技术手段。
通过MATLAB的仿真研究，可以更加直观地揭示GPS信号的本质规律，为实际应用提供有力指导。
二、GPS信号原理及特性全球定位系统（GPS）是一种基于卫星的无线电导航系统，它利用一组在地球轨道上运行的卫星来提供全球范围内的定位和时间服务。
每个GPS卫星都不断地向地球表面发射射频信号，这些信号被地面上的接收器接收并处理，从而确定接收器的三维位置和速度，以及精确的时间信息。
GPS卫星发射的信号是L波段的射频信号，分为两个频段：L1（142 MHz）和L2（160 MHz）。
每个频段都包含两种类型的信号：C/A码（粗捕获码）和P码（精密码）。
C/A码是对公众开放的，用于民用和商业应用，而P码则用于军事和特定的高精度应用。
基于GPS软件接收机的SINS／GPS紧组合系统仿真
能差 等缺 点 ； 超 紧 组合 实 现 了 双方 向辅 助 ， 即 ＧＰ Ｓ
的输 出与 Ｓ Ｉ ＮＳ进 行 组 合 滤 波 ， 修整 Ｓ Ｉ ＮＳ的 累 积 误差 ， Ｓ Ｉ ＮＳ的输 出 也 用 于辅 助 ＧＰ Ｓ的 捕 获 和跟 踪， 解决 了接 收环路 动 态范 围窄 和抗干 扰性 能差 等
收 稿 日期 ：２ ０ １ ２ － ０ ６ － ２ ９
后输 出 ， 而 系统校 正 分 为 输 出校 正 与 反馈 校 正 ， 即
开环 校 正与 闭环校 正 。
资 助 项 目 ：国 家 航 空 科 学 基 金 （ 批 准号 ： ２ ０ １ ０ ０ ８ ５ １ ０ １ ７ ）
联 系人 ： 胡雯婷 Ｅ － ｍａ ｉ ｌ ： ８ ０ ４ ９ ２ １ ２ ７ ８ ＠ｑ ｑ ． ｃ ｏ ｍ
第３ ８卷 第 ５期
２ ０ １ ３年 １ Ｏ月
全 球 定 位 系 统
Ｇ Ｎ ＳＳ Ｗ ｏｒ ｌ ｄ ｏｆ Ｃｈｉ ｎａ
Ｖｏ １ ． ３８， ＮＯ ．５ Ｏｃ ｔ ｏｂ ｅ ｒ ， ２ ０１ ３
基 于 ＧＰ Ｓ软 件 接 收 机 的 Ｓ Ｉ ＮＳ ／ ＧＰ Ｓ 紧 组合 系统 仿 真
性 能差 等 缺点无 法避 免 ， 使得组 合 系统 的可靠 性主 要取决于 Ｓ Ｉ ＮＳ ； 紧组 合是 基 于 伪距 和伪 距 率观 测 的组 合 ， 当可见 卫 星数少 于 ４颗 时也能完 成组 合滤
波， 但 并不 能解 决 Ｇ Ｐ Ｓ的 动 态 范 围 窄 和抗 干扰 性
在 图 １中 ， ＧＰ Ｓ模 拟 信 号 和 惯 导 原 始 数 据 分
： 软件部分
GPS仿真模拟
声源定位（vb程序）Public X, Y As DoublePublic A1, B1, C1, p1, Q1, R1, A2, B2, C2, p2, Q2, R2, A3, B3, C3, p3, Q3, R3, d1, d2, d3 As DoublePrivate Sub Command1_Click()Dim x1, y1, x2, y2, x3, y3, x4, y4, c, t1, t2, t3, t4c = 2980 '声速值x1 = 0 '接收传感器坐标值y1 = 450x2 = 300y2 = 450x3 = 0y3 = 0x4 = 300y4 = 0t1 = Text1.Text '输入测得的t1,t2,t3,t4t2 = Text2.Textt3 = Text3.Textt4 = Text4.Textd1 = c * (t1 - t3) '3个距离差d2 = c * (t2 - t3)d3 = c * (t4 - t3)A1 = x2 * (x1 ^ 2 + y1 ^ 2 - d1 ^ 2) - x1 * (x2 ^ 2 + y2 ^ 2 - d2 ^ 2) '用1,2点求声源坐标B1 = y2 * (x1 ^ 2 + y1 ^ 2 - d1 ^ 2) - y1 * (x2 ^ 2 + y2 ^ 2 - d2 ^ 2)C1 = d1 * (x2 ^ 2 + y2 ^ 2 - d2 ^ 2) - d2 * (x1 ^ 2 + y1 ^ 2 - d1 ^ 2)cosp1 = A1 / Sqr(A1 ^ 2 + B1 ^ 2)sinp1 = B1 / Sqr(A1 ^ 2 + B1 ^ 2)cosQ1_p1 = C1 / Sqr(A1 ^ 2 + B1 ^ 2)sinQ1_p1 = Sqr(1 - cosQ1_p1 ^ 2)cosQ1 = cosQ1_p1 * cosp1 - sinQ1_p1 * sinp1sinQ1 = sinQ1_p1 * cosp1 + sinp1 * cosQ1_p1R1 = (x1 ^ 2 + y1 ^ 2 - d1 ^ 2) / (2 * x1 * cosQ1 + 2 * y1 * sinQ1 + 2 * d1)A2 = x2 * (x4 ^ 2 + y4 ^ 2 - d3 ^ 2) - x4 * (x2 ^ 2 + y2 ^ 2 - d2 ^ 2) '用2,4点求声源坐标B2 = y2 * (x4 ^ 2 + y4 ^ 2 - d3 ^ 2) - y4 * (x2 ^ 2 + y2 ^ 2 - d2 ^ 2)C2 = d3 * (x2 ^ 2 + y2 ^ 2 - d2 ^ 2) - d2 * (x4 ^ 2 + y4 ^ 2 - d3 ^ 2)cosp2 = A2 / Sqr(A2 ^ 2 + B2 ^ 2)sinp2 = B2 / Sqr(A2 ^ 2 + B2 ^ 2)cosQ2_p2 = C2 / Sqr(A2 ^ 2 + B2 ^ 2)sinQ2_p2 = Sqr(1 - cosQ2_p2 ^ 2)cosQ2 = cosQ2_p2 * cosp2 - sinQ2_p2 * sinp2sinQ2 = sinQ2_p2 * cosp2 + sinp2 * cosQ2_p2R2 = (x4 ^ 2 + y4 ^ 2 - d3 ^ 2) / (2 * x4 * cosQ2 + 2 * y4 * sinQ2 + 2 * d3)A3 = x4 * (x1 ^ 2 + y1 ^ 2 - d1 ^ 2) - x1 * (x4 ^ 2 + y4 ^ 2 - d3 ^ 2) '用1,4点求声源坐标B3 = y4 * (x1 ^ 2 + y1 ^ 2 - d1 ^ 2) - y1 * (x4 ^ 2 + y4 ^ 2 - d3 ^ 2)C3 = d1 * (x4 ^ 2 + y4 ^ 2 - d3 ^ 2) - d4 * (x1 ^ 2 + y1 ^ 2 - d1 ^ 2)cosp3 = A3 / Sqr(A3 ^ 2 + B3 ^ 2)sinp3 = B3 / Sqr(A3 ^ 2 + B3 ^ 2)cosQ3_p3 = C3 / Sqr(A3 ^ 2 + B3 ^ 2)sinQ3_p3 = Sqr(1 - cosQ3_p3 ^ 2)cosQ3 = cosQ3_p3 * cosp3 - sinQ3_p3 * sinp3sinQ3 = sinQ3_p3 * cosp3 + sinp3 * cosQ3_p3R3 = (x1 ^ 2 + y1 ^ 2 - d1 ^ 2) / (2 * x1 * cosQ3 + 2 * y1 * sinQ3 + 2 * d1)X = (R1 * cosQ1 + R2 * cosQ2 + R3 * cosQ3) / 3 '取平均求最合适声源坐标点Y = (R1 * sinQ1 + R2 * sinQ2 + R3 * sinQ3) / 3Text5.Text = X '输出结果Text6.Text = YEnd SubPrivate Sub Command2_Click() '清除Text1.Text = ""Text2.Text = ""Text3.Text = ""Text4.Text = ""End SubGPS模拟(vc程序)#include <iostream>#include <cmath>#include <iomanip>using namespace std;#define c 2982#define min 0.0000000001void GPS(double p[10] , double q[10] , double t[10] , double p1 , double q1 , int num){double a1,a2,a3,a4,b1,b2,mult,H=0,I=0;int i;while(1){a1 = a2 = a3 = a4 = b1 = b2 = 0;for(i = 0; i < 10; i ++){a1 = a1-(3*pow((p[i]-p1),2)+pow((q[i]/1000-q1),2)-pow(c*t[i]/1000000,2));a2 = a2+2*(p[i]/1000-p1)*(q1-q[i]/1000);a3 = a2;a4 = a4-(3*pow((q[i]-q1),2)+pow((p[i]/1000-p1),2)-pow(c*t[i]/1000000,2));b1=b1+(pow((p[i]/1000-p1),2)+pow((q[i]/1000-q1),2)-pow(c*t[i]/1000000,2))*(p1-p[i]/1000);b2=b2+(pow((p[i]/1000-p1),2)+pow((q[i]/1000-q1),2)-pow(c*t[i]/1000000,2))*(q1-q[i]/1000);} //end for;mult = a3/a1;a4 -= a2*mult;b2 -= b1*mult;I = b2/a4;H = (b1-a2*I)/a1;p1 += H;q1 += I;if((pow(H,2)+pow(I,2)) < min)break;} // end while;p1 = p1*1000;q1 = q1*1000;cout <<"点P"<< num <<"定位于: P = "<< fixed << setprecision(8)<< p1;cout <<"(mm) , Q = "<< fiped << setprecision(8)<< q1 <<"(mm)"<<endl;} // end GPS;int main(){double t1[10] , t2[10] , t3[10] , p[10] , q[10];cout << "请输入时间数据(um)与坐标数据(mm):" << endl;int i;for(i = 0; i <= 9; i ++){cout << "第" << i+1 << "组:" << endl;cin >> t1[i] >> t2[i] >> t3[i] >> p[i] >> q[i];} // end for;GPS(p,q,t1,0,0.45,1);GPS(p,q,t2,0.3,0.45,2);GPS(p,q,t3,0,0,3);system("pause");return 0;} // end main;。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
GPS系统仿真实验报告姓名：***学号：**********指导教师：***一、实验目的：分析扩频码的构成原理，基于MATLAB 产生GPS 信号，并进行调制传输，接收解调，最终建立完整的GPS 系统并对接收的结果进行分析。
二、实验原理：扩频通信：扩频通信的基本理论依据是信息论中的香农公式)1(log 2NS W C += 其中,C 是信道容量(bit/s),砰是信道带宽,S 是信号的平均功率,N 是噪声功率,香农公式表明了在给定信噪比的情况下,只要采用某种编码系统就能以任意小的差错概率,以接近于C 的传输速率来传送信息。
扩频序列采用的是伪随机序列,伪随机序列具有尖锐的自相关性,近似为零的互相关性这一特点。
基于这一特性,可以通过检测本地伪随机序列与原伪随机序列的互 相关峰值的大小来捕获伪随机序列。
伪随机序列包括m 序列、Gold 序列、M 序列和组合序列等。
GPS 信号中所使用的是一种组合Gold 序列(C/A 码/P 码)。
三、实验内容：1、GPS 信号产生：GPS 卫星同时发射 L1、L2 两路扩频信号，载频分别为 1575.42Mz 、1227.6MHz ，L1 信号将导航电文分别用 C/A 码和 P 码进行扩频。
GPS 卫星的导航电文(D 码),是一组二进制的数码序列,它与编码脉冲相对应,以数字通讯方式传递给用户。
当数码取“0”时,对应的编码脉冲状态取“+1”而码值取“1”,对应的编码脉冲取“-1”,即反像编码。
L2波段上只传送P 码,P 码是加密码,普通用户无法接受和使用。
C/A 码与P 码采用的也是反像编码方式，因此,当C/A 码或P 码与D 码进行模二和(异或运算)时,表现出来的是码元异或,而波形相乘。
2、C/A 码结构：C/A 码是一种Gold 组合码,由两个10级反馈移位寄存器组合产生,可以表示为)()()(021τi N t G t G t G +⊕=其中0τ为码元对应的时间：1/1023ms ，i N 为相位偏置的码元数。
C/A 码的产生如下图所示：10311)(t t t G ++=109863221)(t t t t t t t G ++++++=3、P 码：P 码是由两组各有两个12级反馈移位寄存器结合产生的,而每一组反馈移位寄存器又是由两个12级的m 序列发生器组合而成的。
因此,P 码仿真模型需用4个m 序列发生模块PN sequence Generator 产生,PN1和PN2组成1X 子码发生器;PN3和PN4组成2X 子码发生器。
又因为C/A 码与P 码采用的均是反像编码方式，所以同一组内及两组间的m 序列模二和(异或运算)时,表现出来的也是码元异或,而波形相乘1X 后接一个移位寄存器使P 码延时一个基码，仿真电路如下图所示。
P 码仿真频率设为l0.23MHz,仿真时间设为0.1ms ，子码a X 1系数为[1100101000001],b X 1的系数为[1111100100111],a X 2的系数为[1111110111011],b X 2的系数为[1001100011101]。
4、系统框架及实现：GPS 信号是数据码经扩频码.载波两次变频调制而成的,由卫星负责信号的发射;通过有噪信道的传输进入接收端;接收端负责信号的接收、捕获和跟踪,GPS 仿真系统结构图如图所示。
5、系统基本框架 6、GPS 信号的捕获与跟踪：信号捕获的目的是为了粗略的确定C/A 码在数据块中的起始位置、载波频率和码相位。
信号跟踪的目的是为了获得精确的载波频率和码相位,并且当载波频率和码相位出现一定变化时,能够使本地载波频率和本地C/A 码相位跟踪上这种变化。
信号跟踪分为载波频率/相位跟踪和码相位跟踪。
捕获和跟踪的实质都是一种对发射信号参数的估计。
捕获是在不具备接收信号参数准确值的信息情况下,对接收载波频率和发送码字及相位的粗略估计,跟踪的实质也是一种对发送信号参数的估计,是在具备对接收载波频率和发送码字相位的粗略估计的信息的前提下,使用更好的估计手段,对发送信号的载波频率/相位和发送码相位的更精确估计。
如果接收到的信号和本地复现信号的相位完全对齐,那么在相关过程中将会出现最大的峰值;但如果相位不完全对齐,只要偏移在一个码片内,仍然可以得到一个较大的峰值,所以捕获过程只能给出载波频率和码相位的一个粗略估计。
捕获/跟踪过程中最重要的是实现两个同步:1.码元同步2.载波同步。
在跟踪环节,码元同步负责信号解扩,载波同步负责信号解调。
此外,由于是软件仿真,两者还肩负着信号的捕获工作。
接收信号进行的处理框图和C/A 码的捕获跟踪实现如下图所示：s(t ) 数据源 编码器 m 序列发生器发射机 射频 振荡器{a n } d (t )c (t ) cos(2πf 0t +ϕ0) 中频 滤波器VCOm 序列 发生器射频 滤波器R (t )r (t )v (t ) 至数据 检测器2cos[2π(f 0+f IF +)t +](b)(a)接收后进行的函数处理C/A码的捕获和跟踪实现框图码元同步：当捕获到卫星信号后,本地码的频率和相位与接收到的码元存在较大的误差,因此必须进行调整,以达到本地码与卫星信号C/A码同频同相，这个调整过程由码跟踪环来完成,称之为相关接收或解扩。
四、实验结果：用m文件实现GPS系统的信号发送以及接收部分的仿真，产生两种伪随机序列后分别对两路信号进行调制，然后对C/A码进行捕获和跟踪，对得到的子帧进行识别再计算星历数据和伪距以及卫星的位置和用户位置。
得到的波形如下图所示,其中图2是svn=10的C/A码序列在码相位和载波频率二位平面上搜索的相关峰值：图2 svn=10的C/A 码序列在码相位和载波频率二位平面上搜索的相关峰值0500100015002000250030003500400045005000GPS = 10 max at50码采样点相关峰GPS = 10 Freq =1250000Hz频率搜索点相关峰1234567891011GPS = 10carrier pullin at Freq =1250000Hz频率牵引搜索点相关峰的平方（归一化）050010001500200025003000-10载波环鉴相器输出时间t(单位：ms)鉴相器输出50010001500200025003000-8载波NCO 频率控制字的改变时间t(单位：ms)频率控制字的改变值五、实验总结：本次实验首先对GPS 的信号进行了仿真，GPS 卫星同时发射 L1、L2 两路扩频信号，载频分别为 1575.42Mz 、1227.6MHz ，L1 信号将导航电文分别用 C/A 码和 P 码进行扩频。
GPS 卫星的导航电文(D 码),是一组二进制的数码序列,它与编码脉冲相对应,以数字通讯方式传递给用户。
当数码取“0”时,对应的编码脉冲状态取“+1”而码值取“1”,对应的编码脉冲取“-1”,即反像编码。
L2波段上只传送P 码,P 码是加密码,普通用户无法接受和使用。
C/A 码与P 码采用的也是反像编码方式，因此,当C/A 码或P 码与D 码进行模二和(异或运算)时,表现出来的是码元异或,而波形相乘。
实验中用MATLAB 进行了仿真。
在进行GPS 系统的仿真中，P 码的生成较C/A 码复杂。
在实验过程中进行整体的仿真时没有考虑。
实际中GPS 的信道主要受到大尺度路径的损耗、电离层、对流层、多径效应、码环鉴相器输出时间t(单位：ms)鉴相器输出码NCO 频率控制字的改变时间t(单位：ms)频率控制字的改变值0100200300400500600700800900伪距P R (单位：m )-11L 1-m (单位：m )载波相位测量值多普勒效应等的影响。
电离层分布于地球大气层的顶部，该层中的传播速度与频率有关；对流层对信号传播路径有很大的影响；无线信道由于多径时延扩展，若收发机处于移动状态，各信号支路由于多普勒效应产生频域弥散，在实际中可采用Jace移动信道模型。
实验中只用了加性白噪声来模拟这些影响，与实际的情况存在一些差距。
本次实验采用了较为理想的情况下的系统收发模型，并对其进行了仿真得到了相关的波形和频谱。
GPS是以扩频通信为基础，主要应用扩频通信的抗干扰性能，通过不同的扩频码加以区分。
GPS卫星同时发射两路扩频信号，一路扩频信号将导航电文分别用C/A码和P码进行扩频，另一路扩频信号信号仅使用P码对导航电文进行扩频。
六、实验部分代码：1、GPS信号仿真部分：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)+reg2 a(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)+reg2 a(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';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,8,9,7]; % 生成一个周期的伪码序列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;%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));%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('调制后数据');2、GPS系统接收部分：function digCA=digitizg(svnum,chip_rate,fs) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CA=codegen1(svnum);ts=1/fs;index=ones(1,fs/1000);for n=1:fs/1000index(1,n)=ceil(ts*n/(1/chip_rate));index(end)=1023;enddigCA=CA(index);%产生2.8s的输入信号clear all;fs=5e6; % *** sampling freqts=1/fs; % *** sampling timen=fs/1000; % *** data pt in 1 mssvnum=input('enter satellite number= ');fl=1575.42e6;%RF frequencycodedelay=10; % code delay in code chip widthnn=2900*n; % *** total number of pts% nn=10*n;fc=1.25e6; % *** center freq without Doppler% fd=4812;% *** Doppler frequencyfd=0;d=fd/fl;phase0=0;chip_rate = 1.023e6*(1+d); %gold code clock rate in Hz.code2=codegen(svnum,nn,chip_rate,codedelay,fs);b=[0:nn-1];bpsk=code2.*sin(2*pi*(fc+fd)*b*ts+phase0);% %产生白噪声% noise_std=0;SNRdb=-10;noise_std=sqrt(0.5/(10^(SNRdb/10)));noise=noise_std*randn(1,length(bpsk));% data=bpsk+noise;%加白噪声data=bpsk;save data.mat datafunction[c1,c2]=calculatePLLCoef(beta,damp,K);dT=1e-3;wn=2*beta/(damp+1/(4*damp));%自然频率b1=(1/K)*((8*damp*wn*dT)/(4+4*damp*wn*dT+(wn*dT)^2)); b2=(1/K)*((4*(wn*dT)^2)/(4+4*damp*wn*dT+(wn*dT)^2));c1=b1+b2;c2=b1;clear all;f0=1.25e6;fl=1575.42e6;fc_lo=1.24e6;%IF freq minus max Doppler offsetfc_hi=1.26e6;%IF freq plus max Doppler offsetstep=1000;%size of one Doppler binfs=5e6; % *** sampling freqts=1/fs; % *** sampling timen=fs/1000; % *** data pt in 1 msT=2;%thresholdchip_rate=1.023e6;load data.mat datasvnum=input('enter satellite number= ');CA=digitizg(svnum,chip_rate,fs); % digitize C/A code%*********FFT ACQUISITION***********C=zeros(length(fc_lo:step:fc_hi),5000);i=0;clear j;% xf=fft(data1(1:n));for fc=fc_lo:step:fc_hii=i+1;t=(0:(n-1))/fs; %time vectorI_comp=cos(2*pi*fc*t).*data(1:n);Q_comp=sin(2*pi*fc*t).*data(1:n);X=fft(I_comp+j*Q_comp);F_ca=conj(fft(CA));C(i,:)=ifft(X.*F_ca);% lc=CA.*exp(j*2*pi*fc*t);% lcf=fft(lc);% C(i,:)=ifft(xf.*conj(lcf));end[amp crw]=max(max(abs(C'))); % find highest peak[amp ccn]=max(max(abs(C)));ph=ccn;fr=f0+(crw-11)*step;subplot(2,1,1)plot(abs(C(crw,1:5000)),'b')title(['GPS = ' num2str(svnum) ' ' 'max at' num2str(ph)])xlabel('码采样点','FontSize',9.0)ylabel('相关峰','FontSize',9.0)legend('基于FFT的相关值')subplot(2,1,2)plot(abs(C(:,ccn)),'*')title(['GPS = ' num2str(svnum) ' ' 'Freq =' num2str(fr) 'Hz'])xlabel('频率搜索点','FontSize',9.0)ylabel('相关峰','FontSize',9.0)legend('基于FFT的相关值',2)formatphformat long efrfiguremesh(abs(C))%%%%%%%%%%%%%%**********pull in**********8doppler_index=0;D=zeros(11,1); t=(0:(n-1))/fs; %time vectorfor fc=(fr-500):100:(fr+500) %search+/-250hz with 50hz increments doppler_index=doppler_index+1;for m=1:3Icomp=sin(2*pi*fc*t).*data((ph:5000+ph-1)+5000*(m-1));Qcomp=cos(2*pi*fc*t).*data((ph:5000+ph-1)+5000*(m-1));I=sum(Icomp.*CA);Q=sum(Qcomp.*CA);D(doppler_index)=D(doppler_index)+I.^2+Q.^2;endendD=D/max(D);a=find(D==max(D));f=fr+(a-6)*100;figureplot(D,'x')title(['GPS = ' num2str(svnum) 'carrier pullin at ' 'Freq =' num2str(f) 'Hz']) xlabel('频率牵引搜索点','FontSize',9.0)ylabel('相关峰的平方（归一化）','FontSize',9.0)legend('3次累加',2)format long eFclear all;close all;clc;tic;%%%%%%%%%%%%%%%%%%%initialize*********************load data.mat datasvnum=1;fl=1575.42e6;c=299792458.0;ph=50;%捕获得到的码相位data=data(ph:5000*2800+ph-1);a=pow2(32)/5e6;g=pow2(32);b=length(data);%载波NCO和码NCO均采用32位的% NCOFCWinit=a*1.023e6*(1+4800/fl);%码NCO初始频率控制字NCOFCWinit=a*1.023e6*(1+0/fl)NCOacc=zeros(1,5001);codeclk=zeros(1,5001);%码钟Qcomp_lo=zeros(1,5000); Icomp_lo=zeros(1,5000);M=a*1.25e6;% CNCOFCWinit=a*4800;%载波NCO初始频率控制字CNCOFCWinit=a*0;CNCOacc=zeros(1,5001);NCOFCW=NCOFCWinit;CNCOFCW=CNCOFCWinit;y(1)=0;x(1)=0;j=0;code=zeros(1,5001);prompt=zeros(1,5000);early=zeros(1,5000); late=zeros(1,5000);ca=codegen1(svnum);inda=1;% pll filter cofficients载波环环路滤波器系数[p1,p2]=calculatePLLCoef(50,0.707,(2*pi)/a);%code loop filter coefficients码环环路滤波器系数[c1,c2]=calculatePLLCoef(1,0.707,1/(2*a));% [c1,c2]=calculatePLLCoef(1,0.707,50);for m=1:2800%code tracking%产生1ms的本地码和本地载波for i=1:5000% i=(m-1)*5000+n;NCOacc(i+1)=NCOacc(i)+NCOFCW;%%% if (NCOacc(i+1)>hex2dec('FFFFFFFF'))if (NCOacc(i+1)>4294967295)%%% NCOacc(i+1) = NCOacc(i+1) - hex2dec('100000000');NCOacc(i+1) = NCOacc(i+1) -4294967296 ;end%%% if (NCOacc(i+1)>=hex2dec('80000000'))if (NCOacc(i+1)>=2147483648)codeclk(i+1) = 0;elsecodeclk(i+1) = 1;endif ((codeclk(i+1)-codeclk(i))>0) %码钟上升沿j=j+1;if j>1023j=j-1023;endcode(i+1)=ca(j);%本地码生成elsecode(i+1)=code(i);end%generate I/Q(pnco=phase offset)Icomp_lo(i)=sin(2*pi*(CNCOacc(i)+i*M)/g);%产生载波Qcomp_lo(i)=cos(2*pi*(CNCOacc(i)+i*M)/g);CNCOacc(i+1)=CNCOacc(i)+CNCOFCW;NCOacc(1)=NCOacc(5001);CNCOacc(1)=CNCOacc(5001);code(1)=code(5001);codeclk(1)=codeclk(5001);enddata1=data((m-1)*5000+1:m*5000);%read data%*****generate C/A codesprompt=code(2:5001);%即时码early=[prompt(4999:5000) prompt(1:4998)];%chip_delay=2;超前码late=[prompt(3:5000) prompt(1:2)];%滞后码Icomp=data1.*Icomp_lo;Qcomp=data1.*Qcomp_lo;%In_phaseIE(m)=sum(early.*Icomp)/5000;% earlyIP(m)=sum(prompt.*Icomp)/5000;%promptIL(m)=sum(late.*Icomp)/5000;%late%QuadratureQE(m)=sum(early.*Qcomp)/5000;%earlyQP(m)=sum(prompt.*Qcomp)/5000;%promptQL(m)=sum(late.*Qcomp)/5000;%lateEE(m)=IE(m).^2+QE(m).^2;%earlyPP(m)=IP(m).^2+QP(m).^2;%promptLL(m)=IL(m).^2+QL(m).^2;%late%costas loop for carrier trackingd2(m)=atan(QP(m)./IP(m));%atan discriminator% if m==30% [p1,p2]=calculatePLLCoef(10,0.707,2*pi/a);% % [p1,p2]=calculatePLLCoef(20,0.707,400*pi);% end%loop filterif m<2y(m)=0;else%filtered discriminator outputy(m)=y(m-1)+p1*d2(m)-p2*d2(m-1);endCNCOFCW=CNCOFCWinit+y(m);E_L(m)=(LL(m)-EE(m));%超前减去滞后包络if m<2x(m)=0;elsex(m)=x(m-1)+c1*E_L(m)-c2*E_L(m-1);endif m>=2000PR1(inda)=NCOacc(i+1)*c/(g*1.023e6);L1_m(inda)=CNCOacc(i)*c/(fl*g);inda=inda+1;end%调整码环NCO的频率NCOFCW=NCOFCWinit+x(m);endpr1=PR1(240:550);L1=L1_m(240:550);C1sm=pr1;wmin=0.0;wvar=0.005;w=1-wvar;C2sm=pr1;%begin observation time loopfor b=2:length(pr1)w1=1/b;C1sm(b)=w1*pr1(b)+(1-w1)*(C1sm(b-1)+L1(b)-L1(b-1));C2sm(b)=w*pr1(b)+(1-w)*(C2sm(b-1)+L1(b)-L1(b-1));w=w-wvar;if w<wminw=wmin;endendt1=[1:311];y3=(0.91428)*t1+6.6091;y4=(0.91311)*t1+6.5899;y1=(0.91282+8.6066e-007)*t1+6.6159-2.824e-005;y2=(0.91569+4.2217e-006)*t1+2049.3+0.029 809;subplot(2,1,1)plot(d2),grid ontitle(['载波环鉴相器输出'])xlabel('时间t(单位：ms)','FontSize',9.0)ylabel('鉴相器输出','FontSize',9.0)subplot(2,1,2)plot(y),grid ontitle(['载波NCO频率控制字的改变'])xlabel('时间t(单位：ms)','FontSize',9.0)ylabel('频率控制字的改变值','FontSize',9.0)figure(2)subplot(2,1,1)plot(E_L),grid ontitle(['码环鉴相器输出'])xlabel('时间t(单位：ms)','FontSize',9.0) ylabel('鉴相器输出','FontSize',9.0) subplot(2,1,2)plot(x),grid ontitle(['码NCO频率控制字的改变'])xlabel('时间t(单位：ms)','FontSize',9.0) ylabel('频率控制字的改变值','FontSize',9.0) figure(3)subplot(2,1,1)plot(PR1);hold ontitle(['伪距'])ylabel('PR(单位：m)','FontSize',9.0) subplot(2,1,2)plot(L1_m),hold onylabel('L1-m(单位：m)','FontSize',9.0)title('载波相位测量值')figure(4)subplot(2,1,1)plot(C1sm),title(['平滑伪距，w=1/k'])ylabel('PR(单位：m)','FontSize',9.0) subplot(2,1,2)plot(C1sm-y3)title(['残差,w=1/k'])ylabel('(单位：m)','FontSize',9.0)text(240,-0.08,'标准差：0.02818m')。