#b-navbar { height:0px; display:none; visibility:hidden; }

ページ

2014年10月24日金曜日

Firefoxが起動しない!原因はハードウェアアクセラレーション


セーフモードだと起動するが、設定をクリアしても、アンインストール→再インストールしてもダメ。毎回クラッシュレポータ (ご迷惑をおかけして申し訳ありません) だったが、どうやら原因はハードウェアアクセラレーションだったようだ。


>> 一部のグラフィックカードとドライバの組み合わせによっては、ハードウェアアクセラレーションを利用すると、Firefox がクラッシュするか、ページ上のテキストやオブジェクトの表示に問題が起こることがあります。
  1. メニューボタン New Fx Menu をクリックし、オプション を選択します。
  2. 詳細 パネルを選択し、一般 タブを選択します。
  3. ハードウェアアクセラレーション機能を使用する (可能な場合) オプションをクリックしてチェックを外してください。
  4. メニューボタン New Fx Menu をクリックし、終了 Close 29 をクリックします。
  5. Firefox を再び起動します。

これで解決した。

当方: Windows7 64bit Firefox 33.0




2014年9月30日火曜日

Unix思想

Unix思想


  • モジュール化の原則 : クリーンなインターフェイスで結合される単純な部品を作れ。

"ぶざまな姿をさらさずに複雑なソフトウェアを書く唯一の方法は、全体としての複雑さの度合いを下げることだ""つまり、適切に定義されたインターフェイスで結び付けられた単純な部品からシステムを作り上げるのだ。 こうすれば、ほとんどの問題は局所化されるし、全体を壊さずに部品だけを改良することも不可能ではなくなる。"--『The Art of UNIX Programming』より

  • 明確性の原則 : 巧妙になるより明確であれ。

  • 組み立て部品の原則 : 他のプログラムと組み合わせられるように作れ。

  • 分離の原則 : メカニズムからポリシーを切り離せ。エンジンからインターフェイスを切り離せ。

  • 単純性の原則 : 単純になるように設計せよ。複雑な部分を追加するのは、どうしても必要なときだけに制限せよ。

  • 倹約の原則 : 他のものでは代えられないことが明確に実証されない限り、大きなプログラムを書くな。

  • 透明性の原則 : デバッグや調査が簡単になるように、わかりやすさを目指して設計せよ。

  • 安定性の原則 : 安定性は、透明性と単純性から生まれる。


すべては変化する
仕様が固まることは無い
技術も常に進化する
こだわるな
最初から全部やろうとしない
どこからやるか
「何が一番やばいですか?」
最も困っているところから
お金、個人情報、……
新機能開発から
バグ修正のところから



迷ったらシンプルな方を選ぶ
シンプルさは信頼性の前提である(Dijkstra)
simple と easy は異なる(Rich Hicky の講演より)
エレガント
"Elegance is a combination of power and simplicity."
エレガンスは、力と単純性が結合して生まれる。
"Elegant code is not only correct but visibly, transparently correct."
エレガントなコードは、ただ正しいだけではなく目に見える透明な形で正しい。
Keep It Simple and Small
Keep It Small Stupid!



http://twada.herokuapp.com/presentations/wewlc/wewlc.html より

2014年8月4日月曜日

R言語:nnet:plot.nn:線が出ないときの対処方法

ニューラルネットワークの解析結果を視覚化するplot.nnは素晴らしい。



参考:http://hosho.ees.hokudai.ac.jp/~kubo/ce/NeuralNetwork.html
  • nnet() で作った neural network の例 (作図は久保先生の自作関数 plot.nn())
    http://hosho.ees.hokudai.ac.jp/~kubo/log/2007/img07/plot.nn.txt



ただ、線が出ないときがある。ニューロンだけ表示されて、あとが真っ白という感じ。





11行目
col.w = function(w) ifelse(w > 0, "#ff400040", "#0000ff40"),

このカラーコードが8桁になっているためであろう。


そこで、
色をここから選んで持ってきて
http://html-color-codes.info/japanese/




修正した例

col.w = function(w) ifelse(w > 0, "#0000FF", "#FF0040"),





線が出た








keywords: R nnet grahical graph draw figures bug display line show 

余談ですが、線を太くしたいときは、10行目の w * 3 を w * 5 とか大きい数字に変えるとええで


2014年6月25日水曜日

OCR用のフォント、バーコード数字フォント




バーコードに使われている数字のフォントです。

有償版なら2万円ぐらい
http://www.flashbackj.com/ocr-b/
http://www.ricoh.co.jp/font/  など
エプソンLPユーザーなら
http://www.epson.jp/dl_soft/readme/7038.htm など


<無料>
非商用なら
≫Link:フォント http://ansuz.sooke.bc.ca/fonts-jp.php

海外のサイト(特に条件なし)OCR-A
http://sourceforge.net/projects/ocr-a-font/files/OCR-A/1.0/
OCR-AとB
http://sourceforge.jp/projects/tsukurimashou/releases/56948
上記のバックアップ
http://hp.vector.co.jp/authors/VA023120/data/OCRB.ttf
http://hp.vector.co.jp/authors/VA023120/data/OCRA.ttf

時系列データのシンプルな予測公式

近傍での1次式を仮定するなら、

  • \(f \left(n+1\right) =2f \left(n\right) -f \left(n-1\right) \)


で、次の値を近似できる。
今日と昨日の値から、明日のデータ値を予測する感じです。


2次式を仮定するなら、
今日と昨日とおとといの値から、明日を予測する感じになる。

  • \(f \left(n+1\right) =3f \left(n\right) -3f \left(n-1\right) +f \left(n-2\right) \)


上記で、次の値を近似できるものの・・・精度は下がる気がする。(結局一番最初の式がシンプルでベストなことがほとんどかもしれない)

以降、同様に順次漸化式を適用すれば


  • \(f \left(n+1\right) =4f \left(n\right) -6f \left(n-1\right) +4f \left(n-2\right) -f \left(n-3\right) \)
  • \(f \left(n+1\right) =5f \left(n\right) -10f \left(n-1\right) +10f \left(n-2\right) -5f \left(n-3\right) +f \left(n-4\right) \)

係数は二項定理っぽくなりますが、あまり実用性はなさそう。


超単純なモデルによる外挿なので
誤差が多いことや長期予測には向かないことに注意しつつ
サクッと1点外挿したいときにどうぞ。

離散データをもっと正確にするには、
回帰分析、重回帰分析、ARIMAなどを検討しましょう。



2014年6月22日日曜日

PHP:単回帰分析(係数と重決定R2値)(1次関数~n次関数)

まず1次関数での近似を行った。

出典:http://PolynomialRegression.drque.net/
最新版は上記から取得してください。


PolynomialRegression.php のソース
class PolynomialRegression
{
  private $xPowers;
  private $xyPowers;
  private $numberOfCoefficient;
  private $forcedValue;

  public function __construct( $numberOfCoefficient = 3 )
  {
    $this->numberOfCoefficient = $numberOfCoefficient;
    $this->reset();

  } // __construct

  public function reset()
  {
    $this->forcedValue = array();
    $this->xPowers = array();
    $this->xyPowers = array();

    $squares = ( $this->numberOfCoefficient - 1 ) * 2;

    // Initialize power arrays.
    for ( $index = 0; $index <= $squares; ++$index )
    {
      $this->xPowers[ $index ] = 0;
      $this->xyPowers[ $index ] = 0;
    }

  } // reset

  public function setDegree( $numberOfCoefficient )
  {
    $this->numberOfCoefficient = $numberOfCoefficient;

  } // setDegree

  public function setNumberOfCoefficient( $numberOfCoefficient )
  {
    $this->numberOfCoefficient = $numberOfCoefficient;

  } // setNumberOfCoefficient

  public function getNumberOfCoefficient( $numberOfCoefficient )
  {
    return $this->numberOfCoefficient;

  } // getnumberOfCoefficient

  public function setForcedCoefficient( $coefficient, $value )
  {
    $this->forcedValue[ $coefficient ] = $value;

  } // setForcedCoefficient

  public function getForcedCoefficient( $coefficient, $value )
  {
    $result = null;
    if ( isset( $this->forcedValue[ $coefficient ] ) )
      $result = $this->forcedValue[ $coefficient ];

    return $result;

  } // getForcedCoefficient

  public function addData( $x, $y )
  {
    $squares = ( $this->numberOfCoefficient - 1 ) * 2;

    // Remove the effect of the forced coefficient from this value.
    foreach ( $this->forcedValue as $coefficient => $value )
    {
      $sub = bcpow( $x, $coefficient );
      $sub = bcmul( $sub, $value );
      $y = bcsub( $y, $sub );
    }

    // Accumulate new data to power sums.
    for ( $index = 0; $index <= $squares; ++$index )
    {
      $this->xPowers[ $index ] =
        bcadd( $this->xPowers[ $index ], bcpow( $x, $index ) );

      $this->xyPowers[ $index ] =
        bcadd
        (
          $this->xyPowers[ $index ],
          bcmul( $y, bcpow( $x, $index ) )
        );
    }

  } // addData


