% A test script for ao/whiten2D
% 
% L. Ferraioli 18-11-08
% 
% $Id: test_ao_whiten2D.m,v 1.6 2009/04/20 14:31:14 luigi Exp $
% 

%%

clear all

%% CSD Noise models

% Mdc2 noise model
% mod = ltpnoise();

load mpsd.mat % load mpsd.mat first column is f then csd11, csd12 and csd22

f = mpsd(:,1); % frequency vector
csd11 = mpsd(:,2);
csd12 = mpsd(:,3);
csd22 = mpsd(:,4);
fs = 10; % sampling frequency

% Building AOs
ltp_csd11 = ao(fsdata(f,csd11,fs));
ltp_csd12 = ao(fsdata(f,csd12,fs));
ltp_csd21 = ao(fsdata(f,conj(csd12),fs));
ltp_csd22 = ao(fsdata(f,csd22,fs));

%% 

iplot(ltp_csd11,ltp_csd12,ltp_csd22)

%% Colored Noise generation

% Make white noise

a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 10000));
a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 10000));
a = [a1 a2];

% Generation from white noise

pl = plist(...
    'csd11', ltp_csd11,...
    'csd12', ltp_csd12,...
    'csd21', ltp_csd21,...
    'csd22', ltp_csd22,...
    'MaxIter', 70, ...
    'PoleType', 2, ...
    'MinOrder', 20, ...
    'MaxOrder', 30, ...
    'Weights', 2, ...
    'Plot', false,...
    'MSEVARTOL', 1e-2,...
    'FITTOL', 1e-3,...
    'UseSym', 0,...
    'Disp', false);

ac = noisegen2D(a, pl);

%% Checking generated noise

acxx1 = ac(1).psd;
acxx2 = ac(2).psd;
acxx12 = ac.cpsd;

iplot(abs(acxx1),abs(ltp_csd11))
iplot(abs(acxx12),abs(ltp_csd12))
iplot(abs(acxx2),abs(ltp_csd22))

%% Noise whitening

pl = plist(...
    'csd11', ltp_csd11, ...
    'csd12', ltp_csd12, ...
    'csd21', ltp_csd21, ...
    'csd22', ltp_csd22, ...
    'MaxIter', 75, ...
    'PoleType', 3, ...
    'MinOrder', 20, ...
    'MaxOrder', 40, ...
    'Weights', 2, ...
    'Plot', false,...
    'MSEVARTOL', 1e-2,...
    'FITTOL', 1e-3,...
    'UseSym', 0,...
    'Disp', false);

aw = whiten2D(ac, pl);

%% Checking whitened noise

awxx1 = aw(1).psd;
awxx2 = aw(2).psd;

%%
iplot(awxx1,awxx2)

%% Checking

aw1 = aw(1);
aw2 = aw(2);
awxx1 = awxx(1,1);
awxx2 = awxx(2,2);

mw1 = aw1.mean; % 0 expected
vw1 = aw1.var; % 1 expected

mw2 = aw2.mean; % 0 expected
vw2 = aw2.var; % 1 expected

mwxx1 = awxx1.mean; % 0.2 expected
vwxx1 = awxx1.var;

mwxx2 = awxx2.mean; % 0.2 expected
vwxx2 = awxx2.var;

% correlation coefficients
% diagonal terms are aw1 and aw2 variances
% cross terms are data series covariances, a value close to 0 means no
% correlation between the data series
crc = corrcoef([aw1.data.y.' aw2.data.y.']);


% END