Mercurial > hg > ltpda

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % TEST csd2tf 1dim identification
3 %
4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 % HISTORY: 23-04-2009 L Ferraioli
6 % Creation
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 %% VERSION
9
10 '$Id: test1_csd2tf.m,v 1.1 2009/04/23 10:11:26 luigi Exp $';
11
12 %% Clear
13
14 clear all
15
16 %% Loading spectra
17
18 load ..\m-toolbox\test\mpsd.mat % load mpsd.mat first column is f then psd1, csd and psd2
19
20 f = mpsd(:,1);
21 psd = mpsd(:,2);
22
23 %% identification in z-domain
24
25 % identification params
26 params = struct('Nmaxiter',50, 'minorder',10,...
27 'maxorder',25, 'spolesopt',2, 'weightparam',2, 'plot',0,...
28 'ctp','chival','lrscond',3,'msevar',2,...
29 'fs',10,'dterm',0, 'spy',0, 'fullauto',1,...
30 'extweights', []);
31
32 ostruct = utils.math.csd2tf(psd,f,params);
33
34 %% Model response
35
36 res = ostruct.res;
37 poles = ostruct.poles;
38 dterm = ostruct.dterm;
39
40 pfparams.type = 'disc';
41 pfparams.freq = f;
42 pfparams.fs = 10;
43 pfparams.res = res;
44 pfparams.pol = poles;
45 pfparams.dterm = dterm;
46
47 % response of the model filter
48 pfr = utils.math.pfresp(pfparams);
49 y = pfr.resp;
50
51 %% Checking result
52
53 figure()
54 loglog(f,psd,'k')
55 hold on
56 loglog(f,y.*conj(y),'r')
57 xlabel('Frequency [Hz]')
58 ylabel('psd')
59 legend('Original', 'Model')
60
61 figure()
62 loglog(f,abs(y),'k')
63 hold on
64 xlabel('Frequency [Hz]')
65 ylabel('Amplitude')
66 legend('TF Model')
67
68 figure()
69 semilogx(f,(180/pi).*unwrap(angle(y)),'k')
70 hold on
71 xlabel('Frequency [Hz]')
72 ylabel('Phase Angle [Deg]')
73 legend('TF Model')
74
75 figure()
76 semilogy(ostruct.mse)
77 xlabel('Iteration')
78 ylabel('MSE')
79 legend('Mean Squared Error Progression')