FIRと血中濃度

この項はオーディオとは全く関係ありません。オーディオ関係の方はお戻りください。

目次へ

インパルスレスポンス(IR)

オーディオの測定でインパルスレスポンス(以下IR)をいうのを測定することがあります。インパルスとは極めて短時間の入力で非常に広い周波数帯域を持ちます。IRはインパルスの入力があった時に応答(出力)がどうなるかということです。 IRを測定すると何がわかるのかというと全て(非線形部分を除く)が分かります。IRすなわち伝達関数です。周波数特性、位相特性、任意の入力があった時に出力がどのように変化するかなどが分かります。

IRの利用例としてオーディオのエフェクトでエコー・リバーブを考えてみます。教会やコンサートホール、風呂場、普通の部屋でも残響があります。部屋の中で音が発生した時に壁で反射していろいろなルートで耳まで音が到達します。その距離の違いによって遅れがきまり、反射する壁の材質で反射音の音質が変わります。残響音は反射を繰り返しながら小さくなりますが、計算上はいつまでたってもゼロにはなりません。すなわちIIR(無限インパルス応答)となります。

さてデジタルエコーで実際の部屋の残響を再現する方法ですが、大きくわけて2つあります。

1.モデリング

部屋の大きさ、音源の位置、マイクの位置から反射の計算式を求めてエコーをかけるもの。計算式によって反射を繰り返し小さくなりながら長時間続く残響を作り出すことができる。

2.実測値をFIRとして使う

実際の教会やホールでIRを測定してその結果をFIR(有限インパルス応答)のフィルタ係数として使う。

2.の場合には限りなく小さくなりながら無限に続く残響をどこかで切り取らなければなりません。残響がほとんどゼロとみなせるくらいになった時点で切り取れば問題ありません。どんなに残響が豊かなコンサートホールでも前日の演奏の残響がまだ残っているなんてことはありませんから。例えば2秒の残響をサンプリング周波数96kHz、32bitで測定すると192000サンプル、750kbyteのデータになります。

1.のモデリングによるIIR(無限インパルス応答)と2.のFIR(有限インパルス応答)を主観で比較します。

1.モデリング
モデリングに数学的知識が必要
モデリングに計算が必要だが実際にフィルタをかけるときの計算量は少ない。
エレガント。

2.実測値をFIRとして使う
測定データを使うだけなので数学的知識不要
モデリング不要だが実際にフィルタをかけるときの計算量は多い。
野暮。力ずく。

血中濃度

さて薬物を投与した時の血中濃度を考えます。IRに相当するものは「ボーラス投与した時の血中濃度の推移」です。もしその薬物の血中濃度に線形性があるのならばIR(ボーラス投与した時の血中濃度の推移)さえわかれば任意の投与に対する血中濃度が計算できます。

ここから先は投与と血中濃度に線形性がある場合だけを考えます。
コンパートメントモデルはどうしたですって?

「そんなの関係ねぇ! そんなの関係ねぇ!」
先に上げたエコーの作り方と同じことです。

それでいいのか?

いいんです! 

その理由を述べます。
そもそもなぜモデル化しなければならないかというと、そうすると計算が「楽」だからです。
薬物の血中濃度は後半では半減期で表せるように指数的に減っていきますがそう簡単にはゼロになりません。すなわちIIR(無限インパルス応答)となります。これをFIR(有限インパルス応答)で正確に計算しようとすると無限の長さのサンプルが必要になります。それで普通はフィードバック/フィードフォワードをかけた(出力結果が入力に戻る)IIRとして表します。そのためにコンパートメントモデルが使われます。ちなみに工学系の方はコンパートメントモデルではなくLCR(コイル・コンデンサ・抵抗)による等価回路にしてシミュレーターで計算するのが普通みたいです。

ここで有限のサンプル(FIR)でIIRの代用をすることを考えます。もしFIRで良いのならモデル化なんかしないでボーラス投与時の血中濃度の推移から任意の入力(投与)にたいする血中濃度の予測ができます。実際に測定した血中濃度を適当に曲線あてはめして後半は半減期で減少するようにすればOKです。実際の測定データのほうがモデリングした結果より重要です。モデリングと実測値が乖離した時に実測値のほうが間違っているとするか、モデリングの方を修正するかは言うまでもありません。

無限に続くIIRを有限のFIRで代用するには条件があります。

1.ある程度の有限の時間内に結果がゼロとみなせるくらいに小さくなる。

または

2.実際に予測したい時間よりもFIRのサンプルのほうが長い。