  public function getCoefficients( $numberOfCoefficient = -1 )
  {
    // If no number of coefficients specified, use standard.
    if ( $numberOfCoefficient == -1 )
      $numberOfCoefficient = $this->numberOfCoefficient;

    $matrix = array();
    for ( $row = 0; $row < $numberOfCoefficient; ++$row )
    {
      $matrix[ $row ] = array();
      for ( $column = 0; $column < $numberOfCoefficient; ++$column )
        $matrix[ $row ][ $column ] =
          $this->xPowers[ $row + $column ];
    }

    // Create augmented matrix by adding X*Y powers.
    for ( $row = 0; $row < $numberOfCoefficient; ++$row )
      $matrix[ $row ][ $numberOfCoefficient ] = $this->xyPowers[ $row ];

    foreach ( $this->forcedValue as $coefficient => $value )
    {
      for ( $index = 0; $index < $numberOfCoefficient; ++$index )
      {
        $matrix[ $index ][ $coefficient ] = "0";
        $matrix[ $coefficient ][ $index ] = "0";
      }

      $matrix[ $coefficient ][ $coefficient ] = "1";
      $matrix[ $coefficient ][ $numberOfCoefficient ]      = $value;
    }

    // Determine number of rows in matrix.
    $rows = count( $matrix );

    // Initialize done.
    $isDone = array();
    for ( $column = 0; $column < $rows; ++$column )
      $isDone[ $column ] = false;

    $order = array();
    for ( $column = 0; $column < $rows; ++$column )
    {
      // Find a row to work with.
      // A row that has a term in this column, and has not yet been
      // reduced.
      $activeRow = 0;
      while ( ( ( 0 == $matrix[ $activeRow ][ $column ] )
             || ( $isDone[ $activeRow ] ) )
           && ( $activeRow < $rows ) )
      {
        ++$activeRow;
      }

      // Do we have a term in this row?
      if ( $activeRow < $rows )
      {
        // Remember the order.
        $order[ $column ] = $activeRow;

        // Normalize row--results in the first term being 1.
        $firstTerm = $matrix[ $activeRow ][ $column ];
        for ( $subColumn = $column; $subColumn <= $rows; ++$subColumn )
          $matrix[ $activeRow ][ $subColumn ] =
            bcdiv( $matrix[ $activeRow ][ $subColumn ], $firstTerm );

        // This row is finished.
        $isDone[ $activeRow ] = true;

        // Subtract the active row from all rows that are not finished.
        for ( $row = 0; $row < $rows; ++$row )
          if ( ( ! $isDone[ $row ] )
            && ( 0 != $matrix[ $row ][ $column ] ) )
          {
             // Get first term in row.
             $firstTerm = $matrix[ $row ][ $column ];
             for ( $subColumn = $column; $subColumn <= $rows; ++$subColumn )
             {
               $accumulator = bcmul( $firstTerm, $matrix[ $activeRow ][ $subColumn ] );
               $matrix[ $row ][ $subColumn ] =
                 bcsub( $matrix[ $row ][ $subColumn ], $accumulator );
             }
          }
      }
    }

    // Reset done.
    for ( $row = 0; $row < $rows; ++$row )
     $isDone[ $row ] = false;

    $coefficients = array();

    for ( $column = ( $rows - 1 ); $column >= 0; --$column )
    {
      // The active row is based on order.
      $activeRow = $order[ $column ];

      // The active row is now finished.
      $isDone[ $activeRow ] = true;

      // For all rows not finished...
      for ( $row = 0; $row < $rows; ++$row )
        if ( ! $isDone[ $row ] )
        {
          $firstTerm = $matrix[ $row ][ $column ];

          // Back substitution.
          for ( $subColumn = $column; $subColumn <= $rows; ++$subColumn )
          {
            $accumulator =
              bcmul( $firstTerm, $matrix[ $activeRow ][ $subColumn ] );
            $matrix[ $row ][ $subColumn ] =
              bcsub( $matrix[ $row ][ $subColumn ], $accumulator );
          }
        }

      // Save this coefficient for the return.
      $coefficients[ $column ] = $matrix[ $activeRow ][ $rows ];
    }

    // Coefficients are stored backward, so sort them.
    ksort( $coefficients );

    // Return the coefficients.
    return $coefficients;

  } // getCoefficients


  static public function interpolate( $coefficients, $x )
  {
    $numberOfCoefficient = count( $coefficients );

    $y = 0;
    for ( $coefficentIndex = 0; $coefficentIndex < $numberOfCoefficient; ++$coefficentIndex )
    {
      // y += coefficients[ coefficentIndex ] * x^coefficentIndex
      $y =
        bcadd
        (
          $y,
          bcmul
          (
            $coefficients[ $coefficentIndex ],
            bcpow( $x, $coefficentIndex )
          )
        );
    }

    return floatval( $y );

  } // interpolate

} // Class

テスト実行

データは、配列で指定する。 中の array(1,0),・・・が、array(Xの値,Yの値(目的変数))の順である。

require_once( 'inc_PolynomialRegression.php' ); //class読込

//テスト配列定義

 $data = array (
  array(1,0),
  array(2,4),
  array(3,3),
  array(4,2),
  array(5,5),
  array(6,3),
  array(7,8)
 );
 

//print_r ( $data ) ;

// Precision digits in BC math. 
  bcscale( 10 ); 

  // Start a regression class of order 2--linear regression. 
  $PolynomialRegression = new PolynomialRegression( 2 ); //変数の数、次数+1

  // Add all the data to the regression analysis. 
  foreach ( $data as $dataPoint ) 
    $PolynomialRegression->addData( $dataPoint[ 0 ], $dataPoint[ 1 ] ); 

  // Get coefficients for the polynomial. 
  $coefficients = $PolynomialRegression->getCoefficients(); 

  // 
  // Get average of Y-data. 
  // 
  $Y_Average = 0.0; 
  foreach ( $data as $dataPoint ) 
    $Y_Average += $dataPoint[ 1 ]; 

  $Y_Average /= count( $data ); 

  // 
  // Calculate R Squared. 
  // 

  $Y_MeanSum  = 0.0; 
  $Y_ErrorSum = 0.0; 
  foreach ( $data as $dataPoint ) 
  { 
    $x = $dataPoint[ 0 ]; 
    $y = $dataPoint[ 1 ]; 
    $error  = $y; 
    $error -= $PolynomialRegression->interpolate( $coefficients, $x ); 
    $Y_ErrorSum += $error * $error; 

    $error  = $y; 
    $error -= $Y_Average; 
    $Y_MeanSum += $error * $error; 
  } 

 $R_Squared = 1.0 - ( $Y_ErrorSum / $Y_MeanSum ); 
 
 
 
 
 

// Print slope and intercept of linear regression. 
// 四捨五入 4桁目

$para_a = round( $coefficients[ 1 ], 4 );
$para_b = round( $coefficients[ 0 ], 4 );
$para_R = round(  $R_Squared ,4 );



//結果出力
print "$para_a,$para_b,$para_R" ;





出力結果
0.8571,0.1429,0.5455

よって
y=0.8571 x + 0.1429
重決定 R2 は 0.5455


※X,Yの入れ替えに注意



エクセルとの比較:OK





二次関数で近似するときは

$PolynomialRegression = new PolynomialRegression( 3 );//2→3にする

$para_a = round( $coefficients[ 2 ], 4 );//増やす
$para_b = round( $coefficients[ 1 ], 4 );
$para_c = round( $coefficients[ 0 ], 4 );
$para_R = round(  $R_Squared ,4 );

出力:0.0952,0.0952,1.2857,0.5657
エクセルと比較:OK

/*=========================================================================*/ /* Name: PolynomialRegression.php */ /* Uses: Calculates and returns coefficients for polynomial regression. */ /* Date: 06/01/2009 */ /* Author: Andrew Que (http://www.DrQue.net/) */ /* Revisions: */ /* 0.8 - 06/01/2009- QUE - Creation. */ /* 0.9 - 06/14/2012- QUE - */ /* + Bug fix: removed notice causes by uninitialized variable. */ /* + Converted naming convention. */ /* + Fix spelling errors (or the ones I found). */ /* + Changed to row-echelon method for solving matrix which is much */ /* faster than the determinant method. */ /* 0.91 - 05/17/2013- QUE - */ /* = Changed name to Polynonial regression as this is more fitting to */ /* to the function. */ /* 0.92 - 12/28/2013- QUE - */ /* + Added forced offset. */ /* 1.00 - 12/29/2013 - QUE - */ /* + Forced offset changed to allow any term to be forced. */ /* Unit complete. Correlation coefficient (r-squared) has been */ /* implemented externally in the demos. */ /* 1.1 - 2014/05/05 - QUE - */ /* + 'interpolate' is now static as it does not need an instance to */ /* operate. Useful if coefficients have been calculated elsewhere. */ /* - Deprecated 'setDegree' function. This is the wrong terminology for */ /* what the function does. It actually sets the number of */ /* coefficients for the polynomial. The degree of the polynomial is */ /* the number of coefficients less one. Made the identical function */ /* 'setNumberOfCoefficient' to replace it. */ /* + Added getter functions for anything that has a set function. */ /* */ /* This project is maintained at: */ /* http://PolynomialRegression.drque.net/ */ /* */ /* ----------------------------------------------------------------------- */ /* */ /* Polynomial regression class. */ /* Copyright (C) 2009, 2012-2014 Andrew Que */ /* */ /* This program is free software: you can redistribute it and/or modify */ /* it under the terms of the GNU General Public License as published by */ /* the Free Software Foundation, either version 3 of the License, or */ /* (at your option) any later version. */ /* */ /* This program is distributed in the hope that it will be useful, */ /* but WITHOUT ANY WARRANTY; without even the implied warranty of */ /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */ /* GNU General Public License for more details. */ /* */ /* You should have received a copy of the GNU General Public License */ /* along with this program. If not, see . */ /* */ /* ----------------------------------------------------------------------- */ /* */ /* (C) Copyright 2009, 2012-2014 */ /* Andrew Que */ /* > */ /*=========================================================================*/ /** * Polynomial regression. * *

* Used for calculating polynomial regression coefficients. Useful for * linear and non-linear regression, and polynomial curve fitting. * * @package PolynomialRegression * @author Andrew Que ({@link http://www.DrQue.net/}) * @link http://PolynomialRegression.drque.net/ Project home page. * @copyright Copyright (c) 2009, 2012-2014, Andrew Que * @license http://opensource.org/licenses/gpl-license.php GNU Public License * @version 1.1 */ /** * Used for calculating polynomial regression coefficients and interpolation using * those coefficients. Useful for linear and non-linear regression, and polynomial * curve fitting. * * Note: Requires BC math to be compiled into PHP. Higher-degree polynomials end up * with very large/small numbers, requiring an arbitrary precision arithmetic. Make sure * to set "bcscale" as coefficients will likely have decimal values. * * Quick example of using this unit to calculate linear regression (1st degree polynomial): * * * $regression = new PolynomialRegression( 2 ); * // ... * $regression->addData( $x, $y ); * // ... * $coefficients = $regression->getCoefficients(); * // ... * $y = $regression->interpolate( $coefficients, $x ); * * * * @package PolynomialRegression * @link http://PolynomialRegression.drque.net/ Project home page.

PHPで相関係数(行列)を計算する

参考:http://lib.stat.cmu.edu/multi/pca.c


/** 
* Creates a correlation matrix from a given data matrix.  
* 
* @see http://lib.stat.cmu.edu/multi/pca.c 
* 
* @author   Paul Meagher 
* @version  0.2 
* @modified Mar 23, 2008 
*/ 

   
     
class CorrelationMatrix { 
     
    private $eps    = 0.005; 

    public $nrows   = 0; 
    public $ncols   = 0;     
     
    public $means   = array(); 
    public $stddevs = array();     
    public $cormat  = array();     
     
