Hi William,


On Aug 1, 2012, at 5:34 AM, William Katsak wrote:

> Thanks for the insight. I would definitely be interested in your
> program for measuring the overhead.

	Ah great. So here is what to do:
0) IUse "traceroute 8.8.8.8" to figure out the IP address of the first hop on the other side of your modem (but see below)

1) from my MacBooks terminal I edit and run ping_sweeper_02.sh (inlined below) which talks around 6 hours (I just run this overnight)
	assign the first outside hop's IP address to TARGET

#! /bin/bash
# ideally one would shuffle all size*repetition to better average out transient issues
#
# USEAGE
# 1) use "traceroute 8.8.8.8" and find the network hops just outside of your modem
# 	typically for ATM the modem to DSLAM link will introduce an additional >10ms latency 
#	so the hop directly after the first big latency jump should be close to the DSLAM, use 
#	this hops IP address should be used as TARGET below
# 2) optional confirm this IPs reliability by running "ping -c 100" against this IP
#	the min and max of the ping times should stay close to the average and the stddev
#	 should be small as well, otherwise PINGSPERSIZE might need to be set to >=100
#

TECH=ATM			# CABLE, ATM, ... just to name the log file
DATESTR=`date +%Y%m%d_%H%M%S`	# to allow multiple sequential records
TARGET=96.34.97.78	# use traceroute something far to get the first hop out of the home net
LOG=ping_sweep_${TECH}_${DATESTR}.txt


PINGSPERSIZE=50		# 480 packet sizes take 13.33 hours at 100 repetitions, so stick to 50
SWEEPMINSIZE=16		# on macosx 10.7 64bit, sizes < 16 do not have a timestamp...
SWEEPMAXSIZE=496	# initial size plus 10 full ATM cells (of 48 bytes payload each)
n_SWEEPS=`expr ${SWEEPMAXSIZE} - ${SWEEPMINSIZE}`

# do it
echo "Sweeping ${n_SWEEPS} * ${PINGSPERSIZE} pings to ${TARGET} might take a while..."
echo "ping -c ${PINGSPERSIZE} -g ${SWEEPMINSIZE} -G ${SWEEPMAXSIZE} ${TARGET} > ${LOG}"
ping -c ${PINGSPERSIZE} -g ${SWEEPMINSIZE} -G ${SWEEPMAXSIZE} ${TARGET} > ${LOG} &

# show some progress
tail -f ${LOG}

echo "Done... ($0)"

3) run the following in matlab or octave (tested with octave 3.6) (inlined and attached) as tc_stab_parameter_guide_02('filly_qualified_path_to_the_ping_log', your_uplink_rate_in_Kbit, your_downlink_rate_in_Kbit) and then wait a bit (it might take a minute or two)

function [ output_args ] = tc_stab_parameter_guide_02( sweep_fqn, up_Kbit, down_Kbit )
%TC_STAB_PARAMETER_GUIDE Summary of this function goes here
%   try to read in the result from a ping sweep run
%	sweep_fqn: the log file of the ping sweep against the first hop after
%		the DSL link
%	up_Kbit: the uplink rate in Kilobits per second
%	down_Kbit: the downlink rate in Kilobits per second
%
% TODO:
%	find whether the carrier is ATM quantized
%	if yes:
%		what is the RTT step (try to deduce the combined up and down rates from this)
%		try to figure out the overhead for each packet
%	make sure all sizes are filled (use NaN for emty ones???)
%Thoughts:
%	maybe require to give the nominal up and down rates, to estimate the
%		RTT stepsize
%	ask about IPv4 or IPv6 (what about tunneling?)
%	the sweep should be taken directly connected to the modem to reduce
%		non-ATM routing delays

dbstop if error;
if ~(isoctave)
	timestamps.(mfilename).start = tic;
else
	tic();
end
disp(['Starting: ', mfilename]);

% control options
show_mean = 0;		% the means are noisier than the medians
show_median = 1;	% the mean is the way to go
show_min = 1;		% the min should be the best measure, but in the ATM test sweep is too variable
show_max = 0;		% only useful for debugging
show_sem = 0;		% give some estimate of the variance
show_ci = 1;		% show the confidence interval of the mean
ci_alpha = 0.05;	% alpha for confidence interval calculation
use_measure = 'median';
use_processed_results = 1;


