#!/usr/bin/perl
#
use POSIX qw(floor);
my @tau;
my @adev;
my @mdev;
my @tdev;

# avar($tau,\@data) fills global arrays tau,adev,mdev,tdev
sub avar {
my $j;
my $k;
my $i;
my $s;
my $v;
my $u;
my @D;

my ($tau0,@y) = @_;
my $n = $#y;
my $N = floor(1+log(($n-1.0)/3.0)/log(2.0));
my $m = 1;

for ($j=0; $j<$N; $j++) {
print $j,"\n";
	$tau[$j] = $m*$tau0;
	$s = 0;
	$v = 0;
	for ($k=0;$k<$m;$k++) {
		$s = $s + $y[$k];
		$v = $v + $y[$k+$m];
		}
	$D[0] = ($v-$s)/$m;

	# $D(k)
	for ($k=1;$k<$n+1-2*$m;$k++) {
		$D[$k] = $D[$k-1]+($y[$k-1]-2*$y[$k+$m-1]+$y[$k+2*$m-1])/$m;
		}
	
	# ADEV
	$s = 0;
	for ($k=0;$k<$n+1-2*$m;$k++) {
		$v = $D[$k]*$D[$k];
		$s = $s + $v;
		}
	$adev[$j] = sqrt($s/($n+1-2*$m)/2);

	# MDEV
	$v = 0;
	for ($z=0;$z<$n+2-(3*$m);$z++) {
		$s = 0;
		for ($k=0;$k<$m;$k++) {
			$s = $s+$D[$z+$k];
			}
		$v = $v + $s * $s;
		}
	$v = $v/($n+2-(3*$m));
	$mdev[$j] = sqrt($v/2)/$m;

	# TDEV
	$tdev[$j] = $mdev[$j]*$m*$tau0/sqrt(3.0);

	$m = $m*2;

	}

}

# Add data points from stdin
while ($reading=<STDIN>) {
	if (substr($reading,0,1) =~ "[0-9]") {
		($x,$y) = split(/\s/,$reading);
		push(@xarray,$x);
		push(@yarray,$y);
		}
	}
print "Done reading...\n";
avar(300,@yarray);

print "\tTau\t\tADEV\t\tMDEV\t\tTDEV\n";
for ($q;$q<$#tau;$q++) {
	printf "\t%i\t%5.4E\t%5.4E\t%5.4E\n",
		$tau[$q],$adev[$q],$mdev[$q],$tdev[$q];
	}
