Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/xspec.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 bc767aaa99a8
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % XSPEC performs cross-spectral analysis of various forms.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: XSPEC performs cross-spectral analysis of various forms.
5 % The function is a helper function for various higher level
6 % functions. It is meant to be called from other functions
7 % (e.g., tfe).
8 %
9 % CALL: b = xspec(a, pl, method, iALGO, iVER, invars);
10 %
11 % INPUTS: a - vector of 2 input AOs
12 % pl - input parameter list
13 % method - one of
14 % 'cpsd' - compute cross-spectral density
15 % 'tfe' - estimate transfer function between inputs
16 % 'mscohere' - estimate magnitude-squared cross-coherence
17 % 'cohere' - estimate complex cross-coherence
18 % iALGO - ALGONAME from the calling higher level script
19 % iVER - VERSION from the calling higher level script
20 % invars - invars variable from the calling higher level script
21 %
22 % OUTPUTS: b - output AO
23 %
24 % VERSION: $Id: xspec.m,v 1.61 2011/05/22 22:09:44 mauro Exp $
25 %
26 % HISTORY: 30-05-2007 M Hewitson
27 % Creation
28 %
29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30
31 function varargout = xspec(varargin)
32
33 import utils.const.*
34
35 % unpack inputs
36 as = varargin{1};
37 pl = varargin{2};
38 method = varargin{3};
39 mi = varargin{4};
40 invars = varargin{5};
41
42 VERSION = '$Id: xspec.m,v 1.61 2011/05/22 22:09:44 mauro Exp $';
43
44 % Set the method version string in the minfo object
45 mi.setMversion([VERSION '-->' mi.mversion]);
46
47 %----------------- Select all AOs with time-series data
48 tsao = [];
49 for ll = 1:numel(as)
50 if isa(as(ll).data, 'tsdata')
51 tsao = [tsao as(ll)];
52 else
53 warning('### xspec requires tsdata (time-series) inputs. Skipping AO %s. \nREMARK: The output doesn''t contain this AO', invars{ll});
54 end
55 end
56 if numel(tsao) ~= 2
57 error('### XSPEC needs two time-series AOs.');
58 end
59
60 %----------------- Gather the input history objects
61 inhists = [tsao(:).hist];
62
63 %----------------- Check the time range
64 time_range = mfind(pl, 'split', 'times');
65 for ll = 1:numel(tsao)
66 if ~isempty(time_range)
67 switch class(time_range)
68 case 'double'
69 tsao(ll) = split(tsao(ll), plist(...
70 'times', time_range));
71 case 'timespan'
72 tsao(ll) = split(tsao(ll), plist(...
73 'timespan', time_range));
74 case 'time'
75 tsao(ll) = split(tsao(ll), plist(...
76 'start_time', time_range(1), ...
77 'end_time', time_range(2)));
78 case 'cell'
79 tsao(ll) = split(tsao(ll), plist(...
80 'start_time', time_range{1}, ...
81 'end_time', time_range{2}));
82 otherwise
83 end
84 end
85 if tsao(ll).len <= 0
86 error('### The object is empty! Please revise your settings ...');
87 end
88 end
89
90 %----------------- Resample all AOs
91 copies = zeros(size(tsao));
92
93 fsmin = ao.findFsMin(tsao);
94 fspl = plist(param('fsout', fsmin));
95 for ll = 1:numel(tsao)
96 % Check Fs
97 if tsao(ll).data.fs ~= fsmin
98 utils.helper.msg(msg.PROC1, 'resampling AO %s to %f Hz', tsao(ll).name, fsmin);
99 % Make a deep copy so we don't
100 % affect the original input data
101 tsao(ll) = copy(tsao(ll), 1);
102 copies(ll) = 1;
103 resample(tsao(ll), fspl);
104 end
105 end
106
107 %----------------- Truncate all vectors
108
109 % Get shortest vector
110 lmin = ao.findShortestVector(tsao);
111 nsecs = lmin / fsmin;
112 for ll = 1:numel(tsao)
113 if len(tsao(ll)) ~= lmin
114 utils.helper.msg(msg.PROC2, 'truncating AO %s to %d secs', tsao(ll).name, nsecs);
115 % do we already have a copy?
116 if ~copies(ll)
117 % Make a deep copy so we don't
118 % affect the original input data
119 tsao(ll) = copy(tsao(ll), 1);
120 copies(ll) = 1;
121 end
122 tsao(ll).select(1:lmin);
123 end
124 end
125
126 %----------------- check input parameters
127
128 usepl = utils.helper.process_spectral_options(pl, 'lin', lmin, fsmin);
129
130 % Loop over input AOs
131 utils.helper.msg(msg.PROC1, 'computing %s(%s -> %s)', method, invars{1}, invars{2});
132
133 % -------- Make Xspec estimate
134
135 % Compute xspec using welch and always scale to PSD
136 [txy, f, info, dev] = welch(tsao(1), tsao(2), method, usepl.pset('Scale', 'PSD'));
137
138 % Keep the data shape of the first input AO
139 if size(tsao(1).data.y,1) == 1
140 txy = txy.';
141 f = f.';
142 dev = dev.';
143 end
144
145 % create new output fsdata
146 fsd = fsdata(f, txy, fsmin);
147 fsd.setEnbw(info.enbw);
148 fsd.setXunits('Hz');
149
150 switch lower(method)
151 case 'tfe'
152 fsd.setYunits(tsao(2).data.yunits / tsao(1).data.yunits);
153 case 'cpsd'
154 fsd.setYunits(tsao(2).data.yunits * tsao(1).data.yunits / unit('Hz'));
155 case {'mscohere','cohere'}
156 fsd.setYunits(unit());
157 otherwise
158 error(['### Unknown method:' method]);
159 end
160 fsd.setNavs(info.navs);
161
162 % set object t0
163 if eq(tsao(1).t0, tsao(2).t0)
164 fsd.setT0(tsao(1).t0);
165 end
166
167 % make output analysis object
168 bs = ao(fsd);
169 % add variance
170 bs.data.dy = dev;
171
172 % simplify the units in the case of cpsd calculation
173 if strcmp(method, 'cpsd')
174 bs.simplifyYunits(...
175 plist('prefixes', false, 'exceptions', 'Hz'));
176 end
177
178 %----------- Add history
179
180 % set name
181 bs.name = sprintf('%s(%s->%s)', upper(method), invars{1}, invars{2});
182
183 % we need to get the input histories in the same order tsao the inputs
184 % to this function call, not in the order of the input to xspec;
185 % otherwise the resulting matrix on a 'create from history' will be
186 % mirrored.
187 bs.addHistory(mi, usepl, [invars(:)], inhists);
188
189 % Propagate 'plotinfo'
190 if ~isempty(tsao(1).plotinfo) || ~isempty(tsao(2).plotinfo)
191 bs.plotinfo = combine(tsao(1).plotinfo, tsao(2).plotinfo);
192 end
193
194 % Set output
195 varargout{1} = bs;
196 end
197