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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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