ltpda: m-toolbox/test/MDC2/test_mdc2_noisegen_accuracy

author	Daniele Nicolodi <nicolodi@science.unitn.it>
date	Wed, 23 Nov 2011 19:22:13 +0100
parents
children

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

0

f0afece42f48 Import.