dop_change_interval.m

计算dop（dillution of precision）参数的程序

%% DOP  Change Rate Plot
% analyze dop recalculation interval by investigating optmum dop change
% rate.
% Input constellation simulation PRN and LOS matrix, smapled each 6s,
% of a static receiver in N39 E72 height=0 with 0 elevation mask angle. PRNs and
% LOSs are grouped by 12 channel infomation, so 0 means that channel is not
% assigned due to last SV simulated has fallen and a new one has not been
% assigned to that channel.
% The resulting DOP change rate could be the slowest or the longest due to
% static and 0 degree mask angle, which will give a upper bound of DOP
% calculation interval. Surely more frequently calculated the more
% sensitive to opitmum DOP change, but certainly there will be more power
% used so to have a worst case estimation of the interval is very
% meanningful.
% 2009.02.16 Nick @ETS

%% read in the pre-generated and formatted CSV file
global opt_mindop_sv_combination;
global try_times;

PRN_LOS = csvread('PRN_LOS_12chan_2Days.csv');
len = length(PRN_LOS(:,1))/12;
prn = zeros(12, len);
los = zeros(12, 3, len);
for i = 1:len
    prn(:, i) = PRN_LOS((i-1)*12+1:i*12,1);
    los(:,:,i) = PRN_LOS((i-1)*12+1:i*12,2:4);
end
opt_mindop_sv_combination = zeros(8, 8, len);
%% calculate dop for each snapshot
N = 12;
optim_dops = ones(8-4+1, len);
% simulated SV constellaiton
for try_times = 1:len
    
    % get constellation LOS matrix
    visible_los = los(:,:,try_times);

    for M = 4:8
        % find the optimum combination with minimum PDOP for each M
        optim_dops(M-3, try_times) = opti_pdop_sel_no_savefile( M, visible_los, N );
    end
    
    try_times
end
%% plot 
%by time
dop_choose_7 = optim_dops(4,:);
do_choose_8 = optim_dops(5,:);
for i = 1:len
    if optim_dops(4,i) == 100
        optim_dops(4,i) = 0;
    end
    if optim_dops(5,i) == 100
        optim_dops(5,i) = 0;
    end
end

index = 1:len;
for i = 1:8-4+1
    figure(i);
    plot(index, optim_dops(i,:),'.');
end

%by statistical
mean_dop = zeros(1, 8-4+1);
max_dop = zeros(1, 8-4+1);
std_dop = zeros(1, 8-4+1);
for i = 1:8-4+1
    x = 0:0.01:3;
    figure(40+i);
    hist(optim_dops(i,:),x)
    mean_dop(i) = mean(optim_dops(i,:));
    max_dop(i) = max(optim_dops(i,:));
    std_dop(i) = std(optim_dops(i,:));
end
%% Analyze the interval
% analyze the dop change by optimum SV combination chage and see how rapid
% this happend.

%map channel number of selected SV to prn number
opt_mindop_sv_combination_prn = zeros(5,8,len);
for i = 1:len
    for j = 4:8
        for k = 1:j
            if opt_mindop_sv_combination(j,k,i) ~= 0
                opt_mindop_sv_combination_prn(j-3,k,i) = prn(opt_mindop_sv_combination(j,k,i),i);
            end
        end
    end
end
%%
%calculate prn combination of selected SV change rate
interval = zeros(5, len);
cnt1 = 1;
cnt2 = 1;
cnt3 = 1;
change_4_1 = 0;
change_4_2 = 0;
change_5_1 = 0;
change_5_2 = 0;
change_6_1 = 0;
change_6_2 = 0;
for i = 1:len-1
    if sum(opt_mindop_sv_combination_prn(1,:,i) ~= opt_mindop_sv_combination_prn(1,:,i+1))~= 0
        change_4_1 = change_4_2;
        change_4_2 = i;
        interval(1, cnt1) = (change_4_2 - change_4_1)*6;%6 seconds betweend two samples
        cnt1 = cnt1 + 1;
    end
    if sum(opt_mindop_sv_combination_prn(2,:,i) ~= opt_mindop_sv_combination_prn(2,:,i+1)) ~= 0
        change_5_1 = change_5_2;
        change_5_2 = i;
        interval(2, cnt2) = (change_5_2 - change_5_1)*6;%6 seconds betweend two samples
        cnt2 = cnt2 + 1;
    end
    if sum(opt_mindop_sv_combination_prn(3,:,i) ~= opt_mindop_sv_combination_prn(3,:,i+1)) ~= 0
        change_6_1 = change_6_2;
        change_6_2 = i;
        interval(3, cnt3) = (change_6_2 - change_6_1)*6;%6 seconds betweend two samples
        cnt3 = cnt3 + 1;
    end
end
%%
mark = zeros(1,len);
temp = 0;
for i = 1:len
    if interval(1,i) ~= 0
        temp = temp + interval(1,i)/6;
        mark(temp) = 3;
    end
end
figure(1);
hold on
plot(index, mark,'r.');

mark = zeros(1,len);
temp = 0;
for i = 1:len
    if interval(2,i) ~= 0
        temp = temp + interval(2,i)/6;
        mark(temp) = 3;
    end
end
figure(2);
hold on
plot(index, mark,'r.');

