inclinationerr.m

经典的多点法形状复原算法的逐次2点法和逐次3点法的程序。本程序对当存在随机误差和系统误差时产生的评价误差进行了分析。

function inclinationerr;
% evaluation of random error and systematic error including Inclination
% two point and Inclination three point method

NN=[10,20,30,50,75,100,150,200,300,400, 500,600,700,800,900,1000];
M=[];

for kk=1:16
%-------------------------------------
N=NN(kk);
s=1;
d=1;

a=1;
err1=a.*randn(1,N);
err2=a.*randn(1,N);
err3=a.*randn(1,N);
b2=2.*a;
b3=4.*a;

nn=1:N;

tmpsum = 0;
for ii=1:N;
    E2b(ii)= tmpsum;
    tmpsum = tmpsum + b2;
end

tmpsum = 0;
for ii=1:N;
    E2w(ii)= tmpsum;
    tmpsum = tmpsum + (err2(ii)-err1(ii));
end

E3b(1)=0;
E3b(2)=0;
for n=3:N;
    tmpsum = 0;
    for ii=1:n-2;
        tmpsum = tmpsum + n-1-ii;
    end
    E3b(n)= b3.*tmpsum;
end

E3w(1)=0;
E3w(2)=0;
for n=3:N;
    tmpsum = 0;
    for ii=1:n-2;
        tmpsum = tmpsum + (n-1-ii).*(err1(ii)+err3(ii)-2.*err2(ii));
    end
    E3w(n)= tmpsum;
end

% least square fitting to x axis
p_fx = polyfit(nn, E2b, 1);
p_gx = polyfit(nn, zeros(1,N), 1);
E2b_fit = E2b + polyval(p_gx-p_fx, nn);
U_2b=norm(E2b_fit)./sqrt(N)./a;

p_fx = polyfit(nn, E2w, 1);
p_gx = polyfit(nn, zeros(1,N), 1);
E2w_fit = E2w + polyval(p_gx-p_fx, nn);
U_2w=norm(E2w_fit)./sqrt(N)./a;

p_fx = polyfit(nn, E3b, 1);
p_gx = polyfit(nn, zeros(1,N), 1);
E3b_fit = E3b + polyval(p_gx-p_fx, nn);
U_3b=norm(E3b_fit)./sqrt(N)./a;

p_fx = polyfit(nn, E3w, 1);
p_gx = polyfit(nn, zeros(1,N), 1);
E3w_fit = E3w + polyval(p_gx-p_fx, nn);
U_3w=norm(E3w_fit)./sqrt(N)./a;

% figure(21);
% subplot(4,1,1);
% plot(nn,E2w_fit,'r');
% subplot(4,1,2);
% plot(nn,E2b_fit,'r');
% subplot(4,1,3);
% plot(nn,E3w_fit,'r');
% subplot(4,1,4);
% plot(nn,E3b_fit,'r');

M=[M;U_2b,U_2w,U_3w,U_3b];
%M=[M;NN(kk), U_2b,U_2w,U_3w,U_3b];
end
% end for kk
%-----------------------------------------------

% display the results
figure(32);
loglog([NN',NN',NN'], M(:,2:4), '-*');
%semilogy(xx2,pe1,'b',xx2,pe2,'m',xx2,pe3,'r',xx2,pe4,'g');
%plot(xx1,pe1,'--r',xx1,pe2,'--m',xx2,pe3,'--b',xx2,pe4,'--y');
title('MSE error from up and down: sys of 3P, random of 3P, random of 2P, sys of 2P (x axis)')

%
% M=[NN',M];
% dlmwrite('inclinationerr.csv',M,'delimiter',',','precision','%.6f');