add CDATA + role in the signal_processing module
[scilab.git] / scilab / modules / signal_processing / help / en_US / levin.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <refentry version="5.0-subset Scilab" xml:id="levin" xml:lang="en"
3           xmlns="http://docbook.org/ns/docbook"
4           xmlns:xlink="http://www.w3.org/1999/xlink"
5           xmlns:svg="http://www.w3.org/2000/svg"
6           xmlns:ns3="http://www.w3.org/1999/xhtml"
7           xmlns:mml="http://www.w3.org/1998/Math/MathML"
8           xmlns:db="http://docbook.org/ns/docbook">
9   <info>
10     <pubdate>$LastChangedDate: 2008-03-26 09:50:39 +0100 (Wed, 26 Mar 2008)
11     $</pubdate>
12   </info>
13
14   <refnamediv>
15     <refname>levin</refname>
16
17     <refpurpose>Toeplitz system solver by Levinson algorithm
18     (multidimensional)</refpurpose>
19   </refnamediv>
20
21   <refsynopsisdiv>
22     <title>Calling Sequence</title>
23
24     <synopsis>[la,sig]=levin(n,cov)</synopsis>
25   </refsynopsisdiv>
26
27   <refsection>
28     <title>Parameters</title>
29
30     <variablelist>
31       <varlistentry>
32         <term>n</term>
33
34         <listitem>
35           <para>A scalar with integer value: the maximum order of the
36           filter</para>
37         </listitem>
38       </varlistentry>
39
40       <varlistentry>
41         <term>cov</term>
42
43         <listitem>
44           <para>A <literal>(nlag*d) x d</literal> matrix. It contains the
45           <literal>Rk</literal> (<literal>d x d</literal> matrices for a
46           <literal>d</literal>-dimensional process) stored in the following
47           way :</para>
48
49           <para><inlinemediaobject>
50               <imageobject>
51                 <imagedata>
52                   <mml:math>
53                   <mml:semantics> <mml:mfenced mml:close=""
54                   mml:open=""> <mml:mtable> <mml:mtr> <mml:mtd> <mml:msub>
55                   <mml:mi>R</mml:mi> <mml:mn>0</mml:mn> </mml:msub> </mml:mtd>
56                   </mml:mtr> <mml:mtr> <mml:mtd> <mml:msub> <mml:mi>R</mml:mi>
57                   <mml:mn>1</mml:mn> </mml:msub> </mml:mtd> </mml:mtr>
58                   <mml:mtr> <mml:mtd> <mml:msub> <mml:mi>R</mml:mi>
59                   <mml:mn>2</mml:mn> </mml:msub> </mml:mtd> </mml:mtr>
60                   <mml:mtr> <mml:mtd> <mml:mo mml:stretchy="false">⋮</mml:mo>
61                   </mml:mtd> </mml:mtr> <mml:mtr> <mml:mtd> <mml:msub>
62                   <mml:mi>R</mml:mi> <mml:mi
63                   mml:fontstyle="italic">nlags</mml:mi> </mml:msub> </mml:mtd>
64                   </mml:mtr> </mml:mtable> </mml:mfenced> <mml:annotation
65                   mml:encoding="StarMath 5.0"> left ( matrix{R_0 ## R_1 ##
66                   R_{2} ## dotsvert ## R_{nlags} } right )</mml:annotation>
67                   </mml:semantics></mml:math>
68                 </imagedata>
69               </imageobject>
70             </inlinemediaobject></para>
71         </listitem>
72       </varlistentry>
73
74       <varlistentry>
75         <term>la</term>
76
77         <listitem>
78           <para>A list, the successively calculated Levinson polynomials
79           (degree 1 to <literal>n</literal>), with coefficients
80           <literal>Ak</literal></para>
81         </listitem>
82       </varlistentry>
83
84       <varlistentry>
85         <term>sig</term>
86
87         <listitem>
88           <para>A list, the successive mean-square errors.</para>
89         </listitem>
90       </varlistentry>
91     </variablelist>
92   </refsection>
93
94   <refsection>
95     <title>Description</title>
96
97     <para>function which solves recursively on n the following Toeplitz system
98     (normal equations)</para>
99
100     <para><inlinemediaobject>
101         <imageobject>
102           <imagedata>
103             <mml:math>
104             <mml:semantics> <mml:mrow> <mml:mrow> <mml:mfenced
105             mml:close="" mml:open=""> <mml:mtable> <mml:mtr> <mml:mtd>
106             <mml:mi>I</mml:mi> </mml:mtd> <mml:mtd> <mml:mrow> <mml:mo
107             mml:stretchy="false">−</mml:mo> <mml:msub> <mml:mi>A</mml:mi>
108             <mml:mn>1</mml:mn> </mml:msub> </mml:mrow> </mml:mtd> <mml:mtd>
109             <mml:mo mml:stretchy="false">⋯</mml:mo> </mml:mtd> <mml:mtd>
110             <mml:mrow> <mml:mo mml:stretchy="false">−</mml:mo> <mml:msub>
111             <mml:mi>A</mml:mi> <mml:mi>n</mml:mi> </mml:msub> </mml:mrow>
112             </mml:mtd> </mml:mtr> </mml:mtable> </mml:mfenced> <mml:mo
113             mml:stretchy="false">∗</mml:mo> <mml:mfenced mml:close=""
114             mml:open=""> <mml:mtable> <mml:mtr> <mml:mtd> <mml:msub>
115             <mml:mi>R</mml:mi> <mml:mn>1</mml:mn> </mml:msub> </mml:mtd>
116             <mml:mtd> <mml:msub> <mml:mi>R</mml:mi> <mml:mn>2</mml:mn>
117             </mml:msub> </mml:mtd> <mml:mtd> <mml:mo
118             mml:stretchy="false">⋯</mml:mo> </mml:mtd> <mml:mtd> <mml:msub>
119             <mml:mi>R</mml:mi> <mml:mi>n</mml:mi> </mml:msub> </mml:mtd>
120             </mml:mtr> <mml:mtr> <mml:mtd> <mml:msub> <mml:mi>R</mml:mi>
121             <mml:mn>0</mml:mn> </mml:msub> </mml:mtd> <mml:mtd> <mml:msub>
122             <mml:mi>R</mml:mi> <mml:mn>1</mml:mn> </mml:msub> </mml:mtd>
123             <mml:mtd> <mml:mo mml:stretchy="false">⋯</mml:mo> </mml:mtd>
124             <mml:mtd> <mml:msub> <mml:mi>R</mml:mi> <mml:mrow>
125             <mml:mi>n</mml:mi> <mml:mo mml:stretchy="false">−</mml:mo>
126             <mml:mn>1</mml:mn> </mml:mrow> </mml:msub> </mml:mtd> </mml:mtr>
127             <mml:mtr> <mml:mtd> <mml:msub> <mml:mi>R</mml:mi> <mml:mrow>
128             <mml:mo mml:stretchy="false">−</mml:mo> <mml:mn>1</mml:mn>
129             </mml:mrow> </mml:msub> </mml:mtd> <mml:mtd> <mml:msub>
130             <mml:mi>R</mml:mi> <mml:mn>0</mml:mn> </mml:msub> </mml:mtd>
131             <mml:mtd> <mml:mo mml:stretchy="false">⋯</mml:mo> </mml:mtd>
132             <mml:mtd> <mml:msub> <mml:mi>R</mml:mi> <mml:mrow>
133             <mml:mi>n</mml:mi> <mml:mo mml:stretchy="false">−</mml:mo>
134             <mml:mn>2</mml:mn> </mml:mrow> </mml:msub> </mml:mtd> </mml:mtr>
135             <mml:mtr> <mml:mtd> <mml:mo mml:stretchy="false">⋮</mml:mo>
136             </mml:mtd> <mml:mtd> <mml:mo mml:stretchy="false">⋮</mml:mo>
137             </mml:mtd> <mml:mtd> <mml:mo mml:stretchy="false">⋯</mml:mo>
138             </mml:mtd> <mml:mtd> <mml:mo mml:stretchy="false">⋮</mml:mo>
139             </mml:mtd> </mml:mtr> <mml:mtr> <mml:mtd> <mml:msub>
140             <mml:mi>R</mml:mi> <mml:mrow> <mml:mn>2</mml:mn> <mml:mo
141             mml:stretchy="false">−</mml:mo> <mml:mi>n</mml:mi> </mml:mrow>
142             </mml:msub> </mml:mtd> <mml:mtd> <mml:msub> <mml:mi>R</mml:mi>
143             <mml:mrow> <mml:mn>3</mml:mn> <mml:mo
144             mml:stretchy="false">−</mml:mo> <mml:mi>n</mml:mi> </mml:mrow>
145             </mml:msub> </mml:mtd> <mml:mtd> <mml:mo
146             mml:stretchy="false">⋯</mml:mo> </mml:mtd> <mml:mtd> <mml:msub>
147             <mml:mi>R</mml:mi> <mml:mn>1</mml:mn> </mml:msub> </mml:mtd>
148             </mml:mtr> <mml:mtr> <mml:mtd> <mml:msub> <mml:mi>R</mml:mi>
149             <mml:mrow> <mml:mn>1</mml:mn> <mml:mo
150             mml:stretchy="false">−</mml:mo> <mml:mi>n</mml:mi> </mml:mrow>
151             </mml:msub> </mml:mtd> <mml:mtd> <mml:msub> <mml:mi>R</mml:mi>
152             <mml:mrow> <mml:mn>2</mml:mn> <mml:mo
153             mml:stretchy="false">−</mml:mo> <mml:mi>n</mml:mi> </mml:mrow>
154             </mml:msub> </mml:mtd> <mml:mtd> <mml:mo
155             mml:stretchy="false">⋯</mml:mo> </mml:mtd> <mml:mtd> <mml:msub>
156             <mml:mi>R</mml:mi> <mml:mn>0</mml:mn> </mml:msub> </mml:mtd>
157             </mml:mtr> </mml:mtable> </mml:mfenced> </mml:mrow> <mml:mo
158             mml:stretchy="false">=</mml:mo> <mml:mn>0</mml:mn> </mml:mrow>
159             <mml:annotation mml:encoding="StarMath 5.0">left ( matrix{I # -A_1
160             # dotsaxis #-A_n} right )* left ( matrix{R_1 # R_2 # dotsaxis #R_n
161             ## R_0 # R_1 # dotsaxis #R_{n-1}## R_{-1} # R_0 # dotsaxis
162             #R_{n-2}## dotsvert # dotsvert# dotsaxis # dotsvert## R_{2-n} #
163             R_{3-n} # dotsaxis #R_{1}## R_{1-n} # R_{2-n} # dotsaxis #R_{0}
164             }right ) = 0</mml:annotation> </mml:semantics></mml:math>
165           </imagedata>
166         </imageobject>
167       </inlinemediaobject></para>
168
169     <para>where {<literal>Rk;k=1:nlag</literal>} is the sequence of
170     <literal>nlag</literal> empirical covariances</para>
171   </refsection>
172
173   <refsection>
174     <title>Examples</title>
175
176     <programlisting role="example"><![CDATA[ 
177 //We use the 'levin' macro for solving the normal equations 
178 //on two examples: a one-dimensional and a two-dimensional process.
179 //We need the covariance sequence of the stochastic process.
180 //This example may usefully be compared with the results from 
181 //the 'phc' macro (see the corresponding help and example in it)
182 //
183 //
184 //1) A one-dimensional process
185 //   -------------------------
186 //
187 //We generate the process defined by two sinusoids (1Hz and 2 Hz) 
188 //in additive Gaussian noise (this is the observed process); 
189 //the simulated process is sampled at 10 Hz (step 0.1 in t, underafter).
190
191 t1=0:.1:100;rand('normal');
192 y1=sin(2*%pi*t1)+sin(2*%pi*2*t1);y1=y1+rand(y1);plot(t1,y1);
193
194 //covariance of y1
195
196 nlag=128;
197 c1=corr(y1,nlag);
198 c1=c1';//c1 needs to be given columnwise (see the section PARAMETERS of this help)
199
200 //compute the filter for a maximum order of n=10
201 //la is a list-type variable each element of which 
202 //containing the filters of order ranging from 1 to n; (try varying n)
203 //in the d-dimensional case this is a matrix polynomial (square, d X d)
204 //sig gives, the same way, the mean-square error
205
206 n=15;
207 [la1,sig1]=levin(n,c1);
208
209 //verify that the roots of 'la' contain the 
210 //frequency spectrum of the observed process y
211 //(remember that y is sampled -in our example 
212 //at 10Hz (T=0.1s) so that we need to retrieve 
213 //the original frequencies (1Hz and 2 Hz) through 
214 //the log and correct scaling by the frequency sampling)
215 //we verify this for each filter order
216
217 for i=1:n, s1=roots(la1(i));s1=log(s1)/2/%pi/.1;
218
219 //now we get the estimated poles (sorted, positive ones only !)
220
221 s1=sort(imag(s1));s1=s1(1:i/2);end;
222
223 //the last two frequencies are the ones really present in the observed 
224 //process ---> the others are "artifacts" coming from the used model size.
225 //This is related to the rather difficult problem of order estimation.
226 //
227 //2) A 2-dimensional process 
228 //   -----------------------
229 //(4 frequencies 1, 2, 3, and 4 Hz, sampled at 0.1 Hz :
230 //   |y_1|        y_1=sin(2*Pi*t)+sin(2*Pi*2*t)+Gaussian noise
231 // y=|   | with : 
232 //   |y_2|        y_2=sin(2*Pi*3*t)+sin(2*Pi*4*t)+Gaussian noise
233
234 d=2;dt=0.1;
235 nlag=64;
236 t2=0:2*%pi*dt:100;
237 y2=[sin(t2)+sin(2*t2)+rand(t2);sin(3*t2)+sin(4*t2)+rand(t2)];
238 c2=[];
239 for j=1:2, for k=1:2, c2=[c2;corr(y2(k,:),y2(j,:),nlag)];end;end;
240 c2=matrix(c2,2,128);cov=[];
241 for j=1:64,cov=[cov;c2(:,(j-1)*d+1:j*d)];end;//covar. columnwise
242 c2=cov;
243
244 //in the multidimensional case, we have to compute the 
245 //roots of the determinant of the matrix polynomial 
246 //(easy in the 2-dimensional case but tricky if d>=3 !). 
247 //We just do that here for the maximum desired 
248 //filter order (n); mp is the matrix polynomial of degree n
249
250 [la2,sig2]=levin(n,c2);
251 mp=la2(n);determinant=mp(1,1)*mp(2,2)-mp(1,2)*mp(2,1);
252 s2=roots(determinant);s2=log(s2)/2/%pi/0.1;//same trick as above for 1D process
253 s2=sort(imag(s2));s2=s2(1:d*n/2);//just the positive ones !
254
255 //There the order estimation problem is seen to be much more difficult !
256 //many artifacts ! The 4 frequencies are in the estimated spectrum 
257 //but beneath many non relevant others.
258  ]]></programlisting>
259   </refsection>
260
261   <refsection>
262     <title>See Also</title>
263
264     <simplelist type="inline">
265       <member><link linkend="phc">phc</link></member>
266     </simplelist>
267   </refsection>
268
269   <refsection>
270     <title>Authors</title>
271
272     <para>G. Le Vey</para>
273   </refsection>
274 </refentry>