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