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