  function CorrelationMatrix($data) { 
                 
    // num rows 
    $this->nrows = count($data);    
     
    // num cols 
    $this->ncols = count($data[0]);  
     
    // Determine mean of column vectors of input data matrix  
    for ($j=0; $j < $this->ncols; $j++) { 
      $this->means[$j] = 0.0; 
      for($i=0; $i < $this->nrows; $i++)  
        $this->means[$j] += $data[$i][$j]; 
      $this->means[$j] /= $this->nrows; 
    } 
           
    // Determine standard deviations of column vectors of data matrix.  
    for ($j=0; $j < $this->ncols; $j++) { 
      $this->stddevs[$j] = 0.0; 
      for ($i=0; $i < $this->nrows; $i++)  
        $this->stddevs[$j] += (($data[$i][$j]-$this->means[$j])*($data[$i][$j]-$this->means[$j])); 
      $this->stddevs[$j] /= $this->nrows; 
      $this->stddevs[$j] = sqrt($this->stddevs[$j]); 
      // The following in an inelegant but usual way to handle 
      // near-zero stddev values, which would cause a zero- 
      // divide error.  
      if ($this->stddevs[$j] <= $this->eps)  
        $this->stddevs[$j] = 1.0; 
    } 
     
    // Center and reduce the column vectors.  
    for ($i=0; $i < $this->nrows; $i++) { 
      for ($j=0; $j < $this->ncols; $j++) { 
        $data[$i][$j] -= $this->means[$j]; 
        $x = sqrt($this->nrows); 
        $x *= $this->stddevs[$j]; 
        $data[$i][$j] /= $x; 
      } 
    } 
   
    // Calculate the m * m correlation matrix.  
    for ($j1=0; $j1 < $this->ncols-1; $j1++) { 
      $this->cormat[$j1][$j1] = 1.0; 
      for ($j2=$j1+1; $j2 < $this->ncols; $j2++) { 
        $this->cormat[$j1][$j2] = 0.0; 
        for ($i=0; $i < $this->nrows; $i++) 
          $this->cormat[$j1][$j2] += ( $data[$i][$j1] * $data[$i][$j2]); 
        $this->cormat[$j2][$j1] = $this->cormat[$j1][$j2]; 
      } 
    } 
     
    $this->cormat[$this->ncols-1][$this->ncols-1] = 1.0; 
     
  } 
   
