Mercurial > hg > ltpda

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