PHPで線スペクトル解析
前の記事でフーリエ級数を実装したので、今回は前の記事で使ったプログラムを
編集して線スペクトル解析を実装しました。
フーリエ変換についてはここで説明するより、
東北大学の解説サイトが詳しく書かれているので、そちらを読んでもらうのが良いと思います。
<? class FourierSeries{ //////////////// // $n:分割数 // $L: 周期 //////////////// function __construct($n,$L){ $this->n=$n; $this->L=$L/2; //フーリエ係数を計算 for($i=0;$i<=$n;$i++){ $this->a[$i]=$this->coefficientA($i); $this->b[$i]=$this->coefficientB($i); } } //$funcの$fromから$toまで積分.$nは分割数 function integration($n,$from,$to,$func,$arg){ $dx=($to-$from)/$n; $sum=0; for($i=0;$i<$n;$i++){ $x=$from+$i*$dx; // i番目の区間のx $fx=$this->$func($x,$arg); // 積分する関数 $sum+=$fx*$dx; } return $sum; } function fA($x,$n){ return $this->f($x)*cos($n*$x*M_PI/$this->L); } function fB($x,$n){ return $this->f($x)*sin($n*$x*M_PI/$this->L); } //フーリエ係数a function coefficientA($n){ return 1/$this->L*$this->integration(1000,-$this->L,$this->L,fA,$n); } //フーリエ係数b function coefficientB($n){ return 1/$this->L*$this->integration(1000,-$this->L,$this->L,fB,$n); } //フーリエ級数 public function series($x){ $sum=0; for($n=1;$n<=$this->n;$n++){ $sum+=$this->a[$n]*cos($n*$x*M_PI/$this->L)+$this->b[$n]*sin($n*$x*M_PI/$this->L); } return $this->a[0]/2+$sum; } function ImAbs($re,$im){ return sqrt(pow($re,2)+pow($im,2)); } //線スペクトル public function spectum($x){ for($n=1;$n<=$this->n;$n++){ echo $n." ".$this->ImAbs($this->a[$n],$this->b[$n])."\n"; } } } ?>
使い方
<? class hourAnalysis extends FourierSeries{ function f($x){ return sin(2*$x)+6*sin(4*$x)+2*cos(2*$x)+3*cos(3*$x); } } $hour=new hourAnalysis(10,2*M_PI); $hour->spectum(0); ?>
これを実行し、結果をgnuplotで表示すると
となり、f($x)の周波数成分が確認できます。