diff options
Diffstat (limited to 'macros/durbinlevinson.sci')
| -rw-r--r-- | macros/durbinlevinson.sci | 111 |
1 files changed, 86 insertions, 25 deletions
diff --git a/macros/durbinlevinson.sci b/macros/durbinlevinson.sci index 74dbb46..9fb647b 100644 --- a/macros/durbinlevinson.sci +++ b/macros/durbinlevinson.sci @@ -1,27 +1,88 @@ -function y= durbinlevinson(C, varargin) -// Perform one step of the Durbin-Levinson algorithm.. +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 -//This is an Octave function. +// 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. - rhs=argn(2); - if(rhs<1 | rhs>3) - error("Wrong number of input arguments"); - end - select(rhs) - case 1 then - y=callOctave("durbinlevinson",C); - case 2 then - y=callOctave("durbinlevinson",C, varargin(1)); - case 3 then - y=callOctave("durbinlevinson",C, varargin(1), varargin(2)); - end -endfunction
\ No newline at end of file +//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]); |