作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,matlab项目合作可私信。
个人主页:Matlab科研工作室
个人信条:格物致知。
更多Matlab仿真内容点击
由于目标信号的发射时间未知,无源定位技术大多利用TDOA(到达时间差)进行目标定位.本文将求解GPS单点定位的Bancroft算法
function [eX,eb] = algebraicGPSequations(satPos, tVals)
% Calculates the position of a receiver given the time delay of arrival to
% at least 4 satellites in 3D (or at least 3 satellites in 2D).
% INPUTS: satPos: is Cartesian location of satellites (in 1D, 2D, or 3D).
% Each row is the position of one satellite.
% tVals: is the arrival time of the signal = distance from
% satellite to reciever + offset b.
% OUTPUTS: eX: estimated position of the user.
% eb: estimate of the offset time b.
%
% Implements the solution of Stephen Bancroft, "An Algebraic Solution of the GPS
% Equations," in IEEE Transactions on Aerospace and Electronic Systems,
% vol. AES-21, no. 1, pp. 56-59, Jan. 1985, doi: 10.1109/TAES.1985.310538.
% https://ieeexplore.ieee.org/document/4104017
%
% pdf at
% https://cas.tudelft.nl/Education/courses/ee4c03/assignments/indoor_loc/Bancroft85loc.pdf
% Aaron T. Becker, University of Houston
% Version 1.0.2: Feb 23, 2021 initial submission
% Version 1.0.3: Feb 05, 2023 updated by adjusting the documentation
%% Initialization
if nargin < 1
disp('No arguments provided, so this generates an example problem.')
% Default values: The node (satellite) positions are randomized. The transmitter
% is randomly positioned outside of the area covered by the nodes.
nSat = 4; % number of nodes (satellites)
dim = 3; %TODO: fails with dim = 1. Not sure why. Works for 2 and 3 dimensions
satPos = 50*rand(nSat,dim)-25; % position of nodes
if dim == 2
th = 2*pi*rand; %random angle
Xpos = (50+(50*rand)).*[cos(th) sin(th)];%user location (unknown to satellites)
else
Xpos = 150*rand(1,dim)-75;
end
b = 15+15*rand; % timing offset (unknown to satellites)
noiseLvl = 0; % Ideal signal + Noise.
tVals = sqrt(sum((Xpos-satPos).^2,2))+b + noiseLvl*rand(nSat,1); % Eqn (2)
fprintf("Running algebraicGPSequations with nSat = %d, in %dD.\n Randomly generated satellite positions.\n Actual user position is:",nSat,dim );
disp(Xpos)
else
nSat = length(tVals);
b = NaN; % unknown to satellites
Xpos = NaN(1,size(satPos,2)); % unknown to satellites
end
if nSat < size(satPos,2) +1
disp(['Error: system is underdetermined with ',num2str(nSat),' satellites and ', num2str(size(satPos,2)), ' dimensions'])
end
%% calculation of user location and time offset
A = [satPos, tVals]; % Eqn (5)
i0 = ones(nSat,1); % Eqn (6)
r = ones(nSat,1); % Eqn (7)
for i = 1:nSat
r(i) = minkowskiSum(A(i,:),A(i,:))/2; % Eqn (8)
end
B = pinv(A); % same as inv(A'*A)*A'; % Eqn (9)
u = B*i0; % Eqn (10)
v = B*r; % Eqn (11)
E = minkowskiSum(u,u); % Eqn (12)
F = minkowskiSum(u,v) - 1; % Eqn (13)
G = minkowskiSum(v,v); % Eqn (14)
lam1 = (-2*F+sqrt((2*F)^2-4*E*G))/(2*E); % Eqn (15) +solution to quadratic
lam2 = (-2*F-sqrt((2*F)^2-4*E*G))/(2*E); % Eqn (15) -solution to quadratic
y1 = lam1*u+v; % Eqn (16) +
y2 = lam2*u+v; % Eqn (16) -
Tx1 = y1(1:end-1)'; % Eqn (17)
Tx2 = y2(1:end-1)'; % Eqn (17)
b1 = -y1(end); % Eqn (17)
b2 = -y2(end); % Eqn (17)
%determine which position answer is right by calculating the average error
tvals1 = sqrt(sum((Tx1-satPos).^2,2))+b1;
tvals2 = sqrt(sum((Tx2-satPos).^2,2))+b2;
% e1 = mean(abs(tvals1-tVals));
% e2 = mean(abs(tvals2-tVals));
e1 = std(tvals1-tVals);
e2 = std(tvals2-tVals);
if e1 < e2
legent = {'Satellite positions','X position','estimate 1 (correct)','estimate 2'};
eX = Tx1;
eb = b1;
else
legent = {'Satellite positions','X position','estimate 1','estimate 2 (correct)'};
eX = Tx2;
eb = b2;
end
%% Plot the process
cplot = @(r,x0,y0,linet) plot(x0 + r*cos(linspace(0,2*pi)),y0 + r*sin(linspace(0,2*pi)),linet);
if size(satPos,2) == 1
figure(1); clf; hold on;
p(1) = plot(satPos,zeros(size(satPos)),'bx');
hold on
for i = 1:length(satPos)
cplot(tVals(i)-b,satPos(i),0,'g-');
end
p(2) = plot(Xpos(1),0, 'xg');
set(p(2),'color', [0,0.5,0])
xlabel('x');ylabel('y')
axis equal
if abs(max(tVals)-min(tVals) - (max(satPos)-min(satPos))) < 1e-10
[minV,minInd] = min(satPos);
[maxV,maxInd] = max(satPos);
legent = {'Satellite positions','X position'};
if tVals(maxInd) < tVals(minInd)
str = sprintf("The receiver position is unknown, but > %.2f",maxV);
a = axis;
reg = patch([max(satPos),max(satPos),a(2),a(2)],[a(4),a(3),a(3),a(4)],[0.5,0.5,1]);
uistack(reg,'bottom')
else
str = sprintf("The receiver position is unknown, but < %.2f",minV);
a = axis;
reg = patch([min(satPos),min(satPos),a(1),a(1)],[a(4),a(3),a(3),a(4)],[0.5,0.5,1]);
uistack(reg,'bottom')
end
title(sprintf('Act X = [%.2f], b=%.2f \n %s',Xpos(1),b,str))
else
p(3) = plot(Tx1(1),0, 'ob');
p(4) = plot(Tx2(1),0, 'or');
for i = 1:length(satPos)
cplot(tVals(i)-b,satPos(i),0,'g-');
cplot(tVals(i)-b1,satPos(i),0,'b--');
cplot(tVals(i)-b2,satPos(i),0,'r--');
end
title(sprintf('Act X = [%.2f], b=%.2f \n Est X = [%.2f], b=%.2f',Xpos(1),b,eX(1),eb))
end
legend(p,legent);
axis tight
disp(['Act X = [',num2str(Xpos),'], b=',num2str(b)])
disp(['Est X = [',num2str(eX), '], b=',num2str(eb)])
disp('std(error), b, Tx for two solutions:')
disp([e1,b1,Tx1])
disp([e2,b2,Tx2])
elseif size(satPos,2) == 2
figure(1); clf; hold on;
p(1) = plot(satPos(:,1),satPos(:,2),'bx');
hold on
for i = 1:length(satPos)
cplot(tVals(i)-b,satPos(i,1),satPos(i,2),'g-');
cplot(tVals(i)-b1,satPos(i,1),satPos(i,2),'b--');
cplot(tVals(i)-b2,satPos(i,1),satPos(i,2),'r--');
end
p(2) = plot(Xpos(1),Xpos(2), 'x');
set(p(2),'color', [0,0.5,0])
p(3) = plot(Tx1(1),Tx1(2), 'ob');
p(4) = plot(Tx2(1),Tx2(2), 'or');
legend(p,legent);
xlabel('x');ylabel('y')
title(sprintf('Act X = [%.2f,%.2f], b=%.2f \n Est X = [%.2f,%.2f], b=%.2f',Xpos(1),Xpos(2),b,eX(1),eX(2),eb))
axis equal
elseif size(satPos,2) == 3
figure(1); clf; hold on;
p(1) = plot3(satPos(:,1),satPos(:,2),satPos(:,3),'bx');
hold on
[X,Y,Z] = sphere;
lightGrey = 0.8*[1 1 1];
for i = 1:length(satPos)
radius = tVals(i)-b;
surf(satPos(i,1)+radius*X,satPos(i,2)+radius*Y,satPos(i,3)+radius*Z,'facecolor','g','facealpha',0.02,'EdgeColor',lightGrey,'edgealpha',0.2)
end
p(2) = plot3(Xpos(1),Xpos(2),Xpos(3), 'x');
set(p(2),'color', [0,0.5,0])
p(3) = plot3(Tx1(1),Tx1(2),Tx1(3), 'ob');
p(4) = plot3(Tx2(1),Tx2(2),Tx2(3), 'or');
legend(p,legent);
xlabel('x');ylabel('y');zlabel('z')
title(sprintf('Act X = [%.2f,%.2f,%.2f], b=%.2f \n Est X = [%.2f,%.2f,%.2f], b=%.2f',Xpos(1),Xpos(2),Xpos(3),b,eX(1),eX(2),eX(3),eb))
axis equal
view([1 1 0.75]) % adjust the viewing angle
elseif size(satPos,2) > 2
disp(['plot is only supported for 1D & 2D, this is ',num2str(size(satPos,2)),'D'])
disp(['Act X = [',num2str(Xpos),'], b=',num2str(b)])
disp(['Est X = [',num2str(eX), '], b=',num2str(eb)])
disp('std(error), b, Tx for two solutions:')
disp([e1,b1,Tx1])
disp([e2,b2,Tx2])
end
function mSum = minkowskiSum(a,b)
mSum = sum(a(1:end-1).*b(1:end-1)) - a(end)*b(end); % Eqn (4)
end
end
[1]林云松, 孙卓振, and 彭良福. "基于Bancroft算法的多点定位TOA-LS估计." 电子学报 40.003(2012):621-624.
[2]林云松, 孙卓振, 彭良福. 基于Bancroft算法的多点定位TOA-LS估计[J]. 电子学报, 2012(03):207-210.
️部分理论引用网络文献,若有侵权联系博主删除
️ 关注我领取海量matlab电子书和数学建模资料
文章浏览阅读2.1k次。原文链接先说说编解码问题编码转换时,通常需要以unicode作为中间编码,即先将其他编码的字符串解码(decode)成unicode,再从unicode编码(encode)成另一种编码。 Eg:str1.decode('gb2312') #将gb2312编码的字符串转换成unicode编码str2.encode('gb2312') #将unicode编码..._python中encode在什么模块
文章浏览阅读949次,点赞21次,收藏15次。本文介绍了Java中的数据输入流(DataInputStream)和数据输出流(DataOutputStream)的使用方法。
文章浏览阅读111次。ie无法兼容_ie 浏览器 newdate
文章浏览阅读239次。这篇文章把 Docker 和 K8s 的关系给大家做了一个解答,希望还在迟疑自己现有的知识储备能不能直接学 K8s 的,赶紧行动起来,K8s 是典型的入门有点难,后面越用越香。
文章浏览阅读561次。ADI中文手册获取方法_adi 如何查看数据手册
文章浏览阅读1k次,点赞4次,收藏3次。React 获取接口数据实现分页效果以拼多多接口为例实现思路加载前 加载动画加载后 判断有内容的时候 无内容的时候用到的知识点1、动画效果(用在加载前,加载之后就隐藏或关闭,用开关效果即可)2、axios请求3、map渲染页面4、分页插件(antd)代码实现import React, { Component } from 'react';//引入axiosimport axios from 'axios';//引入antd插件import { Pagination }_react 分页
文章浏览阅读449次,点赞9次,收藏7次。这个变量与验签过程中的SignatureVerificationFilter::PUT_MESSAGE这个宏是对应的,SignatureVerificationFilter::PUT_MESSAGE,如果在签名过程中putMessage设置为true,则在验签过程中需要添加SignatureVerificationFilter::PUT_MESSAGE。项目中使用到了CryPtopp库进行RSA签名与验签,但是在使用过程中反复提示无效的数字签名。否则就会出现文章开头出现的数字签名无效。_cryptopp 签名
文章浏览阅读848次。新闻稿是新闻从业者经常使用的一种文体,它的格式与内容都有着一定的规范。本文将从新闻稿的格式和范文两个方面进行介绍,以帮助读者更好地了解新闻稿的写作_新闻稿时间应该放在什么位置
文章浏览阅读1.7k次。Java中的转换器设计模式 在这篇文章中,我们将讨论 Java / J2EE项目中最常用的 Converter Design Pattern。由于Java8 功能不仅提供了相应类型之间的通用双向转换方式,而且还提供了转换相同类型对象集合的常用方法,从而将样板代码减少到绝对最小值。我们使用Java8 功能编写了..._java转换器模式
文章浏览阅读150次。1,kubectl run创建pods[root@master ~]# kubectl run nginx-deploy --image=nginx:1.14-alpine --port=80 --replicas=1[root@master ~]# kubectl get podsNAME READY STATUS REST...
文章浏览阅读128次。PAT菜鸡进化史_乙级_1003“答案正确”是自动判题系统给出的最令人欢喜的回复。本题属于 PAT 的“答案正确”大派送 —— 只要读入的字符串满足下列条件,系统就输出“答案正确”,否则输出“答案错误”。得到“答案正确”的条件是: 1. 字符串中必须仅有 P、 A、 T这三种字符,不可以包含其它字符; 2. 任意形如 xPATx 的字符串都可以获得“答案正确”,其中 x 或者是空字符串,或..._1003 pat乙级 最优
文章浏览阅读5.6k次。CH340与Android串口通信为何要将CH340的ATD+Eclipse上的安卓工程移植到AndroidStudio移植的具体步骤CH340串口通信驱动函数通信过程中重难点还存在的问题为何要将CH340的ATD+Eclipse上的安卓工程移植到AndroidStudio为了在这个工程基础上进行改动,验证串口的数据和配置串口的参数,我首先在Eclipse上配置了安卓开发环境,注意在配置环境是..._340串口小板 安卓给安卓发指令