📄 visualcompute.m
字号:
% _________________________
end % if-statement output is MDE or success rates
if ismember(out,[1 2 3 5 6]) % MDBs, MDEs or BNRs
if sys == 1
d = kron(c,D);
dQ = d'*iQy;
PN = N*inv(NQN)*N'*iQy;
if SYSTEM == 4 && ~strcmp(scenario,'gf')
lm = size(M,1);
A = [[M;M2] blkdiag(N,N2) ];
A = A(1:lm,:);
MQN = [MQN MQN2];
else
A = [ M N ];
end
if strcmp(scenario,'sr')
PM = 0;
else
PM = M*iMQM*M'*iQy;
end
im = size(iQy);
if strcmp(output,'outlier')
Q2 = -iMQM * MQN * Q1;
Q3 = iMQM - iMQM * MQN * Q2';
PMN = A * [Q3 Q2 ; Q2' Q1] * A' * iQy;
else
PMN = eye(im);
end
else
d = kron(c,D2);
dQ = d'*iQy2;
PN = N2*inv(NQN2)*N2'*iQy2;
if ~strcmp(scenario,'gf')
lm = size(M,1);
A = [[M;M2] blkdiag(N,N2) ];
A = A(lm+1:end,:);
MQN = [MQN MQN2];
else
A = [ M2 N2 ];
end
if strcmp(scenario,'sr')
PM = 0;
else
PM = M2*iMQM*M2'*iQy2;
end
im = size(iQy2);
if strcmp(output,'outlier')
Q2 = iMQM * MQN * Q1;
Q3 = iMQM - iMQM * MQN * Q2';
PMN = A * [Q3 Q2 ; Q2' Q1] * A' * iQy2;
else
PMN = eye(im);
end
end
% div : denominator in eq. (3.13) and (3.14)
div = sum( slip * (dQ*(eye(im) - (1-slip/no_epochs) * PM - ...
(slip/no_epochs) * PMN))'.*d );
% -----------------------
% minimal detectable bias
% -----------------------
MDB = sqrt ( lam0 ./ div ); % eq.(3.13) or (3.14)
Qhat = [];
if out == 1
if strcmp(output,'slip')
MDB = ml * MDB; % multiplication with wavelength:
% cycle slip MDB in units of range
end
OUT(p) = max(MDB);
titel = ' minimal detectable biases ';
% _________________________
% -------------------
% bias-to-noise ratio
% -------------------
elseif out == 3
% div2: term in parentheses in eq. (3.15) and (3.16)
div2 = slip * (dQ*(eye(size(PN)) - (slip/no_epochs) * PN))'.*d;
BNR = (max(MDB)^2) * div2 - lam0 ; % (3.15) or (3.16)
if abs(BNR)<1e-20, BNR = 0; end
OUT(p) = sqrt(BNR);
titel = ' bias-to-noise ratios (baseline)';
% _________________________
% -------------------------
% minimal detectable effect
% -------------------------
elseif ismember(out,[2 5 6]) % MDE / biased success rates
if SYSTEM == 4
if ~strcmp(scenario,'gf')
M = [M ; M2];
N = blkdiag(N,N2);
iQy = blkdiag(iQy,iQy2);
else
M = blkdiag(M,M2);
N = blkdiag(N,N2);
iQy = blkdiag(iQy,iQy2);
end
if sys == 1
d = [d;zeros(length(c)*m2,size(d,2))];
else
d = [zeros(length(c)*m,size(d,2));d];
end
end
if strcmp(scenario,'gf')||strcmp(scenario,'rr')
Ik = eye(no_epochs);
ek = ones(no_epochs,1);
Qab = -Qahat * kron(ek', MQN'*iMQM ); % eq.(3.24)
iMQMk = kron(Ik,iMQM);
Qb = iMQMk - iMQMk * kron(ek, MQN) * Qab; %eq.(3.25)
else
Qab = -Qahat * MQN'*iMQM ; % eq.(3.24)
Qb = (iMQM - iMQM * MQN * Qab); % eq.(3.25)
end
Qhat = [Qb Qab';Qab Qahat]; % eq.(3.21)
nr = size(N,2); % number of unknown ambiguities
nt = size(Qhat,1); % number of unknowns
nabx=[];
lmdb = length(MDB);
mdetmp = zeros(lmdb,1);
for jdx = 1:lmdb
cnab = d(:,jdx) * MDB(jdx);
% eq.(3.26)
if strcmp(scenario,'gf')||strcmp(scenario,'rr')
nabx(:,jdx) = Qhat * [ kron(cl , M'*iQy*cnab); slip * N'*iQy*cnab ];
else
nabx(:,jdx) = slip*Qhat * [M N]' * iQy * cnab;
end
mdetmp(jdx) = sqrt (sum(nabx(1:3,jdx).^2));
end
jdx2 = find(mdetmp==max(mdetmp));
if length(jdx2)>1, jdx2=jdx2(1); end
if out == 2 % MDE on baseline
OUT(p)= mdetmp(jdx);
titel = ' minimal detectable effects (on baseline) ';
else
bias = iL*Z'*nabx(nt-nr+1:nt,:);
bsr = [];
for si = 1:size(bias,2)
bsr(si) = prod ( normcdf((1-2*bias(:,si))./(2*sqrt(DQz'))) ...
+ normcdf((1+2*bias(:,si))./(2*sqrt(DQz')))-1 );% eq.(3.33)
end
BSR(p) = min(bsr); % biased success rate
if out == 5
OUT(p) = BSR(p);
titel = ' bias affected success rates ';
elseif out == 6
MDE(p) = mdetmp(jdx);
end
end
end
% _________________________
end % if-statement output is 'MDB', 'MDE' or 'BNR'
% ---------------------
% DOPs
% ---------------------
if ismember(out,[7 8])
if no_sv+no_sv2 >= 4;
Amat = [LOS ones(no_sv,1)];
if SYSTEM == 4 %eq.(3.36)
Amat = [Amat zeros(no_sv,1);...
LOS2 zeros(no_sv2,1) ones(no_sv2,1)];
end
qx = inv(Amat' * Amat);
OUT(p) = sqrt(trace(qx(1:out-4,1:out-4))); % PDOP: out=7, GDOP: out=8
if OUT(p)>10, OUT(p)=10; end % eq.(3.34) or (3.35)
titel = ' dilution of precision ';
else
OUT(p) = 10; % maximum value for plotting purposes
end;
end
end % if-statement output is 'number of satellites' or misc.
% _________________________
% waitbar
if widx > widx2*wend/100
waitbar(widx/wend,wh);
widx2 = widx2 + 1;
end
end % loop over grid cells
if lp==1
if out == 6
MDEt(t) = MDE(p);
SRt(t) = SR(p);
BSRt(t) = BSR(p);
else
OUTt(t) = OUT(p);
end
if out~=5, NO_SV(t)= OUT2(p); else NO_SV(t) = OUT(p); end
end
end % loop over time stamps
delete(wh);
if out == 6
if lp>1
MDE = reshape (MDE,ll,lp)';
SR = reshape (SR,ll,lp)';
BSR = reshape (BSR,ll,lp)';
else
MDE = MDEt;
SR = SRt;
BSR = BSRt;
end
else
OUT = reshape (OUT,ll,lp)';
end
% plot OUTPUT
if length(time)==1
if out~= 9, NO_SV = reshape (OUT2,ll,lp)'; else NO_SV = OUT; end
if get(findobj(gcf,'Tag','saveout'),'Value')
filename = get(findobj(gcf,'Tag','filename'),'String');
if out == 6
save(filename,'MDE','SR','BSR','NO_SV','xmin','xmax','ymin','ymax', ...
'starttime')
return
else
save(filename,'OUT','NO_SV','xmin','xmax','ymin','ymax','starttime')
end
end
inter = (ymax-ymin)/6; if inter>1, inter=round(inter); end
figure(figNumber)
axes(findobj(gcf,'Tag','ColorBarAxis'));
cla reset;
set(gca,'Tag','ColorBarAxis');
axes(findobj(gcf,'Tag','PictureAxis'));
cla reset;
set(gca,'Tag','PictureAxis');
imagesc([xmin xmax],[ymax ymin],OUT);
hold on
grid on
load coast;
plot (long,lat,'k-');
title(titel,'color',[1 1 1]);
set(gca,'YDir','normal','Color',[1 1 1], 'XColor',[1 1 1], 'YColor',[1 1 1], ...
'Xgrid','on','Ygrid','on', 'Xtick',xmin:inter:xmax,'Ytick',ymin:inter:ymax, ...
'Xlim',[xmin xmax],'Ylim',[ymin ymax],...
'Box','on','Visible','on','Tag','PictureAxis');
h=colorbar (findobj(gcf,'Tag','ColorBarAxis'));
set(h,'Tag','ColorBarAxis','Color',[1 1 1], 'XColor',[1 1 1], 'YColor',[1 1 1]);
zoom on
% also plot in separate window
if get(findobj(gcf,'Tag','plotnew'),'Value')
figure
imagesc([xmin xmax],[ymax ymin],OUT);
hold on
grid on
load coast;
plot (long,lat,'k-');
title(titel);
set(gca,'YDir','normal','Color',[1 1 1], 'XColor',[0 0 0], 'YColor',[0 0 0], ...
'Xgrid','on','Ygrid','on', 'Xtick',xmin:inter:xmax,'Ytick',ymin:inter:ymax, ...
'Xlim',[xmin xmax],'Ylim',[ymin ymax],...
'Box','on','Visible','on');
colorbar
end
else
if out == 9, NO_SV = OUTt; end
if get(findobj(gcf,'Tag','saveout'),'Value')
filename = get(findobj(gcf,'Tag','filename'),'String');
if out == 6
save(filename,'MDE','SR','BSR','NO_SV','xmin','ymin', ...
'starttime','endtime')
return
else
save(filename,'OUTt','NO_SV','xmin','ymin','starttime','endtime')
end
end
xmin = 1; xmax = max(time);
ymin = min(OUTt); ymax= max(OUTt);
figure
x=xmin:xmax;
st=stairs(x,NO_SV,'g');
set(st,'Linewidth',2)
set(gca,'Ytick',min(NO_SV):1:max(NO_SV),'Ylim',[min(NO_SV)-1 max(NO_SV)+1]);
ds = max(NO_SV) - min(NO_SV) + 2;
ylabel('number of satellites (green)');
if out ~= 9
set(gca,'YaxisLocation','right','xlim',[xmin xmax],'xtick',[],'xticklabel',[]);
hold on
h=axes('position',get(gca,'position'));
hold on
set(h,'color','none')
set(h,'xgrid','off','ygrid','off');
set(h,'YaxisLocation','left')
plot(x,OUTt,'r','Linewidth',2);
yl=get(gca,'Ylim');
dy=yl(2)-yl(1);
yy = 0;
if ismember(out,[4 5])
yl(2)=1;
dy = yl(2)-yl(1);
yy = 0.05*dy/ds;
end
if dy<1e-2, yl(1)=yl(1)-1; yy = 0.01; dy = dy+1; end
set(gca,'xlim',[xmin xmax],'Ytick',yl(1):dy/ds:yl(2),'Ylim',[yl(1) yl(2)+yy]);
ylabel('output (red)');
end
grid on
xlabel('starting epoch ');
title(titel);
end
watchoff(figNumber)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -