* Bug #11670 to 11676 fixed - Add examples in the signal processing help pages.
[scilab.git] / scilab / modules / signal_processing / help / en_US / wiener.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org"  version="5.0-subset Scilab" xml:lang="en" xml:id="wiener">
3     <refnamediv>
4         <refname>wiener</refname>
5         <refpurpose>  Wiener estimate</refpurpose>
6     </refnamediv>
7     <refsynopsisdiv>
8         <title>Calling Sequence</title>
9         <synopsis>[xs,ps,xf,pf]=wiener(y,x0,p0,f,g,h,q,r)</synopsis>
10     </refsynopsisdiv>
11     <refsection>
12         <title>Arguments</title>
13         <variablelist>
14             <varlistentry>
15                 <term>f, g, h</term>
16                 <listitem>
17                     <para>
18                         system matrices in the interval <literal>[t0,tf]</literal>
19                     </para>
20                     <variablelist>
21                         <varlistentry>
22                             <term>f</term>
23                             <listitem>
24                                 <para>
25                                     =<literal>[f0,f1,...,ff]</literal>, and <literal>fk</literal> is a nxn matrix
26                                 </para>
27                             </listitem>
28                         </varlistentry>
29                         <varlistentry>
30                             <term>g</term>
31                             <listitem>
32                                 <para>
33                                     =<literal>[g0,g1,...,gf]</literal>, and <literal>gk</literal> is a nxn matrix
34                                 </para>
35                             </listitem>
36                         </varlistentry>
37                         <varlistentry>
38                             <term>h</term>
39                             <listitem>
40                                 <para>
41                                     =<literal>[h0,h1,...,hf]</literal>, and <literal>hk</literal> is a mxn matrix
42                                 </para>
43                             </listitem>
44                         </varlistentry>
45                     </variablelist>
46                 </listitem>
47             </varlistentry>
48             <varlistentry>
49                 <term>q, r</term>
50                 <listitem>
51                     <para>covariance matrices of dynamics and observation noise</para>
52                     <variablelist>
53                         <varlistentry>
54                             <term>q</term>
55                             <listitem>
56                                 <para>
57                                     =<literal>[q0,q1,...,qf]</literal>, and <literal>qk</literal> is a nxn matrix
58                                 </para>
59                             </listitem>
60                         </varlistentry>
61                         <varlistentry>
62                             <term>r</term>
63                             <listitem>
64                                 <para>
65                                     =<literal>[r0,r1,...,rf]</literal>, and <literal>gk</literal> is a mxm matrix
66                                 </para>
67                             </listitem>
68                         </varlistentry>
69                     </variablelist>
70                 </listitem>
71             </varlistentry>
72             <varlistentry>
73                 <term>x0, p0</term>
74                 <listitem>
75                     <para>initial state estimate and error variance</para>
76                 </listitem>
77             </varlistentry>
78             <varlistentry>
79                 <term>y</term>
80                 <listitem>
81                     <para>
82                         observations in the interval <literal>[t0,tf]</literal>. <literal>y=[y0,y1,...,yf]</literal>, and <literal>yk</literal> is a column m-vector
83                     </para>
84                 </listitem>
85             </varlistentry>
86             <varlistentry>
87                 <term>xs</term>
88                 <listitem>
89                     <para>
90                         Smoothed state estimate <literal>xs= [xs0,xs1,...,xsf]</literal>, and <literal>xsk</literal> is a column n-vector
91                     </para>
92                 </listitem>
93             </varlistentry>
94             <varlistentry>
95                 <term>ps</term>
96                 <listitem>
97                     <para>
98                         Error covariance of smoothed estimate <literal>ps=[p0,p1,...,pf]</literal>, and <literal>pk</literal> is a nxn matrix
99                     </para>
100                 </listitem>
101             </varlistentry>
102             <varlistentry>
103                 <term>xf</term>
104                 <listitem>
105                     <para>
106                         Filtered state estimate <literal>xf= [xf0,xf1,...,xff]</literal>, and <literal>xfk</literal> is a column n-vector
107                     </para>
108                 </listitem>
109             </varlistentry>
110             <varlistentry>
111                 <term>pf</term>
112                 <listitem>
113                     <para>
114                         Error covariance of filtered estimate <literal>pf=[p0,p1,...,pf]</literal>, and <literal>pk</literal> is a nxn matrix
115                     </para>
116                 </listitem>
117             </varlistentry>
118         </variablelist>
119     </refsection>
120     <refsection>
121         <title>Description</title>
122         <para>
123             function which gives the Wiener estimate using
124             the forward-backward Kalman filter formulation
125         </para>
126     </refsection>
127     <refsection>
128         <title>Sample</title>
129         <scilab:image>
130             m0=[10 10]';
131             p0=[100 0;0 100];
132             f=[1.1 50.1;0 0.8];
133             g=[1 0;0 1];
134             h=[1 0;0 1];
135             [hi,hj]=size(h);
136             q=[.01 0;0 0.01];
137             r=20*eye(2,2);
138             rand("seed",66);
139             rand("normal");
140             p0c=chol(p0);
141             x0=m0+p0c'*rand(ones(m0));
142             y=h*x0+chol(r)'*rand(ones(1:hi))';
143             yt=y;
144             x=x0;
145             ft=[f];
146             gt=[g];
147             ht=[h];
148             qt=[q];
149             rt=[r];
150             n=10;
151             for k=1:n
152             [x1,y]=system(x0,f,g,h,q,r);
153             x=[x x1];
154             yt=[yt y];
155             x0=x1;
156             ft=[ft f];
157             gt=[gt g];
158             ht=[ht h];
159             qt=[qt q];
160             rt=[rt r];
161             end
162             [xs,ps,xf,pf]=wiener(yt,m0,p0,ft,gt,ht,qt,rt);
163             a=min([x(1,:)-2*sqrt(ps(1,1:2:2*(n+1))),xf(1,:),xs(1,:)]);
164             b=max([x(1,:)+2*sqrt(ps(1,1:2:2*(n+1))),xf(1,:),xs(1,:)]);
165             c=min([x(2,:)-2*sqrt(ps(2,2:2:2*(n+1))),xf(2,:),xs(2,:)]);
166             d=max([x(2,:)+2*sqrt(ps(2,2:2:2*(n+1))),xf(2,:),xs(2,:)]);
167             xmargin=max([abs(a),abs(b)]);
168             ymargin=max([abs(c),abs(d)]);
169             a=-0.1*xmargin+a;
170             b=.1*xmargin+b;
171             c=-0.1*ymargin+c;
172             d=.1*ymargin+d;
173             scf();
174             plot([a a b],[d c c]);
175             plot2d(x(1,:)',x(2,:)',[2],"000")
176             plot2d(xf(1,:)',xf(2,:)',[2],"000")
177             plot2d(xs(1,:)',xs(2,:)',[2],"000")
178             plot2d(xs(1,:)',xs(2,:)',[-2],"000")
179             plot2d(xf(1,:)',xf(2,:)',[-3],"000")
180             plot2d(x(1,:)',x(2,:)',[-4],"000")
181             
182         </scilab:image>
183     </refsection>
184     <refsection>
185         <title>Examples</title>
186         <programlisting role="example"><![CDATA[ 
187 // initialize state statistics (mean and er. variance)
188 m0=[10 10]';
189 p0=[100 0;0 100];
190 // create system
191 f=[1.1 50.1;0 0.8];
192 g=[1 0;0 1];
193 h=[1 0;0 1];
194 [hi,hj]=size(h);
195 // noise statistics
196 q=[.01 0;0 0.01];
197 r=20*eye(2,2);
198 // initialize system process
199 rand("seed",66);
200 rand("normal");
201 p0c=chol(p0);
202 x0=m0+p0c'*rand(ones(m0));
203 y=h*x0+chol(r)'*rand(ones(1:hi))';
204 yt=y;
205 // initialize plotted variables
206 x=x0;
207 // loop
208 ft=[f];
209 gt=[g];
210 ht=[h];
211 qt=[q];
212 rt=[r];
213 n=10;
214 for k=1:n
215     // generate the state and observation 
216     // at time k (i.e. xk and yk)
217     [x1,y]=system(x0,f,g,h,q,r);
218     x=[x x1];
219     yt=[yt y];
220     x0=x1;
221     ft=[ft f];
222     gt=[gt g];
223     ht=[ht h];
224     qt=[qt q];
225     rt=[rt r];
226 end
227 // get the wiener filter estimate
228 [xs,ps,xf,pf]=wiener(yt,m0,p0,ft,gt,ht,qt,rt);
229 // plot result
230 a=min([x(1,:)-2*sqrt(ps(1,1:2:2*(n+1))),xf(1,:),xs(1,:)]);
231 b=max([x(1,:)+2*sqrt(ps(1,1:2:2*(n+1))),xf(1,:),xs(1,:)]);
232 c=min([x(2,:)-2*sqrt(ps(2,2:2:2*(n+1))),xf(2,:),xs(2,:)]);
233 d=max([x(2,:)+2*sqrt(ps(2,2:2:2*(n+1))),xf(2,:),xs(2,:)]);
234 xmargin=max([abs(a),abs(b)]);
235 ymargin=max([abs(c),abs(d)]);
236 a=-0.1*xmargin+a;
237 b=.1*xmargin+b;
238 c=-0.1*ymargin+c;
239 d=.1*ymargin+d;
240 // plot frame, real state (x), and estimates (xf, and xs)
241 scf();
242 plot([a a b],[d c c]);
243 plot2d(x(1,:)',x(2,:)',[2],"000")
244 plot2d(xf(1,:)',xf(2,:)',[2],"000")
245 plot2d(xs(1,:)',xs(2,:)',[2],"000")
246 // mark data points (* for real data, o for estimates)
247 plot2d(xs(1,:)',xs(2,:)',[-2],"000")
248 plot2d(xf(1,:)',xf(2,:)',[-3],"000")
249 plot2d(x(1,:)',x(2,:)',[-4],"000")
250  ]]></programlisting>
251     </refsection>
252     
253 </refentry>