Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/linSubtract.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 % LINSUBTRACT subtracts a linear contribution from an input ao.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: LINSUBTRACT subtracts a linear contribution from an input ao.
5 % The methods assumes the input data to be synchronous. The
6 % user selects a filter to be applied to the data before
7 % fitting and a time segment where the fit is performed.
8 %
9 % CALL: c = linSubtract(a,b1,b2,b3,...,bN, pl)
10 %
11 % INPUTS: a - AO from where subtract linear contributions
12 % b - AOs with noise contributions
13 % pl - parameter list (see below)
14 %
15 % OUTPUTs: c - output AO with contributions subtracted (tsdata)
16 %
17 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'linSubtract')">Parameters Description</a>
18 %
19 % EXAMPLES:
20 %
21 % 1) Given the data (d):
22 %
23 % d = a + c1*b1 + c2*(b2+b3)^2
24 %
25 % where (bs) are noisy contributions added to a signal (a). To recover (a)
26 % in the [1 1e3] segment applying a [5e-2 0.1] 2nd order bandpass
27 % filter to the data, the call to the function would be
28 %
29 % pl = plist('type','bandpass',...
30 % 'fc',[5e-2 0.1],...
31 % 'order',2,...
32 % 'times',[1 1e3],...
33 % 'coupling',{{'n(1)'},{'(n(2) + n(3)).^2'}});
34 %
35 % a = linSubtract(d,b1,b2,b3, pl)
36 %
37 %
38 % VERSION: $Id: linSubtract.m,v 1.13 2011/04/08 08:56:15 hewitson Exp $
39 %
40 % TODO: option for parallel and serial subtraction
41 % handling errors
42 %
43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
44
45 function varargout = linSubtract(varargin)
46
47 % Check if this is a call for parameters
48 if utils.helper.isinfocall(varargin{:})
49 varargout{1} = getInfo(varargin{3});
50 return
51 end
52
53 import utils.const.*
54 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
55
56 % Method can not be used as a modifier
57 if nargout == 0
58 error('### tfe cannot be used as a modifier. Please give an output variable.');
59 end
60
61 % Collect input variable names
62 in_names = cell(size(varargin));
63 for ii = 1:nargin,in_names{ii} = inputname(ii);end
64
65 % Collect all AOs and plists
66 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
67 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
68
69 % Decide on a deep copy or a modify
70 bs = copy(as, nargout);
71
72 % Combine plists
73 pl = parse(pl, getDefaultPlist);
74
75
76 % get parameters
77 fc = find(pl,'fc');
78 if isempty(fc)
79 error('### Please define a cut-off frequency ''fc''');
80 end
81 times = find(pl,'times');
82 cp = find(pl,'coupling');
83 if isempty(cp)
84 error('### Please define a coupling model ''coupling''')
85 end
86 order = find(pl,'order');
87 type = find(pl,'type');
88
89 % split in time
90 if ~isempty(times)
91 cs = split(bs,plist('times',times));
92 else
93 cs = bs;
94 end
95
96 s = cs(1);
97 for i = 2:length(bs)
98 n(i-1) = cs(i);
99 end
100
101 subt = ao();
102 % Loop noise sources
103 for i=1:length(cp)
104 % coupling
105 nterm = ao();
106 if numel(cp{i}) == 1
107 nterm = eval([char(cp{i}) ';']);
108 else
109 nn = numel(cp{i});
110 for j =1:nn
111 nterm(j) = eval([char(cp{i}{j}) ';']);
112 end
113 end
114 % bandpass filter
115 fbp = miir(plist('type',type,'fc',fc,'order',order,'fs',s.fs));
116 sbp = filtfilt(s,fbp);
117 nterm_bp = filtfilt(nterm,fbp);
118 % linear fit
119 c = lscov(nterm_bp,sbp);
120 sn = lincom(nterm,c);
121 % subtract
122 s = s - sn;
123 end
124
125 % new tsdata
126 fsd = tsdata(s.x,s.y,s.fs);
127 % make output analysis object
128 cs = ao(fsd);
129 % set name
130 cs.name = sprintf('linSubtract(%s)', ao_invars{1});
131 % t0
132 if ~isempty(times)
133 cs.setT0(times(1));
134 else
135 cs.setT0(bs(1).t0);
136 end
137 % Propagate 'plotinfo'
138 plotinfo = [as(:).plotinfo];
139 if ~isempty(plotinfo)
140 cs.plotinfo = combine(plotinfo);
141 end
142 % Add history
143 cs.addHistory(getInfo('None'), pl, [ao_invars(:)], [bs(:).hist]);
144 % Set output
145 varargout{1} = cs;
146
147
148 end
149
150
151 %--------------------------------------------------------------------------
152 % Get Info Object
153 %--------------------------------------------------------------------------
154 function ii = getInfo(varargin)
155 if nargin == 1 && strcmpi(varargin{1}, 'None')
156 sets = {};
157 pl = [];
158 else
159 sets = {'Default'};
160 pl = getDefaultPlist;
161 end
162 % Build info object
163 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: linSubtract.m,v 1.13 2011/04/08 08:56:15 hewitson Exp $', sets, pl);
164 end
165
166
167 %--------------------------------------------------------------------------
168 % Get Default Plist
169 %--------------------------------------------------------------------------
170 function plout = getDefaultPlist()
171 persistent pl;
172 if exist('pl', 'var')==0 || isempty(pl)
173 pl = buildplist();
174 end
175 plout = pl;
176 end
177
178 function pl = buildplist()
179
180 pl = plist();
181
182 % Type
183 p = param({'type', 'Sets the type of filter used to fit the data.'}, {1, {'bandpass', 'bandreject', 'highpass', 'lowpass'}, paramValue.SINGLE});
184 pl.append(p);
185
186 % fc
187 p = param({'fc', 'Frequency cut-off of the filter.'}, paramValue.EMPTY_DOUBLE);
188 pl.append(p);
189
190 % Order
191 p = param({'order', 'Order of the filter.'}, {1, {2}, paramValue.OPTIONAL});
192 pl.append(p);
193
194 % Times
195 p = param({'times', 'A set of times where the fit+subtraction is applied.'}, paramValue.EMPTY_DOUBLE);
196 pl.append(p);
197
198 % Coupling
199 p = param({'coupling', ['A cell-array defining the model of the noise<br>'...
200 'terms to be subtracted. In the cell expression<br>'...
201 '''s'' stands for the input ao and ''n(i)'' for the N<br>' ...
202 'N noise contributions.']}, {1, {'{}'}, paramValue.OPTIONAL});
203 pl.append(p);
204
205 end
206 % PARAMETERS:
207 % 'type' - Sets the type of filter used to fit the data
208 % (help miir).
209 % 'fc' - Frequency cut-off of the filter (help miir)'
210 % 'order' - Order of the filter (help miir).
211 % 'times' - Sets split times where the subtraction applies
212 % (help split).
213 % 'coupling' - A cell-array defining the model of the noise
214 % terms to be subtracted. In the cell expression
215 % 's' stands for the input ao and 'n(i)' for the N
216 % N noise contributions.