mark = zeros(1,len);
temp = 0;
for i = 1:len
    if interval(3,i) ~= 0
        temp = temp + interval(3,i)/6;
        mark(temp) = 3;
    end
end
figure(3);
hold on
plot(index, mark,'r.');

%%
interval_in_second_choose4 = interval(1, 2:227);
min_intvl_4 = min(interval_in_second_choose4)
max_intvl_4 = max(interval_in_second_choose4)
mean_intvl_4 = mean(interval_in_second_choose4)
std_intvl_4 = std(interval_in_second_choose4)

interval_in_second_choose5 = interval(2, 2:253);
min_intvl_5 = min(interval_in_second_choose5)
max_intvl_5 = max(interval_in_second_choose5)
mean_intvl_5 = mean(interval_in_second_choose5)
std_intvl_5 = std(interval_in_second_choose5)

interval_in_second_choose6 = interval(3, 2:230);
min_intvl_6 = min(interval_in_second_choose6)
max_intvl_6 = max(interval_in_second_choose6)
mean_intvl_6 = mean(interval_in_second_choose6)
std_intvl_6 = std(interval_in_second_choose6)
%%

x = 1:12:451;
figure(61);
n=histc(interval_in_second_choose4,x);
bar(n);
figure(62);
n=histc(interval_in_second_choose5,x);
bar(n);
figure(63);
n=histc(interval_in_second_choose6,x);
bar(n);
%%

proposed_recalculate_interval = 60;%calculate every proposed_recalculate_interval seconds
violation_4 = 0;
for i = 1:length(interval_in_second_choose4)
    if interval_in_second_choose4(i) <= proposed_recalculate_interval
        violation_4 = violation_4 + 1;
    end
end
violation_4 = violation_4/i;%get the percentage

violation_5 = 0;
for i = 1:length(interval_in_second_choose5)
    if interval_in_second_choose5(i) <= proposed_recalculate_interval
        violation_5 = violation_5 + 1;
    end
end
violation_5 = violation_5/i;%get the percentage

violation_6 = 0;
for i = 1:length(interval_in_second_choose6)
    if interval_in_second_choose6(i) <= proposed_recalculate_interval
        violation_6 = violation_6 + 1;
    end
end
violation_6 = violation_6/i;%get the percentage
%%
% evaluate dop degrade caused by using the same SVs for proposed_recalculate_interval time, 
% and not changing to the currently optimum ones
proposed_recalculate_interval = 60;
actual_dops = ones(8-4+1, len);
index0 = proposed_recalculate_interval/6;
for i = 1:len/index0
    for k = 1:index0
        for m = 4:6
            sel_los = zeros(m, 3);
            for j = 1:m
                index_los = find(prn(:,(i-1)*index0+k)==opt_mindop_sv_combination_prn(m-3,j,(i-1)*index0+1));
                if ~isempty(index_los)
                    sel_los(j,:) = los( index_los, :, (i-1)*index0+k );
                else
                    sel_los(j,:) = zeros(1,3);
                end
            end
            dop_of_los = sqrt(trace(inv(sel_los'*sel_los)));
            actual_dops(m-3,(i-1)*index0+k) = dop_of_los;
        end
    end
end

index = 1:len;
actual_dops(1,2320)=0;
for i = 1:6-4+1
    figure(100+i);
    temp = actual_dops(i,:)-optim_dops(i,:);
    mean(temp)
    plot(index, temp,'.');
end
%% use quasi and my method to process the simulated constellation data
% and compare the result
% Choose 6 from visible SVs

M = 5;
N = 12;
quasi_optim_dops = zeros(1,try_times);
mynew_optim_dops = zeros(1,try_times);
for try_times = 1:len
    if try_times == 445
        try_times
    end
    visible_los = zeros(1,3);
     for i = 1:N
         if los(i,:,try_times) ~= [ 1 0 0 ]
            visible_los = [visible_los;los(i,:,try_times)];
         end
     end
    visible_los(1,:) = [];
    N_visible = length(visible_los(:,1));
    
    % calculate the quasi opti sets
    quasi_sel_flag = quasi( M, visible_los, N_visible );
    H = zeros(1,3); % matrix to hold selected LOS vectors
    counter = 1;
    for i = 1:N_visible
        if quasi_sel_flag(i) == 1
            H = [ H; visible_los(i,:) ];
            counter = counter + 1;
        end
    end
    quasi_optim_dops(try_times) = sqrt(trace(inv(H'*H)));
    % calculate my new method dop results
    temp = quasi_new_method_optmized_for_simulation( visible_los, N_visible);
    mynew_optim_dops(try_times) = temp(M);
end

%analyze
quasi_dop_ratio = quasi_optim_dops./optim_dops(3,:);
my_dop_ratio = mynew_optim_dops./optim_dops(3,:);



%plot choose 6 from N
x = 1:0.001:1.5;
figure(9);
hist(quasi_dop_ratio,x)
figure(49);
hist(my_dop_ratio,x)

mean_rate(1,:) = mean(quasi_dop_ratio);
mean_rate(2,:) = mean(my_dop_ratio);
max_rate(1,:) = max(quasi_dop_ratio);
max_rate(2,:) = max(my_dop_ratio);