%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TEST matrix/MultiChannelNoise
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HISTORY:     23-04-2009 L Ferraioli
%                 Creation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VERSION

'$Id: test_matrix_fromCSD.m,v 1.3 2009/04/23 16:03:55 luigi Exp $';

%% Loading spectra

load ..\m-toolbox\test\mpsd.mat % load mpsd.mat first column is f then psd1, csd and psd2

f = mpsd(:,1);
psd = mpsd(:,2);
fs = 10;

% 1dim model
mod1D = ao(plist('xvals', f, 'yvals', psd, 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH1'));
mod1D.setName;

% 2dim model
csd11 = ao(plist('xvals', f, 'yvals', mpsd(:,2), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH1 PSD'));
csd11.setName;
csd12 = ao(plist('xvals', f, 'yvals', mpsd(:,3), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH12 CSD'));
csd12.setName;
csd21 = ao(plist('xvals', f, 'yvals', conj(mpsd(:,3)), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH21 CSD'));
csd21.setName;
csd22 = ao(plist('xvals', f, 'yvals', mpsd(:,4), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH2 PSD'));
csd22.setName;

mod2D = [csd11 csd12;csd21 csd22];

%%
%--------------------------------------------------------------------------
% 1 Dimensional Noisegen Filter
%--------------------------------------------------------------------------

%% 1 dim noisegen filter in s domain - symbolic precision All pass

% plist for noise generation
plns = plist(...
  'csd',       mod1D, ...
  'targetobj',  'parfrac', ...
  'Nsecs',     1e4, ...
  'fs',        fs, ...
  'MaxIter',   70, ...
  'PoleType',  3, ...
  'MinOrder',  15, ...
  'MaxOrder',  30, ...
  'Weights',   2, ...
  'Plot',      true,...
  'MSEVARTOL', 1e-2,...
  'FITTOL',    1e-3,...
  'UseSym',    'on');

na = matrix(plns);



%% 1 dim noisegen filter in z domain - symbolic precision All pass

% plist for noise generation
plns = plist(...
  'csd',       mod1D, ...
  'targetobj',  'miir', ...
  'Nsecs',     1e4, ...
  'fs',        fs, ...
  'MaxIter',   70, ...
  'PoleType',  3, ...
  'MinOrder',  10, ...
  'MaxOrder',  30, ...
  'Weights',   2, ...
  'Plot',      true,...
  'MSEVARTOL', 1e-2,...
  'FITTOL',    1e-3,...
  'UseSym',    'on');

na = matrix(plns);

%%
%--------------------------------------------------------------------------
% 2 Dimensional Noisegen Filter
%--------------------------------------------------------------------------


%% 2 dim noisegen filter in s domain - symbolic precision All pass

% plist for noise generation
plns = plist(...
  'csd',       mod2D, ...
  'targetobj',  'parfrac', ...
  'MaxIter',   70, ...
  'PoleType',  3, ...
  'MinOrder',  25, ...
  'MaxOrder',  40, ...
  'Weights',   2, ...
  'Plot',      false,...
  'MSEVARTOL', 1e-2,...
  'FITTOL',    1e-4,...
  'UseSym',    'on');

na = matrix(plns);

%% Check response

na11 = resp(na.objs(1,1),plist('f',f,'bank','parallel'));
na12 = resp(na.objs(1,2),plist('f',f,'bank','parallel'));
na21 = resp(na.objs(2,1),plist('f',f,'bank','parallel'));
na22 = resp(na.objs(2,2),plist('f',f,'bank','parallel'));

mna = matrix([na11 na12;na21 na22]);

ECSD = mna*conj(transpose(mna));

iplot(csd11,ECSD.objs(1,1))
iplot(csd12,ECSD.objs(1,2))
iplot(csd22,ECSD.objs(2,2))


%% 2 dim noisegen filter in z domain - symbolic precision All pass

% plist for noise generation
plns = plist(...
  'csd',       mod2D, ...
  'targetobj',  'miir', ...
  'fs',        fs, ...
  'MaxIter',   70, ...
  'PoleType',  3, ...
  'MinOrder',  25, ...
  'MaxOrder',  40, ...
  'Weights',   2, ...
  'Plot',      false,...
  'MSEVARTOL', 1e-1,...
  'FITTOL',    1e-4,...
  'UseSym',    'on');

na = matrix(plns);

%% Check response

na11 = resp(na.objs(1,1).filters,plist('f',f,'bank','parallel'));
na12 = resp(na.objs(1,2).filters,plist('f',f,'bank','parallel'));
na21 = resp(na.objs(2,1).filters,plist('f',f,'bank','parallel'));
na22 = resp(na.objs(2,2).filters,plist('f',f,'bank','parallel'));

mna = matrix([na11 na12;na21 na22]);

ECSD = mna*conj(transpose(mna));

iplot(csd11,ECSD.objs(1,1))
iplot(csd12,ECSD.objs(1,2))
iplot(csd22,ECSD.objs(2,2))

%%
%--------------------------------------------------------------------------
% 2 Dimensional Noise Generation
%--------------------------------------------------------------------------

plng = plist('model',na,...
  'Nsecs',2e5,...
  'fs',fs);

no = matrix(plng);

%% Check psd

no1 = no.objs(1);
no2 = no.objs(2);

no1xx = no1.psd;
no2xx = no2.psd;

iplot(no1xx,csd11,no2xx,csd22)

%% Check cross-coherence

coh = csd12./sqrt(csd11.*csd22);
rcoh = real(coh);
icoh = imag(coh);

no12xx = cpsd(no1,no2,plist('nfft',1e4));
no1xx = no1.psd(plist('nfft',1e4));
no2xx = no2.psd(plist('nfft',1e4));

ecoh = no12xx./(sqrt(no1xx.*no2xx));
recoh = real(ecoh);
iecoh = imag(ecoh);

iplot(recoh,rcoh,plist('Yscales',{'All','lin'}))
iplot(iecoh,icoh,plist('Yscales',{'All','lin'}))