diff --git a/allan.m b/allan.m new file mode 100644 index 0000000..43f4348 --- /dev/null +++ b/allan.m @@ -0,0 +1,576 @@ +function [retval, s, errorb, tau] = allan(data,tau,name,verbose) +% ALLAN Compute the Allan deviation for a set of time-domain frequency data +% [RETVAL, S, ERRORB, TAU] = ALLAN(DATA,TAU,NAME,VERBOSE) +% +% Inputs: +% DATA should be a structure and have the following fields: +% DATA.freq or DATA.phase +% A vector of fractional frequency measurements (df/f) in +% DATA.freq *or* phase offset data (seconds) in DATA.phase . +% If frequency data is not present, it will be generated by +% differentiating the phase data. +% If both fields are present, then DATA.freq will be used. +% Note: for general-purpose calculations of Allan deviation, +% (i.e. a two-sample variance) use DATA.freq . +% +% DATA.rate or DATA.time +% The sampling rate in Hertz (DATA.rate) or a vector of +% timestamps for each measurement in seconds (DATA.time). +% DATA.rate is used if both fields are present. +% If DATA.rate == 0, then the timestamps are used. +% +% DATA.units (optional) +% The units for the data. If present, the string DATA.units +% is added to the plot y-axis label. +% +% TAU is an array of tau values for computing Allan deviation. +% TAU values must be divisible by 1/DATA.rate (data points cannot be +% grouped in fractional quantities!) and invalid values are ignored. +% Leave empty to use default values. +% NAME is an optional label that is added to the plot titles. +% VERBOSE sets the level of status messages: +% 0 = silent & no data plots; +% 1 = status messages & minimum plots; +% 2 = all messages and plots (default) +% +% Outputs: +% 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). +% 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. +% 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). +% +% Example: +% +% To compute the Allan deviation for the data in the variable "lt": +% >> lt +% lt = +% freq: [1x86400 double] +% rate: 0.5 +% +% Use: +% +% >> ad = allan(lt,[2 10 100],'lt data',1); +% +% 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. +% You can also use the default settings, which are usually a good starting point: +% +% >> ad = allan(lt); +% +% +% Notes: +% This function calculates the standard Allan deviation (ADEV), *not* the +% overlapping ADEV. Use "allan_overlap.m" for overlapping ADEV. +% The calculation is performed using fractional frequency data. If only +% phase data is provided, frequency data is generated by differentiating +% the phase data. +% No pre-processing of the data is performed, except to remove any +% initial offset (i.e., starting gap) in the time record. +% 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 +% time-stamped data, only tau values greater than the maximum gap between +% samples and less than half the total time are used. +% 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, +% 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 +% TAUBIN to 1 (search for "#TAUBIN"). +% You can choose between loglog and semilog plotting of results by +% commenting in/out the appropriate line. Search for "#PLOTLOG". +% I recommend installing "dsplot.m", which improves the performance of +% 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. +% 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 +% given in IEEE standard 1139-1999, Annex C. +% The author welcomes other validation results, see contact info below. +% +% For more information, see: +% [1] W. J. Riley, "The Calculation of Time Domain Frequency Stability," +% Available on the web: +% http://www.ieee-uffc.org/frequency_control/teaching.asp?name=paper1ht +% +% +% M.A. Hopcroft +% mhopeng at gmail dot com +% +% I welcome your comments and feedback! +% +% MH Mar2014 +% v2.24 fix bug related to generating freq data from phase with timestamps +% (thanks to S. David-Grignot for finding the bug) +% MH Oct2010 +% v2.22 tau truncation to integer groups; tau sort +% plotting bugfix +% v2.20 sychronize updates across allan, allan_overlap, allan_modified +% v2.16 add TAU as output, fixed unusual error with dsplot v1.1 +% v2.14 update plotting behaviour, default tau values +% + +versionstr = 'allan v2.24'; + +% MH Jun2010 +% v2.12 bugfix for rate data row/col orientation +% add DATA.units for plotting +% use dsplot.m for plotting +% +% MH MAR2010 +% v2.1 minor interface and bugfixes +% update data consistency check +% +% MH FEB2010 +% v2.0 Consistent code behaviour for all "allan_x.m" functions: +% accept phase data +% verbose levels +% +% +% MH JAN2010 +% v1.84 code cleanup +% v1.82 typos in comments and code cleanup +% tau bin plotting changed for performance improvement +% v1.8 Performance improvements: +% vectorize code for rate data +% logical indexing for irregular rate data +% MH APR2008 +% v1.62 loglog plot option +% v1.61 improve error handling, plotting +% fix bug in regular data calc for high-rate data +% fix bug in timestamp data calc for large starting gap +% (thanks to C. B. Ruiz for identifying these bugs) +% uses timestamps for DATA.rate=0 +% progress indicator for large timestamp data processing +% MH JUN2007 +% v1.54 Improve data plotting and optional bin plotting +% MH FEB2007 +% v1.5 use difference from median for plotting +% added MAD calculation for outlier detection +% MH JAN2007 +% v1.48 plotting typos fixes +% MH DEC2006 +% v1.46 hack to plot error bars +% v1.44 further validation (Riley 1000-pt) +% plot mean and std +% MH NOV2006 +% v1.42 typo fix comments +% v1.4 fix irregular rate algorithm +% irregular algorithm rejects tau less than max gap in time data +% validate both algorithms using test data from NBS Monograph 140 +% v1.3 fix time calc if data.time not present +% add error bars (not possible due to bug in MATLAB R14SP3) +% remove offset calculation +% v1.24 improve feedback +% MH SEP2006 +% v1.22 updated comments +% v1.2 errors and warnings +% v1.1 handle irregular interval data +%#ok<*AGROW> + +% defaults +if nargin < 4, verbose=2; end +if nargin < 3, name=''; end +if nargin < 2 || isempty(tau), tau=2.^(-10:10); end + +% plot "tau bins"? #TAUBIN +TAUBIN = 0; % set 0 or 1 % WARNING: this has a significant impact on performance + +% Formatting for plots +FontName = 'Arial'; +FontSize = 14; +plotlinewidth=2; + +if verbose >= 1, fprintf(1,'allan: %s\n\n',versionstr); end + +%% Data consistency checks +if ~(isfield(data,'phase') || isfield(data,'freq')) + error('Either ''phase'' or ''freq'' must be present in DATA. See help file for details. [con0]'); +end +if isfield(data,'time') + if isfield(data,'phase') && (length(data.phase) ~= length(data.time)) + 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]'); + else + error('The time and phase vectors are not the same length. See help for details. [con1]'); + end + end + if isfield(data,'phase') && (any(isnan(data.phase)) || any(isinf(data.phase))) + error('The phase vector contains invalid elements (NaN/Inf). [con3]'); + end + if isfield(data,'freq') && (any(isnan(data.freq)) || any(isinf(data.freq))) + error('The freq vector contains invalid elements (NaN/Inf). [con4]'); + end + if isfield(data,'time') && (any(isnan(data.time)) || any(isinf(data.time))) + error('The time vector contains invalid elements (NaN/Inf). [con5]'); + end +end + +% sort tau vector +tau=sort(tau); + + +%% Basic statistical tests on the data set +if ~isfield(data,'freq') + if isfield(data,'rate') && data.rate ~= 0 + data.freq=diff(data.phase).*data.rate; + elseif isfield(data,'time') + data.freq=diff(data.phase)./diff(data.time); + end + 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 +end +if size(data.freq,2) > size(data.freq,1), data.freq=data.freq'; end % ensure columns + +s.numpoints=length(data.freq); +s.max=max(data.freq); +s.min=min(data.freq); +s.mean=mean(data.freq); +s.median=median(data.freq); +if isfield(data,'time') + 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); +elseif isfield(data,'rate') && data.rate ~= 0; + s.linear=polyfit((1/data.rate:1/data.rate:length(data.freq)/data.rate)',data.freq,1); +else + error('Either "time" or "rate" must be present in DATA. Type "help allan" for details. [err1]'); +end +s.std=std(data.freq); + +if verbose >= 2 + fprintf(1,'allan: input data statistics:\n'); + disp(s); +end + + +% center at median for plotting +medianfreq=data.freq-s.median; +sm=[]; sme=[]; + +% Screen for outliers using 5x Median Absolute Deviation (MAD) criteria +s.MAD = median(abs(medianfreq)/0.6745); +if verbose >= 2 + fprintf(1, 'allan: 5x MAD value for outlier detection: %g\n',5*s.MAD); +end +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'); +end + + +%%%% +% There are two cases, either using timestamps or fixed sample rate: + +%% Fixed Sample Rate Data +% If there is a regular interval between measurements, calculation is much +% easier/faster +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 + + % string for plot title + name=[name ' (' num2str(data.rate) ' Hz)']; + + % what is the time interval between data points? + tmstep = 1/data.rate; + + % Is there time data? Just for curiosity/plotting, does not impact calculation + if isfield(data,'time') + % adjust time data to remove any starting gap; first time step + % should not be zero for comparison with freq data + dtime=data.time-data.time(1)+mean(diff(data.time)); + if verbose >= 2 + fprintf(1,'allan: End of timestamp data: %g sec.\n',dtime(end)); + 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))); + end + end + else + % create time axis data using rate (for plotting only) + dtime=(tmstep:tmstep:length(data.freq)*tmstep)'; % column oriented + end + + % check the range of tau values and truncate if necessary + % find halfway point of time record + halftime = round(tmstep*length(data.freq)/2); + % truncate tau to appropriate values + 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 + + % save the freq data for the loop + dfreq=data.freq; + % find the number of data points in each tau group + m = data.rate.*tau; + % only integer values allowed (no fractional groups of points) + %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 + %m = m(m-round(m)<1e-8); % change to round(m) for integer test v2.22 + m = m(m==round(m)); + %m=round(m); + + if verbose >= 1, fprintf(1,'allan: calculating Allan deviation...\n '); end + + % calculate the Allan deviation for each value of tau + k=0; tic; + for i = tau + if verbose >= 2, fprintf(1,'%g ',i); end + k=k+1; + + % truncate frequency set to an even multiple of this tau value + freq=dfreq(1:end-rem(length(dfreq),m(k))); + % group the data into tau-length groups or bins + f = reshape(freq,m(k),[]); % Vectorize! + % find average in each "tau group", y_k (each colummn of f) + fa=mean(f,1); + % first finite difference + fd=diff(fa); + % calculate two-sample variance for this tau + M=length(fa); + sm(k)=sqrt(0.5/(M-1)*(sum(fd.^2))); + + % estimate error bars + sme(k)=sm(k)/sqrt(M+1); + + if TAUBIN == 1 + % save the binning points for plotting + fs(k,1:length(freq)/m(k))=m(k):m(k):length(freq); fval{k}=mean(f,1); + end + + end % repeat for each value of tau + + if verbose >= 2, fprintf(1,'\n'); end + calctime=toc; if verbose >= 2, fprintf(1,'allan: Elapsed time for calculation: %e seconds\n',calctime); end + + + +%% Irregular data (timestamp) +elseif isfield(data,'time') + % the interval between measurements is irregular + % so we must group the data by time + if verbose >= 1, fprintf(1, 'allan: irregular rate data (no fixed sample rate)\n'); end + + % string for plot title + name=[name ' (timestamp)']; + + % adjust time to remove any initial offset or zero + dtime=data.time-data.time(1)+mean(diff(data.time)); + %dtime=data.time; + % where is the maximum gap in time record? + gap_pos=find(diff(dtime)==max(diff(dtime))); + % what is average data spacing? + avg_gap = mean(diff(dtime)); + + if verbose >= 2 + 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'); + 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'); + pause; + end + + if verbose >= 1 + 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); + if max(diff(dtime)) ~= 1/mean(diff(dtime)) + fprintf(1, ' Max. gap: %g sec at position %d\n',max(diff(dtime)),gap_pos(1)); + end + if max(diff(dtime)) > 5*avg_gap + fprintf(1, ' WARNING: Max. gap in time record is suspiciously large (>5x the average interval).\n'); + end + end + + + % find halfway point + halftime = fix(dtime(end)/2); + % truncate tau to appropriate values + tau = tau(tau >= max(diff(dtime)) & tau <= halftime); + if isempty(tau) + error('allan: ERROR: no appropriate tau values (> %g s, < %g s)\n',max(diff(dtime)),halftime); + end + + % save the freq data for the loop + dfreq=data.freq; + dtime=dtime(1:length(dfreq)); + + if verbose >= 1, fprintf(1,'allan: calculating Allan deviation...\n'); end + + k=0; tic; + for i = tau + if verbose >= 2, fprintf(1,'%d ',i); end + + k=k+1; fa=[]; %f=[]; + km=0; + + % truncate data set to an even multiple of this tau value + freq=dfreq(dtime <= dtime(end)-rem(dtime(end),i)); + time=dtime(dtime <= dtime(end)-rem(dtime(end),i)); + %freq=dfreq; + %time=dtime; + + % break up the data into groups of tau length in sec + while i*km < time(end) + km=km+1; + + % progress bar + if verbose >= 2 + if rem(km,100)==0, fprintf(1,'.'); end + if rem(km,1000)==0, fprintf(1,'%g/%g\n',km,round(time(end)/i)); end + end + + f = freq(i*(km-1) < time & time <= i*km); + f = f(~isnan(f)); % make sure values are valid + + if ~isempty(f) + fa(km)=mean(f); + else + fa(km)=0; + end + + if TAUBIN == 1 % WARNING: this has a significant impact on performance + % save the binning points for plotting + %if find(time <= i*km) > 0 + fs(k,km)=max(time(time <= i*km)); + %else + if isempty(fs(k,km)) + fs(k,km)=0; + end + fval{k}=fa; + end % save tau bin plot points + + end + + if verbose >= 2, fprintf(1,'\n'); end + + % first finite difference of the averaged results + fd=diff(fa); + % calculate Allan deviation for this tau + M=length(fa); + sm(k)=sqrt(0.5/(M-1)*(sum(fd.^2))); + + % estimate error bars + sme(k)=sm(k)/sqrt(M+1); + + + end + + if verbose == 2, fprintf(1,'\n'); end + calctime=toc; if verbose >= 2, fprintf(1,'allan: Elapsed time for calculation: %e seconds\n',calctime); end + + +else + error('allan: WARNING: no DATA.rate or DATA.time! Type "help allan" for more information. [err2]'); +end + + +%%%%%%%% +%% Plotting + +if verbose >= 2 % show all data + + % 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 + try + % dsplot makes a new figure + hd=dsplot(dtime,medianfreq); + catch ME + figure; + 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)); + end + 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 + if verbose >= 2, fprintf(1,' (Message: %s)\n',ME.message); end + end + set(hd,'Marker','.','LineStyle','none','Color','b'); % equivalent to '.-' + hold on; + + % show center (0) + plot(xlim,[0 0],':k'); + % show 5x Median Absolute Deviation (MAD) values + hm=plot(xlim,[5*s.MAD 5*s.MAD],'-r'); + plot(xlim,[-5*s.MAD -5*s.MAD],'-r'); + % show linear fit line + hf=plot(xlim,polyval(s.linear,xlim)-s.median,'-g'); + title(['Data: ' name],'FontSize',FontSize+2,'FontName',FontName); + %set(get(gca,'Title'),'Interpreter','none'); + xlabel('Time [sec]','FontSize',FontSize,'FontName',FontName); + if isfield(data,'units') + ylabel(['data - median(data) [' data.units ']'],'FontSize',FontSize,'FontName',FontName); + else + ylabel('freq - median(freq)','FontSize',FontSize,'FontName',FontName); + end + 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)); + % tighten up + xlim([dtime(1) dtime(end)]); + + + % Optional tau bin (y_k samples) plot + if TAUBIN == 1 + % plot the tau divisions on the data plot + rfs=size(fs,1); + colororder=get(gca,'ColorOrder'); + axis tight; kc=2; + %ap=axis; + for j=1:rfs + 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 + for b=1:find(fs(j,:), 1, 'last' ); + % plot the tau division boundaries + %plot([fs(j,b) fs(j,b)],[ap(3)*1.1 ap(4)*1.1],'-','Color',colororder(kc,:)); + % plot tau group y values + if b == 1 + plot([dtime(1) fs(j,b)],[fval{j}(b)-s.median fval{j}(b)-s.median],'-','Color',colororder(kc,:),'LineWidth',4); + else + plot([fs(j,b-1) fs(j,b)],[fval{j}(b)-s.median fval{j}(b)-s.median],'-','Color',colororder(kc,:),'LineWidth',4); + end + end + end + axis auto + end % End optional bin plot + +end % end plot raw data + + +if verbose >= 1 % show ADEV results + + % plot Allan deviation results + if ~isempty(sm) + figure + + % Choose loglog or semilogx plot here #PLOTLOG + %semilogx(tau,sm,'.-b','LineWidth',plotlinewidth,'MarkerSize',24); + loglog(tau,sm,'.-b','LineWidth',plotlinewidth,'MarkerSize',24); + + % 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 + %errorbar(tau,sm,sme,'.-b'); set(gca,'XScale','log'); + % this is a hack to approximate the error bars + hold on; plot([tau; tau],[sm+sme; sm-sme],'-k','LineWidth',max(plotlinewidth-1,2)); + + grid on; + title(['Allan Deviation: ' name],'FontSize',FontSize+2,'FontName',FontName); + %set(get(gca,'Title'),'Interpreter','none'); + xlabel('\tau [sec]','FontSize',FontSize,'FontName',FontName); + if isfield(data,'units') + ylabel(['\sigma_y(\tau) [' data.units ']'],'FontSize',FontSize,'FontName',FontName); + else + ylabel('\sigma_y(\tau)','FontSize',FontSize,'FontName',FontName); + end + set(gca,'FontSize',FontSize,'FontName',FontName); + % expand the x axis a little bit so that the errors bars look nice + adax = axis; + axis([adax(1)*0.9 adax(2)*1.1 adax(3) adax(4)]); + + % display the minimum value + fprintf(1,'allan: Minimum ADEV value: %g at tau = %g seconds\n',min(sm),tau(sm==min(sm))); + + elseif verbose >= 1 + 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'); + fprintf(1,'Type "help allan" for more information.\n\n'); + end + +end % end plot ADEV data + +retval = sm; +errorb = sme; + +return diff --git a/allan_modified.m b/allan_modified.m new file mode 100644 index 0000000..df6c83b --- /dev/null +++ b/allan_modified.m @@ -0,0 +1,561 @@ +function [retval, s, errorb, tau] = allan_modified(data,tau,name,verbose) +% ALLAN_MODIFIED Compute the modified Allan deviation for a set of +% time-domain frequency data +% [RETVAL, S, ERRORB, TAU] = ALLAN_MODIFIED(DATA,TAU,NAME,VERBOSE) +% +% Inputs: +% DATA should be a struct and have the following fields: +% DATA.freq or DATA.phase +% A vector of fractional frequency measurements (df/f) in +% DATA.freq *or* phase offset data (seconds) in DATA.phase +% If phase data is not present, it will be generated by +% integrating the fractional frequency data. +% If both fields are present, then DATA.phase will be used. +% +% DATA.rate or DATA.time +% The sampling rate in Hertz (DATA.rate) or a vector of +% timestamps for each measurement in seconds (DATA.time). +% DATA.rate is used if both fields are present. +% If DATA.rate == 0, then the timestamps are used. +% +% TAU is an array of tau values for computing Allan deviation. +% TAU values must be divisible by 1/DATA.rate (data points cannot be +% grouped in fractional quantities!). Invalid values are ignored. +% NAME is an optional label that is added to the plot titles. +% VERBOSE sets the level of status messages: +% 0 = silent & no data plots; 1 = status messages; 2 = all messages +% +% Outputs: +% 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). +% 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. +% 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). +% +% Example: +% +% To compute the modified Allan deviation for the data in the variable "lt": +% >> lt +% lt = +% freq: [1x86400 double] +% rate: 0.5 +% +% Use: +% +% >> adm = allan_modified(lt,[2 10 100],'lt data',1); +% +% 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. +% You can also use the default settings, which are usually a good starting point: +% +% >> adm = allan_modified(lt); +% +% +% Notes: +% This function calculates the modifed Allan deviation (MDEV). +% The calculation is performed using phase data. If only frequency data is +% provided, phase data is generated by integrating the frequency data. +% No pre-processing of the data is performed. +% 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 +% time-stamped data, only tau values greater than the maximum gap between +% samples and less than half the total time are used. +% 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, +% then compare with time-stamp-based. Often the differences are insignificant. +% When phase data is generated by integrating time-stamped frequency data, +% the final data point is dropped, because there is no timestamp for it. +% This will create a [usually small] difference between the results from +% analyzing the same data set with timestamp data and analyzing with a +% fixed sample rate. See note in the code near line 350. +% You can choose between loglog and semilog plotting of results by +% commenting in/out the appropriate line. Search for "#PLOTLOG". +% 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 +% given in IEEE standard 1139-1999, Annex C. +% The author welcomes other validation results, see contact info below. +% +% For more information, see: +% [1] W. J. Riley, "The Calculation of Time Domain Frequency Stability," +% Available on the web: +% http://www.ieee-uffc.org/frequency_control/teaching.asp?name=paper1ht +% +% +% M.A. Hopcroft +% mhopeng at gmail dot com +% +% I welcome your comments and feedback! +% +% MH Mar2014 +% v1.24 fix bug related to generating freq data from phase with timestamps +% (thanks to S. David-Grignot for finding the bug) +% MH Oct2010 +% v1.22 tau truncation to integer groups; tau sort +% plotting bugfix +% v1.20 update to match allan.m (dsplot.m, columns) +% discard tau values with timestamp irregularities +% + +versionstr = 'allan_modified v1.24'; + +% MH MAR2010 +% v1.1 bugfixes for irregular sample rates +% update consistency check +% +% MH FEB2010 +% v1.0 based on allan_overlap v2.0 +% + +%#ok<*AGROW> + + +% defaults +if nargin < 4, verbose = 2; end +if nargin < 3, name=''; end +if nargin < 2 || isempty(tau), tau=2.^(-10:10); end +if isfield(data,'rate') && isempty(data.rate), data.rate=0; end % v1.1 + +% Formatting for plots +FontName = 'Arial'; +FontSize = 14; +plotlinewidth=2; + +if verbose >= 1, fprintf(1,'allan_modified: %s\n\n',versionstr); end + +%% Data consistency checks +if ~(isfield(data,'phase') || isfield(data,'freq')) + error('Either ''phase'' or ''freq'' must be present in DATA. See help file for details. [con0]'); +end +if isfield(data,'time') + if isfield(data,'phase') && (length(data.phase) ~= length(data.time)) + 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]'); + else + error('The time and phase vectors are not the same length. See help for details. [con1]'); + end + end + if isfield(data,'phase') && (any(isnan(data.phase)) || any(isinf(data.phase))) + error('The phase vector contains invalid elements (NaN/Inf). [con3]'); + end + if isfield(data,'freq') && (any(isnan(data.freq)) || any(isinf(data.freq))) + error('The freq vector contains invalid elements (NaN/Inf). [con4]'); + end + if isfield(data,'time') && (any(isnan(data.time)) || any(isinf(data.time))) + error('The time vector contains invalid elements (NaN/Inf). [con5]'); + end +end + +% sort tau vector +tau=sort(tau); + + +%% Basic statistical tests on the data set +if ~isfield(data,'freq') + if isfield(data,'rate') && data.rate ~= 0 + data.freq=diff(data.phase).*data.rate; + elseif isfield(data,'time') + data.freq=diff(data.phase)./diff(data.time); + end + if verbose >= 1, fprintf(1,'allan_modified: Fractional frequency data generated from phase data (M=%g).\n',length(data.freq)); end +end +if size(data.freq,2) > size(data.freq,1), data.freq=data.freq'; end % ensure columns + +s.numpoints=length(data.freq); +s.max=max(data.freq); +s.min=min(data.freq); +s.mean=mean(data.freq); +s.median=median(data.freq); +if isfield(data,'time') + 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); +elseif isfield(data,'rate') && data.rate ~= 0; + s.linear=polyfit((1/data.rate:1/data.rate:length(data.freq)/data.rate)',data.freq,1); +else + error('Either "time" or "rate" must be present in DATA. Type "help allan_modified" for details. [err1]'); +end +s.std=std(data.freq); + +if verbose >= 2 + fprintf(1,'allan_modified: fractional frequency data statistics:\n'); + disp(s); +end + +% scale to median for plotting +medianfreq=data.freq-s.median; +sm=[]; sme=[]; + +% Screen for outliers using 5x Median Absolute Deviation (MAD) criteria +MAD = median(abs(medianfreq)/0.6745); +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'); +end + +%%%% +% There are two cases, either using timestamps or rate: + +%% Fixed Sample Rate Data +% If there is a regular interval between measurements, calculation is much +% easier/faster +if isfield(data,'rate') && data.rate > 0 % if data rate was given + if verbose >= 1 + fprintf(1, 'allan_modified: regular data '); + if isfield(data,'freq') + fprintf(1, '(%g freq data points @ %g Hz)\n',length(data.freq),data.rate); + elseif isfield(data,'phase') + fprintf(1, '(%g phase data points @ %g Hz)\n',length(data.phase),data.rate); + else + error('\n phase or freq data missing [err10]'); + end + end + + % string for plot title + name=[name ' (' num2str(data.rate) ' Hz)']; + + % what is the time interval between data points? + tmstep = 1/data.rate; + + % Is there time data? Just for curiosity/plotting, does not impact calculation + if isfield(data,'time') + % adjust time data to remove any starting gap; first time step + % should not be zero for comparison with freq data + dtime=data.time-data.time(1)+mean(diff(data.time)); + dtime=dtime(1:length(medianfreq)); % equalize the data vector lengths for plotting (v1.1) + if verbose >= 2 + fprintf(1,'allan_modified: End of timestamp data: %g sec.\n',dtime(end)); + 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))); + end + end + else + % create time axis data using rate (for plotting only) + dtime=(tmstep:tmstep:length(data.freq)*tmstep); + end + + + % is phase data present? If not, generate it + if ~isfield(data,'phase') + nfreq=data.freq-s.mean; + dphase=zeros(1,length(nfreq)+1); + 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 + else + dphase=data.phase; + end + + + % check the range of tau values and truncate if necessary + % find halfway point of time record + halftime = round(tmstep*length(data.freq)/2); + % truncate tau to appropriate values + 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 + + % find the number of data points in each tau group + % number of samples + N=length(dphase); + m = data.rate.*tau; + % only integer values allowed (no fractional groups of points) + %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 + %m = m(m-round(m)<1e-8); % change to round(m) for integer test v1.22 + m = m(m==round(m)); + %m=round(m); + + if verbose >= 1, fprintf(1,'allan_modified: calculating modified Allan deviation...\n '); end + + + % calculate the modified Allan deviation for each value of tau + k=0; tic; + for i = tau + k=k+1; + pa=[]; + if verbose >= 2, fprintf(1,'%d ',i); end + + mphase = dphase; + + % calculate overlapping "phase averages" (x_k) + for p=1:m(k) + + % truncate frequency set length to an even multiple of this tau value + mphase=mphase(1:end-rem(length(mphase),m(k))); + % group phase values + mp=reshape(mphase,m(k),[]); + % find average in each "tau group" (each column of mp) + pa(p,:)=mean(mp,1); + % shift data vector by -1 and repeat + mphase=circshift(dphase,(size(dphase)>1)*-p); + + end + + % create "modified" y_k freq values + mfreq=diff(pa,1,2)./i; + mfreq=reshape(mfreq,1,[]); + + % calculate modified frequency differences + mfreqd=reshape(mfreq,m(k),[]); % Vectorize! + mfreqd=diff(mfreqd,1,2); + mfreqd=reshape(mfreqd,1,[]); + + + % 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))); + + % estimate error bars + sme(k)=sm(k)/sqrt(N-3*m(k)+1); + + + end % repeat for each value of tau + + 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 + + + +%% Irregular data (timestamp) +elseif isfield(data,'time') + % the interval between measurements is irregular + % so we must group the data by time + if verbose >= 1, fprintf(1, 'allan_modified: irregular rate data (no fixed sample rate)\n'); end + + % string for plot title + name=[name ' (timestamp)']; + + % adjust time to remove any initial offset + dtime=data.time-data.time(1)+mean(diff(data.time)); + %dtime=data.time-data.time(1); + % where is the maximum gap in time record? + gap_pos=find(diff(dtime)==max(diff(dtime))); + % what is average data spacing? + avg_gap = mean(diff(dtime)); + + if verbose >= 2 + 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'); + 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'); + pause; + end + + if verbose >= 1 + 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); + 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)); + end + if max(diff(dtime)) > 5*avg_gap + fprintf(1, ' WARNING: Max. gap in time record is suspiciously large (>5x the average interval).\n'); + end + end + + % is phase data present? If not, generate it + if ~isfield(data,'phase') + nfreq=data.freq-s.mean; + % NOTE: uncommenting the following two lines will artificially + % allow the code to give equivalent results for validation data + % sets using fixed rate data and timestamped data by adding a + % "phantom" data point for frequency integration. Use of this + % phantom point can skew the results of calculations on real data. + %nfreq(end+1)=0; % phantom freq point, with average value + %dtime(end+1)=dtime(end)+avg_gap; % phantom average time step + dphase=zeros(1,length(nfreq)); + 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 + else + dphase=data.phase; + end + + % find halfway point + halftime = fix(dtime(end)/2); + % truncate tau to appropriate values + tau = tau(tau >= max(diff(dtime)) & tau <= halftime); + if isempty(tau) + error('allan_modified: ERROR: no appropriate tau values (> %g s, < %g s)\n',max(diff(dtime)),halftime); + end + +% % save the freq data for the loop +% dfreq=data.freq; + + % number of samples + N=length(dphase); + m=round(tau./mean(diff(dtime))); + + if verbose >= 1, fprintf(1,'allan_modified: calculating modified Allan deviation...\n'); end + + k=0; tic; + for i = tau + + k=k+1; pa=[]; + + mphase = dphase; time = dtime; + + if verbose >= 2, fprintf(1,'%d ',i); end + + % calculate overlapping "phase averages" (x_k) + %for j = 1:i + for j = 1:m(k) % (v1.1) + km=0; + %fprintf(1,'j: %d ',j); + + % (v1.1) truncating not correct for overlapping samples + % truncate data set to an even multiple of this tau value + %mphase = mphase(time <= time(end)-rem(time(end),i)); + %time = time(time <= time(end)-rem(time(end),i)); + + + % break up the data into overlapping groups of tau length + while i*km < time(end) + km=km+1; + + % progress bar + if verbose >= 2 + if rem(km,100)==0, fprintf(1,'.'); end + if rem(km,1000)==0, fprintf(1,'%g/%g\n',km,round(time(end)/i)); end + end + + mp = mphase(i*(km-1) < (time) & (time) <= i*km); + + if ~isempty(mp) + pa(j,km)=mean(mp); + else + pa(j,km)=0; + end + + end + + % shift data vector by -1 and repeat + mphase=circshift(dphase,(size(mphase)>1)*-j); + mphase(end-j+1:end)=[]; + time=circshift(dtime,(size(time)>1)*-j); + time(end-j+1:end)=[]; + time=time-time(1)+avg_gap; % remove time offset + + end + + % create "modified" y_k freq values + mfreq=diff(pa,1,2)./i; + mfreq=reshape(mfreq,1,[]); + + % calculate modified frequency differences + mfreqd=reshape(mfreq,m(k),[]); % Vectorize! + mfreqd=diff(mfreqd,1,2); + mfreqd=reshape(mfreqd,1,[]); + + % calculate two-sample variance for this tau + % only the first N-3*m(k)+1 samples are valid + 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))); + + % estimate error bars + sme(k)=sm(k)/sqrt(N); + + if verbose >= 2, fprintf(1,'\n'); end + else + if verbose >=2, fprintf(1,' tau=%g dropped due to timestamp irregularities\n',tau(k)); end + sm(k)=0; sme(k)=0; + end + + + end + + 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 + + % remove any points that were dropped + tau(sm==0)=[]; + sm(sm==0)=[]; + sme(sme==0)=[]; + + % modify time vector for plotting + dtime=dtime(1:length(medianfreq)); + +else + error('allan_modified: WARNING: no DATA.rate or DATA.time! Type "help allan_modified" for more information. [err2]'); +end + + +%%%%%%%% +%% Plotting + +if verbose >= 2 % show all data + + % 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 + try + % dsplot makes a new figure + hd=dsplot(dtime,medianfreq); + catch ME + figure; + 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 + if verbose >= 2, fprintf(1,' (Message: %s)\n',ME.message); end + end + set(hd,'Marker','.','LineStyle','none','Color','b'); % equivalent to '.-' + hold on; + + fx = xlim; + % plot([fx(1) fx(2)],[s.median s.median],'-k'); + plot([fx(1) fx(2)],[0 0],':k'); + % show 5x Median Absolute deviation (MAD) values + hm=plot([fx(1) fx(2)],[5*MAD 5*MAD],'-r'); + plot([fx(1) fx(2)],[-5*MAD -5*MAD],'-r'); + % show linear fit line + hf=plot(xlim,polyval(s.linear,xlim)-s.median,'-g'); + title(['Data: ' name],'FontSize',FontSize+2,'FontName',FontName); + %set(get(gca,'Title'),'Interpreter','none'); + xlabel('Time [sec]','FontSize',FontSize,'FontName',FontName); + if isfield(data,'units') + ylabel(['data - median(data) [' data.units ']'],'FontSize',FontSize,'FontName',FontName); + else + ylabel('freq - median(freq)','FontSize',FontSize,'FontName',FontName); + end + 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)); + % tighten up + xlim([dtime(1) dtime(end)]); + + +end % end plot raw data + + +if verbose >= 1 % show analysis results + + % plot Allan deviation results + if ~isempty(sm) + figure + + % Choose loglog or semilogx plot here #PLOTLOG + %semilogx(tau,sm,'.-b','LineWidth',plotlinewidth,'MarkerSize',24); + loglog(tau,sm,'.-b','LineWidth',plotlinewidth,'MarkerSize',24); + + % 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 + %errorbar(tau,sm,sme,'.-b'); set(gca,'XScale','log'); + % this is a hack to approximate the error bars + hold on; plot([tau; tau],[sm+sme; sm-sme],'-k','LineWidth',max(plotlinewidth-1,2)); + + grid on; + title(['Modified Allan Deviation: ' name],'FontSize',FontSize+2,'FontName',FontName); + %set(get(gca,'Title'),'Interpreter','none'); + xlabel('\tau [sec]','FontSize',FontSize,'FontName',FontName); + ylabel('Modified \sigma_y(\tau)','FontSize',FontSize,'FontName',FontName); + set(gca,'FontSize',FontSize,'FontName',FontName); + % expand the x axis a little bit so that the errors bars look nice + adax = axis; + axis([adax(1)*0.9 adax(2)*1.1 adax(3) adax(4)]); + + % display the minimum value + fprintf(1,'allan: Minimum modified ADEV value: %g at tau = %g seconds\n',min(sm),tau(sm==min(sm))); + + elseif verbose >= 1 + 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'); + fprintf(1,'Type "help allan_modified" for more information.\n\n'); + end + +end % end plot analysis + +retval = sm; +errorb = sme; + +return diff --git a/allan_overlap.m b/allan_overlap.m new file mode 100644 index 0000000..55689b8 --- /dev/null +++ b/allan_overlap.m @@ -0,0 +1,547 @@ +function [retval, s, errorb, tau] = allan_overlap(data,tau,name,verbose) +% ALLAN_OVERLAP Compute the overlapping Allan deviation for a set of +% time-domain frequency data +% [RETVAL, S, ERRORB, TAU] = ALLAN_OVERLAP(DATA,TAU,NAME,VERBOSE) +% +% Inputs: +% DATA should be a struct and have the following fields: +% DATA.freq or DATA.phase +% A vector of fractional frequency measurements (df/f) in +% DATA.freq *or* phase offset data (seconds) in DATA.phase +% If phase data is not present, it will be generated by +% integrating the fractional frequency data. +% If both fields are present, then DATA.phase will be used. +% +% DATA.rate or DATA.time +% The sampling rate in Hertz (DATA.rate) or a vector of +% timestamps for each measurement in seconds (DATA.time). +% DATA.rate is used if both fields are present. +% If DATA.rate == 0, then the timestamps are used. +% +% TAU is an array of tau values for computing Allan deviation. +% TAU values must be divisible by 1/DATA.rate (data points cannot be +% grouped in fractional quantities!). Invalid values are ignored. +% NAME is an optional label that is added to the plot titles. +% VERBOSE sets the level of status messages: +% 0 = silent & no data plots; 1 = status messages; 2 = all messages +% +% Outputs: +% 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). +% 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. +% 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). +% +% Example: +% +% To compute the overlapping Allan deviation for the data in the variable "lt": +% >> lt +% lt = +% freq: [1x86400 double] +% rate: 0.5 +% +% Use: +% +% >> ado = allan_overlap(lt,[2 10 100],'lt data',1); +% +% The Allan deviation will be computed and plotted at tau = 2,10,100 seconds. +% 1-sigma confidence intervals will be indicated by vertical lines. +% You can also use the default settings, which are usually a good starting point: +% +% >> ado = allan_overlap(lt); +% +% +% Notes: +% This function calculates the overlapping Allan deviation (ADEV), *not* the +% standard ADEV. Use "allan.m" for standard ADEV. +% The calculation is performed using phase data. If only frequency data is +% provided, phase data is generated by integrating the frequency data. +% However, the timestamp-based calculation is performed using frequency +% data. Phase data is differentiated to generate frequency data if necessary. +% No pre-processing of the data is performed, except to remove any +% initial offset in the time record. +% 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 +% time-stamped data, only tau values greater than the maximum gap between +% samples and less than half the total time are used. +% 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, +% then compare with time-stamp-based. Often the differences are insignificant. +% 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 +% overlapping ADEV; a more accurate (and usually smaller uncertainty) +% value can be determined from chi-squared statistics, but that is not +% implemented in this version. +% You can choose between loglog and semilog plotting of results by +% commenting in/out the appropriate line. Search for "#PLOTLOG". +% 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 +% given in IEEE standard 1139-1999, Annex C. +% The author welcomes other validation results, see contact info below. +% +% For more information, see: +% [1] W. J. Riley, "Addendum to a test suite for the calculation of time domain +% frequency stability," presented at IEEE Frequency Control Symposium, +% 1996. +% Available on the web: +% http://www.ieee-uffc.org/frequency_control/teaching.asp?name=paper1ht +% +% +% M.A. Hopcroft +% mhopeng at gmail dot com +% +% I welcome your comments and feedback! +% +% MH Mar2014 +% v2.24 fix bug related to generating freq data from phase with timestamps +% (thanks to S. David-Grignot for finding the bug) +% MH Oct2010 +% v2.22 tau truncation to integer groups; tau sort +% plotting bugfix +% v2.20 update to match allan.m (dsplot.m, columns) +% discard tau values with timestamp irregularities + +versionstr = 'allan_overlap v2.24'; + +% +% MH MAR2010 +% v2.1 bugfixes for irregular sample rates +% (thanks to Ryad Ben-El-Kezadri for feedback and testing) +% handle empty rate field +% fix integer comparisons for fractional sample rates +% update consistency check +% +% MH FEB2010 +% v2.0 use phase data for calculation- much faster +% Consistent code behaviour for all "allan_x.m" functions: +% accept phase data +% verbose levels +% +% MH JAN2010 +% v1.0 based on allan v1.84 +% + +%#ok<*AGROW> + + +% defaults +if nargin < 4, verbose = 2; end +if nargin < 3, name=''; end +if nargin < 2 || isempty(tau), tau=2.^(-10:10); end +if isfield(data,'rate') && isempty(data.rate), data.rate=0; end % v2.1 + +% Formatting for plots +FontName = 'Arial'; +FontSize = 14; +plotlinewidth=2; + +if verbose >= 1, fprintf(1,'allan_overlap: %s\n\n',versionstr); end + +%% Data consistency checks v2.1 +if ~(isfield(data,'phase') || isfield(data,'freq')) + error('Either ''phase'' or ''freq'' must be present in DATA. See help file for details. [con0]'); +end +if isfield(data,'time') + if isfield(data,'phase') && (length(data.phase) ~= length(data.time)) + 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]'); + else + error('The time and phase vectors are not the same length. See help for details. [con1]'); + end + end + if isfield(data,'phase') && (any(isnan(data.phase)) || any(isinf(data.phase))) + error('The phase vector contains invalid elements (NaN/Inf). [con3]'); + end + if isfield(data,'freq') && (any(isnan(data.freq)) || any(isinf(data.freq))) + error('The freq vector contains invalid elements (NaN/Inf). [con4]'); + end + if isfield(data,'time') && (any(isnan(data.time)) || any(isinf(data.time))) + error('The time vector contains invalid elements (NaN/Inf). [con5]'); + end +end + +% sort tau vector +tau=sort(tau); + +%% Basic statistical tests on the data set +if ~isfield(data,'freq') + if isfield(data,'rate') && data.rate ~= 0 + data.freq=diff(data.phase).*data.rate; + elseif isfield(data,'time') + data.freq=diff(data.phase)./diff(data.time); + end + if verbose >= 1, fprintf(1,'allan_overlap: Fractional frequency data generated from phase data (M=%g).\n',length(data.freq)); end +end +if size(data.freq,2) > size(data.freq,1), data.freq=data.freq'; end % ensure columns + +s.numpoints=length(data.freq); +s.max=max(data.freq); +s.min=min(data.freq); +s.mean=mean(data.freq); +s.median=median(data.freq); +if isfield(data,'time') + 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); +elseif isfield(data,'rate') && data.rate ~= 0; + s.linear=polyfit((1/data.rate:1/data.rate:length(data.freq)/data.rate)',data.freq,1); +else + error('Either "time" or "rate" must be present in DATA. Type "help allan_overlap" for details. [err1]'); +end +s.std=std(data.freq); + +if verbose >= 2 + fprintf(1,'allan_overlap: fractional frequency data statistics:\n'); + disp(s); +end + + +% scale to median for plotting +medianfreq=data.freq-s.median; +sm=[]; sme=[]; + +% Screen for outliers using 5x Median Absolute Deviation (MAD) criteria +MAD = median(abs(medianfreq)/0.6745); +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'); +end + +%%%% +% There are four cases, freq or phase data, using timestamps or rate: + +%% Fixed Sample Rate Data +% If there is a regular interval between measurements, calculation is much +% easier/faster +if isfield(data,'rate') && data.rate > 0 % if data rate was given + if verbose >= 1 + fprintf(1, 'allan_overlap: regular data '); + if isfield(data,'freq') + fprintf(1, '(%g freq data points @ %g Hz)\n',length(data.freq),data.rate); + elseif isfield(data,'phase') + fprintf(1, '(%g phase data points @ %g Hz)\n',length(data.phase),data.rate); + else + error('\n phase or freq data missing [err10]'); + end + end + + % string for plot title + name=[name ' (' num2str(data.rate) ' Hz)']; + + % what is the time interval between data points? + tmstep = 1/data.rate; + + % Is there time data? Just for curiosity/plotting, does not impact calculation + if isfield(data,'time') + % adjust time data to remove any starting gap; first time step + % should not be zero for comparison with freq data + dtime=data.time-data.time(1)+mean(diff(data.time)); + dtime=dtime(1:length(medianfreq)); % equalize the data vector lengths for plotting (v2.1) + if verbose >= 2 + fprintf(1,'allan_overlap: End of timestamp data: %g sec.\n',dtime(end)); + 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))); + end + end + else + % create time axis data using rate (for plotting only) + dtime=(tmstep:tmstep:length(data.freq)*tmstep); + end + + + % is phase data present? If not, generate it + if ~isfield(data,'phase') + nfreq=data.freq-s.mean; + dphase=zeros(1,length(nfreq)+1); + 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 + else + dphase=data.phase; + end + + % check the range of tau values and truncate if necessary + % find halfway point of time record + halftime = round(tmstep*length(data.freq)/2); + % truncate tau to appropriate values + 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 + + % number of samples + N=length(dphase); + % number of samples per tau period + m = data.rate.*tau; + % only integer values allowed for m (no fractional groups of points) + %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 + %m = m(m-round(m)<1e-8); % change to round(m) for integer test v2.22 + m = m(m==round(m)); + %m=round(m); + %fprintf(1,'m: %.50f\n',m) + + if verbose >= 1, fprintf(1,'allan_overlap: calculating overlapping Allan deviation...\n '); end + + % calculate the Allan deviation for each value of tau + k=0; tic; + for i = tau + k=k+1; + if verbose >= 2, fprintf(1,'%d ',i); end + + + % pad phase data set length to an even multiple of this tau value + mphase=zeros(ceil(length(dphase)./m(k))*m(k),1); + mphase(1:N)=dphase; + % group phase values + mp=reshape(mphase,m(k),[]); + % compute second differences of phase values (x_k+m - x_k) + md1=diff(mp,1,2); + md2=diff(md1,1,2); + md1=reshape(md2,1,[]); + + % compute overlapping ADEV from phase values + % 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)); + + % estimate error bars + sme(k)=sm(k)/sqrt(N-2*m(k)); + + + end % repeat for each value of tau + + 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 + + + +%% Irregular data, no fixed interval +elseif isfield(data,'time') + % the interval between measurements is irregular + % so we must group the data by time + if verbose >= 1, fprintf(1, 'allan_overlap: irregular rate data (no fixed sample rate)\n'); end + + + % string for plot title + name=[name ' (timestamp)']; + + + % adjust time to remove any starting offset + dtime=data.time-data.time(1)+mean(diff(data.time)); + + % save the freq data for the loop + dfreq=data.freq; + dtime=dtime(1:length(dfreq)); + + dfdtime=diff(dtime); % only need to do this once (v2.1) + % where is the maximum gap in time record? + gap_pos=find(dfdtime==max(dfdtime)); + % what is average data spacing? + avg_gap = mean(dfdtime); + s.avg_rate = 1/avg_gap; % save avg rate for user (v2.1) + + if verbose >= 2 + 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'); + 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'); + pause; + end + + if verbose >= 1 + 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); + 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)); + end + if max(diff(dtime)) > 5*avg_gap + fprintf(1, ' WARNING: Max. gap in time record is suspiciously large (>5x the average interval).\n'); + end + end + + + % find halfway point + halftime = fix(dtime(end)/2); + % truncate tau to appropriate values + tau = tau(tau >= max(dfdtime) & tau <= halftime); + if isempty(tau) + error('allan_overlap: ERROR: no appropriate tau values (> %g s, < %g s)\n',max(dfdtime),halftime); + end + + + % number of samples + M=length(dfreq); + % number of samples per tau period + m=round(tau./avg_gap); + + if verbose >= 1, fprintf(1,'allan_overlap: calculating overlapping Allan deviation...\n'); end + + k=0; tic; + for i = tau + k=k+1; + fa=[]; + + if verbose >= 2, fprintf(1,'%d ',i); end + + freq = dfreq; time = dtime; + + + % compute overlapping samples (y_k) for this tau + %for j = 1:i + for j = 1:m(k) % (v2.1) + km=0; + %fprintf(1,'j: %d ',j); + + % (v2.1) truncating not correct for overlapping samples + % truncate data set to an even multiple of this tau value + %freq = freq(time <= time(end)-rem(time(end),i)); + %time = time(time <= time(end)-rem(time(end),i)); + + % break up the data into overlapping groups of tau length + while i*km <= time(end) + km=km+1; + %i*km + + % progress bar + if verbose >= 2 + if rem(km,100)==0, fprintf(1,'.'); end + if rem(km,1000)==0, fprintf(1,'%g/%g\n',km,round(time(end)/i)); end + end + + f = freq(i*(km-1) < (time) & (time) <= i*km); + + if ~isempty(f) + fa(j,km)=mean(f); + else + fa(j,km)=0; + end + + end + %fa + + % shift data vector by -1 and repeat + freq=circshift(dfreq,(size(freq)>1)*-j); + freq(end-j+1:end)=[]; + time=circshift(dtime,(size(time)>1)*-j); + time(end-j+1:end)=[]; + time=time-time(1)+avg_gap; % remove time offset + + end + + % compute second differences of fractional frequency values (y_k+m - y_k) + fd1=diff(fa,1,2); + fd1=reshape(fd1,1,[]); + % compute overlapping ADEV from fractional frequency values + % only the first M-2*m(k)+1 samples are valid + 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)); + + % estimate error bars + sme(k)=sm(k)/sqrt(M+1); + + if verbose >= 2, fprintf(1,'\n'); end + + else + if verbose >=2, fprintf(1,' tau=%g dropped due to timestamp irregularities\n',tau(k)); end + sm(k)=0; sme(k)=0; + end + + + end + + 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 + + % remove any points that were dropped + tau(sm==0)=[]; + sm(sm==0)=[]; + sme(sme==0)=[]; + + + +else + error('allan_overlap: WARNING: no DATA.rate or DATA.time! Type "help allan" for more information. [err2]'); +end + + +%%%%%%%% +%% Plotting + +if verbose >= 2 % show all data + + % 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 + try + % dsplot makes a new figure + hd=dsplot(dtime,medianfreq); + catch ME + figure; + 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 + if verbose >= 2, fprintf(1,' (Message: %s)\n',ME.message); end + end + set(hd,'Marker','.','LineStyle','none','Color','b'); % equivalent to '.-' + hold on; + + fx = xlim; + % plot([fx(1) fx(2)],[s.median s.median],'-k'); + plot([fx(1) fx(2)],[0 0],':k'); + % show 5x Median Absolute deviation (MAD) values + hm=plot([fx(1) fx(2)],[5*MAD 5*MAD],'-r'); + plot([fx(1) fx(2)],[-5*MAD -5*MAD],'-r'); + % show linear fit line + hf=plot(xlim,polyval(s.linear,xlim)-s.median,'-g'); + title(['Data: ' name],'FontSize',FontSize+2,'FontName','Arial'); + %set(get(gca,'Title'),'Interpreter','none'); + xlabel('Time [sec]','FontSize',FontSize,'FontName',FontName); + if isfield(data,'units') + ylabel(['data - median(data) [' data.units ']'],'FontSize',FontSize,'FontName',FontName); + else + ylabel('freq - median(freq)','FontSize',FontSize,'FontName',FontName); + end + 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)); + % tighten up + xlim([dtime(1) dtime(end)]); + + +end % end plot raw data + + +if verbose >= 1 % show analysis results + + % plot Allan deviation results + if ~isempty(sm) + figure + + % Choose loglog or semilogx plot here #PLOTLOG + %semilogx(tau,sm,'.-b','LineWidth',plotlinewidth,'MarkerSize',24); + loglog(tau,sm,'.-b','LineWidth',plotlinewidth,'MarkerSize',24); + + % 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 + %errorbar(tau,sm,sme,'.-b'); set(gca,'XScale','log'); + % this is a hack to approximate the error bars + hold on; plot([tau; tau],[sm+sme; sm-sme],'-k','LineWidth',max(plotlinewidth-1,2)); + + grid on; + title(['Overlapping Allan Deviation: ' name],'FontSize',FontSize+2,'FontName',FontName); + %set(get(gca,'Title'),'Interpreter','none'); + xlabel('\tau [sec]','FontSize',FontSize,'FontName','Arial'); + ylabel(' Overlapping \sigma_y(\tau)','FontSize',FontSize,'FontName',FontName); + set(gca,'FontSize',FontSize,'FontName',FontName); + % expand the x axis a little bit so that the errors bars look nice + adax = axis; + axis([adax(1)*0.9 adax(2)*1.1 adax(3) adax(4)]); + + % display the minimum value + fprintf(1,'allan: Minimum overlapping ADEV value: %g at tau = %g seconds\n',min(sm),tau(sm==min(sm))); + + elseif verbose >= 1 + 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'); + fprintf(1,'Type "help allan_overlap" for more information.\n\n'); + end + +end % end plot analysis + +retval = sm; +errorb = sme; + +return diff --git a/allanplot.m b/allanplot.m new file mode 100755 index 0000000..85f0471 --- /dev/null +++ b/allanplot.m @@ -0,0 +1,12 @@ +#!/usr/bin/octave-cli --persist + +filename = argv(){1}; +col = eval(argv(){2}); +mult = eval(argv(){3}); + +data.freq = load(filename)(:,col).*mult; +data.rate = 1; + +ad = allan(data, 2.^(0:12), strsplit(filename, '/'){end}, 1); +input("Press to continue..."); +exit diff --git a/dsplot.m b/dsplot.m new file mode 100644 index 0000000..b1d5544 --- /dev/null +++ b/dsplot.m @@ -0,0 +1,319 @@ +function hL = dsplot(x, y, numPoints) + +%DSPLOT Create down sampled plot. +% This function creates a down sampled plot to improve the speed of +% exploration (zoom, pan). +% +% DSPLOT(X, Y) plots Y versus X by downsampling if there are large number +% of elements. X and Y needs to obey the following: +% 1. X must be a monotonically increasing vector. +% 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. +% +% DSPLOT(Y) plots the columns of Y versus their index. +% +% 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. +% +% DSPLOT(X, Y, NUMPOINTS) or DSPLOT(Y, [], NUMPOINTS) specifies the +% number of points (roughly) to display on the screen. The default is +% 50000 points (~390 kB doubles). NUMPOINTS can be a number greater than +% 500. +% +% 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. +% If the signal is stochastic or has a lot of sharp changes, there will +% be more points on plotted on the screen. +% +% The figure title (name) will indicate whether the plot shown is +% downsampled or is the true representation. +% +% 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 +% the base workspace via a menu. +% +% Run the following examples and zoom/pan to see the performance. +% +% Example 1: (with small details) +% x = linspace(0, 2*pi, 1000000); +% y1 = sin(x)+.02*cos(200*x)+0.001*sin(2000*x)+0.0001*cos(20000*x); +% dsplot(x,y1);title('Down Sampled'); +% % compare with +% figure;plot(x,y1);title('Normal Plot'); +% +% Example 2: (with outlier points) +% x = linspace(0, 2*pi, 1000000); +% 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)); +% y1([300000, 700000, 700001, 900000]) = [0, 1, -2, 0.5]; +% y2(300000:500000) = y2(300000:500000) + 1; +% y2(500001:600000) = y2(500001:600000) - 1; +% y2(800000) = 0; +% dsplot(x, [y1;y2]);title('Down Sampled'); +% % compare with +% figure;plot(x, [y1;y2]);title('Normal Plot'); +% +% See also PLOT. + +% Version: +% v1.0 - first version (Aug 1, 2007) +% 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 +% item for saving out the original data back to the base +% workspace. (Aug 10, 2007) +% +% Jiro Doke +% August 1, 2007 + +debugMode = false; + +%-------------------------------------------------------------------------- +% Error checking +error(nargchk(1, 3, nargin, 'struct')); +if nargin < 3 + % Number of points to show on the screen. It's quite possible that more + % points will be displayed if there are outlier points + numPoints = 50000; % ~390 kB for doubles +end +if nargin == 1 || isempty(y) + noXVar = true; + y = x; + x = []; +else + noXVar = false; +end +myErrorCheck; +%-------------------------------------------------------------------------- + +if size(x, 2) > 1 % it's a row vector -> transpose + x = x'; + y = y'; + varTranspose = true; +else + varTranspose = false; +end + +% Number of lines +numSignals = size(y, 2); + +% 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 +% orientation. +if numSignals > size(y, 1) + s = input(sprintf('Are you sure you want to plot %d lines? (y/n) ', ... + numSignals), 's'); + if ~strcmpi(s, 'y') + disp('Canceled. You may want to transpose the matrix.'); + if nargout == 1 + hL = []; + end + return; + end +end + +% Attempt to find outliers. Use a running average technique +filterWidth = ceil(min([50, length(x)/10])); % max window size of 50 +a = y - filter(ones(filterWidth,1)/filterWidth, 1, y); +[iOutliers, jOutliers] = find(abs(a - repmat(mean(a), size(a, 1), 1)) > ... + repmat(4 * std(a), size(a, 1), 1)); +clear a; + +% Always create new figure because it messes around with zoom, pan, +% datacursors. +hFig = figure; +figName = ''; + +% Create template plot using NaNs +hLine = plot(NaN(2, numSignals), NaN(2, numSignals)); +set(hLine, 'tag', 'dsplot_lines'); + +% Define CreateFcn for the figure +set(hFig, 'CreateFcn', @mycreatefcn); +mycreatefcn(); + +% Create menu for exporting data +hMenu = uimenu(hFig, 'Label', 'Data'); +uimenu(hMenu, ... + 'Label' , 'Export data to workspace.', ... + 'Callback', @myExportFcn); + +% Update lines +updateLines([min(x), max(x)]); + +% Deal with output argument +if nargout == 1 + hL = hLine; +end + +%-------------------------------------------------------------------------- + function myExportFcn(varargin) + % 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 + % can get back the data. + + % Determine the variable name + allVarNames = evalin('base', 'who'); + newVarName = genvarname('dsplotData', allVarNames); + + % X + if ~noXVar + if varTranspose + dat.x = x'; + else + dat.x = x; + end + end + + % Y + if varTranspose + dat.y = y'; + else + dat.y = y; + end + + assignin('base', newVarName, dat); + + msgbox(sprintf('Data saved to the base workspace as ''%s''.', ... + newVarName), 'Saved', 'modal'); + + end + +%-------------------------------------------------------------------------- + function mycreatefcn(varargin) + % 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 + % the figure. + + if nargin > 0 + hFig = varargin{1}; + end + hLine = findobj(hFig, 'type', 'axes'); + hLine(strmatch('legend', get(hLine, 'tag'))) = []; + hLine = get(hLine, 'Children'); + + % Create Zoom, Pan, Datacursor objects + hZoom = zoom(hFig); + hPan = pan(hFig); + hDc = datacursormode(hFig); + set(hZoom, 'ActionPostCallback', @mypostcallback); + set(hPan , 'ActionPostCallback', @mypostcallback); + set(hDc , 'UpdateFcn' , @myDCupdatefcn); + + end + +%-------------------------------------------------------------------------- + function mypostcallback(obj, evd) %#ok + % This callback that gets called when the mouse is released after + % zooming or panning. + + % single or double-click + switch get(hFig, 'SelectionType') + case {'normal', 'alt'} + updateLines(xlim(evd.Axes)); + + case 'open' + updateLines([min(x), max(x)]); + + end + + end + +%-------------------------------------------------------------------------- + function updateLines(rng) + % This helper function is for determining the points to plot on the + % screen based on which portion is visible in the current limits. + + % find indeces inside the range + id = find(x >= rng(1) & x <= rng(2)); + + % if there are more points than we want + if length(id) > numPoints / numSignals + + % see how many outlier points are in this range + blah = iOutliers > id(1) & iOutliers < id(end); + + % determine indeces of points to plot. + idid = round(linspace(id(1), id(end), round(numPoints/numSignals)))'; + + x2 = cell(numSignals, 1); + y2 = x2; + for iSignals = 1:numSignals + % add outlier points + ididid = unique([idid; iOutliers(blah & jOutliers == iSignals)]); + x2{iSignals} = x(ididid); + y2{iSignals} = y(ididid, iSignals); + end + + if debugMode + figName = ['downsampled - ', sprintf('%d, ', cellfun('length', y2))]; + else + figName = 'downsampled'; + end + + else % no need to down sample + figName = 'true'; + + x2 = repmat({x(id)}, numSignals, 1); + y2 = mat2cell(y(id, :), length(id), ones(1, numSignals))'; + + end + + % Update plot + set(hLine, {'xdata', 'ydata'} , [x2, y2]); + set(hFig, 'Name', figName); + + end + +%-------------------------------------------------------------------------- + function txt = myDCupdatefcn(empt, event_obj) %#ok + % This function displays appropriate data cursor message based on the + % display type + + pos = get(event_obj,'Position'); + switch figName + case 'true' + txt = {['X: ',num2str(pos(1))],... + ['Y: ',num2str(pos(2))]}; + otherwise + txt = {['X: ',num2str(pos(1))],... + ['Y: ',num2str(pos(2))], ... + 'Warning: Downsampled', ... + 'May not be accurate'}; + end + end + +%-------------------------------------------------------------------------- + function myErrorCheck + % Do some error checking on the input arguments. + + if ~isa(numPoints, 'double') || numel(numPoints) > 1 || numPoints < 500 + error('Third argument must be a scalar greater than 500'); + end + if ~isnumeric(x) || ~isnumeric(y) + error('Arguments must be numeric'); + end + if length(size(x)) > 2 || length(size(y)) > 2 + error('Only 2-D data accepted'); + end + + % If only one input, create index vector X + if isempty(x) + if ismember(1, size(y)) + x = reshape(1:numel(y), size(y)); + else + x = (1:size(y, 1))'; + end + end + + if ~ismember(1, size(x)) + error('First argument has to be a vector'); + end + 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'); + end + if any(diff(x) <= 0) + error('The first argument has to be a monotonically increasing vector'); + end + end + +end \ No newline at end of file diff --git a/pt100.m b/pt100.m new file mode 100644 index 0000000..aa895de --- /dev/null +++ b/pt100.m @@ -0,0 +1,5 @@ +function T = pt100(R) + +pt100 = @(R) (R/100.-1)/0.003850+273.15; + +T = pt100(R); diff --git a/res2temp16941.m b/res2temp16941.m new file mode 100644 index 0000000..d3fe985 --- /dev/null +++ b/res2temp16941.m @@ -0,0 +1,5 @@ +function T = res2temp16941(R) + +res2temp16941 = @(R) 10.^(2.9486 * (log10(1000./R)).^2 + 4.5862 * log10(1000./R) + 2.266); + +T = res2temp16941(R); diff --git a/res2temp16943.m b/res2temp16943.m new file mode 100644 index 0000000..cfc71b1 --- /dev/null +++ b/res2temp16943.m @@ -0,0 +1,5 @@ +function T = res2temp16943(R) + +res2temp16943 = @(R) 10.^(3.4738 * (log10(1000./R)).^2 + 5.1198 * log10(1000./R) + 2.3681); + +T = res2temp16943(R); diff --git a/res2temp16944.m b/res2temp16944.m new file mode 100644 index 0000000..ca82056 --- /dev/null +++ b/res2temp16944.m @@ -0,0 +1,5 @@ +function T = res2temp16944(R) + +res2temp16944 = @(R) 10.^(3.3674 * (log10(1000./R)).^2 + 5.2874 * log10(1000./R) + 2.5165); + +T = res2temp16944(R); diff --git a/res2temp16945.m b/res2temp16945.m new file mode 100644 index 0000000..98461b1 --- /dev/null +++ b/res2temp16945.m @@ -0,0 +1,5 @@ +function T = res2temp16945(R) + +res2temp16945 = @(R) 10.^(3.2497 * (log10(1000./R)).^2 + 5.1777 * log10(1000./R) + 2.499); + +T = res2temp16945(R); diff --git a/res2temp16947.m b/res2temp16947.m new file mode 100644 index 0000000..880380f --- /dev/null +++ b/res2temp16947.m @@ -0,0 +1,5 @@ +function T = res2temp16947(R) + +res2temp16947 = @(R) 10.^(3.4597 * (log10(1000./R)).^2 + 5.2422 * log10(1000./R) + 2.4169); + +T = res2temp16947(R); diff --git a/res2temp625.m b/res2temp625.m new file mode 100644 index 0000000..9be4607 --- /dev/null +++ b/res2temp625.m @@ -0,0 +1,5 @@ +function T = res2temp625(R) + +res2temp625 = @(R) 0.333548856582638109 + 11.7361551595386118 * (1000./R) + -31.32988932320903987 * (1000./R).^2 + 262.878643524833024 * (1000./R).^3 + -704.163538021035492 * (1000./R).^4 + 1056.6040485650301 * (1000./R).^5 + -307.057196729816496 * (1000./R).^6; + +T = res2temp625(R); diff --git a/res2temp627.m b/res2temp627.m new file mode 100644 index 0000000..ad02c46 --- /dev/null +++ b/res2temp627.m @@ -0,0 +1,5 @@ +function T = res2temp627(R) + +res2temp627 = @(R) 0.399341181655472610 + 10.8420092277810909 * (1000./R) + -26.4597939187660813 * (1000./R).^2 + 245.9828566655493379 * (1000./R).^3 + -668.069876596331596 * (1000./R).^4 + 1001.69882618263364 * (1000./R) .^5 + -267.272089680656791 * (1000./R).^6; + +T = res2temp627(R); diff --git a/res2temp628.m b/res2temp628.m new file mode 100644 index 0000000..49d38cb --- /dev/null +++ b/res2temp628.m @@ -0,0 +1,5 @@ +function T = res2temp628(R) + +res2temp628 = @(R) 0.463200932294057566 + 13.5049710820894688 * (1000./R) + -30.5191222755238414 * (1000./R).^2 + 231.098593852017075* (1000./R).^3 + -550.122691885568202 * (1000./R).^4 + 806.038547554984689 * (1000./R).^5 + -198.510489917360246 * (1000./R).^6; + +T = res2temp628(R);