Mercurial > hg > ltpda
view m-toolbox/test/MDC2/test_mdc2_noisegen_accuracy_v2.m @ 17:7afc99ec5f04 database-connection-manager
Update ao_model_retrieve_in_timespan
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
% 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'});