% A test script for ao/noisegen2D
% 
% DESCRIPTION: Run noisegen2D with test data and test procedure accuracy
% 
% L. Ferraioli 10-11-08
% 
% $Id: test_ao_noisegen2D.m,v 1.9 2010/05/03 19:04:45 luigi Exp $
% 

%% General use variables and vectors

userdir = 'C:\Users\Luigi'; % You should set your own dir

f = logspace(-6,log10(5),300);
fs = 10;

Nt = fs*10000;

%% Make white noise

a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e5));
a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e5));
a = [a1 a2];

%% CSD Noise models

fundir = fullfile(userdir,'ltp_data_analysis\MDCs\MDC1_UTN');
cf = cd;
cd(fundir);
[TF,CSD] = mdc1_tf_models(plist('f',f,'fs',fs));
cd(cf);

%% Input data psd

a1xx = a1.psd;
a2xx = a2.psd;

%% Making csd of input data

pl  = plist('Nfft', Nt);
iCSD = cpsd(a,pl);

%% plot input data psd

iplot(a1xx,a2xx)

%% Plot Model Tfs

iplot(TF(1,1),TF(1,2),TF(2,1),TF(2,2))

%% Plot Model CSD

iplot(CSD(1,1),CSD(1,2),CSD(2,1),CSD(2,2))

%% Noise generation

pl = plist(...
    'csd11', CSD(1,1), ...
    'csd12', CSD(1,2), ...
    'csd21', CSD(2,1), ...
    'csd22', CSD(2,2), ...
    'MaxIter', 80, ...
    'PoleType', 2, ...
    'MinOrder', 15, ...
    'MaxOrder', 45, ...
    'Weights', 3, ...
    'FITTOL', 1e-3,...
    'MSEVARTOL', 1e-1,...
    'UseSym', 0,...
    'Plot', false,...
    'Disp', false);

ac = noisegen2D(a, pl);

%% Checking results and starting data

% iplot(a)
iplot(ac)

%% Making cross-spectrum

% acxx = ac.cpsd;
plpsd = plist('navs',2,'order',1,'olap',50);
acxx1 = ac(1).psd(plpsd);
acxx2 = ac(2).psd(plpsd);

iplot(acxx1,CSD(1,1),acxx2,CSD(2,2))

%% Plotting spectra

% iplot(acxx);

iplot(acxx(1,1),CSD(1,1))
iplot(abs(acxx(1,2)),abs(CSD(1,2)))
iplot(acxx(2,2),CSD(2,2))

%%

m1=mean(acxx(2,2).data.y(end-10,end)) % calculate average on the tail of channel 2
m2=mean(CSD(2,2).data.y(end-5,end))% calculate average on the tail of channel 2
m1/m2 % verify that the ratio is near 1

%% 
% ************************************************************************
% Some more analysis for testing the accuracy of noise generation procedure
% ************************************************************************

%% Extracting filters from data

Filt11 = find(ac(1).procinfo,'FILT11');
Filt12 = find(ac(1).procinfo,'FILT12');
Filt21 = find(ac(2).procinfo,'FILT21');
Filt22 = find(ac(2).procinfo,'FILT22');

%% Calculating filters responses

tr11 = resp(Filt11,plist('f',f));
rFilt11 = tr11(1);
for ii = 2:numel(tr11)
  rFilt11 = rFilt11 + tr11(ii);
end
rFilt11.setName('rFilt11', 'internal');

tr12 = resp(Filt12,plist('f',f));
rFilt12 = tr12(1);
for ii = 2:numel(tr12)
  rFilt12 = rFilt12 + tr12(ii);
end
rFilt12.setName('rFilt12', 'internal');

tr21 = resp(Filt21,plist('f',f));
rFilt21 = tr21(1);
for ii = 2:numel(tr21)
  rFilt21 = rFilt21 + tr21(ii);
end
rFilt21.setName('rFilt21', 'internal');

tr22 = resp(Filt22,plist('f',f));
rFilt22 = tr22(1);
for ii = 2:numel(tr22)
  rFilt22 = rFilt22 + tr22(ii);
end
rFilt22.setName('rFilt22', 'internal');

%% Obtaining transfer functions

% calculating transfer functions from data
pl  = plist('Nfft', Nt);
etf11 = tfe(a(1),ac(1),pl);
etf12 = tfe(a(2),ac(1),pl);
etf21 = tfe(a(1),ac(2),pl);
etf22 = tfe(a(2),ac(2),pl);

%% Comparing Filters Responses with estimated TFs (e-TFs)

% Comparing filters responses and calculated TFs
pl = plist('Legends', {'Filter Response','e-TF'});
iplot(rFilt11,etf11,pl)
iplot(rFilt12,etf12,pl)
iplot(rFilt21,etf21,pl)
iplot(rFilt22,etf22,pl)

%% Comparing starting TFs (s-TFs) with estimated TFs

pl = plist('Legends', {'s-TF','e-TF'});
iplot(TF(1,1),etf11/sqrt(fs))
iplot(TF(1,2),etf12/sqrt(fs))
iplot(TF(2,1),etf21/sqrt(fs))
iplot(TF(2,2),etf22/sqrt(fs))

%% Building CSD from estimated TFs (e-TFs)

% Output CSD is obtained as
% eCSD = [TF11 TF12;TF21 TF22]*iCSD*[TF11 TF12;TF21 TF22]'
eCSD = [etf(1,3) etf(2,3);etf(1,4) etf(2,4)]*iCSD*[etf(1,3) etf(2,3);etf(1,4) etf(2,4)]';

ecsd11 = eCSD(1,1);
ecsd12 = eCSD(1,2);
ecsd22 = eCSD(2,2);

