0){continue;}
list($ufx[$i],$ufy[$i],$ufz[$i])=explode(",",$incsvl);
$i++;
}
for($starti=0;$startifft($fx);
$wy = $fft->fft($fy);
$wz = $fft->fft($fz);//fft変換
for($i = 0; $i < 512; $i++){
if(128<$i){//周波数を割り出す
$f=$hz_i*(512-$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]->setReal(($wx[$i]->getReal())*$filter);//フィルタをかける
$wx[$i]->setImag(($wx[$i]->getImag())*$filter);
$wy[$i]->setReal(($wy[$i]->getReal())*$filter);
$wy[$i]->setImag(($wy[$i]->getImag())*$filter);
$wz[$i]->setReal(($wz[$i]->getReal())*$filter);
$wz[$i]->setImag(($wz[$i]->getImag())*$filter);
}
$wx = $fft->ifft($wx);
$wy = $fft->ifft($wy);
$wz = $fft->ifft($wz);//逆変換
//合成
$synthe=array();
for($i = 0; $i < 512; $i++){
$synthe[]=sqrt(
pow($wx[$i]->getReal(),2)+
pow($wy[$i]->getReal(),2)+
pow($wz[$i]->getReal(),2));
}
rsort($synthe);//強い順に並べる
$a=$synthe[30];//0.3秒の所を取り出す
$kshindo=floor(round(2*log10($a)+0.94,2)*10)/10;//計測震度計算
if($kshindo<0){$kshindo=0;}//0未満もある筈だが先例に倣って0以上にする
echo ($starti/100)." ".sprintf("%0.1f",$kshindo)."
\n";
//グラフ描く
$linex=$endi*(XIMGSIZE-1)/30000-1;
$liney=YIMGSIZE-$kshindo*(YIMGSIZE-1)/7-1;
imageline($im,$lastx,$lasty,$linex,$liney,$blk);
$lasty=$liney;
$lastx=$linex;
}
//画像出力
imagepng($im,"AA06EA01.png");
imagedestroy($im);
?>