line source
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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))
