function varargout = ltpda_parfit_derivative(varargin)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DESCRIPTION: LTPDA_PARFIT_DERIVATIVE performs the numeric time derivative using the method% of 5 points parabolic fit [1]. The function can perform the second% derivative, the first derivative and the point estimation. % %%% CALL: Deriv = ltpda_parfit_derivative(a, ... , pl)% %% INPUTS: % - Aos containing the data to be differentiated.% - pl: is a parameter list containing the options for the% calculation. Parameters are:% - 'ORDER' whose accepeted options are:% - 'ZERO' perform the point estimation using the% coefficients vector d0 = [-3 12 17 12 -3]./35% - 'FIRST' perform the first derivative using the% couefficients vector d1 = [-2 -1 0 1 2]./(10*T)% - 'SECOND' perform the second derivative using the% coefficients vector d2 = [2 -1 -2 -1 2]./(7*T^2)% % NOTE1: T is the sampling period% NOTE2: The default option for 'ORDER' is 'SECOND'% %%% OUTPUTS: % - Deriv is the resulting Ao after the operation defined by the% input parameters%% REFERENCES:% [1] L. Carbone et al., Physical Review D 75, 042001 (2007)%% VERSION: $Id: ltpda_parfit_derivative.m,v 1.1 2008/04/24 16:13:12 luigi Exp $% % HISTORY: 18-03-2008 L Ferraioli% Creation%% The following call returns a parameter list object that contains the% default parameter values:%% >> pl = ltpda_parfit_derivative('Params')%% The following call returns a string that contains the routine CVS version:%% >> version = ltpda_parfit_derivative('Version')%% The following call returns a string that contains the routine category:%% >> category = ltpda_parfit_derivative('Category')%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Standard history variablesALGONAME = mfilename;VERSION = '$Id: ltpda_parfit_derivative.m,v 1.1 2008/04/24 16:13:12 luigi Exp $';CATEGORY = 'Signal Processing';%% Check if this is a call for parameters, the CVS version string% or the method categoryif nargin == 1 && ischar(varargin{1}) in = char(varargin{1}); if strcmpi(in, 'Params') varargout{1} = getDefaultPL(); return elseif strcmpi(in, 'Version') varargout{1} = VERSION; return elseif strcmpi(in, 'Category') varargout{1} = CATEGORY; return endend%% Collect input ao's, plist's and ao variable namesin_names = {};for ii = 1:nargin in_names{end+1} = inputname(ii); % Collect the inputs namesend[as, pl, invars] = collect_inputs(varargin, in_names);% produce one parameter listif isa(pl, 'plist') pl = combine(pl, getDefaultPL());else pl = getDefaultPL();end%% Initialize outputsbs = [];%% Go through analysis objects% Assigning coefficients values based on the input optionsswitch find(pl, 'ORDER') case 'ZERO' Coeffs = [-3 12 17 12 -3]./35; case 'FIRST' Coeffs = [-2 -1 0 1 2]./10; case 'SECOND' Coeffs = [2 -1 -2 -1 2]./7;endfor jj = 1:numel(as) a = as(jj); phi = a.data.y; fs = get(a, 'fs'); T = 1/fs; %% Building vectors for calculation phi_temp = [phi(1);phi(1);phi(1);phi(1);phi;phi(end);phi(end);phi(end);phi(end)]; % Switching between the input options switch find(pl, 'ORDER') case 'ZERO' Deriv = (Coeffs(1)*phi_temp(1:end-4) + Coeffs(2)*phi_temp(2:end-3) + Coeffs(3)*phi_temp(3:end-2) + Coeffs(4)*phi_temp(4:end-1) + Coeffs(5)*phi_temp(5:end)); Deriv = Deriv(3:end-2); case 'FIRST' Deriv = (1/T).*(Coeffs(1)*phi_temp(1:end-4) + Coeffs(2)*phi_temp(2:end-3) + Coeffs(3)*phi_temp(3:end-2) + Coeffs(4)*phi_temp(4:end-1) + Coeffs(5)*phi_temp(5:end)); Deriv = Deriv(3:end-2); case 'SECOND' Deriv = (1/T^2).*(Coeffs(1)*phi_temp(1:end-4) + Coeffs(2)*phi_temp(2:end-3) + Coeffs(3)*phi_temp(3:end-2) + Coeffs(4)*phi_temp(4:end-1) + Coeffs(5)*phi_temp(5:end)); Deriv = Deriv(3:end-2); end %% Building output Aos % Creating new time series data Deriv = tsdata(Deriv,fs); Deriv = set(Deriv, 'name', ... sprintf('ltpda_parfit_derivative(%s)', a.data.name), ... 'xunits', a.data.xunits, 'yunits', 'm/s^2'); % make new history objects h = history(ALGONAME, VERSION, pl, a.hist); h = set(h, 'invars', invars); % make output analysis objects b = ao(Deriv, h); clear Deriv % name, mfilename, description for these objects b = setnh(b, 'name', sprintf('ltpda_parfit_derivative(%s)', invars{jj}),... 'description', find(pl,'description')); %% Output the data bs = [bs b]; clear a bend%% Reshape the ouput to the same size of the inputvarargout{1} = reshape(bs, size(as));%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUNCTION: getDefaultPL%% DESCRIPTION: Get default params%% HISTORY: 01-02-2008 M Hueller% Creation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function plo = getDefaultPL()plo = plist('NAME', 'ltpda_parfit_derivative', 'ORDER', 'SECOND', ... 'DESCRIPTION','Parabolic Fit Derivative Method');% END