Mercurial > hg > ltpda
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 |