--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/lisovfit.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,337 @@
+% 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]