Mercurial > hg > ltpda

comparison m-toolbox/test/test_ao_whiten2D_2.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/whiten2D
2 %
3 % L. Ferraioli 18-11-08
4 %
5 % $Id: test_ao_whiten2D_2.m,v 1.3 2009/08/26 18:27:06 luigi Exp $
6 %
7
8 %%
9
10 clear all
11
12 %% CSD Noise models and TF
13
14 [CSD,TF] = get_2D_test_obj_tf_psd();
15
16 % get noise coherence
17 coh = CSD(1,2)./sqrt(abs(CSD(1,1).*CSD(2,2)));
18 rcoh = real(coh);
19 rcoh.setName;
20 icoh = imag(coh);
21 icoh.setName;
22
23 %% Colored Noise generation
24
25 % Make white noise
26
27 a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 5e5));
28 a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 5e5));
29
30 % Generation from white noise
31
32 pl = plist(...
33 'csd11', CSD(1,1),...
34 'csd12', CSD(1,2),...
35 'csd21', CSD(2,1),...
36 'csd22', CSD(2,2),...
37 'MaxIter', 70, ...
38 'PoleType', 2, ...
39 'MinOrder', 25, ...
40 'MaxOrder', 35, ...
41 'Weights', 2, ...
42 'Plot', false,...
43 'MSEVARTOL', 1e-2,...
44 'FITTOL', 1e-4,...
45 'UseSym', 0,...
46 'Disp', false);
47
48 [ac1,ac2] = noisegen2D(a1,a2, pl);
49
50 ac1.setYunits('m');
51 ac2.setYunits('m');
52
53 %% Checking generated noise
54
55 acxx1 = ac1.psd(plist('Nfft',1e6));
56 acxx2 = ac2.psd(plist('Nfft',1e6));
57
58 acxy = cohere(ac1,ac2,plist('Nfft',5e4));
59 racxy = real(acxy);
60 racxy.setName;
61 iacxy = imag(acxy);
62 iacxy.setName;
63
64 iplot(abs(acxx1),abs(CSD(1,1)))
65 iplot(abs(acxx2),abs(CSD(2,2)))
66
67 iplot(racxy,rcoh,plist('Yscales',{'All','lin'}))
68 iplot(iacxy,icoh,plist('Yscales',{'All','lin'}))
69
70 %%
71 iplot(abs(acxx1),abs(acxx2))
72
73 %% Get Coloring filters
74
75 cf11 = ac1.procinfo.find('Filt11');
76 cf12 = ac1.procinfo.find('Filt12');
77 cf21 = ac2.procinfo.find('FILT21');
78 cf22 = ac2.procinfo.find('FILT22');
79
80 %% Get Theoretical CSD from coloring filters
81
82 rc11 = resp(cf11,plist('f',CSD(1,1).x,'bank','parallel'));
83 rc12 = resp(cf12,plist('f',CSD(1,1).x,'bank','parallel'));
84 rc21 = resp(cf21,plist('f',CSD(1,1).x,'bank','parallel'));
85 rc22 = resp(cf22,plist('f',CSD(1,1).x,'bank','parallel'));
86
87
88 psd11 = rc11.*conj(rc11)+rc12.*conj(rc12);
89 psd22 = rc21.*conj(rc21)+rc22.*conj(rc22);
90 psd12 = rc11.*conj(rc21)+rc12.*conj(rc22);
91 psd21 = conj(psd12);
92 coher = psd12./sqrt(abs(psd11).*abs(psd22));
93 rcoher = real(coher);
94 rcoher.setName;
95 icoher = imag(coher);
96 icoher.setName;
97
98 % compare with data
99 iplot(racxy,rcoher,plist('Yscales',{'All','lin'}))
100 iplot(iacxy,icoher,plist('Yscales',{'All','lin'}))
101
102 %% Noise whitening
103
104 pl = plist(...
105 'csd11', CSD(1,1),...
106 'csd12', CSD(1,2),...
107 'csd21', CSD(2,1),...
108 'csd22', CSD(2,2),...
109 'MaxIter', 75, ...
110 'PoleType', 3, ...
111 'MinOrder', 25, ...
112 'MaxOrder', 40, ...
113 'Weights', 2, ...
114 'Plot', false,...
115 'MSEVARTOL', 1e-2,...
116 'FITTOL', 1e-3,...
117 'UseSym', 0,...
118 'Disp', false);
119
120 [aw1,aw2] = whiten2D(ac1,ac2, pl);
121
122 aw1.setYunits('m');
123 aw2.setYunits('m');
124
125 %% Checking whitened noise
126
127 awxx1 = aw1.psd(plist('Nfft',1e6));
128 awxx2 = aw2.psd(plist('Nfft',1e6));
129 %%
130 iplot(awxx1,awxx2)
131
132 %% Get coherence
133
134 awxy = cohere(aw1,aw2,plist('Nfft',6e4));
135 rawxy = real(awxy);
136 rawxy.setName;
137 iawxy = imag(awxy);
138 iawxy.setName;
139
140 %% Check whitened data coherence
141 racxy.setName('a_xy');
142 iacxy.setName('a_xy');
143 rawxy.setName('wa_xy');
144 iawxy.setName('wa_xy');
145 iplot(racxy.split(plist('split_type','frequencies','frequencies',[4e-4 5])),rcoh.split(plist('split_type','frequencies','frequencies',[1e-4 5])),rawxy.split(plist('split_type','frequencies','frequencies',[4e-4 5])),plist('Yscales',{'All','lin'},'Ylabels',{'All','Re(coh)'}))
146 iplot(iacxy.split(plist('split_type','frequencies','frequencies',[4e-4 5])),icoh.split(plist('split_type','frequencies','frequencies',[1e-4 5])),iawxy.split(plist('split_type','frequencies','frequencies',[4e-4 5])),plist('Yscales',{'All','lin'},'Ylabels',{'All','Im(coh)'}))
147
148 %% Get whitening filters
149
150 wf11 = aw1.procinfo.find('Filt11');
151 wf12 = aw1.procinfo.find('Filt12');
152 wf21 = aw2.procinfo.find('FILT21');
153 wf22 = aw2.procinfo.find('FILT22');
154
155 % filters respose
156 rw11 = resp(wf11,plist('f',CSD(1,1).x,'bank','parallel'));
157 rw12 = resp(wf12,plist('f',CSD(1,1).x,'bank','parallel'));
158 rw21 = resp(wf21,plist('f',CSD(1,1).x,'bank','parallel'));
159 rw22 = resp(wf22,plist('f',CSD(1,1).x,'bank','parallel'));
160
161 %% Build expected csd and coherence of whitened data
162
163 wp11 = conj(rw11).*psd11+psd12.*conj(rw12);
164 wp12 = psd11.*conj(rw21)+psd12.*conj(rw22);
165 wp21 = psd21.*conj(rw11)+psd22.*conj(rw12);
166 wp22 = psd21.*conj(rw21)+psd22.*conj(rw22);
167
168 u11 = rw11.*wp11+rw12.*wp21;
169 u11.setName;
170 u12 = rw11.*wp12+rw12.*wp22;
171 u12.setName;
172 u21 = rw21.*wp11+rw22.*wp21;
173 u21.setName;
174 u22 = rw21.*wp12+rw22.*wp22;
175 u22.setName;
176 wcoher = u12./sqrt(abs(u11).*abs(u22));
177 rwcoher = real(wcoher);
178 rwcoher.setName;
179 iwcoher = imag(wcoher);
180 iwcoher.setName;
181
182 % compare with data
183 iplot(rawxy,rwcoher,plist('Yscales',{'All','lin'}))
184 iplot(iawxy,iwcoher,plist('Yscales',{'All','lin'}))
185
186 %% Coherence of starting data
187
188 axy = cohere(a1,a2,plist('Nfft',1e4));
189 raxy = real(axy);
190 raxy.setName;
191 iaxy = imag(axy);
192 iaxy.setName;
193
194 % compare with whitened
195 iplot(rawxy,raxy,plist('Yscales',{'All','lin'}))
196 iplot(iawxy,iaxy,plist('Yscales',{'All','lin'}))
197
198 %% Comparing Theoretical data
199
200 iplot(rcoher,rwcoher,plist('Yscales',{'All','lin'}))
201 iplot(icoher,iwcoher,plist('Yscales',{'All','lin'}))
202
203 iplot(rcoh,rwcoher,plist('Yscales',{'All','lin'}))
204 iplot(icoh,iwcoher,plist('Yscales',{'All','lin'}))
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224 % END