gps信号和干扰的matlab仿真代码

x1=10.23y1=6.24z1=17.8 %假设的卫星位置参数

Code1=fGenerateNavigationData(x1,y1,z1)%产生导航电文,通过调用子函数fGenerateNavigationData实现

%将产生的Code1中的0找出转换为-1

index1=find(Code1==0)

Code1(index1)=-ones(1,length(index1))

SvNum=12 %设定卫星编号为12

Code2=zeros(1,1) %定义Code2的初值为0

Temp=fGenerateCAcode3(SvNum)%将编号为SvNum的卫星通过调用子函数fGenerateCAcode3生成C/A码

%将段键Temp中的0找出并转换为-1

index1=find(Temp==0)

Temp(index1)=-ones(1,length(index1))

Temp=[Temp Temp]

%生成Code2

for i=1:length(Code1)

Code2=[Code2 Code1(1,i)*Temp]

end

Code2=Code2(2:length(Code2))

%每位数据通过正弦波来调制

SinWave=sin([0:2*pi/8:2*pi*3/8])

SinWave=single(SinWave)

GPSsignals=zeros(1,1)

SinWave=[SinWave SinWave]

for i=1:length(Code2)

GPSsignals=[GPSsignals Code2(1,i)*SinWave]

end

whos

GPSsignals=GPSsignals(2:length(GPSsignals))

figure(2)

plot(GPSsignals)axis([100,4000,-1.5 1.5])

figure(3)

plot(10*log10(abs(fft(GPSsignals))))

产生C/A码的方法

clc

clear all

%产生C/A码的方法一

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%,%%%%%%%

k1=2k2=6delay=5 %定义参数k1、k2以及延迟

Reg=-ones(1,10) %定义寄存器1～10的初始值都为-1

%通过循环产生长度为1024的MLS序列以及G2序列

for j=1:1023

MLS(j)=Reg(10) %将士0号寄存器的输出作为MLS输出

modulo=Reg(2)*Reg(3)*Reg(6)*Reg(8)*Reg(9)*Reg(10)

Reg(2:10)=Reg(1:9)

Reg(1)=modulo

g2(j)=Reg(k1)*Reg(k2) %将参数k1和k2所代表的寄存器模2相加后作为G2输出

%g2即为输出

end

%将G2与MLS进行延迟检验

if MLS == g2([delay:1023 1:delay-1])

disp('OK')

else

disp('指族not match')

end

%在G2序列中找出-1并转换为0,找出1并转换为1

ind1=find(g2==1)

ind2=find(g2==1)

g2(ind1)=ones(1,length(ind1))

g2(ind2)=zeros(1,length(ind2))

temp=g2(1:120)

x(1)=0

Show(1)=temp(1)

P=2

%下面的循环是为了将结果显示成握逗巧方波形式

for i=2:length(temp)

if temp(i)==temp(i-1)

x(P)=i-1

Show(P)=temp(i-1)

x(P+1)=i-1+0.01

Show(P+1)=temp(i)

P=P+2

else

Show(P)=temp(i)

x(P)=i

P=P+1

end

end

%画出仿真结果图

plot(x,Show+1)

axis([0 length(x)-60 -0.1 1.1])

grid on

子程序：

%子函数fGenerateNavigationData

function y=fGenerateNavigationData(x1,y1,z1)

%将传进的参数转换为十六进制数

%x=x1

%y=y1

%z=z1

x=10.23

y=6.24

z=17.8

str1=num2hex(x)

str2=num2hex(y)

str3=num2hex(z)

Table1=[0 0 0 0

0 0 0 1

0 0 1 0

0 0 1 1

0 1 0 0

0 1 0 1

0 1 1 0

0 1 1 1

1 0 0 0

1 0 0 1

1 0 1 0

1 0 1 1

1 1 0 0

1 1 0 1

1 1 1 0

1 1 1 1]

TotalStr=[str1 str2 str3]

DataCode=zeros(1,1)

l=length(TotalStr)

%将TotalStr中的数转换为ASCⅡ码表中的数值

for i=1:l

temp=int8(TotalStr(i))

if (temp>58)

temp=temp-96+10

else

temp=temp-47

end

DataCode=[DataCode Table1(temp,:)]

end

y=DataCode(2:length(DataCode))

figure(1)

subplot(2,1,1)

plot(y)axis([0 195 -0.5 1.5])title('导航电文数据（0 1）')

