--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/MDC2/test_mdc2_noisegen_accuracy_v2.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,323 @@
+% A test script for to check the accuracy of the noise generation procedure
+% for the MDC2
+% 
+% 
+% L. Ferraioli 04-02-2009
+% 
+% $Id: test_mdc2_noisegen_accuracy_v2.m,v 1.3 2009/04/20 14:34:12 luigi Exp $
+% 
+
+%% General use variables and vectors
+
+fs = 1;
+f = logspace(-6,log10(fs/2),300);
+f = f.';
+Nsecs = 1e6; % number of seconds
+Nfft1 = 1e5;
+plsp1 = plist('Nfft',Nfft1,'order',2);
+Nfft2 = 5e3;
+plsp2 = plist('Nfft',Nfft2);
+
+%% MDC2 Models
+
+b = ao(plist('built-in','mdc2r2_fd_ltpnoise','f',f));
+
+%%
+
+iplot(b)
+
+%%
+
+CSD = [b(1) b(2);conj(b(2)) b(3)];
+MSC = b(2).*conj(b(2))./(b(1).*b(3));
+
+%% some ploting
+
+iplot(CSD(1,1),CSD(2,2),MSC)
+iplot(CSD(1,1),CSD(2,2),CSD(1,2),CSD(2,1))
+
+%% Make white noise
+
+a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
+a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
+
+%% Noise generation
+
+plng = plist(...
+    'csd11', CSD(1,1), ...
+    'csd12', CSD(1,2), ...
+    'csd21', CSD(2,1), ...
+    'csd22', CSD(2,2), ...
+    'MaxIter', 80, ...
+    'PoleType', 2, ...
+    'MinOrder', 35, ...
+    'MaxOrder', 40, ...
+    'Weights', 2, ...
+    'Plot', false,...
+    'MSEVARTOL', 1e-2,...
+    'FITTOL', 1e-5,...
+    'UseSym', 0,...
+    'Disp', false);
+
+[ac1,ac2] = noisegen2D(a1,a2, plng);
+
+% extracting filters
+h11 = ac1.procinfo.find('Filt11');
+h12 = ac1.procinfo.find('Filt12');
+h21 = ac2.procinfo.find('Filt21');
+h22 = ac2.procinfo.find('Filt22');
+
+%% Find stationary filter state
+
+Nrun = 20;
+
+for jj = 1:Nrun % filter state evolution
+  disp(num2str(jj))
+  ac11 = filter(a1,h11);
+  h11 = ac11.procinfo.find('Filters'); % prooagate histout
+  ac12 = filter(a2,h12);
+  h12 = ac12.procinfo.find('Filters');
+  ac21 = filter(a1,h21);
+  h21 = ac21.procinfo.find('Filters');
+  ac22 = filter(a2,h22);
+  h22 = ac22.procinfo.find('Filters');
+end
+
+% re-define ac1 and ac2
+ac1 = ac11 + ac12;
+ac1.setName;
+ac2 = ac21 + ac22;
+ac2.setName;
+
+%% spectra calculation
+
+acxx1 = psd(ac1,plsp1);
+acxx2 = psd(ac2,plsp1);
+acxx12 = cohere(ac1,ac2,plsp2);
+
+% %%
+% iplot(CSD11,acxx1,CSD22,acxx2)
+% iplot(MSC,acxx12)
+% 
+% %%
+
+% get dofs
+dof1 = getdof(acxx1,plist('method','psd'));
+dof2 = getdof(acxx2,plist('method','psd'));
+dof12 = getdof(acxx12,plist('method','mscohere'));
+
+% Interpolate theorethical data
+CSD11 = interp(CSD(1,1),plist('vertices',acxx1.data.x));
+CSD22 = interp(CSD(2,2),plist('vertices',acxx2.data.x));
+MSC = interp(MSC,plist('vertices',acxx12.data.x));
+
+% getting variance
+[lw1,up1,sig1] = confint(acxx1,plist('method','psd'));
+[lw2,up2,sig2] = confint(acxx2,plist('method','psd'));
+[lw12,up12,sig12] = confint(acxx12,plist('method','mscohere'));
+
+% set limits of reliable regions
+[a,b] = size(acxx12.x);
+la = round(a*0.0212);
+
+% chi square
+k1 = acxx1 - CSD11;
+k1 = k1.^2;
+k1 = k1./sig1;
+k1 = split(k1,plist('split_type','samples','samples',[la a]));
+chi1 = sum(k1)./length(k1.data.x);
+
+k2 = acxx2 - CSD22;
+k2 = k2.^2;
+k2 = k2./sig2;
+k2 = split(k2,plist('split_type','samples','samples',[la a]));
+chi2 = sum(k2)./length(k2.data.x);
+
+k12 = acxx12 - MSC;
+k12= k12.^2;
+k12 = k12./sig12;
+k12 = split(k12,plist('split_type','samples','samples',[la a]));
+chi12 = sum(k12)./length(k12.data.x);
+
+proc1 = acxx1.procinfo;
+proc2 = acxx2.procinfo;
+proc12 = acxx12.procinfo;
+
+%% 
+
+iplot(CSD11,acxx1,lw1,up1)
+iplot(CSD22,acxx2,lw2,up2)
+iplot(MSC,acxx12,lw12,up12)
+
+%% Running the loop
+
+N = 100;
+
+for ii = 2:N
+  disp(num2str(ii))
+  
+  a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
+  a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
+
+  x11 = filter(a1,h11);
+  x12 = filter(a2,h12);
+  x21 = filter(a1,h21);
+  x22 = filter(a2,h22);
+
+  x1 = x11 + x12;
+  x2 = x21 + x22;
+  
+  % sectra
+  xx1 = psd(x1,plsp1);
+  xx2 = psd(x2,plsp1); 
+  xx12 = cohere(x1,x2,plsp2);
+%   xx12 = xx12(1,2);
+  
+  % sum spectra
+  acxx1 = acxx1 + xx1;
+  acxx2 = acxx2 + xx2;
+  acxx12 = acxx12 + xx12;
+  
+  % getting variance
+  [lw,up,sig1] = confint(xx1,plist('method','psd'));
+  [lw,up,sig2] = confint(xx2,plist('method','psd'));
+  [lw,up,sig12] = confint(xx12,plist('method','mscohere'));
+
+  % chi square
+  k1 = xx1 - CSD11;
+  k1 = k1.^2;
+  k1 = k1./sig1;
+  k1 = split(k1,plist('split_type','samples','samples',[la a]));
+  xc1 = sum(k1)./length(k1.data.x);
+
+  k2 = xx2 - CSD22;
+  k2 = k2.^2;
+  k2 = k2./sig2;
+  k2 = split(k2,plist('split_type','samples','samples',[la a]));
+  xc2 = sum(k2)./length(k2.data.x);
+
+  k12 = xx12 - MSC;
+  k12= k12.^2;
+  k12 = k12./sig12;
+  k12 = split(k12,plist('split_type','samples','samples',[la a]));
+  xc12 = sum(k12)./length(k12.data.x);
+  
+  % adding chisquare
+  chi1 = chi1 + xc1;
+  chi2 = chi2 + xc2;
+  chi12 = chi12 + xc12;
+  
+end
+
+%% Do Averages
+
+% do average
+acxx1 = acxx1./N;
+acxx1.setName;
+acxx2 = acxx2./N;
+acxx2.setName;
+acxx12 = acxx12./N;
+acxx12.setName;
+
+% average chi square
+chi1 = chi1./N;
+chi2 = chi2./N;
+chi12 = chi12./N;
+
+% dofs for the averages
+adof1 = dof1.*N;
+adof2 = dof2.*N;
+adof12 = dof12.*N;
+
+% get confidence bounds for the 'true' spactra
+[lw1,up1,sig1] = confint(acxx1,plist('method','psd','dof',adof1));
+[lw2,up2,sig2] = confint(acxx2,plist('method','psd','dof',adof2));
+[lw12,up12,sig12] = confint(acxx12,plist('method','mscohere','dof',adof12));
+
+%% plotting
+
+iplot(CSD(1,1),acxx1)
+iplot(CSD(2,2),acxx2)
+iplot(MSC,acxx12)
+
+iplot(CSD(1,1),lw1,up1)
+iplot(CSD(2,2),lw2,up2)
+iplot(MSC,lw12,up12)
+
+%% plotting conflevels 1
+
+x = lw1.data.x;
+y1 = lw1.data.y;
+y2 = up1.data.y;
+
+figure
+y = [y1 (y2-y1)]; % y1 and y2 are columns
+ha = area(x, y);
+set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
+set(ha(2), 'FaceColor', 'r')
+set(ha, 'LineStyle', 'none')
+grid on
+
+% plot the line edges
+hold on 
+hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
+hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
+hd = plot(CSD(1,1).data.x, CSD(1,1).data.y);
+
+set(gca,'xscale','log','yscale','log');
+% set(gca,'Layer','top')
+ylabel('Power Spectrum [m^{2} / Hz]');
+xlabel('Frequency [Hz]');
+legend([hd hb],{'Model CSD(1,1)','95% Conf. level'});
+
+%% plotting conflevels 2
+
+x = lw2.data.x;
+y1 = lw2.data.y;
+y2 = up2.data.y;
+
+figure
+y = [y1 (y2-y1)]; % y1 and y2 are columns
+ha = area(x, y);
+set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
+set(ha(2), 'FaceColor', 'r')
+set(ha, 'LineStyle', 'none')
+grid on
+
+% plot the line edges
+hold on 
+hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
+hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
+hd = plot(CSD(2,2).data.x, CSD(2,2).data.y);
+
+set(gca,'xscale','log','yscale','log');
+% set(gca,'Layer','top')
+ylabel('Power Spectrum [m^{2} / Hz]');
+xlabel('Frequency [Hz]');
+legend([hd hb],{'Model CSD(2,2)','95% Conf. level'});
+
+%% plotting conflevels 3
+
+x = lw12.data.x;
+y1 = lw12.data.y;
+y2 = up12.data.y;
+
+figure
+y = [y1 (y2-y1)]; % y1 and y2 are columns
+ha = area(x, y);
+set(ha(1), 'FaceColor', 'none') % this makes the bottom area invisible
+set(ha(2), 'FaceColor', 'r')
+set(ha, 'LineStyle', 'none')
+grid on
+
+% plot the line edges
+hold on 
+hb = plot(x, y1, 'LineWidth', 1, 'Color', 'r');
+hc = plot(x, y2, 'LineWidth', 1, 'Color', 'r');
+hd = plot(MSC.data.x, MSC.data.y);
+
+set(gca,'yscale','log');
+% set(gca,'Layer','top')
+ylabel('Magnitude Square Coherence');
+xlabel('Frequency [Hz]');
+legend([hd hb],{'Model MSC','95% Conf. level'});