if (nargin == 0)
	sweep_fqn = fullfile(pwd, 'ping_sweep_ATM.txt');	% was Bridged, LLC/SNAP RFC-1483/2684 connection (overhead 32 bytes - 14 = 18)
%  	sweep_fqn = fullfile(pwd, 'ping_sweep_CABLE_20120426_230227.txt');
%  	sweep_fqn = fullfile(pwd, 'ping_sweep_CABLE_20120801_001235.txt');
	up_Kbit = 512;
	down_Kbit = 3008;
end
if (nargin == 1)
	up_Kbit = 512;
	down_Kbit = 3008;
end
if (nargin == 2)
	down_Kbit = 3008;
end

%ATM
quantum.byte = 48;	% ATM packets are always 53 bytes, 48 thereof payload
quantum.bit = quantum.byte * 8;
ATM_cell.byte = 53;
ATM_cell.bit = ATM_cell.byte * 8;

% CONFIGURE THE NEXT TWO VALUES
line.down.Kbit = down_Kbit;
line.up.Kbit = up_Kbit;
% DONE
line.down.bit = line.down.Kbit * 1024;
line.up.bit = line.up.Kbit * 1024;


% known packet size offsets in bytes
offsets.IPv4 = 20;		% assume no IPv4 options are used, IP 6 would be 40bytes?
offsets.IPv6 = 40;		% not used yet...
offsets.ICMP = 8;		% ICMP header
offsets.ethernet = 14;	% ethernet header
% unknown offsets is what we need to figure out to feed tc-stab...


[sweep_dir, sweep_name] = fileparts(sweep_fqn);
cur_parsed_data_mat = [sweep_fqn(1:end-4), '.mat'];

if (use_processed_results && ~isempty(dir(cur_parsed_data_mat)))
	disp(['Loading processed ping data from ', cur_parsed_data_mat]);
	load(cur_parsed_data_mat, 'ping');
else
	% read in the result from a ping sweep
	disp(['Processing ping data from ', sweep_fqn]);
	ping = parse_ping_output(sweep_fqn);
	save(cur_parsed_data_mat, 'ping');
end


% analyze the data
min_ping_size = min(ping.data(:, ping.cols.size)) - offsets.ICMP;
disp(['Minimum size of ping payload used: ', num2str(min_ping_size), ' bytes.']);
known_overhead = offsets.ICMP + min_ping_size + offsets.IPv4;
ping.data(:, ping.cols.size) = ping.data(:, ping.cols.size) + known_overhead;	% we know we used ping so add the 8 bytes already (since these will not account for overhead)
size_list = unique(ping.data(:, ping.cols.size));

per_size.header = {'size', 'mean', 'median', 'min', 'max', 'std', 'n', 'sem', 'ci'};
per_size.cols = get_column_name_indices(per_size.header);
per_size.data = zeros([length(size_list), length(per_size.header)]);
per_size.data(:, per_size.cols.size) = size_list;
for i_size = 1 : length(size_list)
	cur_size = size_list(i_size);
	cur_size_idx = find(ping.data(:, ping.cols.size) == cur_size);
	per_size.data(i_size, per_size.cols.mean) = mean(ping.data(cur_size_idx, ping.cols.time));
	per_size.data(i_size, per_size.cols.median) = median(ping.data(cur_size_idx, ping.cols.time));
	per_size.data(i_size, per_size.cols.min) = min(ping.data(cur_size_idx, ping.cols.time));
	per_size.data(i_size, per_size.cols.max) = max(ping.data(cur_size_idx, ping.cols.time));
	per_size.data(i_size, per_size.cols.std) = std(ping.data(cur_size_idx, ping.cols.time), 0);
	per_size.data(i_size, per_size.cols.n) = length(cur_size_idx);
	per_size.data(i_size, per_size.cols.sem) = per_size.data(i_size, per_size.cols.std) / sqrt(length(cur_size_idx));	
	per_size.data(i_size, per_size.cols.ci) = calc_cihw(per_size.data(i_size, per_size.cols.std), per_size.data(i_size, per_size.cols.n), ci_alpha);	
