Home > classes > @ao > diff.m

diff

PURPOSE ^

DIFF differentiates the data in AO.

SYNOPSIS ^

function varargout = diff(varargin)

DESCRIPTION ^

 DIFF differentiates the data in AO.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: DIFF differentiates the data in AO. The result is a data
              series the same length as the input series.

 CALL:        >> b = diff(a,pl)

 INPUTS:      pl   - a parameter list
              a    - input analysis object

 OUTPUTS:
              b    - output analysis object containing the differentiated data.

 PARAMETERS:  method - the method to use: [default: '3POINT']
                       '2POINT' - 2 point derivative computed as
                                  [y(i+1)-y(i)]./[x(i+1)-x(i)].
                       '3POINT' - 3 point derivative. Compute derivative
                                  at i as [y(i+1)-y(i-1)] / [x(i+1)-x(i-1)].
                                  For i==1, the output is computed as
                                  [y(2)-y(1)]/[x(2)-x(1)]. The last sample
                                  is computed as [y(N)-y(N-1)]/[x(N)-x(N-1)].
                       '5POINT' - 5 point derivative. Compute derivative dx
                                  at i as
                                  [-y(i+2)+8*y(i+1)-8*y(i-1)+y(i-2)] /
                                       [3*(x(i+2)-x(i-2))].
                                  For i==1, the output is computed as
                                  [y(2)-y(1)]/[x(2)-x(1)]. The last sample
                                  is computed as [y(N)-y(N-1)]/[x(N)-x(N-1)].
                       'ORDER2' - Compute derivative using a 2nd order
                                  method.
                 'ORDER2SMOOTH' - Compute derivative using a 2nd order
                                  method with a parabolic fit to 5
                                  consecutive samples.
                       'filter' - applies an IIR filter built from a
                                  single pole at the chosen frequency. The
                                  filter is applied forwards and backwards
                                  (filtfilt) to achieve the desired f^2
                                  response. This only works for time-series
                                  AOs. For this method, you can specify the
                                  pole frequency with an additional parameter
                                  'f0' [default: 1/Nsecs]

 M-FILE INFO: Get information about this methods by calling
              >> ao.getInfo('diff')

              Get information about a specified set-plist by calling:
              >> ao.getInfo('diff', 'None')

 VERSION:     $Id: diff.m,v 1.10 2008/09/05 11:05:29 ingo Exp $

 HISTORY:     03-07-08 M Hewitson
                 Creation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % DIFF differentiates the data in AO.
