Mercurial > hg > ltpda

comparison m-toolbox/test/utils/test_2dim_vdfit.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 Fitting procedure in 2dim z-domain VDFIT
3 % At the end you have 4 stable transfer functions in partial fractions
4 %
5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 % HISTORY: 02-10-2008 L Ferraioli
7 % Creation
8 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
9 %% VERSION
10
11 '$Id: test_2dim_vdfit.m,v 1.1 2009/04/23 10:11:26 luigi Exp $';
12
13 %% Clear
14
15 clear all
16
17 %% Loading spectra
18
19 load ..\m-toolbox\test\mpsd.mat % load mpsd.mat first column is f then psd1, csd and psd2
20
21 f = mpsd(:,1);
22 csd11 = mpsd(:,2);
23 csd12 = mpsd(:,3);
24 csd21 = [];
25 csd22 = mpsd(:,4);
26
27 fs = 10;
28
29
30 %% Eigendecomposition
31
32 [tf11,tf12,tf21,tf22] = utils.math.eigcsd(csd11,csd12,csd21,csd22,'USESYM',0,'DIG',50,'OTP','TF');
33
34 %% Constructing vector
35
36 f1 = [tf11 tf21];
37 f2 = [tf12 tf22];
38
39 %% VDFIT
40 tic;
41
42 clear mlr1 mlr2
43
44 N = 32; %Order of approximation
45
46
47 % Max Iteration
48 Nmaxiter = 40;
49
50 % mlr1 = zeros(Nmaxiter,1);
51 % mlr2 = zeros(Nmaxiter,1);
52
53 % Fitting params
54 fitin.stable = 0;
55 fitin.dterm = 0;
56 fitin.plot = 0;
57 fitin.fs = fs;
58
59 weight = utils.math.wfun(f1,2);
60 % weight = 1./abs(f1);
61
62 pparams = struct('spolesopt',3, 'type','DISC', 'pamp', 0.98);
63 poles1 = utils.math.startpoles(N,f,pparams);
64 for hh = 1:Nmaxiter
65 [res1,poles1,dterm1,mresp1,rdl1,sqe1] = utils.math.vdfit(f1,f,poles1,weight,fitin); % Fitting
66 disp(['Iter' num2str(hh)])
67
68 % Stop condition checking
69 mlr1(hh) = sqe1(:,1);
70 mlr2(hh) = sqe1(:,2);
71 % [ext1,msg1] = utils.math.stopfit(f1(:,1),rdl1(:,1),mlr1,'lrsrmse',2,15);
72 % [ext2,msg2] = utils.math.stopfit(f1(:,2),rdl1(:,2),mlr2,'lrsrmse',2,15);
73
74 % order = length(poles1);
75 % mlr1(hh) = mean(log10(abs(f1(:,1)))-log10(abs(rdl1(:,1))));
76 % mlr2(hh) = mean(log10(abs(f1(:,2)))-log10(abs(rdl1(:,2))));
77 % [ext1,msg1] = utils.math.stopfit(f1(:,1),rdl1(:,1),mlr1,order,'lrs',1,0.0001,2);
78 % [ext2,msg2] = utils.math.stopfit(f1(:,2),rdl1(:,2),mlr2,order,'lrs',1,0.0001,2);
79
80 % if ext1 && ext2
81 % disp(msg1)
82 % break
83 % end
84
85 end
86 elpstime = toc;
87
88 % plotting squared error
89 figure()
90 semilogy(mlr1,'-ok')
91 hold on
92 grid on
93 semilogy(mlr2,'-or')
94
95 % % plotting RMS error variation
96 % figure()
97 % semilogy(abs(diff(sqrt(mlr1))),'-ok')
98 % hold on
99 % grid on
100 % semilogy(abs(diff(sqrt(mlr2))),'-or')
101
102 % plotting squared error variation
103 figure()
104 semilogy(abs(diff(mlr1)./mlr1(1,end-1)),'-ok')
105 hold on
106 grid on
107 semilogy(abs(diff(mlr2)./mlr2(1,end-1)),'-or')
108
109 %% All Passing
110
111 % [np1,resp1] = pfallpz(res1,poles1,dterm1,mresp1,f,fs);
112 [nr1,np1,nd1,resp1] = utils.math.pfallpsymz(res1,poles1,dterm1,mresp1,f,fs);
113
114
115 figure()
116 subplot(2,1,1);
117 p1 = loglog(f,abs(mresp1),'k');
118 hold on
119 p2 = loglog(f,abs(resp1),'r');
120 xlabel('Frequency [Hz]')
121 ylabel('Amplitude')
122 legend([p1(1) p2(1)],'VDFIT','Stabilized')
123 hold off
124
125 subplot(2,1,2);
126 p4 = semilogx(f,(180/pi).*unwrap(angle(mresp1)),'k');
127 hold on
128 p5 = semilogx(f,(180/pi).*unwrap(angle(resp1)),'r');
129 xlabel('Frequency [Hz]')
130 ylabel('Phase [Deg]')
131 legend([p4(1) p5(1)],'VDFIT', 'Stabilized')
132 hold off
133
134 %% VDFIT
135
136 clear mlr3 mlr4 ext1 ext2
137 N = 24; %Order of approximation
138
139
140 % Max Iteration
141 Nmaxiter = 70;
142
143 % Fitting params
144 fitin.stable = 0;
145 fitin.dterm = 0;
146 fitin.plot = 1;
147 fitin.fs = fs;
148
149 weight = utils.math.wfun(f2,2);
150
151 pparams = struct('spolesopt',2, 'type','DISC', 'pamp', 0.98);
152 poles2 = utils.math.startpoles(N,f,pparams);
153
154 for hh = 1:Nmaxiter
155 [res2,poles2,dterm2,mresp2,rdl2,rmse2] = utils.math.vdfit(f2,f,poles2,weight,fitin); % Fitting
156 disp(['Iter' num2str(hh)])
157
158 % Stop condition checking
159 mlr3(hh) = rmse2(:,1);
160 mlr4(hh) = rmse2(:,2);
161 [ext1,msg1] = utils.math.stopfit(f2(:,1),rdl2(:,1),mlr3,'lrsrmse',2,15);
162 [ext2,msg2] = utils.math.stopfit(f2(:,2),rdl2(:,2),mlr4,'lrsrmse',2,15);
163
164 if ext1 && ext2
165 disp(msg1)
166 break
167 end
168
169 end
170
171 % plotting RMS error
172 figure()
173 semilogy(mlr3,'-ok')
174 hold on
175 grid on
176 semilogy(mlr4,'-or')
177
178 %%
179 % params = struct('spolesopt',2, 'Nmaxiter',250, 'minorder',30,...
180 % 'maxorder',60, 'weightparam',1, 'plot',0,...
181 % 'ctp','lrs','lrscond',1,'lrsvarcond',1,'nrmsecond',2,...
182 % 'stabfit',0,'dterm',0,'spy',0);
183 %
184 % [res2,poles2,dterm2,mresp2,rdl2] = autodfit(f2,f,fs,params);