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
|
function [a, v, ref] = levinson(bcf, p)
//Levinson- Durbin Recurssion Algorithm
//calling syntax
//a = levinson(r)
//a = levinson(r,n)
//[a,e] = levinson(r,n)
//[a,e,k] = levinson(r,n)
// where
// a is the coefficients of a length(r)-1 order autoregressive linear process
//e is the prediction error when order is n
// k is a column vector containing the reflection coefficients of length n
//Example :
//a = [1 0.1 -0.8]; //Estimate the coefficients of an autoregressive process given by x(n) = 0.1x(n-1) - 0.8x(n-2) + w(n)
//
//v = 0.4;
//w = sqrt(v)*rand(15000,1,"normal");
//x = filter(1,a,w);
//
//[r,lg] = xcorr(x,'biased');
//r(lg<0) = [];
//
//ar = levinson(r,length(a)-1)
// // Output :---
// ar =
//
// 1. 0.0983843 - 0.7929775
//
funcprot(0);
nargin = argn(2);
nargout = argn(1);
[rows columns] = size(bcf)
if ( nargin<1 )
error("Wrong input argument ");
elseif( ~isvector(bcf) | length(bcf)<2 )
error( "levinson: arg 1 (bcf) must be vector of length >1\n");
elseif ( nargin>1 & ( ~isscalar(p) | fix(p)~=p ) )
error( "levinson: arg 2 (p) must be integer >0\n");
else
if ((nargin == 1)|(p>=length(bcf))) p = length(bcf) - 1; end
if( columns >1 ) bcf=bcf(:); end
if nargout < 3 & p < 100
// ## direct solution [O(p^3), but no loops so slightly faster for small p]
// ## Kay & Marple Eqn (2.39)
R = toeplitz(bcf(1:p), conj(bcf(1:p)));
a = R \ -bcf(2:p+1);
a = [ 1, a.' ];
v = real( a*conj(bcf(1:p+1)) );
else
// ## durbin-levinson [O(p^2), so significantly faster for large p]
// ## Kay & Marple Eqns (2.42-2.46)
ref = zeros(p,1);
g = -bcf(2)/bcf(1);
a = [ g ];
v = real( ( 1 - g*conj(g)) * bcf(1) );
ref(1) = g;
for t = 2 : p
g = -(bcf(t+1) + a * bcf(t:-1:2)) / v;
a = [ a+g*conj(a(t-1:-1:1)), g ];
v = v * ( 1 - real(g*conj(g)) ) ;
ref(t) = g;
end
a = [1, a];
end
end
endfunction
|
