%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TEST matrix/fromCSD
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HISTORY:     23-04-2009 L Ferraioli
%                 Creation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VERSION '$Id: test_matrix_fromCSD.m,v 1.3 2009/11/06 17:00:30 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 - double precision All pass

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

na = matrix(plns);

%% 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 - double 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',    'off');

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 - double precision All pass

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

na = matrix(plns);

%% 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',  20, ...
  'MaxOrder',  30, ...
  'Weights',   2, ...
  'Plot',      true,...
  'MSEVARTOL', 1e-2,...
  'FITTOL',    1e-3,...
  'UseSym',    'on');

na = matrix(plns);

%% Check response

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

est_csd11 = na11.*conj(na11) + na12.*conj(na12);
est_csd11.setName;
est_csd12 = na11.*conj(na21) + na12.*conj(na22);
est_csd12.setName;
est_csd22 = na12.*conj(na12) + na22.*conj(na22);
est_csd22.setName;

iplot(csd11,est_csd11)
iplot(csd12,est_csd12)
iplot(csd22,est_csd22)


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

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

na = matrix(plns);

%% 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',  35, ...
  'Weights',   2, ...
  'Plot',      false,...
  'MSEVARTOL', 1e-1,...
  'FITTOL',    1e-3,...
  '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,(2/fs).*ECSD.objs(1,1))
iplot(csd12,(2/fs).*ECSD.objs(1,2))
iplot(csd22,(2/fs).*ECSD.objs(2,2))
























































