Mercurial > hg > ltpda

view m-toolbox/test/test_ao_fromCSD.m @ 52:daf4eab1a51e database-connection-manager tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix. Default password should be [] not an empty string
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 07 Dec 2011 17:29:47 +0100
parents f0afece42f48
children
line wrap: on
line source

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% TEST ao/fromCSD
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HISTORY:     23-04-2009 L Ferraioli
%                 Creation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VERSION

'$Id: test_ao_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];


%% Generate 1D Noise

% plist for noise generation

plns = plist(...
  'csd',       mod1D, ...
  'Nsecs',     1e4, ...
  'fs',        fs, ...
  'xunits',    's', ...
  'yunits',    '', ...
  'MaxIter',   50, ...
  'PoleType',  3, ...
  'MinOrder',  7, ...
  'MaxOrder',  35, ...
  'Weights',   2, ...
  'Plot',      false,...
  'MSEVARTOL', 1e-2,...
  'FITTOL',    1e-2);

na = ao(plns);

%% Compare with Noisegen1D

aw = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e4));

pl = plist(...
    'model', mod1D, ...
    'MaxIter', 50, ...
    'PoleType', 3, ...
    'MinOrder', 7, ...
    'MaxOrder', 35, ...
    'Weights', 2, ...
    'Plot', false,...
    'Disp', false,...
    'MSEVARTOL', 1e-2,...
    'FITTOL', 1e-2);

nat = noisegen1D(aw, pl);

%% Make psd

naxx = na.lpsd;
natxx = nat.lpsd;


iplot(naxx,natxx,mod1D)

%% 2 dim noise generation

% plist for noise generation

plns = plist(...
  'csd',       mod2D, ...
  'Nsecs',     1e4, ...
  'fs',        fs, ...
  'xunits',    's', ...
  'yunits',    '', ...
  'MaxIter',   70, ...
  'PoleType',  3, ...
  'MinOrder',  20, ...
  'MaxOrder',  40, ...
  'Weights',   2, ...
  'Plot',      false,...
  'MSEVARTOL', 1e-2,...
  'FITTOL',    1e-2);

na = ao(plns);

%% Compare with Noisegen2D

aw1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e4));
aw2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e4));

pl = plist(...
    'csd11', csd11, ...
    'csd12', csd12, ...
    'csd21', csd21, ...
    'csd22', csd22, ...
    'MaxIter', 70, ...
    'PoleType', 3, ...
    'MinOrder', 20, ...
    'MaxOrder', 40, ...
    'Weights', 2, ...
    'Plot', false,...
    'Disp', false,...
    'MSEVARTOL', 1e-2,...
    'FITTOL', 1e-2);

[nat1,nat2] = noisegen2D(aw1,aw2, pl);

%% Make psd

naxx1 = na(1).lpsd;
natxx1 = nat1.lpsd;
naxx2 = na(2).lpsd;
natxx2 = nat2.lpsd;


iplot(naxx1,natxx1,csd11)
iplot(naxx2,natxx2,csd22)


%% 3 Dim noise generation

% sampling frequency
fs = 1; % Hz
% frequency vector
f = logspace(-6,log10(0.5),300);
f = f.';

%%% Set the 3 channel model
h11 = pzmodel(plist('GAIN', [0.01], 'POLES', [pz(9.9999999999999995e-007,NaN) pz(0.01,NaN) pz(0.02,NaN)], 'ZEROS', [pz(0.0001,NaN) pz(0.002,NaN)]));
h12 = pzmodel(plist('GAIN', [0.0001], 'POLES', [pz(0.10000000000000001,NaN) pz(0.0030000000000000001,NaN) pz(0.0050000000000000001,NaN) pz(1.0000000000000001e-005,NaN)], 'ZEROS', [pz(0.00050000000000000001,NaN) pz(0.001,NaN)]));
h13 = pzmodel(plist('GAIN', [0.00001], 'POLES', [pz(5.0000000000000004e-006,NaN) pz(0.02,NaN)], 'ZEROS', pz(0.00040000000000000002,NaN)));

h21 = pzmodel(plist('GAIN', [0.00001], 'POLES', [pz(5.0000000000000004e-006,NaN) pz(0.10000000000000001,NaN) pz(0.050000000000000003,NaN)], 'ZEROS', [pz(0.002,NaN) pz(0.00020000000000000001,NaN)]));
h22 = pzmodel(plist('GAIN', [0.001], 'POLES', [pz(0.01,NaN) pz(9.9999999999999995e-008,NaN) pz(0.002,NaN) pz(0.001,NaN)], 'ZEROS', [pz(1.0000000000000001e-005,NaN) pz(2.0000000000000002e-005,NaN) pz(5.0000000000000002e-005,NaN) pz(0.20000000000000001,NaN)]));
h23 = pzmodel(plist('GAIN', [0.0001], 'POLES', [pz(0.01,NaN) pz(0.029999999999999999,NaN) pz(0.10000000000000001,NaN) pz(5.0000000000000002e-005,NaN)], 'ZEROS', [pz(0.0001,NaN) pz(0.0050000000000000001,NaN)]));

