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); ?>