読者です 読者をやめる 読者になる 読者になる

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

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

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)の周波数成分が確認できます。

PHPでフーリエ級数展開

フーリエ級数展開PHPで実装してみました。

一般の周期関数に対するフーリエ級数展開です。

<?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;     
			$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;
	}
}
?>

使い方

<?php
class TestCls extends FourierSeries{
        //周期関数
        function f($x){        
		return sin($x/2); //周期4π
	}
}
$N=100;
$from=0;
$to=4*M_PI;
	
$obj=new TestCls(10,4*M_PI);
for($x=$from;$x<$to;$x+=($to-$from)/$N){
	echo $obj->series($x)."\n";		
}
?>

NVelocityのテンプレートファイルを指定するときは相対パス

C#のテンプレートエンジンとしてCastle ProjectのNVelocityを使っていますが

ResourceNotFoundExceptionエラーが出てしまうので試行錯誤したところ、VelocityEngine.GetTemplateでテンプレートファイルのパスの指定を絶対パスから相対パスにすることで無事動きました。

 

Template template = velocity.GetTemplate( Environment.CurrentDirectory+ @"\template\template.vm"); //NG

                                                             

Template template = velocity.GetTemplate(  @"template\template.vm"); //OK