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