add a unit test for lindquist based on the commit 909a8330d3fb367b12653522556ef26606f...
[scilab.git] / scilab / modules / signal_processing / tests / unit_tests / lindquist.tst
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) 2013 - Scilab Enterprises - Charlotte Hecquet
4 // Copyright (C) 2013 - Scilab Enterprises - Sylvestre Ledru
5 //
6 //  This file is distributed under the same license as the Scilab package.
7 // =============================================================================
8
9 //Generate signal
10 x=%pi/10:%pi/10:102.4*%pi;
11 y=[1; 1] * sin(x) + [sin(2*x); sin(1.9*x)];
12
13 //Compute correlations
14 c=[];
15 for j=1:2
16    for k=1:2
17      c=[c;corr(y(k,:),y(j,:),64)];
18    end
19 end
20 c=matrix(c,2,128);
21
22 //Find H,F,G with 6 states
23 hk=hank(20,20,c);
24 [H,F,G]=phc(hk,2,6);
25
26 //Solve Riccati equation
27 R0=c(1:2,1:2);
28 [P,R,T]=lindquist(100,H,F,G,R0);
29
30 assert_checkalmostequal(P,[52.953064387089377,-170.091664767326478;390.114441205931939,-1126.45220755477089;-2285.86604463640651,4579.41083686864113;1084.73413409944055,-1564.67420653525346;-5116.06921833499018,8743.68286689257002;-1976.40587294372335,3503.44251127613916]);
31
32 assert_checkalmostequal([-704.208398385009673,1053.5647942803439;1053.5647942827934,-1981.29618548816984], R);
33
34 assert_checkalmostequal(T, [-0.72821126450254769,-0.45461826616982365;-1.28940233637036883,-1.25076781365101763;0.25073933783035585,2.01963206598920664;3.13684440927242791,0.29283797163618436;-1.64165514197064755,1.68395444469620714;-1.61905018245414611,3.26102925315640091]);