0001 function varargout = ltpda_simulate_discrete(varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
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
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
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
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,{}))
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
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
0100 nsamples = 0;
0101 end
0102
0103
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
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
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
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
0158
0159 hi_new = history('simulate_discrete','version', syst);
0160
0161
0162 h = history(hi_new);
0163
0164
0165 varargout = ao(ts,h) ;
0166 end