Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/getdof.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 % GETDOF Calculates degrees of freedom for psd, lpsd, cohere and lcohere
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: GETDOF Input psd or mscohere (magnitude square coherence)
5 % estimated with the WOSA (Welch's Overlapped Segment Averaging Method) and
6 % return degrees of freedom of the estimator.
7 %
8 % CALL: dof = getdof(a,pl)
9 %
10 % INPUTS:
11 % a - input analysis objects containing power spectral
12 % densities or magnitude squared coherence.
13 % pl - input parameter list
14 %
15 % OUTPUTS:
16 % dof - cdata AO with degrees of freedom for the
17 % corresponding estimator. If the estimators are lpsd
18 % or lcohere then dof number of elements is the same of
19 % the spectral estimator
20 %
21 %
22 % If the last input argument is a parameter list (plist).
23 % The following parameters are recognised.
24 %
25 %
26 %
27 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'getdof')">Parameters Description</a>
28 %
29 % VERSION: $Id: getdof.m,v 1.19 2011/05/23 20:36:47 mauro Exp $
30 %
31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32
33 function varargout = getdof(varargin)
34
35 %%% Check if this is a call for parameters
36 if utils.helper.isinfocall(varargin{:})
37 varargout{1} = getInfo(varargin{3});
38 return
39 end
40
41 import utils.const.*
42 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
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
49 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
50 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
51
52 %%% avoid multiple AO at input
53 if numel(as) > 1
54 error('!!! Too many input AOs, GETDOF can process only one AO per time !!!')
55 end
56
57 %%% check that fsdata is input
58 if ~isa(as.data, 'fsdata')
59 error('!!! Non-fsdata input, GETDOF can process only fsdata !!!')
60 end
61
62 %%% avoid input modification
63 if nargout == 0
64 error('!!! GETDOF cannot be used as a modifier. Please give an output variable !!!');
65 end
66
67 %%% Parse plists
68 pl = parse(pl, getDefaultPlist());
69
70 %%% Find parameters
71 mtd = lower(find(pl, 'method'));
72 if ~ischar(mtd)
73 error('!!! Method must be a string !!!')
74 end
75 mtd = lower(mtd);
76
77 Ntot = find(pl,'DataLength');
78
79 %%% switching over methods
80 switch mtd
81 case 'psd'
82 % get hist
83 hst = as.hist;
84 % get nodes
85 [n,a, nodes] = getNodes(hst);
86 % get plists from nodes
87 pls = [nodes(:).pl];
88 if numel(pls) > 1
89 plss = pls(1);
90 for ii = 2:numel(pls)-1
91 plss = parse(plss, pls(ii));
92 end
93 end
94 % get number of averages
95 navs = as.data.navs;
96 % get window object
97 w = find(plss, 'WIN');
98 % percentage of overlap
99 olap = find(plss, 'OLAP')./100;
100 % number of bins in each fft
101 nfft = find(plss, 'NFFT');
102 psll = find(plss, 'psll');
103
104 % Normalize window data in order to be square integrable to 1
105 if ischar(w)
106 switch lower(w)
107 case 'kaiser'
108 Win = specwin(w, nfft, psll);
109 otherwise
110 Win = specwin(w, nfft);
111 end
112 else
113 Win = w;
114 end
115 win = Win.win ./ sqrt(Win.ws2);
116
117 % Calculates total number of data in the original time-series
118 if isempty(Ntot)
119 Ntot = ceil(navs*(nfft-olap*nfft)+olap*nfft);
120 end
121
122 if navs == 1
123 dofs = round(2*navs);
124 else
125 [R,n] = utils.math.overlapCorr(win,Ntot,navs);
126 dof = 2*navs/(2*R*navs+1);
127 dofs = round(dof);
128 end
129
130 case 'lpsd'
131 % get hist
132 hst = as.hist;
133 % get nodes
134 [n,a, nodes] = getNodes(hst);
135 % get plists from nodes
136 pls = [nodes(:).pl];
137 if numel(pls) > 1
138 plss = pls(1);
139 for ii = 2:numel(pls)-1
140 plss = parse(plss, pls(ii));
141 end
142 end
143 % get window used
144 uwin = find(plss, 'WIN');
145
146 % extract number of frequencies bins
147 nf = length(as.x);
148
149 % dft length for each bin
150 if ~isempty(as.procinfo)
151 L = as.procinfo.find('L');
152 else
153 error('### The AO doesn''t have any procinfo with the key ''L''');
154 end
155
156 % set original data length as the length of the first window
157 if isempty(Ntot)
158 nx = L(1);
159 else
160 nx = Ntot;
161 end
162
163 % windows overlap
164 olap = find(plss, 'OLAP')./100;
165 psll = find(plss, 'psll');
166
167 dofs = ones(nf,1);
168 for jj = 1:nf
169 l = L(jj);
170 % compute window
171 if ischar(uwin)
172 switch lower(uwin)
173 case 'kaiser'
174 w = specwin(uwin, l, psll);
175 otherwise
176 w = specwin(uwin, l);
177 end
178 else
179 w = uwin;
180 end
181 % Normalize window data in order to be square integrable to 1
182 owin = w.win;
183 owin = owin./sqrt(w.ws2);
184
185 % Compute the number of averages we want here
186 segLen = l;
187 nData = nx;
188 ovfact = 1 / (1 - olap);
189 davg = (((nData - segLen)) * ovfact) / segLen + 1;
190 navg = round(davg);
191
192 if navg == 1
193 dof = 2*navg;
194 else
195 [R,n] = utils.math.overlapCorr(owin,nx,navg);
196 dof = 2*navg/(2*R*navg+1);
197 end
198
199 % storing c and dof
200 dofs(jj) = dof;
201
202 end % for jj=1:nf
203
204 case 'mscohere'
205 % get hist
206 hst = as.hist;
207 % get nodes
208 [n,a, nodes] = getNodes(hst);
209 % get plists from nodes
210 pls = [nodes(:).pl];
211 if numel(pls) > 1
212 plss = pls(1);
213 for ii = 2:numel(pls)-1
214 plss = parse(plss, pls(ii));
215 end
216 end
217 % get number of averages
218 navs = as.data.navs;
219 % get window object
220 w = find(plss,'WIN');
221 % percentage of overlap
222 olap = find(plss,'OLAP')./100;
223 % number of bins in each fft
224 nfft = find(plss,'NFFT');
225 psll = find(plss, 'psll');
226
227 % Normalize window data in order to be square integrable to 1
228 if ischar(w)
229 switch lower(w)
230 case 'kaiser'
231 Win = specwin(w, nfft, psll);
232 otherwise
233 Win = specwin(w, nfft);
234 end
235 else
236 Win = w;
237 end
238 win = Win.win ./ sqrt(Win.ws2);
239
240 % Calculates total number of data in the original time-series
241 if isempty(Ntot)
242 Ntot = ceil(navs*(nfft-olap*nfft)+olap*nfft);
243 end
244
245 if navs == 1
246 dofs = round(2*navs);
247 else
248 [R,n] = utils.math.overlapCorr(win,Ntot,navs);
249 dof = 2*navs/(2*R*navs+1);
250 dofs = round(dof);
251 end
252
253 case 'mslcohere'
254 % get hist
255 hst = as.hist;
256 % get nodes
257 [n,a, nodes] = getNodes(hst);
258 % get plists from nodes
259 pls = [nodes(:).pl];
260 if numel(pls) > 1
261 plss = pls(1);
262 for ii = 2:numel(pls)-1
263 plss = parse(plss, pls(ii));
264 end
265 end
266 % get window used
267 uwin = find(plss,'WIN');
268
269 % extract number of frequencies bins
270 nf = length(as.x);
271
272 % dft length for each bin
273 if ~isempty(as.procinfo)
274 L = as.procinfo.find('L');
275 else
276 error('### The AO doesn''t have any procinfo with the key ''L''');
277 end
278
279 % set original data length as the length of the first window
280 if isempty(Ntot)
281 nx = L(1);
282 else
283 nx = Ntot;
284 end
285
286 % windows overlap
287 olap = find(plss, 'OLAP')./100;
288 psll = find(plss, 'psll');
289
290 dofs = ones(nf, 1);
291 for jj = 1:nf
292 l = L(jj);
293 % compute window
294 if ischar(uwin)
295 switch lower(uwin)
296 case 'kaiser'
297 w = specwin(uwin, l, psll);
298 otherwise
299 w = specwin(uwin, l);
300 end
301 else
302 w = uwin;
303 end
304 % Normalize window data in order to be square integrable to 1
305 owin = w.win ./ sqrt(w.ws2);
306
307 % Compute the number of averages we want here
308 segLen = l;
309 nData = nx;
310 ovfact = 1 / (1 - olap);
311 davg = (((nData - segLen)) * ovfact) / segLen + 1;
312 navg = round(davg);
313
314 if navg == 1
315 dof = 2*navg;
316 else
317 [R,n] = utils.math.overlapCorr(owin,nx,navg);
318 dof = 2*navg/(2*R*navg+1);
319 end
320
321 % storing c and dof
322 dofs(jj) = dof;
323
324 end % for jj=1:nf
325
326 end %switch mtd
327
328 % Output data
329
330 % dof
331 ddof = cdata();
332 ddof.setY(dofs);
333 odof = ao(ddof);
334 % Set output AO name
335 odof.name = sprintf('dof(%s)', ao_invars{:});
336 % Add history
337 odof.addHistory(getInfo('None'), pl, [ao_invars(:)], [as.hist]);
338
339 % output
340 if nargout == 1
341 varargout{1} = odof;
342 else
343 error('!!! getdof can have only one output')
344 end
345
346 end
347
348 %--------------------------------------------------------------------------
349 % Get Info Object
350 %--------------------------------------------------------------------------
351 function ii = getInfo(varargin)
352 if nargin == 1 && strcmpi(varargin{1}, 'None')
353 sets = {};
354 pl = [];
355 else
356 sets = {'Default'};
357 pl = getDefaultPlist();
358 end
359 % Build info object
360 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: getdof.m,v 1.19 2011/05/23 20:36:47 mauro Exp $', sets, pl);
361 ii.setModifier(false);
362 ii.setOutmin(1);
363 ii.setOutmax(1);
364 end
365
366 %--------------------------------------------------------------------------
367 % Get Default Plist
368 %--------------------------------------------------------------------------
369 function plout = getDefaultPlist()
370 persistent pl;
371 if ~exist('pl', 'var') || isempty(pl)
372 pl = buildplist();
373 end
374 plout = pl;
375 end
376
377 function pl = buildplist()
378
379 pl = plist();
380
381 p = param({'method', ['Set the desired method. Supported values are<ul>'...
382 '<li>''psd'' power spectrum calculated with ao/psd, whatever the scale</li>'...
383 '<li>''lpsd'' power spectrum calculated with ao/lpsd, whatever the scale</li>'...
384 '<li>''mscohere'' magnitude square coherence spectrum calculated with ao/cohere</li>'...
385 '<li>''mslcohere'' magnitude square coherence spectrum calculated with ao/lcohere</li>']}, ...
386 {1, {'psd', 'lpsd', 'mscohere', 'mslcohere'}, paramValue.OPTIONAL});
387 pl.append(p);
388
389 p = param({'DataLength',['Data length of the time series.'...
390 'It is better to input for more stable calculation.'...
391 'Leave it empty if you do not know its value.']},...
392 paramValue.EMPTY_DOUBLE);
393 pl.append(p);
394
395
396 end
397
398
399
400 % END
401
402