h31 = pzmodel(plist('GAIN', [0.001], 'POLES', [pz(0.01,NaN) pz(0.029999999999999999,NaN) pz(0.10000000000000001,NaN) pz(0.00050000000000000001,NaN)], 'ZEROS', [pz(0.0001,NaN) pz(0.0050000000000000001,NaN)]));
h32 = pzmodel(plist('GAIN', [0.00001], 'POLES', [pz(0.01,NaN) pz(0.029999999999999999,NaN) pz(0.10000000000000001,NaN) pz(0.00050000000000000001,NaN) pz(0.00029999999999999997,NaN)], 'ZEROS', [pz(0.002,NaN) pz(0.059999999999999998,NaN)]));
h33 = pzmodel(plist('GAIN', [0.05], 'POLES', [pz(0.0001,NaN) pz(9.9999999999999995e-008,NaN) pz(0.002,NaN)], 'ZEROS', [pz(3.0000000000000001e-005,NaN) pz(5.0000000000000002e-005,NaN) pz(0.10000000000000001,NaN)]));

%%% TFs resps
% filters response
plresp = plist('f',f);

rh11 = resp(h11,plresp);
rh12 = resp(h12,plresp);
rh13 = resp(h13,plresp);

rh21 = resp(h21,plresp);
rh22 = resp(h22,plresp);
rh23 = resp(h23,plresp);

rh31 = resp(h31,plresp);
rh32 = resp(h32,plresp);
rh33 = resp(h33,plresp);

%%% Spectra calculation with Kay style method
% Note that the definition of cross-spectral matrix follows that of s M
% Kay, Modern Spectral Estimation

% csd11 
G11 = rh11.y.*conj(rh11.y) + rh12.y.*conj(rh12.y) + rh13.y.*conj(rh13.y);
G11 = ao(plist('xvals', f, 'yvals', G11, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
G11.setName;

% csd12
G12 = conj(rh11.y).*rh21.y + conj(rh12.y).*rh22.y + conj(rh13.y).*rh23.y;
G12 = ao(plist('xvals', f, 'yvals', G12, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
G12.setName;

% csd13
G13 = conj(rh11.y).*rh31.y + conj(rh12.y).*rh32.y + conj(rh13.y).*rh33.y;
G13 = ao(plist('xvals', f, 'yvals', G13, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
G13.setName;

% csd21
G21 = conj(G12);
G21.setName;

% csd22
G22 = rh21.y.*conj(rh21.y) + rh22.y.*conj(rh22.y) + rh23.y.*conj(rh23.y);
G22 = ao(plist('xvals', f, 'yvals', G22, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
G22.setName;

% csd23
G23 = conj(rh21.y).*rh31.y + conj(rh22.y).*rh32.y + conj(rh23.y).*rh33.y;
G23 = ao(plist('xvals', f, 'yvals', G23, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
G23.setName;

% csd31
G31 = conj(G13);
G31.setName;

% csd32
G32 = conj(G23);
G32.setName;

% csd33
G33 = conj(rh31.y).*rh31.y + conj(rh32.y).*rh32.y + conj(rh33.y).*rh33.y;
G33 = ao(plist('xvals', f, 'yvals', G33, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
G33.setName;

%%% Build CSD matrix
mod3D = [G11 G12 G13;G21 G22 G23;G31 G32 G33];

%% noise generation

% plist for noise generation

plns = plist(...
  'csd',       mod3D, ...
  'Nsecs',     1e4, ...
  'fs',        fs, ...
  'xunits',    's', ...
  'yunits',    '', ...
  'MaxIter',   90, ...
  'PoleType',  3, ...
  'MinOrder',  20, ...
  'MaxOrder',  50, ...
  'Weights',   2, ...
  'Plot',      false,...
  'MSEVARTOL', 1e-2,...
  'FITTOL',    1e-2);

na = ao(plns);

%% Make psd

naxx1 = na(1).psd;

naxx2 = na(2).psd;

naxx3 = na(3).psd;



iplot(naxx1,G11)
iplot(naxx2,G22)
iplot(naxx3,G33)