diff m-toolbox/test/utils/test1_csd2tf.m @ 0:f0afece42f48

Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/utils/test1_csd2tf.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,79 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% TEST csd2tf 1dim identification
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% HISTORY:     23-04-2009 L Ferraioli
+%                 Creation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% VERSION
+
+'$Id: test1_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);
+psd = mpsd(:,2);
+
+%% identification in z-domain
+
+% identification params
+params = struct('Nmaxiter',50, 'minorder',10,...
+    'maxorder',25, 'spolesopt',2, 'weightparam',2, 'plot',0,...
+    'ctp','chival','lrscond',3,'msevar',2,...
+    'fs',10,'dterm',0, 'spy',0, 'fullauto',1,...
+    'extweights', []);
+
+ostruct = utils.math.csd2tf(psd,f,params);
+
+%% Model response
+
+res = ostruct.res;
+poles = ostruct.poles;
+dterm = ostruct.dterm;
+
+pfparams.type = 'disc';
+pfparams.freq = f;
+pfparams.fs = 10;
+pfparams.res = res;
+pfparams.pol = poles;
+pfparams.dterm = dterm;
+
+% response of the model filter
+pfr = utils.math.pfresp(pfparams);
+y = pfr.resp;
+
+%% Checking result
+
+figure()
+loglog(f,psd,'k')
+hold on
+loglog(f,y.*conj(y),'r')
+xlabel('Frequency [Hz]')
+ylabel('psd')
+legend('Original', 'Model')
+
+figure()
+loglog(f,abs(y),'k')
+hold on
+xlabel('Frequency [Hz]')
+ylabel('Amplitude')
+legend('TF Model')
+
+figure()
+semilogx(f,(180/pi).*unwrap(angle(y)),'k')
+hold on
+xlabel('Frequency [Hz]')
+ylabel('Phase Angle [Deg]')
+legend('TF Model')
+
+figure()
+semilogy(ostruct.mse)
+xlabel('Iteration')
+ylabel('MSE')
+legend('Mean Squared Error Progression')