JP7323911B2 - 同期引込周波数帯域演算装置、同期引込周波数帯域演算方法およびプログラム - Google Patents

Description

本発明は、発振器などの注入同期性能を把握する同期引込周波数帯域演算装置、同期引込周波数帯域演算方法およびプログラムに関する。

発振器が取り付けられた集積回路などの電子回路においては、回路に取り付けられた発振器の注入同期性能の精確な把握が非常に難しいという問題があった。
注入同期とは、自励発振系に外部信号を強制注入すると、発振系が外部信号に同期する普遍的な物理現象である。これを利用する技術は真空管時代に端を発し、現在のミリ波等の高周波数帯での利用、省電力設計、回路の微細化の要請から、さらに進展している。例えば、無線通信信号の低位相雑音化や無線電力伝送の安定化のための要素技術として研究が進められている。

注入同期を応用する上では、発振器の同期性能を精確に評価することが基本的な課題である。同期性能は、注入同期が成立する外部信号の周波数帯域（同期引込周波数帯域）で評価される。
同期引込周波数帯域を得る手法については古くから膨大な研究がある。代表的な従来手法としては、以下の３つの手法（Ａ），（Ｂ），（Ｃ）がある。

（Ａ）位相方程式
位相方程式は古くから知られるアドラーの方程式の一般化であり、位相感受関数と呼ばれる発振器の内部特性から同期引周波数帯域を得ることを可能とする。位相感受関数は、振動子に微小な注入インパルスが与えられた時に発生する僅かな位相のシフト量として、実験的にあるいは数値計算で得ることができる。位相方程式は、発振器の位相感受関数Ｚを用いて、以下の式で与えられる。

