bmarechal / octave_scripts

Download as
Browse Code ยป

Commit 1a0e88f0ce6e736bf1f59e7b02cc394896768bff

Authored by bmarechal
1 parent 183c460e4e
Exists in master

replace 4-spaces by tab

Showing 8 changed files with 1688 additions and 1680 deletions Inline Diff

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