bmarechal / octave_scripts

Download as
Browse Code ยป

Commit b3b851a23e6321841a2d789d31d3d48b06557378

Authored by bmarechal
1 parent 3124193a36
Exists in master

add cov test scripts

Showing 2 changed files with 616 additions and 0 deletions Inline Diff

File was created 1 function [retval, s, errorb, tau] = allan_cov(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 dfreq2=data.freq2;
300 % find the number of data points in each tau group
301 m = data.rate.*tau;
302 % only integer values allowed (no fractional groups of points)
303 %tau = tau(m-round(m)<1e-8); % numerical precision issues (v2.1)
304 tau = tau(m==round(m)); % The round() test is only correct for values < 2^53
305 %m = m(m-round(m)<1e-8); % change to round(m) for integer test v2.22
306 m = m(m==round(m));
307 %m=round(m);
308
309 if verbose >= 1, fprintf(1,'allan: calculating Allan deviation...\n '); end
310
311 % calculate the Allan deviation for each value of tau
312 k=0; tic;
313 for i = tau
314 if verbose >= 2, fprintf(1,'%g ',i); end
315 k=k+1;
316
317 % truncate frequency set to an even multiple of this tau value
318 freq=dfreq(1:end-rem(length(dfreq),m(k)));
319 freq2=dfreq2(1:end-rem(length(dfreq2),m(k)));
320 % group the data into tau-length groups or bins
321 f = reshape(freq,m(k),[]); % Vectorize!
322 f2 = reshape(freq2,m(k),[]); % Vectorize!
323 % find average in each "tau group", y_k (each colummn of f)
324 fa=mean(f,1);
325 fa2=mean(f2,1);
326 % first finite difference
327 fd=diff(fa);
328 fd2=diff(fa2);
329 % calculate two-sample variance for this tau
330 M=length(fa);
331 sm(k)=sqrt(0.5/(M-1)*(sum(fd.*fd2)));
332
333 % estimate error bars
334 sme(k)=sm(k)/sqrt(M+1);
335
336 if TAUBIN == 1
337 % save the binning points for plotting
338 fs(k,1:length(freq)/m(k))=m(k):m(k):length(freq); fval{k}=mean(f,1);
339 end
340
341 end % repeat for each value of tau
342
343 if verbose >= 2, fprintf(1,'\n'); end
344 calctime=toc; if verbose >= 2, fprintf(1,'allan: Elapsed time for calculation: %e seconds\n',calctime); end
345
346
347
348 %% Irregular data (timestamp)
349 elseif isfield(data,'time')
350 % the interval between measurements is irregular
351 % so we must group the data by time
352 if verbose >= 1, fprintf(1, 'allan: irregular rate data (no fixed sample rate)\n'); end
353
354 % string for plot title
355 name=[name ' (timestamp)'];
356
357 % adjust time to remove any initial offset or zero
358 dtime=data.time-data.time(1)+mean(diff(data.time));
359 %dtime=data.time;
360 % where is the maximum gap in time record?
361 gap_pos=find(diff(dtime)==max(diff(dtime)));
362 % what is average data spacing?
363 avg_gap = mean(diff(dtime));
364
365 if verbose >= 2
366 fprintf(1, 'allan: WARNING: irregular timestamp data (no fixed sample rate).\n');
367 fprintf(1, ' Calculation time may be long and the results subject to interpretation.\n');
368 fprintf(1, ' You are advised to estimate using an average sample rate (%g Hz) instead of timestamps.\n',1/avg_gap);
369 fprintf(1, ' Continue at your own risk! (press any key to continue)\n');
370 pause;
371 end
372
373 if verbose >= 1
374 fprintf(1, 'allan: End of timestamp data: %g sec\n',dtime(end));
375 fprintf(1, ' Average rate: %g Hz (%g sec/measurement)\n',1/avg_gap,avg_gap);
376 if max(diff(dtime)) ~= 1/mean(diff(dtime))
377 fprintf(1, ' Max. gap: %g sec at position %d\n',max(diff(dtime)),gap_pos(1));
378 end
379 if max(diff(dtime)) > 5*avg_gap
380 fprintf(1, ' WARNING: Max. gap in time record is suspiciously large (>5x the average interval).\n');
381 end
382 end
383
384
385 % find halfway point
386 halftime = fix(dtime(end)/2);
387 % truncate tau to appropriate values
388 tau = tau(tau >= max(diff(dtime)) & tau <= halftime);
389 if isempty(tau)
390 error('allan: ERROR: no appropriate tau values (> %g s, < %g s)\n',max(diff(dtime)),halftime);
391 end
392
393 % save the freq data for the loop
394 dfreq=data.freq;
395 dtime=dtime(1:length(dfreq));
396
397 if verbose >= 1, fprintf(1,'allan: calculating Allan deviation...\n'); end
398
399 k=0; tic;
400 for i = tau
401 if verbose >= 2, fprintf(1,'%d ',i); end
402
403 k=k+1; fa=[]; %f=[];
404 km=0;
405
406 % truncate data set to an even multiple of this tau value
407 freq=dfreq(dtime <= dtime(end)-rem(dtime(end),i));
408 time=dtime(dtime <= dtime(end)-rem(dtime(end),i));
409 %freq=dfreq;
410 %time=dtime;
411
412 % break up the data into groups of tau length in sec
413 while i*km < time(end)
414 km=km+1;
415
416 % progress bar
417 if verbose >= 2
418 if rem(km,100)==0, fprintf(1,'.'); end
419 if rem(km,1000)==0, fprintf(1,'%g/%g\n',km,round(time(end)/i)); end
420 end
421
422 f = freq(i*(km-1) < time & time <= i*km);
423 f = f(~isnan(f)); % make sure values are valid
424
425 if ~isempty(f)
426 fa(km)=mean(f);
427 else
428 fa(km)=0;
429 end
430
431 if TAUBIN == 1 % WARNING: this has a significant impact on performance
432 % save the binning points for plotting
433 %if find(time <= i*km) > 0
434 fs(k,km)=max(time(time <= i*km));
435 %else
436 if isempty(fs(k,km))
437 fs(k,km)=0;
438 end
439 fval{k}=fa;
440 end % save tau bin plot points
441
442 end
443
444 if verbose >= 2, fprintf(1,'\n'); end
445
446 % first finite difference of the averaged results
447 fd=diff(fa);
448 % calculate Allan deviation for this tau
449 M=length(fa);
450 sm(k)=sqrt(0.5/(M-1)*(sum(fd.^2)));
451
452 % estimate error bars
453 sme(k)=sm(k)/sqrt(M+1);
454
455
456 end
457
458 if verbose == 2, fprintf(1,'\n'); end
459 calctime=toc; if verbose >= 2, fprintf(1,'allan: Elapsed time for calculation: %e seconds\n',calctime); end
460
461
462 else
463 error('allan: WARNING: no DATA.rate or DATA.time! Type "help allan" for more information. [err2]');
464 end
465
466
467 %%%%%%%%
468 %% Plotting
469
470 if verbose >= 2 % show all data
471
472 % plot the frequency data, centered on median
473 if size(dtime,2) > size(dtime,1), dtime=dtime'; end % this should not be necessary, but dsplot 1.1 is a little bit brittle
474 try
475 % dsplot makes a new figure
476 hd=dsplot(dtime,medianfreq);
477 catch ME
478 figure;
479 if length(dtime) ~= length(medianfreq)
480 fprintf(1,'allan: Warning: length of time axis (%d) is not equal to data array (%d)\n',length(dtime),length(medianfreq));
481 end
482 hd=plot(dtime,medianfreq);
483 if verbose >= 1, fprintf(1,'allan: Note: Install dsplot.m for improved plotting of large data sets (File Exchange File ID: #15850).\n'); end
484 if verbose >= 2, fprintf(1,' (Message: %s)\n',ME.message); end
485 end
486 set(hd,'Marker','.','LineStyle','none','Color','b'); % equivalent to '.-'
487 hold on;
488
489 % show center (0)
490 plot(xlim,[0 0],':k');
491 % show 5x Median Absolute Deviation (MAD) values
492 hm=plot(xlim,[5*s.MAD 5*s.MAD],'-r');
493 plot(xlim,[-5*s.MAD -5*s.MAD],'-r');
494 % show linear fit line
495 hf=plot(xlim,polyval(s.linear,xlim)-s.median,'-g');
496 title(['Data: ' name],'FontSize',FontSize+2,'FontName',FontName);
497 %set(get(gca,'Title'),'Interpreter','none');
498 xlabel('Time [sec]','FontSize',FontSize,'FontName',FontName);
499 if isfield(data,'units')
500 ylabel(['data - median(data) [' data.units ']'],'FontSize',FontSize,'FontName',FontName);
501 else
502 ylabel('freq - median(freq)','FontSize',FontSize,'FontName',FontName);
503 end
504 set(gca,'FontSize',FontSize,'FontName',FontName);
File was created 1 #!/usr/bin/octave-cli --persist
2
3 filename = argv(){1};
4 col1 = eval(argv(){2});
5 col2 = eval(argv(){3});
6 mult1 = eval(argv(){4});
7 mult2 = eval(argv(){5});
8
9 if length(col1) == length(mult1)
10 figure
11 hold all
12 grid on
13 cc = 'bkcgmry';
14 for i = [1:length(col1)]
15 data.freq = load(filename)(:,col1(i)).*mult1(i);
16 data.freq2 = load(filename)(:,col2(i)).*mult2(i);
17 if eval(argv(){end-1}) == 1
18 printf('\ndata1 drift removed\n\n')
19 data.freq = detrend(data.freq);
20 end
21 if eval(argv(){end}) == 1
22 printf('\ndata2 drift removed\n\n')
23 data.freq2 = detrend(data.freq2);
24 end
25 data.rate = 1;
26 [ad, S, err, tau] = allan_cov(data, 2.^[0:nextpow2(length(data.freq))-3]./data.rate, strcat(strsplit(filename, '/'){end}, num2str(i)), 0);
27 loglogerr(tau, ad, err, strcat(cc(mod(i, length(cc))), '-s'))
28 leg{i} = strcat(filename, ' cov col', num2str(col1(i)), ' col', num2str(col2(i)));
29 axis(10.^ceil(log10([tau(1), tau(end)])))
30 hold on
31 end
32 legend(leg)
33 input("Press to continue...");