end

figure('Name', sweep_name);
hold on;
legend_str = {};
if (show_mean)
	% means
	legend_str{end + 1} = 'mean';
	plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.mean), 'Color', [0 1 0 ]);
	if  (show_sem)
		legend_str{end + 1} = '+sem';
		legend_str{end + 1} = '-sem';
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.mean) - per_size.data(:, per_size.cols.sem), 'Color', [0 0.66 0]);
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.mean) + per_size.data(:, per_size.cols.sem), 'Color', [0 0.66 0]);
	end
	if  (show_ci)
		legend_str{end + 1} = '+ci';
		legend_str{end + 1} = '-ci';
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.mean) - per_size.data(:, per_size.cols.ci), 'Color', [0 0.37 0]);
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.mean) + per_size.data(:, per_size.cols.ci), 'Color', [0 0.37 0]);
	end

end
if(show_median)
	% median +- standard error of the mean, confidence interval would be
	% better
	legend_str{end + 1} = 'median';
	plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.median), 'Color', [1 0 0]);
	if (show_sem)
		legend_str{end + 1} = '+sem';
		legend_str{end + 1} = '-sem';
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.median) - per_size.data(:, per_size.cols.sem), 'Color', [0.66 0 0]);
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.median) + per_size.data(:, per_size.cols.sem), 'Color', [0.66 0 0]);
	end
	if (show_ci)
		legend_str{end + 1} = '+ci';
		legend_str{end + 1} = '-ci';
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.median) - per_size.data(:, per_size.cols.ci), 'Color', [0.37 0 0]);
		plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.median) + per_size.data(:, per_size.cols.ci), 'Color', [0.37 0 0]);
	end
end
if(show_min)
	% minimum, should be cleanest, but for the test data set looks quite sad...
	legend_str{end + 1} = 'min';
	plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.min), 'Color', [0 0 1]);
end
if(show_max)
	% minimum, should be cleanest, but for the test data set looks quite sad...
	legend_str{end + 1} = 'max';
	plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.max), 'Color', [0 0 0.66]);
end

title(['If this plot shows a (noisy) step function with a stepping of ', num2str(quantum.byte), ' bytes then the data carrier is quantised, make sure to use tc-stab']);
xlabel('Approximate packet size [bytes]');
ylabel('ICMP round trip times (ping RTT) [ms]');
legend(legend_str);
hold off;

% looking at the different plot we see many side peaks
%figure;
%plot(per_size.data(1:end - 1, strmatch('size', per_size.header , 'exact')), diff(per_size.data(:, strmatch('median', per_size.header , 'exact'))), 'Color', [1 0 0]);


% potentially clean up the data, by interpolating values with large sem
% from the neighbours?


% if this is ATM based the expectancy is there are RTT time quantums of 48
% bytes (the payload of an ATM package) find this
%x = fminsearch(@(x)sin(x^2), x0);


% estimate the RTT step size
% at ADSL down 3008kbit/sec up 512kbit/sec we expect, this oes not include
% processing time
expected_RTT_quantum_ms = (ATM_cell.bit / line.down.bit + ATM_cell.bit / line.up.bit ) * 1000;	% this estimate is rather a lower bound for fastpath , so search for best fits
disp(['lower bound estimate of one ATM cell RTT is ', num2str(expected_RTT_quantum_ms), ' ms.']);
% lets search from expected_RTT_quantum_ms to 1.5 * expected_RTT_quantum_ms
% in steps of expected_RTT_quantum_ms / 100

% to allow for interleaved ATM set ups increase the search space up to 32
% times best fastpath RTT estimate
RTT_quantum_list = (expected_RTT_quantum_ms : expected_RTT_quantum_ms / 100 : 32 * expected_RTT_quantum_ms);
quantum_list = (1 : 1 : quantum.byte);

