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'});