Mercurial > hg > ltpda
view m-toolbox/test/test_ao_confint.m @ 18:947e2ff4b1b9 database-connection-manager
Update plist.FROM_REPOSITORY_PLIST and plist.TO_REPOSITORY_PLIST
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
% 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'));