differences = zeros([length(RTT_quantum_list) length(quantum_list)]);
cumulative_differences = differences;
all_stairs = zeros([length(RTT_quantum_list) length(quantum_list) length(per_size.data(:, per_size.cols.(use_measure)))]);
for i_RTT_quant = 1 : length(RTT_quantum_list)
	cur_RTT_quant = RTT_quantum_list(i_RTT_quant);
	for i_quant = 1 : quantum.byte
		[differences(i_RTT_quant, i_quant), cumulative_differences(i_RTT_quant, i_quant), all_stairs(i_RTT_quant, i_quant, :)] = get_difference_between_data_and_stair( per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.(use_measure)), ...
																	quantum_list(i_quant), quantum.byte, 0, cur_RTT_quant );
	end
end

% for the test DSL set the best x_offset is 21.
[min_cum_diff, min_cum_diff_idx] = min(cumulative_differences(:));
[min_cum_diff_row_idx, min_cum_diff_col_idx] = ind2sub(size(cumulative_differences),min_cum_diff_idx);
best_difference = differences(min_cum_diff_row_idx, min_cum_diff_col_idx);
disp(['remaining ATM cell length after ICMP header is ', num2str(quantum_list(min_cum_diff_col_idx)), ' bytes.']);
disp(['ICMP RTT of a single ATM packet is ', num2str(RTT_quantum_list(min_cum_diff_col_idx)), ' ms.']);

figure('Name', 'Comparing ping data with');
hold on
legend_str = {'ping_data', 'fitted_stair'};
plot(per_size.data(:, per_size.cols.size), per_size.data(:, per_size.cols.(use_measure)), 'Color', [1 0 0]);
plot(per_size.data(:, per_size.cols.size), squeeze(all_stairs(min_cum_diff_row_idx, min_cum_diff_col_idx, :)) + best_difference, 'Color', [0 1 0]);
title(['Estimated RTT per quantum: ', num2str(RTT_quantum_list(min_cum_diff_col_idx)), ' ms; ICMP header offset in quantum ', num2str(quantum_list(min_cum_diff_col_idx)), ' bytes']);
xlabel('Approximate packet size [bytes]');
ylabel('ICMP round trip times (ping RTT) [ms]');
if (isoctave)
	legend(legend_str);
else
	legend(legend_str, 'Interpreter', 'none');
end
hold off


% find cumulative_differences minimum, get the differences for that matrix
% and plot the corresponding stair with the real data
% next extrapolate to get the y-intercept


% as first approximation use the ATM cell offset and known offsets (ICMP
% IPv4 min_ping_size) to estimate the number of cells used for per paket
% overhead
% this assumes that no ATM related overhead is >= ATM cell size
% -1 to account for matlab 1 based indices
n_bytes_overhead_2nd_cell = quantum.byte - (quantum_list(min_cum_diff_col_idx) - 1);	% just assume we can not fit all overhead into one cell...
pre_IP_overhead = quantum.byte + (n_bytes_overhead_2nd_cell - known_overhead);	% ths is the one we are after in the end
disp(' ');
disp(['Estimated overhead preceeding the IP header: ', num2str(pre_IP_overhead), ' bytes']);



% use http://ace-host.stuart.id.au/russell/files/tc/tc-atm/ to present the
% most likely ATM setup for a given overhead and present a recommendation
% for the stab invocation
display_protocol_stack_information(pre_IP_overhead);



% now turn this into tc-stab recommendations:
disp(['Add the following to both the ingress and egress root qdiscs:']);
disp(' ');
disp(['A) Assuming the router connects over ethernet to the DSL-modem:']);
disp(['stab mtu 2048 tsize 128 overhead ', num2str(pre_IP_overhead - offsets.ethernet), ' linklayer atm']);
disp(' ');
disp(['B) Assuming the router connects via PPP and non-ethernet to the modem:']);
disp(['stab mtu 2048 tsize 128 overhead ', num2str(pre_IP_overhead), ' linklayer atm']);

disp(' ');
if ~(isoctave)
	timestamps.(mfilename).end = toc(timestamps.(mfilename).start);
	disp([mfilename, ' took: ', num2str(timestamps.(mfilename).end), ' seconds.']);
else
	toc
end

disp('Done...');

return
end


function [ ping_data ] = parse_ping_output( ping_log_fqn )
%PARSE_PING_OUTPUT read the putput of a ping run/sweep
% for further processng

