周波数応答：pade関数の作成

MATLABにはpade関数が用意されています。
Scilabにも同じような関数がないか探しました。
しかし、残念ながら見つけることが出来ませんでした（しっかり探せばあるのかも？）。
そこで、独自にpade関数を作成しました。

仕様

＜概要＞
むだ時間と次数を渡すと、パデ近似手法で得られる有理式の分子と分母の各次数の係数を返す。

コーリングシーケンス	パラメータ
[num,den]=pade(t,n)	入力	t	むだ時間[Sec]
		n	次数
	戻り値	num	分子の各次数の係数
		den	分母の各次数の係数

＜実行例：遅れ時間=0.1Sec、次数=4＞

コンソール画面

-->getf('pade.sci')        //←pade.sciファイルに定義(pade)した関数の使用宣言
-->[num,den]=pade(0.1,4)   //←pade関数の実行
 den  =
 
    16800000.  
    840000.    
    18000.     
    200.       
    1.         
 num  =
 
    16800000.  
  - 840000.    
    18000.     
  - 200.       
    1.

上記結果を伝達関数に書き換えると次式になります。

（注）
MATLABでは次数の高い方から返します。
Scilabでは伝達関数の表示等、次数の低い方から表示します。
よって、今回作成したpade関数もScilabに合わせて、次数の低い方から返すようにしました。

プログラム

前節の最後で導き出した、下記の連立方程式を解くプログラムを作成します。

筆者は、Scilabのプログラミングに関しては初心者です。
非効率・冗長なコーディングが多々あると思いますがご勘弁を（って言い訳をしておく）。

pade.sci

// D:むだ時間(Sec)
// N:パデ近似の次数
//
// p:分子の係数
// q:分母の係数
function [p, q]=pade(D, N)
  getf('fact.sci'); //階乗計算用関数
  num = N*2;

  //e^(-Dx) のラプラス展開
  k=1;
  for i=1:num+1 do
    a(i) = k * 1/fact(i-1) * D^(i-1);
    k=k * -1;
  end
  
  //連立方程式の作成
  tbl = zeros(num, num+1);
  for k=1: num do
    for i=1: k do
      if(k+1-i <=N) then
        tbl(k, k+1-i) = a(i);
      end
    end
    tbl(k, num+1) = a(i+1);
  end
  for i=1:N do
    tbl(i, N+i)=-1;
  end

  //掃き出し法による連立方程式の解の計算
  for k = 1: num do
    p = tbl(k,k);
    for i = k: num+1 do
      tbl(k, i) = tbl(k, i) / p;
    end

    for i = 1: num do
      if(i ~= k) then
        w = tbl(i, k); 
        for j = k: num+1 do
          tbl(i, j) = tbl(i, j) - w * tbl(k, j);
        end
      end
    end
  end

  //分子・分母の計算
  dq = abs(1 / tbl(N, num+1));
  dp = abs(1 / tbl(N*2, num+1));
  q(1) = 0 - 1 / tbl(N, num+1);
  p(1) = 0 - 1 / tbl(num, num+1);
  for i=1: N do
    q(i+1) = 0 - tbl(i, num+1) * dq;
    p(i+1) = 0 - tbl(N+i, num+1) * dp;
  end
endfunction

＜階乗を求める関数＞

fact.sci

function f=fact(n)
  if(n <= 1) then
    f=1;
  else
    f = n * fact( n-1 );
  end
endfunction

筆者は、ここに提示したプログラムの著作権を主張しません。
自由にコピペしてお使いください。
ただし、その事によって生じた損害・不利益等の責任を、筆者は一切持ちません