% A test script for utils.math.eigcsd
% 
% L. Ferraioli 18-11-08
% 
% $Id: test_eigcsd.m,v 1.2 2009/06/10 16:12:06 luigi Exp $
% 

%%

clear all

%% CSD Noise models and TF

[CSD,TF] = get_2D_test_obj_tf_psd();

% get noise coherence
% coh = CSD(1,2)./sqrt(abs(CSD(1,1).*CSD(2,2)));
% rcoh = real(coh);
% rcoh.setName;
% icoh = imag(coh);
% icoh.setName;

%% Get Noisegen TFs numerically

[tf11,tf12,tf21,tf22] = utils.math.eigcsd(CSD(1,1).y,CSD(1,2).y,CSD(2,1).y,CSD(2,2).y,'USESYM',0,'DIG',50,'OTP','TF');

% Build AOs
h11 = ao(plist('xvals', TF(1,1).x, 'yvals', tf11, 'fs', TF(1,1).fs, 'dtype', 'fsdata'));
h11.setName;
h12 = ao(plist('xvals', TF(1,1).x, 'yvals', tf12, 'fs', TF(1,1).fs, 'dtype', 'fsdata'));
h12.setName;
h21 = ao(plist('xvals', TF(1,1).x, 'yvals', tf21, 'fs', TF(1,1).fs, 'dtype', 'fsdata'));
h21.setName;
h22 = ao(plist('xvals', TF(1,1).x, 'yvals', tf22, 'fs', TF(1,1).fs, 'dtype', 'fsdata'));
h22.setName;


%% Checking generated TFs

iplot(abs(abs(TF(1,1))-abs(h11)))
iplot(abs(abs(TF(1,2))-abs(h12)))
iplot(abs(abs(TF(2,1))-abs(h21)))
iplot(abs(abs(TF(2,2))-abs(h22)))

%% Calculate psd

xx11 = h11.*conj(h11)+h12.*conj(h12);
xx12 = h11.*conj(h21)+h12.*conj(h22);
xx22 = h22.*conj(h22)+h21.*conj(h21);
xx21 = conj(xx12);
xx11.setName;
xx12.setName;
xx21.setName;
xx22.setName;

%% Compare

iplot(CSD(1,1),xx11)
iplot(CSD(1,2),xx12)
iplot(CSD(2,1),xx21)
iplot(CSD(2,2),xx22)

%% Get Noisegen TFs symbolically

[tf11,tf12,tf21,tf22] = utils.math.eigcsd(CSD(1,1).y,CSD(1,2).y,CSD(2,1).y,CSD(2,2).y,'USESYM',1,'DIG',50,'OTP','TF');

% Build AOs
h11 = ao(plist('xvals', TF(1,1).x, 'yvals', tf11, 'fs', TF(1,1).fs, 'dtype', 'fsdata'));
h11.setName;
h12 = ao(plist('xvals', TF(1,1).x, 'yvals', tf12, 'fs', TF(1,1).fs, 'dtype', 'fsdata'));
h12.setName;
h21 = ao(plist('xvals', TF(1,1).x, 'yvals', tf21, 'fs', TF(1,1).fs, 'dtype', 'fsdata'));
h21.setName;
h22 = ao(plist('xvals', TF(1,1).x, 'yvals', tf22, 'fs', TF(1,1).fs, 'dtype', 'fsdata'));
h22.setName;


%% Checking generated TFs

iplot(abs(abs(TF(1,1))-abs(h11)))
iplot(abs(abs(TF(1,2))-abs(h12)))
iplot(abs(abs(TF(2,1))-abs(h21)))
iplot(abs(abs(TF(2,2))-abs(h22)))

%% Calculate psd

xx11 = h11.*conj(h11)+h12.*conj(h12);
xx12 = h11.*conj(h21)+h12.*conj(h22);
xx22 = h22.*conj(h22)+h21.*conj(h21);
xx21 = conj(xx12);
xx11.setName;
xx12.setName;
xx21.setName;
xx22.setName;

%% Compare

iplot(CSD(1,1),xx11)
iplot(CSD(1,2),xx12)
iplot(CSD(2,1),xx21)
iplot(CSD(2,2),xx22)