《GPS系统仿真》doc版.docx

GPS系统仿真实验报告姓名：陈珂学号：3150404002指导教师：高瑜翔1、 实验目的：分析扩频码的构成原理，基于MATLAB产生GPS信号，并进行调制传输，接收解调，最终建立完整的GPS系统并对接收的结果进行分析。2、 实验原理：扩频通信：扩频通信的基本理论依据是信息论中的香农公式其中,C是信道容量(bit/s),砰是信道带宽,S是信号的平均功率,N是噪声功率,香农公式表明了在给定信噪比的情况下,只要采用某种编码系统就能以任意小的差错概率,以接近于C的传输速率来传送信息。扩频序列采用的是伪随机序列,伪随机序列具有尖锐的自相关性,近似为零的互相关性这一特点。基于这一特性,可以通过检测本地伪随机序列与原伪随机序列的互相关峰值的大小来捕获伪随机序列。伪随机序列包括m序列、Gold序列、M序列和组合序列等。GPS信号中所使用的是一种组合Gold序列(C/A码/P码)。3、 实验内容：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级反馈移位寄存器组合产生,可以表示为其中为码元对应的时间：1/1023ms，为相位偏置的码元数。C/A码的产生如下图所示：3、 P码：P码是由两组各有两个12级反馈移位寄存器结合产生的,而每一组反馈移位寄存器又是由两个12级的m序列发生器组合而成的。因此,P码仿真模型需用4个m序列发生模块PN sequence Generator产生,PN1和PN2组成子码发生器;PN3和PN4组成子码发生器。又因为C/A码与P码采用的均是反像编码方式，所以同一组内及两组间的m序列模二和(异或运算)时,表现出来的也是码元异或,而波形相乘后接一个移位寄存器使P码延时一个基码，仿真电路如下图所示。P码仿真频率设为l0.23MHz,仿真时间设为0.1ms，子码系数为1100101000001,的系数为1111100100111,的系数为1111110111011,的系数为1001100011101。4、 系统框架及实现：GPS信号是数据码经扩频码.载波两次变频调制而成的,由卫星负责信号的发射;通过有噪信道的传输进入接收端;接收端负责信号的接收、捕获和跟踪,GPS仿真系统结构图如图所示。系统基本框架5、 GPS信号的捕获与跟踪：信号捕获的目的是为了粗略的确定C/A码在数据块中的起始位置、载波频率和码相位。信号跟踪的目的是为了获得精确的载波频率和码相位,并且当载波频率和码相位出现一定变化时,能够使本地载波频率和本地C/A码相位跟踪上这种变化。信号跟踪分为载波频率/相位跟踪和码相位跟踪。捕获和跟踪的实质都是一种对发射信号参数的估计。捕获是在不具备接收信号参数准确值的信息情况下,对接收载波频率和发送码字及相位的粗略估计,跟踪的实质也是一种对发送信号参数的估计,是在具备对接收载波频率和发送码字相位的粗略估计的信息的前提下,使用更好的估计手段,对发送信号的载波频率/相位和发送码相位的更精确估计。如果接收到的信号和本地复现信号的相位完全对齐,那么在相关过程中将会出现最大的峰值;但如果相位不完全对齐,只要偏移在一个码片内,仍然可以得到一个较大的峰值,所以捕获过程只能给出载波频率和码相位的一个粗略估计。捕获/跟踪过程中最重要的是实现两个同步:1.码元同步2.载波同步。在跟踪环节,码元同步负责信号解扩,载波同步负责信号解调。此外,由于是软件仿真,两者还肩负着信号的捕获工作。接收信号进行的处理框图和C/A码的捕获跟踪实现如下图所示：接收后进行的函数处理C/A码的捕获和跟踪实现框图码元同步：当捕获到卫星信号后,本地码的频率和相位与接收到的码元存在较大的误差,因此必须进行调整,以达到本地码与卫星信号C/A码同频同相，这个调整过程由码跟踪环来完成,称之为相关接收或解扩。4、 实验结果：用m文件实现GPS系统的信号发送以及接收部分的仿真，产生两种伪随机序列后分别对两路信号进行调制，然后对C/A码进行捕获和跟踪，对得到的子帧进行识别再计算星历数据和伪距以及卫星的位置和用户位置。得到的波形如下图所示,其中图2是svn=10的C/A码序列在码相位和载波频率二位平面上搜索的相关峰值：图2 svn=10的C/A码序列在码相位和载波频率二位平面上搜索的相关峰值5、 实验总结： 本次实验首先对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的信道主要受到大尺度路径的损耗、电离层、对流层、多径效应、多普勒效应等的影响。电离层分布于地球大气层的顶部，该层中的传播速度与频率有关；对流层对信号传播路径有很大的影响；无线信道由于多径时延扩展，若收发机处于移动状态，各信号支路由于多普勒效应产生频域弥散，在实际中可采用Jace移动信道模型。实验中只用了加性白噪声来模拟这些影响，与实际的情况存在一些差距。本次实验采用了较为理想的情况下的系统收发模型，并对其进行了仿真得到了相关的波形和频谱。GPS是以扩频通信为基础，主要应用扩频通信的抗干扰性能，通过不同的扩频码加以区分。GPS卫星同时发射两路扩频信号，一路扩频信号将导航电文分别用C/A码和P码进行扩频，另一路扩频信号信号仅使用P码对导航电文进行扩频。6、 实验部分代码：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; else z1b=mod(N-C1*y1a),4093); x1bWork=1; x1b=(N-z1b-C1*y1a)/4093; end m=mod(N,(C1+37); y2a=(N-m)/(C1+37); if (m=C1) dv=m-C1; else dv=0; end z2a=mod(m-dv),4092); x2a=mod(m-dv)-z2a)/4092),3750); z2b=mod(m-dv),4093); if (m=C2) x2b=3748; else x2b=(m-z2b)/4093; end %各移位寄存器的状态 for i=1:z1a slave1a=mod(reg1a(6)+reg1a(8)+reg1a(11)+reg1a(12),2); reg1a(2:12)= reg1a(1:11); reg1a(1)=slave1a; end for i=1:z1b slave1b=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; end for i=1:z2a slave2a=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; end for i=1:z2b slave2b=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=4091 rx1a=1; end if z1b=4092 rx1b=1; end if z2a=4091 rx2a=1; x2aWork=0; end if z2b=4092 rx2b=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=1 if rx1b=0 slave1b=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=1 reg1b=0 0 1 0 1 0 1 0 1 0 1 0; rx1b=0; end end else if x1bWork=0 end end if reg1b=0 1 0 1 0 1 0 1 0 1 0 0 rx1b=1; x1bcou=x1bcou+1; if x1bcou=3749 x1bwork=0; x1bcou=0; end end %寄存器x1a的移位函数 if rx1a=0 slave1a=mod(reg1a(6)+reg1a(8)+reg1a(11)+reg1a(12),2); reg1a(2:12)=reg1a(1:11); reg1a(1)=slave1a; else if rx1a=1 reg1a=0 0 0 1 0 0 1 0 0 1 0 0; rx1a=0; end end if reg1a=0 0 1 0 0 1 0 0 1 0 0 0 rx1a=1; x1acou=x1acou+1; if x1acou=3750 x1bwork=1; x1acou=0; end end %寄存器x2b的移位函数 if x2bWork=1 if rx2b=0 slave2b=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=1 x2bout=reg2b(11); reg2b=0 0 1 0 1 0 1 0 1 0 1 0; rx2b=0; end end else if x2bWork=0 reg2b=0 0 1 0 1 0 1 0 1 0 1 0; x2bWork=1; rx2b=0; end end if reg2b=0 1 0 1 0 1 0 1 0 1 0 0 rx2b=1; x2bcou=x2bcou+1; if x2bcou=3749 x2bWork=0; x2bcou=0; end end %寄存器x2a的移位函数 if x2aWork=1 if rx2a=0 slave2a=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=1 reg2a=1 0 1 0 0 1 0 0 1 0 0 1; rx2a=0; end end else if x2aWork=0if rx2a=1 cou37=cou37+1; if cou37=37 rx2a=0; x2awork=1; cou37=0; end end end end if reg2a=0 1 0 0 1 0 0 1 0 0 1 1 rx2a=1; x2acou=x2acou+1; if x2acou=3750 x2awork=0; x2acou=0; end end endfor i=1:NumberPCode p(i)= mod( x1(i)+x2(i),2); endp=p;pcode=p;C/A码产生的子程序CAgenerate：function cacode=CAgenerate(number) if (number37) disp(输入参数必须在1 37之间取值); return endCACode=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:1023 CACode(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:c if y(i)0.5 datacode(i)=0; else datacode(i)=1; endendx(1)=0;show(1)=datacode(1);p=2;for i=1:c for j=1:100 x(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:100 x(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:c if datacode(i)=1 datacodek(i-1)*50+1:i*50)=ones(1,50); else datacodek(i-1)*50+1:i*50)=zeros(1,50); endendfor i=1:c for j=1:50 addr=rem(i-1)*50+j),10); if addr=0 addr=10; end kuopindata(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:100 x(p)=i-1+j*0.01; show(p)=kuopindata(i); p=p+1; endend subplot(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:c if y(i)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; else codeclk(i+1) = 1; end if (codeclk(i+1)-codeclk(i)0) %码钟上升沿 j=j+1; if j1023 j=j-1023; end code(i+1)=ca(j);%本地码生成 else code(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); end data1=data(m-1)*5000+1:m*5000);%read data %*generate C/A codes prompt=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_phase IE(m)=sum(early.*Icomp)/5000;% early IP(m)=sum(prompt.*Icomp)/5000;%prompt IL(m)=sum(late.*Icomp)/5000;%late %Quadrature QE(m)=sum(early.*Qcomp)/5000;%early QP(m)=sum(prompt.*Qcomp)/5000;%prompt QL(m)=sum(late.*Qcomp)/5000;%late EE(m)=IE(m).2+QE(m).2;%early PP(m)=IP(m).2+QP(m).2;%prompt LL(m)=IL(m).2+QL(m).2;%late %costas loop for carrier tracking d2(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 filter if m2 y(m)=0; else %filtered discriminator output y