Rem 数値的逆ラプラス変換filt06(t) Hosonoのアルゴリズム----------------------
Rem atstp:分母 btstp:分子---ラプラス関数値の算出に必要なパラメータをセル参照
Rem 作例 F(s) = ht(8次多項式)/gt(8次多項式) : 有理関数
Public Function filt06(time As Double, atstp As Variant, btstp As Variant, nts As Integer) As Double
Dim ss(500) As Complex
Dim x1L(500) As Complex
Dim res(500) As Double
Dim aa As Double
Dim n As Integer
Dim cnt As Double
Rem パラメータ関連・複素数変数
Dim resigma As Double
Dim ared(50, 2) As Double
Dim bred(50, 2) As Double
Dim pi As Double
Dim a01(50) As Complex
Dim b01(50) As Complex
Dim ht(50) As Complex
Dim gt(50) As Complex
Dim ht01 As Complex
Dim ht02 As Complex
Dim ht03 As Complex
Dim ht04 As Complex
Dim gt01 As Complex
Dim gt02 As Complex
Dim gt03 As Complex
Dim gt04 As Complex
Rem 各種定数
pi = 3.14159265358979
Rem パラメータの取り込み・複素数変換
Rem 指定範囲から順番にデータ用セル値を取り出す（次数の高い順a0-a8,b0-b8）
For Each myCell In atstp
ared(cnt, 0) = myCell.Value
a01(cnt) = ToComplex(ared(cnt, 0), 0)
cnt = cnt + 1
Next
cnt = 0
For Each myCell In btstp
bred(cnt, 0) = myCell.Value
b01(cnt) = ToComplex(bred(cnt, 0), 0)
cnt = cnt + 1
Next
Rem 高速ラプラス変換（細野法）
If time > 0 Then
aa = 8
For n = 1 To nts + 10
ss(n) = ToComplex(aa / time, (n - 0.5) / time * pi)
Rem ラプラス変換関数 Cadd Csub Cmul Cdiv ToComplex F(ss(n))算出
gt(7) = ss(n)
gt(6) = Cmul(gt(7), ss(n))
gt(5) = Cmul(gt(6), ss(n))
gt(4) = Cmul(gt(5), ss(n))
gt(3) = Cmul(gt(4), ss(n))
gt(2) = Cmul(gt(3), ss(n))
gt(1) = Cmul(gt(2), ss(n))
gt(0) = Cmul(gt(1), ss(n))
gt01 = a01(8)
gt01 = Cadd(gt01, Cmul(gt(7), a01(7)))
gt01 = Cadd(gt01, Cmul(gt(6), a01(6)))
gt01 = Cadd(gt01, Cmul(gt(5), a01(5)))
gt01 = Cadd(gt01, Cmul(gt(4), a01(4)))
gt01 = Cadd(gt01, Cmul(gt(3), a01(3)))
gt01 = Cadd(gt01, Cmul(gt(2), a01(2)))
gt01 = Cadd(gt01, Cmul(gt(1), a01(1)))
gt01 = Cadd(gt01, Cmul(gt(0), a01(0)))
ht(7) = ss(n)
ht(6) = Cmul(ht(7), ss(n))
ht(5) = Cmul(ht(6), ss(n))
ht(4) = Cmul(ht(5), ss(n))
ht(3) = Cmul(ht(4), ss(n))
ht(2) = Cmul(ht(3), ss(n))
ht(1) = Cmul(ht(2), ss(n))
ht(0) = Cmul(ht(1), ss(n))
ht01 = b01(8)
ht01 = Cadd(ht01, Cmul(ht(7), b01(7)))
ht01 = Cadd(ht01, Cmul(ht(6), b01(6)))
ht01 = Cadd(ht01, Cmul(ht(5), b01(5)))
ht01 = Cadd(ht01, Cmul(ht(4), b01(4)))
ht01 = Cadd(ht01, Cmul(ht(3), b01(3)))
ht01 = Cadd(ht01, Cmul(ht(2), b01(2)))
ht01 = Cadd(ht01, Cmul(ht(1), b01(1)))
ht01 = Cadd(ht01, Cmul(ht(0), b01(0)))
x1L(n) = Cdiv(ht01, gt01)
Rem ラプラス変換関数 Cadd Csub Cmul Cdiv ToComplex
res(n) = Imz(x1L(n)) * (-1) ^ n
Next
Rem nts項分加算
resigma = 0
For n = 1 To nts
resigma = resigma + res(n)
Next
Rem オイラー変換-10項分追加
resigma = resigma + (res(nts + 1) * 1023 + res(nts + 2) * 1013 + res(nts + 3) * 968 + res(nts + 4) * 848 + res(nts + 5) * 638 ) / 2 ^ 10
resigma = resigma + (res(nts + 6) * 386 + res(nts + 7) * 176 + res(nts + 8) * 56 + res(nts + 9) * 11 + res(nts + 10)) / 2 ^ 10
Rem 変換値出力
filt06 = resigma * Exp(aa) / time
Else
Rem t=0の関数値：初期値の出力
filt06 = 0
End If
End Function
