-
Notifications
You must be signed in to change notification settings - Fork 16
/
groupv_fit.m
69 lines (54 loc) · 1.2 KB
/
groupv_fit.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
function [groupv offset] = groupv_fit(dist,time,snr,mingroupv,maxgroupv)
% function to apply weighted least square polynomial fitting to the data
%
err_tol = 2; % number of err std
x=dist(:);
y=time(:);
w = snr(:);
N=1;
if ~exist('mingroupv')
mingroupv = 2;
end
if ~exist('maxgroupv')
maxgroupv = 5;
end
if length(x)~=length(y) || length(x)~=length(w)
disp('The input vectors should have same length!');
return
end
mat = zeros(length(x),N+1);
for i=1:length(x)
for j=1:N+1
mat(i,j) = x(i)^(N+1-j);
end
end
W = diag(w);
para = (mat'*W*mat)\(mat'*W*y);
% Iterative to remove significant outliers.
% first time
err = mat*para - y;
err(find(w==0))=0;
stderr = std(err);
errind = find(abs(err) > stderr*err_tol);
new_w = w;
new_w(errind) = 0;
W = diag(new_w);
para = (mat'*W*mat)\(mat'*W*y);
% second time
err = mat*para - y;
err(find(new_w==0))=0;
stderr = std(err);
errind = find(abs(err) > stderr*(err_tol));
new_w(errind) = 0;
W = diag(new_w);
para = (mat'*W*mat)\(mat'*W*y);
groupv = 1/para(1);
offset = para(2);
if groupv > maxgroupv
groupv = maxgroupv;
offset = sum((y - x/groupv).*w)/sum(w);
elseif groupv < mingroupv;
groupv = mingroupv;
offset = sum((y - x/groupv).*w)/sum(w);
end
end