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