cur_sweep_fd = fopen(ping_log_fqn, 'r');
ping_data.header = {'size', 'icmp_seq', 'ttl', 'time'};
ping_data.cols = get_column_name_indices(ping_data.header);

ping_data.data = zeros([1, length(ping_data.header)]);
cur_data_lines = 0;

% skip the first line
% PING netblock-75-79-143-1.dslextreme.com (75.79.143.1): (16 ... 1000)
% data bytes
header_line = fgetl(cur_sweep_fd);

while ~feof(cur_sweep_fd)
	cur_line = fgetl(cur_sweep_fd);
	[first_element, remainder] = strtok(cur_line);
	first_element_as_number = str2num(first_element);
	if isempty(first_element);
		% skip empty lines
		continue;
	end
	if strmatch('---', first_element)
		%we reached the end of sweeps
		break;
	end
	% now read in the data
	% 30 bytes from 75.79.143.1: icmp_seq=339 ttl=63 time=14.771 ms
	if ~isempty(first_element_as_number)
		cur_data_lines = cur_data_lines + 1;
		% size of the ICMP package
		ping_data.data(cur_data_lines, ping_data.cols.size) = first_element_as_number;
		% now process the remainder
		while ~isempty(remainder)
			[next_item, remainder] = strtok(remainder);
			equality_pos = strfind(next_item, '=');
			% data items are name+value pairs
			if ~isempty(equality_pos);
				cur_key = next_item(1: equality_pos - 1);
				cur_value = str2num(next_item(equality_pos + 1: end));
				switch cur_key
					% busybox ping and macosx ping return different key names
					case {'seq', 'icmp_seq'}
						ping_data.data(cur_data_lines, ping_data.cols.icmp_seq) = cur_value;
					case 'ttl'
						ping_data.data(cur_data_lines, ping_data.cols.ttl) = cur_value;
					case 'time'
						ping_data.data(cur_data_lines, ping_data.cols.time) = cur_value;
				end
			end
		end
		
	end
	
end

% clean up
fclose(cur_sweep_fd);
return
end



function [ difference , cumulative_difference, stair_y ] = get_difference_between_data_and_stair( data_x, data_y, x_size, stair_x_step_size, y_offset, stair_y_step_size )
% x_size is the flat part of the first stair, that is quantum minus the
% offset
debug = 0;
difference = [];
x_start_val = min(data_x);
x_end_val = max(data_x);

% construct stair
stair_x = data_x;
proto_stair_y = zeros([data_x(end) 1]);
% make sure the x_size values do not exceed the step size...
if (x_size > stair_x_step_size)
	if mod(x_size, stair_x_step_size) == 0
		x_size = stair_x_step_size;
	else
		x_size = mod(x_size, stair_x_step_size);
	end
end

stair_y_step_idx = (x_start_val + x_size : stair_x_step_size : x_end_val);
proto_stair_y(stair_y_step_idx) = stair_y_step_size;
stair_y = cumsum(proto_stair_y);
stair_y = stair_y(x_start_val:x_end_val) + y_offset;
if (debug)
	figure
	hold on;
	title(['x offset used: ', num2str(x_size), ' with quantum ', num2str(stair_x_step_size)]);
	plot(data_x, data_y, 'Color', [0 1 0]);
	plot(stair_x, stair_y, 'Color', [1 0 0]);
	hold off;
end
difference = sum(abs(data_y - stair_y)) / length(data_y);
cumulative_difference = sum(abs(data_y - (stair_y + difference)));

return
end

% function [ stair ] = build_stair(x_vector, x_size, stair_x_step_size, y_offset, stair_y_step_size )
% stair = [];
% 
% return
% end


function [columnnames_struct, n_fields] = get_column_name_indices(name_list)
% return a structure with each field for each member if the name_list cell
% array, giving the position in the name_list, then the columnnames_struct
% can serve as to address the columns, so the functions assitgning values
% to the columns do not have to care too much about the positions, and it
% becomes easy to add fields.
n_fields = length(name_list);
for i_col = 1 : length(name_list)
	cur_name = name_list{i_col};
	columnnames_struct.(cur_name) = i_col;