［数１］式において、Δωは発振器の自然角周波数ωと外部信号の角周波数Ωからの離調(Δω≡ω-Ω）、fは外部信号である。
ここで、同期引込周波数帯域は［数１］式の右辺のφZ(Ωt+ψ)・f(Ωt)dtを計算することで得ることができる。

一方で、この［数１］式で示す位相方程式は、次の［数２］式に示すより精確な位相の時間発展式（以下、「Winfree」モデルと呼ぶ）を近似して得られたものであり、［数２］式に示す同期引込周波数帯域と比較して精度が低いという問題点がある。

（Ｂ）直接シミュレーション
この手法は、発振器の時間発展式を数値的に直接解き、同期引込周波数帯域を求めるものである。
この直接シミュレーションの問題点は、計算コストが非常に高いことである。例えば、小規模の電子回路を計算する場合でも、ワークステーションを使用して数時間の計算時間が必要である。さらに、心臓の１拍をシミュレートする場合、地球シミュレータと称されるスーパーコンピュータを使って、数時間の計算が必要である。また、直接シミュレーションの場合、得られるデータ点そのものは信用できるが、それぞれのデータ点をつないで得られる同期引込周波数帯域の品質保証はない。さらに、同期系を記述する方程式がなければ、直接シミュレーションの手法を適用することは不可能である。

（Ｃ）Pikovskyの方法
Pikovskyの方法は、［数３］式の位相の時間発展式を用いて、同期引込周波数帯域を得るものである。

［数３］式における結合関数Ｑは、［数３］式のＺ・ｆに相当する。
Pikovskyの方法は、この結合関数を同期の「時系列データ」より推定・補間するものである。したがって、その精度の保証は数値的な推定・補間のコストに依存し、時系列データの選定と、その精度に依存する。このため、Pikovskyの方法は、同期引込周波数帯域の算出結果の精度を保証することはできない。
なお、特許文献１には、発振器を有する注入同期系に注入する入力信号の最適波形を演算する手法についての一例が記載されている。

特開２０１５－１４６５３４号公報

上述したように、従来の手法では、いずれも問題があった。すなわち、位相方程式によって同期引込周波数帯域を得る手法（Ａ）では、精度が低いまたは保証されないという問題があった。また、直接シミュレーションによって同期引込周波数帯域を得る手法（Ｂ）では、計算コストが高く、かつ発振器を記述する式が得られていないと適用不可能であるという問題があった。さらに、Pikovskyの方法によって同期引込周波数帯域を得る手法（Ｃ）では、精度の保証が不可能であった。

本発明は、少ない計算コストで精度の高い同期引込周波数帯域を得ることができる同期引込周波数帯域演算装置、同期引込周波数帯域演算方法およびプログラムを提供することを目的とする。

また本発明の同期引込周波数帯域演算装置は、発振器を備える注入同期系の同期引込周波数帯域を演算する同期引込周波数帯域演算装置において、発振器の出力を得る入力部と、入力部に入力した発振器出力に基づいて、発振器の同期引込周波数帯域を演算する演算部と、演算部が演算した同期引込周波数帯域を出力する出力部とを備える。
演算部は、発振器の発振信号の同期引込周波数帯域が所定の精度となる入力強度限界としてのリミットサイクルを、発振器の時間発展方程式を近似したWinfreeモデルの式をPoincare-Lindstedt法による数学的理論を用いて解くことで得、得られたリミットサイクルと、入力無しのリミットサイクルの距離である乖離度を算出すると共に、算出した各入力強度における乖離度から平均二乗誤差を求め、求めた乖離度の平均二乗誤差が急激に増加する箇所が存在しない領域である乖離度の線形領域を入力強度限界として取得する演算を行う。そして、演算部で取得した線形領域の入力強度限界を、発振器の同期引込周波数帯域が精確に得られる入力強度限界とする。

また、本発明の同期引込周波数帯域演算方法は、発振器を備える注入同期系の同期引込周波数帯域を演算する同期引込周波数帯域演算方法において、発振器の発振信号の同期引込周波数帯域が所定の精度となる入力強度限界としてのリミットサイクルを、発振器の時間発展方程式を近似したWinfreeモデルの式をPoincare-Lindstedt法による数学的理論を用いて解くことで得、得られたリミットサイクルと、入力無しのリミットサイクルの距離である乖離度を算出する乖離度算出処理と、乖離度算出処理で算出した各入力強度における乖離度から平均二乗誤差を求める平均二乗誤差算出処理と、平均二乗誤差算出処理で得た乖離度の平均二乗誤差が急激に増加する箇所が存在しない領域である乖離度の線形領域を入力強度限界としてを取得する入力強度限界取得処理と、を含む。

さらに、本発明のプログラムは、上記同期引込周波数帯域演算方法の各処理を実行する手順をコンピュータに実行させるプログラムである。

本発明によると、従来よりも非常に少ない計算コストで、非常に精度の高い同期引込周波数帯域を得ることができる。すなわち、位相方程式により同期引込周波数帯域を得る場合よりも高い精度で同期引込周波数帯域の算出ができると共に、直接シミュレーション法よりも大幅に少ない計算で同期引込周波数帯域の算出を行うことができる。

本発明の第１の実施の形態例による同期引込周波数帯域演算装置の例を示す構成図である。 本発明の第１の実施の形態例によるロックレンジ取得処理の流れを示すフローチャートである。 本発明の第１の実施の形態例による入力強度限界算出処理の流れを示すフローチャートである。 本発明の第１の実施の形態例による入力強度とロックレンジの乖離度の例を示す特性図である。 本発明の第１の実施の形態例による入力強度と平均二乗誤差の例を示す特性図である。 本発明の第１の実施の形態例によるロックレンジおよびリミットサイクルの算出処理の流れを示すフローチャートである。 本発明の第１の実施の形態例によるロックレンジ乖離度の算出処理の例を示すフローチャートである。 本発明の第１の実施の形態例による平均二乗誤差の算出処理の例を示すフローチャートである。 本発明の第１の実施の形態例によるＰＬ法でのロックレンジ算出処理の例を示すフローチャートである。 本発明の第１の実施の形態例を適用した回路構成例を示す回路図である。 本発明の第１の実施の形態例による同期引込周波数帯域演算精度と、従来手法（位相方程式）とを比較した特性図である。 本発明の第２の実施の形態例の適用例であるHodgkin-Huxleyモデルの例を示す図である。 本発明の第２の実施の形態例によるHodgkin-Huxleyモデルのさらに具体的な例を示す図である。 本発明の第２の実施の形態例による同期引込周波数帯域演算精度と、従来手法（位相方程式）とを比較した特性図である。

＜本発明の原理＞
以下、本発明の実施の形態例を順に説明する。
最初に、本発明の具体的な実施の形態例を説明する前に、本発明の原理について説明する。
先に説明したように、従来から知られている位相方程式による手法（Ａ）は、Winfreモデル（［数２］式）を近似して解き、同期引込周波数帯域を得るものである。したがって、位相方程式による手法よりも、Winfreeモデルを直接解く方がより精確な同期引込周波数帯域が得られる。ここで本発明では、Poincare-Lindstedt法（以下、「ＰＬ法」と称する）という古くから知られる数学的理論を用いて、Winfreeモデルを解くことを特徴の１つとする。

一方で、Winfreeモデルは，本来の発振系の時間発展方程式を近似して得られた式である。したがって、元の式を直接計算して得られる手法（Ｂ）の直接シミュレーションで得られる同期引込周波数帯域とは誤差が生じてしまう。
ここで本発明では、同期引込周波数帯域を得る際の条件の設定により、ＰＬ法に生じる誤差を数%以内にとどめるようにしている。

まず、本発明に適用されるＰＬ法について説明する。
ＰＬ法は、先に説明したWinfreeモデル（［数２］式）に従い、同期しているときの位相θと周波数Ω（すなわち同期可能な周波数Ω）を精度良く与える手法である。
まず、［数２］式において、次の［数４］式に示すように、変数変換τ=Ωtを行う。［数４］式では、入力信号の初期位相φを導入した。

そして、次の［数５］式に示すように、位相θと周波数Ωを入力強度εで展開する。

このとき、発振器は同期しているという条件、すなわち各θ_nは周期的という条件から、次の［数６］式が成り立つ。

［数５］式を［数４］式に代入することで、次の［数７］式が得られる。

ここで、［数７］式の両辺でε＾n (n=0,1,2,…)の係数が一致するという条件、および［数６］式から，θ＾((n))とΩ＾((n))は、以下の［数８］式から順次得ることができる。

そして、［数８］式を［数５］式に代入して得たΩから、同期引込周波数帯域が得られる。

次に、本発明での、同期引込周波数帯域を得る際の条件の設定について説明する。
本発明では、ＰＬ法で得られた同期引込周波数帯域の精度数％が保証される入力強度限界ε_ＬＲは、リミットサイクルと呼ばれる周期軌道の乖離度d(ε)が線形的に振る舞う入力強度限界ε_ｌｉｎから決まることを利用している。

ここで、乖離度ｄ（ε）は外部入力信号無しのリミットサイクルと外部入力信号強度εにおける同期引込周波数帯域の端（最大・最小）の外力周波数におけるリミットサイクルの距離である。
乖離度ｄ（ε）を得るためには、まず、それぞれのリミットサイクル上に位相ψ＝０の原点Ｏ_１，Ｏ_２を任意に取り、リミットサイクル上の点を、Ｆ_Ｏ１（ψ），Ｇ_Ｏ２（ψ）と表す。
このとき、２つのリミットサイクルの平均距離を、次の［数９］式で示す。

乖離度d(ε)は、次の［数１０］式で定義される。

ＰＬ法で得られた同期引込周波数帯域が、リミットサイクルの乖離度と関係があることは、伝統的な同期理論から説明可能である。

ここで、本発明の同期引込周波数帯域を得るアルゴリズムの概略を説明する。
最初に、いくつかの入力強度{ε_1<ε_2<…<ε_n}で、直接シミュレーション手法により乖離度d(ε)を得る乖離度算出処理を行う。
次に、乖離度d(ε)を、[0,ε_ｍ](m=1,…,n)の範囲で線形フィッティングし、平均二乗誤差E(ε_ｍ)を計算する平均二乗誤差算出処理を行う。

その後、平均二乗誤差E(ε_ｍ)(m=1,…,n)が急激に増加するｍが存在するかを判定する。
ここまでの処理で、平均二乗誤差E(ε_ｍ) が急激に増加するｍが存在しない場合には、ε_ｎ＋１＞ε_ｎとなるε_ｎ＋１を、｛ε_１＜ε_２＜…＜ε_ｎ｝に追加して、ｎ＋１→ｎと定義し直す処理を行う。
そして、平均二乗誤差E(ε_ｍ)(m=1,…,n)が急激に増加するｍが存在する場合には、ε_ｍが乖離度d(ε)の線形領域の入力強度限界ε_ｌｉｎとなる。このようにして、入力強度限界取得処理が行われる。
このようにして得られた入力強度限界ε_ｌｉｎが、ＰＬ法で得られた同期引込周波数帯域の精度が数％に保証される入力強度限界ε_ＬＲと一致するようになる。

＜同期引込周波数帯域演算装置の構成例＞
以下、ここまで説明したアルゴリズムを適用した、本発明の実施の形態例の同期引込周波数帯域演算装置の詳細を順に説明する。
図１は、本発明の第１の実施の形態例に係る同期引込周波数帯域演算装置の構成例を示した図である。図１に示すように、同期引込周波数帯域演算装置１０は、注入同期系２０および高周波変換器３０に接続されている。
注入同期系２０は、信号源２１、入力信号生成器２２および発振器２３から構成されている。

信号源２１は、発振器２３に注入する入力信号の元になる所定周波数の交流信号を発生する。入力信号生成器２２は、例えばデジタル／アナログ変換器または入力信号の波形（電圧）を生成するハードウェアである。発振器２３は、同期引込を行わせるための通常用いられる発振器であり、例えばＥ級発振器などを含む。
高周波変換器３０は、発振器２３から供給される発信信号の周波数を無線通信用の高周波数に変換し、変換した高周波信号をアンテナ３１に供給するためのハードウェアである。

同期引込周波数帯域演算装置１０は、バス１１に接続されるＲＡＭ１２、ＲＯＭ１３、演算部１５、制御部１６、入出力インターフェース（Ｉ／Ｏ）１７ａ、１７ｂ、入力部１８および表示部１９より構成される。また、バス１１には不揮発性記憶装置（データベース）１４が接続されている。このデータベース１４には、発振器の情報５０が書き込まれる。具体的には、発振器の情報５０として、演算部１５で算出される入力強度限界ε_linが書き込まれる。また、データベース１４に書き込まれる発振器の情報（テーブル）５０として、位相感受関数Ｚ（θ）や、自然角周波数ωがある。なお、ＲＯＭ１３の容量が十分にある場合には、ＲＯＭ１３をデータベース１４として、ＲＯＭ１３に発振器の情報５０を記憶してもよい。

演算部１５は、主に、ＰＬ法によるロックレンジの算出を行うとともに、線形領域算出のための直接シミュレーションを行っている。制御部１６は、例えばＣＰＵ（Central Processing Unit）等の演算制御装置で構成され、同期引込周波数帯域演算装置１０を構成する各部の制御を行う。
同期引込周波数帯域演算装置１０の演算部１５で演算した結果は、入出力インターフェース１７ａを介して注入同期系２０の入力信号生成器２２に供給される。また、注入同期系２０の発振器２３の出力信号は、入出力インターフェース１７ｂを介して同期引込周波数帯域演算装置１０に供給され、演算部１５におけるロックレンジの算出および線形領域算出に利用される。

また、入力部１８からは、発振器の自然角周波数ω、発振器の位相感受関数Ｚ（θ）（以上はデータベースに情報がない場合）、所望の外部入力信号波形u(τ)（例えばsinτ）、初期の｛ε₁＜ε₂＜ε_３，・・・＜ε_n｝および所望のロックレンジの精度Δ_＊が入力され、演算部１５の演算に利用される。表示部１９には、演算部１５で算出されたロックレンジのデータ（ε、Δω_left、Δω_right）が表示される。

次に、図２～図９を参照して、本実施の形態例における注入同期引込周波数帯域を高精度に計算する方法について説明する。

＜ロックレンジの取得処理＞
最初に図２を参照して、ロックレンジ（ＬＲ）の取得処理（メイン）について説明する。図２に示す処理は、ロックレンジの左端Δω_leftと右端Δω_rightの一方についてのみ示しているが、この処理はロックレンジの左端Δω_leftと右端Δω_rightの両側で行われる処理である。なお、図面では、ロックレンジはＬＲと記載する。

ロックレンジの取得処理が開始されると（ステップＳ１０）、所望されるロックレンジの精度Δ_＊が入力される（ステップＳ１１）。このロックレンジの精度Δ_＊は、あくまでも同期引込周波数帯域の計算実行を指示するユーザが所望する精度値であり、ユーザが任意に入力する値である。例えば、精度Δ_＊を数％以下にしたい場合には、「数％」が所望Δ_＊の精度になる。

続いて、発振器２３の自然角周波数ωと発振器の位相感受関数Ｚ（θ）が取得される（ステップＳ１２）。自然角周波数ωと位相感受関数Ｚ（θ）は発振器２３の固有の数値であり、ロックレンジを計算しようとしている発振器自体が具備するものである。

ステップＳ１１とステップＳ１２で、ユーザにより所望の精度Δ_＊が入力され、発振器２３の自然角周波数ωと位相感受関数Ｚ（θ）が取得されると、次に、入力強度限界ε_linをデータベース１４から読み出す（ステップＳ１３）。ここで、入力強度限界ε_linは、リミットサイクル（ＬＣ）の乖離度ｄ（ε）が線形的に振る舞う入力強度限界である。この入力強度限界ε_linの値は、先に説明したＰＬ法から得られるロックレンジの精度、例えば数％が保証される入力強度限界ε_LRとほぼ等しい値となる。なお、入力強度εとは、外部入力信号ｆの振幅のことであり、振幅と波形を分離して表現するためにｆ(Ωt)をεu(Ωt)と表すこともある。また、図面では、リミットサイクルはＬＣと称する。

リミットサイクル（ＬＣ）の乖離度ｄ（ε）とは、ロックレンジの右端または左端の外部入力信号ｆ(Ωt)の周波数のＬＣと、外部信号がないときのリミットサイクルとの距離を表す。なお、ＰＬ法で得られる同期引込周波数帯域がリミットサイクルの乖離度ｄ（ε）と関係があることは、伝統的な同期理論から説明することができる。

ステップＳ１３で、入力強度限界ε_linをデータベース１４から読み出すことができるのは、当然のことながら、入力強度限界ε_linがデータベース１４に保存されている場合に限られる。つまり、発振器２３に対して線形領域が予め分かっている場合には、入力強度限界ε_linがデータベース１４に保存されているので、入力強度限界ε_linをデータベース１４から読み出すことが可能である。そうでない場合、つまり発振器２３に対して線形領域が予め分かっていない場合には、コストをかけても入力強度限界ε_linを計算で求めなければならない。

そこで、次のステップでは、ステップＳ１３においてデータベース１４から入力強度限界ε_linが読み出すことができたか否かが判断される（ステップＳ１４）。ステップＳ１４で入力強度限界ε_linがデータベース１４から読み出すことができたと判断された場合（ステップＳ１４のＹｅｓ）には、所望の外部入力信号波形u(τ)がユーザによって入力される（ステップＳ１５）。そして、入力した外部入力波形を振幅により正規化する（ステップＳ１６）。この正規化は入力信号波形u（τ）をu（τ）の最大値で割り算することによって行われる。

次に、入力強度限界ε_linをよりも小さい外部信号の入力強度（振幅に相当）εに対して、ＰＬ法によりロックレンジを算出する（ステップＳ１７）。このステップＳ１７の処理の詳細については図９で後述する。
そして、外部入力信号の入力強度εと、発振器２３の角周波数と外部入力信号の角周波数差Δωの組み合わせ（ε、Δω）を出力して（ステップＳ１８）、ロックレンジ取得処理（メイン）を終了する（ステップＳ１９）。

ステップＳ１４で、入力強度限界ε_linがデータベース１４から読み出されなかったと判断された場合には（ステップＳ１４のＮｏ）、外部入力信号の入力強度ε｛ε₁＜ε₂＜ε_３，・・・＜ε_n｝の中から、ユーザが適当な入力強度εを選択する（ステップＳ２０）。このとき、選択する作業を行うユーザに、該当する技術についてのスキルがある場合には、スキルに基づいて適切な入力強度εを選択するのが好ましいが、ここで選択する入力強度εは十分小さな値であれば適切でなくても問題ない。

そして、所望の外部入力信号波形u(τ)をサイン波として（u(τ)＝sinτ）入力強度限界ε_linの算出処理を行い（ステップＳ２１）、算出された入力強度限界ε_linをデータベース１４に書き込む（ステップＳ２２）。これ以降、ステップＳ１５～Ｓ１８までの処理が行われる。なおステップＳ２１の処理の詳細は、図３で後述する。

＜入力強度限界ε_linの算出処理＞
次に、図３のフローチャートを参照して、ステップＳ２１における入力強度限界ε_linの算出処理について説明する。
入力強度限界ε_linの算出処理が開始されると（ステップＳ２５）、外部入力信号の入力強度εをｎ段階に分割し、分割した入力強度ε_１～ε_ｎ｛ε₁＜ε₂＜ε_３，・・・＜ε_n｝に対してロックレンジを、ＬＣ算出処理によって求める（ステップＳ２６）。このステップＳ２６の処理の詳細は、図６で後述する。

次に、ステップＳ２６の算出処理で得られた入力強度εと出力電圧Vout(t)からリミットサイクル（ＬＣ）の乖離度ｄ(ε)を算出する（ステップＳ２７）。ここで図４を参照してリミットサイクルの乖離度ｄ(ε)について説明すると、横軸は入力強度ε₁、ε₂、ε₃、ε₄、ε₅を示し、縦軸はリミットサイクルの乖離度ｄ（ε）を示す。図４に示すように、×印で示した乖離度ｄ（ε）は、入力強度ε₃とε₄の間に、大きなギャップ（開き）が認められる。すなわち、入力強度ε₁～ε₃までは、実線で示した直線上に載っているが、入力強度ε₄、ε₅は、直線上にない。この図４から、入力強度限界ε_linが入力強度ε₃に近いことが分かる。
なお、ステップＳ２７の処理の詳細は、図７で後述される。

ステップＳ２７の処理の後、乖離度ｄ(ε_ｍ)(m=１,２・・・ｎ）について、平均二乗誤差E(ε_ｍ）を計算する（ステップＳ２８）。ここで、平均二乗誤差E(ε_ｍ）について、図５を参照して説明する。図４と同様に、図５の横軸も入力強度ε₁、ε₂、ε₃、ε₄、ε₅であり、縦軸は平均二乗誤差Ｅ（ε）である。図５に示すように、入力強度ε₃とε₄の間に平均二乗誤差Ｅ（ε）に大きなギャップがみられる。このことは、入力強度ε₁～ε₃までが線形領域であり、入力強度ε₄、ε₅は非線形領域なっていることを示しており、これによって入力強度限界ε_linは、入力強度ε₃に近いことが理解できる。このステップＳ２８の算出処理の詳細は、図８で後述される。

次に、ステップＳ２８で求めた平均二乗誤差E(ε_m）が、一つ前の段階の平均二乗誤差E(ε_ｍ－１）から所定の倍率（例えば２倍：E(ε_ｍ）=２E(ε_ｍ－１））になるｍが存在するか否かが判断される（ステップＳ２９）。ステップＳ２９で、平均二乗誤差E(ε_m）が一つ前の平均二乗誤差E(ε_ｍ－１）に対して所定倍率（例えば２倍）になっていると判断された場合（ステップＳ２９のＹｅｓ）には、平均二乗誤差がE(ε_ｍ）になったときの入力強度ε_mよりも一つ前の段階の入力強度ε_m-1を入力強度限界ε_linとし（ステップＳ３０）、入力強度限界ε_linの算出処理を終了する（ステップＳ３１）。

ステップＳ２９で、平均二乗誤差E(ε_ｍ）が一つ前の平均二乗誤差E(ε_ｍ－１）に対して所定倍率（例えば、２倍）になっていないと判断された場合（ステップＳ２９のＮｏ）には、ＰＬ法によるロックレンジ算出処理がなされる（ステップＳ３２）。なお、ステップＳ３２の処理の詳細は、図９で後述する。

そして、ステップＳ２６のロックレンジ、リミットサイクル算出処理により得られた（ε,Δω_sim）と、ステップＳ３２のＰＬ法によるＬＲ算出処理で得られた（ε,Δω_PL）に基づいて、ロックレンジの精度Δ(ε_m) (m=１,２・・・ｎ）を算出する（ステップＳ３３）。
次に、ステップＳ３３の算出処理の結果、ロックレンジ（ＬＲ）の精度Δ(ε_ｍ) (m=１,２・・・ｎ）が所望の精度Δ_＊（例えば、１％）より大きいか否かが判断される（ステップＳ３４）。

ステップＳ３４で、ロックレンジの精度Δ(ε_m）が所望の精度Δ_＊（例えば、１％）より大きいｍが存在する判断された場合（ステップＳ３４のＹｅｓ）には、ステップＳ３０に移行して、入力強度ε_ｍの一つ前の段階の入力強度ε_ｍ－１を入力強度限界ε_linとする。ステップＳ３４でこのような判断をしているのは、ロックレンジの精度Δ(ε）の線形領域が広がっていたとしても、精度Δ_＊（例えば、１％）以上に大きくなってしまったら、意味がなくなるからである。

ステップＳ３４で、ロックレンジの精度Δ(ε_m）が所望の精度Δ_＊（例えば、１％）より大きくなるｍが存在しないと判断された場合（ステップＳ３４のＮｏ）には、入力強度εに入力強度ε_nよりも大きな入力強度となるε_n+1｛ε₁＜ε₂＜ε_３，・・・＜ε_n＜ε_n＋１｝を追加し、ｎ＋１→ｎと定義し直して、ステップＳ２６のロックレンジ、リミットサイクル算出処理に戻る（ステップＳ３５）。以上が、入力強度限界ε_linの算出処理である。

＜ロックレンジ、リミットサイクル算出処理＞
次に、図６を参照して、図３のステップＳ２６のロックレンジ（ＬＲ）、リミットサイクル（ＬＣ）算出処理の詳細について説明する。
処理が開始されると（ステップＳ４０）、まず、［数１１］式に示す位相方程式から、発信器と外部入力信号の角周波数差（離調）Δωを算出する（ステップＳ４１）。

ステップＳ４１で、ロックレンジの左端を算出する場合には、ΔωがΔω_leftとなり、ロックレンジの右端を算出する場合には、ΔωがΔω_rightとなる。
次に、離調Δωをステップ幅ｈと置き、その最小の値ｈ_stopをΔω×10^-5に設定する（ステップＳ４２）。この設定では、ステップ幅ｈを当初の値から、1/2、1/4、・・・と小さくしていくようにし、ｈ_stopになったら、その段階でそれ以上小さくすることを止める。

次に、外部入力信号の角周波数ΩをΩ＝ω＋Δωに設定する（ステップＳ４３）。そして、角周波数Ωの外部入力信号を注入し、直接シミュレーションにより時間とともに変わる出力電圧Vout(t)を測定してメモリに記憶する（ステップＳ４４）。
ステップＳ４４で時間とともに変わる出力電圧Vout(t)を測定した後に、周波数カウンターを用いて周波数の時間変化を測定し、これをメモリに記憶する（ステップＳ４５）。そして、平均角周波数ω～（オメガチルダー）を計算する（ステップＳ４６）。次に、計算した平均角周波数ω～と外部信号の角周波数Ωとを比較して、同期が成立しているか否かを判断する（ステップＳ４７）。ここで、計算した平均角周波数ω～が外部信号の角周波数Ωに限りなく近くなったとき、つまりω～／Ωが限りなく「１」になった時に同期が成立したと判断する。

ステップＳ４７で、同期が成立したと判断された場合（ステップＳ４７のＹｅｓ）には、発信器と外部入力信号の角周波数差（離調）Δω（ステップ幅ｈ）がその最小の値ｈ_stopより小さいか否かが判断される（ステップＳ４８）。ステップＳ４８で角周波数差（離調）Δωに相当するステップ幅ｈが最小値ｈ_stopより小さいと判断された場合（ステップＳ４８のＹｅｓ）には、角周波数差（離調）Δωとその時の出力電圧Vout(t)をメモリに格納し（ステップＳ４９）、処理を終了して（ステップＳ５０）、図３のステップＳ２７の処理に移行する。

ステップＳ４７で、平均角周波数ω～と外部信号の角周波数Ωの比が「１」に近い値にならず、同期が成立したと判断されなかった場合（ステップＳ４７のＮｏ）には、ステップ幅ｈをｈ／２とすることにより、発信器と外部入力信号の角周波数差（離調）Δωを小さくして（ステップＳ５１）、再度ステップＳ４３からの処理を行う。
また、ステップＳ４８で、角周波数差（離調）Δωに相当するステップ幅ｈが最小値ｈ_stopより小さくなっていないと判断された場合（ステップＳ４８のＮｏ）には、角周波数差（離調）Δωを大きくして、再度ステップＳ４３からの処理を行う。以上が、ロックレンジ、リミットサイクル算出処理である。

＜ロックレンジ乖離度ｄ（ε）の算出処理＞
次に、図７を参照して、図３のステップＳ２７におけるリミットサイクル（ＬＣ）乖離度ｄ（ε）の算出処理について説明する。処理が開始されると（ステップＳ５５）、外部入力信号無しの場合の出力電圧の時間変化V_OUT ⁽⁰⁾(t)と、外部入力信号有りの場合の出力電圧の時間変化V_OUT ^(ε)(t)を求める。そして、V_OUT ⁽⁰⁾(t)をF(θ)、V_OUT ^(ε)(t)をG(θ)と置く(ステップＳ５６)。なお、θは外部入力信号の各周波数Ωに時間tを掛けたΩtを置き換えたもの、すなわちθ＝Ωtである。

ここで、リミットサイクル（ＬＣ）とは、直接シミュレーションにより得られる電圧等の時間変化から作られる、相空間（状態空間）上の閉軌道をいう。例えば、回路シミュレータ等により、数値的に得られた複数の出力電圧V_OUT1、V_OUT2の時間変化を使ってリミットサイクルを得ることができる。また、得られる出力電圧V_OUTが１つだけの場合は、時間の遅れをτとして、（V_OUT1(t)、V_OUT2(t)）＝（V_OUT(t)、V_OUT(t+τ)）の時間変化を使ってリミットサイクルを便宜的に得ることができる。

次に、位相φの初期値ｉを「０」に設定する(ステップＳ５７)。この初期値「ｉ＝０」は、０～２πまでの位相φをＮ等分したときの０番目の位相φの番号である。そして、位相φ＝ｉ／Ｎ＊２πとして（ステップＳ５８）、以後の処理を行う。
最初に、外部入力信号無しの場合の出力電圧の時間変化F(θ)から位相をφだけずらした、外部入力信号有りの場合の出力電圧の時間変化G(θ+φ)の差の絶対値を、相空間（状態空間）上の閉軌道で一周させて（０～２π）積分した値ｄ_p（φ）を求める（ステップＳ５９）。

そして、位相φの０から２πまでをＮ等分して、ステップＳ５８、Ｓ５９の計算を繰り返す（ステップＳ６０）。すなわち、ステップＳ６０で「ｉ＜Ｎ」であれば（ステップＳ６０のＹｅｓ）、ステップＳ５８、Ｓ５９を繰り返し、ステップＳ６０で「ｉ≧Ｎ」になったと判断された場合（ステップＳ６０のＮｏ）には、ステップＳ５９で計算したｄ_ｐ（φ）（０~２π）の中の最小値を求めてｄ（＝min d_p(φ)）とする（ステップＳ６１）。そして、その時の入力強度εとｄのセット（ε、ｄ）をメモリに格納して（ステップＳ６２）、処理を終了する（ステップＳ６３）。この（ε、ｄ）がリミットサイクルで求めた乖離度ｄ（ε）に相当する。

＜平均二乗誤差Ｅ(ε_m)の算出＞
次に、図８を参照して、図３のステップＳ２８における平均二乗誤差Ｅ(ε_m)の算出処理について説明する。処理を開始して（ステップＳ６５）、初期値ｍを１にセットする（ステップＳ６６）。このｍ＝１は、入力強度εと乖離度ｄ（ε）を｛ε₁＜ε₂＜・・・＜ε_m｝のｍ個とったときの第１番目のデータを意味している。

そして、乖離度のデータ｛（ε_1、ｄ(ε_１)）、（ε_２、ｄ(ε_２)）、・・（ε_ｍ、ｄ(ε_ｍ)）｝を求め、原点（０,０）を通ると仮定して線形フィッティングを行って、１次関数d_lin（ε）を求める（ステップＳ６７）。このステップＳ６７の線形フィッティングを図に示したものが、既に説明した図４である。

次に、d（ε_ｉ）とｄ_lin（ε_ｉ）の差分をとり、ｉ=１～ｍ（ｍ個）の二乗平均値E(ε_m)を求める（ステップＳ６８）。そして、ステップＳ６８の計算をｍ＝ｎになるまで繰り返す(ステップＳ６９)。すなわち、ステップＳ６９で、ｍ≦ｎであれば（ステップＳ６９のＮｏ）、ステップＳ６７に戻り、ｍ＞ｎになったときに、そのときの平均二乗誤差Ｅ（ε_m）をメモリに格納して（ステップＳ７０）、処理を終了する(ステップＳ７１)。以上が、平均二乗誤差Ｅ（ε_m）の算出処理である。

＜ＰＬ法ロックレンジ算出処理＞
次に、図９を参照して、図２のステップＳ１７および図３のステップＳ３２のＰＬ法ロックレンジ算出処理について説明する。ＰＬ法ロックレンジ算出処理は、［数１２］式で示すWinfreeモデルの微分方程式の解を求める方法である。つまり、Winfreeモデルの微分方程式を解くために、Poincare-Lindstedt法に基づき、発信器の角周波数Ωと位相θを入力強度εのべき乗の展開式（［数１３］式）で表す方法である。ＰＬ法ＬＲ算出処理では、まず［数１３］式で示した角周波数Ωと位相θを入力強度εのべき乗の式において、入力強度εと入力強度εの２乗（ε²）の係数を求める。

処理が開始されると（ステップＳ７５）、まず、位相φの初期値ｉが「０」に設定される(ステップＳ７６)。この初期値「ｉ＝０」は、０～２πまでの位相φをＮ等分したときの０番目の位相φのことである。そして、位相φ＝ｉ／Ｎ＊２πとして（ステップＳ７７）、以後の処理を行う。

既に説明したように、発信器の角周波数Ωを入力強度εのべき乗の展開式（［数１３］式）で表したときの入力強度εの係数を順にΩ⁽⁰⁾(φ)、Ω⁽¹⁾(φ)、Ω⁽²⁾(φ)、・・・とすると、０次の係数はΩ⁽⁰⁾(φ)＝ωになり、１次の係数Ω⁽¹⁾(φ)は、［数１４］式によって求めることができる（ステップＳ７８）。

また、位相θ(τ)も［数１３］式に示すように、入力強度εのべき乗の展開式で表すことができるが、そのときの０次の係数はθ⁽⁰⁾(τ)＝τとなる。なお、τは外部入力信号の角周波数Ωに時間ｔを掛けたΩtを置き換えたもの、すなわちτ＝Ωtである。

位相θ(τ)を入力強度εのべき乗で展開した［数１３］式におけるεの係数は、［数１５］式で求められる（ステップＳ７９）。

さらに、発信器の角周波数Ω(φ)の入力強度εの２乗（ε^２）の係数Ω⁽²⁾(φ)は、［数１６］式で求めることができる（ステップＳ８０）。

以上のステップＳ７８～Ｓ８０の計算がｉ＝０からｉ＝Ｎまで繰り返され、位相φが０～２πの間でＮ分割されたそれぞれの値に対して行われる（ステップＳ８１）。すなわち、ステップＳ８１で、ｉ＜Ｎ（φ＜２π）であれば（ステップＳ８１のＹｅｓ）、ｉがＮ（φ＝２π）になるまで、ステップＳ７７～Ｓ８０の処理が繰り返される。そして、ステップＳ８１でｉ≧Ｎになったと判断された場合に（ステップＳ８１のＮｏ）、次の入力強度限界ε_linを求める処理に移行する。

ここでも、入力強度εを求めるための「ｊ」が「０」に設定される(ステップＳ８２)。ここで「ｊ＝０」は、０から入力強度限界ε_linまでの入力強度εをＭ段階で分割したときの初期値を示している。また、入力強度ε_jは、入力強度限界ε_linを使って、ε_j＝ｊ／Ｍ＊ε_linで表される（ステップＳ８３）。この式は、入力強度限界ε_linをＭ分割したときのｊ（ｊ＝０～Ｍ）番目の入力強度ε_jを示す式である。

そして、ステップＳ７８で求めたΩ⁽¹⁾(φ)と、ステップＳ８０で求めたΩ⁽²⁾(φ)を使って、外部入力信号Ω（φ）を近似式で求める（ステップＳ８４）。このステップＳ８４の式は、発信器と外部入力信号の角周波数差（離調）Δω（＝Ω（φ）－ω）が入力強度εの２次式で表されることを意味しており、１次係数がΩ⁽¹⁾(φ)、２次係数がΩ⁽²⁾(φ)である。

次に、φ（０～２π）間の最小値の外部信号の角周波数Ω（φ）（minΩ（φ））と発振器の角周波数ωとの差を求め、ロックレンジの左端Δω_leftを算出する。さらに、φ（０～２π）間の最大値の外部信号の角周波数Ω（φ）（maxΩ（φ））と発振器の角周波数ωとの差を求め、ロックレンジの右端Δω_rightを算出する（ステップＳ８５）。

以上のステップＳ８３～Ｓ８５の計算をｊ＝０からｊ＝Ｍまで繰り返し、入力強度εの計算を０から入力強度限界ε_linの間でＭ分割されたそれぞれの値に対して行う（ステップＳ８６）。すなわち、ステップＳ８６で、ｊ＜Ｍ（ε_j＜ε_lin）であれば（ステップＳ８８のＹｅｓ）、ｊがＭ（ε_j＝ε_lin）になるまで、ステップＳ８３～Ｓ８５の処理が繰り返される。そして、ステップＳ８６でj≧Ｍになったと判断された場合には（ステップＳ８８のＮｏ）、（ε、Δω）を記憶部に格納して（ステップＳ８７）処理を終了する（ステップＳ８８）。そして、図２のステップＳ１８に移行する。

＜発振器の構成例＞
図１０は、図１に示す発振器２３の具体的な構成例を示す図である。ここでは、発振器２３として、リングオシレータとした例である。
信号源（パルスジェネレータ）２１からのパルスが入力信号生成器（バイアスサーキット）２２に供給され、入力信号生成器２２で生成された信号が、発振器（リングオシレータ）２３の入力端子２３ａに供給される。
発振器２３は、３段に直列接続されたインバータ２３ｂ、２３ｃ、２３ｄとスイッチング素子２３ｆとを備え、３段目のインバータ２３ｄの出力信号が、１段目のインバータ２３ｂの入力に戻されると共に、発振器２３の出力として出力端子２３ｅから出力される。

入力端子２３ａに得られる信号は、スイッチング素子２３ｆのゲートに供給される。スイッチング素子２３ｆのドレインは、２段目のインバータ２３ｃと３段目のインバータ２３ｄの間に接続され、スイッチング素子２３ｆのソースは、１段目のインバータ２３ｂの入力側に接続される。
本実施の形態例では、このように注入同期発振器として構成された発振器２３の同期引込周波数帯域を、高精度に算出することができる。

＜同期引込周波数帯域演算精度の例＞
図１１は、本実施の形態例による同期引込周波数帯域演算装置１０で、図１０に示す発振器（リングオシレータ）２３の同期引込周波数帯域を演算した場合を示す。
この図１１において、従来方法と記載された特性は、位相方程式で求めた特性であり、本発明と記載された特性は、本実施の形態例による同期引込周波数帯域演算装置１０で求めた特性である。

図１１の（Ａ）は、ロックレンジの左端におけるリミットサイクル（ＬＣ）の乖離度ｄ_left（ε）を示し、図１１の（Ｂ）は、ロックレンジの右端におけるリミットサイクル（ＬＣ）の乖離度ｄ_right（ε）を示す。
図１１の（Ｃ）は、ロックレンジ（ＬＲ）の左端における精度Δ_left（ε）を示し、図１１の（Ｄ）は、ロックレンジ（ＬＲ）の右端における精度Δ_right（ε）を示す。図１１の（Ｃ），（Ｄ）に示す△印の特性は、位相方程式（従来方法）で求めた精度であり、○印の特性は、本実施の形態例により求めた精度である。これらの精度の数値はパーセントである。
図１１において、横軸は外部入力信号の入力強度（電圧振幅ε［Ｖ］）である。

なお、図１１の（Ａ）および（Ｃ）は、外部入力信号の角周波数を、入力強度毎に注入同期が成立する最小のものにした場合であり、図１１の（Ｂ）および（Ｄ）は、外部入力信号の角周波数を、入力強度毎に注入同期が成立する最大のものにした場合である。
ここでは、図１１の（Ａ）および（Ｂ）に示すデータ点（□印）がリミットサイクル（ＬＣ）の乖離度であり、図１１の（Ｃ）および（Ｄ）に示すデータ点（△印）が従来方法（位相方程式）により得られたロックレンジ（ＬＲ）の精度、データ点（○印）が本実施の形態例（本発明）により得られたロックレンジ（ＬＲ）の精度である。

図１１の（Ａ）および（Ｂ）から判るように、入力強度ε＝０．０９０［Ｖ］およびε＝０．１０５［Ｖ］が、ロックレンジ（ＬＲ）の左端と右端における、それぞれのリミットサイクル（ＬＣ）の乖離度が線形的に振る舞う入力強度限界ε_linである。本実施の形態例の場合、図１１の（Ｃ）および（Ｄ）に示すように、ロックレンジ（ＬＲ）の精度が従来方法（位相方程式）よりも優れており、入力強度がε≦ε_linであれば、予め設定した精度１％よりも良くなることがわかる。

次に、同期引込周波数帯域の算出にかかる時間（ＣＰＵ時間）を、直接シミュレーション法と位相方程式（従来手法）と本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合で比較する。
直接シミュレーション法では、Synopsys社の回路シミュレータHSPICEをインテル社製の2.6GHz×12コアXeonプロセッサE5-2630v2（ただし、デル社のPRECISION T5610を使用）で実行した。
位相方程式（従来手法）および本実施の形態例における同期引込周波数帯域演算装置１０（本発明）では、C言語で実装したプログラムをインテル社製の3.6GHz×8コアXeonプロセッサCore i7-4790（ただし、デル社のOPTIPLEX9020を使用）で実行した。

同期引込周波数帯域の算出にかかる時間（ＣＰＵ時間）は、
直接シミュレーション法の場合は約６時間であり、
位相方程式（従来手法）の場合は約１μ秒であり、
本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合は約２～３μ秒である。

本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合の計算時間は、直接シミュレーション法の場合よりも非常に短く、従来手法である位相方程式の場合と比べても同等程度であることがわかる。

このように本実施の形態例の同期引込周波数帯域演算装置１０によると、従来よりも非常に少ない計算コストで、非常に精度の高い同期引込周波数帯域を得ることができるようになる。すなわち、位相方程式により同期引込周波数帯域を得る場合よりも高い精度で算出ができると共に、直接シミュレーション法よりも大幅に少ない計算で行える。

＜本発明の第２の実施の形態例＞
次に、本発明の第２の実施の形態例を、図１２～図１４を参照して説明する。
本発明の第２の実施の形態例では、生体の電気活動を算出する同期引込周波数帯域演算装置としたものである。同期引込周波数帯域演算装置としての構成や演算処理は、第１の実施の形態例で説明した同期引込周波数帯域演算装置１０と同じである。

図１２は、生体の神経細胞の電気活動を表す代表的モデルである、Hodgkin-Huxleyモデルを示す。図１２の左側に示すように、細胞体１０１は、核１０２、樹状突起１０３、シナプス１０４、軸索丘１０５、ミエリン鞘１０６、および軸索抹消１０７を備える。
ここで、図１２の右側に示すように、神経細胞の電気活動が伝導する箇所である神経線維１１０に、電位計１４０が接続された微細電極１２０および電極１３０を装着して、生体の神経細胞の電気活動による電位を計測する。つまり、生体が備える神経細胞を発振器として、このHodgkin-Huxleyモデルによる同期引込周波数帯域を、同期引込周波数帯域演算装置１０が演算する。

このHodgkin-Huxleyモデルは、図１２に示すような神経細胞に見られる活動電位について、現象論的に記述したｖ，ｍ，ｈ，ｎの4つの変数からなる連立方程式（［数１７］式）である。

［数１７］式の連立方程式において、各変数の意味は以下の通りである。
ｖ：細胞膜の電位（膜電位）
ｍ：Ｎａ^＋チャネルの活性化ゲートが開いている割合
ｈ：Ｎａ^＋チャネルの不活性化ゲートが開いている割合
ｎ：Ｋ^＋チャネルの活性化ゲートが開いている割合

図１３は、活性化ゲートや不活性化ゲートに生じるナトリウムイオン（Ｎａ^＋）やカリウムイオン（Ｋ^＋）を示すものである。

連立方程式でのｍ，ｈ，ｎの時間変化は、それぞれのゲートが開く速度定数αと閉じる速度定数βからなる方程式に従う。速度定数α，βは膜電位ｖに依存し、その依存性は、イカ巨大軸索を用いた実験値により便宜的に以下の［数１８］式で表される。

Hodgkin-Huxleyモデルでの、その他の定数は以下の通りである。
Ｃ_ｍ：細胞膜のコンデンサーとしての容量
Ｉ_ｄ：細胞内に注入する電流の直流成分
Ｉ（ωｔ）：細胞内に注入する電流の交流成分
また、各チャネルを流れる電流は、以下の式で定義される。

Ｉ_Ｎａ：Ｎａ^＋チャネル全体を流れる電流（３つの活性化ゲートと１つの不活性ゲートによって開閉される）
ｇ_Ｎａ：Ｎａ^＋のイオン電流の最大コンダクタンス
ｖ_Ｎａ：Ｎａ^＋のイオン電流がゼロになる平衡電位
Ｉ_Ｋ：Ｋ^＋チャネル全体を流れる電流（４つの活性化ゲートによって開閉される）
ｇ_Ｋ：Ｋ^＋のイオン電流の最大コンダクタンス
ｖ_Ｋ：Ｋ^＋のイオン電流がゼロになる平衡電位
Ｉ_Ｌ：その他のイオン電流（リーク電流と呼ばれる）
ｇ_Ｌ：リーク電流の最大コンダクタンス
ｖ_Ｌ：リーク電流がゼロになる平衡電位

このようなHodgkin-Huxleyモデルを適用して、第１の実施の形態例で説明した処理で、同期引込周波数帯域演算装置１０は、生体の同期引込周波数帯域を算出することができる。
図１４は、本実施の形態例による同期引込周波数帯域演算装置１０で、図１２に示す神経細胞の同期引込周波数帯域を演算した場合を示す。
この図１４において、従来方法と記載された特性は、位相方程式で求めた特性であり、本発明と記載された特性は、本実施の形態例で求めた特性である。

図１４の（Ａ）は、ロックレンジの左端におけるリミットサイクル（ＬＣ）の乖離度ｄ_left（ε）を示し、図１４の（Ｂ）は、ロックレンジの右端におけるリミットサイクル（ＬＣ）の乖離度ｄ_right（ε）を示す。
図１４の（Ｃ）は、ロックレンジ（ＬＲ）の左端における精度Δ_left（ε）を示し、図１１の（Ｄ）は、ロックレンジ（ＬＲ）の右端における精度Δ_right（ε）を示す。図１４の（Ｃ），（Ｄ）に示す△印の特性は、位相方程式（従来方法）で求めた精度であり、○印の特性は、本実施の形態例により求めた精度である。これらの精度の数値はパーセントである。
図１４において、横軸は外部入力信号の入力強度（電流振幅ε［Ａ］）である。

図１４の（Ａ）および（Ｃ）は、外部入力信号の角周波数を、入力強度毎に注入同期が成立する最小のものにした場合であり、図１４の（Ｂ）および（Ｄ）は、外部入力信号の角周波数を、入力強度毎に注入同期が成立する最大のものにした場合である。
ここでは、図１４の（Ａ）および（Ｂ）に示すデータ点（□印）がリミットサイクル（ＬＣ）の乖離度であり、図１４の（Ｃ）および（Ｄ）に示すデータ点（△印）が従来方法（位相方程式）により得られたロックレンジ（ＬＲ）の精度、データ点（○印）が本実施の形態例（本発明）により得られたロックレンジ（ＬＲ）の精度である。

図１４の（Ａ）および（Ｂ）から判るように、入力強度ε＝１．０［Ａ］が、ロックレンジ（ＬＲ）の左端と右端における、それぞれのリミットサイクル（ＬＣ）の乖離度が線形的に振る舞う入力強度限界ε_linである。本実施の形態例の場合、図１４の（Ｃ）および（Ｄ）に示すように、ロックレンジ（ＬＲ）の精度が従来方法（位相方程式）よりも優れており、入力強度がε≦ε_linであれば、予め設定した精度１％よりも良くなることがわかる。

次に、同期引込周波数帯域の算出にかかる時間（ＣＰＵ時間）を、直接シミュレーション法と位相方程式（従来手法）と本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合で比較する。
直接シミュレーション法、位相方程式（従来手法）および本実施の形態例における同期引込周波数帯域演算装置１０（本発明）では、いずれもC言語で実装したプログラムをインテル社製の3.6GHz×8コアXeonプロセッサCore i7-4790（ただし、デル社のOPTIPLEX9020を使用）で実行した。

同期引込周波数帯域の算出にかかる時間（ＣＰＵ時間）は、
直接シミュレーション法の場合は約２時間であり、
位相方程式（従来手法）の場合は約１μ秒であり、
本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合は約２～３μ秒である。

本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合の計算時間は、直接シミュレーション法の場合よりも非常に短く、従来手法である位相方程式の場合と比べても同等程度であることがわかる。

このように本実施の形態例の同期引込周波数帯域演算装置１０においても、従来よりも非常に少ない計算コストで、非常に精度の高い同期引込周波数帯域を得ることができるようになる。すなわち、位相方程式により同期引込周波数帯域を得る場合よりも高い精度で算出ができると共に、直接シミュレーション法よりも大幅に少ない計算で行える。

＜変形例＞
なお、上述した各実施の形態例は、本発明の好適な例を示すものであり、本発明の同期引込周波数帯域演算装置は、様々な発振器を備える注入同期系の同期引込周波数帯域の演算に用いることができる。また、同期引込周波数帯域を演算する際に、上述した実施の形態例では、求めた乖離度の平均二乗誤差の変化から乖離度の線形領域となる入力強度限界を取得する処理まで、演算で求めるようにした。これに対して、線形領域か否かの判断は、演算結果の表示出力結果からユーザが判断し、その判断した結果を同期引込周波数帯域演算装置に入力させるようにしてもよい。

また、図１に示す同期引込周波数帯域演算装置１０は、専用の演算装置で構成する他に、本発明の同期引込周波数帯域演算処理の各手順を実行するプログラムを作成して、そのプログラムのコンピュータ装置での実行で、同期引込周波数帯域演算装置として機能するようにしてもよい。

１０…同期引込周波数帯域演算装置、１１…バス、１２…ＲＡＭ、１３…ＲＯＭ、１４…データベース部、１５…演算部、１６…制御部、１７ａ，１７ｂ…入出力インターフェース（Ｉ／Ｏ）、１８…入力部、１９…表示部、２０…注入同期系、２１…信号源、２２…入力信号生成器、２３…発振器、３０…高周波変換器、３１…アンテナ、５０…テーブル、１０１…細胞体、１０２…核、１０３…樹状突起、１０４…シナプス、１０５…軸索丘、１０６…ミエリン鞘、１０７…軸索抹消、１１０…神経線維、１２０…微細電極、１３０…電極、１４０…電位計

Claims (5)

1. 発振器を備える注入同期系の同期引込周波数帯域を演算する同期引込周波数帯域演算装置において、
   前記発振器の出力を得る入力部と、
   前記入力部に入力した発振器出力に基づいて、前記発振器の同期引込周波数帯域を演算する演算部と、
   前記演算部が演算した同期引込周波数帯域を出力する出力部とを備え、
   前記演算部は、前記発振器の発振信号の同期引込周波数帯域が所定の精度となる入力強度限界としてのリミットサイクルを、発振器の時間発展方程式を近似したWinfreeモデルの式をPoincare-Lindstedt法を用いて解くことで得、得られたリミットサイクルと、入力無しのリミットサイクルの距離である乖離度を算出すると共に、算出した乖離度の平均二乗誤差を求め、求めた乖離度の平均二乗誤差が急激に増加する箇所が存在しない領域である乖離度の線形領域を入力強度限界として取得する演算を行い、
   前記演算部で取得した線形領域の入力強度限界を、前記発振器の同期引込周波数帯域の最小および最大とする
   同期引込周波数帯域演算装置。
2. 前記発振器は、リングオシレータである
   前記出力部は、前記リングオシレータの同期引込周波数帯域を出力する
   請求項１に記載の同期引込周波数帯域演算装置。
3. 前記発振器は、生体の神経細胞であり、
   前記出力部は、前記神経細胞の同期引込周波数帯域を出力する
   請求項１に記載の同期引込周波数帯域演算装置。
4. 発振器を備える注入同期系の同期引込周波数帯域を演算する同期引込周波数帯域演算方法において、
   前記発振器の発振信号の同期引込周波数帯域が所定の精度となる入力強度限界としてのリミットサイクルを、発振器の時間発展方程式を近似したWinfreeモデルの式をPoincare-Lindstedt法による数学的理論を用いて解くことで得、得られたリミットサイクルと、入力無しのリミットサイクルの距離である乖離度を算出する乖離度算出処理と、
   前記乖離度算出処理で算出した乖離度の平均二乗誤差を求める平均二乗誤差算出処理と、
   前記平均二乗誤差算出処理で得た乖離度の平均二乗誤差が急激に増加する箇所が存在しない領域である乖離度の線形領域を入力強度限界として取得する入力強度限界取得処理と、を含み
   前記入力強度限界取得処理で得た線形領域の入力強度限界を、前記発振器の同期引込周波数帯域の最小および最大とする
   同期引込周波数帯域演算方法。
5. 発振器を備える注入同期系の同期引込周波数帯域を演算するプログラムであり、
   前記発振器の発振信号の同期引込周波数帯域が所定の精度となる入力強度限界としてのリミットサイクルを、発振器の時間発展方程式を近似したWinfreeモデルの式をPoincare-Lindstedt法による数学的理論を用いて解くことで得、得られたリミットサイクルと、入力無しのリミットサイクルの距離である乖離度を算出する乖離度算出手順と、
   前記乖離度算出手順で算出した乖離度の平均二乗誤差を求める平均二乗誤差算出手順と、
   前記平均二乗誤差算出手順で得た乖離度の平均二乗誤差が急激に増加する箇所が存在しない領域である乖離度の線形領域を入力強度限界として取得する入力強度限界取得手順と、をコンピュータに実行させる
   プログラム。

Images