Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/lisovfit.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 % LISOVFIT uses LISO to fit a pole/zero model to the input frequency-series.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: LISOVFIT uses LISO to fit a pole/zero model to the input
5 % frequency-series.
6 %
7 % CALL: >> pzm = lisovfit(a,pl)
8 %
9 % INPUTS: pl - a parameter list
10 % a - input analysis object
11 %
12 % OUTPUTS:
13 % pzm - the fitted pzmodel.
14 %
15 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'lisovfit')">Parameters Description</a>
16 %
17 % VERSION: $Id: lisovfit.m,v 1.14 2011/04/08 08:56:14 hewitson Exp $
18 %
19 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
20
21 function varargout = lisovfit(varargin)
22
23 % Check if this is a call for parameters
24 if utils.helper.isinfocall(varargin{:})
25 varargout{1} = getInfo(varargin{3});
26 return
27 end
28
29 import utils.const.*
30 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
31
32 % Collect input variable names
33 in_names = cell(size(varargin));
34 for ii = 1:nargin,in_names{ii} = inputname(ii);end
35
36 % Collect all AOs and plists
37 [bs, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
38 [pl, pl_invars] = utils.helper.collect_objects(varargin(:), 'plist', in_names);
39
40 % combine plists
41 pl = parse(pl, getDefaultPlist());
42
43 % Extract parameters
44 pzm0 = find(pl, 'PZM0');
45 pzml = find(pl, 'PZML');
46 pzmu = find(pl, 'PZMU');
47 delay = find(pl, 'delay');
48 f1 = find(pl, 'f1');
49 f2 = find(pl, 'f2');
50 nf = find(pl, 'nf');
51 method = find(pl, 'method');
52 np1 = find(pl, 'np1');
53 np2 = find(pl, 'np2');
54
55
56 % Check inputs
57 if isempty(f1) || isempty(f2) || isempty(nf)
58 error('### Specify the full frequency range with f1, f2, and nf');
59 end
60 if f1 > f2
61 error('### The starting frequency should be less than the end frequency');
62 end
63 if numel(delay) ~= 3
64 error('### The delay must be specified as a 3 element numerical vector: [limit start upper]');
65 end
66 if ~(delay(1) < delay(2))
67 error('### The lower limit for the delay must be less than the starting guess.');
68 end
69 if ~(delay(2) < delay(3))
70 error('### The upper limit for the delay must be greater than the starting guess.');
71 end
72 if numel(pzm0.poles) ~= numel(pzml.poles)
73 error('### The starting, lower, and upper models must have the same number of poles');
74 end
75 if numel(pzm0.poles) ~= numel(pzmu.poles)
76 error('### The starting, lower, and upper models must have the same number of poles');
77 end
78 if numel(pzm0.zeros) ~= numel(pzml.zeros)
79 error('### The starting, lower, and upper models must have the same number of zeros');
80 end
81 if numel(pzm0.zeros) ~= numel(pzmu.zeros)
82 error('### The starting, lower, and upper models must have the same number of zeros');
83 end
84
85 % Loop over input AOs
86 pzms = [];
87 pzmls = [];
88 pzmus = [];
89 for j=1:numel(bs)
90 if isa(bs(j).data, 'fsdata')
91 % Generate temp file names
92 outfile = [tempname '.fil'];
93 datafile = [tempname '.dat'];
94 % Write LISO fit file
95 switch lower(method)
96 case 'fit'
97 if isreal(bs(j).data.getY)
98 writeLISOfitFile(pzm0, pzml, pzmu, delay, f1, f2, nf, outfile, datafile, false);
99 else
100 writeLISOfitFile(pzm0, pzml, pzmu, delay, f1, f2, nf, outfile, datafile, true);
101 end
102 case 'vfit'
103 writeLISOvfitFile(np1,np2,outfile, datafile);
104 otherwise
105 error('### method should be ''vfit'' or ''fit''.');
106 end
107 % Export AO data file
108 export(bs(j), datafile);
109 % call fil
110 utils.bin.fil(outfile);
111 % get fitted model
112 pzmfit = pzmodel(outfile);
113 % pzmodel returns 3 models, set name for the first one
114 pzmfit(1).name = sprintf('fit(%s)', pzm0.name);
115 % Set units
116 [ounits, iunits] = factor(bs(j).data.yunits);
117 pzmfit(1).setIunits(iunits);
118 pzmfit(1).setOunits(ounits);
119 pzmfit(2).setIunits(iunits);
120 pzmfit(2).setOunits(ounits);
121 pzmfit(3).setIunits(iunits);
122 pzmfit(3).setOunits(ounits);
123 % add history
124 pzmfit(1).addHistory(getInfo('None'), pl, ao_invars(j), bs(j).hist);
125 pzmfit(2).addHistory(getInfo('None'), pl, ao_invars(j), bs(j).hist);
126 pzmfit(3).addHistory(getInfo('None'), pl, ao_invars(j), bs(j).hist);
127 % add to output
128 pzms = [pzms pzmfit(1)];
129 pzmls = [pzmls pzmfit(2)];
130 pzmus = [pzmus pzmfit(3)];
131 else
132 error('### unknown data type.');
133 end
134 end
135
136 % Set outputs
137 if nargout == 1
138 varargout{1} = pzms;
139 elseif nargout == 3
140 varargout{1} = pzms;
141 varargout{2} = pzmls;
142 varargout{3} = pzmus;
143 else
144 error('### Incorrect number of output arguments.');
145 end
146 end
147
148 %----------------------------------------------------------
149 % Write a liso vector fitting file
150 %
151 function writeLISOvfitFile(np1, np2, outfile, datafile)
152 fd = fopen(outfile, 'w+');
153
154 fprintf(fd, '# Temporary LISO fitting file \n');
155
156 % Write fit command
157 fprintf(fd, 'vfit %s reim rel %d %d \n', datafile,np1,np2);
158 fprintf(fd, '\n');
159
160 % Close file
161 fclose(fd);
162 end
163
164 %----------------------------------------------------------
165 % Write a liso fit file
166 %
167 function writeLISOfitFile(pzm0, pzml, pzmu, delay, f1, f2, nf, outfile, datafile, complexData)
168 fd = fopen(outfile, 'w+');
169
170 fprintf(fd, '# Temporary LISO fitting file from pzmodel: %s\n\n', pzm0.name);
171
172 % first output poles
173 Np = numel(pzm0.poles);
174 for k=1:Np
175 pole = pzm0.poles(k);
176 lpole = pzml.poles(k);
177 upole = pzmu.poles(k);
178 fprintf(fd, '# POLE %d\n', k);
179 % write pole start
180 if isnan(pole.q)
181 if ~isnan(lpole.q) || ~isnan(upole.q)
182 error('### Poles in start, lower, and upper models must be of the same for (real, complex)');
183 end
184 fprintf(fd, 'pole %g\n', pole.f);
185 else
186 fprintf(fd, 'pole %g %g\n', pole.f, pole.q);
187 end
188 % write param range
189 % param pole2:f 1.0746e-06 0.00010746
190 fprintf(fd, 'param pole%d:f %g %g\n', k-1, lpole.f, upole.f);
191 if ~isnan(pole.q)
192 fprintf(fd, 'param pole%d:q %g %g\n', k-1, lpole.q, upole.q);
193 end
194 fprintf(fd, '\n');
195 end
196
197 % then output zeros
198 Nz = numel(pzm0.zeros);
199 for k=1:Nz
200 zero = pzm0.zeros(k);
201 lzero = pzml.zeros(k);
202 uzero = pzmu.zeros(k);
203 fprintf(fd, '# ZERO %d\n', k);
204 % write pole start
205 if isnan(zero.q)
206 if ~isnan(lzero.q) || ~isnan(uzero.q)
207 error('### Zeros in start, lower, and upper models must be of the same for (real, complex)');
208 end
209 fprintf(fd, 'zero %g\n', zero.f);
210 else
211 fprintf(fd, 'zero %g %g\n', zero.f, zero.q);
212 end
213 % write param range
214 % param pole2:f 1.0746e-06 0.00010746
215 fprintf(fd, 'param zero%d:f %g %g\n', k-1, lzero.f, uzero.f);
216 if ~isnan(zero.q)
217 fprintf(fd, 'param zero%d:q %g %g\n', k-1, lzero.q, uzero.q);
218 end
219 fprintf(fd, '\n');
220 end
221
222 % Write delay out
223 if ~isempty(delay)
224 fprintf(fd, 'delay %g\n', delay(2));
225 fprintf(fd, 'param delay %g %g\n', delay(1), delay(3));
226 end
227 fprintf(fd, '\n');
228
229 % Write factor
230 fprintf(fd, 'factor %g\n', pzm0.gain);
231 fprintf(fd, 'param factor %g %g\n', pzml.gain, pzmu.gain);
232 fprintf(fd, '\n');
233
234 % Write frequency command
235 fprintf(fd, 'freq lin %g %g %g\n', f1, f2, nf);
236 fprintf(fd, '\n');
237
238 % Write fit command
239 if complexData
240 fprintf(fd, 'fit %s reim semi\n', datafile);
241 else
242 fprintf(fd, 'fit %s abs rel\n', datafile);
243 end
244 fprintf(fd, 'rewrite samebetter\n');
245 fprintf(fd, '\n');
246
247 % Close file
248 fclose(fd);
249 end
250
251 %--------------------------------------------------------------------------
252 % Get Info Object
253 %--------------------------------------------------------------------------
254 function ii = getInfo(varargin)
255 if nargin == 1 && strcmpi(varargin{1}, 'None')
256 sets = {};
257 pls = [];
258 else
259 sets = {'Default'};
260 pls = getDefaultPlist;
261 end
262 % Build info object
263 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: lisovfit.m,v 1.14 2011/04/08 08:56:14 hewitson Exp $', sets, pls);
264 ii.setModifier(false);
265 end
266
267 %--------------------------------------------------------------------------
268 % Get Default Plist
269 %--------------------------------------------------------------------------
270 function plout = getDefaultPlist()
271 persistent pl;
272 if exist('pl', 'var')==0 || isempty(pl)
273 pl = buildplist();
274 end
275 plout = pl;
276 end
277
278 function pl = buildplist()
279
280 pl = plist();
281
282 % PZM0
283 p = param({'PZM0', 'A pzmodel describing the starting guess.'}, {1, {pzmodel}, paramValue.OPTIONAL});
284 pl.append(p);
285
286 % PZML
287 p = param({'PZMU', 'A pzmodel describing the upper-bound.'}, {1, {pzmodel}, paramValue.OPTIONAL});
288 pl.append(p);
289
290 % PZMU
291 p = param({'PZML', 'A pzmodel describing the lower-bound.'}, {1, {pzmodel}, paramValue.OPTIONAL});
292 pl.append(p);
293
294 % Delay
295 p = param({'Delay', 'A 3-element numeric vector describing<br>a time-delay to include in the fit.'}, {1, {[0 1 10]}, paramValue.OPTIONAL});
296 pl.append(p);
297
298 % F1
299 p = param({'F1', 'A start freqeuency to fit over'}, {1, {0}, paramValue.OPTIONAL});
300 pl.append(p);
301
302 % F2
303 p = param({'F2', 'The end freqeuency to fit over'}, {1, {1}, paramValue.OPTIONAL});
304 pl.append(p);
305
306 % Nf
307 p = param({'NF','The number of frequencies to include in the fit.'}, {1, {100}, paramValue.OPTIONAL});
308 pl.append(p);
309
310 % method
311 p = param({'method', 'The fitting method.'}, {1, {'fit', 'vfit'}, paramValue.SINGLE});
312 pl.append(p);
313
314 % NP1
315 p = param({'NP1', 'The minimum number of poles for vfit.'}, {1, {2}, paramValue.OPTIONAL});
316 pl.append(p);
317
318 % NP2
319 p = param({'NP2', 'The maximum number of poles for vfit.'}, {1, {10}, paramValue.OPTIONAL});
320 pl.append(p);
321
322
323
324 end
325
326 % PARAMETERS:
327 % 'PZM0' - a pzmodel describing the starting guess.
328 % 'PZML' - a pzmodel describing the lower-bound.
329 % 'PZMU' - a pzmodel describing the upper-bound.
330 % 'DELAY' - a 3-element numeric vector describing a time-delay
331 % to include in the fit: [lower start upper]
332 % 'f1' - start freqeuency to fit over [default: taken from input AO]
333 % 'f2' - end freqeuency to fit over [default: taken from input AO]
334 % 'nf' - number of frequencies to fit [default: taken from input AO]
335 % 'method'- 'fit' for fitting or 'vfit' for vector fitting [default: fit]
336 % 'np1' - lower num. poles for vfit [default: 2]
337 % 'np2' - upper num. poles for vfit [default: 10]