end
return
end




function [ci_halfwidth_vector] = calc_cihw(std_vector, n, alpha)
%calc_ci : calculate the half width of the confidence interval (for 1 - alpha)
%	the t_value lookup depends on alpha and the samplesize n; the relevant
%	calculation of the degree of freedom is performed inside calc_t_val.
%	ci_halfwidth = t_val(alpha, n-1) * std / sqrt(n)
%	Each groups CI ranges from mean - ci_halfwidth to mean - ci_halfwidth, so
%	the calling function has to perform this calculation...
%
% INPUTS:
%	std_vector: vector containing the standard deviations of all requested
%		groups
%	n: number of samples in each group, if the groups have different
%		samplesizes, specify each groups samplesize in a vector
%	alpha: the desired maximal uncertainty/error in the range of [0, 1]
% OUTPUT:
%	ci_halfwidth_vector: vector containing the confidence intervals half width
%		for each group

% calc_t_val return one sided t-values, for the desired two sidedness one has
% to half the alpha for the table lookup
cur_alpha = alpha / 2;

% if n is scalar use same n for all elements of std_vec
if isscalar(n)
	t_ci = calc_t_val(cur_alpha, n);
	ci_halfwidth_vector = std_vector * t_ci / sqrt(n);
	% if n is a vector, prepare a matching vector of t_ci values
elseif isvector(n)
	t_ci_vector = n;
	% this is probably ugly, but calc_t_val only accepts scalars.
	for i_pos = 1 : length(n)
		t_ci_vector(i_pos) = calc_t_val(cur_alpha, n(i_pos));
	end
	ci_halfwidth_vector = std_vector .* t_ci_vector ./ sqrt(n);
end

return
end

%-----------------------------------------------------------------------------
function [t_val] = calc_t_val(alpha, n)
% the t value for the given alpha and n
% so call with the n of the sample, not with degres of freedom
% see http://mathworld.wolfram.com/Studentst-Distribution.html for formulas
% return values follow Bortz, Statistik fuer Sozialwissenschaftler, Springer 
% 1999, table D page 775. That is it returns one sided t-values.
% primary author S. Moeller

% TODO:
%	sidedness of t-value???

% basic error checking
if nargin < 2
	error('alpha and n have to be specified...');
end

% probabilty of error
tmp_alpha = alpha ;%/ 2;
if (tmp_alpha < 0) || (tmp_alpha > 1)
	msgbox('alpha has to be taken from [0, 1]...');
	t_val = NaN;
	return
end
if tmp_alpha == 0
	t_val = -Inf;
	return
elseif tmp_alpha ==1
	t_val = Inf;
	return
end
% degree of freedom
df = n - 1;
if df < 1
	%msgbox('The n has to be >= 2 (=> df >= 1)...');
% 	disp('The n has to be >= 2 (=> df >= 1)...');
	t_val = NaN;
	return
end


% only calculate each (alpha, df) combination once, store the results
persistent t_val_array;
% create the t_val_array
if ~iscell(t_val_array)
	t_val_array = {[NaN;NaN]};
end
% search for the (alpha, df) tupel, avoid calculation if already stored
if iscell(t_val_array)
	% cell array of 2d arrays containing alpha / t_val pairs
	if df <= length(t_val_array)
		% test whether the required alpha, t_val tupel exists
		if ~isempty(t_val_array{df})
			% search for alpha
			tmp_array = t_val_array{df};
			alpha_index = find(tmp_array(1,:) == tmp_alpha);
			if any(alpha_index)
				t_val = tmp_array(2, alpha_index);
				return
			end
		end
	else
		% grow t_val_array to length of n
		missing_cols = df - length(t_val_array);
		for i_missing_cols = 1: missing_cols
			t_val_array{end + 1} = [NaN;NaN];
		end
	end
end

% check the sign
cdf_sign = 1;
if (1 - tmp_alpha) == 0.5
	t_val = t_cdf;
elseif (1 - tmp_alpha) < 0.5 % the t-cdf is point symmetric around (0, 0.5)
	cdf_sign = -1;
	tmp_alpha = 1 - tmp_alpha; % this will be undone later