% ecsd11 = etf(1,3).*conj(etf(1,3))+etf(2,3).*conj(etf(2,3));
% ecsd12 = etf(1,3).*conj(etf(1,4))+etf(2,3).*conj(etf(2,4));
% ecsd22 = etf(2,4).*conj(etf(2,4))+etf(1,4).*conj(etf(1,4));

%% Comparing original CSD with e-TFs CSD

pl = plist('Legends', {'Original CSD','e-TF CSD'});
iplot(abs(CSD(1,1)),ecsd11/(fs),pl)
iplot(abs(CSD(1,2)),ecsd12/(fs),pl)
iplot(abs(CSD(2,2)),ecsd22/(fs),pl)

%% Filtering data separately

% This operation is performed internally to the noisegen2D. Output data are
% then obtained by b1 = b11 + b12 and b2 = b21 + b22
b11 = filter(a1,plist('filter',Filt11,'bank','parallel'));
b12 = filter(a2,plist('filter',Filt12,'bank','parallel'));
b21 = filter(a1,plist('filter',Filt21,'bank','parallel'));
b22 = filter(a2,plist('filter',Filt22,'bank','parallel'));

%% Extracting transfer functions from separately filtered data

pl  = plist('Nfft', Nt);
etf11 = tfe(a1,b11,pl);
etf12 = tfe(a2,b12,pl);
etf21 = tfe(a1,b21,pl);
etf22 = tfe(a2,b22,pl);

%% Comparing separately-estimated TFs (se-TFs) with filter responses

pl = plist('Legends', {'Filter Response','se-TF'});
iplot(rFilt11,etf11,pl)
iplot(rFilt12,etf12,pl)
iplot(rFilt21,etf21,pl)
iplot(rFilt22,etf22,pl)

%% Building CSD from se-TFs

% Output CSD is obtained as
% eCSD = [TF11 TF12;TF21 TF22]*iCSD*[TF11 TF12;TF21 TF22]'
seCSD = [etf11 etf12;etf21 etf22]*iCSD*[etf11 etf12;etf21 etf22]';

secsd11 = seCSD(1,1);
secsd12 = seCSD(1,2);
secsd22 = seCSD(2,2);

% secsd11 = etf11(1,2).*conj(etf11(1,2))+etf12(1,2).*conj(etf12(1,2));
% secsd12 = etf11(1,2).*conj(etf21(1,2))+etf12(1,2).*conj(etf22(1,2));
% secsd22 = etf22(1,2).*conj(etf22(1,2))+etf21(1,2).*conj(etf21(1,2));

%% Comparing original CSD with e-TFs CSD

pl = plist('Legends', {'Original CSD','se-TF CSD'});
iplot(CSD(1,1),secsd11/(fs),pl)
iplot(CSD(1,2),secsd12/(fs),pl)
iplot(CSD(2,2),secsd22/(fs),pl)

%% Comparing filters with TFs obtained by eigendecomposition

% This function output transfer functions as they are obtained by the
% eigendecomposition process. i.e. before the fitting process

icsd11 = CSD(1,1).data.y*fs/2;
icsd12 = CSD(1,2).data.y*fs/2;
icsd21 = CSD(2,1).data.y*fs/2;
icsd22 = CSD(2,2).data.y*fs/2;

[tf11,tf12,tf21,tf22] = utils.math.eigcsd(icsd11,icsd12,icsd21,icsd22,'USESYM',0,'DIG',50,'OTP','TF');

% Making AOs
eigtf11 = ao(fsdata(f,tf11,fs));
eigtf12 = ao(fsdata(f,tf12,fs));
eigtf21 = ao(fsdata(f,tf21,fs));
eigtf22 = ao(fsdata(f,tf22,fs));

%% Comparing eig-TFs with output filters

% Compare TFs before and after the fitting process

pl = plist('Legends', {'eig-TF','Filter Response'});
iplot(eigtf11,rFilt11,pl)
iplot(eigtf12,rFilt12,pl)
iplot(eigtf21,rFilt21,pl)
iplot(eigtf22,rFilt22,pl)

%% Phase difference between eig-TFs and output filters

% checking that phase differences between TFs are preserved by the fitting
% process
eigph1 = unwrap(angle(eigtf11)-angle(eigtf21));
filtph1 = unwrap(angle(rFilt11)-angle(rFilt21));
eigph2 = unwrap(angle(eigtf22)-angle(eigtf12));
filtph2 = unwrap(angle(rFilt22)-angle(rFilt12));

pl = plist('Legends',{'eig-TF \Delta\phi','Filter \Delta\phi'},'YScales',{'All','lin'});
iplot(eigph1,filtph1+2*pi,pl)
iplot(eigph2,filtph2+2*pi,pl)

%% Comparing eig-TFs with se-TFs

% Compare eigendecomposition results with separately estimated TFs (se-TFs)
pl = plist('Legends', {'eig-TF','se-TF'});
iplot(eigtf11,etf11(1,2),pl)
iplot(eigtf12,etf12(1,2),pl)
iplot(eigtf21,etf21(1,2),pl)
iplot(eigtf22,etf22(1,2),pl)

%% Phase difference between eig-TFs and se-TFs

% checking that phase differences between TFs are preserved by the fitting
% process also for the filtered data (se-TFs)
eigph1 = unwrap(angle(eigtf11)-angle(eigtf21));
filtph1 = unwrap(angle(etf11(1,2))-angle(etf21(1,2)));
eigph2 = unwrap(angle(eigtf22)-angle(eigtf12));
filtph2 = unwrap(angle(etf22(1,2))-angle(etf12(1,2)));

pl = plist('Legends',{'eig-TF \Delta\phi','se-TF \Delta\phi'},'YScales',{'All','lin'});
iplot(eigph1,filtph1+2*pi,pl)
iplot(eigph2,filtph2,pl)

% END