0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % A test script for to check the accuracy of the noise generation procedure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 % for the MDC2
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 % L. Ferraioli 04-02-2009
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 % $Id: test_mdc2_noisegen_accuracy_v2.m,v 1.3 2009/04/20 14:34:12 luigi Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 %% General use variables and vectors
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 fs = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 f = logspace(-6,log10(fs/2),300);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 f = f.';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 Nsecs = 1e6; % number of seconds
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 Nfft1 = 1e5;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 plsp1 = plist('Nfft',Nfft1,'order',2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 Nfft2 = 5e3;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 plsp2 = plist('Nfft',Nfft2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 %% MDC2 Models
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 b = ao(plist('built-in','mdc2r2_fd_ltpnoise','f',f));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 %%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 iplot(b)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 %%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 CSD = [b(1) b(2);conj(b(2)) b(3)];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 MSC = b(2).*conj(b(2))./(b(1).*b(3));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 %% some ploting
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 iplot(CSD(1,1),CSD(2,2),MSC)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 iplot(CSD(1,1),CSD(2,2),CSD(1,2),CSD(2,1))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 %% Make white noise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 %% Noise generation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46 plng = plist(...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 'csd11', CSD(1,1), ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 'csd12', CSD(1,2), ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 'csd21', CSD(2,1), ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 'csd22', CSD(2,2), ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 'MaxIter', 80, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 'PoleType', 2, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 'MinOrder', 35, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 'MaxOrder', 40, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 'Weights', 2, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 'Plot', false,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 'MSEVARTOL', 1e-2,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 'FITTOL', 1e-5,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 'UseSym', 0,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 'Disp', false);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 [ac1,ac2] = noisegen2D(a1,a2, plng);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 % extracting filters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 h11 = ac1.procinfo.find('Filt11');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 h12 = ac1.procinfo.find('Filt12');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 h21 = ac2.procinfo.find('Filt21');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 h22 = ac2.procinfo.find('Filt22');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 %% Find stationary filter state
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 Nrun = 20;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 for jj = 1:Nrun % filter state evolution
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 disp(num2str(jj))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 ac11 = filter(a1,h11);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 h11 = ac11.procinfo.find('Filters'); % prooagate histout
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 ac12 = filter(a2,h12);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79 h12 = ac12.procinfo.find('Filters');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 ac21 = filter(a1,h21);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 h21 = ac21.procinfo.find('Filters');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 ac22 = filter(a2,h22);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 h22 = ac22.procinfo.find('Filters');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 % re-define ac1 and ac2
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87 ac1 = ac11 + ac12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 ac1.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 ac2 = ac21 + ac22;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 ac2.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92 %% spectra calculation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94 acxx1 = psd(ac1,plsp1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 acxx2 = psd(ac2,plsp1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 acxx12 = cohere(ac1,ac2,plsp2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 % %%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 % iplot(CSD11,acxx1,CSD22,acxx2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100 % iplot(MSC,acxx12)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102 % %%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 % get dofs
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105 dof1 = getdof(acxx1,plist('method','psd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106 dof2 = getdof(acxx2,plist('method','psd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107 dof12 = getdof(acxx12,plist('method','mscohere'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109 % Interpolate theorethical data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 CSD11 = interp(CSD(1,1),plist('vertices',acxx1.data.x));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 CSD22 = interp(CSD(2,2),plist('vertices',acxx2.data.x));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112 MSC = interp(MSC,plist('vertices',acxx12.data.x));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114 % getting variance
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115 [lw1,up1,sig1] = confint(acxx1,plist('method','psd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116 [lw2,up2,sig2] = confint(acxx2,plist('method','psd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
117 [lw12,up12,sig12] = confint(acxx12,plist('method','mscohere'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
118
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
119 % set limits of reliable regions
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
120 [a,b] = size(acxx12.x);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
121 la = round(a*0.0212);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
122
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
123 % chi square
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
124 k1 = acxx1 - CSD11;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
125 k1 = k1.^2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
126 k1 = k1./sig1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
127 k1 = split(k1,plist('split_type','samples','samples',[la a]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
128 chi1 = sum(k1)./length(k1.data.x);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
129
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
130 k2 = acxx2 - CSD22;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
131 k2 = k2.^2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
132 k2 = k2./sig2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
133 k2 = split(k2,plist('split_type','samples','samples',[la a]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
134 chi2 = sum(k2)./length(k2.data.x);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
135
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
136 k12 = acxx12 - MSC;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
137 k12= k12.^2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
138 k12 = k12./sig12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
139 k12 = split(k12,plist('split_type','samples','samples',[la a]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
140 chi12 = sum(k12)./length(k12.data.x);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
141
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
142 proc1 = acxx1.procinfo;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
143 proc2 = acxx2.procinfo;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
144 proc12 = acxx12.procinfo;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
145
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
146 %%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
147
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
148 iplot(CSD11,acxx1,lw1,up1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
149 iplot(CSD22,acxx2,lw2,up2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
150 iplot(MSC,acxx12,lw12,up12)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
151
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
152 %% Running the loop
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
153
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
154 N = 100;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
155
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
156 for ii = 2:N
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
157 disp(num2str(ii))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
158
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
159 a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
160 a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
161
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
162 x11 = filter(a1,h11);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
163 x12 = filter(a2,h12);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
164 x21 = filter(a1,h21);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
165 x22 = filter(a2,h22);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
166
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
167 x1 = x11 + x12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
168 x2 = x21 + x22;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
169
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
170 % sectra
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
171 xx1 = psd(x1,plsp1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
172 xx2 = psd(x2,plsp1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
173 xx12 = cohere(x1,x2,plsp2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
174 % xx12 = xx12(1,2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
175
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
176 % sum spectra
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
177 acxx1 = acxx1 + xx1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
178 acxx2 = acxx2 + xx2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
179 acxx12 = acxx12 + xx12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
180
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
181 % getting variance
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
182 [lw,up,sig1] = confint(xx1,plist('method','psd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
183 [lw,up,sig2] = confint(xx2,plist('method','psd'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
184 [lw,up,sig12] = confint(xx12,plist('method','mscohere'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
185
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
186 % chi square
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
187 k1 = xx1 - CSD11;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
188 k1 = k1.^2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
189 k1 = k1./sig1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
190 k1 = split(k1,plist('split_type','samples','samples',[la a]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
191 xc1 = sum(k1)./length(k1.data.x);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
192
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
193 k2 = xx2 - CSD22;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
194 k2 = k2.^2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
195 k2 = k2./sig2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
196 k2 = split(k2,plist('split_type','samples','samples',[la a]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
197 xc2 = sum(k2)./length(k2.data.x);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
198
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
199 k12 = xx12 - MSC;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
200 k12= k12.^2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
201 k12 = k12./sig12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
202 k12 = split(k12,plist('split_type','samples','samples',[la a]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
203 xc12 = sum(k12)./length(k12.data.x);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
204
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
205 % adding chisquare
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
206 chi1 = chi1 + xc1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
207 chi2 = chi2 + xc2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
208 chi12 = chi12 + xc12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
209
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
210 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
211
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
212 %% Do Averages
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
213
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
214 % do average
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
215 acxx1 = acxx1./N;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
216 acxx1.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
217 acxx2 = acxx2./N;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
218 acxx2.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
219 acxx12 = acxx12./N;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
220 acxx12.setName;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
221
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
222 % average chi square
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
223 chi1 = chi1./N;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
224 chi2 = chi2./N;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
225 chi12 = chi12./N;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
226
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
227 % dofs for the averages
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
228 adof1 = dof1.*N;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
229 adof2 = dof2.*N;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
230 adof12 = dof12.*N;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
231
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
232 % get confidence bounds for the 'true' spactra
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
233 [lw1,up1,sig1] = confint(acxx1,plist('method','psd','dof',adof1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
234 [lw2,up2,sig2] = confint(acxx2,plist('method','psd','dof',adof2));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
235 [lw12,up12,sig12] = confint(acxx12,plist('method','mscohere','dof',adof12));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
236
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
237 %% plotting
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
238
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
239 iplot(CSD(1,1),acxx1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
240 iplot(CSD(2,2),acxx2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
241 iplot(MSC,acxx12)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
242
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
243 iplot(CSD(1,1),lw1,up1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
244 iplot(CSD(2,2),lw2,up2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
245 iplot(MSC,lw12,up12)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
246
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
247 %% plotting conflevels 1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
248
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
249 x = lw1.data.x;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
250 y1 = lw1.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
251 y2 = up1.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
252
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
253 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
254 y = [y1 (y2-y1)]; % y1 and y2 are columns
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
255 ha = area(x, y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
256 set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
257 set(ha(2), 'FaceColor', 'r')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
258 set(ha, 'LineStyle', 'none')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
259 grid on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
260
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
261 % plot the line edges
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
262 hold on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
263 hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
264 hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
265 hd = plot(CSD(1,1).data.x, CSD(1,1).data.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
266
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
267 set(gca,'xscale','log','yscale','log');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
268 % set(gca,'Layer','top')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
269 ylabel('Power Spectrum [m^{2} / Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
270 xlabel('Frequency [Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
271 legend([hd hb],{'Model CSD(1,1)','95% Conf. level'});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
272
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
273 %% plotting conflevels 2
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
274
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
275 x = lw2.data.x;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
276 y1 = lw2.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
277 y2 = up2.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
278
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
279 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
280 y = [y1 (y2-y1)]; % y1 and y2 are columns
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
281 ha = area(x, y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
282 set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
283 set(ha(2), 'FaceColor', 'r')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
284 set(ha, 'LineStyle', 'none')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
285 grid on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
286
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
287 % plot the line edges
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
288 hold on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
289 hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
290 hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
291 hd = plot(CSD(2,2).data.x, CSD(2,2).data.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
292
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
293 set(gca,'xscale','log','yscale','log');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
294 % set(gca,'Layer','top')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
295 ylabel('Power Spectrum [m^{2} / Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
296 xlabel('Frequency [Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
297 legend([hd hb],{'Model CSD(2,2)','95% Conf. level'});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
298
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
299 %% plotting conflevels 3
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
300
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
301 x = lw12.data.x;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
302 y1 = lw12.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
303 y2 = up12.data.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
304
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
305 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
306 y = [y1 (y2-y1)]; % y1 and y2 are columns
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
307 ha = area(x, y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
308 set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
309 set(ha(2), 'FaceColor', 'r')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
310 set(ha, 'LineStyle', 'none')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
311 grid on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
312
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
313 % plot the line edges
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
314 hold on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
315 hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
316 hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
317 hd = plot(MSC.data.x, MSC.data.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
318
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
319 set(gca,'yscale','log');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
320 % set(gca,'Layer','top')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
321 ylabel('Magnitude Square Coherence');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
322 xlabel('Frequency [Hz]');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
323 legend([hd hb],{'Model MSC','95% Conf. level'});
|