Mercurial > hg > ltpda

comparison m-toolbox/test/test_ao_noisegen2D.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 % A test script for ao/noisegen2D
2 %
3 % DESCRIPTION: Run noisegen2D with test data and test procedure accuracy
4 %
5 % L. Ferraioli 10-11-08
6 %
7 % $Id: test_ao_noisegen2D.m,v 1.9 2010/05/03 19:04:45 luigi Exp $
8 %
9
10 %% General use variables and vectors
11
12 userdir = 'C:\Users\Luigi'; % You should set your own dir
13
14 f = logspace(-6,log10(5),300);
15 fs = 10;
16
17 Nt = fs*10000;
18
19 %% Make white noise
20
21 a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e5));
22 a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e5));
23 a = [a1 a2];
24
25 %% CSD Noise models
26
27 fundir = fullfile(userdir,'ltp_data_analysis\MDCs\MDC1_UTN');
28 cf = cd;
29 cd(fundir);
30 [TF,CSD] = mdc1_tf_models(plist('f',f,'fs',fs));
31 cd(cf);
32
33 %% Input data psd
34
35 a1xx = a1.psd;
36 a2xx = a2.psd;
37
38 %% Making csd of input data
39
40 pl = plist('Nfft', Nt);
41 iCSD = cpsd(a,pl);
42
43 %% plot input data psd
44
45 iplot(a1xx,a2xx)
46
47 %% Plot Model Tfs
48
49 iplot(TF(1,1),TF(1,2),TF(2,1),TF(2,2))
50
51 %% Plot Model CSD
52
53 iplot(CSD(1,1),CSD(1,2),CSD(2,1),CSD(2,2))
54
55 %% Noise generation
56
57 pl = plist(...
58 'csd11', CSD(1,1), ...
59 'csd12', CSD(1,2), ...
60 'csd21', CSD(2,1), ...
61 'csd22', CSD(2,2), ...
62 'MaxIter', 80, ...
63 'PoleType', 2, ...
64 'MinOrder', 15, ...
65 'MaxOrder', 45, ...
66 'Weights', 3, ...
67 'FITTOL', 1e-3,...
68 'MSEVARTOL', 1e-1,...
69 'UseSym', 0,...
70 'Plot', false,...
71 'Disp', false);
72
73 ac = noisegen2D(a, pl);
74
75 %% Checking results and starting data
76
77 % iplot(a)
78 iplot(ac)
79
80 %% Making cross-spectrum
81
82 % acxx = ac.cpsd;
83 plpsd = plist('navs',2,'order',1,'olap',50);
84 acxx1 = ac(1).psd(plpsd);
85 acxx2 = ac(2).psd(plpsd);
86
87 iplot(acxx1,CSD(1,1),acxx2,CSD(2,2))
88
89 %% Plotting spectra
90
91 % iplot(acxx);
92
93 iplot(acxx(1,1),CSD(1,1))
94 iplot(abs(acxx(1,2)),abs(CSD(1,2)))
95 iplot(acxx(2,2),CSD(2,2))
96
97 %%
98
99 m1=mean(acxx(2,2).data.y(end-10,end)) % calculate average on the tail of channel 2
100 m2=mean(CSD(2,2).data.y(end-5,end))% calculate average on the tail of channel 2
101 m1/m2 % verify that the ratio is near 1
102
103 %%
104 % ************************************************************************
105 % Some more analysis for testing the accuracy of noise generation procedure
106 % ************************************************************************
107
108 %% Extracting filters from data
109
110 Filt11 = find(ac(1).procinfo,'FILT11');
111 Filt12 = find(ac(1).procinfo,'FILT12');
112 Filt21 = find(ac(2).procinfo,'FILT21');
113 Filt22 = find(ac(2).procinfo,'FILT22');
114
115 %% Calculating filters responses
116
117 tr11 = resp(Filt11,plist('f',f));
118 rFilt11 = tr11(1);
119 for ii = 2:numel(tr11)
120 rFilt11 = rFilt11 + tr11(ii);
121 end
122 rFilt11.setName('rFilt11', 'internal');
123
124 tr12 = resp(Filt12,plist('f',f));
125 rFilt12 = tr12(1);
126 for ii = 2:numel(tr12)
127 rFilt12 = rFilt12 + tr12(ii);
128 end
129 rFilt12.setName('rFilt12', 'internal');
130
131 tr21 = resp(Filt21,plist('f',f));
132 rFilt21 = tr21(1);
133 for ii = 2:numel(tr21)
134 rFilt21 = rFilt21 + tr21(ii);
135 end
136 rFilt21.setName('rFilt21', 'internal');
137
138 tr22 = resp(Filt22,plist('f',f));
139 rFilt22 = tr22(1);
140 for ii = 2:numel(tr22)
141 rFilt22 = rFilt22 + tr22(ii);
142 end
143 rFilt22.setName('rFilt22', 'internal');
144
145 %% Obtaining transfer functions
146
147 % calculating transfer functions from data
148 pl = plist('Nfft', Nt);
149 etf11 = tfe(a(1),ac(1),pl);
150 etf12 = tfe(a(2),ac(1),pl);
151 etf21 = tfe(a(1),ac(2),pl);
152 etf22 = tfe(a(2),ac(2),pl);
153
154 %% Comparing Filters Responses with estimated TFs (e-TFs)
155
156 % Comparing filters responses and calculated TFs
157 pl = plist('Legends', {'Filter Response','e-TF'});
158 iplot(rFilt11,etf11,pl)
159 iplot(rFilt12,etf12,pl)
160 iplot(rFilt21,etf21,pl)
161 iplot(rFilt22,etf22,pl)
162
163 %% Comparing starting TFs (s-TFs) with estimated TFs
164
165 pl = plist('Legends', {'s-TF','e-TF'});
166 iplot(TF(1,1),etf11/sqrt(fs))
167 iplot(TF(1,2),etf12/sqrt(fs))
168 iplot(TF(2,1),etf21/sqrt(fs))
169 iplot(TF(2,2),etf22/sqrt(fs))
170
171 %% Building CSD from estimated TFs (e-TFs)
172
173 % Output CSD is obtained as
174 % eCSD = [TF11 TF12;TF21 TF22]*iCSD*[TF11 TF12;TF21 TF22]'
175 eCSD = [etf(1,3) etf(2,3);etf(1,4) etf(2,4)]*iCSD*[etf(1,3) etf(2,3);etf(1,4) etf(2,4)]';
176
177 ecsd11 = eCSD(1,1);
178 ecsd12 = eCSD(1,2);
179 ecsd22 = eCSD(2,2);
180
181 % ecsd11 = etf(1,3).*conj(etf(1,3))+etf(2,3).*conj(etf(2,3));
182 % ecsd12 = etf(1,3).*conj(etf(1,4))+etf(2,3).*conj(etf(2,4));
183 % ecsd22 = etf(2,4).*conj(etf(2,4))+etf(1,4).*conj(etf(1,4));
184
185 %% Comparing original CSD with e-TFs CSD
186
187 pl = plist('Legends', {'Original CSD','e-TF CSD'});
188 iplot(abs(CSD(1,1)),ecsd11/(fs),pl)
189 iplot(abs(CSD(1,2)),ecsd12/(fs),pl)
190 iplot(abs(CSD(2,2)),ecsd22/(fs),pl)
191
192 %% Filtering data separately
193
194 % This operation is performed internally to the noisegen2D. Output data are
195 % then obtained by b1 = b11 + b12 and b2 = b21 + b22
196 b11 = filter(a1,plist('filter',Filt11,'bank','parallel'));
197 b12 = filter(a2,plist('filter',Filt12,'bank','parallel'));
198 b21 = filter(a1,plist('filter',Filt21,'bank','parallel'));
199 b22 = filter(a2,plist('filter',Filt22,'bank','parallel'));
200
201 %% Extracting transfer functions from separately filtered data
202
203 pl = plist('Nfft', Nt);
204 etf11 = tfe(a1,b11,pl);
205 etf12 = tfe(a2,b12,pl);
206 etf21 = tfe(a1,b21,pl);
207 etf22 = tfe(a2,b22,pl);
208
209 %% Comparing separately-estimated TFs (se-TFs) with filter responses
210
211 pl = plist('Legends', {'Filter Response','se-TF'});
212 iplot(rFilt11,etf11,pl)
213 iplot(rFilt12,etf12,pl)
214 iplot(rFilt21,etf21,pl)
215 iplot(rFilt22,etf22,pl)
216
217 %% Building CSD from se-TFs
218
219 % Output CSD is obtained as
220 % eCSD = [TF11 TF12;TF21 TF22]*iCSD*[TF11 TF12;TF21 TF22]'
221 seCSD = [etf11 etf12;etf21 etf22]*iCSD*[etf11 etf12;etf21 etf22]';
222
223 secsd11 = seCSD(1,1);
224 secsd12 = seCSD(1,2);
225 secsd22 = seCSD(2,2);
226
227 % secsd11 = etf11(1,2).*conj(etf11(1,2))+etf12(1,2).*conj(etf12(1,2));
228 % secsd12 = etf11(1,2).*conj(etf21(1,2))+etf12(1,2).*conj(etf22(1,2));
229 % secsd22 = etf22(1,2).*conj(etf22(1,2))+etf21(1,2).*conj(etf21(1,2));
230
231 %% Comparing original CSD with e-TFs CSD
232
233 pl = plist('Legends', {'Original CSD','se-TF CSD'});
234 iplot(CSD(1,1),secsd11/(fs),pl)
235 iplot(CSD(1,2),secsd12/(fs),pl)
236 iplot(CSD(2,2),secsd22/(fs),pl)
237
238 %% Comparing filters with TFs obtained by eigendecomposition
239
240 % This function output transfer functions as they are obtained by the
241 % eigendecomposition process. i.e. before the fitting process
242
243 icsd11 = CSD(1,1).data.y*fs/2;
244 icsd12 = CSD(1,2).data.y*fs/2;
245 icsd21 = CSD(2,1).data.y*fs/2;
246 icsd22 = CSD(2,2).data.y*fs/2;
247
248 [tf11,tf12,tf21,tf22] = utils.math.eigcsd(icsd11,icsd12,icsd21,icsd22,'USESYM',0,'DIG',50,'OTP','TF');
249
250 % Making AOs
251 eigtf11 = ao(fsdata(f,tf11,fs));
252 eigtf12 = ao(fsdata(f,tf12,fs));
253 eigtf21 = ao(fsdata(f,tf21,fs));
254 eigtf22 = ao(fsdata(f,tf22,fs));
255
256 %% Comparing eig-TFs with output filters
257
258 % Compare TFs before and after the fitting process
259
260 pl = plist('Legends', {'eig-TF','Filter Response'});
261 iplot(eigtf11,rFilt11,pl)
262 iplot(eigtf12,rFilt12,pl)
263 iplot(eigtf21,rFilt21,pl)
264 iplot(eigtf22,rFilt22,pl)
265
266 %% Phase difference between eig-TFs and output filters
267
268 % checking that phase differences between TFs are preserved by the fitting
269 % process
270 eigph1 = unwrap(angle(eigtf11)-angle(eigtf21));
271 filtph1 = unwrap(angle(rFilt11)-angle(rFilt21));
272 eigph2 = unwrap(angle(eigtf22)-angle(eigtf12));
273 filtph2 = unwrap(angle(rFilt22)-angle(rFilt12));
274
275 pl = plist('Legends',{'eig-TF \Delta\phi','Filter \Delta\phi'},'YScales',{'All','lin'});
276 iplot(eigph1,filtph1+2*pi,pl)
277 iplot(eigph2,filtph2+2*pi,pl)
278
279 %% Comparing eig-TFs with se-TFs
280
281 % Compare eigendecomposition results with separately estimated TFs (se-TFs)
282 pl = plist('Legends', {'eig-TF','se-TF'});
283 iplot(eigtf11,etf11(1,2),pl)
284 iplot(eigtf12,etf12(1,2),pl)
285 iplot(eigtf21,etf21(1,2),pl)
286 iplot(eigtf22,etf22(1,2),pl)
287
288 %% Phase difference between eig-TFs and se-TFs
289
290 % checking that phase differences between TFs are preserved by the fitting
291 % process also for the filtered data (se-TFs)
292 eigph1 = unwrap(angle(eigtf11)-angle(eigtf21));
293 filtph1 = unwrap(angle(etf11(1,2))-angle(etf21(1,2)));
294 eigph2 = unwrap(angle(eigtf22)-angle(eigtf12));
295 filtph2 = unwrap(angle(etf22(1,2))-angle(etf12(1,2)));
296
297 pl = plist('Legends',{'eig-TF \Delta\phi','se-TF \Delta\phi'},'YScales',{'All','lin'});
298 iplot(eigph1,filtph1+2*pi,pl)
299 iplot(eigph2,filtph2,pl)
300
301 % END