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