--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_ao_confint.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,281 @@
+% Test script for ao/confint
+% 
+% L Ferraioli 20-03-09
+% 
+% $Id: test_ao_confint.m,v 1.4 2011/04/29 14:03:08 luigi Exp $
+
+%% params
+
+fs = 1; % Hz
+
+%% Set the model
+
+d = [1 -1.1 0.5];
+h11 = miir([1 -0.5],d,fs);
+h12 = miir([0 -0.5],d,fs);
+h21 = miir([0 0.4],d,fs);
+h22 = miir([1 -0.6],d,fs);
+
+%% Theoretical sectra
+
+f = logspace(-5,log10(0.5),300);
+f = f.';
+
+plresp = plist('f',f);
+
+rh11 = resp(h11,plresp);
+rh12 = resp(h12,plresp);
+rh21 = resp(h21,plresp);
+rh22 = resp(h22,plresp);
+
+% psd11
+G11 = rh11.y.*conj(rh11.y) + rh12.y.*conj(rh12.y);
+G11 = ao(plist('xvals', f, 'yvals', G11, 'fs', fs, 'type', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+G11.setName;
+
+% psd 22
+G22 = rh21.y.*conj(rh21.y) + rh22.y.*conj(rh22.y);
+G22 = ao(plist('xvals', f, 'yvals', G22, 'fs', fs, 'type', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+G22.setName;
+
+% cpsd
+G12 = conj(rh11.y).*rh21.y + conj(rh12.y).*rh22.y;
+G12 = ao(plist('xvals', f, 'yvals', G12, 'fs', fs, 'type', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+G12.setName;
+
+% mscohere
+sK12 = (G12.y.*conj(G12.y))./(G11.y.*G22.y);
+sK12 = ao(plist('xvals', f, 'yvals', sK12, 'fs', fs, 'type', 'fsdata'));
+sK12.setName;
+
+%% Noise generation
+
+w1 = ao(plist('tsfcn','randn(size(t))','fs',1,'nsecs',1e5,'yunits','m'));
+w2 = ao(plist('tsfcn','randn(size(t))','fs',1,'nsecs',1e5,'yunits','m'));
+
+x11 = filter(w1,h11);
+x12 = filter(w2,h12);
+x21 = filter(w1,h21);
+x22 = filter(w2,h22);
+
+x1 = x11 + x12;
+x2 = x21 + x22;
+
+%% mslcohere
+
+sk12log = lcohere(x1,x2,plist('type','MS'));
+if numel(sk12log)>1
+  sk12log = sk12log(1,2);
+end
+
+% iplot(sK12,sk12log,plist('Yscales',{'All','linear'}))
+
+out = confint(sk12log,plist('method','mslcohere'));
+lc = out.getObjectAtIndex(1);
+uc = out.getObjectAtIndex(2);
+var = out.getObjectAtIndex(3);
+
+iplot(sK12,sk12log,lc,uc,plist('Yscales',{'All','linear'}))
+
+dof = getdof(sk12log,plist('method','mslcohere'));
+
+%% shaded plot
+
+x = lc.data.x;
+y1 = lc.data.y;
+y2 = uc.data.y;
+mod = sK12;
+
+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(mod.data.x, mod.data.y);
+hf = plot(sk12log.data.x,sk12log.data.y,'Color', 'g');
+
+set(gca,'xscale','log','yscale','lin');
+% set(gca,'Layer','top')
+ylabel('Magnitude Squared Coherence');
+xlabel('Frequency [Hz]');
+legend([hd hf ha(2)],{'Model MSC','Sample MSC','95% Conf. level'});
+
+%% lpsd
+
+g11lpsd = lpsd(x1);
+
+% iplot(2.*G11,g11lpsd)
+
+outlpsd = confint(g11lpsd,plist('method','lpsd'));
+lclpsd = outlpsd.getObjectAtIndex(1);
+uclpsd = outlpsd.getObjectAtIndex(2);
+varlpsd = outlpsd.getObjectAtIndex(3);
+
+iplot(2.*G11,g11lpsd,lclpsd,uclpsd)
+
+doflpsd = getdof(g11lpsd,plist('method','lpsd'));
+
+%% shaded plot
+
+x = lclpsd.data.x;
+y1 = lclpsd.data.y;
+y2 = uclpsd.data.y;
+mod = 2.*G11;
+
+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(mod.data.x, mod.data.y);
+hf = plot(g11lpsd.data.x,g11lpsd.data.y,'Color', 'g');
+
+set(gca,'xscale','log','yscale','log');
+% set(gca,'Layer','top')
+ylabel('Power Spectral Density [m^{2} / Hz]');
+xlabel('Frequency [Hz]');
+legend([hd hf ha(2)],{'Model Spectrum','Sample Spectrum','95% Conf. level'});
+
+%% mscohere
+
+sk12lin = cohere(x1,x2,plist('Nfft',1e3,'type','MS'));
+if numel(sk12lin)>1
+  sk12lin = sk12lin(1,2);
+end
+
+iplot(sK12,sk12lin,plist('Yscales',{'All','linear'}))
+
+outlin = confint(sk12lin,plist('method','mscohere'));
+lclin = outlin.getObjectAtIndex(1);
+uclin = outlin.getObjectAtIndex(2);
+varlin = outlin.getObjectAtIndex(3);
+
+iplot(sK12,sk12lin,lclin,uclin,plist('Yscales',{'All','linear'}))
+
+doflin = getdof(sk12lin,plist('method','mscohere'));
+
+%% shaded plot
+
+x = lclin.data.x;
+y1 = lclin.data.y;
+y2 = uclin.data.y;
+mod = sK12;
+
+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(mod.data.x, mod.data.y);
+hf = plot(sk12lin.data.x,sk12lin.data.y,'Color', 'g');
+
+set(gca,'xscale','log','yscale','lin');
+% set(gca,'Layer','top')
+ylabel('Magnitude Squared Coherence');
+xlabel('Frequency [Hz]');
+legend([hd hf ha(2)],{'Model MSC','Sample MSC','95% Conf. level'});
+
+%% psd
+
+g11psd = psd(x1,plist('Nfft',1e4));
+
+iplot(2.*G11,g11psd)
+
+outpsd = confint(g11psd,plist('method','psd'));
+lcpsd = outpsd.getObjectAtIndex(1);
+ucpsd = outpsd.getObjectAtIndex(2);
+varpsd = outpsd.getObjectAtIndex(3);
+
+iplot(2.*G11,g11psd,lcpsd,ucpsd)
+
+dofpsd = getdof(g11psd,plist('method','psd'));
+
+%% shaded plot
+
+x = lcpsd.data.x;
+y1 = lcpsd.data.y;
+y2 = ucpsd.data.y;
+mod = 2.*G11;
+
+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(mod.data.x, mod.data.y);
+hf = plot(g11psd.data.x,g11psd.data.y,'Color', 'g');
+
+set(gca,'xscale','log','yscale','log');
+% set(gca,'Layer','top')
+ylabel('Power Spectral Density [m^{2} / Hz]');
+xlabel('Frequency [Hz]');
+legend([hd hf ha(2)],{'Model Spectrum','Sample Spectrum','95% Conf. level'});
+
+%% Shaded plot 2
+
+x = [lcpsd.data.x; flipud(lcpsd.data.x)];
+y = [ucpsd.data.y; flipud(lcpsd.data.y)];
+mod = 2.*G11;
+
+figure
+ha = fill(x,y,'r');
+set(ha,'EdgeColor','r');
+
+grid on
+hold on
+
+hd = plot(mod.data.x, mod.data.y);
+hf = plot(g11psd.data.x,g11psd.data.y,'Color','g');
+
+set(gca,'xscale','log','yscale','log');
+% set(gca,'Layer','top')
+ylabel('Power Spectral Density [m^{2} / Hz]');
+xlabel('Frequency [Hz]');
+legend([hd hf ha],{'Model Spectrum','Sample Spectrum','95% Conf. level'});
+
+%% Test it is working also with processed data
+
+g11psd = psd(x1,plist('Nfft',1e4));
+
+g11psd.setName('psd');
+
+plsp = plist('samples',[5 numel(g11psd.x)]);
+g11psds = split(g11psd,plsp);
+
+outpsd = confint(g11psds,plist('method','psd'));
+lcpsd = outpsd.getObjectAtIndex(1);
+ucpsd = outpsd.getObjectAtIndex(2);
+varpsd = outpsd.getObjectAtIndex(3);
+
+iplot(2.*G11,g11psds,lcpsd,ucpsd)
+
+dofpsd = getdof(g11psd,plist('method','psd'));
+