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]);