1.の場合、レミフェンタニルなど半減期が十分短い薬剤は例えばどんな量を投与しても数時間後にはゼロで良いでしょう。そうしても誤差は十分に小さいはずです。

2.の場合、例えば初回投与から血中濃度の予測が必要なのが20時間だと考えれば(大抵の手術はもっと短いはず)20時間後までのデータがあればその時点でゼロでなくてもかまわないことになります。もっとも20時間後にゼロでないところでいきなりぶった切ると20時間以降の予測値はガタガタになります。ちなみにオーディオ関係でサンプルを途中でぶった切った時には窓関数をかけて端を徐々にゼロにします。

実際に計算してみよう!

まずどれくらいのサンプルが必要か考えてみます。
サンプリング間隔ですが10秒くらいあれば良いでしょう。一秒刻みが必要になるとは考えられません。長さは20時間あれば大抵の手術には間に合います。そうすると7200ポイント必要になります。ここは切り良く8192としましょう。

計算量が多くなり過ぎないか? オーディオの世界では普通のパソコンでサンプリング周波数96kで6chにtap数8192のフィルタをかけても余裕です。それに比べたら計算量は10秒に一回なので576万分の1にしかなりません。

次に実際のIRのデータが必要です。ここではシミュレーターを使ってみます。実測値ではないですがご勘弁を。ついでに血中濃度ではなく仮想の効果部位濃度で計算してみます。投与量との関係に線形性があれば血中濃度でも効果部位濃度でも同じことです。

適当に計算してGNU Octaveに入力します。
試しにIRのデーターを作ってみます。
広島大学麻酔科のExcel_PkPdを使わせてもらいました。
フェンタニルを1mcg/kgボーラス投与した時の結果をoctaveで変数f1に取り込みます。
f1をプロットしてみます。横軸の単位は時間です。

k1

10秒間隔で8192サンプルなので23時間弱のデータです。120秒間隔で作ったデータをoctaveのinterp1関数で直線補間して10秒間隔としたのであまり良いデータではありませんが、実用上十分です。
これがIR、すなわちFIRフィルタの係数となります。

次に入力データを作ります。

a=[2 zeros(1,6*30-1), 1, zeros(1,6*30-1), 1, zeros(1,6*120)];

最初に2mcg/kg、30分と1時間後に1mcg/kgを投与したデータになります。
これをconv関数で畳み込みます。

c=conv(a,f1);

一瞬で計算が終わります。結果をプロットしでみます。横軸は時間です。

k2

もっともらしいグラフが描けました。
このグラフの終わりの方を示します。

k3

firフィルタが23時間弱しかないので初回投与から22時間後以降で値が不正確になってます。

次にレミフェンタニルを計算してみます。
同じようにExcel_PkPdをつかって1mcg/kg投与時の値を計算してoctaveで変数r1に取り込みます。すぐに値が小さくなるので10秒間隔で2.5時間くらいのデータにしました。横軸の単位は時間です。

k4

今度は持続投与をしてみます。1mcg/kg/minの投与量はサンプリングが10秒なので1/6となります。0.5γ15分、0.25γ75分、0.5γ30分、0.25γ30分、0.125γ30分、0.25γ30分のデータを作ります。

s1=1/6; #1γ 1mcg/kg/60sec
a=[ones(1,6*15)*s1*0.5,ones(1,6*75)*s1*0.25,ones(1,6*30)*s1*0.5,ones(1,6*30)*s1*0.25,ones(1,6*30)*s1*0.125,ones(1,6*30)*s1*0.25,zeros(1,6*120)];

aをプロットします。横軸の単位は時間、縦軸の単位はγです。

k5

次にconv関数で畳み込みます。

c=conv(r1,a);

結果をプロットします。

k6

今度もそれらしいグラフが描けました。

この方法でIR(ボーラス投与の血中濃度や効果部位濃度の推移)がわかればどんな薬剤でも計算できます(線形性があれば)。サンプリング時間を10分おきにすれば8192サンプルで56日間のデータになるので内服薬の血中濃度の予測なんかにも使えます。

ちなみにIRの逆関数を計算して任意の出力を得られるような入力を計算することができるか?考えてみると、計算は簡単にできます。ただし、残念ながらマイナスの値が必要になります。オーディオ信号なら交流なのでマイナスの値でも構いませんが、薬剤投与にマイナスの値は不可能です。結局無理な出力(例えば一瞬で効果部位濃度をゼロにする)を得ることはできません。

2014年3月22日

目次へ

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

Time limit is exhausted. Please reload CAPTCHA.