Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/detrend.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 % DETREND detrends the input analysis object using a polynomial of degree N.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: DETREND detrends the input analysis object using a
5 % polynomial of degree N.
6 %
7 % CALL: b = detrend(a1,a2,a3,..., pl)
8 %
9 % INPUTS: aN - a list of analysis objects
10 % pl - a parameter list
11 %
12 % OUTPUTS: b - array of analysis objects
13 %
14 % If the last input argument is a parameter list (plist) it is used.
15 % The following parameters are recognised.
16 %
17 % NOTE: detrend uses two possible algorithms. A fast C-code implementation
18 % is used for orders less than 11. If the order is higher than 10, then a
19 % MATLAB code is used which is typically much slower. You can also force
20 % the use of the MATLAB code using a plist option. When interpreting the
21 % resulting coefficients, you must be clear which algorithm was used. For
22 % the C-code algorithm, the coefficients are scaled from the original by
23 %
24 % z = 2.*ii/(n-1)-1;
25 %
26 % such that the functional form (for a three-coefficient example) that is
27 % subtracted from the data can be recovered with:
28 %
29 % fitted_c = c.y(3).*z.^2 + c.y(2).*z + c.y(1);
30 %
31 %
32 % The procinfo field of the output AOs is filled with the following key/value
33 % pairs:
34 %
35 % 'COEFFS' - pest object with coefficients describing the subtracted trend function
36 % Note that they correspond to physical trend coefficents
37 % only if using the (slower) option 'M-FILE ONLY'. In the
38 % case of using the C-code algorith, see the note above.
39 %
40 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'detrend')">Parameters Description</a>
41 %
42 % VERSION: $Id: detrend.m,v 1.28 2011/04/08 08:56:11 hewitson Exp $
43 %
44 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45
46 function varargout = detrend(varargin)
47
48 % Check if this is a call for parameters
49 if utils.helper.isinfocall(varargin{:})
50 varargout{1} = getInfo(varargin{3});
51 return
52 end
53
54 import utils.const.*
55 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
56
57 % Collect input variable names
58 in_names = cell(size(varargin));
59 for ii = 1:nargin,in_names{ii} = inputname(ii);end
60
61 % Collect all AOs
62 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
63 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
64
65 % Decide on a deep copy or a modify
66 bs = copy(as, nargout);
67
68 % Combine plists
69 pl = parse(pl, getDefaultPlist);
70
71 % Unpack parameter list
72 % Leaving the "N" for backwards compatibility
73 order = find(pl, 'N', find(pl, 'order'));
74
75 % Loop over analysis objects
76 for jj = 1:numel(bs)
77 % check data
78 if ~isa(bs(jj).data, 'tsdata')
79 warning('!!! This method only works with time-series at the moment. Skipping %s', ao_invars{jj});
80 else
81 % detrend with polynomial
82 if find(pl, 'M-FILE ONLY')
83 if order >= 0
84 [y,c] = polydetrend(bs(jj).x, bs(jj).y, order);
85 else
86 y = bs(jj).y;
87 c = [];
88 end
89 else
90 try
91 [y,c] = ltpda_polyreg(bs(jj).y, order);
92 catch
93 warning('!!! failed to execture ltpda_polyreg.mex. Using m-file call.');
94 if order >= 0
95 [y,c] = polydetrend(bs(jj).x, bs(jj).y, order);
96 else
97 y = bs(jj).y;
98 c = [];
99 end
100 end
101 end
102 % set the values in the y and procinfo fields
103 bs(jj).data.setY(y);
104 % build a pest object with the coefficients
105 units = unit.initObjectWithSize(1, order + 1);
106 names = [];
107 mdl = [];
108
109 if order >= 0
110 for kk = 1:order+1
111 ll = order - kk + 1;
112 units(kk) = bs(jj).yunits / simplify((unit(bs(jj).xunits)).^(ll));
113 names{kk}= ['P' num2str(ll)];
114 mdl = [mdl '+P' num2str(ll) '*t.^' num2str(ll) ' '];
115 end
116
117 coeff = pest(plist(...
118 'names', names, ...
119 'y', c, ...
120 'dy', zeros(size(c)), ...
121 'yunits', units ...
122 ));
123 bs(jj).procinfo = plist('coeffs', coeff.setModels(smodel(mdl)));
124 end
125 % add name
126 bs(jj).name = sprintf('detrend(%s)', bs(jj).name);
127 % add history
128 bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), bs(jj).hist);
129 % Clear the errors since they don't make sense anymore
130 clearErrors(bs(jj));
131 end
132 end
133
134
135 % Set output
136 if nargout == numel(bs)
137 % List of outputs
138 for ii = 1:numel(bs)
139 varargout{ii} = bs(ii);
140 end
141 else
142 % Single output
143 varargout{1} = bs;
144 end
145 end
146
147 %--------------------------------------------------------------------------
148 % Get Info Object
149 %--------------------------------------------------------------------------
150 function ii = getInfo(varargin)
151
152 if nargin == 1 && strcmpi(varargin{1}, 'None')
153 sets = {};
154 pl = [];
155 else
156 sets = {'Default'};
157 pl = getDefaultPlist;
158 end
159 % Build info object
160 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);
161 end
162
163 %--------------------------------------------------------------------------
164 % Get Default Plist
165 %--------------------------------------------------------------------------
166
167 function plout = getDefaultPlist()
168 persistent pl;
169 if exist('pl', 'var')==0 || isempty(pl)
170 pl = buildplist();
171 end
172 plout = pl;
173 end
174
175 function pl = buildplist()
176 pl = plist();
177
178 % order
179 p = param({'order',['The order of detrending:<ul>', ...
180 '<li>-1 - no detrending</li>', ...
181 '<li>0 - subtract mean</li>', ...
182 '<li>1 - subtract linear fit</li>', ...
183 '<li>N - subtract fit of polynomial, order N</li></ul>']}, {3, {-1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9}, paramValue.SINGLE});
184 pl.append(p);
185
186 % m-file only
187 p = param({'M-FILE ONLY','Using SLOWER m-file call'}, paramValue.FALSE_TRUE);
188 pl.append(p);
189
190 end
191
192 function [y,p] = polydetrend(varargin)
193 % POLYDETREND detrends the input data vector with a polynomial.
194 %
195 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
196 %
197 % DESCRIPTION: POLYDETREND detrends the input data vector with a
198 % polynomial.
199 %
200 % CALL: [y,p] = polydetrend(x,order);
201 %
202 % INPUTS: x - vector of x values
203 % order - order of polynomial to fit and subtract
204 %
205 % OUTPUTS: y - detrended data
206 % p - coefficients fitted to the trend function
207 %
208 % VERSION: $Id: detrend.m,v 1.28 2011/04/08 08:56:11 hewitson Exp $
209 %
210 % HISTORY: 30-05-07 M Hewitson
211 % Creation
212 %
213 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
214
215 if nargin == 2
216 x = varargin{1};
217 order = varargin{2};
218 t = [1:length(x)].';
219 elseif nargin == 3
220 t = varargin{1};
221 x = varargin{2};
222 order = varargin{3};
223 else
224 error('### incorrect inputs.');
225 end
226
227 % fit polynomial
228 p = polyfit(t, x, order);
229
230 % make polynomial series
231 py = polyval(p, t);
232
233 % detrend
234 y = x - py;
235 end
236 % END
237