Mercurial > hg > ltpda

comparison m-toolbox/classes/@matrix/mchNoisegen.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 % MCHNOISEGEN Generates multichannel noise data series given a model
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % FUNCTION: mchNoisegen generates multichannel noise data series given a
5 % mutichannel noise generating filter.
6 %
7 % DESCRIPTION:
8 %
9 %
10 % CALL: out = mchNoisegen(mod, pl)
11 %
12 % INPUT:
13 % mod: is a matrix containing a multichannel noise generating
14 % filter aiming to generate colored noise from unitary variance
15 % independent white data series. Each element of the multichannel
16 % filter must be a parallel bank of filters as that produced by
17 % matrix/mchNoisegenFilter or ao/Noisegen2D. The filter matrix must
18 % be square.
19 %
20 % OUTPUT:
21 % out: is a matrix containg a multichannel colored noise time
22 % series which csd matrix is defined by mod'*mod.
23 %
24 % HISTORY: 19-10-2009 L Ferraioli
25 % Creation
26 %
27 % ------------------------------------------------------------------------
28 %
29 % <a href="matlab:utils.helper.displayMethodInfo('matrix', 'mchNoisegen')">Parameters Description</a>
30 %
31 % VERSION: $Id: mchNoisegen.m,v 1.6 2011/04/08 08:56:31 hewitson Exp $
32 %
33 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34 function varargout = mchNoisegen(varargin)
35
36 % Check if this is a call for parameters
37 if utils.helper.isinfocall(varargin{:})
38 varargout{1} = getInfo(varargin{3});
39 return
40 end
41
42 import utils.const.*
43 utils.helper.msg(msg.OMNAME, 'running %s/%s', mfilename('class'), mfilename);
44
45 % Collect input variable names
46 in_names = cell(size(varargin));
47 for ii = 1:nargin,in_names{ii} = inputname(ii);end
48
49 % Collect all ltpdauoh objects
50 [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
51 [pl, invars] = utils.helper.collect_objects(varargin(:), 'plist');
52
53 inhists = mtxs.hist;
54
55 % combine plists
56 pl = parse(pl, getDefaultPlist());
57 pl.getSetRandState();
58
59 outm = copy(mtxs,nargout);
60
61 % Get parameters and set params for fit
62 Nsecs = find(pl, 'nsecs');
63 fs = find(pl, 'fs');
64 yunit = find(pl, 'yunits');
65
66 % total number of data in the series
67 Ntot = round(Nsecs*fs);
68
69 % chose case for input filter
70 if isa(outm,'matrix') && isa(outm.objs(1),'filterbank')
71 % discrete system
72 sys = 'z';
73 mfil = outm.objs;
74 [nn,mm] = size(mfil);
75 if nn~=mm
76 error('!!! Filter matrix must be square')
77 end
78 % check if filter is a matrix built with noisegen2D
79 fromnsg2D = false;
80 if nn == 2
81 nn11 = numel(mfil(1,1).filters);
82 nn21 = numel(mfil(2,1).filters);
83 if nn11 == nn21
84 % get poles out of the filters
85 pl11 = zeros(nn11,1);
86 pl21 = zeros(nn21,1);
87 for ii=1:nn11
88 pl11(ii) = -1*mfil(1,1).filters(ii).b(2);
89 pl21(ii) = -1*mfil(2,1).filters(ii).b(2);
90 end
91 % check if poles are equal that means they are produced with
92 % noisegen2D
93 if all((pl11-pl21)<eps)
94 fromnsg2D = true;
95 end
96 end
97 end
98 elseif isa(outm,'matrix') && isa(outm.objs(1),'parfrac')
99 % continuous system
100 sys = 's';
101 mfil = outm.objs;
102 [nn,mm] = size(mfil);
103 if nn~=mm
104 error('!!! Filter matrix must be square')
105 end
106 else
107 error('!!! Input filter must be a ''matrix'' of ''filterbank'' or ''parfrac'' objects.')
108 end
109
110 % init output
111 out(nn,1) = ao;
112
113 % switch between input filter type
114 switch sys
115 case 'z' % discrete input filters
116
117 % init output data
118 o = zeros(nn,Ntot);
119
120 for zz=1:nn % moving along system dimension
121
122 % extract residues and poles from input objects
123 % run over matrix dimensions
124 res = [];
125 pls = [];
126 filtsz = [];
127 for jj=1:nn % collect filters coefficients along the columns zz
128 bfil = mfil(jj,zz).filters;
129 filtsz = [filtsz; numel(bfil)];
130 for kk=1:numel(bfil)
131 num = bfil(kk).a;
132 den = bfil(kk).b;
133 res = [res; num(1)];
134 pls = [pls; -1*den(2)];
135 end
136 end
137
138 % rescaling residues to get the correct result for univariate psd
139 res = res.*sqrt(fs/2);
140
141
142 % Nrs = numel(res);
143 % % get covariance matrix
144 % R = zeros(Nrs);
145 % for aa=1:Nrs
146 % for bb=1:Nrs
147 % R(aa,bb) = (res(aa)*conj(res(bb)))/(1-pls(aa)*conj(pls(bb)));
148 % end
149 % end
150 %
151 % % avoiding problems caused by roundoff errors
152 % HH = triu(R,0); % get upper triangular part of R
153 % HH1 = triu(R,1); % get upper triangular part of R above principal diagonal
154 % HH2 = HH1'; % do transpose conjugate
155 % R = HH + HH2; % reconstruct R in order to be really hermitian
156 %
157 % % get singular value decomposition of R
158 % [U,S,V] = svd(R,0);
159 %
160 % % conversion matrix
161 % A = V*sqrt(S);
162 %
163 % % generate unitary variance gaussian random noise
164 % %ns = randn(Nrs,Ntot);
165 % ns = randn(Nrs,1);
166 %
167 % % get correlated starting data points
168 % cns = A*ns;
169 %
170 % % need to correct for roundoff errors
171 % % cleaning up results for numerical approximations
172 % idx = imag(pls(:,1))==0;
173 % cns(idx) = real(cns(idx));
174 %
175 % % cleaning up results for numerical roundoff errors
176 % % states associated to complex conjugate poles must be complex
177 % % conjugate
178 % idxi = imag(pls(:,1))~=0;
179 % icns = cns(idxi);
180 % for jj = 1:2:numel(icns)
181 % icns(jj+1,1) = conj(icns(jj,1));
182 % end
183 % cns(idxi) = icns;
184
185 cns = getinitz(res,pls,filtsz,fromnsg2D);
186
187 y = zeros(sum(filtsz),2);
188 rnoise = zeros(sum(filtsz),1);
189 rns = randn(1,Ntot);
190 %rns = utils.math.blwhitenoise(Ntot,fs,1/Nsecs,fs/2);
191 %rns = rns.'; % willing to work with raw
192
193 y(:,1) = cns;
194
195 % start recurrence
196 for xx = 2:Ntot+1
197 rnoise(:,1) = rns(xx-1);
198 y(:,2) = pls.*y(:,1) + res.*rnoise;
199 idxr = 0;
200 for kk=1:nn
201 o(kk,xx-1) = o(kk,xx-1) + sum(y(idxr+1:idxr+filtsz(kk),2));
202 idxr = idxr+filtsz(kk);
203 end
204 y(:,1) = y(:,2);
205 end
206
207 end
208
209 clear rns
210
211 % build output ao
212 for dd=1:nn
213 out(dd,1) = ao(tsdata(o(dd,:),fs));
214 out(dd,1).setYunits(unit(yunit));
215 end
216
217 outm = matrix(out);
218
219 case 's' % continuous input filters
220
221 o = zeros(nn,Ntot);
222
223 T = 1/fs; % sampling period
224
225 for zz=1:nn % moving along system dimension
226
227 % extract residues and poles from input objects
228 % run over matrix dimensions
229 res = [];
230 pls = [];
231 filtsz = [];
232 for jj=1:nn % collect filters coefficients along the columns zz
233 bfil = mfil(jj,zz);
234 filtsz = [filtsz; numel(bfil.res)];
235 res = [res; bfil.res];
236 pls = [pls; bfil.poles];
237 end
238
239 % rescaling residues to get the correct result for univariate psd
240 res = res.*sqrt(fs/2);
241
242 Nrs = numel(res);
243
244 % get covariance matrix for innovation
245 Ri = zeros(Nrs);
246 for aa=1:Nrs
247 for bb=1:Nrs
248 Ri(aa,bb) = (res(aa)*conj(res(bb)))*(exp((pls(aa) + conj(pls(bb)))*T)-1)/(pls(aa) + conj(pls(bb)));
249 end
250 end
251
252 % avoiding problems caused by roundoff errors
253 HH = triu(Ri,0); % get upper triangular part of R
254 HH1 = triu(Ri,1); % get upper triangular part of R above principal diagonal
255 HH2 = HH1'; % do transpose conjugate
256 Ri = HH + HH2; % reconstruct R in order to be really hermitian
257
258 % get singular value decomposition of R
259 [Ui,Si,Vi] = svd(Ri,0);
260
261 % conversion matrix for innovation
262 Ai = Vi*sqrt(Si);
263
264 % get covariance matrix for initial state
265 Rx = zeros(Nrs);
266 for aa=1:Nrs
267 for bb=1:Nrs
268 Rx(aa,bb) = -1*(res(aa)*conj(res(bb)))/(pls(aa) + conj(pls(bb)));
269 end
270 end
271
272 % avoiding problems caused by roundoff errors
273 HH = triu(Rx,0); % get upper triangular part of R
274 HH1 = triu(Rx,1); % get upper triangular part of R above principal diagonal
275 HH2 = HH1'; % do transpose conjugate
276 Rx = HH + HH2; % reconstruct R in order to be really hermitian
277
278 % get singular value decomposition of R
279 [Ux,Sx,Vx] = svd(Rx,0);
280
281 % conversion matrix for initial state
282 Ax = Vx*sqrt(Sx);
283
284 % generate unitary variance gaussian random noise
285 %ns = randn(Nrs,Ntot);
286 ns = randn(Nrs,1);
287
288 % get correlated starting data points
289 cns = Ax*ns;
290
291 % need to correct for roundoff errors
292 % cleaning up results for numerical approximations
293 idx = imag(pls(:,1))==0;
294 cns(idx) = real(cns(idx));
295
296 % cleaning up results for numerical roundoff errors
297 % states associated to complex conjugate poles must be complex
298 % conjugate
299 idxi = imag(pls(:,1))~=0;
300 icns = cns(idxi);
301 for jj = 1:2:numel(icns)
302 icns(jj+1,1) = conj(icns(jj,1));
303 end
304 cns(idxi) = icns;
305
306 y = zeros(sum(filtsz),2);
307 rnoise = zeros(sum(filtsz),1);
308 rns = randn(1,Ntot);
309 %rns = utils.math.blwhitenoise(Ntot,fs,1/Nsecs,fs/2);
310 %rns = rns.'; % willing to work with raw
311
312 y(:,1) = cns;
313
314 % start recurrence
315 for xx = 2:Ntot+1
316 % innov = Ai*randn(sum(filtsz),1);
317 rnoise(:,1) = rns(xx-1);
318 innov = Ai*rnoise;
319 % need to correct for roundoff errors
320 % cleaning up results for numerical approximations
321 innov(idx) = real(innov(idx));
322
323 % cleaning up results for numerical roundoff errors
324 % states associated to complex conjugate poles must be complex
325 % conjugate
326 iinnov = innov(idxi);
327 for jj = 1:2:numel(iinnov)
328 iinnov(jj+1,1) = conj(iinnov(jj,1));
329 end
330 innov(idxi) = iinnov;
331
332 y(:,2) = diag(exp(pls.*T))*y(:,1) + innov;
333
334 idxr = 0;
335 for kk=1:nn
336 o(kk,xx-1) = o(kk,xx-1) + sum(y(idxr+1:idxr+filtsz(kk),2));
337 idxr = idxr+filtsz(kk);
338 end
339 y(:,1) = y(:,2);
340 end
341
342 end
343
344 % clear rns
345
346 % build output ao
347 for dd=1:nn
348 out(dd,1) = ao(tsdata(o(dd,:),fs));
349 out(dd,1).setYunits(unit(yunit));
350 end
351
352 outm = matrix(out);
353
354 end
355
356 % set output
357 varargout{1} = outm;
358
359
360 end
361
362 %--------------------------------------------------------------------------
363 % Get Info Object
364 %--------------------------------------------------------------------------
365 function ii = getInfo(varargin)
366 if nargin == 1 && strcmpi(varargin{1}, 'None')
367 sets = {};
368 pl = [];
369 else
370 sets = {'Default'};
371 pl = getDefaultPlist;
372 end
373 % Build info object
374 ii = minfo(mfilename, 'matrix', 'ltpda', utils.const.categories.sigproc, '$Id: mchNoisegen.m,v 1.6 2011/04/08 08:56:31 hewitson Exp $', sets, pl);
375 ii.setArgsmin(1);
376 ii.setOutmin(1);
377 ii.setOutmax(1);
378 end
379
380 %--------------------------------------------------------------------------
381 % Get Default Plist
382 %--------------------------------------------------------------------------
383 function plout = getDefaultPlist()
384 persistent pl;
385 if exist('pl', 'var')==0 || isempty(pl)
386 pl = buildplist();
387 end
388 plout = pl;
389 end
390
391 function pl = buildplist()
392
393 pl = plist();
394
395 % Nsecs
396 p = param({'nsecs', 'Number of seconds in the desired noise data series.'}, {1, {1}, paramValue.OPTIONAL});
397 pl.append(p);
398
399 % Fs
400 p = param({'fs', 'The sampling frequency of the noise data series.'}, {1, {1}, paramValue.OPTIONAL});
401 pl.append(p);
402
403 % Yunits
404 p = param({'yunits','Unit on Y axis.'}, paramValue.STRING_VALUE(''));
405 pl.append(p);
406
407 end
408
409 %--------------------------------------------------------------------------
410 % Local function
411 % Estimate init values for the z domain case
412 %--------------------------------------------------------------------------
413 function cns = getinitz(res,pls,filtsz,fromnsg2D)
414
415 if fromnsg2D % init in the case of filters produced with noisegen2D
416 % divide in 2 the problem
417 res1 = res(1:filtsz(1));
418 res2 = res(filtsz(1)+1:filtsz(2));
419 pls1 = pls(1:filtsz(1));
420 pls2 = pls(filtsz(1)+1:filtsz(2));
421
422 %%% the first series %%%
423 Nrs = numel(res1);
424 % get covariance matrix
425 R = zeros(Nrs);
426 for aa=1:Nrs
427 for bb=1:Nrs
428 R(aa,bb) = (res1(aa)*conj(res1(bb)))/(1-pls1(aa)*conj(pls1(bb)));
429 end
430 end
431
432 % avoiding problems caused by roundoff errors
433 HH = triu(R,0); % get upper triangular part of R
434 HH1 = triu(R,1); % get upper triangular part of R above principal diagonal
435 HH2 = HH1'; % do transpose conjugate
436 R = HH + HH2; % reconstruct R in order to be really hermitian
437
438 % get singular value decomposition of R
439 [U,S,V] = svd(R,0);
440
441 % conversion matrix
442 A = V*sqrt(S);
443
444 % generate unitary variance gaussian random noise
445 %ns = randn(Nrs,Ntot);
446 ns = randn(Nrs,1);
447
448 % get correlated starting data points
449 cns1 = A*ns;
450
451 % need to correct for roundoff errors
452 % cleaning up results for numerical approximations
453 idx = imag(pls1(:,1))==0;
454 cns1(idx) = real(cns1(idx));
455
456 % cleaning up results for numerical roundoff errors
457 % states associated to complex conjugate poles must be complex
458 % conjugate
459 idxi = imag(pls1(:,1))~=0;
460 icns1 = cns1(idxi);
461 for jj = 1:2:numel(icns1)
462 icns1(jj+1,1) = conj(icns1(jj,1));
463 end
464 cns1(idxi) = icns1;
465
466 %%% the second series %%%
467 Nrs = numel(res2);
468 % get covariance matrix
469 R = zeros(Nrs);
470 for aa=1:Nrs
471 for bb=1:Nrs
472 R(aa,bb) = (res2(aa)*conj(res2(bb)))/(1-pls2(aa)*conj(pls2(bb)));
473 end
474 end
475
476 % avoiding problems caused by roundoff errors
477 HH = triu(R,0); % get upper triangular part of R
478 HH1 = triu(R,1); % get upper triangular part of R above principal diagonal
479 HH2 = HH1'; % do transpose conjugate
480 R = HH + HH2; % reconstruct R in order to be really hermitian
481
482 % get singular value decomposition of R
483 [U,S,V] = svd(R,0);
484
485 % conversion matrix
486 A = V*sqrt(S);
487
488 % generate unitary variance gaussian random noise
489 %ns = randn(Nrs,Ntot);
490 ns = randn(Nrs,1);
491
492 % get correlated starting data points
493 cns2 = A*ns;
494
495 % need to correct for roundoff errors
496 % cleaning up results for numerical approximations
497 idx = imag(pls2(:,1))==0;
498 cns2(idx) = real(cns2(idx));
499
500 % cleaning up results for numerical roundoff errors
501 % states associated to complex conjugate poles must be complex
502 % conjugate
503 idxi = imag(pls2(:,1))~=0;
504 icns2 = cns2(idxi);
505 for jj = 1:2:numel(icns2)
506 icns2(jj+1,1) = conj(icns2(jj,1));
507 end
508 cns2(idxi) = icns2;
509
510 %%% combine results %%%
511 cns = [cns1;cns2];
512
513 else % init in the case of filters produced with matrix constructor
514
515 Nrs = numel(res);
516 % get covariance matrix
517 R = zeros(Nrs);
518 for aa=1:Nrs
519 for bb=1:Nrs
520 R(aa,bb) = (res(aa)*conj(res(bb)))/(1-pls(aa)*conj(pls(bb)));
521 end
522 end
523
524 % avoiding problems caused by roundoff errors
525 HH = triu(R,0); % get upper triangular part of R
526 HH1 = triu(R,1); % get upper triangular part of R above principal diagonal
527 HH2 = HH1'; % do transpose conjugate
528 R = HH + HH2; % reconstruct R in order to be really hermitian
529
530 % get singular value decomposition of R
531 [U,S,V] = svd(R,0);
532
533 % conversion matrix
534 A = V*sqrt(S);
535
536 % generate unitary variance gaussian random noise
537 %ns = randn(Nrs,Ntot);
538 ns = randn(Nrs,1);
539
540 % get correlated starting data points
541 cns = A*ns;
542
543 % need to correct for roundoff errors
544 % cleaning up results for numerical approximations
545 idx = imag(pls(:,1))==0;
546 cns(idx) = real(cns(idx));
547
548 % cleaning up results for numerical roundoff errors
549 % states associated to complex conjugate poles must be complex
550 % conjugate
551 idxi = imag(pls(:,1))~=0;
552 icns = cns(idxi);
553 for jj = 1:2:numel(icns)
554 icns(jj+1,1) = conj(icns(jj,1));
555 end
556 cns(idxi) = icns;
557
558 end
559
560 end