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