#!/usr/bin/perl
# 以下は変更する必要なし
#-----------------------------------------------------------------------------
#
# degs2allan (version  $Id: eodv,v 1.4 2007/05/22 02:22:33 Exp $) 
#
# Allan分散の算出
# 
# Usage:
# > degs2allan.pl hogehoge.csv 
#
# options:
#         特に無し
#
# example:
# degs2allan.pl hogehoge.csv                                #(1)
#
#-----------------------------------------------------------------------------
# Libraries
my @omg;
my @MDEV;
my @ADEV;

# tau[sec]
my $tau0= 0.1;
my @tau;
my @ave;
my @D;

my $n	= 	0;
my $m	=	1;
my $s	=	0;
my $v	=	0;

#-----------------------------------------------------------------------------
# Code begin
#-----------------------------------------------------------------------------
($opt) = @ARGV; 
if ($opt eq "-h")
{
    &show_usage();
    exit;
}

# オプションがない場合
if ($#ARGV == -1)
{
	print "Error:ファイル名を指定してください\n";
}

# オプションが１つある場合
if ($#ARGV ==  0)
{
	# 対象ファイルのオープン
	open(IN,"<$ARGV[0]") || die "Error: can't open $ARGV[0]\n";
		@line = <IN>;
	close(IN);
	
	$cnt = 0;
	
	# ファイル読込み(角速度読込み)
	foreach (@line)
	{
		# 角度[deg/s]
		$omg[$cnt] = (split(/\,/,$_))[0];
		
		# インクリメント処理
		$cnt++;
	}
	
	my $MAX	=	int(1+log10(($cnt - 1.0)/3.0)/log10(2.0));
	
	$m = 1;

	# Allan分散算出
	for($i=0;$i<$MAX;$i++)
	{
		# 1.tau計算
		$tau[$i] = $m*$tau0;
			
		# 2.差分値
		$s = 0;
		$v = 0;
		for($k=0;$k<$m;$k++)
		{
			$s = $s + $omg[$k];
			$v = $v + $omg[$k+$m];
		}
		$D[0] = ($v - $s)/$m;
			
		# 3.D算出
		for($k=1;$k<$cnt+1-2*$m;$k++)
		{
			$D[$k] = $D[$k-1] + ($omg[$k-1]-2*$omg[$k+$m-1]+$omg[$k+2*$m-1])/$m;
		}

		
		# 4.ADEV
		$s = 0;
		for($k=0;$k<$cnt+1-2*$m;$k++)
		{
			$v = $D[$k]**2;
			$s = $s + $v;
		}
		$ADEV[$i] = sqrt($s/($cnt+1-2*$m)/2);
		
		# 5.MDEV算出
		$v = 0;
		for($j=0;$j<$cnt+2-(3*$m);$j++)
		{
			$s = 0;
			for($k=0;$k<$m;$k++)
			{
				$s = $s + $D[$j+$k];
			}
			$v = $v + $s*$s;
		}
		$v = $v/($cnt + 2 -3*$m);
		$MDEV[$i] = sqrt($v/2)/$m;
		
		# 5.m値のインクリメント
		$m = $m*2;
	}
	
	# 出力
	for($j=0;$j<$MAX;$j++)
	{
		$disp = sprintf("%e,%e,%e\n",$tau[$j],$MDEV[$j],$ADEV[$j]);
		print $disp;
	}
}

# オプションが２つある場合
if ($#ARGV == 1)
{

}
exit;

#-----------------------------------------------------------------------------
# Sub code begin
#-----------------------------------------------------------------------------
sub show_usage()
{
  print "Usage:\n";
  print "> csv2allan [-h]\n";
  print "\n";
  print "options:\n";
  print "         -h : 利用方法を提示する\n";
  print "\n";
  print "example:\n";
  print "eodv -h                         #(1) HELP\n";
}

sub log10 {
    my $n = shift;
    return log($n)/log(10);
}
