1、选择公交车辆车牌号码进行实时监控。
2、通过选择开始时间和结束时间进行车辆的轨迹回放。
3、根据开始个结束时间打印报表。
4、针对哈尔滨公交公司内部人员使用。
在下面的仿真的代码中,理想的观测量不是真实数据,而是自生成的正弦波数据,在真实的应用场景中,应该是一系列的参考数据。%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 卡尔曼滤波器在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('滤波后的观测量')
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上面有相关资料可以下载,可以参考一下
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)