bmarechal / octave_scripts

Download as
Browse Code ยป

Commit b197c3fdf011f9972e68675ac4ca59e0cd09ad67

Authored by bmarechal
0 parents
Exists in master

first commit

Showing 14 changed files with 2060 additions and 0 deletions Inline Diff

File was created 1 function [retval, s, errorb, tau] = allan(data,tau,name,verbose)
2 % ALLAN Compute the Allan deviation for a set of time-domain frequency data
3 % [RETVAL, S, ERRORB, TAU] = ALLAN(DATA,TAU,NAME,VERBOSE)
4 %
5 % Inputs:
6 % DATA should be a structure and have the following fields:
7 % DATA.freq or DATA.phase
8 % A vector of fractional frequency measurements (df/f) in
9 % DATA.freq *or* phase offset data (seconds) in DATA.phase .
10 % If frequency data is not present, it will be generated by
11 % differentiating the phase data.
12 % If both fields are present, then DATA.freq will be used.
13 % Note: for general-purpose calculations of Allan deviation,
14 % (i.e. a two-sample variance) use DATA.freq .
15 %
16 % DATA.rate or DATA.time
17 % The sampling rate in Hertz (DATA.rate) or a vector of
18 % timestamps for each measurement in seconds (DATA.time).
19 % DATA.rate is used if both fields are present.
20 % If DATA.rate == 0, then the timestamps are used.
21 %
22 % DATA.units (optional)
23 % The units for the data. If present, the string DATA.units
24 % is added to the plot y-axis label.
25 %
26 % TAU is an array of tau values for computing Allan deviation.
27 % TAU values must be divisible by 1/DATA.rate (data points cannot be
28 % grouped in fractional quantities!) and invalid values are ignored.
29 % Leave empty to use default values.
30 % NAME is an optional label that is added to the plot titles.
31 % VERBOSE sets the level of status messages:
32 % 0 = silent & no data plots;
33 % 1 = status messages & minimum plots;
34 % 2 = all messages and plots (default)
35 %
36 % Outputs:
37 % RETVAL is the array of Allan deviation values at each TAU.
38 % S is an optional output of other statistical measures of the data (mean, std, etc).
39 % ERRORB is an optional output containing the error estimates for a 1-sigma
40 % confidence interval. These values are shown on the figure for each point.
41 % TAU is an optional output containing the array of tau values used in the
42 % calculation (which may be a truncated subset of the input or default values).
43 %
44 % Example:
45 %
46 % To compute the Allan deviation for the data in the variable "lt":
47 % >> lt
48 % lt =
49 % freq: [1x86400 double]
50 % rate: 0.5
51 %
52 % Use:
53 %
54 % >> ad = allan(lt,[2 10 100],'lt data',1);
55 %
56 % The Allan deviation will be computed and plotted at tau = 2,10,100 seconds.
57 % 1-sigma confidence intervals will be indicated by vertical lines at each point.
58 % You can also use the default settings, which are usually a good starting point:
59 %
60 % >> ad = allan(lt);
61 %
62 %
63 % Notes:
64 % This function calculates the standard Allan deviation (ADEV), *not* the
65 % overlapping ADEV. Use "allan_overlap.m" for overlapping ADEV.
66 % The calculation is performed using fractional frequency data. If only
67 % phase data is provided, frequency data is generated by differentiating
68 % the phase data.
69 % No pre-processing of the data is performed, except to remove any
70 % initial offset (i.e., starting gap) in the time record.
71 % For rate-based data, ADEV is computed only for tau values greater than the
72 % minimum time between samples and less than the half the total time. For
73 % time-stamped data, only tau values greater than the maximum gap between
74 % samples and less than half the total time are used.
75 % The calculation for fixed sample rate data is *much* faster than for
76 % time-stamp data. You may wish to run the rate-based calculation first,
77 % then compare with time-stamp-based. Often the differences are insignificant.
78 % To show the "tau bins" (y_k samples) on the data plot, set the variable
79 % TAUBIN to 1 (search for "#TAUBIN").
80 % You can choose between loglog and semilog plotting of results by
81 % commenting in/out the appropriate line. Search for "#PLOTLOG".
82 % I recommend installing "dsplot.m", which improves the performance of
83 % plotting large data sets. Download from File Exchange, File ID: #15850.
84 % allan.m will use dsplot.m if it is present on your MATLAB path.
85 % This function has been validated using the test data from NBS Monograph
86 % 140, the 1000-point test data set given by Riley [1], and the example data
87 % given in IEEE standard 1139-1999, Annex C.
88 % The author welcomes other validation results, see contact info below.
89 %
90 % For more information, see:
91 % [1] W. J. Riley, "The Calculation of Time Domain Frequency Stability,"
92 % Available on the web:
93 % http://www.ieee-uffc.org/frequency_control/teaching.asp?name=paper1ht
94 %
95 %
96 % M.A. Hopcroft
97 % mhopeng at gmail dot com
98 %
99 % I welcome your comments and feedback!
100 %
101 % MH Mar2014
102 % v2.24 fix bug related to generating freq data from phase with timestamps
103 % (thanks to S. David-Grignot for finding the bug)
104 % MH Oct2010
105 % v2.22 tau truncation to integer groups; tau sort
106 % plotting bugfix
107 % v2.20 sychronize updates across allan, allan_overlap, allan_modified
108 % v2.16 add TAU as output, fixed unusual error with dsplot v1.1
109 % v2.14 update plotting behaviour, default tau values
110 %
111
112 versionstr = 'allan v2.24';
113
114 % MH Jun2010
115 % v2.12 bugfix for rate data row/col orientation
116 % add DATA.units for plotting
117 % use dsplot.m for plotting
118 %
119 % MH MAR2010
120 % v2.1 minor interface and bugfixes
121 % update data consistency check
122 %
123 % MH FEB2010
124 % v2.0 Consistent code behaviour for all "allan_x.m" functions:
125 % accept phase data
126 % verbose levels
127 %
128 %
129 % MH JAN2010
130 % v1.84 code cleanup
131 % v1.82 typos in comments and code cleanup
132 % tau bin plotting changed for performance improvement
133 % v1.8 Performance improvements:
134 % vectorize code for rate data
135 % logical indexing for irregular rate data
136 % MH APR2008
137 % v1.62 loglog plot option
138 % v1.61 improve error handling, plotting
139 % fix bug in regular data calc for high-rate data
140 % fix bug in timestamp data calc for large starting gap
141 % (thanks to C. B. Ruiz for identifying these bugs)
142 % uses timestamps for DATA.rate=0
143 % progress indicator for large timestamp data processing
144 % MH JUN2007
145 % v1.54 Improve data plotting and optional bin plotting
146 % MH FEB2007
147 % v1.5 use difference from median for plotting
148 % added MAD calculation for outlier detection
149 % MH JAN2007
150 % v1.48 plotting typos fixes
151 % MH DEC2006
152 % v1.46 hack to plot error bars
153 % v1.44 further validation (Riley 1000-pt)
154 % plot mean and std
155 % MH NOV2006
156 % v1.42 typo fix comments
157 % v1.4 fix irregular rate algorithm
158 % irregular algorithm rejects tau less than max gap in time data
159 % validate both algorithms using test data from NBS Monograph 140
160 % v1.3 fix time calc if data.time not present
161 % add error bars (not possible due to bug in MATLAB R14SP3)
162 % remove offset calculation
163 % v1.24 improve feedback
164 % MH SEP2006
165 % v1.22 updated comments
166 % v1.2 errors and warnings
167 % v1.1 handle irregular interval data
168 %#ok<*AGROW>
169
170 % defaults
171 if nargin < 4, verbose=2; end
172 if nargin < 3, name=''; end
173 if nargin < 2 || isempty(tau), tau=2.^(-10:10); end
174
175 % plot "tau bins"? #TAUBIN
176 TAUBIN = 0; % set 0 or 1 % WARNING: this has a significant impact on performance
177
178 % Formatting for plots
179 FontName = 'Arial';
180 FontSize = 14;
181 plotlinewidth=2;
182
183 if verbose >= 1, fprintf(1,'allan: %s\n\n',versionstr); end
184
185 %% Data consistency checks
186 if ~(isfield(data,'phase') || isfield(data,'freq'))
187 error('Either ''phase'' or ''freq'' must be present in DATA. See help file for details. [con0]');
188 end
189 if isfield(data,'time')
190 if isfield(data,'phase') && (length(data.phase) ~= length(data.time))
191 if isfield(data,'freq') && (length(data.freq) ~= length(data.time))
192 error('The time and freq vectors are not the same length. See help for details. [con2]');
193 else
194 error('The time and phase vectors are not the same length. See help for details. [con1]');
195 end
196 end
197 if isfield(data,'phase') && (any(isnan(data.phase)) || any(isinf(data.phase)))
198 error('The phase vector contains invalid elements (NaN/Inf). [con3]');
199 end
200 if isfield(data,'freq') && (any(isnan(data.freq)) || any(isinf(data.freq)))
201 error('The freq vector contains invalid elements (NaN/Inf). [con4]');
202 end
203 if isfield(data,'time') && (any(isnan(data.time)) || any(isinf(data.time)))
204 error('The time vector contains invalid elements (NaN/Inf). [con5]');
205 end
206 end
207
208 % sort tau vector
209 tau=sort(tau);
210
211
212 %% Basic statistical tests on the data set
213 if ~isfield(data,'freq')
214 if isfield(data,'rate') && data.rate ~= 0
215 data.freq=diff(data.phase).*data.rate;
216 elseif isfield(data,'time')
217 data.freq=diff(data.phase)./diff(data.time);
218 end
219 if verbose >= 1, fprintf(1,'allan: Fractional frequency data generated from phase data (M=%g).\n',length(data.freq)); end
220 data.time(1)=[]; % make time stamps correspond to freq data
221 end
222 if size(data.freq,2) > size(data.freq,1), data.freq=data.freq'; end % ensure columns
223
224 s.numpoints=length(data.freq);
225 s.max=max(data.freq);
226 s.min=min(data.freq);
227 s.mean=mean(data.freq);
228 s.median=median(data.freq);
229 if isfield(data,'time')
230 if size(data.time,2) > size(data.time,1), data.time=data.time'; end % ensure columns
231 s.linear=polyfit(data.time(1:length(data.freq)),data.freq,1);
232 elseif isfield(data,'rate') && data.rate ~= 0;
233 s.linear=polyfit((1/data.rate:1/data.rate:length(data.freq)/data.rate)',data.freq,1);
234 else
235 error('Either "time" or "rate" must be present in DATA. Type "help allan" for details. [err1]');
236 end
237 s.std=std(data.freq);
238
239 if verbose >= 2
240 fprintf(1,'allan: input data statistics:\n');
241 disp(s);
242 end
243
244
245 % center at median for plotting
246 medianfreq=data.freq-s.median;
247 sm=[]; sme=[];
248
249 % Screen for outliers using 5x Median Absolute Deviation (MAD) criteria
250 s.MAD = median(abs(medianfreq)/0.6745);
251 if verbose >= 2
252 fprintf(1, 'allan: 5x MAD value for outlier detection: %g\n',5*s.MAD);
253 end
254 if verbose >= 1 && any(abs(medianfreq) > 5*s.MAD)
255 fprintf(1, 'allan: NOTE: There appear to be outliers in the frequency data. See plot.\n');
256 end
257
258
259 %%%%
260 % There are two cases, either using timestamps or fixed sample rate:
261
262 %% Fixed Sample Rate Data
263 % If there is a regular interval between measurements, calculation is much
264 % easier/faster
265 if isfield(data,'rate') && data.rate > 0 % if data rate was given
266 if verbose >= 1, fprintf(1, 'allan: regular data (%g data points @ %g Hz)\n',length(data.freq),data.rate); end
267
268 % string for plot title
269 name=[name ' (' num2str(data.rate) ' Hz)'];
270
271 % what is the time interval between data points?
272 tmstep = 1/data.rate;
273
274 % Is there time data? Just for curiosity/plotting, does not impact calculation
275 if isfield(data,'time')
276 % adjust time data to remove any starting gap; first time step
277 % should not be zero for comparison with freq data
278 dtime=data.time-data.time(1)+mean(diff(data.time));
279 if verbose >= 2
280 fprintf(1,'allan: End of timestamp data: %g sec.\n',dtime(end));
281 if (data.rate - 1/mean(diff(dtime))) > 1e-6
282 fprintf(1,'allan: NOTE: data.rate (%f Hz) does not match average timestamped sample rate (%f Hz)\n',data.rate,1/mean(diff(dtime)));
283 end
284 end
285 else
286 % create time axis data using rate (for plotting only)
287 dtime=(tmstep:tmstep:length(data.freq)*tmstep)'; % column oriented
288 end
289
290 % check the range of tau values and truncate if necessary
291 % find halfway point of time record
292 halftime = round(tmstep*length(data.freq)/2);
293 % truncate tau to appropriate values
294 tau = tau(tau >= tmstep & tau <= halftime);
295 if verbose >= 2, fprintf(1, 'allan: allowable tau range: %g to %g sec. (1/rate to total_time/2)\n',tmstep,halftime); end
296
297 % save the freq data for the loop
298 dfreq=data.freq;
299 % find the number of data points in each tau group
300 m = data.rate.*tau;
301 % only integer values allowed (no fractional groups of points)
302 %tau = tau(m-round(m)<1e-8); % numerical precision issues (v2.1)
303 tau = tau(m==round(m)); % The round() test is only correct for values < 2^53
304 %m = m(m-round(m)<1e-8); % change to round(m) for integer test v2.22
305 m = m(m==round(m));
306 %m=round(m);
307
308 if verbose >= 1, fprintf(1,'allan: calculating Allan deviation...\n '); end
309
310 % calculate the Allan deviation for each value of tau
311 k=0; tic;
312 for i = tau
313 if verbose >= 2, fprintf(1,'%g ',i); end
314 k=k+1;
315
316 % truncate frequency set to an even multiple of this tau value
317 freq=dfreq(1:end-rem(length(dfreq),m(k)));
318 % group the data into tau-length groups or bins
319 f = reshape(freq,m(k),[]); % Vectorize!
320 % find average in each "tau group", y_k (each colummn of f)
321 fa=mean(f,1);
322 % first finite difference
323 fd=diff(fa);
324 % calculate two-sample variance for this tau
325 M=length(fa);
326 sm(k)=sqrt(0.5/(M-1)*(sum(fd.^2)));
327
328 % estimate error bars
329 sme(k)=sm(k)/sqrt(M+1);
330
331 if TAUBIN == 1
332 % save the binning points for plotting
333 fs(k,1:length(freq)/m(k))=m(k):m(k):length(freq); fval{k}=mean(f,1);
334 end
335
336 end % repeat for each value of tau
337
338 if verbose >= 2, fprintf(1,'\n'); end
339 calctime=toc; if verbose >= 2, fprintf(1,'allan: Elapsed time for calculation: %e seconds\n',calctime); end
340
341
342
343 %% Irregular data (timestamp)
344 elseif isfield(data,'time')
345 % the interval between measurements is irregular
346 % so we must group the data by time
347 if verbose >= 1, fprintf(1, 'allan: irregular rate data (no fixed sample rate)\n'); end
348
349 % string for plot title
350 name=[name ' (timestamp)'];
351
352 % adjust time to remove any initial offset or zero
353 dtime=data.time-data.time(1)+mean(diff(data.time));
354 %dtime=data.time;
355 % where is the maximum gap in time record?
356 gap_pos=find(diff(dtime)==max(diff(dtime)));
357 % what is average data spacing?
358 avg_gap = mean(diff(dtime));
359
360 if verbose >= 2
361 fprintf(1, 'allan: WARNING: irregular timestamp data (no fixed sample rate).\n');
362 fprintf(1, ' Calculation time may be long and the results subject to interpretation.\n');
363 fprintf(1, ' You are advised to estimate using an average sample rate (%g Hz) instead of timestamps.\n',1/avg_gap);
364 fprintf(1, ' Continue at your own risk! (press any key to continue)\n');
365 pause;
366 end
367
368 if verbose >= 1
369 fprintf(1, 'allan: End of timestamp data: %g sec\n',dtime(end));
370 fprintf(1, ' Average rate: %g Hz (%g sec/measurement)\n',1/avg_gap,avg_gap);
371 if max(diff(dtime)) ~= 1/mean(diff(dtime))
372 fprintf(1, ' Max. gap: %g sec at position %d\n',max(diff(dtime)),gap_pos(1));
373 end
374 if max(diff(dtime)) > 5*avg_gap
375 fprintf(1, ' WARNING: Max. gap in time record is suspiciously large (>5x the average interval).\n');
376 end
377 end
378
379
380 % find halfway point
381 halftime = fix(dtime(end)/2);
382 % truncate tau to appropriate values
383 tau = tau(tau >= max(diff(dtime)) & tau <= halftime);
384 if isempty(tau)
385 error('allan: ERROR: no appropriate tau values (> %g s, < %g s)\n',max(diff(dtime)),halftime);
386 end
387
388 % save the freq data for the loop
389 dfreq=data.freq;
390 dtime=dtime(1:length(dfreq));
391
392 if verbose >= 1, fprintf(1,'allan: calculating Allan deviation...\n'); end
393
394 k=0; tic;
395 for i = tau
396 if verbose >= 2, fprintf(1,'%d ',i); end
397
398 k=k+1; fa=[]; %f=[];
399 km=0;
400
401 % truncate data set to an even multiple of this tau value
402 freq=dfreq(dtime <= dtime(end)-rem(dtime(end),i));
403 time=dtime(dtime <= dtime(end)-rem(dtime(end),i));
404 %freq=dfreq;
405 %time=dtime;
406
407 % break up the data into groups of tau length in sec
408 while i*km < time(end)
409 km=km+1;
410
411 % progress bar
412 if verbose >= 2
413 if rem(km,100)==0, fprintf(1,'.'); end
414 if rem(km,1000)==0, fprintf(1,'%g/%g\n',km,round(time(end)/i)); end
415 end
416
417 f = freq(i*(km-1) < time & time <= i*km);
418 f = f(~isnan(f)); % make sure values are valid
419
420 if ~isempty(f)
421 fa(km)=mean(f);
422 else
423 fa(km)=0;
424 end
425
426 if TAUBIN == 1 % WARNING: this has a significant impact on performance
427 % save the binning points for plotting
428 %if find(time <= i*km) > 0
429 fs(k,km)=max(time(time <= i*km));
430 %else
431 if isempty(fs(k,km))
432 fs(k,km)=0;
433 end
434 fval{k}=fa;
435 end % save tau bin plot points
436
437 end
438
439 if verbose >= 2, fprintf(1,'\n'); end
440
441 % first finite difference of the averaged results
442 fd=diff(fa);
443 % calculate Allan deviation for this tau
444 M=length(fa);
445 sm(k)=sqrt(0.5/(M-1)*(sum(fd.^2)));
446
447 % estimate error bars
448 sme(k)=sm(k)/sqrt(M+1);
449
450
451 end
452
453 if verbose == 2, fprintf(1,'\n'); end
454 calctime=toc; if verbose >= 2, fprintf(1,'allan: Elapsed time for calculation: %e seconds\n',calctime); end
455
456
457 else
458 error('allan: WARNING: no DATA.rate or DATA.time! Type "help allan" for more information. [err2]');
459 end
460
461
462 %%%%%%%%
463 %% Plotting
464
465 if verbose >= 2 % show all data
466
467 % plot the frequency data, centered on median
468 if size(dtime,2) > size(dtime,1), dtime=dtime'; end % this should not be necessary, but dsplot 1.1 is a little bit brittle
469 try
470 % dsplot makes a new figure
471 hd=dsplot(dtime,medianfreq);
472 catch ME
473 figure;
474 if length(dtime) ~= length(medianfreq)
475 fprintf(1,'allan: Warning: length of time axis (%d) is not equal to data array (%d)\n',length(dtime),length(medianfreq));
476 end
477 hd=plot(dtime,medianfreq);
478 if verbose >= 1, fprintf(1,'allan: Note: Install dsplot.m for improved plotting of large data sets (File Exchange File ID: #15850).\n'); end
479 if verbose >= 2, fprintf(1,' (Message: %s)\n',ME.message); end
480 end
481 set(hd,'Marker','.','LineStyle','none','Color','b'); % equivalent to '.-'
482 hold on;
483
484 % show center (0)
485 plot(xlim,[0 0],':k');
486 % show 5x Median Absolute Deviation (MAD) values
487 hm=plot(xlim,[5*s.MAD 5*s.MAD],'-r');
488 plot(xlim,[-5*s.MAD -5*s.MAD],'-r');
489 % show linear fit line
490 hf=plot(xlim,polyval(s.linear,xlim)-s.median,'-g');
491 title(['Data: ' name],'FontSize',FontSize+2,'FontName',FontName);
492 %set(get(gca,'Title'),'Interpreter','none');
493 xlabel('Time [sec]','FontSize',FontSize,'FontName',FontName);
494 if isfield(data,'units')
495 ylabel(['data - median(data) [' data.units ']'],'FontSize',FontSize,'FontName',FontName);
496 else
497 ylabel('freq - median(freq)','FontSize',FontSize,'FontName',FontName);
498 end
499 set(gca,'FontSize',FontSize,'FontName',FontName);
File was created 1 function [retval, s, errorb, tau] = allan_modified(data,tau,name,verbose)
2 % ALLAN_MODIFIED Compute the modified Allan deviation for a set of
3 % time-domain frequency data
4 % [RETVAL, S, ERRORB, TAU] = ALLAN_MODIFIED(DATA,TAU,NAME,VERBOSE)
5 %
6 % Inputs:
7 % DATA should be a struct and have the following fields:
8 % DATA.freq or DATA.phase
9 % A vector of fractional frequency measurements (df/f) in
10 % DATA.freq *or* phase offset data (seconds) in DATA.phase
11 % If phase data is not present, it will be generated by
12 % integrating the fractional frequency data.
13 % If both fields are present, then DATA.phase will be used.
14 %
15 % DATA.rate or DATA.time
16 % The sampling rate in Hertz (DATA.rate) or a vector of
17 % timestamps for each measurement in seconds (DATA.time).
18 % DATA.rate is used if both fields are present.
19 % If DATA.rate == 0, then the timestamps are used.
20 %
21 % TAU is an array of tau values for computing Allan deviation.
22 % TAU values must be divisible by 1/DATA.rate (data points cannot be
23 % grouped in fractional quantities!). Invalid values are ignored.
24 % NAME is an optional label that is added to the plot titles.
25 % VERBOSE sets the level of status messages:
26 % 0 = silent & no data plots; 1 = status messages; 2 = all messages
27 %
28 % Outputs:
29 % RETVAL is the array of modified Allan deviation values at each TAU.
30 % S is an optional output of other statistical measures of the data (mean, std, etc).
31 % ERRORB is an optional output containing the error estimates for a
32 % 1-sigma confidence interval. These values are shown on the figure for each point.
33 % TAU is an optional output containing the array of tau values used in the
34 % calculation (which may be a truncated subset of the input or default values).
35 %
36 % Example:
37 %
38 % To compute the modified Allan deviation for the data in the variable "lt":
39 % >> lt
40 % lt =
41 % freq: [1x86400 double]
42 % rate: 0.5
43 %
44 % Use:
45 %
46 % >> adm = allan_modified(lt,[2 10 100],'lt data',1);
47 %
48 % The modified Allan deviation will be computed and plotted at tau = 2,10,100 seconds.
49 % 1-sigma confidence intervals will be indicated by vertical lines at each point.
50 % You can also use the default settings, which are usually a good starting point:
51 %
52 % >> adm = allan_modified(lt);
53 %
54 %
55 % Notes:
56 % This function calculates the modifed Allan deviation (MDEV).
57 % The calculation is performed using phase data. If only frequency data is
58 % provided, phase data is generated by integrating the frequency data.
59 % No pre-processing of the data is performed.
60 % For rate-based data, MDEV is computed only for tau values greater than the
61 % minimum time between samples and less than the half the total time. For
62 % time-stamped data, only tau values greater than the maximum gap between
63 % samples and less than half the total time are used.
64 % The calculation for fixed sample rate data is *much* faster than for
65 % time-stamp data. You may wish to run the rate-based calculation first,
66 % then compare with time-stamp-based. Often the differences are insignificant.
67 % When phase data is generated by integrating time-stamped frequency data,
68 % the final data point is dropped, because there is no timestamp for it.
69 % This will create a [usually small] difference between the results from
70 % analyzing the same data set with timestamp data and analyzing with a
71 % fixed sample rate. See note in the code near line 350.
72 % You can choose between loglog and semilog plotting of results by
73 % commenting in/out the appropriate line. Search for "#PLOTLOG".
74 % This function has been validated using the test data from NBS Monograph
75 % 140, the 1000-point test data set given by Riley [1], and the example data
76 % given in IEEE standard 1139-1999, Annex C.
77 % The author welcomes other validation results, see contact info below.
78 %
79 % For more information, see:
80 % [1] W. J. Riley, "The Calculation of Time Domain Frequency Stability,"
81 % Available on the web:
82 % http://www.ieee-uffc.org/frequency_control/teaching.asp?name=paper1ht
83 %
84 %
85 % M.A. Hopcroft
86 % mhopeng at gmail dot com
87 %
88 % I welcome your comments and feedback!
89 %
90 % MH Mar2014
91 % v1.24 fix bug related to generating freq data from phase with timestamps
92 % (thanks to S. David-Grignot for finding the bug)
93 % MH Oct2010
94 % v1.22 tau truncation to integer groups; tau sort
95 % plotting bugfix
96 % v1.20 update to match allan.m (dsplot.m, columns)
97 % discard tau values with timestamp irregularities
98 %
99
100 versionstr = 'allan_modified v1.24';
101
102 % MH MAR2010
103 % v1.1 bugfixes for irregular sample rates
104 % update consistency check
105 %
106 % MH FEB2010
107 % v1.0 based on allan_overlap v2.0
108 %
109
110 %#ok<*AGROW>
111
112
113 % defaults
114 if nargin < 4, verbose = 2; end
115 if nargin < 3, name=''; end
116 if nargin < 2 || isempty(tau), tau=2.^(-10:10); end
117 if isfield(data,'rate') && isempty(data.rate), data.rate=0; end % v1.1
118
119 % Formatting for plots
120 FontName = 'Arial';
121 FontSize = 14;
122 plotlinewidth=2;
123
124 if verbose >= 1, fprintf(1,'allan_modified: %s\n\n',versionstr); end
125
126 %% Data consistency checks
127 if ~(isfield(data,'phase') || isfield(data,'freq'))
128 error('Either ''phase'' or ''freq'' must be present in DATA. See help file for details. [con0]');
129 end
130 if isfield(data,'time')
131 if isfield(data,'phase') && (length(data.phase) ~= length(data.time))
132 if isfield(data,'freq') && (length(data.freq) ~= length(data.time))
133 error('The time and freq vectors are not the same length. See help for details. [con2]');
134 else
135 error('The time and phase vectors are not the same length. See help for details. [con1]');
136 end
137 end
138 if isfield(data,'phase') && (any(isnan(data.phase)) || any(isinf(data.phase)))
139 error('The phase vector contains invalid elements (NaN/Inf). [con3]');
140 end
141 if isfield(data,'freq') && (any(isnan(data.freq)) || any(isinf(data.freq)))
142 error('The freq vector contains invalid elements (NaN/Inf). [con4]');
143 end
144 if isfield(data,'time') && (any(isnan(data.time)) || any(isinf(data.time)))
145 error('The time vector contains invalid elements (NaN/Inf). [con5]');
146 end
147 end
148
149 % sort tau vector
150 tau=sort(tau);
151
152
153 %% Basic statistical tests on the data set
154 if ~isfield(data,'freq')
155 if isfield(data,'rate') && data.rate ~= 0
156 data.freq=diff(data.phase).*data.rate;
157 elseif isfield(data,'time')
158 data.freq=diff(data.phase)./diff(data.time);
159 end
160 if verbose >= 1, fprintf(1,'allan_modified: Fractional frequency data generated from phase data (M=%g).\n',length(data.freq)); end
161 end
162 if size(data.freq,2) > size(data.freq,1), data.freq=data.freq'; end % ensure columns
163
164 s.numpoints=length(data.freq);
165 s.max=max(data.freq);
166 s.min=min(data.freq);
167 s.mean=mean(data.freq);
168 s.median=median(data.freq);
169 if isfield(data,'time')
170 if size(data.time,2) > size(data.time,1), data.time=data.time'; end % ensure columns
171 s.linear=polyfit(data.time(1:length(data.freq)),data.freq,1);
172 elseif isfield(data,'rate') && data.rate ~= 0;
173 s.linear=polyfit((1/data.rate:1/data.rate:length(data.freq)/data.rate)',data.freq,1);
174 else
175 error('Either "time" or "rate" must be present in DATA. Type "help allan_modified" for details. [err1]');
176 end
177 s.std=std(data.freq);
178
179 if verbose >= 2
180 fprintf(1,'allan_modified: fractional frequency data statistics:\n');
181 disp(s);
182 end
183
184 % scale to median for plotting
185 medianfreq=data.freq-s.median;
186 sm=[]; sme=[];
187
188 % Screen for outliers using 5x Median Absolute Deviation (MAD) criteria
189 MAD = median(abs(medianfreq)/0.6745);
190 if verbose >= 1 && any(abs(medianfreq) > 5*MAD)
191 fprintf(1, 'allan_modified: NOTE: There appear to be outliers in the frequency data. See plot.\n');
192 end
193
194 %%%%
195 % There are two cases, either using timestamps or rate:
196
197 %% Fixed Sample Rate Data
198 % If there is a regular interval between measurements, calculation is much
199 % easier/faster
200 if isfield(data,'rate') && data.rate > 0 % if data rate was given
201 if verbose >= 1
202 fprintf(1, 'allan_modified: regular data ');
203 if isfield(data,'freq')
204 fprintf(1, '(%g freq data points @ %g Hz)\n',length(data.freq),data.rate);
205 elseif isfield(data,'phase')
206 fprintf(1, '(%g phase data points @ %g Hz)\n',length(data.phase),data.rate);
207 else
208 error('\n phase or freq data missing [err10]');
209 end
210 end
211
212 % string for plot title
213 name=[name ' (' num2str(data.rate) ' Hz)'];
214
215 % what is the time interval between data points?
216 tmstep = 1/data.rate;
217
218 % Is there time data? Just for curiosity/plotting, does not impact calculation
219 if isfield(data,'time')
220 % adjust time data to remove any starting gap; first time step
221 % should not be zero for comparison with freq data
222 dtime=data.time-data.time(1)+mean(diff(data.time));
223 dtime=dtime(1:length(medianfreq)); % equalize the data vector lengths for plotting (v1.1)
224 if verbose >= 2
225 fprintf(1,'allan_modified: End of timestamp data: %g sec.\n',dtime(end));
226 if (data.rate - 1/mean(diff(dtime))) > 1e-6
227 fprintf(1,'allan_modified: NOTE: data.rate (%f Hz) does not match average timestamped sample rate (%f Hz)\n',data.rate,1/mean(diff(dtime)));
228 end
229 end
230 else
231 % create time axis data using rate (for plotting only)
232 dtime=(tmstep:tmstep:length(data.freq)*tmstep);
233 end
234
235
236 % is phase data present? If not, generate it
237 if ~isfield(data,'phase')
238 nfreq=data.freq-s.mean;
239 dphase=zeros(1,length(nfreq)+1);
240 dphase(2:end) = cumsum(nfreq).*tmstep;
241 if verbose >= 1, fprintf(1,'allan_modified: phase data generated from fractional frequency data (N=%g).\n',length(dphase)); end
242 else
243 dphase=data.phase;
244 end
245
246
247 % check the range of tau values and truncate if necessary
248 % find halfway point of time record
249 halftime = round(tmstep*length(data.freq)/2);
250 % truncate tau to appropriate values
251 tau = tau(tau >= tmstep & tau <= halftime);
252 if verbose >= 2, fprintf(1, 'allan_modified: allowable tau range: %g to %g sec. (1/rate to total_time/2)\n',tmstep,halftime); end
253
254 % find the number of data points in each tau group
255 % number of samples
256 N=length(dphase);
257 m = data.rate.*tau;
258 % only integer values allowed (no fractional groups of points)
259 %tau = tau(m-round(m)<1e-8); % numerical precision issues (v1.1)
260 tau = tau(m==round(m)); % The round() test is only correct for values < 2^53
261 %m = m(m-round(m)<1e-8); % change to round(m) for integer test v1.22
262 m = m(m==round(m));
263 %m=round(m);
264
265 if verbose >= 1, fprintf(1,'allan_modified: calculating modified Allan deviation...\n '); end
266
267
268 % calculate the modified Allan deviation for each value of tau
269 k=0; tic;
270 for i = tau
271 k=k+1;
272 pa=[];
273 if verbose >= 2, fprintf(1,'%d ',i); end
274
275 mphase = dphase;
276
277 % calculate overlapping "phase averages" (x_k)
278 for p=1:m(k)
279
280 % truncate frequency set length to an even multiple of this tau value
281 mphase=mphase(1:end-rem(length(mphase),m(k)));
282 % group phase values
283 mp=reshape(mphase,m(k),[]);
284 % find average in each "tau group" (each column of mp)
285 pa(p,:)=mean(mp,1);
286 % shift data vector by -1 and repeat
287 mphase=circshift(dphase,(size(dphase)>1)*-p);
288
289 end
290
291 % create "modified" y_k freq values
292 mfreq=diff(pa,1,2)./i;
293 mfreq=reshape(mfreq,1,[]);
294
295 % calculate modified frequency differences
296 mfreqd=reshape(mfreq,m(k),[]); % Vectorize!
297 mfreqd=diff(mfreqd,1,2);
298 mfreqd=reshape(mfreqd,1,[]);
299
300
301 % calculate two-sample variance for this tau
302 sm(k)=sqrt((1/(2*(N-3*m(k)+1)))*(sum(mfreqd(1:N-3*m(k)+1).^2)));
303
304 % estimate error bars
305 sme(k)=sm(k)/sqrt(N-3*m(k)+1);
306
307
308 end % repeat for each value of tau
309
310 if verbose >= 2, fprintf(1,'\n'); end
311 calctime=toc; if verbose >= 2, fprintf(1,'allan_modified: Elapsed time for calculation: %g seconds\n',calctime); end
312
313
314
315 %% Irregular data (timestamp)
316 elseif isfield(data,'time')
317 % the interval between measurements is irregular
318 % so we must group the data by time
319 if verbose >= 1, fprintf(1, 'allan_modified: irregular rate data (no fixed sample rate)\n'); end
320
321 % string for plot title
322 name=[name ' (timestamp)'];
323
324 % adjust time to remove any initial offset
325 dtime=data.time-data.time(1)+mean(diff(data.time));
326 %dtime=data.time-data.time(1);
327 % where is the maximum gap in time record?
328 gap_pos=find(diff(dtime)==max(diff(dtime)));
329 % what is average data spacing?
330 avg_gap = mean(diff(dtime));
331
332 if verbose >= 2
333 fprintf(1, 'allan_modified: WARNING: irregular timestamp data (no fixed sample rate).\n');
334 fprintf(1, ' Calculation time may be long and the results subject to interpretation.\n');
335 fprintf(1, ' You are advised to estimate using an average sample rate (%g Hz) instead of timestamps.\n',1/avg_gap);
336 fprintf(1, ' Continue at your own risk! (press any key to continue)\n');
337 pause;
338 end
339
340 if verbose >= 1
341 fprintf(1, 'allan_modified: End of timestamp data: %g sec\n',dtime(end));
342 fprintf(1, ' Average sample rate: %g Hz (%g sec/measurement)\n',1/avg_gap,avg_gap);
343 if max(diff(dtime)) ~= 1/mean(diff(dtime))
344 fprintf(1, ' Max. gap in time record: %g sec at position %d\n',max(diff(dtime)),gap_pos(1));
345 end
346 if max(diff(dtime)) > 5*avg_gap
347 fprintf(1, ' WARNING: Max. gap in time record is suspiciously large (>5x the average interval).\n');
348 end
349 end
350
351 % is phase data present? If not, generate it
352 if ~isfield(data,'phase')
353 nfreq=data.freq-s.mean;
354 % NOTE: uncommenting the following two lines will artificially
355 % allow the code to give equivalent results for validation data
356 % sets using fixed rate data and timestamped data by adding a
357 % "phantom" data point for frequency integration. Use of this
358 % phantom point can skew the results of calculations on real data.
359 %nfreq(end+1)=0; % phantom freq point, with average value
360 %dtime(end+1)=dtime(end)+avg_gap; % phantom average time step
361 dphase=zeros(1,length(nfreq));
362 dphase(2:end) = cumsum(nfreq(1:end-1)).*diff(dtime);
363 if verbose >= 1, fprintf(1,'allan_modified: Phase data generated from fractional frequency data (N=%g).\n',length(dphase)); end
364 else
365 dphase=data.phase;
366 end
367
368 % find halfway point
369 halftime = fix(dtime(end)/2);
370 % truncate tau to appropriate values
371 tau = tau(tau >= max(diff(dtime)) & tau <= halftime);
372 if isempty(tau)
373 error('allan_modified: ERROR: no appropriate tau values (> %g s, < %g s)\n',max(diff(dtime)),halftime);
374 end
375
376 % % save the freq data for the loop
377 % dfreq=data.freq;
378
379 % number of samples
380 N=length(dphase);
381 m=round(tau./mean(diff(dtime)));
382
383 if verbose >= 1, fprintf(1,'allan_modified: calculating modified Allan deviation...\n'); end
384
385 k=0; tic;
386 for i = tau
387
388 k=k+1; pa=[];
389
390 mphase = dphase; time = dtime;
391
392 if verbose >= 2, fprintf(1,'%d ',i); end
393
394 % calculate overlapping "phase averages" (x_k)
395 %for j = 1:i
396 for j = 1:m(k) % (v1.1)
397 km=0;
398 %fprintf(1,'j: %d ',j);
399
400 % (v1.1) truncating not correct for overlapping samples
401 % truncate data set to an even multiple of this tau value
402 %mphase = mphase(time <= time(end)-rem(time(end),i));
403 %time = time(time <= time(end)-rem(time(end),i));
404
405
406 % break up the data into overlapping groups of tau length
407 while i*km < time(end)
408 km=km+1;
409
410 % progress bar
411 if verbose >= 2
412 if rem(km,100)==0, fprintf(1,'.'); end
413 if rem(km,1000)==0, fprintf(1,'%g/%g\n',km,round(time(end)/i)); end
414 end
415
416 mp = mphase(i*(km-1) < (time) & (time) <= i*km);
417
418 if ~isempty(mp)
419 pa(j,km)=mean(mp);
420 else
421 pa(j,km)=0;
422 end
423
424 end
425
426 % shift data vector by -1 and repeat
427 mphase=circshift(dphase,(size(mphase)>1)*-j);
428 mphase(end-j+1:end)=[];
429 time=circshift(dtime,(size(time)>1)*-j);
430 time(end-j+1:end)=[];
431 time=time-time(1)+avg_gap; % remove time offset
432
433 end
434
435 % create "modified" y_k freq values
436 mfreq=diff(pa,1,2)./i;
437 mfreq=reshape(mfreq,1,[]);
438
439 % calculate modified frequency differences
440 mfreqd=reshape(mfreq,m(k),[]); % Vectorize!
441 mfreqd=diff(mfreqd,1,2);
442 mfreqd=reshape(mfreqd,1,[]);
443
444 % calculate two-sample variance for this tau
445 % only the first N-3*m(k)+1 samples are valid
446 if length(mfreqd) >= N-3*m(k)+1
447 sm(k)=sqrt((1/(2*(N-3*m(k)+1)))*(sum(mfreqd(1:N-3*m(k)+1).^2)));
448
449 % estimate error bars
450 sme(k)=sm(k)/sqrt(N);
451
452 if verbose >= 2, fprintf(1,'\n'); end
453 else
454 if verbose >=2, fprintf(1,' tau=%g dropped due to timestamp irregularities\n',tau(k)); end
455 sm(k)=0; sme(k)=0;
456 end
457
458
459 end
460
461 if verbose >= 2, fprintf(1,'\n'); end
462 calctime=toc; if verbose >= 2, fprintf(1,'allan_modified: Elapsed time for calculation: %g seconds\n',calctime); end
463
464 % remove any points that were dropped
465 tau(sm==0)=[];
466 sm(sm==0)=[];
File was created 1 function [retval, s, errorb, tau] = allan_overlap(data,tau,name,verbose)
2 % ALLAN_OVERLAP Compute the overlapping Allan deviation for a set of
3 % time-domain frequency data
4 % [RETVAL, S, ERRORB, TAU] = ALLAN_OVERLAP(DATA,TAU,NAME,VERBOSE)
5 %
6 % Inputs:
7 % DATA should be a struct and have the following fields:
8 % DATA.freq or DATA.phase
9 % A vector of fractional frequency measurements (df/f) in
10 % DATA.freq *or* phase offset data (seconds) in DATA.phase
11 % If phase data is not present, it will be generated by
12 % integrating the fractional frequency data.
13 % If both fields are present, then DATA.phase will be used.
14 %
15 % DATA.rate or DATA.time
16 % The sampling rate in Hertz (DATA.rate) or a vector of
17 % timestamps for each measurement in seconds (DATA.time).
18 % DATA.rate is used if both fields are present.
19 % If DATA.rate == 0, then the timestamps are used.
20 %
21 % TAU is an array of tau values for computing Allan deviation.
22 % TAU values must be divisible by 1/DATA.rate (data points cannot be
23 % grouped in fractional quantities!). Invalid values are ignored.
24 % NAME is an optional label that is added to the plot titles.
25 % VERBOSE sets the level of status messages:
26 % 0 = silent & no data plots; 1 = status messages; 2 = all messages
27 %
28 % Outputs:
29 % RETVAL is the array of overlapping Allan deviation values at each TAU.
30 % S is an optional output of other statistical measures of the data (mean, std, etc).
31 % ERRORB is an optional output containing the error estimates for a 1-sigma
32 % confidence interval. Error bars are plotted as vertical lines at each point.
33 % TAU is an optional output containing the array of tau values used in the
34 % calculation (which may be a truncated subset of the input or default values).
35 %
36 % Example:
37 %
38 % To compute the overlapping Allan deviation for the data in the variable "lt":
39 % >> lt
40 % lt =
41 % freq: [1x86400 double]
42 % rate: 0.5
43 %
44 % Use:
45 %
46 % >> ado = allan_overlap(lt,[2 10 100],'lt data',1);
47 %
48 % The Allan deviation will be computed and plotted at tau = 2,10,100 seconds.
49 % 1-sigma confidence intervals will be indicated by vertical lines.
50 % You can also use the default settings, which are usually a good starting point:
51 %
52 % >> ado = allan_overlap(lt);
53 %
54 %
55 % Notes:
56 % This function calculates the overlapping Allan deviation (ADEV), *not* the
57 % standard ADEV. Use "allan.m" for standard ADEV.
58 % The calculation is performed using phase data. If only frequency data is
59 % provided, phase data is generated by integrating the frequency data.
60 % However, the timestamp-based calculation is performed using frequency
61 % data. Phase data is differentiated to generate frequency data if necessary.
62 % No pre-processing of the data is performed, except to remove any
63 % initial offset in the time record.
64 % For rate-based data, ADEV is computed only for tau values greater than the
65 % minimum time between samples and less than the half the total time. For
66 % time-stamped data, only tau values greater than the maximum gap between
67 % samples and less than half the total time are used.
68 % The calculation for fixed sample rate data is *much* faster than for
69 % time-stamp data. You may wish to run the rate-based calculation first,
70 % then compare with time-stamp-based. Often the differences are insignificant.
71 % The error bars at each point are calculated using the 1-sigma intervals
72 % based on the size of the data set. This is usually an overestimate for
73 % overlapping ADEV; a more accurate (and usually smaller uncertainty)
74 % value can be determined from chi-squared statistics, but that is not
75 % implemented in this version.
76 % You can choose between loglog and semilog plotting of results by
77 % commenting in/out the appropriate line. Search for "#PLOTLOG".
78 % This function has been validated using the test data from NBS Monograph
79 % 140, the 1000-point test data set given by Riley [1], and the example data
80 % given in IEEE standard 1139-1999, Annex C.
81 % The author welcomes other validation results, see contact info below.
82 %
83 % For more information, see:
84 % [1] W. J. Riley, "Addendum to a test suite for the calculation of time domain
85 % frequency stability," presented at IEEE Frequency Control Symposium,
86 % 1996.
87 % Available on the web:
88 % http://www.ieee-uffc.org/frequency_control/teaching.asp?name=paper1ht
89 %
90 %
91 % M.A. Hopcroft
92 % mhopeng at gmail dot com
93 %
94 % I welcome your comments and feedback!
95 %
96 % MH Mar2014
97 % v2.24 fix bug related to generating freq data from phase with timestamps
98 % (thanks to S. David-Grignot for finding the bug)
99 % MH Oct2010
100 % v2.22 tau truncation to integer groups; tau sort
101 % plotting bugfix
102 % v2.20 update to match allan.m (dsplot.m, columns)
103 % discard tau values with timestamp irregularities
104
105 versionstr = 'allan_overlap v2.24';
106
107 %
108 % MH MAR2010
109 % v2.1 bugfixes for irregular sample rates
110 % (thanks to Ryad Ben-El-Kezadri for feedback and testing)
111 % handle empty rate field
112 % fix integer comparisons for fractional sample rates
113 % update consistency check
114 %
115 % MH FEB2010
116 % v2.0 use phase data for calculation- much faster
117 % Consistent code behaviour for all "allan_x.m" functions:
118 % accept phase data
119 % verbose levels
120 %
121 % MH JAN2010
122 % v1.0 based on allan v1.84
123 %
124
125 %#ok<*AGROW>
126
127
128 % defaults
129 if nargin < 4, verbose = 2; end
130 if nargin < 3, name=''; end
131 if nargin < 2 || isempty(tau), tau=2.^(-10:10); end
132 if isfield(data,'rate') && isempty(data.rate), data.rate=0; end % v2.1
133
134 % Formatting for plots
135 FontName = 'Arial';
136 FontSize = 14;
137 plotlinewidth=2;
138
139 if verbose >= 1, fprintf(1,'allan_overlap: %s\n\n',versionstr); end
140
141 %% Data consistency checks v2.1
142 if ~(isfield(data,'phase') || isfield(data,'freq'))
143 error('Either ''phase'' or ''freq'' must be present in DATA. See help file for details. [con0]');
144 end
145 if isfield(data,'time')
146 if isfield(data,'phase') && (length(data.phase) ~= length(data.time))
147 if isfield(data,'freq') && (length(data.freq) ~= length(data.time))
148 error('The time and freq vectors are not the same length. See help for details. [con2]');
149 else
150 error('The time and phase vectors are not the same length. See help for details. [con1]');
151 end
152 end
153 if isfield(data,'phase') && (any(isnan(data.phase)) || any(isinf(data.phase)))
154 error('The phase vector contains invalid elements (NaN/Inf). [con3]');
155 end
156 if isfield(data,'freq') && (any(isnan(data.freq)) || any(isinf(data.freq)))
157 error('The freq vector contains invalid elements (NaN/Inf). [con4]');
158 end
159 if isfield(data,'time') && (any(isnan(data.time)) || any(isinf(data.time)))
160 error('The time vector contains invalid elements (NaN/Inf). [con5]');
161 end
162 end
163
164 % sort tau vector
165 tau=sort(tau);
166
167 %% Basic statistical tests on the data set
168 if ~isfield(data,'freq')
169 if isfield(data,'rate') && data.rate ~= 0
170 data.freq=diff(data.phase).*data.rate;
171 elseif isfield(data,'time')
172 data.freq=diff(data.phase)./diff(data.time);
173 end
174 if verbose >= 1, fprintf(1,'allan_overlap: Fractional frequency data generated from phase data (M=%g).\n',length(data.freq)); end
175 end
176 if size(data.freq,2) > size(data.freq,1), data.freq=data.freq'; end % ensure columns
177
178 s.numpoints=length(data.freq);
179 s.max=max(data.freq);
180 s.min=min(data.freq);
181 s.mean=mean(data.freq);
182 s.median=median(data.freq);
183 if isfield(data,'time')
184 if size(data.time,2) > size(data.time,1), data.time=data.time'; end % ensure columns
185 s.linear=polyfit(data.time(1:length(data.freq)),data.freq,1);
186 elseif isfield(data,'rate') && data.rate ~= 0;
187 s.linear=polyfit((1/data.rate:1/data.rate:length(data.freq)/data.rate)',data.freq,1);
188 else
189 error('Either "time" or "rate" must be present in DATA. Type "help allan_overlap" for details. [err1]');
190 end
191 s.std=std(data.freq);
192
193 if verbose >= 2
194 fprintf(1,'allan_overlap: fractional frequency data statistics:\n');
195 disp(s);
196 end
197
198
199 % scale to median for plotting
200 medianfreq=data.freq-s.median;
201 sm=[]; sme=[];
202
203 % Screen for outliers using 5x Median Absolute Deviation (MAD) criteria
204 MAD = median(abs(medianfreq)/0.6745);
205 if verbose >= 1 && any(abs(medianfreq) > 5*MAD)
206 fprintf(1, 'allan_overlap: NOTE: There appear to be outliers in the frequency data. See plot.\n');
207 end
208
209 %%%%
210 % There are four cases, freq or phase data, using timestamps or rate:
211
212 %% Fixed Sample Rate Data
213 % If there is a regular interval between measurements, calculation is much
214 % easier/faster
215 if isfield(data,'rate') && data.rate > 0 % if data rate was given
216 if verbose >= 1
217 fprintf(1, 'allan_overlap: regular data ');
218 if isfield(data,'freq')
219 fprintf(1, '(%g freq data points @ %g Hz)\n',length(data.freq),data.rate);
220 elseif isfield(data,'phase')
221 fprintf(1, '(%g phase data points @ %g Hz)\n',length(data.phase),data.rate);
222 else
223 error('\n phase or freq data missing [err10]');
224 end
225 end
226
227 % string for plot title
228 name=[name ' (' num2str(data.rate) ' Hz)'];
229
230 % what is the time interval between data points?
231 tmstep = 1/data.rate;
232
233 % Is there time data? Just for curiosity/plotting, does not impact calculation
234 if isfield(data,'time')
235 % adjust time data to remove any starting gap; first time step
236 % should not be zero for comparison with freq data
237 dtime=data.time-data.time(1)+mean(diff(data.time));
238 dtime=dtime(1:length(medianfreq)); % equalize the data vector lengths for plotting (v2.1)
239 if verbose >= 2
240 fprintf(1,'allan_overlap: End of timestamp data: %g sec.\n',dtime(end));
241 if (data.rate - 1/mean(diff(dtime))) > 1e-6
242 fprintf(1,'allan_overlap: NOTE: data.rate (%f Hz) does not match average timestamped sample rate (%f Hz)\n',data.rate,1/mean(diff(dtime)));
243 end
244 end
245 else
246 % create time axis data using rate (for plotting only)
247 dtime=(tmstep:tmstep:length(data.freq)*tmstep);
248 end
249
250
251 % is phase data present? If not, generate it
252 if ~isfield(data,'phase')
253 nfreq=data.freq-s.mean;
254 dphase=zeros(1,length(nfreq)+1);
255 dphase(2:end) = cumsum(nfreq)./data.rate;
256 if verbose >= 1, fprintf(1,'allan_overlap: phase data generated from fractional frequency data (N=%g).\n',length(dphase)); end
257 else
258 dphase=data.phase;
259 end
260
261 % check the range of tau values and truncate if necessary
262 % find halfway point of time record
263 halftime = round(tmstep*length(data.freq)/2);
264 % truncate tau to appropriate values
265 tau = tau(tau >= tmstep & tau <= halftime);
266 if verbose >= 2, fprintf(1, 'allan_overlap: allowable tau range: %g to %g sec. (1/rate to total_time/2)\n',tmstep,halftime); end
267
268 % number of samples
269 N=length(dphase);
270 % number of samples per tau period
271 m = data.rate.*tau;
272 % only integer values allowed for m (no fractional groups of points)
273 %tau = tau(m-round(m)<1e-8); % numerical precision issues (v2.1)
274 tau = tau(m==round(m)); % The round() test is only correct for values < 2^53
275 %m = m(m-round(m)<1e-8); % change to round(m) for integer test v2.22
276 m = m(m==round(m));
277 %m=round(m);
278 %fprintf(1,'m: %.50f\n',m)
279
280 if verbose >= 1, fprintf(1,'allan_overlap: calculating overlapping Allan deviation...\n '); end
281
282 % calculate the Allan deviation for each value of tau
283 k=0; tic;
284 for i = tau
285 k=k+1;
286 if verbose >= 2, fprintf(1,'%d ',i); end
287
288
289 % pad phase data set length to an even multiple of this tau value
290 mphase=zeros(ceil(length(dphase)./m(k))*m(k),1);
291 mphase(1:N)=dphase;
292 % group phase values
293 mp=reshape(mphase,m(k),[]);
294 % compute second differences of phase values (x_k+m - x_k)
295 md1=diff(mp,1,2);
296 md2=diff(md1,1,2);
297 md1=reshape(md2,1,[]);
298
299 % compute overlapping ADEV from phase values
300 % only the first N-2*m(k) samples are valid
301 sm(k)=sqrt((1/(2*(N-2*m(k))*i^2))*sum(md1(1:N-2*m(k)).^2));
302
303 % estimate error bars
304 sme(k)=sm(k)/sqrt(N-2*m(k));
305
306
307 end % repeat for each value of tau
308
309 if verbose >= 2, fprintf(1,'\n'); end
310 calctime=toc; if verbose >= 2, fprintf(1,'allan_overlap: Elapsed time for calculation: %g seconds\n',calctime); end
311
312
313
314 %% Irregular data, no fixed interval
315 elseif isfield(data,'time')
316 % the interval between measurements is irregular
317 % so we must group the data by time
318 if verbose >= 1, fprintf(1, 'allan_overlap: irregular rate data (no fixed sample rate)\n'); end
319
320
321 % string for plot title
322 name=[name ' (timestamp)'];
323
324
325 % adjust time to remove any starting offset
326 dtime=data.time-data.time(1)+mean(diff(data.time));
327
328 % save the freq data for the loop
329 dfreq=data.freq;
330 dtime=dtime(1:length(dfreq));
331
332 dfdtime=diff(dtime); % only need to do this once (v2.1)
333 % where is the maximum gap in time record?
334 gap_pos=find(dfdtime==max(dfdtime));
335 % what is average data spacing?
336 avg_gap = mean(dfdtime);
337 s.avg_rate = 1/avg_gap; % save avg rate for user (v2.1)
338
339 if verbose >= 2
340 fprintf(1, 'allan_overlap: WARNING: irregular timestamp data (no fixed sample rate).\n');
341 fprintf(1, ' Calculation time may be long and the results subject to interpretation.\n');
342 fprintf(1, ' You are advised to estimate using an average sample rate (%g Hz) instead of timestamps.\n',1/avg_gap);
343 fprintf(1, ' Continue at your own risk! (press any key to continue)\n');
344 pause;
345 end
346
347 if verbose >= 1
348 fprintf(1, 'allan_overlap: End of timestamp data: %g sec\n',dtime(end));
349 fprintf(1, ' Average rate: %g Hz (%g sec/measurement)\n',1/avg_gap,avg_gap);
350 if max(diff(dtime)) ~= 1/mean(diff(dtime))
351 fprintf(1, ' Max. gap in time record: %g sec at position %d\n',max(dfdtime),gap_pos(1));
352 end
353 if max(diff(dtime)) > 5*avg_gap
354 fprintf(1, ' WARNING: Max. gap in time record is suspiciously large (>5x the average interval).\n');
355 end
356 end
357
358
359 % find halfway point
360 halftime = fix(dtime(end)/2);
361 % truncate tau to appropriate values
362 tau = tau(tau >= max(dfdtime) & tau <= halftime);
363 if isempty(tau)
364 error('allan_overlap: ERROR: no appropriate tau values (> %g s, < %g s)\n',max(dfdtime),halftime);
365 end
366
367
368 % number of samples
369 M=length(dfreq);
370 % number of samples per tau period
371 m=round(tau./avg_gap);
372
373 if verbose >= 1, fprintf(1,'allan_overlap: calculating overlapping Allan deviation...\n'); end
374
375 k=0; tic;
376 for i = tau
377 k=k+1;
378 fa=[];
379
380 if verbose >= 2, fprintf(1,'%d ',i); end
381
382 freq = dfreq; time = dtime;
383
384
385 % compute overlapping samples (y_k) for this tau
386 %for j = 1:i
387 for j = 1:m(k) % (v2.1)
388 km=0;
389 %fprintf(1,'j: %d ',j);
390
391 % (v2.1) truncating not correct for overlapping samples
392 % truncate data set to an even multiple of this tau value
393 %freq = freq(time <= time(end)-rem(time(end),i));
394 %time = time(time <= time(end)-rem(time(end),i));
395
396 % break up the data into overlapping groups of tau length
397 while i*km <= time(end)
398 km=km+1;
399 %i*km
400
401 % progress bar
402 if verbose >= 2
403 if rem(km,100)==0, fprintf(1,'.'); end
404 if rem(km,1000)==0, fprintf(1,'%g/%g\n',km,round(time(end)/i)); end
405 end
406
407 f = freq(i*(km-1) < (time) & (time) <= i*km);
408
409 if ~isempty(f)
410 fa(j,km)=mean(f);
411 else
412 fa(j,km)=0;
413 end
414
415 end
416 %fa
417
418 % shift data vector by -1 and repeat
419 freq=circshift(dfreq,(size(freq)>1)*-j);
420 freq(end-j+1:end)=[];
421 time=circshift(dtime,(size(time)>1)*-j);
422 time(end-j+1:end)=[];
423 time=time-time(1)+avg_gap; % remove time offset
424
425 end
426
427 % compute second differences of fractional frequency values (y_k+m - y_k)
428 fd1=diff(fa,1,2);
429 fd1=reshape(fd1,1,[]);
430 % compute overlapping ADEV from fractional frequency values
431 % only the first M-2*m(k)+1 samples are valid
432 if length(fd1) >= M-2*m(k)+1
433 sm(k)=sqrt((1/(2*(M-2*m(k)+1)))*sum(fd1(1:M-2*m(k)+1).^2));
434
435 % estimate error bars
436 sme(k)=sm(k)/sqrt(M+1);
437
438 if verbose >= 2, fprintf(1,'\n'); end
439
440 else
441 if verbose >=2, fprintf(1,' tau=%g dropped due to timestamp irregularities\n',tau(k)); end
442 sm(k)=0; sme(k)=0;
443 end
444
445
446 end
447
448 if verbose >= 2, fprintf(1,'\n'); end
449 calctime=toc; if verbose >= 1, fprintf(1,'allan_overlap: Elapsed time for calculation: %g seconds\n',calctime); end
450
451 % remove any points that were dropped
452 tau(sm==0)=[];
453 sm(sm==0)=[];
454 sme(sme==0)=[];
455
456
File was created 1 #!/usr/bin/octave-cli --persist
2
3 filename = argv(){1};
4 col = eval(argv(){2});
5 mult = eval(argv(){3});
6
7 data.freq = load(filename)(:,col).*mult;
8 data.rate = 1;
9
File was created 1 function hL = dsplot(x, y, numPoints)
2
3 %DSPLOT Create down sampled plot.
4 % This function creates a down sampled plot to improve the speed of
5 % exploration (zoom, pan).
6 %
7 % DSPLOT(X, Y) plots Y versus X by downsampling if there are large number
8 % of elements. X and Y needs to obey the following:
9 % 1. X must be a monotonically increasing vector.
10 % 2. If Y is a vector, it must be the same size as X.
11 % 3. If Y is a matrix, one of the dimensions must line up with X.
12 %
13 % DSPLOT(Y) plots the columns of Y versus their index.
14 %
15 % hLine = DSPLOT(X, Y) returns the handles of the line. Note that the
16 % lines may be downsampled, so they may not represent the full data set.
17 %
18 % DSPLOT(X, Y, NUMPOINTS) or DSPLOT(Y, [], NUMPOINTS) specifies the
19 % number of points (roughly) to display on the screen. The default is
20 % 50000 points (~390 kB doubles). NUMPOINTS can be a number greater than
21 % 500.
22 %
23 % It is very likely that more points will be displayed than specified by
24 % NUMPOINTS, because it will try to plot any outlier points in the range.
25 % If the signal is stochastic or has a lot of sharp changes, there will
26 % be more points on plotted on the screen.
27 %
28 % The figure title (name) will indicate whether the plot shown is
29 % downsampled or is the true representation.
30 %
31 % The figure can be saved as a .fig file, which will include the actual
32 % data. The figure can be reloaded and the actual data can be exported to
33 % the base workspace via a menu.
34 %
35 % Run the following examples and zoom/pan to see the performance.
36 %
37 % Example 1: (with small details)
38 % x = linspace(0, 2*pi, 1000000);
39 % y1 = sin(x)+.02*cos(200*x)+0.001*sin(2000*x)+0.0001*cos(20000*x);
40 % dsplot(x,y1);title('Down Sampled');
41 % % compare with
42 % figure;plot(x,y1);title('Normal Plot');
43 %
44 % Example 2: (with outlier points)
45 % x = linspace(0, 2*pi, 1000000);
46 % y1 = sin(x) + .01*cos(200*x) + 0.001*sin(2000*x);
47 % y2 = sin(x) + 0.3*cos(3*x) + 0.001*randn(size(x));
48 % y1([300000, 700000, 700001, 900000]) = [0, 1, -2, 0.5];
49 % y2(300000:500000) = y2(300000:500000) + 1;
50 % y2(500001:600000) = y2(500001:600000) - 1;
51 % y2(800000) = 0;
52 % dsplot(x, [y1;y2]);title('Down Sampled');
53 % % compare with
54 % figure;plot(x, [y1;y2]);title('Normal Plot');
55 %
56 % See also PLOT.
57
58 % Version:
59 % v1.0 - first version (Aug 1, 2007)
60 % v1.1 - added CreateFcn for the figure so that when the figure is saved
61 % and re-loaded, the zooming and panning works. Also added a menu
62 % item for saving out the original data back to the base
63 % workspace. (Aug 10, 2007)
64 %
65 % Jiro Doke
66 % August 1, 2007
67
68 debugMode = false;
69
70 %--------------------------------------------------------------------------
71 % Error checking
72 error(nargchk(1, 3, nargin, 'struct'));
73 if nargin < 3
74 % Number of points to show on the screen. It's quite possible that more
75 % points will be displayed if there are outlier points
76 numPoints = 50000; % ~390 kB for doubles
77 end
78 if nargin == 1 || isempty(y)
79 noXVar = true;
80 y = x;
81 x = [];
82 else
83 noXVar = false;
84 end
85 myErrorCheck;
86 %--------------------------------------------------------------------------
87
88 if size(x, 2) > 1 % it's a row vector -> transpose
89 x = x';
90 y = y';
91 varTranspose = true;
92 else
93 varTranspose = false;
94 end
95
96 % Number of lines
97 numSignals = size(y, 2);
98
99 % If the number of lines is greater than the number of data points per
100 % line, it's possible that the user may have mistaken the matrix
101 % orientation.
102 if numSignals > size(y, 1)
103 s = input(sprintf('Are you sure you want to plot %d lines? (y/n) ', ...
104 numSignals), 's');
105 if ~strcmpi(s, 'y')
106 disp('Canceled. You may want to transpose the matrix.');
107 if nargout == 1
108 hL = [];
109 end
110 return;
111 end
112 end
113
114 % Attempt to find outliers. Use a running average technique
115 filterWidth = ceil(min([50, length(x)/10])); % max window size of 50
116 a = y - filter(ones(filterWidth,1)/filterWidth, 1, y);
117 [iOutliers, jOutliers] = find(abs(a - repmat(mean(a), size(a, 1), 1)) > ...
118 repmat(4 * std(a), size(a, 1), 1));
119 clear a;
120
121 % Always create new figure because it messes around with zoom, pan,
122 % datacursors.
123 hFig = figure;
124 figName = '';
125
126 % Create template plot using NaNs
127 hLine = plot(NaN(2, numSignals), NaN(2, numSignals));
128 set(hLine, 'tag', 'dsplot_lines');
129
130 % Define CreateFcn for the figure
131 set(hFig, 'CreateFcn', @mycreatefcn);
132 mycreatefcn();
133
134 % Create menu for exporting data
135 hMenu = uimenu(hFig, 'Label', 'Data');
136 uimenu(hMenu, ...
137 'Label' , 'Export data to workspace.', ...
138 'Callback', @myExportFcn);
139
140 % Update lines
141 updateLines([min(x), max(x)]);
142
143 % Deal with output argument
144 if nargout == 1
145 hL = hLine;
146 end
147
148 %--------------------------------------------------------------------------
149 function myExportFcn(varargin)
150 % This callback allows for extracting the actual data from the figure.
151 % This means that if you save this figure and load it back later, you
152 % can get back the data.
153
154 % Determine the variable name
155 allVarNames = evalin('base', 'who');
156 newVarName = genvarname('dsplotData', allVarNames);
157
158 % X
159 if ~noXVar
160 if varTranspose
161 dat.x = x';
162 else
163 dat.x = x;
164 end
165 end
166
167 % Y
168 if varTranspose
169 dat.y = y';
170 else
171 dat.y = y;
172 end
173
174 assignin('base', newVarName, dat);
175
176 msgbox(sprintf('Data saved to the base workspace as ''%s''.', ...
177 newVarName), 'Saved', 'modal');
178
179 end
180
181 %--------------------------------------------------------------------------
182 function mycreatefcn(varargin)
183 % This callback defines the custom zoom/pan functions. It is defined as
184 % the CreateFcn of the figure, so it allows for saving and reloading of
185 % the figure.
186
187 if nargin > 0
188 hFig = varargin{1};
189 end
190 hLine = findobj(hFig, 'type', 'axes');
191 hLine(strmatch('legend', get(hLine, 'tag'))) = [];
192 hLine = get(hLine, 'Children');
193
194 % Create Zoom, Pan, Datacursor objects
195 hZoom = zoom(hFig);
196 hPan = pan(hFig);
197 hDc = datacursormode(hFig);
198 set(hZoom, 'ActionPostCallback', @mypostcallback);
199 set(hPan , 'ActionPostCallback', @mypostcallback);
200 set(hDc , 'UpdateFcn' , @myDCupdatefcn);
201
202 end
203
204 %--------------------------------------------------------------------------
205 function mypostcallback(obj, evd) %#ok
206 % This callback that gets called when the mouse is released after
207 % zooming or panning.
208
209 % single or double-click
210 switch get(hFig, 'SelectionType')
211 case {'normal', 'alt'}
212 updateLines(xlim(evd.Axes));
213
214 case 'open'
215 updateLines([min(x), max(x)]);
216
217 end
218
219 end
220
221 %--------------------------------------------------------------------------
222 function updateLines(rng)
223 % This helper function is for determining the points to plot on the
224 % screen based on which portion is visible in the current limits.
225
226 % find indeces inside the range
227 id = find(x >= rng(1) & x <= rng(2));
228
229 % if there are more points than we want
230 if length(id) > numPoints / numSignals
231
232 % see how many outlier points are in this range
233 blah = iOutliers > id(1) & iOutliers < id(end);
234
235 % determine indeces of points to plot.
236 idid = round(linspace(id(1), id(end), round(numPoints/numSignals)))';
237
238 x2 = cell(numSignals, 1);
239 y2 = x2;
240 for iSignals = 1:numSignals
241 % add outlier points
242 ididid = unique([idid; iOutliers(blah & jOutliers == iSignals)]);
243 x2{iSignals} = x(ididid);
244 y2{iSignals} = y(ididid, iSignals);
245 end
246
247 if debugMode
248 figName = ['downsampled - ', sprintf('%d, ', cellfun('length', y2))];
249 else
250 figName = 'downsampled';
251 end
252
253 else % no need to down sample
254 figName = 'true';
255
256 x2 = repmat({x(id)}, numSignals, 1);
257 y2 = mat2cell(y(id, :), length(id), ones(1, numSignals))';
258
259 end
260
261 % Update plot
262 set(hLine, {'xdata', 'ydata'} , [x2, y2]);
263 set(hFig, 'Name', figName);
264
265 end
266
267 %--------------------------------------------------------------------------
268 function txt = myDCupdatefcn(empt, event_obj) %#ok
269 % This function displays appropriate data cursor message based on the
270 % display type
File was created 1 function T = pt100(R)
2
3 pt100 = @(R) (R/100.-1)/0.003850+273.15;
File was created 1 function T = res2temp16941(R)
2
3 res2temp16941 = @(R) 10.^(2.9486 * (log10(1000./R)).^2 + 4.5862 * log10(1000./R) + 2.266);
File was created 1 function T = res2temp16943(R)
2
3 res2temp16943 = @(R) 10.^(3.4738 * (log10(1000./R)).^2 + 5.1198 * log10(1000./R) + 2.3681);
File was created 1 function T = res2temp16944(R)
2
3 res2temp16944 = @(R) 10.^(3.3674 * (log10(1000./R)).^2 + 5.2874 * log10(1000./R) + 2.5165);
File was created 1 function T = res2temp16945(R)
2
3 res2temp16945 = @(R) 10.^(3.2497 * (log10(1000./R)).^2 + 5.1777 * log10(1000./R) + 2.499);
File was created 1 function T = res2temp16947(R)
2
3 res2temp16947 = @(R) 10.^(3.4597 * (log10(1000./R)).^2 + 5.2422 * log10(1000./R) + 2.4169);
File was created 1 function T = res2temp625(R)
2
3 res2temp625 = @(R) 0.333548856582638109 + 11.7361551595386118 * (1000./R) + -31.32988932320903987 * (1000./R).^2 + 262.878643524833024 * (1000./R).^3 + -704.163538021035492 * (1000./R).^4 + 1056.6040485650301 * (1000./R).^5 + -307.057196729816496 * (1000./R).^6;
File was created 1 function T = res2temp627(R)
2
3 res2temp627 = @(R) 0.399341181655472610 + 10.8420092277810909 * (1000./R) + -26.4597939187660813 * (1000./R).^2 + 245.9828566655493379 * (1000./R).^3 + -668.069876596331596 * (1000./R).^4 + 1001.69882618263364 * (1000./R) .^5 + -267.272089680656791 * (1000./R).^6;
File was created 1 function T = res2temp628(R)
2
3 res2temp628 = @(R) 0.463200932294057566 + 13.5049710820894688 * (1000./R) + -30.5191222755238414 * (1000./R).^2 + 231.098593852017075* (1000./R).^3 + -550.122691885568202 * (1000./R).^4 + 806.038547554984689 * (1000./R).^5 + -198.510489917360246 * (1000./R).^6;