Home > m > timetools > statespacefunctions > ltpda_simulate_discrete.m

ltpda_simulate_discrete

PURPOSE ^

simulates a discrete state space system

SYNOPSIS ^

function varargout = ltpda_simulate_discrete(varargin)

DESCRIPTION ^

simulates a discrete state space system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: simulates a discrete state space system

 CALL: ao_out = simulate_discrete(ao_in, syst)

 INPUTS: ao_in_list - plist of :
                         - ao containing concatenated data stream
                         - plist describing an input noise with a mean 
                             'MEAN' (vector) and a sigma 'SIGMA' (Matrix)
 syst - system data in plist format

 OUTPUTS: ao_out _ output analysis object

 ***** THERE ARE NO DEFAULT PARAMETERS *****

 &&VERSION: $Id: ltpda_simulate_discrete.html,v 1.1 2008/03/01 12:29:32 hewitson Exp $

 HISTORY: 19-02-2008 A Grynagier
 17-02-2008 A Grynagier
 Creation 30-01-2008 M Weyrich

 TO DO: options to be defined
 implement input type of three kinds:
 - AO
 - Constant 
 - gaussian white noise
 use LTPDA functions to generate white noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function varargout = ltpda_simulate_discrete(varargin)
0002 %simulates a discrete state space system
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % DESCRIPTION: simulates a discrete state space system
0006 %
0007 % CALL: ao_out = simulate_discrete(ao_in, syst)
0008 %
0009 % INPUTS: ao_in_list - plist of :
0010 %                         - ao containing concatenated data stream
0011 %                         - plist describing an input noise with a mean
0012 %                             'MEAN' (vector) and a sigma 'SIGMA' (Matrix)
0013 % syst - system data in plist format
0014 %
0015 % OUTPUTS: ao_out _ output analysis object
0016 %
0017 % ***** THERE ARE NO DEFAULT PARAMETERS *****
0018 %
0019 % &&VERSION: $Id: ltpda_simulate_discrete.html,v 1.1 2008/03/01 12:29:32 hewitson Exp $
0020 %
0021 % HISTORY: 19-02-2008 A Grynagier
0022 % 17-02-2008 A Grynagier
0023 % Creation 30-01-2008 M Weyrich
0024 %
0025 % TO DO: options to be defined
0026 % implement input type of three kinds:
0027 % - AO
0028 % - Constant
0029 % - gaussian white noise
0030 % use LTPDA functions to generate white noise
0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0032 
0033 ALGONAME = mfilename;
0034 VERSION = '$Id: ltpda_simulate_discrete.html,v 1.1 2008/03/01 12:29:32 hewitson Exp $';
0035 CATEGORY = 'STATESPACE';
0036 display(['starting ' ALGONAME]);
0037 
0038 if isequal( varargin{1}, 'Version')
0039     varargout = VERSION;
0040     return;
0041 elseif isequal(varargin{1}, 'Params')
0042     varargout = plist();
0043     return;
0044 elseif isequal(varargin{1}, 'Category')
0045     varargout = CATEGORY;
0046     return;
0047 end
0048 syst = varargin{2};
0049 ao_in_list = varargin{1};
0050 
0051 %% retrieve system infos
0052 A = find(syst, 'AMAT');
0053 B = find(syst, 'BMATS');
0054 C = find(syst, 'CMAT');
0055 D = find(syst, 'DMATS');
0056 A = A{1};
0057 C = C{1};
0058 nss = size(A,1);
0059 %check iniside A, Bs, c, Ds there is no symbolic variable
0060 Xini = find(syst, 'XINI');
0061 timestep = find(syst, 'TIMESTEP');
0062 NBINPUTS = find(syst, 'NBINPUTS');
0063 INPUTISUSED = find(syst, 'INPUTISUSED');
0064 INPUTNAMES = find(syst, 'INPUTNAMES');
0065 INPUTSIZES = find(syst, 'INPUTSIZES');
0066 XISOUTPUT = find(syst, 'XISOUTPUT');
0067 YISOUTPUT = find(syst, 'YISOUTPUT');
0068 
0069 
0070 
0071 %% read out AO (time-series data)
0072 ao_in = cell(NBINPUTS,1);
0073 for i=1:NBINPUTS
0074     if INPUTISUSED(i) 
0075         ao_in{i} = find(ao_in_list, INPUTNAMES{i});
0076         if isequal(class(ao_in{i}), 'ao')
0077             if not(isequal(ao_in{i}.data.x,{})) %making sure ao is not empty
0078                 t = ao_in{i}.data.x;
0079                 dt = diff(t);
0080                 if sum(not(isequal(dt,timestep)))>1
0081                     msg = ['error in input number ', i , ' called ', INPUTNAMES{i}, 'because wrong timestep for simulation'];
0082                     error(msg);
0083                 end
0084                 if not( isequal( size(B{i},1), size(D{i},1), size(ao_in{i}.data.y,1) ) )
0085                     msg = ['error in input number ', i , ' called ', INPUTNAMES{i}, 'because wrong input size'];
0086                     error(msg);
0087                 end
0088             end
0089         end
0090     end
0091 end
0092 %% getting correct number of samples
0093 nsamples = inf;
0094 for i=1:NBINPUTS
0095     if isequal(class(ao_in{i}), 'ao')
0096         nsamples = min(nsamples,size(ao_in{i}.data.y,1));
0097     end
0098 end
0099 if nsamples == inf % case there is no input!
0100     nsamples = 0;
0101 end
0102 
0103 %% generating random data
0104 for i=1:NBINPUTS
0105     if isequal(class(ao_in{i}), 'plist')
0106         sigma = find(ao_in{i}, 'SIGMA');
0107         mean = find(ao_in{i}, 'MEAN');
0108         if sigma==0
0109             ao_in{i} = ao(ones(nsamples,INPUTSIZES(i))*diag(mean));
0110         else
0111             ao_in{i} = ao(ones(nsamples,INPUTSIZES(i))*diag(mean) + randn(nsamples,INPUTSIZES(i))*sigma);
0112         end
0113     end
0114 end
0115 
0116 %% simulation first step
0117 x = zeros(nss, nsamples); 
0118 y = zeros(size(C,1), nsamples);
0119 x(:,1) = Xini;
0120 IncrementY = zeros(size(C,1),1);
0121 for i=1:NBINPUTS
0122     if INPUTISUSED(i)
0123         IncrementY = IncrementY +  D{i}* ao_in{i}.data.y(1,:)' ;
0124     end
0125 end
0126 y(:,1) = C* x(:,1) +IncrementY;
0127 
0128 %% simulation loop
0129 for k = 2:nsamples
0130     IncrementX = zeros(nss,1);
0131     IncrementY = zeros(size(C,1),1);
0132     for i=1:NBINPUTS
0133         if INPUTISUSED(i)
0134             IncrementX = IncrementX +  B{i}*ao_in{i}.data.y(k,:)' ;
0135             IncrementY = IncrementY +  D{i}*ao_in{i}.data.y(k,:)' ;
0136         end
0137     end
0138     x(:,k) = A* x(:,k-1) + IncrementX ;
0139     y(:,k) = C* x(:,k) + IncrementY ;
0140 end
0141 
0142 %% data object
0143 out_dat = zeros(0, nsamples);
0144 if XISOUTPUT
0145     out_dat = [out_dat; x];
0146 end
0147 if YISOUTPUT
0148     out_dat = [out_dat; y];
0149 end
0150 
0151 if exist('t','var')
0152     ts = tsdata(t,out_dat) ;
0153 else
0154     ts = tsdata(1:nsamples,out_dat) ;
0155 end
0156 
0157 %% history object, old history from ao_in plus new
0158 % hi_old = ao_in.hist ;
0159 hi_new = history('simulate_discrete','version', syst);
0160 % hi = [hi_old hi_new];
0161 % h = history(hi);
0162 h = history(hi_new);
0163 
0164 %% construct output analysis object
0165 varargout  = ao(ts,h) ;
0166 end

Generated on Tue 26-Feb-2008 10:52:52 by m2html © 2003