Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/quasiSweptSine.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 % QUASISWEPTSING computes a transfer function from swept-sine measurements
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: QUASISWEPTSING computes a transfer function from discrete
5 % swept-sine measurements.
6 %
7 % In order for the calculation to work, you need to give it an array of
8 % start and stop times (or durations), and (optionally) an array of
9 % amplitudes and frequencies of the injected sine-waves. If you don't
10 % specify the frequencies, you must give a time-series of the injected
11 % signal and the algorithm will try to determine the amplitudes and
12 % frequencies from the data.
13 %
14 %
15 % CALL: T = quasiSweptSine(out, pl);
16 %
17 % INPUTS: out - The measured output of the system
18 % PL - parameter list
19 %
20 % OUTPUT: T - the measured transfer function
21 %
22 % The procinfo of the output AOs contains the following fields:
23 %
24 % 'frequencies' - the frequencies used in the DFT estimation.
25 % 'timespans' - an array of timespan objects, one for each sine-wave
26 % segment
27 %
28 %
29 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'quasiSweptSine')">Parameters Description</a>
30 %
31 % VERSION: $Id: quasiSweptSine.m,v 1.10 2011/04/08 08:56:16 hewitson Exp $
32 %
33 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34
35
36 function varargout = quasiSweptSine(varargin)
37
38 % check if this is a call for parameters
39 if utils.helper.isinfocall(varargin{:})
40 varargout{1} = getInfo(varargin{3});
41 return
42 end
43
44 % tell the system we are runing
45 import utils.const.*
46 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
47
48 % collect input variable names
49 in_names = cell(size(varargin));
50 for ii = 1:nargin,in_names{ii} = inputname(ii);end
51
52 % collect all AOs and plists
53 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
54 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
55
56 if nargout == 0
57 error('### quasiSweptSine can not be used as a modifier method. Please give at least one output');
58 end
59
60 % Make copies or handles to inputs
61 bs = copy(as, nargout);
62
63 % combine plists
64 pl = parse(pl, getDefaultPlist());
65
66 % Parameters
67 input = pl.find('input');
68 startTimes = pl.find('Start times');
69 stopTimes = pl.find('Stop times');
70 durations = pl.find('durations');
71 amplitudes = pl.find('amplitudes');
72 frequencies = pl.find('frequencies');
73 % phases = pl.find('phases');
74 inUnits = pl.find('Input Units');
75 win = pl.find('Win');
76 Nerror = pl.find('Nerror');
77
78 if isa(input, 'ao') && ~isa(input.data, 'tsdata')
79 utils.helper.err('quasiSweptSine requires time-series as input');
80 end
81
82
83 %---------------------------------------------------
84 %--------- Convert the times to timespans
85 if isempty(durations) && isempty(stopTimes)
86 utils.helper.err('You need to specify either an array of stop times, or an array of durations');
87 end
88
89
90 %---------------------------------------------------
91 %--------- if we have the input, we use it to determine amplitudes and
92 %--------- frequencies
93
94 if isempty(input) && isempty(amplitudes)
95 utils.helper.err('You need to specify either an input signal, or a full description of the signals including amplitudes and frequecies.');
96 end
97
98 computeFrequecies = false;
99 if isempty(frequencies)
100 computeFrequecies = true;
101 end
102
103
104 % Go through each ao
105 for jj=1:numel(bs)
106
107
108 if ~isa(bs.data, 'tsdata')
109 utils.helper.err('quasiSweptSine requires time-series as input');
110 end
111
112 [timespans, startTimes, stopTimes] = generateTimespans(bs(jj).t0, startTimes, stopTimes, durations);
113
114
115 % Go through each time-span
116 Txx = zeros(size(timespans));
117 dT = zeros(size(timespans));
118 for kk=1:numel(timespans)
119 if isa(input, 'ao')
120 inseg = input.split(plist('timespan', timespans(kk)));
121 else
122 % We don't have an input, so we need to create inputs from the
123 % signal specs
124
125 inseg = ao(plist('waveform', 'sine wave', ...
126 'nsecs', timespans(kk).interval, ...
127 'fs', bs(jj).fs, ...
128 'A', amplitudes(kk), ...
129 'f', frequencies(kk), ...
130 'toff', 0, ...
131 't0', startTimes(kk), ...
132 'yunits', inUnits));
133 end
134 if computeFrequecies
135 b = sineParams(inseg, plist('N', 1));
136 frequencies(kk) = b.y(2);
137 end
138
139 utils.helper.msg(msg.PROC1, 'Computing TF at %g Hz', frequencies(kk));
140
141 % Window
142 w = ao(plist('win', win, 'length', inseg.len));
143
144 % Output
145 outseg = bs(jj).split(plist('timespan', timespans(kk)));
146
147 % Compute DFT
148 fs = outseg.data.fs;
149 N = outseg.len;
150 J = -2*pi*1i.*(0:N-1)/fs;
151 outxx = exp(frequencies(kk)*J)*(w.y.*outseg.data.getY);
152 inxx = exp(frequencies(kk)*J)*(w.y.*inseg.data.getY);
153 Txx(kk) = outxx./inxx;
154
155 % Compute error
156 %
157 nout = noiseAroundLine(outseg, frequencies(kk), Nerror);
158 nin = noiseAroundLine(inseg, frequencies(kk), Nerror);
159 snr1 = abs(outxx)./nout;
160 snr2 = abs(inxx)./nin;
161 dT(kk) = abs(Txx(kk)) * sqrt( (1./snr1)^2 + (1./snr2)^2);
162
163
164 end % End loop over timespans (segments)
165
166 % Make output
167 bs(jj).data = fsdata(frequencies, Txx, bs(jj).data.fs);
168 bs(jj).data.dy = dT;
169 bs(jj).data.setXunits('Hz');
170 bs(jj).data.setYunits(outseg.yunits./inseg.yunits);
171 % Set name
172 bs(jj).name = sprintf('sweptsine(%s,%s)', input.name, ao_invars{jj});
173 % Set procinfo
174 bs(jj).procinfo = plist('frequencies', frequencies, ...
175 'timespan', timespans);
176 % Add history
177 bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), bs(jj).hist);
178
179 end % Loop over input aos
180
181 % set outputs
182 varargout{1} = bs;
183
184 end
185
186 function n = noiseAroundLine(sig, f0, N)
187
188 xx = abs(fft(sig));
189
190 % get the noise floor around the frequency of interest
191 % - we need the bin that is nearest the frequency
192 f = abs(xx.x - f0);
193 [m, mi] = min(f);
194 % xx.x(mi)
195 % iplot(xx)
196 % plot(xx.x(mi), xx.y(mi), 'x')
197
198 M = N-1;
199 % Measure below the line frequency if we can
200 nl = 0;
201 if (mi>1)
202 ls = max(1, mi-2*M);
203 le = max(1, mi-M);
204 nl = xx.y(ls:le);
205 end
206 % Measure above the line frequency if we can
207 nu = 0;
208 if mi<xx.len
209 us = min(xx.len, mi+M);
210 ue = min(xx.len, mi+2*M);
211 nu = xx.y(us:ue);
212 end
213 n = mean([nl;nu]);
214
215 end
216
217 function [ts, startTimes, stopTimes] = generateTimespans(t0, startTimes, stopTimes, durations)
218
219
220 if isempty(durations) && numel(startTimes) ~= numel(stopTimes)
221 utils.helper.err('You need to specify the same number of start and stop times.');
222 end
223 if isempty(stopTimes) && numel(startTimes) ~= numel(durations)
224 utils.helper.err('You need to specify the same number of durations and stop times.');
225 end
226
227
228 ts = timespan.initObjectWithSize(1,numel(startTimes));
229
230 % Convert to time objects
231 if iscell(startTimes)
232 startTimes = time(startTimes);
233 end
234 if isnumeric(startTimes)
235 newStarts = time.initObjectWithSize(1,numel(startTimes));
236 for kk=1:numel(startTimes)
237 newStarts(kk) = t0+startTimes(kk);
238 end
239 startTimes = newStarts;
240 end
241
242 if isempty(durations) && iscell(stopTimes)
243 stopTimes = time(stopTimes);
244 end
245
246 if isempty(durations) && isnumeric(stopTimes)
247 if isnumeric(stopTimes)
248 newStops = time.initObjectWithSize(1,numel(startTimes));
249 for kk=1:numel(stopTimes)
250 newStops(kk) = t0+stopTimes(kk);
251 end
252 stopTimes = newStops;
253 end
254 end
255
256 useDuration = isempty(stopTimes);
257 for kk=1:numel(startTimes)
258 if useDuration
259 ts(kk) = timespan(startTimes(kk), startTimes(kk)+durations(kk));
260 else
261 ts(kk) = timespan(startTimes(kk), stopTimes(kk));
262 end
263 end
264
265 end
266
267
268 % get info object
269 function ii = getInfo(varargin)
270 if nargin == 1 && strcmpi(varargin{1}, 'None')
271 sets = {};
272 pl = [];
273 else
274 sets = {'Default'};
275 pl = getDefaultPlist;
276 end
277 % build info object
278 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.op, '$Id: quasiSweptSine.m,v 1.10 2011/04/08 08:56:16 hewitson Exp $', sets, pl);
279 end
280
281 % get default plist
282 function plout = getDefaultPlist()
283 persistent pl;
284 if exist('pl', 'var')==0 || isempty(pl)
285 pl = buildplist();
286 end
287 plout = pl;
288 end
289
290 function pl = buildplist()
291
292 % default plist for linear fitting
293 pl = plist();
294
295 % Input
296 p = param({'input', 'The input data series.'}, paramValue.EMPTY_DOUBLE);
297 pl.append(p);
298
299 % Start times
300 p = param({'Start Times', 'A cell array of start times, or an array of time objects.'}, ...
301 paramValue.EMPTY_CELL);
302 pl.append(p);
303
304 % Stop times
305 p = param({'Stop Times', 'A cell array of stop times, or an array of time objects.'}, ...
306 paramValue.EMPTY_CELL);
307 pl.append(p);
308
309 % Durations
310 p = param({'Durations', 'An array of durations that can be used instead of the stop times.'}, ...
311 paramValue.EMPTY_DOUBLE);
312 pl.append(p);
313
314 % Amplitudes
315 p = param({'Amplitudes', 'An array of amplitudes.'}, ...
316 paramValue.EMPTY_DOUBLE);
317 pl.append(p);
318
319 % Frequencies
320 p = param({'Frequencies', 'An array of frequencies [Hz].'}, ...
321 paramValue.EMPTY_DOUBLE);
322 pl.append(p);
323
324 % Input units
325 p = param({'Input Units', 'If you don''t give an input signal AO, you can specify the units of the signal that will be constructed internally.'}, ...
326 {1, {'V'}, paramValue.OPTIONAL});
327 pl.append(p);
328
329 % Window
330 p = param({'Win', 'A window to apply to each segment when computing the DFT.'}, ...
331 paramValue.WINDOW);
332 pl.append(p);
333
334 % Error samples
335 p = param({'Nerror', ['The number of samples either side of the line frequency to use to estimate the noise floor.'...
336 'The noise is estimated from <br><br><tt>mean([y(idx-2*M:idx-M);y(idx+M:idx+2M)])</tt><br><br> where <tt>M=N-1</tt> and <tt>idx</tt> is the index '...
337 'of the bin nearest to the frequency of the signal.']}, ...
338 paramValue.DOUBLE_VALUE(5));
339 pl.append(p);
340
341 % % Phases
342 % p = param({'Phases', 'An array of phases [degrees].'}, paramValue.DOUBLE_VALUE(0));
343 % pl.append(p);
344
345 end