end

% init some variables
n_iterations = 0;
delta_t = 1;
last_alpha = 1;
higher_t = 50;
lower_t = 0;
% find a t-value pair around the desired alpha value 
while norm_students_cdf(higher_t, df) < (1 - tmp_alpha);
	lower_t = higher_t;
	higher_t = higher_t * 2;
end

% search the t value for the given alpha...
while (n_iterations < 1000) && (abs(delta_t) >= 0.0001)
	n_iterations = n_iterations + 1;
	% get the test_t (TODO linear interpolation) 
 	% higher_alpha = norm_students_cdf(higher_t, df);
 	% lower_alpha = norm_students_cdf(lower_t, df);
	test_t = lower_t + ((higher_t - lower_t) / 2);
	cur_alpha = norm_students_cdf(test_t, df);
	% just in case we hit the right t spot on...
	if cur_alpha == (1 - tmp_alpha)
		t_crit = test_t;
		break;
	% probably we have to search for the right t
	elseif cur_alpha < (1 - tmp_alpha)
		% test_t is the new lower_t
		lower_t = test_t;
		%higher_t = higher_t;	% this stays as is...
	elseif cur_alpha > (1 - tmp_alpha)
		% 
		%lower_t = lower_t;	% this stays as is...
		higher_t = test_t;
	end
	delta_t = higher_t - lower_t;
	last_alpha = cur_alpha;
end
t_crit = test_t;

% set the return value, correct for negative t values
t_val = t_crit * cdf_sign;
if cdf_sign < 0
	tmp_alpha = 1 - tmp_alpha;
end

% store the alpha, n, t_val tupel in t_val_array
pos = size(t_val_array{df}, 2);
t_val_array{df}(1, (pos + 1)) = tmp_alpha;
t_val_array{df}(2, (pos + 1)) = t_val;

return
end

%-----------------------------------------------------------------------------
function [scaled_cdf] = norm_students_cdf(t, df)
% calculate the cdf of students distribution for a given degree of freedom df,
% and all given values of t, then normalize the result
% the extreme values depend on the values of df!!!

% get min and max by calculating values for extrem t-values (e.g. -10000000,
% 10000000)
extreme_cdf_vals = students_cdf([-10000000, 10000000], df);

tmp_cdf = students_cdf(t, df);

scaled_cdf =	(tmp_cdf - extreme_cdf_vals(1)) /...
				(extreme_cdf_vals(2) - extreme_cdf_vals(1));
return
end

%-----------------------------------------------------------------------------
function [cdf_value_array] = students_cdf(t_value_array, df)
%students_cdf: calc the cumulative density function for a t-distribution
% Calculate the CDF value for each value t of the input array
% see http://mathworld.wolfram.com/Studentst-Distribution.html for formulas
% INPUTS:	t_value_array:	array containing the t values for which to
%							calculate the cdf
%			df:	degree of freedom; equals n - 1 for the t-distribution

cdf_value_array = 0.5 +...
		((betainc(1, 0.5 * df, 0.5) / beta(0.5 * df, 0.5)) - ...
		(betainc((df ./ (df + t_value_array.^2)), 0.5 * df, 0.5) /...
		beta(0.5 * df, 0.5))) .*...
		sign(t_value_array);
	
return
end

%-----------------------------------------------------------------------------
function [t_prob_dist] = students_pf(df, t_arr)
%  calculate the probability function for students t-distribution

t_prob_dist =	(df ./ (df + t_arr.^2)).^((1 + df) / 2) /...
				(sqrt(df) * beta(0.5 * df, 0.5));

% % calculate and scale the cdf by hand...
% cdf = cumsum(t_prob_dist);
% discrete_t_cdf = (cdf - min(cdf)) / (max(cdf) - min(cdf));
% % numericaly get the t-value for the given alpha
% tmp_index = find(discrete_t_cdf > (1 - tmp_alpha));
% t_crit = t(tmp_index(1));
				
return
end

function in = isoctave ()
 persistent inout;

 if isempty(inout),
   inout = exist('OCTAVE_VERSION','builtin') ~= 0;
 end;
 in = inout;

