Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/getinitstate.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % GETINITSTATE
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION:
5 %
6 % Initialize filters for noise generation
7 %
8 %
9 % CALL:
10 %
11 % INPUT:
12 %
13 % - res
14 % - poles
15 % - S0
16 %
17 % OUTPUT:
18 %
19 % Zi
20 %
21 % REFERENCES:
22 %
23 %
24 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
25 % HISTORY: 23-04-2009 L Ferraioli
26 % Creation
27 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28 %
29 % VERSION: $Id: ndeigcsd.m,v 1.2 2009/06/10 16:15:32 luigi Exp $
30 %
31 %
32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 function Zi = getinitstate(res,poles,S0,varargin)
34
35 % default value for the method used
36 mtd = 'svd';
37
38 if ~isempty(varargin)
39 for j=1:length(varargin)
40 if strcmpi(varargin{j},'mtd')
41 mtd = lower(varargin{j+1});
42 end
43 end
44 end
45
46 % get problem dimesionality, it is assumed that res for a given filter are
47 % input as columns and that each different filter has the same number of
48 % residues and poles
49 [rw,cl] = size(res);
50
51 N = numel(res(:,1));
52
53
54 % check size
55 if cl == 1 % one dimensional
56 A = res(:);
57 p = poles(:);
58
59 % Marking complex and real poles
60 % cindex = 1; pole is complex, next conjugate pole is marked with cindex
61 % = 2. cindex = 0; pole is real
62 cindex=zeros(numel(p),1);
63 for mm=1:numel(p)
64 if imag(p(mm))~=0
65 if mm==1
66 cindex(mm)=1;
67 else
68 if cindex(mm-1)==0 || cindex(mm-1)==2
69 cindex(mm)=1; cindex(mm+1)=2;
70 else
71 cindex(mm)=2;
72 end
73 end
74 end
75 end
76
77 NA = numel(A);
78
79 % Build covariance matrix for filter states
80 H = zeros(NA);
81 for aa = 1:NA
82 for bb = 1:NA
83 H(aa,bb) = (p(aa)*conj(p(bb))*A(aa)*conj(A(bb))*S0)/(1-p(aa)*conj(p(bb)));
84 end
85 end
86
87 % avoiding problems caused by roundoff errors
88 HH = triu(H,0); % get upper triangular part of H
89 HH1 = triu(H,1); % get upper triangular part of H above principal diagonal
90 HH2 = HH1'; % do transpose conjugate
91 H = HH + HH2; % reconstruct H in order to be really hermitian
92
93 % switch between methods
94 switch mtd
95 case 'svd'
96 % get decomposition of Hr
97 [U,S,V] = svd(H,0);
98 % output initial states for gaussian data series
99 ZZ = V*(sqrt(diag(S)).*randn(N,1));
100
101 % cleaning up results for numerical approximations
102 idx = imag(poles(:,1))==0;
103 ZZ(idx) = real(ZZ(idx));
104
105 % cleaning up results for numerical roundoff errors
106 % states associated to complex conjugate poles must be complex
107 % conjugate
108 for jj = 1:numel(ZZ)
109 if cindex(jj)==1
110 ZZ(jj+1) = conj(ZZ(jj));
111 end
112 end
113
114 case 'mvnorm'
115
116 ZZ = mvnrnd(zeros(N,1),H,1);
117 % willing to work with columns
118 if size(ZZ,2)>1
119 ZZ = ZZ.';
120 end
121 if imag(ZZ(1))~=0 && imag(p(1))==0
122 % flip
123 ZZ = flipud(ZZ);
124 end
125
126 % cleaning up results for numerical roundoff errors
127 % states associated to complex conjugate poles must be complex
128 % conjugate
129 for jj = 1:numel(ZZ)
130 if cindex(jj)==1
131 ZZ(jj+1) = conj(ZZ(jj));
132 end
133 end
134
135 end
136
137
138
139 else % cl dmensional
140
141 % join residues and poles
142 A = [];
143 p = [];
144 pdim = [];
145
146 for ii=1:cl
147 % Join poles and residues as a single column
148 A = [A; res(:,ii)];
149 p = [p; poles(:,ii)];
150 pdim = [pdim; numel(poles(:,ii))];
151 end
152
153 % Marking complex and real poles
154 % cindex = 1; pole is complex, next conjugate pole is marked with cindex
155 % = 2. cindex = 0; pole is real
156 cindex=zeros(numel(p),1);
157 for mm=1:numel(p)
158 if imag(p(mm))~=0
159 if mm==1
160 cindex(mm)=1;
161 else
162 if cindex(mm-1)==0 || cindex(mm-1)==2
163 cindex(mm)=1; cindex(mm+1)=2;
164 else
165 cindex(mm)=2;
166 end
167 end
168 end
169 end
170
171 % sanity check, search if poles are equals
172 eqpoles = false;
173 if ~all(logical(diff(pdim))) % is executed if the elements of pdim are equal
174 % compare poles series
175 for ii=2:cl
176 if all((poles(:,ii)-poles(:,ii-1))<eps)
177 eqpoles = true;
178 end
179 end
180 end
181
182 if ~eqpoles % poles are different
183
184 NA = numel(A);
185
186 % Build covariance matrix for filter states
187 H = zeros(NA);
188 for aa = 1:NA
189 for bb = 1:NA
190 H(aa,bb) = (p(aa)*conj(p(bb))*A(aa)*conj(A(bb))*S0)/(1-p(aa)*conj(p(bb)));
191 end
192 end
193
194 % avoiding problems caused by roundoff errors
195 HH = triu(H,0); % get upper triangular part of H
196 HH1 = triu(H,1); % get upper triangular part of H above principal diagonal
197 HH2 = HH1'; % do transpose conjugate
198 H = HH + HH2; % reconstruct H in order to be really hermitian
199
200 % get full rank H
201 [U,S,V] = svd(H,0);
202
203 % reducing size
204 Ur = U(1:N,1:N);
205 Sr = S(1:N,1:N);
206 Vr = V(1:N,1:N);
207
208 % New full rank covariance
209 Hr = Vr*Sr*Vr';
210
211 % avoiding problems caused by roundoff errors
212 HH = triu(Hr,0); % get upper triangular part of H
213 HH1 = triu(Hr,1); % get upper triangular part of H above principal diagonal
214 HH2 = HH1'; % do transpose conjugate
215 Hr = HH + HH2; % reconstruct H in order to be really hermitian
216
217 % switch between methods
218 switch mtd
219 case 'svd'
220 % get decomposition of Hr
221 [UU,SS,VV] = svd(Hr,0);
222 % output initial states for gaussian data series
223 ZZ = VV*(sqrt(diag(SS)).*randn(N,1));
224
225 % cleaning up results for numerical roundoff errors
226 idx = imag(p)==0;
227 ZZ(idx) = real(ZZ(idx));
228
229 % cleaning up results for numerical roundoff errors
230 % states associated to complex conjugate poles must be complex
231 % conjugate
232 for jj = 1:numel(ZZ)
233 if cindex(jj)==1
234 ZZ(jj+1) = conj(ZZ(jj));
235 end
236 end
237
238
239 case 'mvnorm'
240
241 ZZ = mvnrnd(zeros(N,1),Hr,1);
242 % willing to work with columns
243 if size(ZZ,2)>1
244 ZZ = ZZ.';
245 end
246 if imag(ZZ(1))~=0 && imag(p(1))==0
247 % flip
248 ZZ = flipud(ZZ);
249 end
250
251 % cleaning up results for numerical roundoff errors
252 % states associated to complex conjugate poles must be complex
253 % conjugate
254 for jj = 1:numel(ZZ)
255 if cindex(jj)==1
256 ZZ(jj+1) = conj(ZZ(jj));
257 end
258 end
259
260 end
261
262 else % poles are in common
263
264 NA = numel(A);
265
266 % Build block diagonal covariance matrix for filter states
267 H = zeros(NA);
268 for ii = 1:numel(pdim)
269 for aa = 1+(ii-1)*pdim(ii):(ii-1)*pdim(ii)+pdim(ii)
270 for bb = 1+(ii-1)*pdim(ii):(ii-1)*pdim(ii)+pdim(ii)
271 H(aa,bb) = (p(aa)*conj(p(bb))*A(aa)*conj(A(bb))*S0)/(1-p(aa)*conj(p(bb)));
272 end
273 end
274 end
275
276 % avoiding problems caused by roundoff errors
277 HH = triu(H,0); % get upper triangular part of H
278 HH1 = triu(H,1); % get upper triangular part of H above principal diagonal
279 HH2 = HH1'; % do transpose conjugate
280 H = HH + HH2; % reconstruct H in order to be really hermitian
281
282 % switch between methods
283 switch mtd
284 case 'svd'
285 % get decomposition of Hr
286 [UU,SS,VV] = svd(H,0);
287 % output initial states for gaussian data series
288 rd = randn(N,1);
289 rdv = [];
290 for jj = 1:numel(pdim);
291 rdv = [rdv; rd];
292 end
293 ZZ = VV*(sqrt(diag(SS)).*rdv);
294
295 % cleaning up results for numerical roundoff errors
296 idx = imag(p)==0;
297 ZZ(idx) = real(ZZ(idx));
298
299 % cleaning up results for numerical roundoff errors
300 % states associated to complex conjugate poles must be complex
301 % conjugate
302 for jj = 1:numel(ZZ)
303 if cindex(jj)==1
304 ZZ(jj+1) = conj(ZZ(jj));
305 end
306 end
307
308
309 case 'mvnorm'
310
311 ZZ = mvnrnd(zeros(N,1),H,1);
312 % willing to work with columns
313 if size(ZZ,2)>1
314 ZZ = ZZ.';
315 end
316 if imag(ZZ(1))~=0 && imag(p(1))==0
317 % flip
318 ZZ = flipud(ZZ);
319 end
320
321 % cleaning up results for numerical roundoff errors
322 % states associated to complex conjugate poles must be complex
323 % conjugate
324 for jj = 1:numel(ZZ)
325 if cindex(jj)==1
326 ZZ(jj+1) = conj(ZZ(jj));
327 end
328 end
329
330 end
331
332 end
333
334 end
335
336 Zi = ZZ;
337
338 end