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