comparison m-toolbox/classes/@matrix/MultiChannelNoise.m @ 0:f0afece42f48

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