Mercurial > hg > ltpda

comparison m-toolbox/test/test_ao_xfit.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 % Tests for xfit
2 %
3 % $Id: test_ao_xfit.m,v 1.7 2011/05/12 07:58:57 mauro Exp $
4 %
5
6 %% Case 1: Fit with function in plist
7 % Fit to a frequency-series
8
9 % Create a frequency-series
10 datapl = plist('fsfcn', '0.01./(0.0001+f) + 5*abs(randn(size(f))) ', 'f1', 1e-5, 'f2', 5, 'nf', 1000, ...
11 'xunits', 'Hz', 'yunits', 'N/Hz');
12 data = ao(datapl);
13 data.setName;
14
15 % Do fit
16 fitpl = plist('Function', 'P(1)./(P(2) + Xdata) + P(3)', ...
17 'P0', [0.1 0.01 1]);
18 params = xfit(data, fitpl);
19
20 % Evaluate model
21 BestModel = eval(params, plist('type','fsdata','xdata',data,'xfield','x'));
22 BestModel.setName;
23
24 % Display results
25 iplot(data,BestModel)
26
27 %% Case 2: Fit with function in plist
28
29 % Create a noisy sine-wave
30 fs = 10;
31 nsecs = 500;
32 datapl = plist('waveform', 'Sine wave', 'f', 0.01, 'A', 0.6, 'fs', fs, 'nsecs', nsecs, ...
33 'xunits', 's', 'yunits', 'm');
34 sw = ao(datapl);
35 noise = ao(plist('tsfcn', '0.01*randn(size(t))', 'fs', fs, 'nsecs', nsecs));
36 data = sw+noise;
37 data.setName;
38
39 % Do fit
40 fitpl = plist('Function', 'P(1).*sin(2*pi*P(2).*Xdata + P(3))', ...
41 'P0', [1 0.01 0]);
42 params = xfit(data, fitpl);
43
44 % Evaluate model
45 BestModel = eval(params, plist('type','tsdata','xdata',data,'xfield','x'));
46 BestModel.setName;
47
48 % Display results
49 iplot(data,BestModel)
50
51 %% Case 3: Fit with smodel
52
53 % Fit an smodel of a straight line to some data
54
55 % Create a noisy straight-line
56 datapl = plist('xyfcn', '2.33 + 0.1*x + 0.01*randn(size(x))', 'x', 0:0.1:10, ...
57 'xunits', 's', 'yunits', 'm');
58 data = ao(datapl);
59 data.setName;
60
61 % Model to fit
62 mdl = smodel('a + b*x');
63 mdl.setXvar('x');
64 mdl.setParams({'a', 'b'}, {1 2});
65
66 % Fit model
67 fitpl = plist('Function', mdl, 'P0', [1 1]);
68 params = xfit(data, fitpl);
69
70 % Evaluate model
71 BestModel = eval(params,plist('xdata',data,'xfield','x'));
72 BestModel.setName;
73
74 % Display results
75 iplot(data,BestModel)
76
77 %% Case 4: Fit with smodel:
78 % Fit a chirp-sine firstly starting from an initial guess (quite close
79 % to the true values) (bad convergency) and secondly by a Monte Carlo
80 % search (good convergency)
81
82 % Create a noisy chirp-sine
83 fs = 10;
84 nsecs = 1000;
85
86 % Model to fit and generate signal
87 mdl = smodel(plist('name', 'chirp', 'expression', 'A.*sin(2*pi*(f + f0.*t).*t + p)', ...
88 'params', {'A','f','f0','p'}, 'xvar', 't', 'xunits', 's', 'yunits', 'm'));
89
90 % signal
91 s = mdl.setValues({10,1e-4,1e-5,0.3});
92 s.setXvals(0:1/fs:nsecs-1/fs);
93 signal = s.eval;
94 signal.setName;
95
96 % noise
97 noise = ao(plist('tsfcn', '1*randn(size(t))', 'fs', fs, 'nsecs', nsecs));
98
99 % data
100 data = signal + noise;
101 data.setName;
102
103 % Fit model from the starting guess
104 fitpl_ig = plist('Function', mdl, 'P0',[8,9e-5,9e-6,0]);
105 params_ig = xfit(data, fitpl_ig);
106
107 % Evaluate model
108 BestModel_ig = eval(params_ig,plist('xdata',data,'xfield','x'));
109 BestModel_ig.setName;
110
111 % Display results
112 iplot(data,BestModel_ig)
113
114 % Fit model by a Monte Carlo search
115 fitpl_mc = plist('Function', mdl, ...
116 'MonteCarlo', 'yes', 'Npoints', 1000, 'LB', [8,9e-5,9e-6,0], 'UB', [11,3e-4,2e-5,2*pi]);
117 params_mc = xfit(data, fitpl_mc);
118
119 % Evaluate model
120 BestModel_mc = eval(params_mc,plist('xdata',data,'xfield','x'));
121 BestModel_mc.setName;
122
123 % Display results
124 iplot(data,BestModel_mc)
125
126 %% Case 5: Fit multichannel with smodel
127
128 % Ch.1 data
129 datapl = plist('xyfcn', '0.1*x + 0.01*randn(size(x))', 'x', 0:0.1:10, 'name', 'channel 1', ...
130 'xunits', 'K', 'yunits', 'Pa');
131 a1 = ao(datapl);
132 % Ch.2 data
133 datapl = plist('xyfcn', '2.5*x + 0.1*sin(2*pi*x) + 0.01*randn(size(x))', 'x', 0:0.1:10, 'name', 'channel 2', ...
134 'xunits', 'K', 'yunits', 'T');
135 a2 = ao(datapl);
136
137 % Model to fit
138 mdl1 = smodel('a*x');
139 mdl1.setXvar('x');
140 mdl1.setParams({'a'}, {1});
141 mdl1.setXunits('K');
142 mdl1.setYunits('Pa');
143
144 mdl2 = smodel('b*x + a*sin(2*pi*x)');
145 mdl2.setXvar('x');
146 mdl2.setParams({'a','b'}, {1,2});
147 mdl2.setXunits('K');
148 mdl2.setYunits('T');
149
150 % Fit model
151 params = xfit(a1,a2, plist('Function', [mdl1,mdl2]));
152
153 % evaluate model
154 b = eval(params, plist('index',1,'xdata',a1,'xfield','x'));
155 b.setName('fit Ch.1');
156 r = a1-b;
157 r.setName('residuals');
158 iplot(a1,b,r)
159
160 b = eval(params, plist('index',2,'xdata',a2,'xfield','x'));
161 b.setName('fit Ch.2');
162 r = a2-b;
163 r.setName('residuals');
164 iplot(a2,b,r)
165
166
167