--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/utils/test2_csd2tf.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,143 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% TEST csd2tf 2dim identification
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% HISTORY:     23-04-2009 L Ferraioli
+%                 Creation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% VERSION
+
+'$Id: test2_csd2tf.m,v 1.1 2009/04/23 10:11:26 luigi Exp $';
+
+%% Clear
+
+clear all
+
+%% Loading spectra
+
+load ..\m-toolbox\test\mpsd.mat % load mpsd.mat first column is f then psd1, csd and psd2
+
+f = mpsd(:,1);
+% psd1 = mpsd(:,2);
+% csd12 = mpsd(:,3);
+% psd2 = mpsd(:,4);
+
+csd = zeros(2,2,numel(f));
+for ii = 1:numel(f)
+  csd(:,:,ii) = [mpsd(ii,2) mpsd(ii,3);conj(mpsd(ii,3)) mpsd(ii,4)];
+end
+
+%% identification in z-domain
+
+% identification params
+params = struct('Nmaxiter',70, 'minorder',20,...
+    'maxorder',40, 'spolesopt',2, 'weightparam',2, 'plot',0,...
+    'ctp','chival','lrscond',4,'msevar',2,...
+    'fs',10,'dterm',0, 'spy',0, 'fullauto',1,...
+    'extweights', []);
+
+ostruct = utils.math.csd2tf(csd,f,params);
+
+%% Model response TF11
+
+pfparams.type = 'disc';
+pfparams.freq = f;
+pfparams.fs = 10;
+pfparams.res = ostruct(1).res(:,1);
+pfparams.pol = ostruct(1).poles(:,1);
+pfparams.dterm = ostruct(1).dterm(:,1);
+
+% response of the model filter
+pfr = utils.math.pfresp(pfparams);
+y11 = pfr.resp;
+
+%% Model response TF12
+
+pfparams.type = 'disc';
+pfparams.freq = f;
+pfparams.fs = 10;
+pfparams.res = ostruct(2).res(:,1);
+pfparams.pol = ostruct(2).poles(:,1);
+pfparams.dterm = ostruct(2).dterm(:,1);
+
+% response of the model filter
+pfr = utils.math.pfresp(pfparams);
+y12 = pfr.resp;
+
+%% Model response TF21
+
+pfparams.type = 'disc';
+pfparams.freq = f;
+pfparams.fs = 10;
+pfparams.res = ostruct(1).res(:,2);
+pfparams.pol = ostruct(1).poles(:,1);
+pfparams.dterm = ostruct(1).dterm(:,2);
+
+% response of the model filter
+pfr = utils.math.pfresp(pfparams);
+y21 = pfr.resp;
+
+%% Model response TF22
+
+pfparams.type = 'disc';
+pfparams.freq = f;
+pfparams.fs = 10;
+pfparams.res = ostruct(2).res(:,2);
+pfparams.pol = ostruct(2).poles(:,1);
+pfparams.dterm = ostruct(2).dterm(:,2);
+
+% response of the model filter
+pfr = utils.math.pfresp(pfparams);
+y22 = pfr.resp;
+
+%% Model PSD calculation
+
+%     psd1(f) = tf11(f)*conj(tf11(f))+tf12(f)*conj(tf12(f))
+%     csd(f) = tf11(f)*conj(tf21(f))+tf12(f)*conj(tf22(f))
+%     psd2(f) = tf21(f)*conj(tf21(f))+tf22(f)*conj(tf22(f))
+
+% mpsd1 = y11.*conj(y11)+y12.*conj(y12);
+% mcsd = y11.*conj(y21)+y12.*conj(y22);
+% mpsd2 = y21.*conj(y21)+y22.*conj(y22);
+
+mpsd1 = conj(y11).*y11+conj(y12).*y12;
+mcsd = conj(y11).*y21+conj(y12).*y22;
+mpsd2 = conj(y21).*y21+conj(y22).*y22;
+
+%% Checking result
+
+figure()
+loglog(f,squeeze(csd(1,1,:)),'k')
+hold on
+loglog(f,mpsd1,'r')
+xlabel('Frequency [Hz]')
+ylabel('psd')
+legend('Original', 'Model')
+
+figure()
+loglog(f,abs(squeeze(csd(1,2,:))),'k')
+hold on
+loglog(f,abs(mcsd),'r')
+xlabel('Frequency [Hz]')
+ylabel('psd')
+legend('Original', 'Model')
+
+figure()
+loglog(f,squeeze(csd(2,2,:)),'k')
+hold on
+loglog(f,mpsd2,'r')
+xlabel('Frequency [Hz]')
+ylabel('psd')
+legend('Original', 'Model')
+
+figure()
+semilogy(ostruct(1).mse)
+xlabel('Iteration')
+ylabel('MSE')
+legend('Mean Squared Error Progression for filters 11 and 21')
+
+figure()
+semilogy(ostruct(2).mse)
+xlabel('Iteration')
+ylabel('MSE')
+legend('Mean Squared Error Progression for filters 12 and 22')