Mercurial > hg > ltpda

comparison m-toolbox/test/utils/test_autodfit.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 script for autodfit
3 %
4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 % HISTORY: 12-09-2008 L Ferraioli
6 % Creation
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 %% VERSION
9
10 '$Id: test_autodfit.m,v 1.1 2009/04/23 10:11:26 luigi Exp $';
11
12 %% Cleaning
13
14 clear all
15
16 %% data
17
18
19
20 % % frequencies
21 % f = logspace(-6,log10(5),300); % frequency vector in Hz
22 % f = f.';
23 % fs = 10;
24
25 % % Continuous rational transfer function
26 % Num = [-2.726e-007; 1.665e-005; 1.303e-007; 8.381e-010]; % Numerator
27 % Den = [1; 0.2189; 0.01922; 0.0007803; 0]; % Denominator
28 % y = freqs(Num,Den,w);
29
30 % % filt1
31 % dRes = [0.2+j*0.003;0.2-j*0.003;0.45+j*0.0007;0.45-j*0.0007];
32 % dPoles = [0.97+j*0.0003;0.97-j*0.0003;0.75+j*0.00005;0.75-j*0.00005];
33 % dDTerms = 0;
34
35 % % filt2
36 % dRes = [0.2+j*0.003;0.2-j*0.003;0.45+j*0.0007;0.45-j*0.0007;12;0.45+j*0.07;0.45-j*0.07];
37 % dPoles = [0.97+j*0.0003;0.97-j*0.0003;0.75+j*0.00005;0.75-j*0.00005;0.1;0.998+j*0.00005;0.998-j*0.00005];
38 % dDTerms = 0;
39
40 % % Model aStefano Tf11
41 % dRes = [2.44554138162509e-011 - 1.79482547894083e-011i;
42 % 2.44554138162509e-011 + 1.79482547894083e-011i;
43 % 2.66402334803101e-009 + 1.1025122049153e-009i;
44 % 2.66402334803101e-009 - 1.1025122049153e-009i;
45 % -7.3560293387644e-009;
46 % -1.82811618589835e-009 - 1.21803627800855e-009i;
47 % -1.82811618589835e-009 + 1.21803627800855e-009i;
48 % 1.16258677367555e-009;
49 % 1.65216557639319e-016;
50 % -1.78092396888606e-016;
51 % -2.80420398962379e-017;
52 % 9.21305973049041e-013 - 8.24686706827269e-014i;
53 % 9.21305973049041e-013 + 8.24686706827269e-014i;
54 % 5.10730060739905e-010 - 3.76571756625722e-011i;
55 % 5.10730060739905e-010 + 3.76571756625722e-011i;
56 % 3.45893698149735e-009;
57 % 3.98139182134446e-014 - 8.25503935419059e-014i;
58 % 3.98139182134446e-014 + 8.25503935419059e-014i;
59 % -1.40595719147164e-011];
60 %
61 % dPoles = [0.843464045655194 - 0.0959986292915475i;
62 % 0.843464045655194 + 0.0959986292915475i;
63 % 0.953187595424927 - 0.0190043625473383i;
64 % 0.953187595424927 + 0.0190043625473383i;
65 % 0.967176277937188;
66 % 0.995012027005247 - 0.00268322602801729i;
67 % 0.995012027005247 + 0.00268322602801729i;
68 % 0.996564761885673;
69 % 0.999999366165445;
70 % 0.999981722418555;
71 % 0.999921882627659;
72 % 0.999624431675213 - 0.000813407848742761i;
73 % 0.999624431675213 + 0.000813407848742761i;
74 % 0.997312006278751 - 0.00265611346834941i;
75 % 0.997312006278751 + 0.00265611346834941i;
76 % 0.990516544257531;
77 % 0.477796923118318 - 0.311064085401834i;
78 % 0.477796923118318 + 0.311064085401834i;
79 % 0];
80 %
81 % dDTerms = 0;
82
83 % % Model Stefano TF12
84 % dRes = [1.44258422208796e-017 + 7.07359428613009e-019i;
85 % 1.44258422208796e-017 - 7.07359428613009e-019i;
86 % -3.4918408053655e-021 - 1.05662874569329e-021i;
87 % -3.4918408053655e-021 + 1.05662874569329e-021i;
88 % -7.61773292876976e-021;
89 % 4.84357724603939e-020 + 2.38824204294595e-019i;
90 % 4.84357724603939e-020 - 2.38824204294595e-019i;
91 % -4.07088520945753e-020 - 2.31474543846105e-019i;
92 % -4.07088520945753e-020 + 2.31474543846105e-019i;
93 % 8.73316588658882e-023;
94 % -5.21840635377469e-020;
95 % 1.8461911504859e-023;
96 % 5.20105247464461e-020;
97 % -4.68960092394415e-022;
98 % -1.44261407664171e-017 + 6.8922564526833e-019i;
99 % -1.44261407664171e-017 - 6.8922564526833e-019i;
100 % 3.13688133935426e-022];
101 %
102 % dPoles = [0.477546340377332 - 0.310830571032376i;
103 % 0.477546340377332 + 0.310830571032376i;
104 % 0.99790715414307 - 0.0028490561287024i;
105 % 0.99790715414307 + 0.0028490561287024i;
106 % 0.998014205354671 ;
107 % 0.999585354543332 - 0.000780408757425194i;
108 % 0.999585354543332 + 0.000780408757425194i;
109 % 0.99966003029931 - 0.000830944038363768i;
110 % 0.99966003029931 + 0.000830944038363768i;
111 % 0.999962770401331 ;
112 % 0.999981881865521 ;
113 % 0.999999365763457 ;
114 % 0.999981706320212 ;
115 % 0.99992421574188 ;
116 % 0.477898460791003 + 0.311001926610074i;
117 % 0.477898460791003 - 0.311001926610074i;
118 % 0];
119 % dDTerms = 0;
120
121 % % Model Stefano Tf21
122 % dRes = [-1.80035241582968e-016 + 1.99543917791863e-015i;
123 % -1.80035241582968e-016 - 1.99543917791863e-015i;
124 % -1.85590889333759e-013 - 1.23844418827409e-014i;
125 % -1.85590889333759e-013 + 1.23844418827409e-014i;
126 % 5.03656596876842e-013 ;
127 % -2.62470963499904e-013 + 2.30024232938878e-012i;
128 % -2.62470963499904e-013 - 2.30024232938878e-012i;
129 % -9.83780507870955e-018 ;
130 % 3.40426735130194e-021 ;
131 % 9.78322351492755e-018 ;
132 % -1.65010934542937e-020 ;
133 % 2.60918565203438e-015 + 1.0546609464659e-015i;
134 % 2.60918565203438e-015 - 1.0546609464659e-015i;
135 % 3.18105585405455e-014 + 2.48839990780042e-013i;
136 % 3.18105585405455e-014 - 2.48839990780042e-013i;
137 % 3.23021641947666e-013 ;
138 % 4.81265000078114e-016 - 3.18269170053848e-017i;
139 % 4.81265000078114e-016 + 3.18269170053848e-017i;
140 % 5.16260024128201e-018];
141 %
142 % dPoles = [0.872004077421604 - 0.110344282822693i;
143 % 0.872004077421604 + 0.110344282822693i;
144 % 0.956884129232757 - 0.0225532091775074i;
145 % 0.956884129232757 + 0.0225532091775074i;
146 % 0.966514825697177 ;
147 % 0.995140550419744 - 0.000755127639524413i;
148 % 0.995140550419744 + 0.000755127639524413i;
149 % 0.999981802194393 ;
150 % 0.99999936576546 ;
151 % 0.999981722418555 ;
152 % 0.999921882627659 ;
153 % 0.999624431675213 + 0.000813407848742761i;
154 % 0.999624431675213 - 0.000813407848742761i;
155 % 0.997312006278751 + 0.00265611346834941i;
156 % 0.997312006278751 - 0.00265611346834941i;
157 % 0.990516544257531 ;
158 % 0.477796923118318 + 0.311064085401834i;
159 % 0.477796923118318 - 0.311064085401834i;
160 % 0];
161 %
162 % dDTerms = 0;
163
164 % % Model Stefano Tf22
165 % dRes = [1.1284521501259e-014;
166 % -3.72133611555879e-014 - 2.08232683444075e-014i;
167 % -3.72133611555879e-014 + 2.08232683444075e-014i;
168 % 9.84930639106637e-014 - 1.46640810672565e-013i;
169 % 9.84930639106637e-014 + 1.46640810672565e-013i;
170 % 2.51323684013671e-014 ;
171 % -5.64078525288305e-014 ;
172 % -1.62476406586366e-014 ;
173 % -1.08424815979566e-011 + 8.32328079357669e-012i;
174 % -1.08424815979566e-011 - 8.32328079357669e-012i;
175 % 2.41831559776112e-011];
176 %
177 % dPoles = [0.988511243978897;
178 % 0.997305870640646 + 0.00211760900132725i;
179 % 0.997305870640646 - 0.00211760900132725i;
180 % 0.999626453270255 + 0.0008125673525946i;
181 % 0.999626453270255 - 0.0008125673525946i;
182 % 0.999999366366222 ;
183 % 0.999981706320212 ;
184 % 0.99992421574188 ;
185 % 0.477898460791003 - 0.311001926610074i;
186 % 0.477898460791003 + 0.311001926610074i;
187 % 0];
188 %
189 % dDTerms = 0;
190
191 % % filt4
192 % dRes = [2.44554138162509e-011;
193 % 2.66402334803101e-009;
194 % -7.3560293387644e-009;
195 % -1.82811618589835e-009;
196 % 1.16258677367555e-009;
197 % 1.65216557639319e-016;
198 % -1.78092396888606e-016;
199 % -2.80420398962379e-017;
200 % 9.21305973049041e-013;
201 % 5.10730060739905e-010;
202 % 3.45893698149735e-009;
203 % 3.98139182134446e-014;
204 % -1.40595719147164e-011];
205 %
206 % dPoles = [0.843464045655194;
207 % 0.953187595424927;
208 % 0.967176277937188;
209 % 0.995012027005247;
210 % 0.996564761885673;
211 % 0.999999366165445;
212 % 0.999981722418555;
213 % 0.999921882627659;
214 % 0.999624431675213;
215 % 0.997312006278751;
216 % 0.990516544257531;
217 % 0.477796923118318;
218 % 0];
219 %
220 % dDTerms = 0;
221 %
222 %
223 % pfparams.type = 'disc';
224 % pfparams.freq = f;
225 % pfparams.fs = fs;
226 % pfparams.res = dRes;
227 % pfparams.pol = dPoles;
228 % pfparams.dterm = dDTerms;
229 %
230 % % response of the model filter
231 % pfr = pfresp(pfparams);
232 % y = pfr.resp;
233
234 % % Model Continuous minphase
235 % r = [1.4871879011802e-005;
236 % 6.10992134443159e-016 - 9.32963644121133e-014i;
237 % 6.10992134443159e-016 + 9.32963644121133e-014i;
238 % 6.5700947750201e-017 - 8.67535202691611e-015i;
239 % 6.5700947750201e-017 + 8.67535202691611e-015i;
240 % 7.60899607753932e-017 - 2.31523533550992e-015i;
241 % 7.60899607753932e-017 + 2.31523533550992e-015i;
242 % 7.81925532239946e-010 - 5.60617358881436e-009i;
243 % 7.81925532239946e-010 + 5.60617358881436e-009i;
244 % -1.36723085989008e-009 + 5.08414219753465e-009i;
245 % -1.36723085989008e-009 - 5.08414219753465e-009i;
246 % -8.98316119086409e-009 - 2.15634331132174e-007i;
247 % -8.98316119086409e-009 + 2.15634331132174e-007i;
248 % 9.22133132976409e-009 - 3.24455764655643e-008i;
249 % 9.22133132976409e-009 + 3.24455764655643e-008i;
250 % 6.27280610849718e-011 - 1.19002054635157e-009i;
251 % 6.27280610849718e-011 + 1.19002054635157e-009i;
252 % 3.28962376794087e-010 - 3.35707070704021e-010i;
253 % 3.28962376794087e-010 + 3.35707070704021e-010i];
254 %
255 % p = [-1878769.47366179;
256 % -9.61587854136845e-006 + 4.02759537557155e-008i;
257 % -9.61587854136845e-006 - 4.02759537557155e-008i;
258 % -3.97325785368196e-005 + 7.21480997398869e-007i;
259 % -3.97325785368196e-005 - 7.21480997398869e-007i;
260 % -0.00016489329503076 + 4.21265809470161e-006i;
261 % -0.00016489329503076 - 4.21265809470161e-006i;
262 % -0.0488243534364908 + 0.00506677723199764i;
263 % -0.0488243534364908 - 0.00506677723199764i;
264 % -0.060222623381404 + 0.0123923197320473i;
265 % -0.060222623381404 - 0.0123923197320473i;
266 % -0.404482556535086 + 0.0153085700658377i;
267 % -0.404482556535086 - 0.0153085700658377i;
268 % -0.955593480749454 + 0.0869402878111827i;
269 % -0.955593480749454 - 0.0869402878111827i;
270 % -2.56698459926376 + 0.591616326379104i;
271 % -2.56698459926376 - 0.591616326379104i;
272 % -4.4135552953192 + 0.652758967843539i;
273 % -4.4135552953192 - 0.652758967843539i];
274 %
275 % d = 0;
276 %
277 % f = logspace(-6,log10(5),300); % frequency vector in Hz
278 % f = f.';
279 % pfparams.type = 'cont';
280 % pfparams.freq = f;
281 % pfparams.res = r;
282 % pfparams.pol = p;
283 % pfparams.dterm = d;
284 %
285 % % response of the model filter
286 % pfr = pfresp(pfparams);
287 % y = pfr.resp;
288
289
290 %% psd
291
292 load ..\m-toolbox\test\mpsd.mat % load mpsd.mat first column is f then psd1, csd and psd2
293
294 f = mpsd(:,1);
295 % psd = mpsd(:,2:4);
296 psd = mpsd(:,2);
297 y = abs(sqrt(psd));
298
299 %% auto search for the order
300
301 fs = 10; % Sampling frequency in Hz
302 % Fitting params
303 params = struct('spolesopt',2, 'Nmaxiter',70, 'minorder',24,...
304 'maxorder',35, 'weightparam',1, 'plot',0,...
305 'ctp','chivar','lrscond',5,'msevar',2,...
306 'stabfit',0,'dterm',0,'spy',1);
307 [res,poles,dterm,mresp,rdl,mse] = utils.math.autodfit(y,f,fs,params);
308
309
310 %% comparison
311 figure()
312 subplot(2,1,1);
313 loglog(f,abs(y),'k')
314 hold on
315 loglog(f,abs(mresp),'r')
316 xlabel('Frequency [Hz]')
317 ylabel('Amplitude')
318 legend('Original', 'CTFIT')
319
320 subplot(2,1,2);
321 semilogx(f,angle(y),'k')
322 hold on
323 semilogx(f,angle(mresp),'r')
324 xlabel('Frequency [Hz]')
325 ylabel('Amplitude')
326 legend('Original', 'CTFIT')
327
328 % plotting mean squared error
329 figure()
330 semilogy(mse,'-ok')
331 % hold on
332 % grid on
333 % semilogy(mlr2,'-or')
334
335 % plotting squared error variation
336 figure()
337 semilogy(abs(diff(mse)./mse(1:end-1)),'-ok')
338 % hold on
339 % grid on
340 % semilogy(abs(diff(mlr2)),'-or')
341