diff options
Diffstat (limited to 'macros/levinson.sci')
| -rw-r--r-- | macros/levinson.sci | 213 |
1 files changed, 144 insertions, 69 deletions
diff --git a/macros/levinson.sci b/macros/levinson.sci index ccb0ab9..63ae994 100644 --- a/macros/levinson.sci +++ b/macros/levinson.sci @@ -1,83 +1,158 @@ -//levinson levinson- durbin algorithm -//calling syntax +//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 +// 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 -// Author Debdeep Dey -function [a, v_f, ref_f] = levinson (acf, p) -if ( argn(2)<1 ) - error("Too few input arguments"); - // elseif( length(acf)<2 ) - // error( "levinson: arg 1 (acf) must be vector of length >1\n"); -elseif ( argn(2)>1 & ( ~isscalar(p) | fix(p)~=p ) ) +//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 +// + +//************************************************************************************************** +//______________________________________________version1 code (not working)_________________________ +//__________________________________________________________________________________________________ +//************************************************************************************************** +//function [a, v_f, ref_f] = levinson (acf, p) + +//if ( argn(2)<1 ) +// error("Too few input arguments"); +// elseif( length(acf)<2 ) +// error( "levinson: arg 1 (acf) must be vector of length >1\n"); +//elseif ( argn(2)>1 & ( ~isscalar(p) | fix(p)~=p ) ) +// error( "levinson: arg 2 (p) must be integer >0\n"); +//elseif (isempty(acf)) +// error("R cannot be empty"); +//else +// if ((argn(2) == 1)|(p>=length(acf))) +// p = length(acf) - 1; +// end +// if( size(acf,1)==1 & size(acf,2)>1 ) then +// acf=acf(:); +// +// end force a column vector +// if size(acf,1)>=1 then handles matrix i/p +// +// acf=acf'; +// a=acf; +// rows=size(acf,1); +// for i=1:rows +// acf_temp=acf(i,:); +// acf_temp=acf_temp(:); +// p=length(acf_temp)-1; +// disp(acf_temp); +// if argn(1) < 3 & p < 100 +// +// // Kay & Marple Eqn (2.39) +// +// R = toeplitz(acf_temp(1:p), conj(acf_temp(1:p))); +// an = R \ -acf_temp(2:p+1); +// an= [ 1, an.' ]; +// v_f(i,:)= real( an*conj(acf_temp(1:p+1)) ); +// a(i,:)=an; +// an=[]; +// +// else +// +// // Kay & Marple Eqns (2.42-2.46) +// +// ref = zeros(p,1); +// g = -acf_temp(2)/acf_temp(1); +// +// an = [ g ]; +// v= real( ( 1 - g*conj(g)) * acf_temp(1) ); +// ref(1) = g; +// for t = 2 : p +// g = -(acf_temp(t+1) + an * acf_temp(t:-1:2)) / v; +// an = [ an+g*conj(an(t-1:-1:1)), g ]; +// v = v * ( 1 - real(g*conj(g)) ) ; +// ref(t) = g; +// end +// v_f(i,:)=v; +// v=[]; +// ref_f(:,i)=ref; +// an = [1, an]; +// a(i,:)=an; +// an=[]; +// end end if +// end end for +// end +//end + +// +//endfunction + + +//************************************************************************************************** +//______________________________________________version2 code ( working)____________________________ +//__________________________________________________________________________________________________ +//************************************************************************************************** + + +function [a, v, ref] = levinson(bcf, p) + + 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"); -elseif (isempty(acf)) - error("R cannot be empty"); -else - if ((argn(2) == 1)|(p>=length(acf))) - p = length(acf) - 1; - end - if( size(acf,1)==1 & size(acf,2)>1 ) then - acf=acf(:); - - end // force a column vector - if size(acf,1)>=1 then // handles matrix i/p - - acf=acf'; - a=acf; - rows=size(acf,1); - for i=1:rows - acf_temp=acf(i,:); - acf_temp=acf_temp(:); - p=length(acf_temp)-1; - //disp(acf_temp); - if argn(1) < 3 & p < 100 - - //// Kay & Marple Eqn (2.39) - - R = toeplitz(acf_temp(1:p), conj(acf_temp(1:p))); - an = R \ -acf_temp(2:p+1); - an= [ 1, an.' ]; - v_f(i,:)= real( an*conj(acf_temp(1:p+1)) ); - a(i,:)=an; - an=[]; - - else - - //// Kay & Marple Eqns (2.42-2.46) - - ref = zeros(p,1); - g = -acf_temp(2)/acf_temp(1); - - an = [ g ]; - v= real( ( 1 - g*conj(g)) * acf_temp(1) ); - ref(1) = g; - for t = 2 : p - g = -(acf_temp(t+1) + an * acf_temp(t:-1:2)) / v; - an = [ an+g*conj(an(t-1:-1:1)), g ]; - v = v * ( 1 - real(g*conj(g)) ) ; - ref(t) = g; - end - v_f(i,:)=v; - v=[]; - ref_f(:,i)=ref; - an = [1, an]; - a(i,:)=an; - an=[]; - end //end if - end //end for - end -end + 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 |