Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/xcorr.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 % XCORR makes cross-correlation estimates of the time-series
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: XCORR makes cross-correlation estimates of the time-series
5 % objects in the input analysis objects. The cross-correlation is
6 % computed using MATLAB's xcorr (>> help xcorr).
7 %
8 % CALL: b = xcorr(a1,a2,pl)
9 %
10 % INPUTS: b - output analysis objects
11 % a1,a2 - input analysis objects (only two)
12 % pl - input parameter list
13 %
14 % The function makes correlation estimates between a1 and a2.
15 %
16 % If only on AO is input, the auto-correlation is computed.
17 %
18 % If the last input argument is a parameter list (plist) it is used.
19 % The following parameters are recognised.
20 %
21 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'xcorr')">Parameters Description</a>
22 %
23 % VERSION: $Id: xcorr.m,v 1.23 2011/04/08 08:56:15 hewitson Exp $
24 %
25 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26
27 function varargout = xcorr(varargin)
28
29 % Check if this is a call for parameters
30 if utils.helper.isinfocall(varargin{:})
31 varargout{1} = getInfo(varargin{3});
32 return
33 end
34
35 import utils.const.*
36 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
37
38 if nargout == 0
39 error('### xcorr cannot be used as a modifier. Please give an output variable.');
40 end
41
42 % Collect input variable names
43 in_names = cell(size(varargin));
44 for ii = 1:nargin,in_names{ii} = inputname(ii);end
45
46 % Collect all AOs and plists
47 [as, invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
48 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
49
50 % combine plists
51 pl = parse(pl, getDefaultPlist());
52
53 na = numel(as);
54 if na > 2
55 error('### XCORR accepts only two AOs to cross-correlate.');
56 end
57
58 %----------------- Keep the history to suppress the history of the
59 %----------------- intermediate steps
60 inhists = [as(:).hist];
61
62 %----------------- Resample all AOs
63 copies = zeros(size(as));
64
65 fsmax = findFsMax(as);
66 fspl = plist('fsout', fsmax);
67 for jj=1:na
68 % check this is a time-series object
69 if ~isa(as(jj).data, 'tsdata')
70 error('### ltpda_xspec requires tsdata (time-series) inputs.');
71 end
72 % Check Fs
73 if as(jj).fs ~= fsmax
74 utils.helper.msg(msg.PROC1, 'resampling AO %s to %f Hz', as(jj).name, fsmax);
75 % Make a deep copy so we don't
76 % affect the original input data
77 as(jj) = copy(as(jj), 1);
78 copies(jj) = 1;
79 as(jj).resample(fspl);
80 end
81 end
82
83 %----------------- Truncate all vectors
84
85 % Get shortest vector
86 utils.helper.msg(msg.PROC1, '*** Truncating all vectors...');
87 lmin = findShortestVector(as);
88 nsecs = lmin / fsmax;
89 for jj=1:na
90 if len(as(jj)) ~= lmin
91 utils.helper.msg(msg.PROC2, 'truncating AO %s to %d secs', as(jj).name, nsecs);
92 % do we already have a copy?
93 if ~copies(jj)
94 % Make a deep copy so we don't
95 % affect the original input data
96 as(jj) = copy(as(jj), 1);
97 copies(jj) = 1;
98 end
99 as(jj).select(1:lmin);
100 end
101 end
102
103 %----------------- check input parameters
104
105 % Maximum lag for Xcorr
106 MaxLag = find(pl, 'MaxLag');
107
108 % Scale for Xcorr
109 scale = find(pl, 'Scale');
110
111 % Loop over input AOs
112 bs = ao;
113
114 % -------- Make Xspec estimate
115
116 % Compute cross-correlation estimates using XCORR
117 if MaxLag == -1
118 MaxLag = len(as(1));
119 end
120 % Use .data.y syntax (rather than .y) to preserve y vector shape
121 [c,lags] = xcorr(as(1).data.y, as(2).data.y, MaxLag, scale);
122
123 % Keep the data shape of the first input AO
124 if size(as(1).y,1) == 1
125 c = c.';
126 end
127
128 % create new output xydata
129 xy = xydata(lags./fsmax, c);
130 xy.setXunits('s');
131 switch scale
132 case {'none', 'biased', 'unbiased'}
133 xy.setYunits(as(1).yunits * as(2).yunits);
134 case 'coeff'
135 xy.setYunits('');
136 otherwise
137 error(['Unsupported scaling option ' scale]);
138 end
139
140
141 %----------- create new output history
142
143 % make output analysis object
144 bs.data = xy;
145 % set name
146 bs.name = sprintf('xcorr(%s->%s)', invars{1}, invars{2});
147 % Propagate 'plotinfo'
148 plotinfo = [as(:).plotinfo];
149 if ~isempty(plotinfo)
150 bs.plotinfo = combine(plotinfo);
151 end
152 % we need to get the input histories in the same order as the inputs
153 % to this function call, not in the order of the input to tfestimate;
154 % otherwise the resulting matrix on a 'create from history' will be
155 % mirrored.
156 bs.addHistory(getInfo('None'), pl, [invars(:)], inhists);
157
158 % Set output
159 varargout{1} = bs;
160 % end
161 end
162
163 %--------------------------------------------------------------------------
164 % Returns the length of the shortest vector in samples
165 function lmin = findShortestVector(as)
166 lmin = 1e20;
167 for jj=1:numel(as)
168 if len(as(jj)) < lmin
169 lmin = len(as(jj));
170 end
171 end
172 end
173 %--------------------------------------------------------------------------
174 % Returns the max Fs of a set of AOs
175 function fs = findFsMax(as)
176 fs = 0;
177 for jj=1:numel(as)
178 a = as(jj);
179 if a.fs > fs
180 fs = a.fs;
181 end
182 end
183 end
184
185 %--------------------------------------------------------------------------
186 % Get Info Object
187 %--------------------------------------------------------------------------
188 function ii = getInfo(varargin)
189 if nargin == 1 && strcmpi(varargin{1}, 'None')
190 sets = {};
191 pl = [];
192 else
193 sets = {'Default'};
194 pl = getDefaultPlist;
195 end
196 % Build info object
197 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: xcorr.m,v 1.23 2011/04/08 08:56:15 hewitson Exp $', sets, pl);
198 ii.setModifier(false);
199 ii.setArgsmin(2);
200 end
201
202 %--------------------------------------------------------------------------
203 % Get Default Plist
204 %--------------------------------------------------------------------------
205 function plout = getDefaultPlist()
206 persistent pl;
207 if exist('pl', 'var')==0 || isempty(pl)
208 pl = buildplist();
209 end
210 plout = pl;
211 end
212
213 function pl = buildplist()
214
215 pl = plist();
216
217 % MaxLag
218 p = param({'MaxLag', 'Compute over a range of lags -MaxLag to MaxLag [default: M-1]'}, {1, {-1}, paramValue.OPTIONAL});
219 pl.append(p);
220
221 % Scale
222 p = param({'Scale', ['normalisation of the correlation. Choose from:<ul>'...
223 '<li>''biased'' - scales the raw cross-correlation by 1/M</li>'...
224 '<li>''unbiased'' - scales the raw correlation by 1/(M-abs(lags))</li>'...
225 '<li>''coeff'' - normalizes the sequence so that the auto-correlations<br>'...
226 'at zero lag are identically 1.0.</li>'...
227 '<li>''none'' - no scaling</li></ul>']}, {1, {'none', 'biased', 'unbiased', 'coeff'}, paramValue.SINGLE});
228 pl.append(p);
229
230 end
231 % PARAMETERS: 'MaxLag' - compute over range of lags -MaxLag to MaxLag [default: M-1]
232 % 'Scale' - normalisation of the correlation. Choose from:
233 % 'biased' - scales the raw cross-correlation by 1/M.
234 % 'unbiased' - scales the raw correlation by 1/(M-abs(lags)).
235 % 'coeff' - normalizes the sequence so that the auto-correlations
236 % at zero lag are identically 1.0.
237 % 'none' - no scaling (this is the default).
238 %
239 % M is the length of longest input vector. If one vector is shorted,
240 % it is zero padded.