Mercurial > hg > ltpda
view m-toolbox/classes/@ao/lisovfit.m @ 16:91f21a0aab35 database-connection-manager
Update utils.jquery
* * *
Update utils.jmysql.getsinfo
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
% LISOVFIT uses LISO to fit a pole/zero model to the input frequency-series. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % DESCRIPTION: LISOVFIT uses LISO to fit a pole/zero model to the input % frequency-series. % % CALL: >> pzm = lisovfit(a,pl) % % INPUTS: pl - a parameter list % a - input analysis object % % OUTPUTS: % pzm - the fitted pzmodel. % % <a href="matlab:utils.helper.displayMethodInfo('ao', 'lisovfit')">Parameters Description</a> % % VERSION: $Id: lisovfit.m,v 1.14 2011/04/08 08:56:14 hewitson Exp $ % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function varargout = lisovfit(varargin) % Check if this is a call for parameters if utils.helper.isinfocall(varargin{:}) varargout{1} = getInfo(varargin{3}); return end import utils.const.* utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename); % Collect input variable names in_names = cell(size(varargin)); for ii = 1:nargin,in_names{ii} = inputname(ii);end % Collect all AOs and plists [bs, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names); [pl, pl_invars] = utils.helper.collect_objects(varargin(:), 'plist', in_names); % combine plists pl = parse(pl, getDefaultPlist()); % Extract parameters pzm0 = find(pl, 'PZM0'); pzml = find(pl, 'PZML'); pzmu = find(pl, 'PZMU'); delay = find(pl, 'delay'); f1 = find(pl, 'f1'); f2 = find(pl, 'f2'); nf = find(pl, 'nf'); method = find(pl, 'method'); np1 = find(pl, 'np1'); np2 = find(pl, 'np2'); % Check inputs if isempty(f1) || isempty(f2) || isempty(nf) error('### Specify the full frequency range with f1, f2, and nf'); end if f1 > f2 error('### The starting frequency should be less than the end frequency'); end if numel(delay) ~= 3 error('### The delay must be specified as a 3 element numerical vector: [limit start upper]'); end if ~(delay(1) < delay(2)) error('### The lower limit for the delay must be less than the starting guess.'); end if ~(delay(2) < delay(3)) error('### The upper limit for the delay must be greater than the starting guess.'); end if numel(pzm0.poles) ~= numel(pzml.poles) error('### The starting, lower, and upper models must have the same number of poles'); end if numel(pzm0.poles) ~= numel(pzmu.poles) error('### The starting, lower, and upper models must have the same number of poles'); end if numel(pzm0.zeros) ~= numel(pzml.zeros) error('### The starting, lower, and upper models must have the same number of zeros'); end if numel(pzm0.zeros) ~= numel(pzmu.zeros) error('### The starting, lower, and upper models must have the same number of zeros'); end % Loop over input AOs pzms = []; pzmls = []; pzmus = []; for j=1:numel(bs) if isa(bs(j).data, 'fsdata') % Generate temp file names outfile = [tempname '.fil']; datafile = [tempname '.dat']; % Write LISO fit file switch lower(method) case 'fit' if isreal(bs(j).data.getY) writeLISOfitFile(pzm0, pzml, pzmu, delay, f1, f2, nf, outfile, datafile, false); else writeLISOfitFile(pzm0, pzml, pzmu, delay, f1, f2, nf, outfile, datafile, true); end case 'vfit' writeLISOvfitFile(np1,np2,outfile, datafile); otherwise error('### method should be ''vfit'' or ''fit''.'); end % Export AO data file export(bs(j), datafile); % call fil utils.bin.fil(outfile); % get fitted model pzmfit = pzmodel(outfile); % pzmodel returns 3 models, set name for the first one pzmfit(1).name = sprintf('fit(%s)', pzm0.name); % Set units [ounits, iunits] = factor(bs(j).data.yunits); pzmfit(1).setIunits(iunits); pzmfit(1).setOunits(ounits); pzmfit(2).setIunits(iunits); pzmfit(2).setOunits(ounits); pzmfit(3).setIunits(iunits); pzmfit(3).setOunits(ounits); % add history pzmfit(1).addHistory(getInfo('None'), pl, ao_invars(j), bs(j).hist); pzmfit(2).addHistory(getInfo('None'), pl, ao_invars(j), bs(j).hist); pzmfit(3).addHistory(getInfo('None'), pl, ao_invars(j), bs(j).hist); % add to output pzms = [pzms pzmfit(1)]; pzmls = [pzmls pzmfit(2)]; pzmus = [pzmus pzmfit(3)]; else error('### unknown data type.'); end end % Set outputs if nargout == 1 varargout{1} = pzms; elseif nargout == 3 varargout{1} = pzms; varargout{2} = pzmls; varargout{3} = pzmus; else error('### Incorrect number of output arguments.'); end end %---------------------------------------------------------- % Write a liso vector fitting file % function writeLISOvfitFile(np1, np2, outfile, datafile) fd = fopen(outfile, 'w+'); fprintf(fd, '# Temporary LISO fitting file \n'); % Write fit command fprintf(fd, 'vfit %s reim rel %d %d \n', datafile,np1,np2); fprintf(fd, '\n'); % Close file fclose(fd); end %---------------------------------------------------------- % Write a liso fit file % function writeLISOfitFile(pzm0, pzml, pzmu, delay, f1, f2, nf, outfile, datafile, complexData) fd = fopen(outfile, 'w+'); fprintf(fd, '# Temporary LISO fitting file from pzmodel: %s\n\n', pzm0.name); % first output poles Np = numel(pzm0.poles); for k=1:Np pole = pzm0.poles(k); lpole = pzml.poles(k); upole = pzmu.poles(k); fprintf(fd, '# POLE %d\n', k); % write pole start if isnan(pole.q) if ~isnan(lpole.q) || ~isnan(upole.q) error('### Poles in start, lower, and upper models must be of the same for (real, complex)'); end fprintf(fd, 'pole %g\n', pole.f); else fprintf(fd, 'pole %g %g\n', pole.f, pole.q); end % write param range % param pole2:f 1.0746e-06 0.00010746 fprintf(fd, 'param pole%d:f %g %g\n', k-1, lpole.f, upole.f); if ~isnan(pole.q) fprintf(fd, 'param pole%d:q %g %g\n', k-1, lpole.q, upole.q); end fprintf(fd, '\n'); end % then output zeros Nz = numel(pzm0.zeros); for k=1:Nz zero = pzm0.zeros(k); lzero = pzml.zeros(k); uzero = pzmu.zeros(k); fprintf(fd, '# ZERO %d\n', k); % write pole start if isnan(zero.q) if ~isnan(lzero.q) || ~isnan(uzero.q) error('### Zeros in start, lower, and upper models must be of the same for (real, complex)'); end fprintf(fd, 'zero %g\n', zero.f); else fprintf(fd, 'zero %g %g\n', zero.f, zero.q); end % write param range % param pole2:f 1.0746e-06 0.00010746 fprintf(fd, 'param zero%d:f %g %g\n', k-1, lzero.f, uzero.f); if ~isnan(zero.q) fprintf(fd, 'param zero%d:q %g %g\n', k-1, lzero.q, uzero.q); end fprintf(fd, '\n'); end % Write delay out if ~isempty(delay) fprintf(fd, 'delay %g\n', delay(2)); fprintf(fd, 'param delay %g %g\n', delay(1), delay(3)); end fprintf(fd, '\n'); % Write factor fprintf(fd, 'factor %g\n', pzm0.gain); fprintf(fd, 'param factor %g %g\n', pzml.gain, pzmu.gain); fprintf(fd, '\n'); % Write frequency command fprintf(fd, 'freq lin %g %g %g\n', f1, f2, nf); fprintf(fd, '\n'); % Write fit command if complexData fprintf(fd, 'fit %s reim semi\n', datafile); else fprintf(fd, 'fit %s abs rel\n', datafile); end fprintf(fd, 'rewrite samebetter\n'); fprintf(fd, '\n'); % Close file fclose(fd); end %-------------------------------------------------------------------------- % Get Info Object %-------------------------------------------------------------------------- function ii = getInfo(varargin) if nargin == 1 && strcmpi(varargin{1}, 'None') sets = {}; pls = []; else sets = {'Default'}; pls = getDefaultPlist; end % Build info object 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); ii.setModifier(false); end %-------------------------------------------------------------------------- % Get Default Plist %-------------------------------------------------------------------------- function plout = getDefaultPlist() persistent pl; if exist('pl', 'var')==0 || isempty(pl) pl = buildplist(); end plout = pl; end function pl = buildplist() pl = plist(); % PZM0 p = param({'PZM0', 'A pzmodel describing the starting guess.'}, {1, {pzmodel}, paramValue.OPTIONAL}); pl.append(p); % PZML p = param({'PZMU', 'A pzmodel describing the upper-bound.'}, {1, {pzmodel}, paramValue.OPTIONAL}); pl.append(p); % PZMU p = param({'PZML', 'A pzmodel describing the lower-bound.'}, {1, {pzmodel}, paramValue.OPTIONAL}); pl.append(p); % Delay p = param({'Delay', 'A 3-element numeric vector describing<br>a time-delay to include in the fit.'}, {1, {[0 1 10]}, paramValue.OPTIONAL}); pl.append(p); % F1 p = param({'F1', 'A start freqeuency to fit over'}, {1, {0}, paramValue.OPTIONAL}); pl.append(p); % F2 p = param({'F2', 'The end freqeuency to fit over'}, {1, {1}, paramValue.OPTIONAL}); pl.append(p); % Nf p = param({'NF','The number of frequencies to include in the fit.'}, {1, {100}, paramValue.OPTIONAL}); pl.append(p); % method p = param({'method', 'The fitting method.'}, {1, {'fit', 'vfit'}, paramValue.SINGLE}); pl.append(p); % NP1 p = param({'NP1', 'The minimum number of poles for vfit.'}, {1, {2}, paramValue.OPTIONAL}); pl.append(p); % NP2 p = param({'NP2', 'The maximum number of poles for vfit.'}, {1, {10}, paramValue.OPTIONAL}); pl.append(p); end % PARAMETERS: % 'PZM0' - a pzmodel describing the starting guess. % 'PZML' - a pzmodel describing the lower-bound. % 'PZMU' - a pzmodel describing the upper-bound. % 'DELAY' - a 3-element numeric vector describing a time-delay % to include in the fit: [lower start upper] % 'f1' - start freqeuency to fit over [default: taken from input AO] % 'f2' - end freqeuency to fit over [default: taken from input AO] % 'nf' - number of frequencies to fit [default: taken from input AO] % 'method'- 'fit' for fitting or 'vfit' for vector fitting [default: fit] % 'np1' - lower num. poles for vfit [default: 2] % 'np2' - upper num. poles for vfit [default: 10]