return;
end


function [] = display_protocol_stack_information(pre_IP_overhead)
% use http://ace-host.stuart.id.au/russell/files/tc/tc-atm/ to present the
% most likely ATM prtocol stack setup for a given overhead so the user can
% compare with his prior knowledge

disp(' ');
disp('According to http://ace-host.stuart.id.au/russell/files/tc/tc-atm/');
disp(['', num2str(pre_IP_overhead), ' bytes overhead indicate']);

switch pre_IP_overhead
	case 8
		disp('Connection: IPoA, VC/Mux RFC-2684');
		disp('Protocol (bytes): ATM AAL5 SAR (8) : 8');
		
	case 16
		disp('Connection: IPoA, LLC/SNAP RFC-2684');
		disp('Protocol (bytes): ATM LLC (3), ATM SNAP (5), ATM AAL5 SAR (8) : 16');
		
	case 24
		disp('Connection: Bridged, VC/Mux RFC-1483/2684');
		disp('Protocol (bytes): Ethernet Header (14), ATM pad (2), ATM AAL5 SAR (8) : 24');

	case 28
		disp('Connection: Bridged, VC/Mux+FCS RFC-1483/2684');
		disp('Protocol (bytes): Ethernet Header (14), Ethernet PAD [8] (0), Ethernet Checksum (4), ATM pad (2), ATM AAL5 SAR (8) : 28');
	
	case 32
		disp('Connection: Bridged, LLC/SNAP RFC-1483/2684');
		disp('Protocol (bytes): Ethernet Header (14), ATM LLC (3), ATM SNAP (5), ATM pad (2), ATM AAL5 SAR (8) : 32');
		disp('OR');
		disp('Connection: PPPoE, VC/Mux RFC-2684');
		disp('Protocol (bytes): PPP (2), PPPoE (6), Ethernet Header (14), ATM pad (2), ATM AAL5 SAR (8) : 32');
	
	case 36
		disp('Connection: Bridged, LLC/SNAP+FCS RFC-1483/2684');
		disp('Protocol (bytes): Ethernet Header (14), Ethernet PAD [8] (0), Ethernet Checksum (4), ATM LLC (3), ATM SNAP (5), ATM pad (2), ATM AAL5 SAR (8) : 36');
		disp('OR');
		disp('Connection: PPPoE, VC/Mux+FCS RFC-2684');
		disp('Protocol (bytes): PPP (2), PPPoE (6), Ethernet Header (14), Ethernet PAD [8] (0), Ethernet Checksum (4), ATM pad (2), ATM AAL5 SAR (8) : 36');
	
	case 10
		disp('Connection: PPPoA, VC/Mux RFC-2364');
		disp('Protocol (bytes): PPP (2), ATM AAL5 SAR (8) : 10');
	
	case 14  
		disp('Connection: PPPoA, LLC RFC-2364');
		disp('Protocol (bytes): PPP (2), ATM LLC (3), ATM LLC-NLPID (1), ATM AAL5 SAR (8) : 14');
	
	case 40
		disp('Connection: PPPoE, LLC/SNAP RFC-2684');
		disp('Protocol (bytes): PPP (2), PPPoE (6), Ethernet Header (14), ATM LLC (3), ATM SNAP (5), ATM pad (2), ATM AAL5 SAR (8) : 40');
	
	case 44
		disp('Connection: PPPoE, LLC/SNAP+FCS RFC-2684');
		disp('Protocol (bytes): PPP (2), PPPoE (6), Ethernet Header (14), Ethernet PAD [8] (0), Ethernet Checksum (4), ATM LLC (3), ATM SNAP (5), ATM pad (2), ATM AAL5 SAR (8) : 44');  
	
	otherwise
		disp('a protocol stack this program does know (yet)...');
end

disp(' ');
return;
end


this should give you an estimate of the required overhead. 
	Now, I only tested it with the one data set I have available so it would be great if you could share the log file from the sweep (plus any information about the protocol stack in use for checking of the results…) Oh, the script is pretty rough but since all is pretty pedestrian it should not require to much time to understand/edit if necessary.

best
	sebastian