Mercurial > hg > ltpda

comparison m-toolbox/test/utils/test2_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 2dim identification
3 %
4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 % HISTORY: 23-04-2009 L Ferraioli
6 % Creation
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 %% VERSION
9
10 '$Id: test2_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 % psd1 = mpsd(:,2);
22 % csd12 = mpsd(:,3);
23 % psd2 = mpsd(:,4);
24
25 csd = zeros(2,2,numel(f));
26 for ii = 1:numel(f)
27 csd(:,:,ii) = [mpsd(ii,2) mpsd(ii,3);conj(mpsd(ii,3)) mpsd(ii,4)];
28 end
29
30 %% identification in z-domain
31
32 % identification params
33 params = struct('Nmaxiter',70, 'minorder',20,...
34 'maxorder',40, 'spolesopt',2, 'weightparam',2, 'plot',0,...
35 'ctp','chival','lrscond',4,'msevar',2,...
36 'fs',10,'dterm',0, 'spy',0, 'fullauto',1,...
37 'extweights', []);
38
39 ostruct = utils.math.csd2tf(csd,f,params);
40
41 %% Model response TF11
42
43 pfparams.type = 'disc';
44 pfparams.freq = f;
45 pfparams.fs = 10;
46 pfparams.res = ostruct(1).res(:,1);
47 pfparams.pol = ostruct(1).poles(:,1);
48 pfparams.dterm = ostruct(1).dterm(:,1);
49
50 % response of the model filter
51 pfr = utils.math.pfresp(pfparams);
52 y11 = pfr.resp;
53
54 %% Model response TF12
55
56 pfparams.type = 'disc';
57 pfparams.freq = f;
58 pfparams.fs = 10;
59 pfparams.res = ostruct(2).res(:,1);
60 pfparams.pol = ostruct(2).poles(:,1);
61 pfparams.dterm = ostruct(2).dterm(:,1);
62
63 % response of the model filter
64 pfr = utils.math.pfresp(pfparams);
65 y12 = pfr.resp;
66
67 %% Model response TF21
68
69 pfparams.type = 'disc';
70 pfparams.freq = f;
71 pfparams.fs = 10;
72 pfparams.res = ostruct(1).res(:,2);
73 pfparams.pol = ostruct(1).poles(:,1);
74 pfparams.dterm = ostruct(1).dterm(:,2);
75
76 % response of the model filter
77 pfr = utils.math.pfresp(pfparams);
78 y21 = pfr.resp;
79
80 %% Model response TF22
81
82 pfparams.type = 'disc';
83 pfparams.freq = f;
84 pfparams.fs = 10;
85 pfparams.res = ostruct(2).res(:,2);
86 pfparams.pol = ostruct(2).poles(:,1);
87 pfparams.dterm = ostruct(2).dterm(:,2);
88
89 % response of the model filter
90 pfr = utils.math.pfresp(pfparams);
91 y22 = pfr.resp;
92
93 %% Model PSD calculation
94
95 % psd1(f) = tf11(f)*conj(tf11(f))+tf12(f)*conj(tf12(f))
96 % csd(f) = tf11(f)*conj(tf21(f))+tf12(f)*conj(tf22(f))
97 % psd2(f) = tf21(f)*conj(tf21(f))+tf22(f)*conj(tf22(f))
98
99 % mpsd1 = y11.*conj(y11)+y12.*conj(y12);
100 % mcsd = y11.*conj(y21)+y12.*conj(y22);
101 % mpsd2 = y21.*conj(y21)+y22.*conj(y22);
102
103 mpsd1 = conj(y11).*y11+conj(y12).*y12;
104 mcsd = conj(y11).*y21+conj(y12).*y22;
105 mpsd2 = conj(y21).*y21+conj(y22).*y22;
106
107 %% Checking result
108
109 figure()
110 loglog(f,squeeze(csd(1,1,:)),'k')
111 hold on
112 loglog(f,mpsd1,'r')
113 xlabel('Frequency [Hz]')
114 ylabel('psd')
115 legend('Original', 'Model')
116
117 figure()
118 loglog(f,abs(squeeze(csd(1,2,:))),'k')
119 hold on
120 loglog(f,abs(mcsd),'r')
121 xlabel('Frequency [Hz]')
122 ylabel('psd')
123 legend('Original', 'Model')
124
125 figure()
126 loglog(f,squeeze(csd(2,2,:)),'k')
127 hold on
128 loglog(f,mpsd2,'r')
129 xlabel('Frequency [Hz]')
130 ylabel('psd')
131 legend('Original', 'Model')
132
133 figure()
134 semilogy(ostruct(1).mse)
135 xlabel('Iteration')
136 ylabel('MSE')
137 legend('Mean Squared Error Progression for filters 11 and 21')
138
139 figure()
140 semilogy(ostruct(2).mse)
141 xlabel('Iteration')
142 ylabel('MSE')
143 legend('Mean Squared Error Progression for filters 12 and 22')