単極性矩形波の周波数帯域

我々の研究室では,電気刺激にパルス幅0.5 msの単極性矩形波を用いています.これの周波数帯域は次のように計算できます.

パルス幅がTで振幅が1/Tの矩形波y(t)のラプラス変換Y(s)は,
\frac{1}{T}\left(\frac{1}{s}-\frac{e^{-Ts}}{s}\right)
になります.

形式的にs=j\omegaを代入して絶対値を求めると,
|Y(j\omega)|=|\frac{\sin^2(T\omega/2)}{T\omega/2}|
になります.つまりsinc関数の絶対値です.

この関数はTが0.5 msのときに,

Graph1

になります.-3 dBの帯域で886 Hz位になります.

筋音の周波数帯域は,加速度では100 Hz以下(変位ではさらに低い)ですから,0.5 msの単極性矩形波は十分フラットな周波数特性を有しています(100 Hzのゲインは0.996).

工学的な説明は上記のとおりですが,電気生理学的には活動電位を1つ発生させればよく,活動電位の持続時間は数ミリ秒ですから,上記より帯域は狭くなるでしょうが,機械的性質を計測する上では問題にはならないと考えられます.活動電位で筋のモデルの力発生要素(f_1(t))がインパルス的に収縮力を発揮することを仮定しています.

fig3

Windkesselモデル

物理情報数学Bで,状態微分方程式の例として図の回路を取り上げました.

この回路図は,血液循環のモデルであるWindkesselモデルとして知られているものです.Windkesselモデルには,素子数によって2要素,3要素,4要素のモデルがあります.2要素はCとRの並列接続で簡単ですが,近似としては粗くなります.3要素が一番よく使われていると思います.4要素は,3要素のモデルでは近似しきれない低周波数での特性も良好に近似できます.

この回路から,状態微分方程式をたてる演習は,3年生のバイオシステムで心循環系の話題のときにももう一度行っています.

教室変更

毎年,2年生のデジタル基礎の講義を,1回だけはPCが設置されている教室で実施します.PCを使って実際にIP addressを確認したり,arp,netstat,tracertやnslookupを使ってインターネットの通信の仕組みを少しでも垣間見てもらうためです.日吉は第5校舎が3.11の震災で取り壊され,教室が減ったためか,教室の空きが少なく,一昨年までの教室を利用できず,別の教室になりました.教室が変われば,PCの設定状況も変わるので,昨日教材作成のためにその教室でPCを使ってみました.設定値が変わるだけで本質的な変更はなかったので助かりました.一昨年には,OSが変わり,IPv6の設定も追加されていて,tracertで追えない部分も増えていて大幅な教材の見直しが必要でした.定員が今年の履修者数に等しいの教室なので,履修者が今年より増えると対応できなくなります.

学科の研究室のネットワークの管理は,学生が行っているところが多いですから,上記のような内容を扱っておくことも必要であろうと考えています.ITCあるいは学科の管理者と会話が成り立たたないということは減ったのではないかと思います.

特性方程式と重根

今日の講義で微分方程式の特性根と一般解の話をしたら,重根をもつ場合について,高校で習う隣接3項漸化式で重根を持つ場合の解き方と同じ考え方で微分方程式も説明できることが分かったと講義の後で教えてくれました.以下にまとめてみます.

隣接3項漸化式で特性方程式が重根をもつ場合

y(k+2)=2ay(k+1)-a^2y(k)
を変形して,差が等比数列の形になるようにする.

y(k+2)-ay(k+1)=a(y(k+1)-ay(k))

右辺は
a^{k+1}(y(1)-ay(0))
となる.

両辺をa^{k+1}で割る.

  \displaystyle\frac{y(k+2)}{a^{k+1}}-\frac{y(k+1)}{a^k}=y(1)-ay(0)
となり,\displaystyle\frac{y(k+1)}{a^k}は,初項y(0)/a^{-1},公差y(1)-ay(0)の等差数列になる.

したがって
  \displaystyle\frac{y(k)}{a^{k-1}}=\frac{y(0)}{a^{-1}} + (k-1)(y(1)-ay(0))
となる.

両辺をa^{k-1}倍して
  y(k)=a^ky(0)+a^{k-1}(k-1)(y(1)-ay(0))
になる.

隣接3項漸化式では,y(k)を求めることになるが,線形差分方程式とみた場合には,基底にa^kka^kを用いればよいことが分かる.

微分方程式

微分方程式の特性根が重根の場合,
  \displaystyle\frac{d^2y(t)}{dt^2} -2a\frac{dy(t)}{dt} + a^2y(t)=0
を変形して
  \displaystyle\frac{d^2y(t)}{dt^2} - a\frac{dy(t)}{dt} = a\left(\frac{dy(t)}{dt} - ay(t)\right)
となる.
これは,
  \displaystyle\frac{d}{dt}\left(\frac{dy(t)}{dt}-ay(t)\right)= a\left(\frac{dy(t)}{dt} - ay(t)\right)
なので,
  \displaystyle\frac{dy(t)}{dt}-ay(t) = e^{at}\left(\frac{dy(0)}{dt} - ay(0)\right)
になる.
e^{at}で両辺を割ると,
  e^{-at}\displaystyle\frac{dy(t)}{dt}-ae^{-at}y(t) = \frac{dy(0)}{dt} -ay(0)
なので,左辺を次のように書き換えて
  \displaystyle\frac{d}{dt}\left(e^{at}y(t)\right) = \frac{dy(0)}{dt} -ay(0)
になる.
積分して,
  e^{-at}y(t) = \left(\displaystyle\frac{dy(0)}{dt} -ay(0)\right)t + c
になり(cは積分定数),
  y(t)=\left(\displaystyle\frac{dy(0)}{dt} -ay(0)\right)te^{at} + ce^{at}
になる,
つまり,基底はe^{at}te^{at}になる.

実際問題として,重根のときの基底を導出する時間はないでしょうから,重根のときの基底を覚えておくとよいでしょう.また,講義で説明した定数変化法を理解しておくとよいと思います.

線形と非線形

月曜日は2時限に物理情報数学B,4時限に生体制御を担当しています.

物理情報数学Bは線形代数です.先週提出の宿題,ランチェスターのモデルの問題は,講義中に例を板書した,α=0.1,β=0.1では,さすがにつまらないです.ノートを右から左に写すだけになってしまいます.自分でαとβの値を設定して,つまり問題を自ら設定できないと評価されません.ランチェスターのモデルでは固有値が簡単な値になるようにαとβを定めることは簡単なことです.

生体制御は,非線形振動を扱いました.BVPモデル,好中球の密度の振動,心循環系の圧反射のモデル,サーカディアンリズムです.私が大学院生のとき,非線形振動だけで1つの講義があり,van der Polの式を数値的に解いて,リミットサイクルを描き,また位相面にメッシュを切って各格子点の速度ベクトルを描くことが課題だったように思います.当時のPCは,8086 8 MHzやV30 10 MHzの16ビットでMS-DOSで動いていました.XYプロッターに出力するプログラムを書いた記憶があります.

7セグメント表示器

デジタル基礎の講義で7セグメント表示器の任意のセグメントを光らせる組み合わせ回路を設計する問題でした.Don’t careがあるカルノーマップの問題です.7セグメント表示器では,10〜15はdont’ careになります.

 

丁度,ネットで
1=2, 2=5, 3=5, 4=4…
という問題がアナログな頭では解けないというのを見かけたところでした.7セグメント表示の文字はやはりデジタルの代名詞なのかと思ったところです.