  function printMeans() {   
    printf("Column Means: 
"); 
    for ($j=0; $j < $this->ncols; $j++)   
      printf("%7.2f", $this->means[$j]);   
    printf("
");   
  } 

  function printStdDevs() { 
    printf("Column Standard Deviations: 
"); 
    for ($j=0; $j < $this->ncols; $j++)  
      printf("%7.2f", $this->stddevs[$j]);  
    printf("
"); 
  } 

  function printCorMat() { 
    printf("Correlation Matrix: 
"); 
    for ($i=0; $i < $this->ncols; $i++) {      
      for ($j=0; $j < $this->ncols; $j++)  
        //echo $this->cormat[$i][$j]." "; 
        printf("%7.4f", $this->cormat[$i][$j]);  
      printf("
"); 
    } 
  } 
       
} 


動作デモ

include "CorrelationMatrix.php";
 
$data[0] = array(1, 2, 3);
$data[1] = array(2, 3, 4);
$data[2] = array(3, 4, 5);
$data[3] = array(4, 6, 8);
$data[4] = array(5, 8, 10);
 
$cm = new CorrelationMatrix($data);
print "
";
$cm->printMeans();
print "
"; $cm->printStdDevs(); print "
"; $cm->printCorMat(); print "
"; print_r ( $cm->cormat ) ; print "
";

デモ出力


Column Means: 
   3.00   4.60   6.00

Column Standard Deviations: 
   1.41   2.15   2.61
Correlation Matrix: 
 1.0000 0.9848 0.9762
 0.9848 1.0000 0.9970
 0.9762 0.9970 1.0000

Array
(
    [0] => Array
        (
            [0] => 1
            [1] => 0.98479824644792
            [2] => 0.97618706018395
        )

    [1] => Array
        (
            [0] => 0.98479824644792
            [1] => 1
            [2] => 0.9969527608178
        )

    [2] => Array
        (
            [0] => 0.97618706018395
            [1] => 0.9969527608178
            [2] => 1
        )

)


10行目の、$cm = new CorrelationMatrix($data); の配列に
結果数値が入って帰ってくる。

2014年6月19日木曜日

PHPで重回帰分析

サイトはこちら:
Class for computing multiple linear regression of the form
http://mnshankar.wordpress.com/2011/05/01/regression-analysis-using-php/
Copyright (c)  2011 Shankar Manamalkav <nshankarあっとまーくufl.edu>

------------------------------------------------------------
ソースはこちら:classライブラリです。includeして使ってください。
・class Lib_Regression
・class Lib_Matrix の2つ
後半に
・テストデータ
・サンプルコード
・サンプルからの出力を掲載。
・エクセルとの比較検証もしました。
------------------------------------------------------------


class Lib_Regression
{

    protected $SSEScalar; //sum of squares due to error
    protected $SSRScalar; //sum of squares due to regression
    protected $SSTOScalar; //Total sum of squares
    protected $RSquare;         //R square
    protected $RSquarePValue;   //R square p falue (significant f)
    protected $F;               //F statistic
    protected $coefficients;    //regression coefficients array
    protected $stderrors;    //standard errror array
    protected $tstats;     //t statistics array
    protected $pvalues;     //p values array
    protected $x = array();
    protected $y = array();

    //implement singleton
    /**
     * An instance of this class (singleton)
     * @var Lib_Regression 
     */
    public function setX($x)
    {
        $this->x = $x;
    }

    public function setY($y)
    {
        $this->y = $y;
    }

    public function getSSE()
    {
        return $this->SSEScalar;
    }

    public function getSSR()
    {
        return $this->SSRScalar;
    }

    public function getSSTO()
    {
        return $this->SSTOScalar;
    }

    public function getRSQUARE()
    {
        return $this->RSquare;
    }
    public function getRSQUAREPValue()
    {
        return $this->RSquarePValue;
    }

    public function getF()
    {
        return $this->F;
    }

    public function getCoefficients()
    {
        return $this->coefficients;
    }

    public function getStandardError()
    {
        return $this->stderrors;
    }

    public function getTStats()
    {
        return $this->tstats;
    }

    public function getPValues()
    {
        return $this->pvalues;
    }

    /**
     * @example $reg->loadCSV('abc.csv',array(0), array(1,2,3));
     * @param type $file
     * @param array $xcolnumbers
     * @param type $ycolnumber 
     */
    public function LoadCSV($file, array $ycol, array $xcol, $hasHeader=true)
    {
        $xarray = array();
        $yarray = array();
        $handle = fopen($file, "r");

        //if first row has headers.. ignore
        if ($hasHeader)
            $data = fgetcsv($handle);
        //get the data into array
        while (($data = fgetcsv($handle)) !== FALSE)
        {
            $rawData[] = array($data);
        }
        $sampleSize = count($rawData);

        $r = 0;
        while ($r < $sampleSize)
        {
            $xarray[] = $this->GetArray($rawData, $xcol, $r, true);
            $yarray[] = $this->GetArray($rawData, $ycol, $r);   //y always has 1 col!
            $r++;
        }
        $this->x = $xarray;
        $this->y = $yarray;
    }

    private function GetArray($rawData, $cols, $r, $incIntercept=false)
    {
        $arrIdx = "";
        $z = 0;
        $arr = array();
        if ($incIntercept)
            //prepend an all 1's column for the intercept
            $arr[]=1;
        foreach ($cols as $key => $val)
        {
            $arrIdx = '$rawData[' . $r . '][0][' . $val . '];';
            eval("\$z = $arrIdx");
            $arr[] = $z;
        }
        return $arr;
    }
    //if true is passed to this method, compute normalized coefficients
    public function Compute($computeNormalized=false)
    {
        if ((count($this->x)==0)||(count($this->y)==0))
                throw new Exception ('Please supply valid X and Y arrays');
        $mx = new Lib_Matrix($this->x);
        $my = new Lib_Matrix($this->y);

        //coefficient(b) = (X'X)-1X'Y 
        $xTx = $mx->Transpose()->Multiply($mx)->Inverse();
        $xTy = $mx->Transpose()->Multiply($my);

        $coeff = $xTx->Multiply($xTy);

        $num_independent = $mx->NumColumns();   //note: intercept is included
        $sample_size = $mx->NumRows();
        $dfTotal = $sample_size - 1;
        $dfModel = $num_independent - 1;
        $dfResidual = $dfTotal - $dfModel;
        //create unit vector..
        for ($ctr = 0; $ctr < $sample_size; $ctr++)
            $u[] = array(1);

        $um = new Lib_Matrix($u);
        //SSR = b(t)X(t)Y - (Y(t)UU(T)Y)/n        
        //MSE = SSE/(df)
        $SSR = $coeff->Transpose()->Multiply($mx->Transpose())->Multiply($my)
                ->Subtract(
                        ($my->Transpose()
                        ->Multiply($um)
                        ->Multiply($um->Transpose())
                        ->Multiply($my)
                        ->ScalarDivide($sample_size))
        );

        $SSE = $my->Transpose()->Multiply($my)->Subtract(
                        $coeff->Transpose()
                                ->Multiply($mx->Transpose())
                                ->Multiply($my)
        );

        $SSTO = $SSR->Add($SSE);
        $this->SSEScalar = $SSE->GetElementAt(0, 0);
        $this->SSRScalar = $SSR->GetElementAt(0, 0);
        $this->SSTOScalar = $SSTO->GetElementAt(0, 0);

        $this->RSquare = $this->SSRScalar / $this->SSTOScalar;


        $this->F = (($this->SSRScalar / ($dfModel)) / ($this->SSEScalar / ($dfResidual)));
	
	$this->RSquarePValue= $this->FishersF($this->F, $dfModel, $dfResidual);
	
        $MSE = $SSE->ScalarDivide($dfResidual);
        //this is a scalar.. get element
        $e = ($MSE->GetElementAt(0, 0));

        $stdErr = $xTx->ScalarMultiply($e);
        for ($i = 0; $i < $num_independent; $i++)
        {
            //get the diagonal elements
            $searray[] = array(sqrt($stdErr->GetElementAt($i, $i)));
            //compute the t-statistic
            $tstat[] = array($coeff->GetElementAt($i, 0) / $searray[$i][0]);
            //compute the student p-value from the t-stat
            $pvalue[] = array($this->Student_PValue($tstat[$i][0], $dfResidual));
        }
        //convert into 1-d vectors and store
        for ($ctr = 0; $ctr < $num_independent; $ctr++)
        {
            $this->coefficients[] = $coeff->GetElementAt($ctr, 0);
            $this->stderrors[] = $searray[$ctr][0];
            $this->tstats[] = $tstat[$ctr][0];
            $this->pvalues[] = $pvalue[$ctr][0];
        }
	/* multiplying the unstandardized coefficient by the ratio of the standard deviations 
         * of the independent variable (x) and dependent variable(y) gives us NORMALIZED coefft.
         */
        if ($computeNormalized)
        {
            //if normalized betas are desired, compute them
            $this->coefficients = $this->ComputeNormalizedBeta($this->coefficients, 
                    $mx->GetInnerArray(), $my->GetInnerArray());
        }
    }

    
    private function ComputeNormalizedBeta($beta, $x, $y)
    {
        $betaNormalized=array();
        $cols = count($x[0]);
        for ($i = 1; $i < $cols; $i++)
        {
            ${'xvar' . $i} = $this->ComputeStdev($x, $i);
        }
        $yvar1 = $this->ComputeStdev($y, 0);
        $betaNormalized[0] = 0;  //Normalized coeff. intercept is not defined
        for ($i = 1; $i < count($beta); $i++)
        {
            $betaNormalized[] = $beta[$i] * ${'xvar' . $i} / $yvar1;
        }
        return $betaNormalized;
    }

    private function ComputeStdev($arr, $column)
    {
        //$arr is a 2d matrix
        //compute stdev of a particular column of a 2d array
        //standardized x and y.
        $rows = count($arr);
        
        for ($i = 0; $i < $rows; $i++)
        {
            $x[] = $arr[$i][$column];
        }
        $stDevX = $this->StdDeviation($x);
        return $stDevX;
    }

    //compute the standard deviation of 1-d array that is passed in
    private function StdDeviation($sample)
    {
        //compute mean
        $mean = array_sum($sample) / sizeof($sample);
        //loop through all and find square of difference betn each and mean
        foreach ($sample as $num)
        {
            $devs[] = pow(($num - $mean), 2);
        }
        //finally, take square root of average of those..
        $std = sqrt(array_sum($devs) / count($devs));
        return $std;
    }
    /**
     * @link http://home.ubalt.edu/ntsbarsh/Business-stat/otherapplets/pvalues.htm#rtdist
     * @param float $t_stat
     * @param float $deg_F
     * @return float 
     */
    private function Student_PValue($t_stat, $deg_F)
    {
        $t_stat = abs($t_stat);
        $mw = $t_stat / sqrt($deg_F);
        $th = atan2($mw, 1);
        if ($deg_F == 1)
            return 1.0 - $th / (M_PI / 2.0);
        $sth = sin($th);
        $cth = cos($th);
        if ($deg_F % 2 == 1)
            return 1.0 - ($th + $sth * $cth * $this->statcom($cth * $cth, 2, $deg_F - 3, -1)) / (M_PI / 2.0);
        else
            return 1.0 - ($sth * $this->statcom($cth * $cth, 1, $deg_F - 3, -1));
    }

    /**
     * @link http://home.ubalt.edu/ntsbarsh/Business-stat/otherapplets/pvalues.htm#rtdist
     * @param float $q
     * @param float $i
     * @param float $j
     * @param float $b
     * @return float 
     */
    private function statcom($q, $i, $j, $b)
    {
        $zz = 1;
        $z = $zz;
        $k = $i;
        while ($k <= $j)
        {
            $zz = $zz * $q * $k / ( $k - $b);
            $z = $z + $zz;
            $k = $k + 2;
        }
        return $z;
    }
    /**
     *Formula to return Fishers F (P Value of RSquare) (Significance F)
     * @param type $f - F statistic computed from regression analysis
     * @param float $n1 - Degrees of freedom1 (regression)
     * @param int $n2 - Degrees of freedom2 (residual)
     * @return float 
     */
    private function  FishersF($f,$n1, $n2)
    {
        $x = $n2/($n1*$f+$n2);
        if ($n1 % 2 == 0)
            return ($this->statcom(1-$x,$n2,$n1+$n2-4,$n2-2)*pow($x,$n2/2));
        if ($n2 % 2 == 0)
            return (1 - $this->statcom($x,$n1,$n1+$n2-4,$n1-2)*pow(1-$x,$n1/2));
        $th = atan(sqrt($n1*$f/$n2));
        $a = $th/(M_PI / 2.0);
        $sth = sin($th);
        $cth = cos($th);
        if ($n2>1)
            $a+=$sth*$cth*($this->statcom($cth*$cth, 2, $n2-3, -1)/(M_PI / 2.0));
        if ($n1==1)
            return 1-$a;
        $c = 4*$this->statcom($sth*$sth, $n2+1, $n1+$n2-4, $n2-2)*$sth*pow($cth,$n2)/M_PI;
        if ($n2==1)
            return (1-$a+$c/2);
        $k=2;
        while ($k<=($n2-1)/2)
        {
            $c = $c*$k/($k-.5);
            $k++;
        }
        return (1-$a+$c);
    }

}





class Lib_Matrix
{

    //global vars
    protected $rows;
    protected $columns;
    protected $MainMatrix = array();

    /**
     * Matrix Constructor
     * 
     * Initialize the Matrix object. Throw an exception if jagged array is passed.
     *
     * @param array $matrix - The array
     */
    function __construct($matrix)
    {
        for ($i = 0; $i < count($matrix); $i++)
        {
            for ($j = 0; $j < count($matrix[$i]); $j++)
                $this->MainMatrix[$i][$j] = $matrix[$i][$j];
        }
        $this->rows = count($this->MainMatrix);
        $this->columns = count($this->MainMatrix[0]);
        if (!$this->isValidMatrix())
        {
            throw new Exception("Invalid matrix");
        }
    }

    /**
     * Is it a valid matrix?
     * 
     * Returns 'False' if it is not a rectangular matrix
     *
     * @return bool
     */
    private function isValidMatrix()
    {
        for ($i = 0; $i < $this->rows; $i++)
        {
            $numCol = count($this->MainMatrix [$i]);
            if ($this->columns != $numCol)
                return false;
        }
        return true;
    }

    /**
     * Display the matrix
     * Formatted display of matrix for debugging.
     */
    public function DisplayMatrix()
    {
        $rows = $this->rows;
        $cols = $this->columns;
        echo "Order of the matrix is ($rows rows X $cols columns)\n";
        for ($r = 0; $r < $rows; $r++)
        {
            for ($c = 0; $c < $cols; $c++)
            {
                echo $this->MainMatrix[$r][$c];
            }
            echo "\n";
        }
    }

    /**
     * Get the inner array stored in matrix object
     * 
     * @return array 
     */
    public function GetInnerArray()
    {
        return $this->MainMatrix;
    }

    /**
     * Number of rows in the matrix
     * @return integer 
     */
    public function NumRows()
    {
        return count($this->MainMatrix);
    }

    /**
     * Number of columns in the matrix
     * @return integer
     */
    public function NumColumns()
    {
        return count($this->MainMatrix[0]);
    }

    /**
     * Return element found at location $row, $col.
     * 
     * @param int $row
     * @param int $col
     * @return object(depends on input)
     */
    public function GetElementAt($row, $col)
    {
        return $this->MainMatrix[$row][$col];
    }

    /**
     * Is this a square matrix?
     * 
     * Determinants and inverses only exist for square matrices!
     * 
     * @return bool 
     */
    public function isSquareMatrix()
    {
        if ($this->rows == $this->columns)
            return true;

        return false;
    }

    /**
     * Subtract matrix2 from matrix object on which this method is called
     * @param Lib_Matrix $matrix2
     * @return Lib_Matrix Note that original matrix is left unchanged
     */
    public function Subtract(Lib_Matrix $matrix2)
    {
        $rows1 = $this->rows;
        $columns1 = $this->columns;

        $rows2 = $matrix2->NumRows();
        $columns2 = $matrix2->NumColumns();

        if (($rows1 != $rows2) || ($columns1 != $columns2))
            throw new Exception('Matrices are not the same size!');

        for ($i = 0; $i < $rows1; $i++)
        {
            for ($j = 0; $j < $columns1; $j++)
            {
                $newMatrix[$i][$j] = $this->MainMatrix[$i][$j] -
                        $matrix2->GetElementAt($i, $j);
            }
        }
        return new Lib_Matrix($newMatrix);
    }

    /**
     * Add matrix2 to matrix object that calls this method.
     * @param Model_Matrix $matrix2
     * @return Lib_Matrix Note that original matrix is left unchanged
     */
    function Add(Lib_Matrix $matrix2)
    {
        $rows1 = $this->rows;
        $rows2 = $matrix2->NumRows();
        $columns1 = $this->columns;
        $columns2 = $matrix2->NumColumns();
        if (($rows1 != $rows2) || ($columns1 != $columns2))
            throw new Exception('Matrices are not the same size!');

        for ($i = 0; $i < $rows1; $i++)
        {
            for ($j = 0; $j < $columns1; $j++)
            {
                $newMatrix[$i][$j] = $this->MainMatrix[$i][$j] +
                        $matrix2->GetElementAt($i, $j);
            }
        }
        return new Lib_Matrix($newMatrix);
    }

    /**
     * Multiply matrix2 into matrix object that calls this method
     * @param Model_Matrix $matrix2
     * @return Lib_Matrix Note that original matrix is left unaltered
     */
    function Multiply(Lib_Matrix $matrix2)
    {
        $sum = 0;
        $rows1 = $this->rows;
        $columns1 = $this->columns;

        $columns2 = $matrix2->NumColumns();
        $rows2 = $matrix2->NumRows();
        if ($columns1 != $rows2)
            throw new Exception('Incompatible matrix types supplied');
        for ($i = 0; $i < $rows1; $i++)
        {
            for ($j = 0; $j < $columns2; $j++)
            {
                $newMatrix[$i][$j] = 0;
                for ($ctr = 0; $ctr < $columns1; $ctr++)
                {
                    $newMatrix[$i][$j] += $this->MainMatrix[$i][$ctr] *
                            $matrix2->GetElementAt($ctr, $j);
                }
            }
        }
        return new Lib_Matrix($newMatrix);
    }

    /**
     * Multiply every element of matrix on which this method is called by the scalar
     * @param object $scalar
     * @return Lib_Matrix 
     */
    public function ScalarMultiply($scalar)
    {
        $rows = $this->rows;
        $columns = $this->columns;

        $newMatrix = array();
        for ($i = 0; $i < $rows; $i++)
        {
            for ($j = 0; $j < $columns; $j++)
            {
                $newMatrix[$i][$j] = $this->MainMatrix[$i][$j] * $scalar;
            }
        }
        return new Lib_Matrix($newMatrix);
    }

    /**
     * Divide every element of matrix on which this method is called by the scalar
     * @param object $scalar
     * @return Lib_Matrix 
     */
    public function ScalarDivide($scalar)
    {
        $rows = $this->rows;
        $columns = $this->columns;

        $newMatrix = array();
        for ($i = 0; $i < $rows; $i++)
        {
            for ($j = 0; $j < $columns; $j++)
            {
                $newMatrix[$i][$j] = $this->MainMatrix[$i][$j] / $scalar;
            }
        }
        return new Lib_Matrix($newMatrix);
    }

    /**
     * Return the sub-matrix after crossing out the $crossx and $crossy row and column respectively
     * Part of determinant expansion by minors method
     * @param int $crossX
     * @param int $crossY
     * @return Lib_Matrix 
     */
    public function GetSubMatrix($crossX, $crossY)
    {
        $rows = $this->rows;
        $columns = $this->columns;

        $newMatrix = array();
        $p = 0; // submatrix row counter
        for ($i = 0; $i < $rows; $i++)
        {
            $q = 0; // submatrix col counter
            if ($crossX != $i)
            {
                for ($j = 0; $j < $columns; $j++)
                {
                    if ($crossY != $j)
                    {
                        $newMatrix[$p][$q] = $this->GetElementAt($i, $j);
                        //$matrix[$i][$j];
                        $q++;
                    }
                }
                $p++;
            }
        }
        return new Lib_Matrix($newMatrix);
    }

    /**
     * Compute the determinant of the square matrix on which this method is called
     * @link http://mathworld.wolfram.com/DeterminantExpansionbyMinors.html
     * @return object(depends on input)
     */
    public function Determinant()
    {
        if (!$this->isSquareMatrix())
            throw new Exception("Not a square matrix!");
        $rows = $this->rows;
        $columns = $this->columns;
        $determinant = 0;
        if ($rows == 1 && $columns == 1)
        {
            return $this->MainMatrix[0][0];
        }
        if ($rows == 2 && $columns == 2)
        {
            $determinant = $this->MainMatrix[0][0] * $this->MainMatrix[1][1] -
                    $this->MainMatrix[0][1] * $this->MainMatrix[1][0];
        } else
        {
            for ($j = 0; $j < $columns; $j++)
            {
                $subMatrix = $this->GetSubMatrix(0, $j);
                if (fmod($j, 2) == 0)
                {
                    $determinant += $this->MainMatrix[0][$j] * $subMatrix->Determinant();
                } else
                {
                    $determinant -= $this->MainMatrix[0][$j] * $subMatrix->Determinant();
                }
            }
        }
        return $determinant;
    }

    /**
     * Compute the transpose of matrix on which this method is called (invert rows and columns)
     * @return Lib_Matrix 
     */
    public function Transpose()
    {
        $rows = $this->rows;
        $columns = $this->columns;
        $newArray = array();
        for ($i = 0; $i < $rows; $i++)
        {
            for ($j = 0; $j < $columns; $j++)
            {
                $newArray[$j][$i] = $this->MainMatrix[$i][$j];
            }
        }
        return new Lib_Matrix($newArray);
    }

    /**
     * Compute the inverse of the matrix on which this method is found (A*A(-1)=I)
     * (cofactor(a))T/(det a)
     * @link http://www.mathwords.com/i/inverse_of_a_matrix.htm
     * @return Lib_Matrix 
     */
    function Inverse()
    {
        if (!$this->isSquareMatrix())
            throw new Exception("Not a square matrix!");
        $rows = $this->rows;
        $columns = $this->columns;

        $newMatrix = array();
        for ($i = 0; $i < $rows; $i++)
        {
            for ($j = 0; $j < $columns; $j++)
            {
                $subMatrix = $this->GetSubMatrix($i, $j);
                if (fmod($i + $j, 2) == 0)
                {
                    $newMatrix[$i][$j] = ($subMatrix->Determinant());
                } else
                {
                    $newMatrix[$i][$j] = -($subMatrix->Determinant());
                }
            }
        }
        $cofactorMatrix = new Lib_Matrix($newMatrix);
        return $cofactorMatrix->Transpose()
                ->ScalarDivide($this->Determinant());
    }
}










----------------------------------------------------------------------
テストデータ



   y   x0   x1   x2
4.5 1 8 2
22.5 1 40.5 24.5
2 1 4.5 0.5
0.5 1 0.5 2
18 1 4.5 4.5
2 1 7 8
32 1 24.5 40.5
4.5 1 4.5 2
40.5 1 32 24.5
2 1 0.5 4.5

----------------------------------------------------------------------
利用サンプル

Y目的変数、X配列:説明変数

      //independent vars.. note 1 has been added to first column for computing
        //intercept
       $x = array(
            array(1, 8, 2,2),
            array(1, 40.5, 24.5,3),
            array(1, 4.5, .5,4),
            array(1, .5, 2,5),
            array(1, 4.5, 4.5,6),
            array(1, 7, 8,7),
            array(1, 24.5, 40.5,8),
            array(1, 4.5, 2,9),
            array(1, 32, 24.5,11),
            array(1, .5, 4.5,22),
        );
        //dependent matrix..note it is a 2d array
        $y = array(array(4.5),
            array(22.5),
            array(2),
            array(.5),
            array(18),
            array(2),
            array(32),
            array(4.5),
            array(40.5),
            array(2)
        );
        
 
$reg = new Lib_Regression();
$reg->setX($x);
$reg->setY($y);
 
//NOTE: passing true to the compute method generates standardized coefficients
$reg->Compute();    //go!
 
var_dump($reg->getSSE());
var_dump($reg->getSSR());
var_dump($reg->getSSTO());
var_dump($reg->getRSQUARE());
var_dump($reg->getF());
var_dump($reg->getRSQUAREPValue());
var_dump($reg->getCoefficients());
var_dump($reg->getStandardError());
var_dump($reg->getTStats());
var_dump($reg->getPValues());






----------------------------------------------------------------------
サンプルからの出力

float(447.80833871098)
float(1448.216661289)
float(1896.025)
float(0.76381728157014)//重決定R2
float(11.319035123602)
float(0.0064027774890735)

array(3) {//係数
[0]=>
float(1.5647582746097)//定数項 
[1]=>
float(0.37870241945753)//x1
[2]=>
float(0.57474832913739)//x2
}//(X0の係数は無しというか、定数項になる)


array(3) {
[0]=>
float(3.4900561380127)
[1]=>
float(0.32797487084168)
[2]=>
float(0.34303873326455)
}

array(3) {
[0]=>
float(0.44834759463231)
[1]=>
float(1.1546690101155)
[2]=>
float(1.6754619038724)
}

array(3) {
[0]=>
float(0.6674511497586)
[1]=>
float(0.28611718855617)
[2]=>
float(0.13775128117929)
}

--------------------------------------------------------------------
エクセルとの結果比較:OK




回帰統計
重相関 R 0.873966
重決定 R2 0.763817
補正 R2 0.553479
標準誤差 7.998289
観測数 10
分散分析表
  自由度 変動 分散 観測された分散比
回帰 3 1448.217 482.7389 11.31904
残差 7 447.8083 63.97262
合計 10 1896.025    
  係数 標準誤差 P-値
切片 1.564758  3.490056  0.448348  0.667451
x0 0 0 65535 #NUM!
x1 0.378702 0.327975 1.154669  0.286117
x2 0.574748 0.343039 1.675462  0.137751


--------------------------------------------------------------------
コメント:

従属変数(X)で、冒頭1列が余計に必要な点に注意がいるが(1で埋めるだけ)、
あとは、スムーズに使える。

変数の数は、適宜増やせる。
処理も速い。

エクセルだと、従属変数が16個までなので、これを使えば、突破できる。
R使ってもいいのだろうが、ウェブシステムに簡易的に組み込むならこっちも検討する価値はある。



/**
 * Copyright (c)  2011 Shankar Manamalkav <nshankarあっとまーくufl.edu>
 * 
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 * 
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 *
 * Class for computing multiple linear regression of the form
 * y=a+b1x1+b2x2+b3x3...
 *
 * @author shankar<nshankarあっとまーくufl.edu>
 * 
 * 
 */


2014年6月10日火曜日

フーリエ解析、フーリエ級数をPHP関数で使う


逆変換もできる関数が揃っているリンクはこちら
http://www.phpkode.com/scripts/item/fast-fourier-transform/

ソースは下記ですが、最新版は上記から探してみてください。


FFT.class.php

<?php
/**
 * @author Michele Andreoli <michi.andreoliあっとまーくgmail.com>
 * @name FFT.class.php
 * @version 0.5 updated 06-07-2010
 * @license http://opensource.org/licenses/gpl-license-php GNU Public License
 * @package FFT
 */
require_once 'Complex.class.php';

/**
 * Class that calculate the FFT and the inverse FFT of a 1D signal
 */
class FFT {
private $dim;
private $p;
private $ind;
private $func;
private $w1;
private $w1i;
private $w2;

/**
 * Constructor for FFT class
 * @param int $dim dimension of the signal
 */
public function __construct($dim) {
$this->dim = $dim;
$this->p = log($this->dim, 2);
}

/**
 * Calculate the FFT of a signal
 * @param array<double> $func input signal
 * @return array<complex> return the DFT for the signal in input
 */
public function fft($func) {
$this->func = $func;

for ($i = 0; $i < $this->dim; $i++)
$this->w1[$i] = new Complex($func[$i], 0);

$w[0] = new Complex(1, 0);
$w[1] = new Complex(cos((-2 * M_PI) / $this->dim), sin((-2 * M_PI) / $this->dim));

for ($i = 2; $i < $this->dim; $i++)
$w[$i] = Complex::Cmul($w[$i-1], $w[1]);

return $this->calculate($w);
}

/**
 * Calculate the inverse FFT of a signal
 * @param array<complex> $func input signal
 * @return array<complex> result of inverse FFT for the signal in input
 */
public function ifft($func) {
$this->func = $func;
$norm = 1 / $this->dim;

for ($i = 0; $i < $this->dim; $i++)
$this->w1[$i] = new Complex($func[$i]->getReal(), $func[$i]->getImag());

$w[0] = new Complex(1, 0);
$w[1] = new Complex(cos((2 * M_PI) / $this->dim), sin((2 * M_PI) / $this->dim));

for ($i = 2; $i < $this->dim; $i++)
$w[$i] = Complex::Cmul($w[$i-1], $w[1]);

$this->w1i = $this->calculate($w);

for ($i = 0; $i < $this->dim; $i++)
$this->w1i[$i] = Complex::RCmul($norm, $this->w1i[$i]);

return $this->w1i;
}

private function calculate($w) {
$k = 1;
$ind[0] = 0;

for ($j = 0; $j < $this->p; $j++) {
for ($i = 0; $i < $k; $i++) {
$ind[$i] *= 2;
$ind[$i+$k] = $ind[$i] + 1;
}
$k *= 2;
}

for ($i = 0; $i < $this->p; $i++) {
$indw = 0;
for ($j = 0; $j < pow(2, $i); $j++) {
$inf = ($this->dim / pow(2, $i)) * $j;
$sup = (($this->dim / pow(2, $i)) * ($j+1)) - 1;
$comp = ($this->dim / pow(2, $i)) / 2;

for ($k = $inf; $k <= floor($inf+(($sup-$inf)/2)); $k++)
     $this->w2[$k] = Complex::Cadd(Complex::Cmul($this->w1[$k], $w[0]), Complex::Cmul($this->w1[$k+$comp], $w[$ind[$indw]]));

$indw++;

for ($k = floor($inf+(($sup-$inf)/2)+1); $k <= $sup; $k++)
     $this->w2[$k] = Complex::Cadd(Complex::Cmul($this->w1[$k], $w[$ind[$indw]]), Complex::Cmul($this->w1[$k-$comp], $w[0]));

$indw++;
}

for($j = 0; $j < $this->dim; $j++)
      $this->w1[$j] = $this->w2[$j];
}

for ($i = 0; $i < $this->dim; $i++)
$this->w1[$i] = $this->w2[$ind[$i]];

return $this->w1;
}

/**
 * Getter for the FFT
 * @return array<complex> get the FFT of the signal
 */
public function getFFT() {
return $this->w1;
}

/**
 * Getter for the inverse FFT
 * @return array<complex> get the inverse FFT of the signal
 */
public function getIFFT() {
return $this->w1i;
}

/**
 * Get the absolute value of the signal
 * @param array<complex> $fft
 * @return array<complex> get the absolute value of FFT
 */
public function getAbsFFT($w) {
for ($i = 0; $i < $this->dim; $i++)
$temp[$i] = Complex::Cabs($w[$i]);

return $temp;
}

/**
 * Getter for the dimension of the signal
 * @return int return the dimension of the signal
 */
public function getDim() {
return $this->dim;
}

/**
 * Convert an array of double into an array of complex
 * @param array<double> $func
 * @return array<complex> return an array of complex
 */
public function doubleToComplex($func) {
for ($i = 0; $i < count($func); $i++)
$aux[$i] = new Complex($func[$i], 0);

return $aux;
}

/**
 * Convert an array of complex into an array of double
 * @param array<complex> $func
 * @return array<double> return an array of double
 */
public function complexToDouble($func) {
for ($i = 0; $i < count($func); $i++)
$aux[$i] = $func[$i]->getReal();

return $aux;
}
}
?>




Complex.class.php
<?php
/**
 * @author Michele Andreoli <michi.andreoliあっとまーくgmail.com>
 * @name Complex.class.php
 * @version 0.3 updated 09-05-2010
 * @license http://opensource.org/licenses/gpl-license-php GNU Public License
 * @package FFT
 */

/**
 * Class that implements a complex data number
 */
class Complex {
private $real;
private $imag;

/**
 * Create a complex datatype
 * @param double $real
 * @param double $imag
 */
public function __construct($real, $imag) {
$this->real = $real;
$this->imag = $imag;
}

/**
 * Return the multiplication of two complex number
 * @param complex $a
 * @param complex $b
 * @return complex
 */
static public function Cmul($a, $b) {
$c = new Complex(0, 0);

$c->setReal($a->getReal() * $b->getReal() - $a->getImag() * $b->getImag());
$c->setImag($a->getImag() * $b->getReal() + $a->getReal() * $b->getImag());

return $c;
}

/**
 * Return the sum of two complex number
 * @param complex $a
 * @param complex $b
 * @return complex
 */
static public function Cadd($a, $b) {
$c = new Complex(0, 0);

$c->setReal($a->getReal() + $b->getReal());
$c->setImag($a->getImag() + $b->getImag());

return $c;
}

/**
 * Return the difference of two complex number
 * @param complex $a
 * @param complex $b
 * @return complex
 */
static public function Csub($a, $b) {
$c = new Complex(0, 0);
      
$c->setReal($a->getReal() - $b->getReal());
$c->setImag($a->getImag() - $b->getImag());
      
        return $c;
}

/**
 * Return the absolute value of a complex number
 * @param complex $z
 * @return double
 */
static public function Cabs($z) {
        $x = abs($z->getReal());
        $y = abs($z->getImag());
      
        if ($x == 0.0)
        $ans = $y;
        else if ($y == 0.0)
        $ans = $x;
        else if ($x > $y) {
        $temp = $y / $x;
        $ans = $x * sqrt(1.0 + $temp * $temp);
        }
        else {
        $temp = $x / $y;
            $ans = $y * sqrt(1.0 + $temp * $temp);
        }
      
        return $ans;
}

/**
 * Return the radix of two complex number
 * @param complex $z
 * @return complex
 */
static public function Csqrt($z) {
if (($z->getReal() == 0.0) && ($z->getImag() == 0.0)) {
$c = new Complex(0, 0);       
return $c;
}
        else {
        $x = abs($z->getReal());
            $y = abs($z->getImag());
          
            if ($x >= $y) {
            $r = $y / $x;
                $w = sqrt($x) * sqrt(0.5 * (1.0 + sqrt(1.0 + $r * $r)));
            }
            else {
            $r = $x / $y;
                $w = sqrt($y) * sqrt(0.5 * ($r + sqrt(1.0 + $r * $r)));
            }
          
            $c = new Complex(0, 0);
          
            if ($z->getReal() >= 0.0) {
            $c->setReal($w);
            $c->setImag($z->getImag() / (2.0 * $w));
            }
            else {
            if ($z->getImag() >= 0)
            $c->setImag($w);
            else
            $c->setImag(-$w);
            $c->setReal($z->getReal() / (2.0 * $c->getImag()));
            }
          
            return $c;
        }
}

/**
 * Return the inverse of a complex number
 * @param complex $z
 * @return complex
 */
static public function Cinv($z) {
$c = new Complex(0, 0);

$c->setReal($z->getReal());
$c->setImag(-$z->getImag());

return $c;
}

/**
 * Return the multiplication of a complex number with a scalar value
 * @param complex $n
 * @param complex $z
 * @return complex
 */
static public function RCmul($n, $z) {
$c = new Complex(0, 0);

$c->setReal($z->getReal() * $n);
$c->setImag(-$z->getImag() * $n);

return $c;
}

/**
* Getters and setters
*/
public function setReal($real) {
$this->real = $real;
}

public function setImag($imag) {
$this->imag = $imag;
}

public function getReal() {
return $this->real;
}

public function getImag() {
return $this->imag;
}
}
?>



サンプルファイル
index.php
<?php
/**
 * @author Michele Andreoli <michi.andreoli@gmail.com>
 * @name index.php
 * @version 0.1 updated 07-05-2010
 * @license http://opensource.org/licenses/gpl-license-php GNU Public License
 * @package FFT
 */
?>

<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
<html>
    <head>
        <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
        <title>Fast Fourier Transform, FFT</title>
    </head>
    <body>
    <?php
    require_once 'FFT.class.php';
    
    // Define a generic function
    $f = array();
    for ($i = 0; $i < 256; $i++) {
     if (($i >= 0 && $i <= 8) || ($i >= 248 && $i <= 255))
    $f[$i] = 1;
     else
    $f[$i] = 0;
}
    
    $fft = new FFT(256);
    
    // Calculate the FFT of the function $f
    $w = $fft->fft($f);
    
    echo "<h1 style=\"font: bold 14px verdana;\">FFT: </h1>";
    for ($i = 0; $i < $fft->getDim(); $i++)
    echo "<p style=\"font: normal 10px verdana;\"><b>".$i."</b>  (".$w[$i]->getReal().", ".$w[$i]->getImag().")</p>";
    
    // Calculate the inverse FFT of the function $w
    $w = $fft->ifft($w);
    
    echo "<br/><h1 style=\"font: bold 14px verdana;\">FFT inverse:</h1>";
    for ($i = 0; $i < $fft->getDim(); $i++)
    echo "<p style=\"font: normal 10px verdana;\"><b>".$i."</b>  (".$w[$i]->getReal().")</p>";
    ?>
    </body>
</html>


サンプルテスト出力
FFT:
0 (17, 0)
1 (16.877376851216, 1.3239409568655E-14)
2 (16.512674319733, 1.8041124150159E-14)
3 (15.915296255567, 4.2299497238218E-14)
4 (15.100595973038, 2.4868995751604E-14)
5 (14.089407730913, 5.7842619582971E-14)
6 (12.9074126744, 6.1672889017927E-14)
7 (11.584360195483, 1.0079437284816E-13)
8 (10.153170387609, 2.8421709430404E-14)
9 (8.6489471020338, 6.9819150461115E-14)
10 (7.1079339256376, 6.8889338677991E-14)
11 (5.5664470995395, 1.1018963519405E-13)
12 (4.0598199297288, 6.3782312764715E-14)
13 (2.6213925918954, 9.886536034287E-14)
14 (1.281579431215, 9.2745255919624E-14)
15 (0.067042973107172, 1.2192330478555E-13)
16 (-1, 1.6168475976272E-14)
17 (-1.9023196470646, 4.4762804574106E-14)
18 (-2.6280462447661, 3.8344327712991E-14)
19 (-3.1708658792552, 5.6926685587655E-14)
20 (-3.5300352907046, 3.0420110874729E-14)
21 (-3.7102186573424, 4.0800696154975E-14)
22 (-3.7211550883243, 3.7359004778637E-14)
23 (-3.5771712024474, 4.0897840669629E-14)
24 (-3.2965582089384, 2.0539125955565E-14)
25 (-2.9008372059234, 2.0969337377608E-14)
26 (-2.4139398281315, 1.8651746813703E-14)
27 (-1.861333797279, 1.8623991238087E-14)
28 (-1.2691242823883, 1.6181500583912E-14)
29 (-0.66316222802768, 1.4557799410397E-14)
30 (-0.068189961048635, 1.4738210651899E-14)
31 (0.49294751478403, 1.4013096238941E-14)
32 (0.99999999999999, 1.5505964705829E-14)
33 (1.435896549321, 1.702284146976E-14)
34 (1.7871995226591, 1.7548462682981E-14)
35 (2.0444158500796, 2.1357915436226E-14)
36 (2.2021584578235, 1.8235413179468E-14)
37 (2.259156292821, 2.4008572907519E-14)
38 (2.2181168286947, 2.5313084961454E-14)
39 (2.0854501403741, 3.3292812950947E-14)
40 (1.8708684117894, 1.5598633495983E-14)
41 (1.5868789244465, 2.4744095661333E-14)
42 (1.2481920207445, 2.3758772726978E-14)
43 (0.87106813048327, 3.486100297323E-14)
44 (0.4726296108078, 2.120525977034E-14)
45 (0.070163833099847, 2.9559688030645E-14)
46 (-0.31955635390749, 2.7713942252205E-14)
47 (-0.68090991947076, 3.5672853559987E-14)
48 (-1, 4.4087703354695E-15)
49 (-1.2651415593768, 1.1664280652468E-14)
50 (-1.4672482570444, 9.1454621653497E-15)
51 (-1.6001065896392, 1.3933298959046E-14)
52 (-1.6605299057846, 6.883382752676E-15)
53 (-1.6483896157951, 8.8817841970013E-15)
54 (-1.5665256474979, 8.5209617139981E-15)
55 (-1.4205427576073, 8.2850393212652E-15)
56 (-1.218503525588, 6.4670491184415E-15)
57 (-0.97053257959026, 6.3976601794025E-15)
58 (-0.68834969871281, 6.8556271770603E-15)
59 (-0.38475179521194, 5.911937606129E-15)
60 (-0.073065326541002, 8.3544282603043E-15)
61 (0.23340862538904, 7.4593109467003E-15)
62 (0.52193455378641, 9.1801566348693E-15)
63 (0.78084698453526, 1.0143795525774E-14)
64 (0.99999999999999, 1.2170819907453E-14)
65 (1.1711451790866, 1.2678226524176E-14)
66 (1.288224448254, 1.4023504579797E-14)
67 (1.347567954397, 1.6556200854723E-14)
68 (1.3479910105092, 1.4530043834782E-14)
69 (1.2907882704857, 1.8804402479589E-14)
70 (1.1796273977999, 1.8651746813703E-14)
71 (1.0203484436007, 2.3717139363555E-14)
72 (0.82067879082865, 1.2434497875802E-14)
73 (0.5898767140905, 1.7194579093882E-14)
74 (0.33831923140421, 1.5903944827755E-14)
75 (0.077051888329222, 2.1038726316647E-14)
76 (-0.18268065064689, 1.2795320358805E-14)
77 (-0.42990386290684, 1.6014967130218E-14)
78 (-0.65436345108323, 1.4058199049316E-14)
79 (-0.84693828436805, 1.6549261960819E-14)
80 (-1, 1.6495537311177E-15)
81 (-1.1077054569259, 3.4486302702419E-15)
82 (-1.1662111434634, 1.4571677198205E-15)
83 (-1.1738019485656, 1.7486012637846E-15)
84 (-1.1309302754798, -6.6613381477509E-16)
85 (-1.0401651618179, -2.2759572004816E-15)
86 (-0.90605472289853, -2.5535129566379E-15)
87 (-0.73490871010854, -4.8433479449272E-15)
88 (-0.5345111359508, -2.2204460492503E-15)
89 (-0.31377563932833, -4.787836793696E-15)
90 (-0.08235844401145, -3.9968028886506E-15)
91 (0.14975468163902, -5.3568260938164E-15)
92 (0.37267018089141, -2.3037127760972E-15)
93 (0.57699973210147, -2.7478019859473E-15)
94 (0.75424763347272, -1.04777297949E-15)
95 (0.89715686611487, -3.9898639947467E-17)
96 (0.99999999999999, 2.1832884103274E-15)
97 (1.0588034996283, 4.0002723356025E-15)
98 (1.0714968217812, 5.4053983511437E-15)
99 (1.0379808584093, 8.7291285311153E-15)
100 (0.96011363307539, 7.6882944455292E-15)
101 (0.84161457218897, 1.1601830607333E-14)
102 (0.68789200394856, 1.2434497875802E-14)
103 (0.50580165410466, 1.6084356069257E-14)
104 (0.30334668360734, 8.1601392309949E-15)
105 (0.089332138671147, 1.1143863609675E-14)
106 (-0.12701153259932, 1.0380585280245E-14)
107 (-0.33642007008792, 1.3322676295502E-14)
108 (-0.52998575303648, 7.8270723236074E-15)
109 (-0.69953032173529, 8.3544282603043E-15)
110 (-0.83794509408247, 6.5919492087119E-15)
111 (-0.93948434240108, 5.8772431366094E-15)
112 (-1, 8.6583076418495E-17)
113 (-1.0171082499518, -1.3392065234541E-15)
114 (-0.9902814070487, -3.2335245592208E-15)
115 (-0.92086161973402, -5.6343818499727E-15)
116 (-0.81199616386893, -5.384581669432E-15)
117 (-0.66849734162249, -8.9928064994638E-15)
118 (-0.49663310430789, -9.4368957093138E-15)
119 (-0.30385736133689, -1.2559397966072E-14)
120 (-0.098491403357184, -6.7723604502135E-15)
121 (0.11063014353736, -1.0061396160665E-14)
122 (0.31453398385415, -9.2703622556201E-15)
123 (0.50448777275728, -1.1324274851177E-14)
124 (0.6723688525762, -7.2164496600635E-15)
125 (0.8110067489246, -7.9936057773011E-15)
126 (0.91448518154279, -5.8841820305133E-15)
127 (0.97839104176565, -4.6074255521944E-15)
128 (0.99999999999997, -2.3037127760972E-15)
129 (0.97839104176564, -1.6375789613221E-15)
130 (0.91448518154278, -1.2212453270877E-15)
131 (0.81100674892459, 8.8817841970013E-16)
132 (0.6723688525762, -4.4408920985006E-16)
133 (0.50448777275728, 2.1094237467878E-15)
134 (0.31453398385414, 2.4980018054066E-15)
135 (0.11063014353735, 5.0931481254679E-15)
136 (-0.098491403357185, -5.5511151231257E-16)
137 (-0.30385736133688, 1.2073675392799E-15)
138 (-0.49663310430789, 7.2164496600635E-16)
139 (-0.6684973416225, 1.720845688169E-15)
140 (-0.81199616386894, -6.1062266354384E-16)
141 (-0.92086161973403, -1.498801083244E-15)
142 (-0.9902814070487, -2.3453461395206E-15)
143 (-1.0171082499518, -4.3645642655576E-15)
144 (-1, -1.2384394644671E-15)
145 (-0.93948434240107, -4.0314973581701E-15)
146 (-0.83794509408246, -4.4547698863084E-15)
147 (-0.69953032173528, -8.076872504148E-15)
148 (-0.52998575303648, -4.6629367034257E-15)
149 (-0.33642007008791, -8.9372953482325E-15)
150 (-0.12701153259932, -8.3821838359199E-15)
151 (0.089332138671149, -1.2281842209916E-14)
152 (0.30334668360734, -3.441691376338E-15)
153 (0.50580165410467, -7.4523720527964E-15)
154 (0.68789200394857, -6.7723604502135E-15)
155 (0.84161457218897, -9.3536289824669E-15)
156 (0.9601136330754, -5.1347814888913E-15)
157 (1.0379808584093, -6.3143934525556E-15)
158 (1.0714968217812, -5.0237591864288E-15)
159 (1.0588034996283, -4.6108949991464E-15)
160 (0.99999999999999, -3.335144798376E-15)
161 (0.89715686611487, -2.6281060661049E-15)
162 (0.75424763347272, -1.9914625504214E-15)
163 (0.57699973210146, -1.3877787807814E-17)
164 (0.37267018089141, -1.3600232051658E-15)
165 (0.14975468163902, 8.0491169285324E-16)
166 (-0.082358444011456, 9.9920072216264E-16)
167 (-0.31377563932834, 1.873501354055E-15)
168 (-0.53451113595081, -1.0547118733939E-15)
169 (-0.73490871010853, -1.346145417358E-15)
170 (-0.90605472289853, -1.942890293094E-15)
171 (-1.0401651618179, -3.5527136788005E-15)
172 (-1.1309302754799, -3.663735981263E-15)
173 (-1.1738019485656, -7.3552275381417E-15)
174 (-1.1662111434634, -8.7568841067309E-15)
175 (-1.1077054569259, -1.4398204850607E-14)
176 (-1, -2.8014101191664E-15)
177 (-0.84693828436803, -9.3050567251396E-15)
178 (-0.65436345108322, -1.0227929614359E-14)
179 (-0.42990386290684, -1.7263968032921E-14)
180 (-0.18268065064689, -1.1157741397483E-14)
181 (0.077051888329233, -1.8374191057546E-14)
182 (0.33831923140422, -1.8485213360009E-14)
183 (0.58987671409051, -2.4216739724636E-14)
184 (0.82067879082866, -1.2184697695261E-14)
185 (1.0203484436007, -1.8027246362351E-14)
186 (1.1796273977999, -1.7347234759768E-14)
187 (1.2907882704857, -2.0969337377608E-14)
188 (1.3479910105092, -1.5834555888716E-14)
189 (1.3475679543971, -1.7638668303732E-14)
190 (1.288224448254, -1.6261297863807E-14)
191 (1.1711451790866, -1.5792922525293E-14)
192 (0.99999999999999, -1.4474532683551E-14)
193 (0.78084698453525, -1.3745081461902E-14)
194 (0.52193455378641, -1.4269835313385E-14)
195 (0.23340862538904, -1.2823075934421E-14)
196 (-0.073065326541002, -1.3697376566313E-14)
197 (-0.38475179521194, -1.3086753902769E-14)
198 (-0.68834969871281, -1.301736496373E-14)
199 (-0.97053257959026, -1.40304434737E-14)
200 (-1.218503525588, -1.3600232051658E-14)
201 (-1.4205427576073, -1.5751289161869E-14)
202 (-1.5665256474979, -1.6292522886374E-14)
203 (-1.6483896157951, -2.18436380095E-14)
204 (-1.6605299057846, -1.6930901125534E-14)
205 (-1.6001065896392, -2.4397150966138E-14)
206 (-1.4672482570444, -2.4771851236949E-14)
207 (-1.2651415593768, -3.4590386110978E-14)
208 (-1, -5.5606267235181E-15)
209 (-0.68090991947075, -1.5938639297275E-14)
210 (-0.31955635390748, -1.4557799410397E-14)
211 (0.070163833099856, -2.4896751327219E-14)
212 (0.4726296108078, -1.2267964422108E-14)
213 (0.87106813048328, -2.1982415887578E-14)
214 (1.2481920207446, -2.1094237467878E-14)
215 (1.5868789244465, -2.9046209881756E-14)
216 (1.8708684117894, -8.4376949871512E-15)
217 (2.0854501403741, -1.5445977830097E-14)
218 (2.2181168286947, -1.4210854715202E-14)
219 (2.259156292821, -1.9345636204093E-14)
220 (2.2021584578235, -1.3239409568655E-14)
221 (2.0444158500796, -1.5737411374062E-14)
222 (1.7871995226591, -1.5480672299617E-14)
223 (1.435896549321, -1.5943843467703E-14)
224 (0.99999999999999, -1.6657821093878E-14)
225 (0.49294751478403, -1.8537255064288E-14)
226 (-0.068189961048636, -1.9852175459079E-14)
227 (-0.66316222802769, -2.152444888992E-14)
228 (-1.2691242823883, -2.36199948489E-14)
229 (-1.861333797279, -2.647881913731E-14)
230 (-2.4139398281315, -3.0087043967342E-14)
231 (-2.9008372059234, -3.9315772859538E-14)
232 (-3.2965582089384, -2.8366198279173E-14)
233 (-3.5771712024474, -4.203581926987E-14)
234 (-3.7211550883243, -4.4242387531312E-14)
235 (-3.7102186573424, -6.6058269965197E-14)
236 (-3.5300352907046, -4.9349413444588E-14)
237 (-3.1708658792552, -7.4690253981657E-14)
238 (-2.6280462447661, -7.7562956057875E-14)
239 (-1.9023196470647, -1.1105699693204E-13)
240 (-1, -1.7320332364321E-14)
241 (0.0670429731072, -5.1077198026661E-14)
242 (1.281579431215, -4.6809778275758E-14)
243 (2.6213925918954, -8.2572837456496E-14)
244 (4.0598199297288, -3.9468428525424E-14)
245 (5.5664470995396, -7.2164496600635E-14)
246 (7.1079339256377, -7.0499162063697E-14)
247 (8.6489471020339, -9.8157593164672E-14)
248 (10.153170387609, -2.0428103653103E-14)
249 (11.584360195483, -4.3479109201883E-14)
250 (12.9074126744, -3.4694469519536E-14)
251 (14.089407730913, -5.2402526762307E-14)
252 (15.100595973038, -1.987299214079E-14)
253 (15.915296255568, -2.076117056049E-14)
254 (16.512674319733, -1.0935696792558E-14)
255 (16.877376851216, -4.4408920985006E-16)

FFT inverse: 逆変換
0 (1)
1 (1)
2 (1)
3 (1)
4 (1)
5 (1)
6 (1)
7 (1)
8 (1)
9 (-1.492017013091E-15)
10 (-1.3672015133964E-15)
11 (-1.6996440096182E-15)
12 (-1.2863480513928E-15)
13 (-1.5587683962791E-15)
14 (-1.5016291471669E-15)
15 (-1.6816325715574E-15)
16 (-7.4297729778648E-16)
17 (-1.0539044870699E-15)
18 (-1.0639905291507E-15)
19 (-1.3505522638021E-15)
20 (-1.0055282904977E-15)
21 (-1.2705716369273E-15)
22 (-1.2692721198654E-15)
23 (-1.4216144240241E-15)
24 (-2.1057182970151E-15)
25 (-1.867667083695E-15)
26 (-1.8321892485993E-15)
27 (-2.0814935382908E-15)
28 (-1.8419014220781E-15)
29 (-2.0550607086268E-15)
30 (-2.0158296086263E-15)
31 (-2.1691698800222E-15)
32 (-9.365943393149E-16)
33 (-1.1401018292578E-15)
34 (-1.1826072002246E-15)
35 (-1.3013645472888E-15)
36 (-1.2983146940507E-15)
37 (-1.4208495467536E-15)
38 (-1.4754605295649E-15)
39 (-1.6164844402858E-15)
40 (-1.5195956075257E-15)
41 (-9.9636106447977E-16)
42 (-1.0732633594539E-15)
43 (-1.2789339263661E-15)
44 (-1.0791084846866E-15)
45 (-1.2951014272876E-15)
46 (-1.2947108099685E-15)
47 (-1.4547926136329E-15)
48 (-7.3163494500363E-16)
49 (-8.7024357800018E-16)
50 (-9.5419792448397E-16)
51 (-1.1567488817289E-15)
52 (-9.1037531451755E-16)
53 (-1.1268657940098E-15)
54 (-1.1827819191168E-15)
55 (-1.3350199557727E-15)
56 (-3.188425109829E-15)
57 (-3.0245646842624E-15)
58 (-3.0252419075903E-15)
59 (-3.3207548064253E-15)
60 (-3.0911596029309E-15)
61 (-3.3091555033119E-15)
62 (-3.1818270138031E-15)
63 (-3.3695938743309E-15)
64 (-1.424191033118E-15)
65 (-1.5637786891793E-15)
66 (-1.5534733103025E-15)
67 (-1.7288728183304E-15)
68 (-1.7184179952312E-15)
69 (-1.7403457982855E-15)
70 (-1.8861274350464E-15)
71 (-1.9371085554654E-15)
72 (-2.0498557684724E-15)
73 (-1.012135590701E-15)
74 (-1.0864084616883E-15)
75 (-1.2528758928463E-15)
76 (-1.1556351459127E-15)
77 (-1.2300250974298E-15)
78 (-1.3418856072278E-15)
79 (-1.4616842235028E-15)
80 (-8.0012558618783E-16)
81 (-8.4112119101412E-16)
82 (-8.4782378723931E-16)
83 (-1.074041317721E-15)
84 (-8.7248324027388E-16)
85 (-1.0535146860055E-15)
86 (-1.0719161805933E-15)
87 (-1.2953142869287E-15)
88 (-1.8314926186071E-15)
89 (-1.6327361994757E-15)
90 (-1.6185793486579E-15)
91 (-1.802471513715E-15)
92 (-1.637156567619E-15)
93 (-1.8436696879157E-15)
94 (-1.775199972167E-15)
95 (-2.0095845587555E-15)
96 (-9.3979483213529E-16)
97 (-1.0479376107947E-15)
98 (-1.1253203833854E-15)
99 (-1.2100307147055E-15)
100 (-1.1236165240294E-15)
101 (-1.2317064737292E-15)
102 (-1.2767354929274E-15)
103 (-1.4702649386853E-15)
104 (-1.4987301307366E-15)
105 (-1.0394094576463E-15)
106 (-1.0783908549738E-15)
107 (-1.2926283889006E-15)
108 (-1.1062930635837E-15)
109 (-1.2623368725199E-15)
110 (-1.3034525271972E-15)
111 (-1.4596606618053E-15)
112 (-7.4400884754136E-16)
113 (-8.495059973224E-16)
114 (-8.7000388844759E-16)
115 (-1.0325677670822E-15)
116 (-8.6864591804354E-16)
117 (-1.0838828941156E-15)
118 (-1.0701171512984E-15)
119 (-1.2450859934057E-15)
120 (-5.6621374255883E-15)
121 (-5.3290705182008E-15)
122 (-5.3290705182008E-15)
123 (-5.4400928206633E-15)
124 (-5.3290705182008E-15)
125 (-5.4400928206633E-15)
126 (-5.4400928206633E-15)
127 (-5.6621374255883E-15)
128 (-2.1094237467878E-15)
129 (-2.2204460492503E-15)
130 (-2.3314683517128E-15)
131 (-2.5535129566379E-15)
132 (-2.2204460492503E-15)
133 (-2.5535129566379E-15)
134 (-2.6645352591004E-15)
135 (-2.6645352591004E-15)
136 (-2.7755575615629E-15)
137 (-1.0664723756337E-15)
138 (-1.0653632260952E-15)
139 (-1.3034231938083E-15)
140 (-1.0922414885677E-15)
141 (-1.2422479088185E-15)
142 (-1.2967615215007E-15)
143 (-1.4375446473715E-15)
144 (-7.0861208260941E-16)
145 (-8.1334397417021E-16)
146 (-8.5575897398995E-16)
147 (-1.0240971889493E-15)
148 (-8.4921747447718E-16)
149 (-1.0175128832528E-15)
150 (-1.0494516646105E-15)
151 (-1.1687972770004E-15)
152 (-1.8376590204648E-15)
153 (-1.596199898963E-15)
154 (-1.6552289120617E-15)
155 (-1.7994157191065E-15)
156 (-1.7257188016859E-15)
157 (-1.8789415438276E-15)
158 (-1.8620591401959E-15)
159 (-2.014658766747E-15)
160 (-9.394548898454E-16)
161 (-1.0099298826948E-15)
162 (-1.0476023046346E-15)
163 (-1.1980375100262E-15)
164 (-1.0470001506969E-15)
165 (-1.1466985809999E-15)
166 (-1.1606235980402E-15)
167 (-1.3938997130671E-15)
168 (-1.4312855457147E-15)
169 (-9.9750083518017E-16)
170 (-1.0086736196821E-15)
171 (-1.254721128795E-15)
172 (-1.0474870480059E-15)
173 (-1.1876502118413E-15)
174 (-1.233255510814E-15)
175 (-1.3821432938111E-15)
176 (-7.8240867335056E-16)
177 (-7.700602400744E-16)
178 (-8.3253046443266E-16)
179 (-9.4779184981862E-16)
180 (-8.6569033938408E-16)
181 (-9.621950387616E-16)
182 (-1.0580314097878E-15)
183 (-1.1308636059932E-15)
184 (-3.0288238280719E-15)
185 (-2.7485950437885E-15)
186 (-2.858940122923E-15)
187 (-3.0075164339381E-15)
188 (-2.9040447300449E-15)
189 (-3.0191157370515E-15)
190 (-3.0354219240978E-15)
191 (-3.1807219709575E-15)
192 (-1.3513665284449E-15)
193 (-1.4338234773087E-15)
194 (-1.4441288561854E-15)
195 (-1.6017962555451E-15)
196 (-1.5012287761817E-15)
197 (-1.59032327559E-15)
198 (-1.6665862437541E-15)
199 (-1.8376497282602E-15)
200 (-1.7249025152532E-15)
201 (-9.6713422438931E-16)
202 (-9.8946529908879E-16)
203 (-1.2180104829718E-15)
204 (-1.0537751409324E-15)
205 (-1.1976128930248E-15)
206 (-1.2517525584894E-15)
207 (-1.3676976630405E-15)
208 (-7.2268527341301E-16)
209 (-7.8395576313589E-16)
210 (-8.2698016624175E-16)
211 (-9.3543817621353E-16)
212 (-8.1746670441286E-16)
213 (-9.5230243502066E-16)
214 (-1.0080419807132E-15)
215 (-1.1560575711521E-15)
216 (-1.774646631364E-15)
217 (-1.4537126631547E-15)
218 (-1.555340638432E-15)
219 (-1.6440911914137E-15)
220 (-1.5675836588305E-15)
221 (-1.6608223246185E-15)
222 (-1.7854055439993E-15)
223 (-1.9112148742389E-15)
224 (-9.0340307119869E-16)
225 (-9.0985586836578E-16)
226 (-9.7433990779345E-16)
227 (-1.0645262338677E-15)
228 (-9.7196072972361E-16)
229 (-9.7470440440547E-16)
230 (-1.0831839902808E-15)
231 (-1.1814883335501E-15)
232 (-1.3235484440738E-15)
233 (-9.2217557726121E-16)
234 (-9.3546210646645E-16)
235 (-1.13585940817E-15)
236 (-9.4987347145686E-16)
237 (-1.0182644144254E-15)
238 (-1.1016264466494E-15)
239 (-1.3011637813802E-15)
240 (-7.6275162708356E-16)
241 (-7.3471406819511E-16)
242 (-6.8760816992127E-16)
243 (-7.4992408814179E-16)
244 (-5.8295316860667E-16)
245 (-4.7124925797664E-16)
246 (-5.6154910747191E-16)
247 (-5.7312029257445E-16)
248 (1)
249 (1)
250 (1)
251 (1)
252 (1)
253 (1)
254 (1)
255 (1)

元に戻っている。



以上です。当方PHP5、Winにて動作確認。
エクセルやR言語などでの出力結果と比較確認してから使うべし。
類似:
http://mitsui725.hatenablog.com/entry/2013/01/02/045859
タグ:
関数、クラス、ライブラリ、離散、FFT、
PHP class libray  lib fourier analytics series transform


/**
 * @author Michele Andreoli <michi.andreoliアットマークgmail.com>
 * @name FFT.class.php
 * @version 0.5 updated 06-07-2010
 * @license http://opensource.org/licenses/gpl-license-php GNU Public License
 * @package FFT
 */


/**
 * @author Michele Andreoli <michi.andreoliアットマークgmail.com>
 * @name Complex.class.php
 * @version 0.3 updated 09-05-2010
 * @license http://opensource.org/licenses/gpl-license-php GNU Public License
 * @package FFT
 */