Mercurial > hg > ltpda

comparison m-toolbox/test/test_ao_noisegen2D_mdc2.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 mdc2 models and test procedure accuracy
4 %
5 % L. Ferraioli 04-02-2009
6 %
7 % $Id: test_ao_noisegen2D_mdc2.m,v 1.3 2009/02/19 17:46:41 luigi Exp $
8 %
9
10 %% General use variables and vectors
11
12 f = logspace(-6,log10(5),300);
13 fs = 10;
14 Nsecs = 1e5; % number of seconds
15 Nfft = 1e5; % number of samples for the fft
16 pls = plist('Nfft', Nfft,'Order',0); % plist for spectra
17
18 %% MDC2 Models
19
20 b = ao(plist('built-in','mdc2r2_fd_ltpnoise','f1',1e-6,'f2',5,'nf',300));
21 CSD = [b(1) b(2);conj(b(2)) b(3)];
22
23 %% Make white noise
24
25 a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs));
26 a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs));
27 a1.setYunits('m');
28 a2.setYunits('m');
29 a = [a1 a2];
30
31 % axx = a.cpsd(pls);
32
33 %%
34 iplot(a)
35
36 %% some ploting
37
38 iplot(axx)
39
40 iplot(CSD(1,1),CSD(1,2),CSD(2,1),CSD(2,2))
41
42 %% Noise generation
43
44 pl = plist(...
45 'csd11', CSD(1,1), ...
46 'csd12', CSD(1,2), ...
47 'csd21', CSD(2,1), ...
48 'csd22', CSD(2,2), ...
49 'MaxIter', 80, ...
50 'PoleType', 2, ...
51 'MinOrder', 35, ...
52 'MaxOrder', 40, ...
53 'Weights', 2, ...
54 'Plot', false,...
55 'FitTolerance', 2,...
56 'RMSEVar', 7,...
57 'UseSym', 0,...
58 'Disp', false);
59
60 ac = noisegen2D(a, pl);
61
62 %% Checking results and starting data
63
64 % iplot(a)
65 iplot(ac)
66
67 %% Making cross-spectrum
68
69 acxx = ac.cpsd(pls);
70 acch = ac.cohere(pls);
71
72 %% Plotting spectra
73
74 % iplot(acxx);
75
76 iplot(abs(acxx(1,1)),abs(CSD(1,1))) % model data need to be multiplied by 2 because acxx is the onesided cpsd
77 iplot(abs(acxx(1,2)),abs(CSD(1,2)))
78 iplot(abs(acxx(2,2)),abs(CSD(2,2)))
79
80 iplot(abs(acxx(1,1)),abs(acxx(1,2)),abs(acxx(2,2)))
81
82 %% Plotting coherence
83
84 iplot(acch(2,1),(abs(CSD(1,2)).^2)./(CSD(1,1).*CSD(2,2)))
85
86 %%
87
88 m1=mean(acxx(2,2).data.y(end-10,end)) % calculate average on the tail of channel 2
89 m2=mean(CSD(2,2).data.y(end-5,end))% calculate average on the tail of channel 2
90 m1/m2 % verify that the ratio is near 1
91
92 %%
93 % ************************************************************************
94 % Some more analysis for testing the accuracy of noise generation procedure
95 % ************************************************************************
96
97 %% Extracting filters from data
98
99 Filt11 = find(ac(1).procinfo,'Filt11');
100 Filt12 = find(ac(1).procinfo,'Filt12');
101 Filt21 = find(ac(2).procinfo,'Filt21');
102 Filt22 = find(ac(2).procinfo,'Filt22');
103
104 %% Calculating filters responses
105
106 tr11 = resp(Filt11,plist('f',f));
107 rFilt11 = tr11(1);
108 for ii = 2:numel(tr11)
109 rFilt11 = rFilt11 + tr11(ii);
110 end
111 rFilt11.setName('rFilt11', 'internal');
112
113 tr12 = resp(Filt12,plist('f',f));
114 rFilt12 = tr12(1);
115 for ii = 2:numel(tr12)
116 rFilt12 = rFilt12 + tr12(ii);
117 end
118 rFilt12.setName('rFilt12', 'internal');
119
120 tr21 = resp(Filt21,plist('f',f));
121 rFilt21 = tr21(1);
122 for ii = 2:numel(tr21)
123 rFilt21 = rFilt21 + tr21(ii);
124 end
125 rFilt21.setName('rFilt21', 'internal');
126
127 tr22 = resp(Filt22,plist('f',f));
128 rFilt22 = tr22(1);
129 for ii = 2:numel(tr22)
130 rFilt22 = rFilt22 + tr22(ii);
131 end
132 rFilt22.setName('rFilt22', 'internal');
133
134 %% Obtaining transfer functions
135
136 % calculating transfer functions from data
137 etf = tfe(a,ac,pls);
138
139 %% Comparing Filters Responses with estimated TFs (e-TFs)
140
141 % Comparing filters responses and calculated TFs
142 pl = plist('Legends', {'Filter Response','e-TF'});
143 iplot(rFilt11,etf(1,3),pl)
144 iplot(rFilt12,etf(2,3),pl)
145 iplot(rFilt21,etf(1,4),pl)
146 iplot(rFilt22,etf(2,4),pl)
147
148 %% Filtering data separately
149
150 % This operation is performed internally to the noisegen2D. Output data are
151 % then obtained by b1 = b11 + b12 and b2 = b21 + b22
152 b11 = filter(a1,plist('filter',Filt11,'bank','parallel'));
153 b12 = filter(a2,plist('filter',Filt12,'bank','parallel'));
154 b21 = filter(a1,plist('filter',Filt21,'bank','parallel'));
155 b22 = filter(a2,plist('filter',Filt22,'bank','parallel'));
156
157 %% Extracting transfer functions from separately filtered data se-TFs
158
159 etf11 = tfe(a1,b11,pls);
160 etf12 = tfe(a2,b12,pls);
161 etf21 = tfe(a1,b21,pls);
162 etf22 = tfe(a2,b22,pls);
163
164 %% Comparing separately-estimated TFs (se-TFs) with filter responses
165
166 pl = plist('Legends', {'Filter Response','se-TF'});
167 iplot(rFilt11,etf11(1,2),pl)
168 iplot(rFilt12,etf12(1,2),pl)
169 iplot(rFilt21,etf21(1,2),pl)
170 iplot(rFilt22,etf22(1,2),pl)
171
172 %% Comparing filters with TFs obtained by eigendecomposition
173
174 % This function output transfer functions as they are obtained by the
175 % eigendecomposition process. i.e. before the fitting process
176
177 icsd11 = CSD(1,1).data.y*fs/2;
178 icsd12 = CSD(1,2).data.y*fs/2;
179 icsd21 = CSD(2,1).data.y*fs/2;
180 icsd22 = CSD(2,2).data.y*fs/2;
181
182 [tf11,tf12,tf21,tf22] = utils.math.eigcsd(icsd11,icsd12,icsd21,icsd22,'USESYM',0,'DIG',50,'OTP','TF');
183
184 % Making AOs
185 eigtf11 = ao(fsdata(f,tf11,fs));
186 eigtf12 = ao(fsdata(f,tf12,fs));
187 eigtf21 = ao(fsdata(f,tf21,fs));
188 eigtf22 = ao(fsdata(f,tf22,fs));
189
190 %% Comparing eig-TFs with output filters
191
192 % Compare TFs before and after the fitting process
193
194 pl = plist('Legends', {'eig-TF','Filter Response'});
195 iplot(eigtf11,rFilt11,pl)
196 iplot(eigtf12,rFilt12,pl)
197 iplot(eigtf21,rFilt21,pl)
198 iplot(eigtf22,rFilt22,pl)
199
200 %% Phase difference between eig-TFs and output filters
201
202 % checking that phase differences between TFs are preserved by the fitting
203 % process
204 eigph1 = unwrap(angle(eigtf11)-angle(eigtf21));
205 eigph1.setYunits('rad')
206 filtph1 = unwrap(angle(rFilt11)-angle(rFilt21));
207 filtph1.setYunits('rad')
208 eigph2 = unwrap(angle(eigtf22)-angle(eigtf12));
209 eigph2.setYunits('rad')
210 filtph2 = unwrap(angle(rFilt22)-angle(rFilt12));
211 filtph2.setYunits('rad')
212
213 pl = plist('Legends',{'eig-TF','Filter'},'YScales',{'All','lin'});
214 iplot(eigph1+2*pi,filtph1,pl)
215 iplot(eigph2,filtph2,pl)
216
217
218 %% Comparing eig-TFs with se-TFs
219
220 % Compare eigendecomposition results with separately estimated TFs (se-TFs)
221 pl = plist('Legends', {'eig-TF','se-TF'});
222 iplot(eigtf11,etf11(1,2),pl)
223 iplot(eigtf12,etf12(1,2),pl)
224 iplot(eigtf21,etf21(1,2),pl)
225 iplot(eigtf22,etf22(1,2),pl)
226
227 %% Phase difference between eig-TFs and se-TFs
228
229 % checking that phase differences between TFs are preserved by the fitting
230 % process also for the filtered data (se-TFs)
231 eigph1 = unwrap(angle(eigtf11)-angle(eigtf21));
232 filtph1 = unwrap(angle(etf11(1,2))-angle(etf21(1,2)));
233 eigph2 = unwrap(angle(eigtf22)-angle(eigtf12));
234 filtph2 = unwrap(angle(etf22(1,2))-angle(etf12(1,2)));
235
236 pl = plist('Legends',{'eig-TF \Delta\phi','se-TF \Delta\phi'},'YScales',{'All','lin'});
237 iplot(eigph1,filtph1,pl)
238 iplot(eigph2,filtph2,pl)
239
240 % END