0002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0003 %
0004 % DESCRIPTION: DIFF differentiates the data in AO. The result is a data
0005 %              series the same length as the input series.
0006 %
0007 % CALL:        >> b = diff(a,pl)
0008 %
0009 % INPUTS:      pl   - a parameter list
0010 %              a    - input analysis object
0011 %
0012 % OUTPUTS:
0013 %              b    - output analysis object containing the differentiated data.
0014 %
0015 % PARAMETERS:  method - the method to use: [default: '3POINT']
0016 %                       '2POINT' - 2 point derivative computed as
0017 %                                  [y(i+1)-y(i)]./[x(i+1)-x(i)].
0018 %                       '3POINT' - 3 point derivative. Compute derivative
0019 %                                  at i as [y(i+1)-y(i-1)] / [x(i+1)-x(i-1)].
0020 %                                  For i==1, the output is computed as
0021 %                                  [y(2)-y(1)]/[x(2)-x(1)]. The last sample
0022 %                                  is computed as [y(N)-y(N-1)]/[x(N)-x(N-1)].
0023 %                       '5POINT' - 5 point derivative. Compute derivative dx
0024 %                                  at i as
0025 %                                  [-y(i+2)+8*y(i+1)-8*y(i-1)+y(i-2)] /
0026 %                                       [3*(x(i+2)-x(i-2))].
0027 %                                  For i==1, the output is computed as
0028 %                                  [y(2)-y(1)]/[x(2)-x(1)]. The last sample
0029 %                                  is computed as [y(N)-y(N-1)]/[x(N)-x(N-1)].
0030 %                       'ORDER2' - Compute derivative using a 2nd order
0031 %                                  method.
0032 %                 'ORDER2SMOOTH' - Compute derivative using a 2nd order
0033 %                                  method with a parabolic fit to 5
0034 %                                  consecutive samples.
0035 %                       'filter' - applies an IIR filter built from a
0036 %                                  single pole at the chosen frequency. The
0037 %                                  filter is applied forwards and backwards
0038 %                                  (filtfilt) to achieve the desired f^2
0039 %                                  response. This only works for time-series
0040 %                                  AOs. For this method, you can specify the
0041 %                                  pole frequency with an additional parameter
0042 %                                  'f0' [default: 1/Nsecs]
0043 %
0044 % M-FILE INFO: Get information about this methods by calling
0045 %              >> ao.getInfo('diff')
0046 %
0047 %              Get information about a specified set-plist by calling:
0048 %              >> ao.getInfo('diff', 'None')
0049 %
0050 % VERSION:     $Id: diff.m,v 1.10 2008/09/05 11:05:29 ingo Exp $
0051 %
0052 % HISTORY:     03-07-08 M Hewitson
0053 %                 Creation
0054 %
0055 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0056 
0057 function varargout = diff(varargin)
0058 
0059   % Check if this is a call for parameters
0060   if utils.helper.isinfocall(varargin{:})
0061     varargout{1} = getInfo(varargin{3});
0062     return
0063   end
0064 
0065   import utils.const.*
0066   utils.helper.msg(msg.MNAME, 'running %s/%s', mfilename('class'), mfilename);
0067   
0068   % Collect input variable names
0069   in_names = cell(size(varargin));
0070   for ii = 1:nargin,in_names{ii} = inputname(ii);end
0071 
0072   % Collect all AOs and plists
0073   [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
0074   pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
0075 
0076   % Decide on a deep copy or a modify
0077   bs = copy(as, nargout);
0078 
0079   % combine plists
0080   pl = combine(pl, getDefaultPlist());
0081 
0082   % Extract method
0083   method = find(pl, 'method');
0084 
0085   for j=1:numel(bs)
0086 
0087     % Diff can't work for cdata objects since we need x data
0088     if isa(bs(j).data, 'cdata')
0089       error('### diff doesn''t work with cdata AOs since we need an x-data vector.');
0090     end
0091 
0092     % Compute derivative with selected method
0093     switch method
0094       case '2POINT'
0095         x      = bs(j).data.getX;
0096         dx     = diff(x);
0097         y      = bs(j).data.getY;
0098         dy     = diff(y);
0099         z      = dy./dx;
0100         bs(j).setY([z; 2*z(end)-z(end-1)], 'internal');
0101       case '3POINT'
0102         x  = bs(j).data.getX;
0103         dx = diff(x);
0104         y  = bs(j).data.getY;
0105         z  = zeros(size(y));
0106         z(2:end-1) = (y(3:end)-y(1:end-2)) ./ (dx(2:end)+dx(1:end-1));
0107         z(1)       = (y(2)-y(1)) ./ (dx(1));
0108         z(end)     = 2*z(end-1)-z(end-2);
0109         bs(j).setY(z, 'internal');
0110       case '5POINT'
0111         x  = bs(j).data.getX;
0112         dx = diff(x);
0113         y  = bs(j).data.getY;
0114         z  = zeros(size(y));
0115         z(1)       = (y(2)-y(1)) ./ (dx(1));
0116         z(2)       = (y(3)-y(1))./(dx(2)+dx(1));
0117         z(3:end-2) = (-y(5:end) + 8.*y(4:end-1) - 8.*y(2:end-3) + y(1:end-4)) ./ (3.*(x(5:end)-x(1:end-4)));
0118         z(end-1)   = 2*z(end-2)-z(end-3);
0119         z(end)     = 2*z(end-1)-z(end-2);
0120         bs(j).setY(z, 'internal');
0121       case 'ORDER2'
0122         x     = bs(j).data.getX;
0123         dx    = diff(x);
0124         y     = bs(j).data.getY;
0125         z     = zeros(size(y));
0126         m     = length(y);
0127         % y'(x1)
0128         z(1) = (1/dx(1)+1/dx(2))*(y(2)-y(1))+...
0129           dx(1)/(dx(1)*dx(2)+dx(2)^2)*(y(1)-y(3));
0130         % y'(xm)
0131         z(m) = (1/dx(m-2)+1/dx(m-1))*(y(m)-y(m-1))+...
0132           dx(m-1)/(dx(m-1)*dx(m-2)+dx(m-2)^2)*(y(m-2)-y(m));
0133         % y'(xi) (i>1 & i<m)
0134         dx1 = repmat(dx(1:m-2),1,1);
0135         dx2 = repmat(dx(2:m-1),1,1);
0136         y1 = y(1:m-2); y2 = y(2:m-1); y3 = y(3:m);
0137         z(2:m-1) = 1./(dx1.*dx2.*(dx1+dx2)).*...
0138           (-dx2.^2.*y1+(dx2.^2-dx1.^2).*y2+dx1.^2.*y3);
0139         bs(j).setY(z, 'internal');
0140       case 'ORDER2SMOOTH'
0141         x  = bs(j).data.getX;
0142         y  = bs(j).data.getY;
0143         dx = diff(x);
0144         m  = length(y);
0145         if max(abs(diff(dx)))>sqrt(eps(max(abs(dx))))
0146           error('### The x-step must be constant for method ''ORDER2SMOOTH''')
0147         elseif m<5
0148           error('### Length of y must be at least 5 for method ''ORDER2SMOOTH''.')
0149         end
0150         h = mean(dx);
0151         z = zeros(size(y));
0152         % y'(x1)
0153         z(1) = sum(y(1:5).*[-54; 13; 40; 27; -26])/70/h;
0154         % y'(x2)
0155         z(2) = sum(y(1:5).*[-34; 3; 20; 17; -6])/70/h;
0156         % y'(x{m-1})
0157         z(m-1) = sum(y(end-4:end).*[6; -17; -20; -3; 34])/70/h;
0158         % y'(xm)
0159         z(m) = sum(y(end-4:end).*[26; -27; -40; -13; 54])/70/h;
0160         % y'(xi) (i>2 & i<(N-1))
0161         Dc = [2 1 0 -1 -2];
0162         tmp = convn(Dc,y)/10/h;
0163         z(3:m-2) = tmp(5:m);
0164         bs(j).setY(z, 'internal');
0165       case 'filter'
0166       otherwise
0167         error('### Unknown method for computing the derivative.');
0168     end
0169 
0170     % add history
0171     bs(j).addHistory(getInfo, pl, ao_invars(j), bs(j).hist);
0172     % name for this object
0173     bs(j).setName(sprintf('diff(%s)', ao_invars{j}), 'internal');
0174   end
0175 
0176   % Set outputs
0177   if nargout > 0
0178     varargout{1} = bs;
0179   end
0180 end
0181 
0182 %--------------------------------------------------------------------------
0183 % Get Info Object
0184 %--------------------------------------------------------------------------
0185 function ii = getInfo(varargin)
0186 
0187   if nargin == 1 && strcmpi(varargin{1}, 'None')
0188     sets = {};
0189     pl   = [];
0190   else
0191     sets = {'Default'};
0192     pl   = getDefaultPlist;
0193   end
0194   % Build info object
0195   ii = minfo(mfilename, 'ao', '', utils.const.categories.sigproc, '$Id: diff.m,v 1.10 2008/09/05 11:05:29 ingo Exp $', sets, pl);
0196 end
0197 
0198 %--------------------------------------------------------------------------
0199 % Get Default Plist
0200 %--------------------------------------------------------------------------
0201 function pl_default = getDefaultPlist()
0202   pl_default = plist('method',  '3POINT');
0203 end
0204

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003