--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/detrend.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,237 @@
+% DETREND detrends the input analysis object using a polynomial of degree N.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: DETREND detrends the input analysis object using a
+%              polynomial of degree N.
+%
+% CALL:        b = detrend(a1,a2,a3,..., pl)
+%
+% INPUTS:      aN  - a list of analysis objects
+%              pl  - a parameter list
+%
+% OUTPUTS:     b - array of analysis objects
+%
+%              If the last input argument is a parameter list (plist) it is used.
+%              The following parameters are recognised.
+%
+% NOTE: detrend uses two possible algorithms. A fast C-code implementation
+% is used for orders less than 11. If the order is higher than 10, then a
+% MATLAB code is used which is typically much slower. You can also force
+% the use of the MATLAB code using a plist option. When interpreting the
+% resulting coefficients, you must be clear which algorithm was used. For
+% the C-code algorithm, the coefficients are scaled from the original by
+% 
+%     z = 2.*ii/(n-1)-1;
+% 
+% such that the functional form (for a three-coefficient example) that is
+% subtracted from the data can be recovered with:
+% 
+% fitted_c = c.y(3).*z.^2 + c.y(2).*z + c.y(1); 
+% 
+% 
+% The procinfo field of the output AOs is filled with the following key/value
+% pairs:
+%
+%    'COEFFS' - pest object with coefficients describing the subtracted trend function
+%               Note that they correspond to physical trend coefficents
+%               only if using the (slower) option 'M-FILE ONLY'. In the
+%               case of using the C-code algorith, see the note above.
+%
+% <a href="matlab:utils.helper.displayMethodInfo('ao', 'detrend')">Parameters Description</a>
+%
+% VERSION:     $Id: detrend.m,v 1.28 2011/04/08 08:56:11 hewitson Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = detrend(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
+  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+  
+  % Decide on a deep copy or a modify
+  bs = copy(as, nargout);
+  
+  % Combine plists
+  pl = parse(pl, getDefaultPlist);
+  
+  % Unpack parameter list
+  % Leaving the "N" for backwards compatibility
+  order = find(pl, 'N', find(pl, 'order'));
+  
+  % Loop over analysis objects
+  for jj = 1:numel(bs)
+    % check data
+    if ~isa(bs(jj).data, 'tsdata')
+      warning('!!! This method only works with time-series at the moment. Skipping %s', ao_invars{jj});
+    else
+      % detrend with polynomial
+      if find(pl, 'M-FILE ONLY')
+        if order >= 0
+          [y,c] = polydetrend(bs(jj).x, bs(jj).y, order);
+        else
+          y = bs(jj).y;
+          c = [];
+        end
+      else
+        try
+          [y,c] = ltpda_polyreg(bs(jj).y, order);
+        catch
+          warning('!!! failed to execture ltpda_polyreg.mex. Using m-file call.');
+          if order >= 0
+            [y,c] = polydetrend(bs(jj).x, bs(jj).y, order);
+          else
+            y = bs(jj).y;
+            c = [];
+          end
+        end
+      end
+      % set the values in the y and procinfo fields
+      bs(jj).data.setY(y);
+      % build a pest object with the coefficients
+      units = unit.initObjectWithSize(1, order + 1);
+      names = [];
+      mdl   = [];
+      
+      if order >= 0
+        for kk = 1:order+1
+          ll = order - kk + 1;
+          units(kk) = bs(jj).yunits / simplify((unit(bs(jj).xunits)).^(ll));
+          names{kk}= ['P' num2str(ll)];
+          mdl = [mdl '+P' num2str(ll) '*t.^' num2str(ll) ' '];
+        end
+        
+        coeff = pest(plist(...
+          'names', names, ...
+          'y', c, ...
+          'dy', zeros(size(c)), ...
+          'yunits', units ...
+          ));
+        bs(jj).procinfo = plist('coeffs', coeff.setModels(smodel(mdl)));
+      end
+      % add name
+      bs(jj).name = sprintf('detrend(%s)', bs(jj).name);
+      % add history
+      bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), bs(jj).hist);
+      % Clear the errors since they don't make sense anymore
+      clearErrors(bs(jj));
+    end
+  end
+  
+  
+  % Set output
+  if nargout == numel(bs)
+    % List of outputs
+    for ii = 1:numel(bs)
+      varargout{ii} = bs(ii);
+    end
+  else
+    % Single output
+    varargout{1} = bs;
+  end
+end
+
+%--------------------------------------------------------------------------
+% Get Info Object
+%--------------------------------------------------------------------------
+function ii = getInfo(varargin)
+  
+  if nargin == 1 && strcmpi(varargin{1}, 'None')
+    sets = {};
+    pl   = [];
+  else
+    sets = {'Default'};
+    pl   = getDefaultPlist;
+  end
+  % Build info object
+  ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: detrend.m,v 1.28 2011/04/08 08:56:11 hewitson Exp $', sets, pl);
+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();
+  
+  % order
+  p = param({'order',['The order of detrending:<ul>', ...
+    '<li>-1 - no detrending</li>', ...
+    '<li>0 - subtract mean</li>', ...
+    '<li>1 - subtract linear fit</li>', ...
+    '<li>N - subtract fit of polynomial, order N</li></ul>']}, {3, {-1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9}, paramValue.SINGLE});
+  pl.append(p);
+  
+  % m-file only
+  p = param({'M-FILE ONLY','Using SLOWER m-file call'}, paramValue.FALSE_TRUE);
+  pl.append(p);
+  
+end
+
+function [y,p] = polydetrend(varargin)
+  % POLYDETREND detrends the input data vector with a polynomial.
+  %
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  %
+  % DESCRIPTION: POLYDETREND detrends the input data vector with a
+  %              polynomial.
+  %
+  % CALL:        [y,p] = polydetrend(x,order);
+  %
+  % INPUTS:      x     - vector of x values
+  %              order - order of polynomial to fit and subtract
+  %
+  % OUTPUTS:     y   - detrended data
+  %              p   - coefficients fitted to the trend function
+  %
+  % VERSION:     $Id: detrend.m,v 1.28 2011/04/08 08:56:11 hewitson Exp $
+  %
+  % HISTORY:     30-05-07 M Hewitson
+  %                 Creation
+  %
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  
+  if nargin == 2
+    x     = varargin{1};
+    order = varargin{2};
+    t     = [1:length(x)].';
+  elseif nargin == 3
+    t     = varargin{1};
+    x     = varargin{2};
+    order = varargin{3};
+  else
+    error('### incorrect inputs.');
+  end
+  
+  % fit polynomial
+  p = polyfit(t, x, order);
+  
+  % make polynomial series
+  py = polyval(p, t);
+  
+  % detrend
+  y = x - py;
+end
+% END
+