Scilab Home page | Wiki | Bug tracker | Forge | Mailing list archives | ATOMS | File exchange
Please login or create an account
Change language to: English - Français - Português - 日本語 -

Please note that the recommended version of Scilab is 6.0.1. This page might be outdated.
See the recommended documentation of this function

Справка Scilab >> Signal Processing > filters > levin


Toeplitz system solver by Levinson algorithm (multidimensional)

Calling Sequence




A scalar with integer value: the maximum order of the filter


A (nlag*d) x d matrix. It contains the Rk (d x d matrices for a d-dimensional process) stored in the following way :


A list, the successively calculated Levinson polynomials (degree 1 to n), with coefficients Ak


A list, the successive mean-square errors.


function which solves recursively on n the following Toeplitz system (normal equations)

where {Rk;k=1:nlag} is the sequence of nlag empirical covariances


//We use the 'levin' macro for solving the normal equations 
//on two examples: a one-dimensional and a two-dimensional process.
//We need the covariance sequence of the stochastic process.
//This example may usefully be compared with the results from 
//the 'phc' macro (see the corresponding help and example in it)
//1) A one-dimensional process
//   -------------------------
//We generate the process defined by two sinusoids (1Hz and 2 Hz) 
//in additive Gaussian noise (this is the observed process); 
//the simulated process is sampled at 10 Hz (step 0.1 in t, underafter).


//covariance of y1

c1=c1';//c1 needs to be given columnwise (see the section PARAMETERS of this help)

//compute the filter for a maximum order of n=10
//la is a list-type variable each element of which 
//containing the filters of order ranging from 1 to n; (try varying n)
//in the d-dimensional case this is a matrix polynomial (square, d X d)
//sig gives, the same way, the mean-square error


//verify that the roots of 'la' contain the 
//frequency spectrum of the observed process y
//(remember that y is sampled -in our example 
//at 10Hz (T=0.1s) so that we need to retrieve 
//the original frequencies (1Hz and 2 Hz) through 
//the log and correct scaling by the frequency sampling)
//we verify this for each filter order

for i=1:n, s1=roots(la1(i));s1=log(s1)/2/%pi/.1;

//now we get the estimated poles (sorted, positive ones only !)


//the last two frequencies are the ones really present in the observed 
//process ---> the others are "artifacts" coming from the used model size.
//This is related to the rather difficult problem of order estimation.
//2) A 2-dimensional process 
//   -----------------------
//(4 frequencies 1, 2, 3, and 4 Hz, sampled at 0.1 Hz :
//   |y_1|        y_1=sin(2*Pi*t)+sin(2*Pi*2*t)+Gaussian noise
// y=|   | with : 
//   |y_2|        y_2=sin(2*Pi*3*t)+sin(2*Pi*4*t)+Gaussian noise

for j=1:2, for k=1:2, c2=[c2;corr(y2(k,:),y2(j,:),nlag)];end;end;
for j=1:64,cov=[cov;c2(:,(j-1)*d+1:j*d)];end;//covar. columnwise

//in the multidimensional case, we have to compute the 
//roots of the determinant of the matrix polynomial 
//(easy in the 2-dimensional case but tricky if d>=3 !). 
//We just do that here for the maximum desired 
//filter order (n); mp is the matrix polynomial of degree n

s2=roots(determinant);s2=log(s2)/2/%pi/0.1;//same trick as above for 1D process
s2=gsort(imag(s2));s2=s2(1:d*n/2);//just the positive ones !

//There the order estimation problem is seen to be much more difficult !
//many artifacts ! The 4 frequencies are in the estimated spectrum 
//but beneath many non relevant others.

See Also

  • phc — Markovian representation
Scilab Enterprises
Copyright (c) 2011-2017 (Scilab Enterprises)
Copyright (c) 1989-2012 (INRIA)
Copyright (c) 1989-2007 (ENPC)
with contributors
Last updated:
Thu Oct 02 14:01:06 CEST 2014