Mercurial > hg > ltpda
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')