<?php 
define("SAMPLEN",512);//サンプリングの個数

wiringPiSetup();//wiringPi初期化
for($i=0;$i<5;$i++){
	pinMode($i, 1);//0から4を出力に設定
}
ledout(0,0);//震度0表示で初期化

$fp=fopen("/dev/ttyACM0",'r+');//端末開く
fwrite($fp,'c');//開始コマンド
$laststr="";//受信文字列の最後の断片
$up=0;$uf=0;//upサンプルの現在位置 uf最大位置
$ref=array(12581177,12577495,14245523);//基準値
$u=array();//測定値をSAMPLEN個保存
$k=array();//3軸の測定値が揃うまで一時保存
$p=array(false,false,false);//3軸の測定値が揃ったかどうかのフラグ
$lt=microtime(true);//最後に震度を計算した時刻
while(1){
	$str=$laststr.fread($fp,512);//読み込み.断片を接続
	$laststr="";
	$lines=explode("\n",$str);
	foreach($lines as $line){
		if($line==""){continue;}
		$line=trim($line);
		if(strlen($line)!=7){$laststr=$line;continue;}//最後の断片を保存
		$v=substr($line,-1,1);//軸
		$s=substr($line,0,6);//測定値
		if($s=="800000"){continue;}//err?
		$k[$v]=hexdec($s);
		$p[$v]=true;//軸フラグセット
		if($p[0] && $p[1] && $p[2]){//揃った
			$k[0]=($k[0]-$ref[0])*980/($ref[2]-$ref[0]);//galに変換
			$k[1]=($k[1]-$ref[1])*980/($ref[2]-$ref[1]);
			$k[2]=($k[2]-($ref[0]+$ref[1])/2)*980/($ref[2]-($ref[0]+$ref[1])/2)-980;
			$up++;//現在位置を進める
			if($up>=SAMPLEN){$up=0;}//超えたら0に戻る
			$u[$up]=$k;//測定値を入れる
			$nt=microtime(true);//nt現在時刻
			if($uf<SAMPLEN){$uf++;$st=$nt-$lt;}//st:SAMPLEN個のサンプリングに要した秒数
			$p=array(false,false,false);//軸フラグクリア
		}
		if($uf>=SAMPLEN && $nt-$lt>1){
			$lt=$nt;//最後に震度計算した時刻更新
			for($i = 0; $i < SAMPLEN-$up; $i++){//ループになってる配列を整形する
				$fx[$i*2]=$u[$i+$up][0];$fx[$i*2+1]=0;
				$fy[$i*2]=$u[$i+$up][1];$fy[$i*2+1]=0;
				$fz[$i*2]=$u[$i+$up][2];$fz[$i*2+1]=0;
			}
			for($i = SAMPLEN-$up; $i < SAMPLEN; $i++){
				$fx[$i*2]=$u[$i+$up-SAMPLEN][0];$fx[$i*2+1]=0;
				$fy[$i*2]=$u[$i+$up-SAMPLEN][1];$fy[$i*2+1]=0;
				$fz[$i*2]=$u[$i+$up-SAMPLEN][2];$fz[$i*2+1]=0;
			}
			$kshindo=calkshindo($fx,$fy,$fz,$st);//震度計算
			echo "kshindo=$kshindo\n";
			if($kshindo<0.5){ledout(0,0);}
			if($kshindo>=0.5 && 1.5>$kshindo){ledout(1,0);}//led出力
			if($kshindo>=1.5 && 2.5>$kshindo){ledout(2,0);}
			if($kshindo>=2.5 && 3.5>$kshindo){ledout(3,0);}
			if($kshindo>=3.5 && 4.5>$kshindo){ledout(4,0);}
			if($kshindo>=4.5 && 5.0>$kshindo){ledout(5,1);}
			if($kshindo>=5.0 && 5.5>$kshindo){ledout(5,2);}
			if($kshindo>=5.5 && 6.0>$kshindo){ledout(6,1);}
			if($kshindo>=6.0 && 6.5>$kshindo){ledout(6,2);}
			if($kshindo>=6.5){ledout(7,0);}
		}
	}
}

