constrained_weighted_NLOS_ireative.m
上传用户:doryuen
上传日期:2013-10-30
资源大小:23k
文件大小:3k
源码类别:

通讯/手机编程

开发平台:

Matlab

  1. % 2007.4.1
  2.     
  3.     clc;
  4.     clear;
  5.     
  6.     n=4;      % the number of base stations
  7.     m=100;    % the number of measurements
  8.     
  9.     delta=60;    % variance of measurement noise
  10.     
  11. %     bx(1:n)=0; by(1:5)=0;   % the coordinate of base staions
  12.     
  13.     radius=1000;
  14.     bx=[0   1.5*radius         0               -1.5*radius        -1.5*radius          0               1.5*radius];
  15.     by=[0   sqrt(3)*radius/2   sqrt(3)*radius   sqrt(3)*radius/2   -sqrt(3)*radius/2  -sqrt(3)*radius  -sqrt(3)*radius/2];
  16.     
  17.     p=[1 0 0;
  18.         0 1 0;
  19.         0 0 0];
  20.     q=[0 0 -1]';
  21.     
  22.        noise(1:m,1:n)=0;
  23.        for a=1:n
  24.            noise(:,a)=(add_noise((0.5+0.5*rand(1,m))*300,0.01,800,1000))';
  25.        end     
  26.     
  27.        dt(1:m,1:n)=0;
  28.        d(1:m,1:n)=0;
  29.        
  30.      for loop=1:m;
  31.         
  32.         txr=1*rand*radius;    tx_angle=rand*2*pi;
  33.         tx(loop)=txr*cos(tx_angle);  ty(loop)=txr*sin(tx_angle);
  34. %         tx(loop)=500;  ty(loop)=800;
  35.         for k=1:n
  36.             dt(loop,k)=sqrt((tx(loop)-bx(k))^2+(ty(loop)-by(k))^2);
  37. %             d(loop,k)=0.95*dt(loop,k)-0*delta*rand+00;         % distance measurement
  38. %             d(loop,k)=dt(loop,k)+exprnd(1,1,1)*100;         % distance measurement
  39.             d(loop,k)=dt(loop,k)+(0.5+0.5*rand)*300;
  40. %             d(loop,k)=dt(loop,k)+noise(loop,k);         % distance measurement
  41.         end
  42.         clear k
  43.             
  44.         w=[]; w(1:n,1)=1;
  45.         
  46.             A(1:n,1:3)=0;
  47.             for k=1:n
  48.                 A(k,:)=[bx(k) by(k) -0.5];
  49.             end
  50.         b(1:n,1)=0; n_f=0; flag=0;
  51.            
  52.             while flag==0
  53.                 n_f=n_f+1;
  54.                 b=0.5*[ (bx(1:n).^2)'+(by(1:n).^2)'-(w(:,n_f).*d(loop,:)').^2];
  55.                 mroot=root_find(b,bx(1:n),by(1:n),p,q);
  56.         
  57.                 if length(mroot)>=1 && mroot(1)~= -1
  58.                     theta=inv(A'*A+mroot*p)*(A'*b-0.5*mroot*q);
  59.                     x(n_f)=theta(1);
  60.                     y(n_f)=theta(2);
  61.                 end
  62.                 est_d=sqrt((x(n_f)-bx(1:n)).^2+(y(n_f)-by(1:n)).^2);
  63. %                 w(:,n_f+1)=(est_d./(w(:,n_f)'.*d(loop,:)))';
  64.                 w(:,n_f+1)=(est_d./(d(loop,:)))';
  65. %                 if max(w(:,n_f+1))>=1
  66. %                     w(:,n_f+1)=w(:,n_f+1)/((1+0.1)*max(w(:,n_f+1)));       %normalized
  67. %                 end                 
  68.                 
  69.                 if n_f<1000 && n_f>100 && abs(mean(diff(x(round(0.9*n_f):n_f))))<5e-2 && abs(mean(diff(y(round(0.9*n_f):n_f))))<5e-2
  70.                     flag=1;
  71.                 end
  72.             end
  73.             
  74.             est_x=mean(x(round(0.9*n_f):n_f));
  75.             est_y=mean(y(round(0.9*n_f):n_f));
  76.         
  77.                 [tx(loop) ty(loop)];
  78.                 [est_x est_y];
  79.                 err_Huber(loop)=sqrt((tx(loop)-est_x)^2+(ty(loop)-est_y)^2);
  80.                 [ex,ey]=TOA_LS(radius,bx(1:n),by(1:n),d(loop,:));
  81.                 [ex,ey];
  82.                 err_LS(loop)=sqrt((tx(loop)-ex)^2+(ty(loop)-ey)^2);
  83.         
  84.         loop
  85.     end
  86.         
  87.     sqrt(sum(err_Huber.^2)/m)
  88.     sqrt(sum(err_LS.^2)/m)
  89.     
  90.