% A test script for ao/noisegen2D
% test ao/noisegen2D with input the model and output the filter
% 
% DESCRIPTION: Run noisegen2D with test data and test procedure accuracy
% 
% L. Ferraioli 10-11-08
% 
% $Id: test_ao_noisegen2D_2.m,v 1.1 2010/05/03 19:04:45 luigi Exp $
% 

%% General use variables and vectors

userdir = 'C:\Users\Luigi'; % You should set your own dir

f = logspace(-6,log10(5),300);
fs = 10;

Nt = fs*10000;


%% CSD Noise models

fundir = fullfile(userdir,'ltp_data_analysis\MDCs\MDC1_UTN');
cf = cd;
cd(fundir);
[TF,CSD] = mdc1_tf_models(plist('f',f,'fs',fs));
cd(cf);


%% Plot Model Tfs

iplot(TF(1,1),TF(1,2),TF(2,1),TF(2,2))

%% Plot Model CSD

iplot(CSD(1,1),CSD(1,2),CSD(2,1),CSD(2,2))

%% Noise generation - 3 inputs

pl = plist(...
    'fs', fs, ...
    'Iunits', '', ...
    'Ounits', 'm', ...
    'MaxIter', 80, ...
    'PoleType', 2, ...
    'MinOrder', 15, ...
    'MaxOrder', 45, ...
    'Weights', 3, ...
    'FITTOL', 1e-3,...
    'MSEVARTOL', 1e-1,...
    'UseSym', 0,...
    'Plot', false,...
    'Disp', false);

fil = noisegen2D(CSD(1,1),CSD(1,2),CSD(2,2), pl);

%% filtering data

a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e5));
a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e5));

a = matrix(a1,a2,plist('shape',[2 1]));

b = filter(a,fil);

%% Checking results and starting data

iplot(b.objs)

%% Make psd

plpsd = plist('navs',8,'order',1,'olap',50);
bxx1 = b.objs(1).psd(plpsd);
bxx2 = b.objs(2).psd(plpsd);

iplot(bxx1,CSD(1,1),bxx2,CSD(2,2))

%% Noise generation - 4 inputs

pl = plist(...
    'fs', fs, ...
    'Iunits', '', ...
    'Ounits', 'm', ...
    'MaxIter', 80, ...
    'PoleType', 2, ...
    'MinOrder', 15, ...
    'MaxOrder', 45, ...
    'Weights', 3, ...
    'FITTOL', 1e-3,...
    'MSEVARTOL', 1e-1,...
    'UseSym', 0,...
    'Plot', false,...
    'Disp', false);

fil = noisegen2D(CSD(1,1),CSD(1,2),CSD(2,1),CSD(2,2), pl);

%% filtering data

a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e5));
a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e5));

a = matrix(a1,a2,plist('shape',[2 1]));

b = filter(a,fil);

%% Checking results and starting data

iplot(b.objs)

%% Make psd

plpsd = plist('navs',8,'order',1,'olap',50);
bxx1 = b.objs(1).psd(plpsd);
bxx2 = b.objs(2).psd(plpsd);

iplot(bxx1,CSD(1,1),bxx2,CSD(2,2))