Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/fpsder.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % FPSDER performs the numeric time derivative
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: FPSDER (Five Points Stencil Derivative) performs the numeric
5 % time derivative using the method of five points stencil.
6 % The function can perform first, second, third and fourth derivetive
7 % of a series of input data.
8 % Given a discrete series of data points, the five-point-stencil method for
9 % the derivative approximation, at a given time t0, is calculated by means
10 % of finite differences between the element at t0 with its four neighbors.
11 % The n-order derivative at a certain time can be approximated by a five
12 % point difference equation:
13 %
14 % d^{n}y[k]
15 % --------- = (1/T^n) * {a*y[k-2] + b*y[k-1] + c*y[k] + d*y[k+1] + e*y[k+2]}
16 % dt^{n}
17 %
18 % It can be demonstrated [1,2] that the five coefficients [a, b, c, d, e] can
19 % be written in terms of only one of them. In fpsder the independent
20 % coefficient is fixed to be the first and is called m. It can be input as
21 % a parameter when the function is called.
22 %
23 %
24 % CALL: Deriv = fpsder(data, params)
25 %
26 % INPUTS:
27 %
28 % - a is a vector containing the data to be differentiated.
29 % - params is a struct with the input parameters:
30 %
31 % - 'ORDER' set the derivative order. Its allowed options are:
32 % - 'ZERO' perform data smoothing using the couefficients
33 % vector d0 = [m -4*m 1+6*m -4*m m].
34 % - 'FIRST' perform the first derivative using the
35 % couefficients vector d1 = [m -(0.5+2*m) 0 (0.5+2*m) m]./T.
36 % Recomended values of m are in the interval [-0.1, 0.1].
37 % - 'SECOND' perform the second derivative using the
38 % coefficients vector d2 = [m 1-4*m 6*m-2 1-4*m m]./(T^2).
39 % Recomended values of m are in the interval [-0.11, 0.3].
40 % - 'THIRD' perform the third derivative using the
41 % coefficients vector d3 = []./(T^3)
42 % - 'FOURTH' perform the third derivative using the
43 % coefficients vector d4 = []./(T^4)
44 %
45 % - 'COEFF' set m coefficient values.
46 % In case of data smoothing: m = -3/35 correspond to the
47 % parabolic fit approximation.
48 % In case of first order derivative: m = -1/5 correspond to the
49 % parabolic fit approximation, m = 1/12 correspond to the
50 % Taylor series approximation.
51 % In case of second order derivative: m = 2/7 corresponds to
52 % the parabolic fit approximation, m = -1/12 corresponds to the
53 % Taylor series approximation and m = 1/4 gives the notch
54 % feature at the Nyquist frequency
55 %
56 % - 'FS' set the data sampling frequency in Hz
57 %
58 % NOTE1: T is the sampling period
59 % NOTE2: The default option for 'ORDER' is 'SECOND'
60 % NOTE3: The default option for 'COEFF' is 2/7
61 % NOTE4: The default option for 'FS' is 10
62 %
63 % OUTPUTS:
64 % - D is a vector containing the resulting data after the
65 % differentiation procedure
66 %
67 % REFERENCES:
68 % [1] L. Ferraioli, M. Hueller and S. Vitale, Discrete derivative
69 % estimation in LISA Pathfinder data reduction,
70 % <a
71 % href="matlab:web('http://www.iop.org/EJ/abstract/0264-9381/26/9/094013/','-browser')">Class. Quantum Grav. 26 (2009) 094013.</a>
72 % [2] L. Ferraioli, M. Hueller and S. Vitale, Discrete derivative
73 % estimation in LISA Pathfinder data reduction,
74 % <a
75 % href="matlab:web('http://arxiv.org/abs/0903.0324v1','-browser')">http://arxiv.org/abs/0903.0324v1</a>
76 %
77 % EXAMPLES:
78 % - Performing the second order derivative of a series of data, m
79 % coefficient is fixed to 2/7 and data sampling frequency is
80 % fixed to 10 Hz.
81 % params = struct('ORDER', 'SECOND', 'COEFF', 2/7, 'FS', 10);
82 % Deriv = fpsder(data, params);
83 %
84 % VERSION: '$Id: fpsder.m,v 1.8 2010/04/09 09:59:59 mauro Exp $';
85 %
86 % HISTORY: 18-03-2008 L Ferraioli
87 % Creation
88 %
89 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
90
91
92 function Deriv = fpsder(a, params)
93
94 % Getting input parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
95 % Collect inputs
96
97 % Default input struct
98 defaultparams = struct('ORDER','SECOND',...
99 'COEFF',2/7,...
100 'FS',10);
101
102 names = {'ORDER','COEFF','FS'};
103
104 % collecting input and default params
105 if ~isempty(params)
106 for jj=1:length(names)
107 if isfield(params, names(jj)) && ~isempty(params.(names{1,jj}))
108 defaultparams.(names{1,jj}) = params.(names{1,jj});
109 end
110 end
111 end
112
113 % values for input variables
114 order = defaultparams.ORDER;
115 m = defaultparams.COEFF;
116 fs = defaultparams.FS;
117
118 % willing to work with columns
119 if size(a,2)>1
120 a = a.';
121 end
122
123 % Assigning coefficients values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
124
125 % Assigning coefficients values based on the input options
126 switch upper(order)
127 case 'ZERO'
128 Coeffs = [m -4*m 1+6*m -4*m m];
129 case 'FIRST'
130 Coeffs = [m -(0.5+2*m) 0 (0.5+2*m) -m];
131 case 'SECOND'
132 Coeffs = [m 1-4*m 6*m-2 1-4*m m];
133 case 'THIRD'
134 Coeffs = [0 0 0 0 0];
135 disp('Not yet implemented, sorry!');
136 case 'FOURTH'
137 Coeffs = [0 0 0 0 0];
138 disp('Not yet implemented, sorry!');
139 otherwise
140 error('### Unknown order %s', order);
141 end
142
143 % Sampling period
144 T = 1/fs;
145
146 % Building vectors for calculation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147
148 % Building the 'extended' vector for calculation
149 % a_temp = [a(1);a(1);a(1);a(1);a;a(end);a(end);a(end);a(end)];
150 a_temp = [2*a(1)-a((4+1):-1:2);a;2*a(end)-a((end-1):-1:end-4)];
151
152 % Switching between the input options differentiate
153 switch upper(order)
154 case 'ZERO'
155 Deriv = (Coeffs(1)*a_temp(1:end-4) + Coeffs(2)*a_temp(2:end-3) + Coeffs(3)*a_temp(3:end-2) + Coeffs(4)*a_temp(4:end-1) + Coeffs(5)*a_temp(5:end));
156 Deriv = Deriv(3:end-2);
157 case 'FIRST'
158 Deriv = (1/T).*(Coeffs(1)*a_temp(1:end-4) + Coeffs(2)*a_temp(2:end-3) + Coeffs(3)*a_temp(3:end-2) + Coeffs(4)*a_temp(4:end-1) + Coeffs(5)*a_temp(5:end));
159 Deriv = Deriv(3:end-2);
160 case 'SECOND'
161 Deriv = (1/T^2).*(Coeffs(1)*a_temp(1:end-4) + Coeffs(2)*a_temp(2:end-3) + Coeffs(3)*a_temp(3:end-2) + Coeffs(4)*a_temp(4:end-1) + Coeffs(5)*a_temp(5:end));
162 Deriv = Deriv(3:end-2);
163 case 'THIRD'
164 Deriv = (1/T^3).*(Coeffs(1)*a_temp(1:end-4) + Coeffs(2)*a_temp(2:end-3) + Coeffs(3)*a_temp(3:end-2) + Coeffs(4)*a_temp(4:end-1) + Coeffs(5)*a_temp(5:end));
165 Deriv = Deriv(3:end-2);
166 case 'FOURTH'
167 Deriv = (1/T^4).*(Coeffs(1)*a_temp(1:end-4) + Coeffs(2)*a_temp(2:end-3) + Coeffs(3)*a_temp(3:end-2) + Coeffs(4)*a_temp(4:end-1) + Coeffs(5)*a_temp(5:end));
168 Deriv = Deriv(3:end-2);
169 otherwise
170 error('### Unknown order %s', order);
171 end
172
173 end