Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/.#psd.m.1.68 @ 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 % PSD makes power spectral density estimates of the time-series objects
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: PSD makes power spectral density estimates of the
5 % time-series objects in the input analysis objects
6 % using the Welch Overlap method. PSD is computed
7 % using a modified version of MATLAB's welch (>> help welch).
8 %
9 % CALL: bs = psd(a1,a2,a3,...,pl)
10 % bs = psd(as,pl)
11 % bs = as.psd(pl)
12 %
13 % INPUTS: aN - input analysis objects
14 % as - input analysis objects array
15 % pl - input parameter list
16 %
17 % OUTPUTS: bs - array of analysis objects, one for each input
18 %
19 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'psd')">Parameters Description</a>
20 %
21 % VERSION: $Id: psd.m,v 1.68 2011/04/27 05:41:08 mauro Exp $
22 %
23 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
24
25
26 function varargout = psd(varargin)
27
28 % Check if this is a call for parameters
29 if utils.helper.isinfocall(varargin{:})
30 varargout{1} = getInfo(varargin{3});
31 return
32 end
33
34 import utils.const.*
35 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
36
37 % Collect input variable names
38 in_names = cell(size(varargin));
39 for ii = 1:nargin,in_names{ii} = inputname(ii);end
40
41 % Collect all AOs
42 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
43
44 % Decide on a deep copy or a modify
45 bs = copy(as, nargout);
46
47 % Apply defaults to plist
48 usepl = applyDefaults(getDefaultPlist, varargin{:});
49
50 inhists = [];
51
52 % Loop over input AOs
53 for jj = 1 : numel(bs)
54 % gather the input history objects
55 inhists = [inhists bs(jj).hist];
56
57 % check input data
58 if isa(bs(jj).data, 'tsdata')
59
60 % Check the time range.
61 time_range = mfind(usepl, 'split', 'times');
62 if ~isempty(time_range)
63 switch class(time_range)
64 case 'double'
65 bs(jj) = split(bs(jj), plist(...
66 'times', time_range));
67 case 'timespan'
68 bs(jj) = split(bs(jj), plist(...
69 'timespan', time_range));
70 case 'time'
71 bs(jj) = split(bs(jj), plist(...
72 'start_time', time_range(1), ...
73 'end_time', time_range(2)));
74 case 'cell'
75 bs(jj) = split(bs(jj), plist(...
76 'start_time', time_range{1}, ...
77 'end_time', time_range{2}));
78 otherwise
79 end
80 end
81
82 % Check the length of the object
83 if bs(jj).len <= 0
84 error('### The object is empty! Please revise your settings ...');
85 end
86
87 % Utility to deal with nfft, win, olap etc
88 use_pl = utils.helper.process_spectral_options(usepl, 'lin', len(bs(jj)), bs(jj).fs);
89
90 % Compute PSD using pwelch
91 [pxx, f, info, dev] = welch(bs(jj), 'psd', use_pl);
92
93 % Keep the data shape of the input AO
94 if size(bs(jj).data.y,1) == 1
95 pxx = pxx.';
96 f = f.';
97 end
98 % create new output fsdata
99 fs = bs(jj).fs;
100 fsd = fsdata(f, pxx, fs);
101 fsd.setXunits('Hz');
102 fsd.setYunits(info.units);
103 fsd.setEnbw(info.enbw);
104 fsd.setNavs(info.navs);
105 fsd.setT0(bs(jj).data.t0);
106 % update AO
107 bs(jj).data = fsd;
108 % add std deviation
109 bs(jj).data.dy = dev;
110 % set name: scaling of spectrum
111 scale = upper(find(use_pl, 'Scale'));
112 bs(jj).name = sprintf('%s(%s)', lower(scale), ao_invars{jj});
113 % Add history
114 bs(jj).addHistory(getInfo('None'), use_pl, ao_invars(jj), inhists(jj));
115 else
116 warning('### Ignoring input AO number %d (%s); it is not a time-series.', jj, bs(jj).name)
117 end
118 end % loop over analysis objects
119
120 % Set output
121 varargout = utils.helper.setoutputs(nargout, bs);
122 end
123
124 %--------------------------------------------------------------------------
125 % Get Info Object
126 %--------------------------------------------------------------------------
127 function ii = getInfo(varargin)
128 if nargin == 1 && strcmpi(varargin{1}, 'None')
129 sets = {};
130 pl = [];
131 else
132 sets = {'Default'};
133 pl = getDefaultPlist();
134 end
135 % Build info object
136 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: psd.m,v 1.68 2011/04/27 05:41:08 mauro Exp $', sets, pl);
137 end
138
139 %--------------------------------------------------------------------------
140 % Get Default Plist
141 %--------------------------------------------------------------------------
142 function plout = getDefaultPlist()
143 persistent pl;
144 if ~exist('pl', 'var') || isempty(pl)
145 pl = buildplist();
146 end
147 plout = pl;
148 end
149
150 function pl = buildplist()
151
152 % General plist for Welch-based, linearly spaced spectral estimators
153 pl = plist.WELCH_PLIST;
154
155 % Scale
156 p = param({'Scale',['The scaling of output. Choose from:<ul>', ...
157 '<li>PSD - Power Spectral Density</li>', ...
158 '<li>ASD - Amplitude (linear) Spectral Density</li>', ...
159 '<li>PS - Power Spectrum</li>', ...
160 '<li>AS - Amplitude (linear) Spectrum</li></ul>']}, {1, {'PSD', 'ASD', 'PS', 'AS'}, paramValue.SINGLE});
161 pl.append(p);
162
163 end
164