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