解析系エンジニアのブログ

数値解析の備忘録がメインですが、オープンソースの配布なども行っています。

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:id:mitsui725:20130123003259p:plain
となり、f($x)の周波数成分が確認できます。