← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to man/signal/levin.cat

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
man/signal/levin.cat
 
1
levin            Scilab Group            Scilab Function              levin
 
2
NAME
 
3
   levin - Toeplitz system solver by Levinson algorithm (multidimensional)
 
4
  
 
5
CALLING SEQUENCE
 
6
 [la,sig]=levin(n,cov)
 
7
PARAMETERS
 
8
 n            : maximum order of the filter
 
9
              
 
10
 cov          : matrix containing the 
 
11
              
 
12
             (d x d matrices for a d-dimensional process). It must be
 
13
              given the following way :
 
14
              
 
15
               |  R0  |
 
16
               |  R1  |
 
17
               |  R2  |
 
18
               |  .   |
 
19
               |  .   |
 
20
               | Rnlag|
 
21
 la           : list-type variable, giving the successively calculated
 
22
              Levinson polynomials (degree 1 to n), with coefficients Ak
 
23
              
 
24
 sig          : list-type variable, giving the successive mean-square
 
25
              errors.
 
26
              
 
27
DESCRIPTION
 
28
   function which solves recursively on n the following Toeplitz system
 
29
  (normal equations)
 
30
  
 
31
                 |R1   R2   . . . Rn  |
 
32
                 |R0   R1   . . . Rn-1|
 
33
                 |R-1  R0   . . . Rn-2|
 
34
                 | .    .   . . .  .  |
 
35
 |I -A1 .. -An|  | .    .   . . .  .  | = 0
 
36
                 | .    .   . . .  .  |
 
37
                 | .    .   . . .  .  |
 
38
                 |R2-n R3-n . . . R1  |
 
39
                 |R1-n R2-n . . . R0  |
 
40
   where {Rk;k=1,nlag} is the sequence of nlag empirical covariances
 
41
  
 
42
AUTHOR
 
43
   G. Le Vey
 
44
  
 
45
EXAMPLE
 
46
 //We use the 'levin' macro for solving the normal equations 
 
47
 //on two examples: a one-dimensional and a two-dimensional process.
 
48
 //We need the covariance sequence of the stochastic process.
 
49
 //This example may usefully be compared with the results from 
 
50
 //the 'phc' macro (see the corresponding help and example in it)
 
51
 //
 
52
 //
 
53
 //1) A one-dimensional process
 
54
 //   -------------------------
 
55
 //
 
56
 //We generate the process defined by two sinusoids (1Hz and 2 Hz) 
 
57
 //in additive Gaussian noise (this is the observed process); 
 
58
 //the simulated process is sampled at 10 Hz (step 0.1 in t, underafter).
 
59
 //
 
60
 t1=0:.1:100;rand('normal');
 
61
 y1=sin(2*%pi*t1)+sin(2*%pi*2*t1);y1=y1+rand(y1);plot(t1,y1);
 
62
 //
 
63
 //covariance of y1
 
64
 //
 
65
 nlag=128;
 
66
 c1=corr(y1,nlag);
 
67
 c1=c1';//c1 needs to be given columnwise (see the section PARAMETERS of this help)
 
68
 //
 
69
 //compute the filter for a maximum order of n=10
 
70
 //la is a list-type variable each element of which 
 
71
 //containing the filters of order ranging from 1 to n; (try varying n)
 
72
 //in the d-dimensional case this is a matrix polynomial (square, d X d)
 
73
 //sig gives, the same way, the mean-square error
 
74
 //
 
75
 n=15;
 
76
 [la1,sig1]=levin(n,c1);
 
77
 //
 
78
 //verify that the roots of 'la' contain the 
 
79
 //frequency spectrum of the observed process y
 
80
 //(remember that y is sampled -in our example 
 
81
 //at 10Hz (T=0.1s) so that we need to retrieve 
 
82
 //the original frequencies (1Hz and 2 Hz) through 
 
83
 //the log and correct scaling by the frequency sampling)
 
84
 //we verify this for each filter order
 
85
 //
 
86
 for i=1:n, s1=roots(la1(i));s1=log(s1)/2/%pi/.1;
 
87
 //
 
88
 //now we get the estimated poles (sorted, positive ones only !)
 
89
 //
 
90
 s1=sort(imag(s1));s1=s1(1:i/2);end;
 
91
 //
 
92
 //the last two frequencies are the ones really present in the observed 
 
93
 //process ---> the others are "artifacts" coming from the used model size.
 
94
 //This is related to the rather difficult problem of order estimation.
 
95
 //
 
96
 //2) A 2-dimensional process 
 
97
 //   -----------------------
 
98
 //(4 frequencies 1, 2, 3, and 4 Hz, sampled at 0.1 Hz :
 
99
 //   |y_1|        y_1=sin(2*Pi*t)+sin(2*Pi*2*t)+Gaussian noise
 
100
 // y=|   | with : 
 
101
 //   |y_2|        y_2=sin(2*Pi*3*t)+sin(2*Pi*4*t)+Gaussian noise
 
102
 //
 
103
 //
 
104
 d=2;dt=0.1;
 
105
 nlag=64;
 
106
 t2=0:2*%pi*dt:100;
 
107
 y2=[sin(t2)+sin(2*t2)+rand(t2);sin(3*t2)+sin(4*t2)+rand(t2)];
 
108
 c2=[];
 
109
 for j=1:2, for k=1:2, c2=[c2;corr(y2(k,:),y2(j,:),nlag)];end;end;
 
110
 c2=matrix(c2,2,128);cov=[];
 
111
 for j=1:64,cov=[cov;c2(:,(j-1)*d+1:j*d)];end;//covar. columnwise
 
112
 c2=cov;
 
113
 //
 
114
 //in the multidimensional case, we have to compute the 
 
115
 //roots of the determinant of the matrix polynomial 
 
116
 //(easy in the 2-dimensional case but tricky if d>=3 !). 
 
117
 //We just do that here for the maximum desired 
 
118
 //filter order (n); mp is the matrix polynomial of degree n
 
119
 //
 
120
 [la2,sig2]=levin(n,c2);
 
121
 mp=la2(n);determinant=mp(1,1)*mp(2,2)-mp(1,2)*mp(2,1);
 
122
 s2=roots(determinant);s2=log(s2)/2/%pi/0.1;//same trick as above for 1D process
 
123
 s2=sort(imag(s2));s2=s2(1:d*n/2);//just the positive ones !
 
124
 //
 
125
 //There the order estimation problem is seen to be much more difficult !
 
126
 //many artifacts ! The 4 frequencies are in the estimated spectrum 
 
127
 //but beneath many non relevant others.
 
128
 //
 
129
SEE ALSO
 
130
   phc  
 
131