数字中频GPS信号的MATLAB仿真

【matlab国外编程代写】GPS信号捕获matlab仿真clear;%指标是signalpower=-185dBW，此时可以进行室内定位；%指标是SNR=-52dB，此时noisepower=-133dBWDopler_bound = 10; %原来是20，为加快运行速度，减少频率槽数Freq_bin = 500;Fc = 1250000;Length_of_data = 5000;Length_of_fft = 4096;NUM = 5;%50010;Signalpower = -177 %I路信号的功率（dBW）Noisepower = -143 %噪声的复功率（dBW）Threshold = 1e5;Times_of_ncoh = 40;%变量初始化successcount=0;%在卫星从地平线升起或降落时，用户速度900m/s时,径向速度差不多(+ -)10K，fid_ca=fopen('CA_5000.txt','r');CA_samp=fscanf(fid_ca,'%d');%5000点fclose(fid_ca);CA_samp=CA_samp';%本地伪码FFT:CA_fftCA_local=Average(Length_of_data, Length_of_fft, CA_samp);%plot(CA_local)CA_fft=fft(CA_local);%plot(imag(CA_fft))Idata=ReadIn( ); %把文本Idata全部读到Idata数组中Qdata=ReadIn( ); %Idata=Idata';Qdata=Qdata';Idata(1:NUM)=[]; %删去前NUM个元素Qdata(1:NUM)=[];AMP=sqrt(2*10^(Signalpower/10)) %输入量为幅度系数，单位相当于（V），%AMP=1e-9，此时signal power=-183dBW DIV_AWGN= 10^(Noisepower/10) %输入量为高斯白噪声的方差。

摘要扩频通信是近几年来迅速发展起来的一种通信技术。

在早期研究这种技术的主要目的是为提高军事通信的和抗干扰性能，因此这种技术的开发和应用一直是处于状态。

扩频技术在军事应用上的最成功例可以以美国和俄国的全球定位系统（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目录引言41GPS理论与其特性51.1GPS系统概述51.2GPS信号构成71.2.1M序列101.2.2C/A码131.2.3P码151.3小结192MATLAB软件192.1MATLAB软件简介202.2MATLAB应用概述212.2.1 MATLAB功能介绍212.2.2 MATLAB使用方法233 GPS卫星导航信号算法与其MATLAB仿真293.1 C/A码仿真代码与其仿真结果293.2 P码的仿真代码与其仿真结果323.3 结果分析与其相关性分析36结论37致38参考文献39附录A 英文原文40附录B 中文翻译41附录C C/A码源代码41附录D P码源代码54引言全球卫星定位系统，简称GPS系统，可在全球围，全天候为用户连续地提供高精度的位置、速度和时间信息。

GPS卫星轨道计算及其MATLAB仿真黎奇，白征东，李帅，陈波波（清华大学地球空间信息研究所，北京 100084）一、程序设计思路1. 读取RINEX文件（注意：文件路径）2. 计算测量日周积秒（测量日的格里历→GPST）3. 按卫星轨道计算步骤计算WGS-84坐标系坐标（内插）4. 按需要将WGS-84坐标系下坐标转换为所需坐标系坐标5. 画图输出二、n 文件说明及读取程序参考时刻oe t 的RINEX 格式的 “”广播星历文件具体如下：（加粗部分为本次轨道化Ω，率i ，弧度/秒4-22）标svacc ，米）收到的卫星信号解，秒）文件名：RinexNreader.m 输 入：文件地址，卫星编号三、计算测量日的周积秒文件名：GCtoGPS.m （其中调用函数：GCtoJD.m）输入：指定公历的年、月、日、时、分、秒文件名：GCtoJD.m输 入：指定公历的 年、月、日四、GPS 卫星轨道计算步骤及计算程序1. 计算卫星运动的平均角速度n平均角速度()03n =经摄动参数n ∆改正后的平均角速度0n n n =+∆3#61-79），n ∆（2#42-60）；14323.98600510/GM m s =⨯ 2. 计算归化时间k t说明：①广播星历是oe t 时刻的，对应的轨道参数也是oe t 时刻的，而观测时间在t 时刻，显然oe t t ＜。

所以，要想获得t 时刻的轨道参数，需要知道t 与oe t 之间的差值即k t 。

以此，按照oe t 时刻轨道参数，外推t 时刻轨道参数。

②k t 的起算时间是星期六/星期日子夜0点，当302400k t s ＞时，604800k t s -；当302400k t s ＜-时，+604800k t s 。

（604800s=1周） =k oe t t t -，且604800302400604800302400k k k k k k t t t t t t =-⎧⎨=+⎩＞ ＜-已知：oe t （1#4-21）3. 计算观测时刻的平近点角k M0k k M M nt =+已知：0M （2#61-79），n （见1），k t （见2） 4. 计算观测时刻的偏近点角k Esin k k k E M e E =+已知：k M （见3），e （3#23-41）方法：迭代解算，设初值0k k E M =，迭代2次基本收敛。

配套毕业设计论文见百度文库请搜索《基于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]);。

数字中频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接收机信号处理所需的数字信号。