ltpda: m-toolbox/classes/@matrix/MultiChannelNoise.m annotate

author	Daniele Nicolodi <nicolodi@science.unitn.it>
date	Wed, 23 Nov 2011 19:22:13 +0100
parents
children

rev	line source
0 f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	1 % MULTICHANNELNOISE
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	3 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	4 % FUNCTION: MultiChannelNoise
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	5 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	6 % DESCRIPTION:
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	7 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	8 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	9 % CALL: a = MultiChannelNoise(a, pl)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	10 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	11 % PARAMETER:
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	12 % pl: plist containing Nsecs, fs
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	13 %
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	15 function a = MultiChannelNoise(a, pli)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	16
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	17 VERSION = '$Id: MultiChannelNoise.m,v 1.10 2010/10/29 16:09:13 ingo Exp $';
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	18
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	19 % get AO info
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	20 iobj = matrix.getInfo('matrix', 'Multichannel Noise');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	21
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	22 % Set the method version string in the minfo object
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	23 iobj.setMversion([VERSION '-->' iobj.mversion]);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	24
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	25 % Add default values
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	26 pl = parse(pli, iobj.plists);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	27
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	28 % Get parameters and set params for fit
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	29 Nsecs = find(pl, 'nsecs');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	30 fs = find(pl, 'fs');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	31 filter = find(pl, 'model');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	32 yunit = find(pl, 'yunits');
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	33
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	34 % total number of data in the series
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	35 Ntot = round(Nsecs*fs);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	36
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	37 % chose case for input filter
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	38 if isa(filter,'matrix') && isa(filter.objs(1),'filterbank')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	39 % discrete system
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	40 sys = 'z';
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	41 mfil = filter.objs;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	42 [nn,mm] = size(mfil);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	43 if nn~=mm
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	44 error('!!! Filter matrix must be square')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	45 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	46 elseif isa(filter,'matrix') && isa(filter.objs(1),'parfrac')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	47 % continuous system
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	48 sys = 's';
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	49 mfil = filter.objs;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	50 [nn,mm] = size(mfil);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	51 if nn~=mm
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	52 error('!!! Filter matrix must be square')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	53 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	54 else
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	55 error('!!! Input filter must be a ''matrix'' of ''filterbank'' or ''parfrac'' objects.')
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	56 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	57
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	58 % init output
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	59 out(nn,1) = ao;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	60
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	61 % switch between input filter type
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	62 switch sys
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	63 case 'z' % discrete input filters
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	64
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	65 o = zeros(nn,Ntot);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	66
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	67 for zz=1:nn % moving along system dimension
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	68
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	69 % extract residues and poles from input objects
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	70 % run over matrix dimensions
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	71 res = [];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	72 pls = [];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	73 filtsz = [];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	74 for jj=1:nn % collect filters coefficients along the columns zz
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	75 bfil = mfil(jj,zz).filters;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	76 filtsz = [filtsz; numel(bfil)];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	77 for kk=1:numel(bfil)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	78 num = bfil(kk).a;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	79 den = bfil(kk).b;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	80 res = [res; num(1)];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	81 pls = [pls; -1*den(2)];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	82 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	83 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	84
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	85 % rescaling residues to get the correct result for univariate psd
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	86 res = res.*sqrt(fs/2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	87
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	88 Nrs = numel(res);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	89 % get covariance matrix
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	90 R = zeros(Nrs);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	91 for aa=1:Nrs
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	92 for bb=1:Nrs
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	93 R(aa,bb) = (res(aa)conj(res(bb)))/(1-pls(aa)conj(pls(bb)));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	94 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	95 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	96
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	97 % avoiding problems caused by roundoff errors
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	98 HH = triu(R,0); % get upper triangular part of R
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	99 HH1 = triu(R,1); % get upper triangular part of R above principal diagonal
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	100 HH2 = HH1'; % do transpose conjugate
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	101 R = HH + HH2; % reconstruct R in order to be really hermitian
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	102
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	103 % get singular value decomposition of R
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	104 [U,S,V] = svd(R,0);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	105
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	106 % conversion matrix
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	107 A = V*sqrt(S);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	108
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	109 % generate unitary variance gaussian random noise
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	110 %ns = randn(Nrs,Ntot);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	111 ns = randn(Nrs,1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	112
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	113 % get correlated starting data points
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	114 cns = A*ns;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	115
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	116 % need to correct for roundoff errors
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	117 % cleaning up results for numerical approximations
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	118 idx = imag(pls(:,1))==0;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	119 cns(idx) = real(cns(idx));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	120
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	121 % cleaning up results for numerical roundoff errors
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	122 % states associated to complex conjugate poles must be complex
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	123 % conjugate
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	124 idxi = imag(pls(:,1))~=0;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	125 icns = cns(idxi);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	126 for jj = 1:2:numel(icns)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	127 icns(jj+1,1) = conj(icns(jj,1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	128 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	129 cns(idxi) = icns;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	130
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	131 y = zeros(sum(filtsz),2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	132 rnoise = zeros(sum(filtsz),1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	133 rns = randn(1,Ntot);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	134 %rns = utils.math.blwhitenoise(Ntot,fs,1/Nsecs,fs/2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	135 %rns = rns.'; % willing to work with raw
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	136
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	137 y(:,1) = cns;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	138
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	139 % start recurrence
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	140 for xx = 2:Ntot+1
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	141 rnoise(:,1) = rns(xx-1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	142 y(:,2) = pls.y(:,1) + res.rnoise;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	143 idxr = 0;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	144 for kk=1:nn
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	145 o(kk,xx-1) = o(kk,xx-1) + sum(y(idxr+1:idxr+filtsz(kk),2));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	146 idxr = idxr+filtsz(kk);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	147 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	148 y(:,1) = y(:,2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	149 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	150
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	151 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	152
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	153 clear rns
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	154
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	155 % build output ao
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	156 for dd=1:nn
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	157 out(dd,1) = ao(tsdata(o(dd,:),fs));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	158 out(dd,1).setYunits(unit(yunit));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	159 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	160
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	161 a = matrix(out);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	162
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	163 case 's' % continuous input filters
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	164
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	165 o = zeros(nn,Ntot);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	166
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	167 T = 1/fs; % sampling period
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	168
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	169 for zz=1:nn % moving along system dimension
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	170
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	171 % extract residues and poles from input objects
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	172 % run over matrix dimensions
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	173 res = [];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	174 pls = [];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	175 filtsz = [];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	176 for jj=1:nn % collect filters coefficients along the columns zz
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	177 bfil = mfil(jj,zz);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	178 filtsz = [filtsz; numel(bfil.res)];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	179 res = [res; bfil.res];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	180 pls = [pls; bfil.poles];
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	181 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	182
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	183 % rescaling residues to get the correct result for univariate psd
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	184 res = res.*sqrt(fs/2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	185
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	186 Nrs = numel(res);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	187
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	188 % get covariance matrix for innovation
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	189 Ri = zeros(Nrs);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	190 for aa=1:Nrs
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	191 for bb=1:Nrs
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	192 Ri(aa,bb) = (res(aa)conj(res(bb)))(exp((pls(aa) + conj(pls(bb)))*T)-1)/(pls(aa) + conj(pls(bb)));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	193 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	194 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	195
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	196 % avoiding problems caused by roundoff errors
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	197 HH = triu(Ri,0); % get upper triangular part of R
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	198 HH1 = triu(Ri,1); % get upper triangular part of R above principal diagonal
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	199 HH2 = HH1'; % do transpose conjugate
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	200 Ri = HH + HH2; % reconstruct R in order to be really hermitian
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	201
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	202 % get singular value decomposition of R
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	203 [Ui,Si,Vi] = svd(Ri,0);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	204
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	205 % conversion matrix for innovation
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	206 Ai = Vi*sqrt(Si);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	207
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	208 % get covariance matrix for initial state
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	209 Rx = zeros(Nrs);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	210 for aa=1:Nrs
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	211 for bb=1:Nrs
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	212 Rx(aa,bb) = -1(res(aa)conj(res(bb)))/(pls(aa) + conj(pls(bb)));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	213 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	214 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	215
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	216 % avoiding problems caused by roundoff errors
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	217 HH = triu(Rx,0); % get upper triangular part of R
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	218 HH1 = triu(Rx,1); % get upper triangular part of R above principal diagonal
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	219 HH2 = HH1'; % do transpose conjugate
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	220 Rx = HH + HH2; % reconstruct R in order to be really hermitian
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	221
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	222 % get singular value decomposition of R
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	223 [Ux,Sx,Vx] = svd(Rx,0);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	224
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	225 % conversion matrix for initial state
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	226 Ax = Vx*sqrt(Sx);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	227
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	228 % generate unitary variance gaussian random noise
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	229 %ns = randn(Nrs,Ntot);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	230 ns = randn(Nrs,1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	231
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	232 % get correlated starting data points
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	233 cns = Ax*ns;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	234
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	235 % need to correct for roundoff errors
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	236 % cleaning up results for numerical approximations
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	237 idx = imag(pls(:,1))==0;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	238 cns(idx) = real(cns(idx));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	239
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	240 % cleaning up results for numerical roundoff errors
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	241 % states associated to complex conjugate poles must be complex
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	242 % conjugate
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	243 idxi = imag(pls(:,1))~=0;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	244 icns = cns(idxi);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	245 for jj = 1:2:numel(icns)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	246 icns(jj+1,1) = conj(icns(jj,1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	247 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	248 cns(idxi) = icns;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	249
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	250 y = zeros(sum(filtsz),2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	251 rnoise = zeros(sum(filtsz),1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	252 rns = randn(1,Ntot);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	253 %rns = utils.math.blwhitenoise(Ntot,fs,1/Nsecs,fs/2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	254 %rns = rns.'; % willing to work with raw
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	255
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	256 y(:,1) = cns;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	257
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	258 % start recurrence
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	259 for xx = 2:Ntot+1
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	260 % innov = Ai*randn(sum(filtsz),1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	261 rnoise(:,1) = rns(xx-1);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	262 innov = Ai*rnoise;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	263 % need to correct for roundoff errors
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	264 % cleaning up results for numerical approximations
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	265 innov(idx) = real(innov(idx));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	266
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	267 % cleaning up results for numerical roundoff errors
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	268 % states associated to complex conjugate poles must be complex
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	269 % conjugate
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	270 iinnov = innov(idxi);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	271 for jj = 1:2:numel(iinnov)
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	272 iinnov(jj+1,1) = conj(iinnov(jj,1));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	273 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	274 innov(idxi) = iinnov;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	275
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	276 y(:,2) = diag(exp(pls.T))y(:,1) + innov;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	277
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	278 idxr = 0;
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	279 for kk=1:nn
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	280 o(kk,xx-1) = o(kk,xx-1) + sum(y(idxr+1:idxr+filtsz(kk),2));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	281 idxr = idxr+filtsz(kk);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	282 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	283 y(:,1) = y(:,2);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	284 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	285
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	286 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	287
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	288 % clear rns
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	289
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	290 % build output ao
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	291 for dd=1:nn
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	292 out(dd,1) = ao(tsdata(o(dd,:),fs));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	293 out(dd,1).setYunits(unit(yunit));
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	294 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	295
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	296 a = matrix(out);
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	297
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	298 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	299
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	300
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	301 end
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	302
f0afece42f48 Import. Daniele Nicolodi <nicolodi@science.unitn.it> parents: diff changeset	303

0

f0afece42f48 Import.