一、仿真要求

要求一：编写两个函数TDOA_CHAN和TDOA_Taylor得到位置的估计。
要求二：用RMSE实现两种算法的性能比较, 得到两种算法的RMSE曲线对比图，横坐标为噪声方差，纵坐标为RMSE。

二、仿真方案的设计、仿真结果及结论

1.仿真方案的设计

TDOA基本定位过程如下

2.仿真结果

3.仿真结论

三、主要仿真代码

1.main函数中的代码如下：

%TDOA定位算法的仿真
clear all;
clc;
%定义四个参与基站的坐标位置
BS1=[0,0];BS2=[500,0];BS3=[500,500];BS4=[0,500];
%BS5=[600,500];
%移动台MS的初始估计位置
MS=[50,100];
std_var=[1e-2,5e-2,1e-1,5e-1,1];%范围
A=[BS1;BS2;BS3;BS4];%矩阵A包含4个初始坐标
number=10000;
for j=1:length(std_var)%循环
    error1=0;error2=0;%初始误差置为0
    std_var1=std_var(j);%令std_var1等于当前数组的值
    for i=1:number%多次循环
        r1=A-ones(4,1)*MS;
        r2=(sum(r1.^2,2)).^(1/2);%计算移动台到各个基站的实际距离
        r=r2(2:end,:)-ones(3,1)*r2(1,:)+std_var1*randn(3,1); %表示从[2,i]开始MS与基站i和基站1的距离差
        sigma=std_var1^2;
        theta1=TDOA_CHAN(A,r,sigma);%调用TDOA_CHAN函数
        theta2=TDOA_Taylor(A,r,sigma);%调用TDOA_Talor函数
        error1=error1+norm(MS-theta1)^2;%移动台MS估计位置与计算的到的距离的平方
        error2=error2+norm(MS-theta2)^2;
    end
    RMSE1(j)=(error1/number)^(1/2); %均方根误差
    RMSE2(j)=(error2/number)^(1/2);
end
% 绘图
semilogx(std_var,RMSE1,'-O',std_var,RMSE2,'-s');% x轴取对数，X轴范围是1e-2到1,Y轴的范围是变动的
xlabel('测量噪声的标准差(m)');
ylabel('RMSE');
legend('TDOA-CHAN','TDOA-Taylor');

2.TDOA_CHAN函数中代码如下：

function theta = TDOA_CHAN(A,p,sigma)
% A是BSs的坐标
% p是距离测量
% sigma是TOA测量的方差
[m,~]=size(A);%size得到A的行列数赋值给[m,~]，~表示占位，就是只要行m的值
k=sum(A.^2,2);%矩阵A每个元素分别平方，得到新的矩阵，再按行求和，为矩阵K
G1=[A(2:end,:)+ones(m-1,1)*A(1,:),p]; %得到Xm1,Ym1,Rm1,的值，m取值[2,i],构建矩阵Ga
h1=1/2*(p.^2-k(2:end,:)+ones(m-1,1)*k(1,:)); %构建矩阵h
Q=diag(ones(m-1,1)*sigma);%构建TDOA的协方差矩阵
% 初始预计值
theta0=inv(G1'*inv(Q)*G1)*G1'*inv(Q)*h1;%通过一次wls算法进行求解
s=A(2:end,:)-ones(m-1,1)*theta0(1:2,:)';
d=sum(s.^2,2);%矩阵s每个元素分别平方，得到新矩阵，再行求和，最为矩阵d
B1=diag(d.^(1/2));
cov1=B1*Q*B1;
% 第一次wls
theta1=inv(G1'*inv(cov1)*G1)*G1'*inv(cov1)*h1;%进行第一次wls计算
cov_theta1=inv(G1'*inv(cov1)*G1);%得到theta1的协方差矩阵
%第二次wls
G2=[1,0;0,1;1,1];%构建Ga'
h2=[(theta1(1,1)-A(1,1))^2;(theta1(2,1)-A(1,2))^2;theta1(3,1)^2];%构建h'
B2=diag([theta1(1,1)-A(1,1),theta1(2,1)-A(1,2),theta1(3,1)]);%构建b'
cov2=4*B2*cov_theta1*B2;%得到误差矢量的协方差矩阵
theta2=inv(G2'*inv(cov2)*G2)*G2'*inv(cov2)*h2;%运用最大似然估计得到
theta=theta2.^(1/2)+[A(1,1);A(1,2)];%得到MS位置的估计值坐标，以及符号
theta=theta';%转换为（x,y）形式,得到(x0,y0)

3.TDOA_Taylor函数代码如下：

function theta = TDOA_Taylor(A,p,sigma)
% A是BSs的坐标
% p是距离测量
% sigma是TOA测量的方差
% 初始预计值
theta0=TDOA_CHAN(A,p,sigma);%调用TDOA_CHAN得到一个初始的估计位置
delta=norm(theta0);%得到矩阵范数
while norm(delta)>1e-2%得到足够小的值(0.01)
    [m,~]=size(A);%size得到A的行列数赋值给[m,~]，~表示占位，就是只要行m的值
    d=sum((A-ones(m,1)*theta0).^2,2);
    R=d.^(1/2);
    G1=ones(m-1,1)*(A(1,1)-theta0(1,1))/R(1,1)-(A(2:m,1)-theta0(1,1))./R(2:m,:);
    G2=ones(m-1,1)*(A(1,2)-theta0(1,2))/R(1,1)-(A(2:m,2)-theta0(1,2))./R(2:m,:);
    G=[G1,G2];%构建Gt
    h=p-(R(2:m,:)-ones(m-1,1)*R(1,:));%构建ht
    Q=diag(ones(m-1,1)*sigma);%TDOA测量值的协方差矩阵
    delta=inv(G'*inv(Q)*G)*G'*inv(Q)*h;% wls加权最小二乘解
    theta0=theta0+delta';%累加,不满足条件继续循环
end
theta=theta0;%输出对应的位置坐标

---

• 以上在matlab环境下进行的TDOA定位算法性能仿真，与我们老师讲的基本是一致的，有疑问欢迎评论，加油！ 🌙
  
  
  

---