bmarechal / octave_scripts

Download as
Browse Code ยป

Commit e5aeec49133ba19a381291ad6dcad1164d621c61

Authored by bmarechal
1 parent 8ca6bc2fe9
Exists in master

use real part

Showing 1 changed file with 1 additions and 1 deletions Inline Diff

function [retval, s, errorb, tau] = allan_cov(data,tau,name,verbose) 1 1 function [retval, s, errorb, tau] = allan_cov(data,tau,name,verbose)
% ALLAN Compute the Allan deviation for a set of time-domain frequency data 2 2 % ALLAN Compute the Allan deviation for a set of time-domain frequency data
% [RETVAL, S, ERRORB, TAU] = ALLAN(DATA,TAU,NAME,VERBOSE) 3 3 % [RETVAL, S, ERRORB, TAU] = ALLAN(DATA,TAU,NAME,VERBOSE)
% 4 4 %
% Inputs: 5 5 % Inputs:
% DATA should be a structure and have the following fields: 6 6 % DATA should be a structure and have the following fields:
% DATA.freq or DATA.phase 7 7 % DATA.freq or DATA.phase
% A vector of fractional frequency measurements (df/f) in 8 8 % A vector of fractional frequency measurements (df/f) in
% DATA.freq *or* phase offset data (seconds) in DATA.phase . 9 9 % DATA.freq *or* phase offset data (seconds) in DATA.phase .
% If frequency data is not present, it will be generated by 10 10 % If frequency data is not present, it will be generated by
% differentiating the phase data. 11 11 % differentiating the phase data.
% If both fields are present, then DATA.freq will be used. 12 12 % If both fields are present, then DATA.freq will be used.
% Note: for general-purpose calculations of Allan deviation, 13 13 % Note: for general-purpose calculations of Allan deviation,
% (i.e. a two-sample variance) use DATA.freq . 14 14 % (i.e. a two-sample variance) use DATA.freq .
% 15 15 %
% DATA.rate or DATA.time 16 16 % DATA.rate or DATA.time
% The sampling rate in Hertz (DATA.rate) or a vector of 17 17 % The sampling rate in Hertz (DATA.rate) or a vector of
% timestamps for each measurement in seconds (DATA.time). 18 18 % timestamps for each measurement in seconds (DATA.time).
% DATA.rate is used if both fields are present. 19 19 % DATA.rate is used if both fields are present.
% If DATA.rate == 0, then the timestamps are used. 20 20 % If DATA.rate == 0, then the timestamps are used.
% 21 21 %
% DATA.units (optional) 22 22 % DATA.units (optional)
% The units for the data. If present, the string DATA.units 23 23 % The units for the data. If present, the string DATA.units
% is added to the plot y-axis label. 24 24 % is added to the plot y-axis label.
% 25 25 %
% TAU is an array of tau values for computing Allan deviation. 26 26 % TAU is an array of tau values for computing Allan deviation.
% TAU values must be divisible by 1/DATA.rate (data points cannot be 27 27 % TAU values must be divisible by 1/DATA.rate (data points cannot be
% grouped in fractional quantities!) and invalid values are ignored. 28 28 % grouped in fractional quantities!) and invalid values are ignored.
% Leave empty to use default values. 29 29 % Leave empty to use default values.
% NAME is an optional label that is added to the plot titles. 30 30 % NAME is an optional label that is added to the plot titles.
% VERBOSE sets the level of status messages: 31 31 % VERBOSE sets the level of status messages:
% 0 = silent & no data plots; 32 32 % 0 = silent & no data plots;
% 1 = status messages & minimum plots; 33 33 % 1 = status messages & minimum plots;
% 2 = all messages and plots (default) 34 34 % 2 = all messages and plots (default)
% 35 35 %
% Outputs: 36 36 % Outputs:
% RETVAL is the array of Allan deviation values at each TAU. 37 37 % RETVAL is the array of Allan deviation values at each TAU.
% S is an optional output of other statistical measures of the data (mean, std, etc). 38 38 % S is an optional output of other statistical measures of the data (mean, std, etc).
% ERRORB is an optional output containing the error estimates for a 1-sigma 39 39 % ERRORB is an optional output containing the error estimates for a 1-sigma
% confidence interval. These values are shown on the figure for each point. 40 40 % confidence interval. These values are shown on the figure for each point.
% TAU is an optional output containing the array of tau values used in the 41 41 % TAU is an optional output containing the array of tau values used in the
% calculation (which may be a truncated subset of the input or default values). 42 42 % calculation (which may be a truncated subset of the input or default values).
% 43 43 %
% Example: 44 44 % Example:
% 45 45 %
% To compute the Allan deviation for the data in the variable "lt": 46 46 % To compute the Allan deviation for the data in the variable "lt":
% >> lt 47 47 % >> lt
% lt = 48 48 % lt =
% freq: [1x86400 double] 49 49 % freq: [1x86400 double]
% rate: 0.5 50 50 % rate: 0.5
% 51 51 %
% Use: 52 52 % Use:
% 53 53 %
% >> ad = allan(lt,[2 10 100],'lt data',1); 54 54 % >> ad = allan(lt,[2 10 100],'lt data',1);
% 55 55 %
% The Allan deviation will be computed and plotted at tau = 2,10,100 seconds. 56 56 % The Allan deviation will be computed and plotted at tau = 2,10,100 seconds.
% 1-sigma confidence intervals will be indicated by vertical lines at each point. 57 57 % 1-sigma confidence intervals will be indicated by vertical lines at each point.
% You can also use the default settings, which are usually a good starting point: 58 58 % You can also use the default settings, which are usually a good starting point:
% 59 59 %
% >> ad = allan(lt); 60 60 % >> ad = allan(lt);
% 61 61 %
% 62 62 %
% Notes: 63 63 % Notes:
% This function calculates the standard Allan deviation (ADEV), *not* the 64 64 % This function calculates the standard Allan deviation (ADEV), *not* the
% overlapping ADEV. Use "allan_overlap.m" for overlapping ADEV. 65 65 % overlapping ADEV. Use "allan_overlap.m" for overlapping ADEV.
% The calculation is performed using fractional frequency data. If only 66 66 % The calculation is performed using fractional frequency data. If only
% phase data is provided, frequency data is generated by differentiating 67 67 % phase data is provided, frequency data is generated by differentiating
% the phase data. 68 68 % the phase data.
% No pre-processing of the data is performed, except to remove any 69 69 % No pre-processing of the data is performed, except to remove any
% initial offset (i.e., starting gap) in the time record. 70 70 % initial offset (i.e., starting gap) in the time record.
% For rate-based data, ADEV is computed only for tau values greater than the 71 71 % For rate-based data, ADEV is computed only for tau values greater than the
% minimum time between samples and less than the half the total time. For 72 72 % minimum time between samples and less than the half the total time. For
% time-stamped data, only tau values greater than the maximum gap between 73 73 % time-stamped data, only tau values greater than the maximum gap between
% samples and less than half the total time are used. 74 74 % samples and less than half the total time are used.
% The calculation for fixed sample rate data is *much* faster than for 75 75 % The calculation for fixed sample rate data is *much* faster than for
% time-stamp data. You may wish to run the rate-based calculation first, 76 76 % time-stamp data. You may wish to run the rate-based calculation first,
% then compare with time-stamp-based. Often the differences are insignificant. 77 77 % then compare with time-stamp-based. Often the differences are insignificant.
% To show the "tau bins" (y_k samples) on the data plot, set the variable 78 78 % To show the "tau bins" (y_k samples) on the data plot, set the variable
% TAUBIN to 1 (search for "#TAUBIN"). 79 79 % TAUBIN to 1 (search for "#TAUBIN").
% You can choose between loglog and semilog plotting of results by 80 80 % You can choose between loglog and semilog plotting of results by
% commenting in/out the appropriate line. Search for "#PLOTLOG". 81 81 % commenting in/out the appropriate line. Search for "#PLOTLOG".
% I recommend installing "dsplot.m", which improves the performance of 82 82 % I recommend installing "dsplot.m", which improves the performance of
% plotting large data sets. Download from File Exchange, File ID: #15850. 83 83 % plotting large data sets. Download from File Exchange, File ID: #15850.
% allan.m will use dsplot.m if it is present on your MATLAB path. 84 84 % allan.m will use dsplot.m if it is present on your MATLAB path.
% This function has been validated using the test data from NBS Monograph 85 85 % This function has been validated using the test data from NBS Monograph
% 140, the 1000-point test data set given by Riley [1], and the example data 86 86 % 140, the 1000-point test data set given by Riley [1], and the example data
% given in IEEE standard 1139-1999, Annex C. 87 87 % given in IEEE standard 1139-1999, Annex C.
% The author welcomes other validation results, see contact info below. 88 88 % The author welcomes other validation results, see contact info below.
% 89 89 %
% For more information, see: 90 90 % For more information, see:
% [1] W. J. Riley, "The Calculation of Time Domain Frequency Stability," 91 91 % [1] W. J. Riley, "The Calculation of Time Domain Frequency Stability,"
% Available on the web: 92 92 % Available on the web:
% http://www.ieee-uffc.org/frequency_control/teaching.asp?name=paper1ht 93 93 % http://www.ieee-uffc.org/frequency_control/teaching.asp?name=paper1ht
% 94 94 %
% 95 95 %
% M.A. Hopcroft 96 96 % M.A. Hopcroft
% mhopeng at gmail dot com 97 97 % mhopeng at gmail dot com
% 98 98 %
% I welcome your comments and feedback! 99 99 % I welcome your comments and feedback!
% 100 100 %
% MH Mar2014 101 101 % MH Mar2014
% v2.24 fix bug related to generating freq data from phase with timestamps 102 102 % v2.24 fix bug related to generating freq data from phase with timestamps
% (thanks to S. David-Grignot for finding the bug) 103 103 % (thanks to S. David-Grignot for finding the bug)
% MH Oct2010 104 104 % MH Oct2010
% v2.22 tau truncation to integer groups; tau sort 105 105 % v2.22 tau truncation to integer groups; tau sort
% plotting bugfix 106 106 % plotting bugfix
% v2.20 sychronize updates across allan, allan_overlap, allan_modified 107 107 % v2.20 sychronize updates across allan, allan_overlap, allan_modified
% v2.16 add TAU as output, fixed unusual error with dsplot v1.1 108 108 % v2.16 add TAU as output, fixed unusual error with dsplot v1.1
% v2.14 update plotting behaviour, default tau values 109 109 % v2.14 update plotting behaviour, default tau values
% 110 110 %
111 111
versionstr = 'allan v2.24'; 112 112 versionstr = 'allan v2.24';
113 113
% MH Jun2010 114 114 % MH Jun2010
% v2.12 bugfix for rate data row/col orientation 115 115 % v2.12 bugfix for rate data row/col orientation
% add DATA.units for plotting 116 116 % add DATA.units for plotting
% use dsplot.m for plotting 117 117 % use dsplot.m for plotting
% 118 118 %
% MH MAR2010 119 119 % MH MAR2010
% v2.1 minor interface and bugfixes 120 120 % v2.1 minor interface and bugfixes
% update data consistency check 121 121 % update data consistency check
% 122 122 %
% MH FEB2010 123 123 % MH FEB2010
% v2.0 Consistent code behaviour for all "allan_x.m" functions: 124 124 % v2.0 Consistent code behaviour for all "allan_x.m" functions:
% accept phase data 125 125 % accept phase data
% verbose levels 126 126 % verbose levels
% 127 127 %
% 128 128 %
% MH JAN2010 129 129 % MH JAN2010
% v1.84 code cleanup 130 130 % v1.84 code cleanup
% v1.82 typos in comments and code cleanup 131 131 % v1.82 typos in comments and code cleanup
% tau bin plotting changed for performance improvement 132 132 % tau bin plotting changed for performance improvement
% v1.8 Performance improvements: 133 133 % v1.8 Performance improvements:
% vectorize code for rate data 134 134 % vectorize code for rate data
% logical indexing for irregular rate data 135 135 % logical indexing for irregular rate data
% MH APR2008 136 136 % MH APR2008
% v1.62 loglog plot option 137 137 % v1.62 loglog plot option
% v1.61 improve error handling, plotting 138 138 % v1.61 improve error handling, plotting
% fix bug in regular data calc for high-rate data 139 139 % fix bug in regular data calc for high-rate data
% fix bug in timestamp data calc for large starting gap 140 140 % fix bug in timestamp data calc for large starting gap
% (thanks to C. B. Ruiz for identifying these bugs) 141 141 % (thanks to C. B. Ruiz for identifying these bugs)
% uses timestamps for DATA.rate=0 142 142 % uses timestamps for DATA.rate=0
% progress indicator for large timestamp data processing 143 143 % progress indicator for large timestamp data processing
% MH JUN2007 144 144 % MH JUN2007
% v1.54 Improve data plotting and optional bin plotting 145 145 % v1.54 Improve data plotting and optional bin plotting
% MH FEB2007 146 146 % MH FEB2007
% v1.5 use difference from median for plotting 147 147 % v1.5 use difference from median for plotting
% added MAD calculation for outlier detection 148 148 % added MAD calculation for outlier detection
% MH JAN2007 149 149 % MH JAN2007
% v1.48 plotting typos fixes 150 150 % v1.48 plotting typos fixes
% MH DEC2006 151 151 % MH DEC2006
% v1.46 hack to plot error bars 152 152 % v1.46 hack to plot error bars
% v1.44 further validation (Riley 1000-pt) 153 153 % v1.44 further validation (Riley 1000-pt)
% plot mean and std 154 154 % plot mean and std
% MH NOV2006 155 155 % MH NOV2006
% v1.42 typo fix comments 156 156 % v1.42 typo fix comments
% v1.4 fix irregular rate algorithm 157 157 % v1.4 fix irregular rate algorithm
% irregular algorithm rejects tau less than max gap in time data 158 158 % irregular algorithm rejects tau less than max gap in time data
% validate both algorithms using test data from NBS Monograph 140 159 159 % validate both algorithms using test data from NBS Monograph 140
% v1.3 fix time calc if data.time not present 160 160 % v1.3 fix time calc if data.time not present
% add error bars (not possible due to bug in MATLAB R14SP3) 161 161 % add error bars (not possible due to bug in MATLAB R14SP3)
% remove offset calculation 162 162 % remove offset calculation
% v1.24 improve feedback 163 163 % v1.24 improve feedback
% MH SEP2006 164 164 % MH SEP2006
% v1.22 updated comments 165 165 % v1.22 updated comments
% v1.2 errors and warnings 166 166 % v1.2 errors and warnings
% v1.1 handle irregular interval data 167 167 % v1.1 handle irregular interval data
%#ok<*AGROW> 168 168 %#ok<*AGROW>
169 169
% defaults 170 170 % defaults
if nargin < 4, verbose=2; end 171 171 if nargin < 4, verbose=2; end
if nargin < 3, name=''; end 172 172 if nargin < 3, name=''; end
if nargin < 2 || isempty(tau), tau=2.^(-10:10); end 173 173 if nargin < 2 || isempty(tau), tau=2.^(-10:10); end
174 174
% plot "tau bins"? #TAUBIN 175 175 % plot "tau bins"? #TAUBIN
TAUBIN = 0; % set 0 or 1 % WARNING: this has a significant impact on performance 176 176 TAUBIN = 0; % set 0 or 1 % WARNING: this has a significant impact on performance
177 177
% Formatting for plots 178 178 % Formatting for plots
FontName = 'Arial'; 179 179 FontName = 'Arial';
FontSize = 14; 180 180 FontSize = 14;
plotlinewidth=2; 181 181 plotlinewidth=2;
182 182
if verbose >= 1, fprintf(1,'allan: %s\n\n',versionstr); end 183 183 if verbose >= 1, fprintf(1,'allan: %s\n\n',versionstr); end
184 184
%% Data consistency checks 185 185 %% Data consistency checks
if ~(isfield(data,'phase') || isfield(data,'freq')) 186 186 if ~(isfield(data,'phase') || isfield(data,'freq'))
error('Either ''phase'' or ''freq'' must be present in DATA. See help file for details. [con0]'); 187 187 error('Either ''phase'' or ''freq'' must be present in DATA. See help file for details. [con0]');
end 188 188 end
if isfield(data,'time') 189 189 if isfield(data,'time')
if isfield(data,'phase') && (length(data.phase) ~= length(data.time)) 190 190 if isfield(data,'phase') && (length(data.phase) ~= length(data.time))
if isfield(data,'freq') && (length(data.freq) ~= length(data.time)) 191 191 if isfield(data,'freq') && (length(data.freq) ~= length(data.time))
error('The time and freq vectors are not the same length. See help for details. [con2]'); 192 192 error('The time and freq vectors are not the same length. See help for details. [con2]');
else 193 193 else
error('The time and phase vectors are not the same length. See help for details. [con1]'); 194 194 error('The time and phase vectors are not the same length. See help for details. [con1]');
end 195 195 end
end 196 196 end
if isfield(data,'phase') && (any(isnan(data.phase)) || any(isinf(data.phase))) 197 197 if isfield(data,'phase') && (any(isnan(data.phase)) || any(isinf(data.phase)))
error('The phase vector contains invalid elements (NaN/Inf). [con3]'); 198 198 error('The phase vector contains invalid elements (NaN/Inf). [con3]');
end 199 199 end
if isfield(data,'freq') && (any(isnan(data.freq)) || any(isinf(data.freq))) 200 200 if isfield(data,'freq') && (any(isnan(data.freq)) || any(isinf(data.freq)))
error('The freq vector contains invalid elements (NaN/Inf). [con4]'); 201 201 error('The freq vector contains invalid elements (NaN/Inf). [con4]');
end 202 202 end
if isfield(data,'time') && (any(isnan(data.time)) || any(isinf(data.time))) 203 203 if isfield(data,'time') && (any(isnan(data.time)) || any(isinf(data.time)))
error('The time vector contains invalid elements (NaN/Inf). [con5]'); 204 204 error('The time vector contains invalid elements (NaN/Inf). [con5]');
end 205 205 end
end 206 206 end
207 207
% sort tau vector 208 208 % sort tau vector
tau=sort(tau); 209 209 tau=sort(tau);
210 210
211 211
%% Basic statistical tests on the data set 212 212 %% Basic statistical tests on the data set
if ~isfield(data,'freq') 213 213 if ~isfield(data,'freq')
if isfield(data,'rate') && data.rate ~= 0 214 214 if isfield(data,'rate') && data.rate ~= 0
data.freq=diff(data.phase).*data.rate; 215 215 data.freq=diff(data.phase).*data.rate;
elseif isfield(data,'time') 216 216 elseif isfield(data,'time')
data.freq=diff(data.phase)./diff(data.time); 217 217 data.freq=diff(data.phase)./diff(data.time);
end 218 218 end
if verbose >= 1, fprintf(1,'allan: Fractional frequency data generated from phase data (M=%g).\n',length(data.freq)); end 219 219 if verbose >= 1, fprintf(1,'allan: Fractional frequency data generated from phase data (M=%g).\n',length(data.freq)); end
data.time(1)=[]; % make time stamps correspond to freq data 220 220 data.time(1)=[]; % make time stamps correspond to freq data
end 221 221 end
if size(data.freq,2) > size(data.freq,1), data.freq=data.freq'; end % ensure columns 222 222 if size(data.freq,2) > size(data.freq,1), data.freq=data.freq'; end % ensure columns
223 223
s.numpoints=length(data.freq); 224 224 s.numpoints=length(data.freq);
s.max=max(data.freq); 225 225 s.max=max(data.freq);
s.min=min(data.freq); 226 226 s.min=min(data.freq);
s.mean=mean(data.freq); 227 227 s.mean=mean(data.freq);
s.median=median(data.freq); 228 228 s.median=median(data.freq);
if isfield(data,'time') 229 229 if isfield(data,'time')
if size(data.time,2) > size(data.time,1), data.time=data.time'; end % ensure columns 230 230 if size(data.time,2) > size(data.time,1), data.time=data.time'; end % ensure columns
s.linear=polyfit(data.time(1:length(data.freq)),data.freq,1); 231 231 s.linear=polyfit(data.time(1:length(data.freq)),data.freq,1);
elseif isfield(data,'rate') && data.rate ~= 0; 232 232 elseif isfield(data,'rate') && data.rate ~= 0;
s.linear=polyfit((1/data.rate:1/data.rate:length(data.freq)/data.rate)',data.freq,1); 233 233 s.linear=polyfit((1/data.rate:1/data.rate:length(data.freq)/data.rate)',data.freq,1);
else 234 234 else
error('Either "time" or "rate" must be present in DATA. Type "help allan" for details. [err1]'); 235 235 error('Either "time" or "rate" must be present in DATA. Type "help allan" for details. [err1]');
end 236 236 end
s.std=std(data.freq); 237 237 s.std=std(data.freq);
238 238
if verbose >= 2 239 239 if verbose >= 2
fprintf(1,'allan: input data statistics:\n'); 240 240 fprintf(1,'allan: input data statistics:\n');
disp(s); 241 241 disp(s);
end 242 242 end
243 243
244 244
% center at median for plotting 245 245 % center at median for plotting
medianfreq=data.freq-s.median; 246 246 medianfreq=data.freq-s.median;
sm=[]; sme=[]; 247 247 sm=[]; sme=[];
248 248
% Screen for outliers using 5x Median Absolute Deviation (MAD) criteria 249 249 % Screen for outliers using 5x Median Absolute Deviation (MAD) criteria
s.MAD = median(abs(medianfreq)/0.6745); 250 250 s.MAD = median(abs(medianfreq)/0.6745);
if verbose >= 2 251 251 if verbose >= 2
fprintf(1, 'allan: 5x MAD value for outlier detection: %g\n',5*s.MAD); 252 252 fprintf(1, 'allan: 5x MAD value for outlier detection: %g\n',5*s.MAD);
end 253 253 end
if verbose >= 1 && any(abs(medianfreq) > 5*s.MAD) 254 254 if verbose >= 1 && any(abs(medianfreq) > 5*s.MAD)
fprintf(1, 'allan: NOTE: There appear to be outliers in the frequency data. See plot.\n'); 255 255 fprintf(1, 'allan: NOTE: There appear to be outliers in the frequency data. See plot.\n');
end 256 256 end
257 257
258 258
%%%% 259 259 %%%%
% There are two cases, either using timestamps or fixed sample rate: 260 260 % There are two cases, either using timestamps or fixed sample rate:
261 261
%% Fixed Sample Rate Data 262 262 %% Fixed Sample Rate Data
% If there is a regular interval between measurements, calculation is much 263 263 % If there is a regular interval between measurements, calculation is much
% easier/faster 264 264 % easier/faster
if isfield(data,'rate') && data.rate > 0 % if data rate was given 265 265 if isfield(data,'rate') && data.rate > 0 % if data rate was given
if verbose >= 1, fprintf(1, 'allan: regular data (%g data points @ %g Hz)\n',length(data.freq),data.rate); end 266 266 if verbose >= 1, fprintf(1, 'allan: regular data (%g data points @ %g Hz)\n',length(data.freq),data.rate); end
267 267
% string for plot title 268 268 % string for plot title
name=[name ' (' num2str(data.rate) ' Hz)']; 269 269 name=[name ' (' num2str(data.rate) ' Hz)'];
270 270
% what is the time interval between data points? 271 271 % what is the time interval between data points?
tmstep = 1/data.rate; 272 272 tmstep = 1/data.rate;
273 273
% Is there time data? Just for curiosity/plotting, does not impact calculation 274 274 % Is there time data? Just for curiosity/plotting, does not impact calculation
if isfield(data,'time') 275 275 if isfield(data,'time')
% adjust time data to remove any starting gap; first time step 276 276 % adjust time data to remove any starting gap; first time step
% should not be zero for comparison with freq data 277 277 % should not be zero for comparison with freq data
dtime=data.time-data.time(1)+mean(diff(data.time)); 278 278 dtime=data.time-data.time(1)+mean(diff(data.time));
if verbose >= 2 279 279 if verbose >= 2
fprintf(1,'allan: End of timestamp data: %g sec.\n',dtime(end)); 280 280 fprintf(1,'allan: End of timestamp data: %g sec.\n',dtime(end));
if (data.rate - 1/mean(diff(dtime))) > 1e-6 281 281 if (data.rate - 1/mean(diff(dtime))) > 1e-6
fprintf(1,'allan: NOTE: data.rate (%f Hz) does not match average timestamped sample rate (%f Hz)\n',data.rate,1/mean(diff(dtime))); 282 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)));
end 283 283 end
end 284 284 end
else 285 285 else
% create time axis data using rate (for plotting only) 286 286 % create time axis data using rate (for plotting only)
dtime=(tmstep:tmstep:length(data.freq)*tmstep)'; % column oriented 287 287 dtime=(tmstep:tmstep:length(data.freq)*tmstep)'; % column oriented
end 288 288 end
289 289
% check the range of tau values and truncate if necessary 290 290 % check the range of tau values and truncate if necessary
% find halfway point of time record 291 291 % find halfway point of time record
halftime = round(tmstep*length(data.freq)/2); 292 292 halftime = round(tmstep*length(data.freq)/2);
% truncate tau to appropriate values 293 293 % truncate tau to appropriate values
tau = tau(tau >= tmstep & tau <= halftime); 294 294 tau = tau(tau >= tmstep & tau <= halftime);
if verbose >= 2, fprintf(1, 'allan: allowable tau range: %g to %g sec. (1/rate to total_time/2)\n',tmstep,halftime); end 295 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 296
% save the freq data for the loop 297 297 % save the freq data for the loop
dfreq=data.freq; 298 298 dfreq=data.freq;
dfreq2=data.freq2; 299 299 dfreq2=data.freq2;
% find the number of data points in each tau group 300 300 % find the number of data points in each tau group
m = data.rate.*tau; 301 301 m = data.rate.*tau;
% only integer values allowed (no fractional groups of points) 302 302 % only integer values allowed (no fractional groups of points)
%tau = tau(m-round(m)<1e-8); % numerical precision issues (v2.1) 303 303 %tau = tau(m-round(m)<1e-8); % numerical precision issues (v2.1)
tau = tau(m==round(m)); % The round() test is only correct for values < 2^53 304 304 tau = tau(m==round(m)); % The round() test is only correct for values < 2^53
%m = m(m-round(m)<1e-8); % change to round(m) for integer test v2.22 305 305 %m = m(m-round(m)<1e-8); % change to round(m) for integer test v2.22
m = m(m==round(m)); 306 306 m = m(m==round(m));
%m=round(m); 307 307 %m=round(m);
308 308
if verbose >= 1, fprintf(1,'allan: calculating Allan deviation...\n '); end 309 309 if verbose >= 1, fprintf(1,'allan: calculating Allan deviation...\n '); end
310 310
% calculate the Allan deviation for each value of tau 311 311 % calculate the Allan deviation for each value of tau
k=0; tic; 312 312 k=0; tic;
for i = tau 313 313 for i = tau
if verbose >= 2, fprintf(1,'%g ',i); end 314 314 if verbose >= 2, fprintf(1,'%g ',i); end
k=k+1; 315 315 k=k+1;
316 316
% truncate frequency set to an even multiple of this tau value 317 317 % truncate frequency set to an even multiple of this tau value
freq=dfreq(1:end-rem(length(dfreq),m(k))); 318 318 freq=dfreq(1:end-rem(length(dfreq),m(k)));
freq2=dfreq2(1:end-rem(length(dfreq2),m(k))); 319 319 freq2=dfreq2(1:end-rem(length(dfreq2),m(k)));
% group the data into tau-length groups or bins 320 320 % group the data into tau-length groups or bins
f = reshape(freq,m(k),[]); % Vectorize! 321 321 f = reshape(freq,m(k),[]); % Vectorize!
f2 = reshape(freq2,m(k),[]); % Vectorize! 322 322 f2 = reshape(freq2,m(k),[]); % Vectorize!
% find average in each "tau group", y_k (each colummn of f) 323 323 % find average in each "tau group", y_k (each colummn of f)
fa=mean(f,1); 324 324 fa=mean(f,1);
fa2=mean(f2,1); 325 325 fa2=mean(f2,1);
% first finite difference 326 326 % first finite difference
fd=diff(fa); 327 327 fd=diff(fa);
fd2=diff(fa2); 328 328 fd2=diff(fa2);
% calculate two-sample variance for this tau 329 329 % calculate two-sample variance for this tau
M=length(fa); 330 330 M=length(fa);
sm(k)=sqrt(0.5/(M-1)*(real(sum(fd.*fd2)))); 331 331 sm(k)=sqrt(0.5/(M-1)*abs(real(sum(fd.*fd2))));
332 332
% estimate error bars 333 333 % estimate error bars
sme(k)=sm(k)/sqrt(M+1); 334 334 sme(k)=sm(k)/sqrt(M+1);
335 335
if TAUBIN == 1 336 336 if TAUBIN == 1
% save the binning points for plotting 337 337 % save the binning points for plotting
fs(k,1:length(freq)/m(k))=m(k):m(k):length(freq); fval{k}=mean(f,1); 338 338 fs(k,1:length(freq)/m(k))=m(k):m(k):length(freq); fval{k}=mean(f,1);
end 339 339 end
340 340
end % repeat for each value of tau 341 341 end % repeat for each value of tau
342 342
if verbose >= 2, fprintf(1,'\n'); end 343 343 if verbose >= 2, fprintf(1,'\n'); end
calctime=toc; if verbose >= 2, fprintf(1,'allan: Elapsed time for calculation: %e seconds\n',calctime); end 344 344 calctime=toc; if verbose >= 2, fprintf(1,'allan: Elapsed time for calculation: %e seconds\n',calctime); end
345 345
346 346
347 347
%% Irregular data (timestamp) 348 348 %% Irregular data (timestamp)
elseif isfield(data,'time') 349 349 elseif isfield(data,'time')
% the interval between measurements is irregular 350 350 % the interval between measurements is irregular
% so we must group the data by time 351 351 % so we must group the data by time
if verbose >= 1, fprintf(1, 'allan: irregular rate data (no fixed sample rate)\n'); end 352 352 if verbose >= 1, fprintf(1, 'allan: irregular rate data (no fixed sample rate)\n'); end
353 353
% string for plot title 354 354 % string for plot title
name=[name ' (timestamp)']; 355 355 name=[name ' (timestamp)'];
356 356
% adjust time to remove any initial offset or zero 357 357 % adjust time to remove any initial offset or zero
dtime=data.time-data.time(1)+mean(diff(data.time)); 358 358 dtime=data.time-data.time(1)+mean(diff(data.time));
%dtime=data.time; 359 359 %dtime=data.time;
% where is the maximum gap in time record? 360 360 % where is the maximum gap in time record?
gap_pos=find(diff(dtime)==max(diff(dtime))); 361 361 gap_pos=find(diff(dtime)==max(diff(dtime)));
% what is average data spacing? 362 362 % what is average data spacing?
avg_gap = mean(diff(dtime)); 363 363 avg_gap = mean(diff(dtime));
364 364
if verbose >= 2 365 365 if verbose >= 2
fprintf(1, 'allan: WARNING: irregular timestamp data (no fixed sample rate).\n'); 366 366 fprintf(1, 'allan: WARNING: irregular timestamp data (no fixed sample rate).\n');
fprintf(1, ' Calculation time may be long and the results subject to interpretation.\n'); 367 367 fprintf(1, ' Calculation time may be long and the results subject to interpretation.\n');
fprintf(1, ' You are advised to estimate using an average sample rate (%g Hz) instead of timestamps.\n',1/avg_gap); 368 368 fprintf(1, ' You are advised to estimate using an average sample rate (%g Hz) instead of timestamps.\n',1/avg_gap);
fprintf(1, ' Continue at your own risk! (press any key to continue)\n'); 369 369 fprintf(1, ' Continue at your own risk! (press any key to continue)\n');
pause; 370 370 pause;
end 371 371 end
372 372
if verbose >= 1 373 373 if verbose >= 1
fprintf(1, 'allan: End of timestamp data: %g sec\n',dtime(end)); 374 374 fprintf(1, 'allan: End of timestamp data: %g sec\n',dtime(end));
fprintf(1, ' Average rate: %g Hz (%g sec/measurement)\n',1/avg_gap,avg_gap); 375 375 fprintf(1, ' Average rate: %g Hz (%g sec/measurement)\n',1/avg_gap,avg_gap);
if max(diff(dtime)) ~= 1/mean(diff(dtime)) 376 376 if max(diff(dtime)) ~= 1/mean(diff(dtime))
fprintf(1, ' Max. gap: %g sec at position %d\n',max(diff(dtime)),gap_pos(1)); 377 377 fprintf(1, ' Max. gap: %g sec at position %d\n',max(diff(dtime)),gap_pos(1));
end 378 378 end
if max(diff(dtime)) > 5*avg_gap 379 379 if max(diff(dtime)) > 5*avg_gap
fprintf(1, ' WARNING: Max. gap in time record is suspiciously large (>5x the average interval).\n'); 380 380 fprintf(1, ' WARNING: Max. gap in time record is suspiciously large (>5x the average interval).\n');
end 381 381 end
end 382 382 end
383 383
384 384
% find halfway point 385 385 % find halfway point
halftime = fix(dtime(end)/2); 386 386 halftime = fix(dtime(end)/2);
% truncate tau to appropriate values 387 387 % truncate tau to appropriate values
tau = tau(tau >= max(diff(dtime)) & tau <= halftime); 388 388 tau = tau(tau >= max(diff(dtime)) & tau <= halftime);
if isempty(tau) 389 389 if isempty(tau)
error('allan: ERROR: no appropriate tau values (> %g s, < %g s)\n',max(diff(dtime)),halftime); 390 390 error('allan: ERROR: no appropriate tau values (> %g s, < %g s)\n',max(diff(dtime)),halftime);
end 391 391 end
392 392
% save the freq data for the loop 393 393 % save the freq data for the loop
dfreq=data.freq; 394 394 dfreq=data.freq;
dtime=dtime(1:length(dfreq)); 395 395 dtime=dtime(1:length(dfreq));
396 396
if verbose >= 1, fprintf(1,'allan: calculating Allan deviation...\n'); end 397 397 if verbose >= 1, fprintf(1,'allan: calculating Allan deviation...\n'); end
398 398
k=0; tic; 399 399 k=0; tic;
for i = tau 400 400 for i = tau
if verbose >= 2, fprintf(1,'%d ',i); end 401 401 if verbose >= 2, fprintf(1,'%d ',i); end
402 402
k=k+1; fa=[]; %f=[]; 403 403 k=k+1; fa=[]; %f=[];
km=0; 404 404 km=0;
405 405
% truncate data set to an even multiple of this tau value 406 406 % truncate data set to an even multiple of this tau value
freq=dfreq(dtime <= dtime(end)-rem(dtime(end),i)); 407 407 freq=dfreq(dtime <= dtime(end)-rem(dtime(end),i));
time=dtime(dtime <= dtime(end)-rem(dtime(end),i)); 408 408 time=dtime(dtime <= dtime(end)-rem(dtime(end),i));
%freq=dfreq; 409 409 %freq=dfreq;
%time=dtime; 410 410 %time=dtime;
411 411
% break up the data into groups of tau length in sec 412 412 % break up the data into groups of tau length in sec
while i*km < time(end) 413 413 while i*km < time(end)
km=km+1; 414 414 km=km+1;
415 415
% progress bar 416 416 % progress bar
if verbose >= 2 417 417 if verbose >= 2
if rem(km,100)==0, fprintf(1,'.'); end 418 418 if rem(km,100)==0, fprintf(1,'.'); end
if rem(km,1000)==0, fprintf(1,'%g/%g\n',km,round(time(end)/i)); end 419 419 if rem(km,1000)==0, fprintf(1,'%g/%g\n',km,round(time(end)/i)); end
end 420 420 end
421 421
f = freq(i*(km-1) < time & time <= i*km); 422 422 f = freq(i*(km-1) < time & time <= i*km);
f = f(~isnan(f)); % make sure values are valid 423 423 f = f(~isnan(f)); % make sure values are valid
424 424
if ~isempty(f) 425 425 if ~isempty(f)
fa(km)=mean(f); 426 426 fa(km)=mean(f);
else 427 427 else
fa(km)=0; 428 428 fa(km)=0;
end 429 429 end
430 430
if TAUBIN == 1 % WARNING: this has a significant impact on performance 431 431 if TAUBIN == 1 % WARNING: this has a significant impact on performance
% save the binning points for plotting 432 432 % save the binning points for plotting
%if find(time <= i*km) > 0 433 433 %if find(time <= i*km) > 0
fs(k,km)=max(time(time <= i*km)); 434 434 fs(k,km)=max(time(time <= i*km));
%else 435 435 %else
if isempty(fs(k,km)) 436 436 if isempty(fs(k,km))
fs(k,km)=0; 437 437 fs(k,km)=0;
end 438 438 end
fval{k}=fa; 439 439 fval{k}=fa;
end % save tau bin plot points 440 440 end % save tau bin plot points
441 441
end 442 442 end
443 443
if verbose >= 2, fprintf(1,'\n'); end 444 444 if verbose >= 2, fprintf(1,'\n'); end
445 445
% first finite difference of the averaged results 446 446 % first finite difference of the averaged results
fd=diff(fa); 447 447 fd=diff(fa);
% calculate Allan deviation for this tau 448 448 % calculate Allan deviation for this tau
M=length(fa); 449 449 M=length(fa);
sm(k)=sqrt(0.5/(M-1)*(sum(fd.^2))); 450 450 sm(k)=sqrt(0.5/(M-1)*(sum(fd.^2)));
451 451
% estimate error bars 452 452 % estimate error bars
sme(k)=sm(k)/sqrt(M+1); 453 453 sme(k)=sm(k)/sqrt(M+1);
454 454
455 455
end 456 456 end
457 457
if verbose == 2, fprintf(1,'\n'); end 458 458 if verbose == 2, fprintf(1,'\n'); end
calctime=toc; if verbose >= 2, fprintf(1,'allan: Elapsed time for calculation: %e seconds\n',calctime); end 459 459 calctime=toc; if verbose >= 2, fprintf(1,'allan: Elapsed time for calculation: %e seconds\n',calctime); end
460 460
461 461
else 462 462 else
error('allan: WARNING: no DATA.rate or DATA.time! Type "help allan" for more information. [err2]'); 463 463 error('allan: WARNING: no DATA.rate or DATA.time! Type "help allan" for more information. [err2]');
end 464 464 end
465 465
466 466
%%%%%%%% 467 467 %%%%%%%%
%% Plotting 468 468 %% Plotting
469 469
if verbose >= 2 % show all data 470 470 if verbose >= 2 % show all data
471 471
% plot the frequency data, centered on median 472 472 % plot the frequency data, centered on median
if size(dtime,2) > size(dtime,1), dtime=dtime'; end % this should not be necessary, but dsplot 1.1 is a little bit brittle 473 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
try 474 474 try
% dsplot makes a new figure 475 475 % dsplot makes a new figure
hd=dsplot(dtime,medianfreq); 476 476 hd=dsplot(dtime,medianfreq);
catch ME 477 477 catch ME
figure; 478 478 figure;
if length(dtime) ~= length(medianfreq) 479 479 if length(dtime) ~= length(medianfreq)
fprintf(1,'allan: Warning: length of time axis (%d) is not equal to data array (%d)\n',length(dtime),length(medianfreq)); 480 480 fprintf(1,'allan: Warning: length of time axis (%d) is not equal to data array (%d)\n',length(dtime),length(medianfreq));
end 481 481 end
hd=plot(dtime,medianfreq); 482 482 hd=plot(dtime,medianfreq);
if verbose >= 1, fprintf(1,'allan: Note: Install dsplot.m for improved plotting of large data sets (File Exchange File ID: #15850).\n'); end 483 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
if verbose >= 2, fprintf(1,' (Message: %s)\n',ME.message); end 484 484 if verbose >= 2, fprintf(1,' (Message: %s)\n',ME.message); end
end 485 485 end
set(hd,'Marker','.','LineStyle','none','Color','b'); % equivalent to '.-' 486 486 set(hd,'Marker','.','LineStyle','none','Color','b'); % equivalent to '.-'
hold on; 487 487 hold on;
488 488
% show center (0) 489 489 % show center (0)
plot(xlim,[0 0],':k'); 490 490 plot(xlim,[0 0],':k');
% show 5x Median Absolute Deviation (MAD) values 491 491 % show 5x Median Absolute Deviation (MAD) values
hm=plot(xlim,[5*s.MAD 5*s.MAD],'-r'); 492 492 hm=plot(xlim,[5*s.MAD 5*s.MAD],'-r');
plot(xlim,[-5*s.MAD -5*s.MAD],'-r'); 493 493 plot(xlim,[-5*s.MAD -5*s.MAD],'-r');
% show linear fit line 494 494 % show linear fit line
hf=plot(xlim,polyval(s.linear,xlim)-s.median,'-g'); 495 495 hf=plot(xlim,polyval(s.linear,xlim)-s.median,'-g');
title(['Data: ' name],'FontSize',FontSize+2,'FontName',FontName); 496 496 title(['Data: ' name],'FontSize',FontSize+2,'FontName',FontName);
%set(get(gca,'Title'),'Interpreter','none'); 497 497 %set(get(gca,'Title'),'Interpreter','none');
xlabel('Time [sec]','FontSize',FontSize,'FontName',FontName); 498 498 xlabel('Time [sec]','FontSize',FontSize,'FontName',FontName);
if isfield(data,'units') 499 499 if isfield(data,'units')
ylabel(['data - median(data) [' data.units ']'],'FontSize',FontSize,'FontName',FontName); 500 500 ylabel(['data - median(data) [' data.units ']'],'FontSize',FontSize,'FontName',FontName);
else 501 501 else
ylabel('freq - median(freq)','FontSize',FontSize,'FontName',FontName); 502 502 ylabel('freq - median(freq)','FontSize',FontSize,'FontName',FontName);
end 503 503 end
set(gca,'FontSize',FontSize,'FontName',FontName); 504 504 set(gca,'FontSize',FontSize,'FontName',FontName);
legend([hd hm hf],{'data (centered on median)','5x MAD outliers',['Linear Fit (' num2str(s.linear(1),'%g') ')']},'FontSize',max(10,FontSize-2)); 505 505 legend([hd hm hf],{'data (centered on median)','5x MAD outliers',['Linear Fit (' num2str(s.linear(1),'%g') ')']},'FontSize',max(10,FontSize-2));