Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/whiten1D.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 % WHITEN1D whitens the input time-series.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: WHITEN1D whitens the input time-series. The filter is built
5 % by fitting to the model provided. If no model is provided, a
6 % fit is made to a spectral-density estimate of the
7 % time-series (made using psd+bin_data or lpsd).
8 % Note: The function assumes that the input model corresponds
9 % to the one-sided psd of the data to be whitened.
10 %
11 % ALGORITHM:
12 % 1) If no model provided, make psd+bin_data or lpsd
13 % of time-series and take it as a model
14 % for the data power spectral density
15 % 2) Fit a set of partial fraction z-domain filters using
16 % utils.math.psd2wf. The fit is automatically stopped when
17 % the accuracy tolerance is reached.
18 % 3) Convert to array of MIIR filters
19 % 4) Assemble into a parallel filterbank object
20 % 5) Filter time-series in parallel
21 %
22 %
23 % CALL: b = whiten1D(a, pl)
24 % [b1,b2,...,bn] = whiten1D(a1,a2,...,an, pl);
25 %
26 % INPUT:
27 % - as are time-series analysis objects or a vector of
28 % analysis objects
29 % - pl is a plist with the input parameters
30 %
31 % OUTPUT:
32 % - bs "whitened" time-series AOs. The whitening filters used
33 % are stored in the objects procinfo field under the
34 % parameter 'Filter'.
35 %
36 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'whiten1D')">Parameters Description</a>
37 %
38 % VERSION: $Id: whiten1D.m,v 1.43 2011/04/08 08:56:12 hewitson Exp $
39 %
40 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
41
42 function varargout = whiten1D(varargin)
43
44 % Check if this is a call for parameters
45 if utils.helper.isinfocall(varargin{:})
46 varargout{1} = getInfo(varargin{3});
47 return
48 end
49
50 import utils.const.*
51 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
52
53 % Collect input variable names
54 in_names = cell(size(varargin));
55 for ii = 1:nargin,in_names{ii} = inputname(ii);end
56
57 % Collect all AOs and plists
58 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
59 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
60
61 % Decide on a deep copy or a modify
62 bs = copy(as, nargout);
63 inhists = [as.hist];
64
65 % combine plists
66 if isempty(pl)
67 model = 'psd';
68 else
69 model = find(pl, 'model');
70 if isempty(model)
71 model = 'psd';
72 end
73 end
74
75 if ischar(model)
76 pl = parse(pl, getDefaultPlist(model));
77 else
78 pl = parse(pl, getDefaultPlist('Default'));
79 end
80 pl.getSetRandState();
81
82 scale = find(pl, 'scaleOut');
83 flim = find(pl, 'flim');
84
85
86 % Loop over input AOs
87 for jj = 1:numel(as)
88 if ~isa(as(jj).data, 'tsdata')
89 utils.helper.msg(msg.IMPORTANT, '%s expects ao/tsdata objects. Skipping AO %s', mfilename, ao_invars{jj});
90 else
91
92 %-------------- Whiten this AO
93
94 % 1) Build whitening filterbank
95 switch class(model)
96 case 'char'
97 % Model is to be evaluated from data
98 in = as(jj);
99 pl.pset('model', model);
100 case 'ao'
101 % Model was provided as fsdata
102 in = model;
103 pl.pset('model', []);
104 end
105 wf = buildWhitener1D(in, pl);
106
107 % 1.5) Scale the date if demanded
108 if (scale)
109 spsd = lpsd(as(jj));
110 freqs = spsd.x;
111 if isempty(flim)
112 error('Please specify a flim field, to know the analysis band.');
113 elseif (flim(2) < flim(1))
114 error('flim should go from the smaller frequency to the bigger frequency. Please reverse them!')
115 else
116 index = find((freqs > flim(1)) & (freqs < flim(2)));
117 end
118
119 v1 = spsd.y(index(1):index(end-1));
120 v2 = spsd.y(index(2):index(end));
121 m = (v1 + v2) /2;
122 p = sum(m.* diff(freqs(index(1):index(end))));
123 end
124
125 % 2) Filter data and scale it if necessary
126 bs(jj).filter(wf);
127 if (scale)
128 bs(jj) = bs(jj) * sqrt(p);
129 end
130
131
132 % 3) Output data
133 % name for this object
134 bs(jj).name = sprintf('whiten1D(%s)', ao_invars{jj});
135 % Collect the filters into procinfo
136 bs(jj).procinfo = combine(plist('Filter', wf.filters), as(jj).procinfo);
137 if(scale)
138 bs(jj).procinfo = combine(plist('ScaleFactor', p, 'Filter', wf.filters), as(jj).procinfo);
139 end
140 % Make sure that the output yunits are empty
141 if ~eq(bs(jj).yunits, unit(''))
142 utils.helper.msg(msg.PROC1, 'Resetting output yunits to empty');
143 bs(jj).setYunits(unit(''));
144 end
145 % add history
146 bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), inhists(jj));
147 % clear errors
148 bs(jj).clearErrors;
149
150
151 end
152 end
153
154
155
156 % Set output
157 if nargout == numel(bs)
158 % List of outputs
159 for ii = 1:numel(bs)
160 varargout{ii} = bs(ii);
161 end
162 else
163 % Single output
164 varargout{1} = bs;
165 end
166 end
167
168 %--------------------------------------------------------------------------
169 % Get Info Object
170 %--------------------------------------------------------------------------
171 function ii = getInfo(varargin)
172 if nargin == 1 && strcmpi(varargin{1}, 'None')
173 sets = {};
174 pl = [];
175 elseif nargin == 1 && ~isempty(varargin{1}) && ischar(varargin{1})
176 sets{1} = varargin{1};
177 pl = getDefaultPlist(sets{1});
178 else
179 sets = SETS();
180 % get plists
181 pl(size(sets)) = plist;
182 for kk = 1:numel(sets)
183 pl(kk) = getDefaultPlist(sets{kk});
184 end
185 end
186 % Build info object
187 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: whiten1D.m,v 1.43 2011/04/08 08:56:12 hewitson Exp $', sets, pl);
188 end
189
190
191 %--------------------------------------------------------------------------
192 % Defintion of Sets
193 %--------------------------------------------------------------------------
194
195 function out = SETS()
196 out = ao.getInfo('buildWhitener1D').sets;
197 end
198
199 %--------------------------------------------------------------------------
200 % Get Default Plist
201 %--------------------------------------------------------------------------
202 function plout = getDefaultPlist(set)
203 persistent pl;
204 persistent lastset;
205 if ~exist('pl', 'var') || isempty(pl) || ~strcmp(lastset, set)
206 pl = buildplist(set);
207 lastset = set;
208 end
209 plout = pl;
210 end
211
212 function pl = buildplist(set)
213
214 pl = plist();
215
216 % Append sets of parameters according to the chosen spectral estimator
217 if ~utils.helper.ismember(lower(SETS), lower(set))
218 error('### Unknown set [%s]', set);
219 else
220 pl = ao.getInfo('buildWhitener1D', lower(set)).plists;
221 end
222
223 switch lower(set)
224 case 'default'
225 % Model
226 p = param({'model', ['A frequency-series AO describing the model<br>'...
227 'response to build the filter from. <br>' ...
228 'As an alternative, the user '...
229 'can choose a model estimation technique:<br>'...
230 '<li>PSD - using <tt>psd</tt> + <tt>bin_data</tt></li>'...
231 '<li>LPSD - using <tt>lpsd</tt></li>']}, paramValue.EMPTY_DOUBLE);
232 pl = combine(plist(p), pl);
233 otherwise
234 end
235
236 p = param({'scaleOut', ['Scale your output by the inband power']},paramValue.FALSE_TRUE);
237 pl = combine(plist(p), pl);
238
239 p = param({'flim', ['Band to calculate the scaling power']},[1e-3 30e-3]);
240 pl = combine(plist(p), pl);
241
242
243 end
244