%将产生的Code1中的0找出转换为-1

index1=find(y==0)

y(index1)=-ones(1,length(index1))

subplot(2,1,2)plot(y)axis([0 195 -1.5 1.5])title('导航电文数据（-1 1）')

自己慢慢理解吧

pudn上面有相关资料可以下载，可以参考一下

具体GPS信号跟踪误差怎么求不知道，如果你能求得该误差信号， 具体的包络算法解释如下

求包络(envelop)的最方便的方法是隐脊用希灶伍渗尔伯特变换， 希尔伯特变换返回的是analytical signal。

对于一个任意的橘锋时序信号x(t)，一个简单的程序就是：

y(t)=hilbert(x(t))             %%%return analytical signal

envelop=abs(y(t))          %%%get the magnitude of the analytical signal,namely envelop

我自己的以前用的一个图附上了， 红色部分是原始信号，绿色部分是包络部分，也就是信号的幅值。 包络线正好将信号包住。

在下面的仿真的代码中，理想的观测量不是真实数据，而是自生成的正弦波数据，在真实的应用场景中，应该是一系列的参考数键凳没据。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 卡尔曼滤波器在INS-GPS组合导航中应用仿真

% Author : lylogn

% Email : lylogn@gmail.com

% Company: BUAA-Dep3

% Time : 2013.01.06

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 参考文献:

% [1]. 邓正隆. 惯导技术, 哈尔滨工业大学粗仔出版社.2006.

clear all

%% 惯性-GPS组合导航模型参数初始化

we = 360/24/60/60*pi/180%地球自转角速度，弧度/s

psi = 10*pi/180 %psi角度 / 弧度

Tge = 0.12

Tgn = 0.10

Tgz = 0.10 %这三个参数的含义详见参考文献

sigma_ge=1

sigma_gn=1

sigma_gz=1

%% 连续空间系统状态方程

% X_dot(t) = A(t)*X(t) + B(t)*W(t)

A=[0 we*sin(psi) -we*cos(psi) 1 0 0 1 0 0

-we*sin(psi) 0 0 0 1 0 0 1 0

we*cos(psi) 0 0 0 0 1 0 0 1

00 0 -1/Tge 0 0 0 0 0

00 0 0 -1/Tgn 0 0 0 0

00 0 0 0 -1/Tgz 0 0 0

00 0 0 0 0 0 0 0

00 0 0 0 0 0 0 0

00 0 0 0 0 0 0 0]%状态转移矩阵

B=[00 0sigma_ge*sqrt(2/Tge) 0 0 0 0 0

00 00 sigma_gn*sqrt(2/Tgn) 0 0 0 0

00 00 0 sigma_gz*sqrt(2/Tgz) 0 0 0]'%输入控制矩阵

%% 转化为离散时间系统状态方程

% X(k+1) = F*X(k) + G*W(k)

T = 0.1

[F,G]=c2d(A,B,T)

H=[10 0 0 0 0 0 0 0

0 -sec(psi) 0 0 0 0 0 0 0]%观测矩稿纳阵

%% 卡尔曼滤波器参数初始化

t=0:T:50-T

length=size(t,2)

y=zeros(2,length)

Q=0.5^2*eye(3)%系统噪声协方差

R=0.25^2*eye(2) %测量噪声协方差

y(1,:)=2*sin(pi*t*0.5)

y(2,:)=2*cos(pi*t*0.5)

Z=y+sqrt(R)*randn(2,length) %生成的含有噪声的假定观测值，2维

X=zeros(9,length) %状态估计值，9维

X(:,1)=[0,0,0,0,0,0,0,0,0]' %状态估计初始值设定

P=eye(9) %状态估计协方差

%% 卡尔曼滤波算法迭代过程

for n=2:length

X(:,n)=F*X(:,n-1)

P=F*P*F'+ G*Q*G'

Kg=P*H'/(H*P*H'+R)

X(:,n)=X(:,n)+Kg*(Z(:,n)-H*X(:,n))

P=(eye(9,9)-Kg*H)*P

end

%% 绘图代码

figure(1)

plot(y(1,:))

hold on

plot(y(2,:))

hold off

title('理想的观测量')

figure(2)

plot(Z(1,:))

hold on

plot(Z(2,:))

hold off

title('带有噪声的观测量')

figure(3)

plot(X(1,:))

hold on

plot(X(2,:))

hold off

title('滤波后的观测量')

---