基于MATLAB实现的无线传感器网络 CAB定位仿真程序 这是无线传感器节点定位CAB算法的仿真程序+使用说明文档.rar

function x_est = CAB(A, r, R) B = A(:, r~=0);%filter the pseudo-circles with radius of 0 r1 = r(r~=0); center = [B A]; radius = [r1; R]; n = length(radius); tolerance = 1e-5; % this parameter is used to cancel the impact of numerical accuracy. intersect_point = []; for i = 1:n-1 for j = i+1:n a1 = center(:,i); a2 = center(:,j); if norm(a1-a2)==0 continue; end r1 = radius(i); r2 = radius(j); pp = intersect(a1, a2, r1, r2); intersect_point = [intersect_point pp]; end end if length(intersect_point)==0 error('all circles do not intersect each other', '%s', 'exceptional case'); end m = length(r); n = size(intersect_point, 2); degree = zeros(n, 1); for i=1:n pp = intersect_point(:,i); for j=1:m dist = norm(pp-A(:,j)); if dist<=R(j)+tolerance & dist>=r(j)-tolerance degree(i) = degree(i)+1; end end end %degree valid_intersect_point = []; if max(degree)==m %the feasible case valid_intersect_point = intersect_point(:, degree==m); x_est = mean(valid_intersect_point, 2); else %TPDS method in random walk case %%select k points with highest degree in the infeasible case %[deg, index] = sort(degree, 'descend'); %k = 3;%assume k=3, that is, the first three points are selected to estimate the position of unknown node, but you can vary k to obtain different estimation performance %x_est = mean(intersect_point(:, index(1:k)), 2); %variant of TPDS method %select all points with highest degree in the feasible case valid_intersect_point = intersect_point(:, degree==max(degree)); x_est = mean(valid_intersect_point, 2); end return function pp = intersect(a1, a2, r1, r2) x1 = a1(1); y1 = a1(2); x2 = a2(1); y2 = a2(2); t = (x1^2+y1^2)-(x2^2+y2^2)-r1^2+r2^2; bx = 2*(x1-x2); by = 2*(y1-y2); d = r1^2-x1^2-y1^2; if bx==0%the exceptional case y = t/by; a = 1; b = -2*x1; c = -d+y^2-2*y*y1; temp = b^2-4*a*c; if temp==0 x = -b/(2*a); pp = [x; y]; elseif temp>0 x = (-b+sqrt(temp))/(2*a); p1 = [x; y]; x = (-b-sqrt(temp))/(2*a); p2 = [x; y]; pp = [p1 p2]; else pp = []; end return; end a = (by/bx)^2+1; b = -2*by*t/bx^2+2*x1*by/bx-2*y1; c = (t/bx)^2-2*x1*t/bx-d; temp = b^2-4*a*c; if temp==0%only one intersection point y = -b/(2*a); x = (t-by*y)/bx; pp = [x; y]; elseif temp>0%two intersection points y = (-b+sqrt(temp))/(2*a); x = (t-by*y)/bx; p1 = [x; y]; y = (-b-sqrt(temp))/(2*a); x = (t-by*y)/bx; p2 = [x; y]; pp = [p1 p2]; else pp = []; %error('exception', '%s', 'there exists exception'); end return;