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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
|
function [newphi, newv] = durbinlevinson (c, oldphi, oldv)
// Perform one step of the Durbin-Levinson algorithm.
//Calling Sequence
// durbinlevinson(c)
// durbinlevinson(c, oldphi)
// durbinlevinson(c, oldphi, oldv)
//Parameters:
//c: The vector c specifies the autocovariances '[gamma_0, ..., gamma_t]' from lag 0 to t.
//oldphi: It specifies the coefficients based on c(t-1).
//oldv: It specifies the corresponding error.
//Description:
//Perform one step of the Durbin-Levinson.
//If 'oldphi' and 'oldv' are omitted, all steps from 1 to t of the algorithm are performed.
//Example:
//durbinlevinson([1, 2, 3], 1, 2)
//ans = [0.5, 0.5]
funcprot(0);
rhs = argn(2);
if (rhs ~= 1 & rhs ~= 3)
error("durbinlevinson: wrong number of input arguments");
end
if (size(c, 'c') > 1)
c = c';
end
newphi = 0;
newv = 0;
if (rhs == 3)
t = length (oldphi) + 1;
if (length (c) < t + 1)
error ("durbinlevinson: arg1 (c) is too small");
end
if (oldv == 0)
error ("durbinlevinson: arg3 (oldv) must be non-zero");
end
if (size(oldphi, 'r') > 1)
oldphi = oldphi';
end
newphi = zeros (1, t);
newphi(1) = (c(t+1) - oldphi * c(2:t)) / oldv;
for i = 2 : t
newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1);
end
newv = (1 - newphi(1)^2) * oldv;
else
tt = length (c)-1;
oldphi = c(2) / c(1);
oldv = (1 - oldphi^2) * c(1);
for t = 2 : tt
newphi = zeros (1, t);
newphi(1) = (c(t+1) - oldphi * c(2:t)) / oldv;
for i = 2 : t
newphi(i) = oldphi(i-1) - newphi(1) * oldphi(t-i+1);
end
newv = (1 - newphi(1)^2) * oldv;
oldv = newv;
oldphi = newphi;
end
end
endfunction
//input validation:
//assert_checkerror("durbinlevinson()", "durbinlevinson: wrong number of input arguments");
//assert_checkerror("durbinlevinson(1, 2, 3, 4)", "Wrong number of input arguments.");
//assert_checkerror("durbinlevinson([1, 2], [1, 2, 3], 5)", "durbinlevinson: arg1 (c) is too small");
//assert_checkerror("durbinlevinson([1, 2, 3, 4], [1, 2], 0)", "durbinlevinson: arg3 (oldv) must be non-zero");
//test mx1 input:
//assert_checkequal(durbinlevinson([2, 4, 3]), durbinlevinson([2, 4, 3], 2, -6));
//assert_checkequal(durbinlevinson([9, 12, 43, 55], 3, -4), [-1.75, 8.25]);
//assert_checkequal(durbinlevinson([2, 5, 3, 5, 7], [2, 6], -4), [5.75, -32.5, -5.5]);
//assert_checkequal(durbinlevinson([6, 5, 8, 4], [1; 2], -1), [17, -33, -15]);
//test 1xm input:
//assert_checkequal(durbinlevinson([1; 4; 7]), [0.6, 1.6]);
//assert_checkequal(durbinlevinson([3; 6; 7], -5, -2), [-18.5, -97.5]);
//assert_checkequal(durbinlevinson([3; 6; 4; 2; 1; 9], [4; 6; 5], 3), [-19, 99, 120, 81]);
//assert_checkequal(durbinlevinson([6; 5; 8; 4], [1, 2], -1), [17, -33, -15]);
|
