From b197c3fdf011f9972e68675ac4ca59e0cd09ad67 Mon Sep 17 00:00:00 2001 From: bma Date: Tue, 26 Sep 2017 20:16:36 +0200 Subject: [PATCH] first commit --- allan.m | 576 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ allan_modified.m | 561 +++++++++++++++++++++++++++++++++++++++++++++++++++++ allan_overlap.m | 547 ++++++++++++++++++++++++++++++++++++++++++++++++++++ allanplot.m | 12 ++ dsplot.m | 319 ++++++++++++++++++++++++++++++ pt100.m | 5 + res2temp16941.m | 5 + res2temp16943.m | 5 + res2temp16944.m | 5 + res2temp16945.m | 5 + res2temp16947.m | 5 + res2temp625.m | 5 + res2temp627.m | 5 + res2temp628.m | 5 + 14 files changed, 2060 insertions(+) create mode 100644 allan.m create mode 100644 allan_modified.m create mode 100644 allan_overlap.m create mode 100755 allanplot.m create mode 100644 dsplot.m create mode 100644 pt100.m create mode 100644 res2temp16941.m create mode 100644 res2temp16943.m create mode 100644 res2temp16944.m create mode 100644 res2temp16945.m create mode 100644 res2temp16947.m create mode 100644 res2temp625.m create mode 100644 res2temp627.m create mode 100644 res2temp628.m 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); -- 2.16.4