function calkshindo($fx,$fy,$fz,$ftime){
	$hz_i=1/$ftime;//1番目の要素の周波数
	$avex=0;$avey=0;$avez=0;	//平均を減算
	for($i = 0; $i < SAMPLEN*2; $i+=2){
		$avex+=$fx[$i];
		$avey+=$fy[$i];
		$avez+=$fz[$i];
	}
	$avex=$avex/SAMPLEN;
	$avey=$avey/SAMPLEN;
	$avez=$avez/SAMPLEN;
	for($i = 0; $i < SAMPLEN*2; $i+=2){
		$fx[$i]-=$avex;
		$fy[$i]-=$avey;
		$fz[$i]-=$avez;
	}

	$wx =Fourier($fx, 1);//fft変換
	$wy =Fourier($fy, 1);
	$wz =Fourier($fz, 1);

	for($i = 0; $i < SAMPLEN; $i++){
		if(SAMPLEN/2<$i){//周波数を割り出す
			$f=$hz_i*(SAMPLEN-$i);//左右対称
		}else{
			$f=$hz_i*$i;
		}
		if($f==0){//直流は消す
			$filter=0;
		}else{
			$y=$f/10;//フィルタを計算
			$filter=sqrt(1/$f);
			$filter=$filter*(1/sqrt(1+ 0.694*pow($y,2) + 0.241*pow($y,4) + 0.0557*pow($y,6) + 0.009664*pow($y,8) + 0.00134*pow($y,10) + 0.000155*pow($y,12)));
			$filter=$filter*(sqrt(1-exp(- pow(($f/0.5),3) )));
		}
		$wx[$i*2]=$wx[$i*2]*$filter;//フィルタをかける
		$wx[$i*2+1]=$wx[$i*2+1]*$filter;
		$wy[$i*2]=$wy[$i*2]*$filter;
		$wy[$i*2+1]=$wy[$i*2+1]*$filter;
		$wz[$i*2]=$wz[$i*2]*$filter;
		$wz[$i*2+1]=$wz[$i*2+1]*$filter;
	}
	$wx =Fourier($wx, -1);//逆フーリエ変換
	$wy =Fourier($wy, -1);
	$wz =Fourier($wz, -1);

	$synthe=array();//合成
	for($i = 0; $i < SAMPLEN*2; $i+=2){
		$synthe[]=sqrt(
		pow($wx[$i],2)+
		pow($wy[$i],2)+
		pow($wz[$i],2));
	}
	rsort($synthe);//強い順に並べる
	$a=$synthe[(int)(0.3*SAMPLEN/$ftime)];//0.3秒の所を取り出す
	$kshindo=floor(round(2*log10($a)+0.94,2)*10)/10;//計測震度計算
	if($kshindo<0){$kshindo=0;}
	return $kshindo;
}

//$out 震度表示
//$kj 無0 弱1 強2
function ledout($out,$kj){
	digitalWrite(0, $out&0b001);
	digitalWrite(1, $out&0b010);
	digitalWrite(2, $out&0b100);
	digitalWrite(3, $kj!=1);
	digitalWrite(4, $kj!=2);
}

//https://gist.github.com/mbijon/1332348
###################################################################
# PHP_Fourier 0.03
# Original Fortran source by Numerical Recipies
# PHP port by Mathew Binkley (binkleym@nukote.com)
###################################################################
 
###################################################################
# Fourier($data, $sign) - Performs the FFT on the *complex*
# array $data
#
# Presumes that count($data) is an integer power of two
#
# $data[even] holds the real portion
# $data[odd] hold the imaginary portion
#
# Example: (1 + 2i) -> $data[0] = 1; $data[1] = 2;
#
# $sign = 1 performs the Fourier Transform
# $sign = -1 performs the Inverse Fourier Transform
#
# Use:
# $fourier_array = Fourier($inputarray, 1);
#
###################################################################
 
function Fourier($input, $isign) {
	#####################################################################
	# We need to shift the array up one because this script is a direct
	# port of the fortran program from NR. Should fix in future.
	#####################################################################
	$data[0] = 0;
	for ($i = 0; $i < count($input); $i++) $data[($i+1)] = $input[$i];

	$n = count($input);

	$j = 1;

	for ($i = 1; $i < $n; $i += 2) {
		if ($j > $i) {
			list($data[($j+0)], $data[($i+0)]) = array($data[($i+0)], $data[($j+0)]);
			list($data[($j+1)], $data[($i+1)]) = array($data[($i+1)], $data[($j+1)]);
		}

		$m = $n >> 1;

		while (($m >= 2) && ($j > $m)) {
			$j -= $m;
			$m = $m >> 1;
		}

		$j += $m;

	}

	$mmax = 2;
	while ($n > $mmax) { # Outer loop executed log2(nn) times
		$istep = $mmax << 1;

		$theta = $isign * 2*pi()/$mmax;

		$wtemp = sin(0.5 * $theta);
		$wpr = -2.0*$wtemp*$wtemp;
		$wpi = sin($theta);
		$wr = 1.0;
		$wi = 0.0;
		for ($m = 1; $m < $mmax; $m += 2) { # Here are the two nested inner loops
			for ($i = $m; $i <= $n; $i+= $istep) {

				$j = $i + $mmax;

				$tempr = $wr * $data[$j] - $wi * $data[($j+1)];
				$tempi = $wr * $data[($j+1)] + $wi * $data[$j];

				$data[$j] = $data[$i] - $tempr;
				$data[($j+1)] = $data[($i+1)] - $tempi;

				$data[$i] += $tempr;
				$data[($i+1)] += $tempi;

			}
			$wtemp = $wr;
			$wr = ($wr * $wpr) - ($wi * $wpi) + $wr;
			$wi = ($wi * $wpr) + ($wtemp * $wpi) + $wi;
		}
		$mmax = $istep;
	}

	for ($i = 1; $i < count($data); $i++) {
		$data[$i] *= sqrt(2/$n); # Normalize the data
		if (abs($data[$i]) < 1E-8) $data[$i] = 0; # Let's round small numbers to zero
		$input[($i-1)] = $data[$i]; # We need to shift array back (see beginning)
	}

	return $input;
}