Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/filtSubtract.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 % FILTSUBTRACT subtracts a frequency dependent noise contribution from an input ao.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: FILTSUBTRACT subtracts a frequency dependent noise contribution from an input ao.
5 % The method computes the transfer function between both input AOs and
6 % fits a miir model to it. The frequency band is applied is set by a
7 % threshold in the coherence that the user defines as an input
8 % parameter.
9 %
10 % CALL: c = filtSubtract(a,b pl)
11 %
12 % INPUTS: a - AO from where subtract linear contributions
13 % b - AOs with noise contributions
14 % pl - parameter list (see below)
15 %
16 % OUTPUTs: c - output AO with contributions subtracted (tsdata)
17 %
18 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'filtSubtract')">Parameters Description</a>
19 %
20 %
21 % VERSION: $Id: filtSubtract.m,v 1.19 2011/05/22 22:26:56 mauro Exp $
22 %
23 % TODO: handling errors
24 % split by coherence function
25 %
26 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
27
28 function varargout = filtSubtract(varargin)
29
30 % Check if this is a call for parameters
31 if utils.helper.isinfocall(varargin{:})
32 varargout{1} = getInfo(varargin{3});
33 return
34 end
35
36 import utils.const.*
37 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
38
39 % Method can not be used as a modifier
40 if nargout == 0
41 error('### filtSubtract cannot be used as a modifier. Please give an output variable.');
42 end
43
44 % Collect input variable names
45 in_names = cell(size(varargin));
46 for ii = 1:nargin,in_names{ii} = inputname(ii);end
47
48 % Collect all AOs and plists
49 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
50 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
51
52 % Only two inputs ao's accepted
53 if numel(as) ~= 2
54 error('### filtSubtract only accepts two inputs AOs.');
55 end
56
57 % Decide on a deep copy or a modify
58 bs = copy(as, nargout);
59
60 % combine plists
61 pl = combine(pl, getDefaultPlist());
62
63 % get parameters
64 times = find(pl,'times');
65 filttimes = find(pl,'times postfilter');
66 fs = find(pl,'fs');
67 filt = find(pl,'filt');
68 fspl = find(pl,'frequencies');
69
70 % resample and consolidate
71 if isempty(fs)
72 bs1r = bs(1);
73 pl = pset(pl,'fs', bs(1).fs);
74 else
75 if fs ~= bs(1).fs
76 bs1r = resample(bs(1), plist('fsout',fs));
77 else
78 bs1r = bs(1);
79 end
80 end
81
82 % consolidate and split
83 c = consolidate(bs1r, bs(2), pl);
84 if ~isempty(times)
85 cs = split(c(1), plist('times', times));
86 cn = split(c(2), plist('times', times));
87 else
88 cs = c(1);
89 cn = c(2);
90 end
91
92 if isempty(filt)
93 % Transfer functions
94 tf = ltfe(cn, cs, pl);
95 % split transfer function to relevant frequencies
96 if ~isempty(fspl)
97 tf_spl = split(tf, plist('frequencies', [fspl(1) fspl(2)]));
98 else
99 tf_spl = tf;
100 end
101
102 % get filter from transfer function
103 fp = zDomainFit(tf_spl, pl);
104 fp.filters.setIunits(cn.yunits);
105 fp.filters.setOunits(cs.yunits);
106 else
107 fp = filt;
108 end
109
110 % get noise contribution
111 cs_cont = filter(cn,fp);
112 cs_cont.simplifyYunits();
113
114 % remove filter transient
115 if ~isempty(filttimes)
116 cs_cont = split(cs_cont, plist('times', filttimes));
117 cs = split(cs, plist('times', filttimes));
118 end
119
120 % subtraction
121 cs_subt = detrend(cs) - detrend(cs_cont);
122
123 % new tsdata
124 fsd = tsdata(cs_subt.x, cs_subt.y, cs_subt.fs);
125 % make output analysis object
126 cs = ao(fsd);
127 % set name
128 cs.name = sprintf('filtSubtract(%s)', ao_invars{1});
129 % set units
130 cs.setYunits(cs_subt.yunits);
131 % t0
132 if ~isempty(times) && ~isempty(filttimes)
133 cs.setT0(bs(1).t0 + times(1) + filttimes(1));
134 elseif isempty(times) && ~isempty(filttimes)
135 cs.setT0(bs(1).t0 + filttimes(1));
136 elseif ~isempty(times) && isempty(filttimes)
137 cs.setT0(bs(1).t0 + times(1));
138 else
139 cs.setT0(bs(1).t0);
140 end
141 % Add procinfo
142 cs.procinfo = plist('filter', fp);
143 % Add history
144 cs.addHistory(getInfo('None'), pl, ao_invars, [bs(:).hist]);
145 % Set output
146 varargout{1} = cs;
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: filtSubtract.m,v 1.19 2011/05/22 22:26:56 mauro 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') || isempty(pl)
173 pl = buildplist();
174 end
175 plout = pl;
176 end
177
178 function pl = buildplist()
179
180 pl = plist();
181
182 % fs
183 p = param({'fs','target sampling frequency to resample the data.'}, paramValue.EMPTY_DOUBLE);
184 pl.append(p);
185
186 % times
187 p = param({'times', 'selects the interval where the subtraction is applied.'},...
188 paramValue.EMPTY_STRING);
189 pl.append(p);
190
191 % times
192 p = param({'times postfilter', 'selects the filter transient intervals to be removed.'},...
193 paramValue.EMPTY_STRING);
194 pl.append(p);
195
196 % filtspl
197 p = param({'frequencies', 'selects the frequency band where the transfer<br>'...
198 'function is fitted.'}, paramValue.EMPTY_STRING);
199 pl.append(p);
200
201 % filt
202 p = param({'filt', ['a miir/mfir object which will be used as a<br>'...
203 'transfer function. If this option is selected<br>'...
204 ' the fit is avoided.']}, paramValue.EMPTY_STRING);
205 pl.append(p);
206
207 end
208