JP 7323911 B2 2023.8.9(57)【 特 許 請 求 の 範 囲 】【請求項１】
発振器を備える注入同期系の同期引込周波数帯域を演算する同期引込周波数帯域演算装置において、
前記発振器の出力を得る入力部と、前記入力部に入力した発振器出力に基づいて、前記発振器の同期引込周波数帯域を演算する演算部と、
前記演算部が演算した同期引込周波数帯域を出力する出力部とを備え、前記演算部は、前記発振器の発振信号の同期引込周波数帯域が所定の精度となる入力強度限界としてのリミットサイクルを、発振器の時間発展方程式を近似したWinfreeモデル

10の式をPoincare-Lindstedt法を用いて解くことで得、得られたリミットサイクルと、入力無しのリミットサイクルの距離である乖離度を算出すると共に、算出した乖離度の平均二乗誤差を求め、求めた乖離度の平均二乗誤差が急激に増加する箇所が存在しない領域である乖離度の線形領域を入力強度限界として取得する演算を行い、前記演算部で取得した線形領域の入力強度限界を、前記発振器の同期引込周波数帯域の最小および最大とする
同期引込周波数帯域演算装置。
【請求項２】
前記発振器は、リングオシレータである
前記出力部は、前記リングオシレータの同期引込周波数帯域を出力する

20

(2)

JP 7323911 B2 2023.8.9

請求項１に記載の同期引込周波数帯域演算装置。
【請求項３】
前記発振器は、生体の神経細胞であり、前記出力部は、前記神経細胞の同期引込周波数帯域を出力する
請求項１に記載の同期引込周波数帯域演算装置。
【請求項４】
発振器を備える注入同期系の同期引込周波数帯域を演算する同期引込周波数帯域演算方法において、
前記発振器の発振信号の同期引込周波数帯域が所定の精度となる入力強度限界としてのリミットサイクルを、発振器の時間発展方程式を近似したWinfreeモデルの式をPoincareLindstedt法による数学的理論を用いて解くことで得、得られたリミットサイクルと、入力無しのリミットサイクルの距離である乖離度を算出する乖離度算出処理と、前記乖離度算出処理で算出した乖離度の平均二乗誤差を求める平均二乗誤差算出処理と、前記平均二乗誤差算出処理で得た乖離度の平均二乗誤差が急激に増加する箇所が存在しない領域である乖離度の線形領域を入力強度限界として取得する入力強度限界取得処理と、を含み
前記入力強度限界取得処理で得た線形領域の入力強度限界を、前記発振器の同期引込周
波数帯域の最小および最大とする
同期引込周波数帯域演算方法。
【請求項５】
発振器を備える注入同期系の同期引込周波数帯域を演算するプログラムであり、前記発振器の発振信号の同期引込周波数帯域が所定の精度となる入力強度限界としてのリミットサイクルを、発振器の時間発展方程式を近似したWinfreeモデルの式をPoincareLindstedt法による数学的理論を用いて解くことで得、得られたリミットサイクルと、入力無しのリミットサイクルの距離である乖離度を算出する乖離度算出手順と、前記乖離度算出手順で算出した乖離度の平均二乗誤差を求める平均二乗誤差算出手順と、前記平均二乗誤差算出手順で得た乖離度の平均二乗誤差が急激に増加する箇所が存在しない領域である乖離度の線形領域を入力強度限界として取得する入力強度限界取得手順と、をコンピュータに実行させるプログラム。
【発明の詳細な説明】【技術分野】【０００１】
本発明は、発振器などの注入同期性能を把握する同期引込周波数帯域演算装置、同期引込周波数帯域演算方法およびプログラムに関する。
【背景技術】【０００２】
発振器が取り付けられた集積回路などの電子回路においては、回路に取り付けられた発振器の注入同期性能の精確な把握が非常に難しいという問題があった。
注入同期とは、自励発振系に外部信号を強制注入すると、発振系が外部信号に同期する
普遍的な物理現象である。これを利用する技術は真空管時代に端を発し、現在のミリ波等の高周波数帯での利用、省電力設計、回路の微細化の要請から、さらに進展している。例
えば、無線通信信号の低位相雑音化や無線電力伝送の安定化のための要素技術として研究が進められている。
【０００３】
注入同期を応用する上では、発振器の同期性能を精確に評価することが基本的な課題である。同期性能は、注入同期が成立する外部信号の周波数帯域（同期引込周波数帯域）で評価される。

10

20

30

40

50

(3)

JP 7323911 B2 2023.8.9

同期引込周波数帯域を得る手法については古くから膨大な研究がある。代表的な従来手法としては、以下の３つの手法（Ａ），（Ｂ），（Ｃ）がある。
【０００４】（Ａ）位相方程式
位相方程式は古くから知られるアドラーの方程式の一般化であり、位相感受関数と呼ばれる発振器の内部特性から同期引周波数帯域を得ることを可能とする。位相感受関数は、振動子に微小な注入インパルスが与えられた時に発生する僅かな位相のシフト量として、実験的にあるいは数値計算で得ることができる。位相方程式は、発振器の位相感受関数Ｚを用いて、以下の式で与えられる。
【０００５】【数１】【０００６】［数１］式において、Δωは発振器の自然角周波数ωと外部信号の角周波数Ωからの離調 (Δ ω ≡ ω -Ω ） 、 fは 外 部 信 号 で あ る 。
こ こ で 、 同 期 引 込 周 波 数 帯 域 は ［ 数 １ ］ 式 の 右 辺 の φ Z(Ω t+ψ )・ f(Ω t)dtを 計 算 す る
ことで得ることができる。
【０００７】
一方で、この［数１］式で示す位相方程式は、次の［数２］式に示すより精確な位相の時 間 発 展 式 （ 以 下 、 「 Winfree」 モ デ ル と 呼 ぶ ） を 近 似 し て 得 ら れ た も の で あ り 、 ［ 数 ２］式に示す同期引込周波数帯域と比較して精度が低いという問題点がある。
【０００８】【数２】【０００９】（Ｂ）直接シミュレーション
この手法は、発振器の時間発展式を数値的に直接解き、同期引込周波数帯域を求めるものである。
この直接シミュレーションの問題点は、計算コストが非常に高いことである。例えば、小規模の電子回路を計算する場合でも、ワークステーションを使用して数時間の計算時間が必要である。さらに、心臓の１拍をシミュレートする場合、地球シミュレータと称されるスーパーコンピュータを使って、数時間の計算が必要である。また、直接シミュレーシ
ョンの場合、得られるデータ点そのものは信用できるが、それぞれのデータ点をつないで得られる同期引込周波数帯域の品質保証はない。さらに、同期系を記述する方程式がなければ、直接シミュレーションの手法を適用することは不可能である。
【００１０】（ Ｃ ） Pikovskyの 方 法
Pikovskyの 方 法 は 、 ［ 数 ３ ］ 式 の 位 相 の 時 間 発 展 式 を 用 い て 、 同 期 引 込 周 波 数 帯 域 を 得るものである。
【００１１】

10

20

30

40

【数３】

【００１２】［数３］式における結合関数Ｑは、［数３］式のＺ・ｆに相当する。
Pikovskyの 方 法 は 、 こ の 結 合 関 数 を 同 期 の 「 時 系 列 デ ー タ 」 よ り 推 定 ・ 補 間 す る も の である。したがって、その精度の保証は数値的な推定・補間のコストに依存し、時系列デー
タ の 選 定 と 、 そ の 精 度 に 依 存 す る 。 こ の た め 、 Pikovskyの 方 法 は 、 同 期 引 込 周 波 数 帯 域 の

50

(4)

JP 7323911 B2 2023.8.9

算出結果の精度を保証することはできない。
なお、特許文献１には、発振器を有する注入同期系に注入する入力信号の最適波形を演算する手法についての一例が記載されている。
【先行技術文献】【特許文献】【００１３】【特許文献１】特開２０１５−１４６５３４号公報【発明の概要】【発明が解決しようとする課題】【００１４】
上述したように、従来の手法では、いずれも問題があった。すなわち、位相方程式によって同期引込周波数帯域を得る手法（Ａ）では、精度が低いまたは保証されないという問題があった。また、直接シミュレーションによって同期引込周波数帯域を得る手法（Ｂ）
では、計算コストが高く、かつ発振器を記述する式が得られていないと適用不可能であると い う 問 題 が あ っ た 。 さ ら に 、 Pikovskyの 方 法 に よ っ て 同 期 引 込 周 波 数 帯 域 を 得 る 手 法 （Ｃ）では、精度の保証が不可能であった。
【００１５】
本発明は、少ない計算コストで精度の高い同期引込周波数帯域を得ることができる同期
引込周波数帯域演算装置、同期引込周波数帯域演算方法およびプログラムを提供することを目的とする。
【課題を解決するための手段】【００１６】
また本発明の同期引込周波数帯域演算装置は、発振器を備える注入同期系の同期引込周
波数帯域を演算する同期引込周波数帯域演算装置において、発振器の出力を得る入力部と、入力部に入力した発振器出力に基づいて、発振器の同期引込周波数帯域を演算する演算部と、演算部が演算した同期引込周波数帯域を出力する出力部とを備える。
演算部は、発振器の発振信号の同期引込周波数帯域が所定の精度となる入力強度限界としてのリミットサイクルを、発振器の時間発展方程式を近似したWinfreeモデルの式をPoi
ncare-Lindstedt法による数学的理論を用いて解くことで得、得られたリミットサイクルと、入力無しのリミットサイクルの距離である乖離度を算出すると共に、算出した各入力強度における乖離度から平均二乗誤差を求め、求めた乖離度の平均二乗誤差が急激に増加する箇所が存在しない領域である乖離度の線形領域を入力強度限界として取得する演算を行う。そして、演算部で取得した線形領域の入力強度限界を、発振器の同期引込周波数帯
域が精確に得られる入力強度限界とする。
【００１７】
また、本発明の同期引込周波数帯域演算方法は、発振器を備える注入同期系の同期引込
周波数帯域を演算する同期引込周波数帯域演算方法において、発振器の発振信号の同期引込周波数帯域が所定の精度となる入力強度限界としてのリミットサイクルを、発振器の時
間発展方程式を近似したWinfreeモデルの式をPoincare-Lindstedt法による数学的理論を用いて解くことで得、得られたリミットサイクルと、入力無しのリミットサイクルの距離である乖離度を算出する乖離度算出処理と、乖離度算出処理で算出した各入力強度における乖離度から平均二乗誤差を求める平均二乗誤差算出処理と、平均二乗誤差算出処理で得た乖離度の平均二乗誤差が急激に増加する箇所が存在しない領域である乖離度の線形領域
を入力強度限界としてを取得する入力強度限界取得処理と、を含む。
【００１８】
さらに、本発明のプログラムは、上記同期引込周波数帯域演算方法の各処理を実行する
手順をコンピュータに実行させるプログラムである。
【発明の効果】【００１９】
本発明によると、従来よりも非常に少ない計算コストで、非常に精度の高い同期引込周

10

20

30

40

50

(5)

JP 7323911 B2 2023.8.9

波数帯域を得ることができる。すなわち、位相方程式により同期引込周波数帯域を得る場
合よりも高い精度で同期引込周波数帯域の算出ができると共に、直接シミュレーション法
よりも大幅に少ない計算で同期引込周波数帯域の算出を行うことができる。
【図面の簡単な説明】【００２０】【図１】本発明の第１の実施の形態例による同期引込周波数帯域演算装置の例を示す構成
図である。
【図２】本発明の第１の実施の形態例によるロックレンジ取得処理の流れを示すフローチ
ャートである。
【図３】本発明の第１の実施の形態例による入力強度限界算出処理の流れを示すフローチ
ャートである。
【図４】本発明の第１の実施の形態例による入力強度とロックレンジの乖離度の例を示す
特性図である。
【図５】本発明の第１の実施の形態例による入力強度と平均二乗誤差の例を示す特性図である。
【図６】本発明の第１の実施の形態例によるロックレンジおよびリミットサイクルの算出
処理の流れを示すフローチャートである。
【図７】本発明の第１の実施の形態例によるロックレンジ乖離度の算出処理の例を示すフローチャートである。
【図８】本発明の第１の実施の形態例による平均二乗誤差の算出処理の例を示すフローチ
ャートである。
【図９】本発明の第１の実施の形態例によるＰＬ法でのロックレンジ算出処理の例を示す
フローチャートである。
【図１０】本発明の第１の実施の形態例を適用した回路構成例を示す回路図である。
【図１１】本発明の第１の実施の形態例による同期引込周波数帯域演算精度と、従来手法（位相方程式）とを比較した特性図である。
【 図 １ ２ 】 本 発 明 の 第 ２ の 実 施 の 形 態 例 の 適 用 例 で あ る Hodgkin-Huxleyモ デ ル の 例 を 示 す図である。
【 図 １ ３ 】 本 発 明 の 第 ２ の 実 施 の 形 態 例 に よ る Hodgkin-Huxleyモ デ ル の さ ら に 具 体 的 な 例を示す図である。
【図１４】本発明の第２の実施の形態例による同期引込周波数帯域演算精度と、従来手法（位相方程式）とを比較した特性図である。
【発明を実施するための形態】【００２１】＜本発明の原理＞
以下、本発明の実施の形態例を順に説明する。
最初に、本発明の具体的な実施の形態例を説明する前に、本発明の原理について説明する。
先 に 説 明 し た よ う に 、 従 来 か ら 知 ら れ て い る 位 相 方 程 式 に よ る 手 法 （ Ａ ） は 、 Winfreモ
デル（［数２］式）を近似して解き、同期引込周波数帯域を得るものである。したがって、 位 相 方 程 式 に よ る 手 法 よ り も 、 Winfreeモ デ ル を 直 接 解 く 方 が よ り 精 確 な 同 期 引 込 周 波
数 帯 域 が 得 ら れ る 。 こ こ で 本 発 明 で は 、 Poincare-Lindstedt法 （ 以 下 、 「 Ｐ Ｌ 法 」 と 称 する ） と い う 古 く か ら 知 ら れ る 数 学 的 理 論 を 用 い て 、 Winfreeモ デ ル を 解 く こ と を 特 徴 の １
つとする。
【００２２】
一 方 で 、 Winfreeモ デ ル は ， 本 来 の 発 振 系 の 時 間 発 展 方 程 式 を 近 似 し て 得 ら れ た 式 で ある。したがって、元の式を直接計算して得られる手法（Ｂ）の直接シミュレーションで得られる同期引込周波数帯域とは誤差が生じてしまう。
ここで本発明では、同期引込周波数帯域を得る際の条件の設定により、ＰＬ法に生じる
誤 差 を 数 %以 内 に と ど め る よ う に し て い る 。

10

20

30

40

50

(6)

JP 7323911 B2 2023.8.9【００２３】
まず、本発明に適用されるＰＬ法について説明する。
Ｐ Ｌ 法 は 、 先 に 説 明 し た Winfreeモ デ ル （ ［ 数 ２ ］ 式 ） に 従 い 、 同 期 し て い る と き の 位
相θと周波数Ω（すなわち同期可能な周波数Ω）を精度良く与える手法である。
ま ず 、 ［ 数 ２ ］ 式 に お い て 、 次 の ［ 数 ４ ］ 式 に 示 す よ う に 、 変 数 変 換 τ =Ω tを 行 う 。 ［数４］式では、入力信号の初期位相φを導入した。
【００２４】【数４】【００２５】
そして、次の［数５］式に示すように、位相θと周波数Ωを入力強度εで展開する。
【００２６】【数５】【００２７】
こ の と き 、 発 振 器 は 同 期 し て い る と い う 条 件 、 す な わ ち 各 θ ̲nは 周 期 的 と い う 条 件 か ら、次の［数６］式が成り立つ。
【００２８】【数６】

20【００２９】［数５］式を［数４］式に代入することで、次の［数７］式が得られる。
【００３０】【数７】

30【００３１】
こ こ で 、 ［ 数 ７ ］ 式 の 両 辺 で ε ＾ n (n=0,1,2,… )の 係 数 が 一 致 す る と い う 条 件 、 お よ び［ 数 ６ ］ 式 か ら ， θ ＾ ((n))と Ω ＾ ((n))は 、 以 下 の ［ 数 ８ ］ 式 か ら 順 次 得 る こ と が で き る。
【００３２】【数８】【００３３】
そして、［数８］式を［数５］式に代入して得たΩから、同期引込周波数帯域が得られる。
【００３４】
次に、本発明での、同期引込周波数帯域を得る際の条件の設定について説明する。

50

(7)

JP 7323911 B2 2023.8.9

本発明では、ＰＬ法で得られた同期引込周波数帯域の精度数％が保証される入力強度限界 ε Ｌ Ｒ は 、 リ ミ ッ ト サ イ ク ル と 呼 ば れ る 周 期 軌 道 の 乖 離 度 d(ε )が 線 形 的 に 振 る 舞 う 入
力強度限界εｌｉｎから決まることを利用している。
【００３５】
ここで、乖離度ｄ（ε）は外部入力信号無しのリミットサイクルと外部入力信号強度εにおける同期引込周波数帯域の端（最大・最小）の外力周波数におけるリミットサイクルの距離である。
乖離度ｄ（ε）を得るためには、まず、それぞれのリミットサイクル上に位相ψ＝０の原点Ｏ１，Ｏ２を任意に取り、リミットサイクル上の点を、ＦＯ１（ψ），ＧＯ２（ψ）と表す。
このとき、２つのリミットサイクルの平均距離を、次の［数９］式で示す。
【００３６】【数９】

10【００３７】
乖 離 度 d(ε )は 、 次 の ［ 数 １ ０ ］ 式 で 定 義 さ れ る 。
【００３８】【数１０】【００３９】
ＰＬ法で得られた同期引込周波数帯域が、リミットサイクルの乖離度と関係があること
は、伝統的な同期理論から説明可能である。
【００４０】
ここで、本発明の同期引込周波数帯域を得るアルゴリズムの概略を説明する。
最 初 に 、 い く つ か の 入 力 強 度 {ε ̲1<ε ̲2<… <ε ̲n}で 、 直 接 シ ミ ュ レ ー シ ョ ン 手 法 に より 乖 離 度 d(ε )を 得 る 乖 離 度 算 出 処 理 を 行 う 。
次 に 、 乖 離 度 d(ε )を 、 [0,ε ｍ ](m=1,… ,n)の 範 囲 で 線 形 フ ィ ッ テ ィ ン グ し 、 平 均 二 乗
誤 差 E(ε ｍ )を 計 算 す る 平 均 二 乗 誤 差 算 出 処 理 を 行 う 。
【００４１】
そ の 後 、 平 均 二 乗 誤 差 E(ε ｍ )(m=1,… ,n)が 急 激 に 増 加 す る ｍ が 存 在 す る か を 判 定 す る。
こ こ ま で の 処 理 で 、 平 均 二 乗 誤 差 E(ε ｍ ) が 急 激 に 増 加 す る ｍ が 存 在 し な い 場 合 に は 、
εｎ＋１＞εｎとなるεｎ＋１を、｛ε１＜ε２＜…＜εｎ｝に追加して、ｎ＋１→ｎと定義し直す処理を行う。
そ し て 、 平 均 二 乗 誤 差 E(ε ｍ )(m=1,… ,n)が 急 激 に 増 加 す る ｍ が 存 在 す る 場 合 に は 、 εｍ が 乖 離 度 d(ε )の 線 形 領 域 の 入 力 強 度 限 界 ε ｌ ｉ ｎ と な る 。 こ の よ う に し て 、 入 力 強 度
限界取得処理が行われる。
このようにして得られた入力強度限界εｌｉｎが、ＰＬ法で得られた同期引込周波数帯
域の精度が数％に保証される入力強度限界εＬＲと一致するようになる。
【００４２】＜同期引込周波数帯域演算装置の構成例＞
以下、ここまで説明したアルゴリズムを適用した、本発明の実施の形態例の同期引込周
波数帯域演算装置の詳細を順に説明する。
図１は、本発明の第１の実施の形態例に係る同期引込周波数帯域演算装置の構成例を示した図である。図１に示すように、同期引込周波数帯域演算装置１０は、注入同期系２０および高周波変換器３０に接続されている。
注入同期系２０は、信号源２１、入力信号生成器２２および発振器２３から構成されて
いる。

30

40

(8)

JP 7323911 B2 2023.8.9【００４３】
信号源２１は、発振器２３に注入する入力信号の元になる所定周波数の交流信号を発生する。入力信号生成器２２は、例えばデジタル／アナログ変換器または入力信号の波形（電圧）を生成するハードウェアである。発振器２３は、同期引込を行わせるための通常用
いられる発振器であり、例えばＥ級発振器などを含む。
高周波変換器３０は、発振器２３から供給される発信信号の周波数を無線通信用の高周波数に変換し、変換した高周波信号をアンテナ３１に供給するためのハードウェアである。
【００４４】
同期引込周波数帯域演算装置１０は、バス１１に接続されるＲＡＭ１２、ＲＯＭ１３、演算部１５、制御部１６、入出力インターフェース（Ｉ／Ｏ）１７ａ、１７ｂ、入力部１８および表示部１９より構成される。また、バス１１には不揮発性記憶装置（データベース）１４が接続されている。このデータベース１４には、発振器の情報５０が書き込まれる。具体的には、発振器の情報５０として、演算部１５で算出される入力強度限界εlinが書き込まれる。また、データベース１４に書き込まれる発振器の情報（テーブル）５０
として、位相感受関数Ｚ（θ）や、自然角周波数ωがある。なお、ＲＯＭ１３の容量が十分にある場合には、ＲＯＭ１３をデータベース１４として、ＲＯＭ１３に発振器の情報５０を記憶してもよい。
【００４５】
演算部１５は、主に、ＰＬ法によるロックレンジの算出を行うとともに、線形領域算出の た め の 直 接 シ ミ ュ レ ー シ ョ ン を 行 っ て い る 。 制 御 部 １ ６ は 、 例 え ば Ｃ Ｐ Ｕ （ Central Pr
ocessing Unit） 等 の 演 算 制 御 装 置 で 構 成 さ れ 、 同 期 引 込 周 波 数 帯 域 演 算 装 置 １ ０ を 構 成する各部の制御を行う。
同期引込周波数帯域演算装置１０の演算部１５で演算した結果は、入出力インターフェ
ース１７ａを介して注入同期系２０の入力信号生成器２２に供給される。また、注入同期
系２０の発振器２３の出力信号は、入出力インターフェース１７ｂを介して同期引込周波
数帯域演算装置１０に供給され、演算部１５におけるロックレンジの算出および線形領域算出に利用される。
【００４６】
また、入力部１８からは、発振器の自然角周波数ω、発振器の位相感受関数Ｚ（θ）（
以 上 は デ ー タ ベ ー ス に 情 報 が な い 場 合 ） 、 所 望 の 外 部 入 力 信 号 波 形 u(τ )（ 例 え ば sinτ ）、初期の｛ε1＜ε2＜ε３，・・・＜εn｝および所望のロックレンジの精度Δ＊が入力され、演算部１５の演算に利用される。表示部１９には、演算部１５で算出されたロック
レンジのデータ（ε、Δωleft、Δωright）が表示される。
【００４７】
次に、図２〜図９を参照して、本実施の形態例における注入同期引込周波数帯域を高精度に計算する方法について説明する。
【００４８】＜ロックレンジの取得処理＞
最初に図２を参照して、ロックレンジ（ＬＲ）の取得処理（メイン）について説明する。図２に示す処理は、ロックレンジの左端Δωleftと右端Δωrightの一方についてのみ示しているが、この処理はロックレンジの左端Δωleftと右端Δωrightの両側で行われる処理である。なお、図面では、ロックレンジはＬＲと記載する。
【００４９】
ロックレンジの取得処理が開始されると（ステップＳ１０）、所望されるロックレンジの精度Δ＊が入力される（ステップＳ１１）。このロックレンジの精度Δ＊は、あくまでも同期引込周波数帯域の計算実行を指示するユーザが所望する精度値であり、ユーザが任意に入力する値である。例えば、精度Δ＊を数％以下にしたい場合には、「数％」が所望
Δ＊の精度になる。
【００５０】

10

20

30

40

50

(9)

JP 7323911 B2 2023.8.9

続いて、発振器２３の自然角周波数ωと発振器の位相感受関数Ｚ（θ）が取得される（ステップＳ１２）。自然角周波数ωと位相感受関数Ｚ（θ）は発振器２３の固有の数値であり、ロックレンジを計算しようとしている発振器自体が具備するものである。
【００５１】
ステップＳ１１とステップＳ１２で、ユーザにより所望の精度Δ＊が入力され、発振器
２３の自然角周波数ωと位相感受関数Ｚ（θ）が取得されると、次に、入力強度限界εlinをデータベース１４から読み出す（ステップＳ１３）。ここで、入力強度限界εlinは、
リミットサイクル（ＬＣ）の乖離度ｄ（ε）が線形的に振る舞う入力強度限界である。この入力強度限界εlinの値は、先に説明したＰＬ法から得られるロックレンジの精度、例
えば数％が保証される入力強度限界εLRとほぼ等しい値となる。なお、入力強度εとは、
外 部 入 力 信 号 ｆ の 振 幅 の こ と で あ り 、 振 幅 と 波 形 を 分 離 し て 表 現 す る た め に ｆ (Ω t)を ε u(Ω t)と 表 す こ と も あ る 。 ま た 、 図 面 で は 、 リ ミ ッ ト サ イ ク ル は Ｌ Ｃ と 称 す る 。
【００５２】
リミットサイクル（ＬＣ）の乖離度ｄ（ε）とは、ロックレンジの右端または左端の外部 入 力 信 号 ｆ (Ω t)の 周 波 数 の Ｌ Ｃ と 、 外 部 信 号 が な い と き の リ ミ ッ ト サ イ ク ル と の 距 離を表す。なお、ＰＬ法で得られる同期引込周波数帯域がリミットサイクルの乖離度ｄ（ε）と関係があることは、伝統的な同期理論から説明することができる。
【００５３】
ステップＳ１３で、入力強度限界εlinをデータベース１４から読み出すことができるのは、当然のことながら、入力強度限界εlinがデータベース１４に保存されている場合に限られる。つまり、発振器２３に対して線形領域が予め分かっている場合には、入力強度限界εlinがデータベース１４に保存されているので、入力強度限界εlinをデータベー
ス１４から読み出すことが可能である。そうでない場合、つまり発振器２３に対して線形領域が予め分かっていない場合には、コストをかけても入力強度限界εlinを計算で求めなければならない。
【００５４】
そこで、次のステップでは、ステップＳ１３においてデータベース１４から入力強度限界εlinが読み出すことができたか否かが判断される（ステップＳ１４）。ステップＳ１
４で入力強度限界εlinがデータベース１４から読み出すことができたと判断された場合（ ス テ ッ プ Ｓ １ ４ の Ｙ ｅ ｓ ） に は 、 所 望 の 外 部 入 力 信 号 波 形 u(τ )が ユ ー ザ に よ っ て 入 力される（ステップＳ１５）。そして、入力した外部入力波形を振幅により正規化する（ステ ッ プ Ｓ １ ６ ） 。 こ の 正 規 化 は 入 力 信 号 波 形 u（ τ ） を u（ τ ） の 最 大 値 で 割 り 算 す る こ とによって行われる。
【００５５】
次に、入力強度限界εlinをよりも小さい外部信号の入力強度（振幅に相当）εに対して、ＰＬ法によりロックレンジを算出する（ステップＳ１７）。このステップＳ１７の処理の詳細については図９で後述する。
そして、外部入力信号の入力強度εと、発振器２３の角周波数と外部入力信号の角周波
数差Δωの組み合わせ（ε、Δω）を出力して（ステップＳ１８）、ロックレンジ取得処理（メイン）を終了する（ステップＳ１９）。
【００５６】
ステップＳ１４で、入力強度限界εlinがデータベース１４から読み出されなかったと判断された場合には（ステップＳ１４のＮｏ）、外部入力信号の入力強度ε｛ε1＜ε2＜
ε３，・・・＜εn｝の中から、ユーザが適当な入力強度εを選択する（ステップＳ２０）。このとき、選択する作業を行うユーザに、該当する技術についてのスキルがある場合には、スキルに基づいて適切な入力強度εを選択するのが好ましいが、ここで選択する入力強度εは十分小さな値であれば適切でなくても問題ない。
【００５７】
そ し て 、 所 望 の 外 部 入 力 信 号 波 形 u(τ )を サ イ ン 波 と し て （ u(τ )＝ sinτ ） 入 力 強 度 限界εlinの算出処理を行い（ステップＳ２１）、算出された入力強度限界εlinをデータベ

10

20

30

40

50

(10)

JP 7323911 B2 2023.8.9

ース１４に書き込む（ステップＳ２２）。これ以降、ステップＳ１５〜Ｓ１８までの処理が行われる。なおステップＳ２１の処理の詳細は、図３で後述する。
【００５８】＜入力強度限界εlinの算出処理＞
次に、図３のフローチャートを参照して、ステップＳ２１における入力強度限界εlinの算出処理について説明する。
入力強度限界εlinの算出処理が開始されると（ステップＳ２５）、外部入力信号の入力強度εをｎ段階に分割し、分割した入力強度ε１〜εｎ｛ε1＜ε2＜ε３，・・・＜εn｝に対してロックレンジを、ＬＣ算出処理によって求める（ステップＳ２６）。このス
テップＳ２６の処理の詳細は、図６で後述する。
【００５９】
次 に 、 ス テ ッ プ Ｓ ２ ６ の 算 出 処 理 で 得 ら れ た 入 力 強 度 ε と 出 力 電 圧 Vout(t)か ら リ ミ ッ
ト サ イ ク ル （ Ｌ Ｃ ） の 乖 離 度 ｄ (ε )を 算 出 す る （ ス テ ッ プ Ｓ ２ ７ ） 。 こ こ で 図 ４ を 参 照 して リ ミ ッ ト サ イ ク ル の 乖 離 度 ｄ (ε )に つ い て 説 明 す る と 、 横 軸 は 入 力 強 度 ε 1 、 ε 2 、 ε 3、ε4、ε5を示し、縦軸はリミットサイクルの乖離度ｄ（ε）を示す。図４に示すように、×印で示した乖離度ｄ（ε）は、入力強度ε3とε4の間に、大きなギャップ（開き）が認められる。すなわち、入力強度ε1〜ε3までは、実線で示した直線上に載っているが、
入力強度ε4、ε5は、直線上にない。この図４から、入力強度限界εlinが入力強度ε3に近いことが分かる。
なお、ステップＳ２７の処理の詳細は、図７で後述される。
【００６０】
ス テ ッ プ Ｓ ２ ７ の 処 理 の 後 、 乖 離 度 ｄ (ε ｍ )(m=１ ,２ ・ ・ ・ ｎ ） に つ い て 、 平 均 二 乗 誤
差 E(ε ｍ ） を 計 算 す る （ ス テ ッ プ Ｓ ２ ８ ） 。 こ こ で 、 平 均 二 乗 誤 差 E(ε ｍ ） に つ い て 、 図
５を参照して説明する。図４と同様に、図５の横軸も入力強度ε1、ε2、ε3、ε4、ε5であり、縦軸は平均二乗誤差Ｅ（ε）である。図５に示すように、入力強度ε3とε4の間に平均二乗誤差Ｅ（ε）に大きなギャップがみられる。このことは、入力強度ε1〜ε3までが線形領域であり、入力強度ε4、ε5は非線形領域なっていることを示しており、これによって入力強度限界εlinは、入力強度ε3に近いことが理解できる。このステップＳ２８の算出処理の詳細は、図８で後述される。
【００６１】
次 に 、 ス テ ッ プ Ｓ ２ ８ で 求 め た 平 均 二 乗 誤 差 E(ε m ） が 、 一 つ 前 の 段 階 の 平 均 二 乗 誤 差 E(ε ｍ − １ ） か ら 所 定 の 倍 率 （ 例 え ば ２ 倍 ： E(ε ｍ ） =２ E(ε ｍ − １ ） ） に な る ｍ が 存 在 す
る か 否 か が 判 断 さ れ る （ ス テ ッ プ Ｓ ２ ９ ） 。 ス テ ッ プ Ｓ ２ ９ で 、 平 均 二 乗 誤 差 E(ε m ） が
一 つ 前 の 平 均 二 乗 誤 差 E(ε ｍ − １ ） に 対 し て 所 定 倍 率 （ 例 え ば ２ 倍 ） に な っ て い る と 判 断さ れ た 場 合 （ ス テ ッ プ Ｓ ２ ９ の Ｙ ｅ ｓ ） に は 、 平 均 二 乗 誤 差 が E(ε ｍ ） に な っ た と き の 入力強度εmよりも一つ前の段階の入力強度εm-1を入力強度限界εlinとし（ステップＳ３０）、入力強度限界εlinの算出処理を終了する（ステップＳ３１）。
【００６２】
ス テ ッ プ Ｓ ２ ９ で 、 平 均 二 乗 誤 差 E(ε ｍ ） が 一 つ 前 の 平 均 二 乗 誤 差 E(ε ｍ − １ ） に 対 して所定倍率（例えば、２倍）になっていないと判断された場合（ステップＳ２９のＮｏ）には、ＰＬ法によるロックレンジ算出処理がなされる（ステップＳ３２）。なお、ステッ
プＳ３２の処理の詳細は、図９で後述する。
【００６３】
そして、ステップＳ２６のロックレンジ、リミットサイクル算出処理により得られた（
ε ,Δ ω s i m ） と 、 ス テ ッ プ Ｓ ３ ２ の Ｐ Ｌ 法 に よ る Ｌ Ｒ 算 出 処 理 で 得 ら れ た （ ε ,Δ ω P L ）に 基 づ い て 、 ロ ッ ク レ ン ジ の 精 度 Δ (ε m ) (m=１ ,２ ・ ・ ・ ｎ ） を 算 出 す る （ ス テ ッ プ Ｓ ３
３）。
次 に 、 ス テ ッ プ Ｓ ３ ３ の 算 出 処 理 の 結 果 、 ロ ッ ク レ ン ジ （ Ｌ Ｒ ） の 精 度 Δ (ε ｍ ) (m=１
,２ ・ ・ ・ ｎ ） が 所 望 の 精 度 Δ ＊ （ 例 え ば 、 １ ％ ） よ り 大 き い か 否 か が 判 断 さ れ る （ ス テ
ップＳ３４）。

10

20

30

40

50

(11)

JP 7323911 B2 2023.8.9【００６４】
ス テ ッ プ Ｓ ３ ４ で 、 ロ ッ ク レ ン ジ の 精 度 Δ (ε m ） が 所 望 の 精 度 Δ ＊ （ 例 え ば 、 １ ％ ） より大きいｍが存在する判断された場合（ステップＳ３４のＹｅｓ）には、ステップＳ３０に移行して、入力強度εｍの一つ前の段階の入力強度εｍ−１を入力強度限界εlinとする 。 ス テ ッ プ Ｓ ３ ４ で こ の よ う な 判 断 を し て い る の は 、 ロ ッ ク レ ン ジ の 精 度 Δ (ε ） の 線
形領域が広がっていたとしても、精度Δ＊（例えば、１％）以上に大きくなってしまったら、意味がなくなるからである。
【００６５】
ス テ ッ プ Ｓ ３ ４ で 、 ロ ッ ク レ ン ジ の 精 度 Δ (ε m ） が 所 望 の 精 度 Δ ＊ （ 例 え ば 、 １ ％ ） よ
り大きくなるｍが存在しないと判断された場合（ステップＳ３４のＮｏ）には、入力強度
εに入力強度εnよりも大きな入力強度となるεn+1｛ε1＜ε2＜ε３，・・・＜εn＜εn＋１｝を追加し、ｎ＋１→ｎと定義し直して、ステップＳ２６のロックレンジ、リミット
サイクル算出処理に戻る（ステップＳ３５）。以上が、入力強度限界εlinの算出処理である。
【００６６】＜ロックレンジ、リミットサイクル算出処理＞
次に、図６を参照して、図３のステップＳ２６のロックレンジ（ＬＲ）、リミットサイ
クル（ＬＣ）算出処理の詳細について説明する。
処理が開始されると（ステップＳ４０）、まず、［数１１］式に示す位相方程式から、発信器と外部入力信号の角周波数差（離調）Δωを算出する（ステップＳ４１）。
【００６７】【数１１】【００６８】
ステップＳ４１で、ロックレンジの左端を算出する場合には、ΔωがΔωleftとなり、ロックレンジの右端を算出する場合には、ΔωがΔωrightとなる。
-5
次 に 、 離 調 Δ ω を ス テ ッ プ 幅 ｈ と 置 き 、 そ の 最 小 の 値 ｈ s t o p を Δ ω × 10 に 設 定 す る （ス テ ッ プ Ｓ ４ ２ ） 。 こ の 設 定 で は 、 ス テ ッ プ 幅 ｈ を 当 初 の 値 か ら 、 1/2、 1/4、 ・ ・ ・ と 小
さくしていくようにし、ｈstopになったら、その段階でそれ以上小さくすることを止める。
【００６９】
次に、外部入力信号の角周波数ΩをΩ＝ω＋Δωに設定する（ステップＳ４３）。そして、角周波数Ωの外部入力信号を注入し、直接シミュレーションにより時間とともに変わる 出 力 電 圧 Vout(t)を 測 定 し て メ モ リ に 記 憶 す る （ ス テ ッ プ Ｓ ４ ４ ） 。
ス テ ッ プ Ｓ ４ ４ で 時 間 と と も に 変 わ る 出 力 電 圧 Vout(t)を 測 定 し た 後 に 、 周 波 数 カ ウ ン
ターを用いて周波数の時間変化を測定し、これをメモリに記憶する（ステップＳ４５）。
そして、平均角周波数ω〜（オメガチルダー）を計算する（ステップＳ４６）。次に、計算した平均角周波数ω〜と外部信号の角周波数Ωとを比較して、同期が成立しているか否かを判断する（ステップＳ４７）。ここで、計算した平均角周波数ω〜が外部信号の角周
波数Ωに限りなく近くなったとき、つまりω〜／Ωが限りなく「１」になった時に同期が成立したと判断する。
【００７０】
ステップＳ４７で、同期が成立したと判断された場合（ステップＳ４７のＹｅｓ）には、発信器と外部入力信号の角周波数差（離調）Δω（ステップ幅ｈ）がその最小の値ｈst
opより小さいか否かが判断される（ステップＳ４８）。ステップＳ４８で角周波数差（離調）Δωに相当するステップ幅ｈが最小値ｈstopより小さいと判断された場合（ステップＳ ４ ８ の Ｙ ｅ ｓ ） に は 、 角 周 波 数 差 （ 離 調 ） Δ ω と そ の 時 の 出 力 電 圧 Vout(t)を メ モ リ に格納し（ステップＳ４９）、処理を終了して（ステップＳ５０）、図３のステップＳ２７の処理に移行する。

10

20

30

40

50

(12)

JP 7323911 B2 2023.8.9【００７１】
ステップＳ４７で、平均角周波数ω〜と外部信号の角周波数Ωの比が「１」に近い値にならず、同期が成立したと判断されなかった場合（ステップＳ４７のＮｏ）には、ステッ
プ幅ｈをｈ／２とすることにより、発信器と外部入力信号の角周波数差（離調）Δωを小さくして（ステップＳ５１）、再度ステップＳ４３からの処理を行う。
また、ステップＳ４８で、角周波数差（離調）Δωに相当するステップ幅ｈが最小値ｈstopより小さくなっていないと判断された場合（ステップＳ４８のＮｏ）には、角周波数差（離調）Δωを大きくして、再度ステップＳ４３からの処理を行う。以上が、ロックレ
ンジ、リミットサイクル算出処理である。
【００７２】＜ロックレンジ乖離度ｄ（ε）の算出処理＞
次に、図７を参照して、図３のステップＳ２７におけるリミットサイクル（ＬＣ）乖離
度ｄ（ε）の算出処理について説明する。処理が開始されると（ステップＳ５５）、外部(0)
入 力 信 号 無 し の 場 合 の 出 力 電 圧 の 時 間 変 化 V O U T (t)と 、 外 部 入 力 信 号 有 り の 場 合 の 出 力
(ε)
(0)
(ε)
電 圧 の 時 間 変 化 VO U T(t)を 求 め る 。 そ し て 、 V O U T (t)を F(θ )、 V O U T(t)を G(θ )と 置 く (ス テ ッ プ Ｓ ５ ６ )。 な お 、 θ は 外 部 入 力 信 号 の 各 周 波 数 Ω に 時 間 tを 掛 け た Ω tを 置
き 換 え た も の 、 す な わ ち θ ＝ Ω tで あ る 。
【００７３】
ここで、リミットサイクル（ＬＣ）とは、直接シミュレーションにより得られる電圧等の時間変化から作られる、相空間（状態空間）上の閉軌道をいう。例えば、回路シミュレ
ー タ 等 に よ り 、 数 値 的 に 得 ら れ た 複 数 の 出 力 電 圧 VO U T 1 、 VO U T 2 の 時 間 変 化 を 使 っ て リ ミ ッ
ト サ イ ク ル を 得 る こ と が で き る 。 ま た 、 得 ら れ る 出 力 電 圧 VO U T が １ つ だ け の 場 合 は 、 時 間の 遅 れ を τ と し て 、 （ V O U T 1 (t)、 V O U T 2 (t)） ＝ （ V O U T (t)、 V O U T (t+τ )） の 時 間 変 化 を 使ってリミットサイクルを便宜的に得ることができる。
【００７４】
次 に 、 位 相 φ の 初 期 値 ｉ を 「 ０ 」 に 設 定 す る (ス テ ッ プ Ｓ ５ ７ )。 こ の 初 期 値 「 ｉ ＝ ０ 」
は、０〜２πまでの位相φをＮ等分したときの０番目の位相φの番号である。そして、位相φ＝ｉ／Ｎ＊２πとして（ステップＳ５８）、以後の処理を行う。
最 初 に 、 外 部 入 力 信 号 無 し の 場 合 の 出 力 電 圧 の 時 間 変 化 F(θ )か ら 位 相 を φ だ け ず ら した 、 外 部 入 力 信 号 有 り の 場 合 の 出 力 電 圧 の 時 間 変 化 G(θ +φ )の 差 の 絶 対 値 を 、 相 空 間 （ 状態空間）上の閉軌道で一周させて（０〜２π）積分した値ｄp（φ）を求める（ステップＳ５９）。
【００７５】
そして、位相φの０から２πまでをＮ等分して、ステップＳ５８、Ｓ５９の計算を繰り返す（ステップＳ６０）。すなわち、ステップＳ６０で「ｉ＜Ｎ」であれば（ステップＳ
６０のＹｅｓ）、ステップＳ５８、Ｓ５９を繰り返し、ステップＳ６０で「ｉ≧Ｎ」になったと判断された場合（ステップＳ６０のＮｏ）には、ステップＳ５９で計算したｄｐ（φ ） （ ０ ˜２ π ） の 中 の 最 小 値 を 求 め て ｄ （ ＝ min d p (φ )） と す る （ ス テ ッ プ Ｓ ６ １ ） 。
そして、その時の入力強度εとｄのセット（ε、ｄ）をメモリに格納して（ステップＳ６２）、処理を終了する（ステップＳ６３）。この（ε、ｄ）がリミットサイクルで求めた乖離度ｄ（ε）に相当する。
【００７６】＜ 平 均 二 乗 誤 差 Ｅ (ε m )の 算 出 ＞
次 に 、 図 ８ を 参 照 し て 、 図 ３ の ス テ ッ プ Ｓ ２ ８ に お け る 平 均 二 乗 誤 差 Ｅ (ε m )の 算 出 処
理について説明する。処理を開始して（ステップＳ６５）、初期値ｍを１にセットする（ステップＳ６６）。このｍ＝１は、入力強度εと乖離度ｄ（ε）を｛ε1＜ε2＜・・・＜
εm｝のｍ個とったときの第１番目のデータを意味している。
【００７７】
そ し て 、 乖 離 度 の デ ー タ ｛ （ ε 1 、 ｄ (ε １ )） 、 （ ε ２ 、 ｄ (ε ２ )） 、 ・ ・ （ ε ｍ 、 ｄ (ε ｍ )） ｝ を 求 め 、 原 点 （ ０ ,０ ） を 通 る と 仮 定 し て 線 形 フ ィ ッ テ ィ ン グ を 行 っ て 、 １ 次 関

10

20

30

40

50

(13)

JP 7323911 B2 2023.8.9

数 dl i n （ ε ） を 求 め る （ ス テ ッ プ Ｓ ６ ７ ） 。 こ の ス テ ッ プ Ｓ ６ ７ の 線 形 フ ィ ッ テ ィ ン グ を図に示したものが、既に説明した図４である。
【００７８】
次 に 、 d（ ε ｉ ） と ｄ l i n （ ε ｉ ） の 差 分 を と り 、 ｉ =１ 〜 ｍ （ ｍ 個 ） の 二 乗 平 均 値 E(ε m)を 求 め る （ ス テ ッ プ Ｓ ６ ８ ） 。 そ し て 、 ス テ ッ プ Ｓ ６ ８ の 計 算 を ｍ ＝ ｎ に な る ま で 繰 り
返 す (ス テ ッ プ Ｓ ６ ９ )。 す な わ ち 、 ス テ ッ プ Ｓ ６ ９ で 、 ｍ ≦ ｎ で あ れ ば （ ス テ ッ プ Ｓ ６ ９のＮｏ）、ステップＳ６７に戻り、ｍ＞ｎになったときに、そのときの平均二乗誤差Ｅ（
ε m ） を メ モ リ に 格 納 し て （ ス テ ッ プ Ｓ ７ ０ ） 、 処 理 を 終 了 す る (ス テ ッ プ Ｓ ７ １ )。 以 上が、平均二乗誤差Ｅ（εm）の算出処理である。
【００７９】＜ＰＬ法ロックレンジ算出処理＞
次に、図９を参照して、図２のステップＳ１７および図３のステップＳ３２のＰＬ法ロ
ックレンジ算出処理について説明する。ＰＬ法ロックレンジ算出処理は、［数１２］式で
示 す Winfreeモ デ ル の 微 分 方 程 式 の 解 を 求 め る 方 法 で あ る 。 つ ま り 、 Winfreeモ デ ル の 微 分
方 程 式 を 解 く た め に 、 Poincare-Lindstedt法 に 基 づ き 、 発 信 器 の 角 周 波 数 Ω と 位 相 θ を 入力強度εのべき乗の展開式（［数１３］式）で表す方法である。ＰＬ法ＬＲ算出処理では、まず［数１３］式で示した角周波数Ωと位相θを入力強度εのべき乗の式において、入2
力強度εと入力強度εの２乗（ε ）の係数を求める。
【００８０】
【数１２】

10

20

【００８１】【数１３】【００８２】
処理が開始されると（ステップＳ７５）、まず、位相φの初期値ｉが「０」に設定される (ス テ ッ プ Ｓ ７ ６ )。 こ の 初 期 値 「 ｉ ＝ ０ 」 は 、 ０ 〜 ２ π ま で の 位 相 φ を Ｎ 等 分 し た と きの０番目の位相φのことである。そして、位相φ＝ｉ／Ｎ＊２πとして（ステップＳ７７）、以後の処理を行う。
【００８３】
既に説明したように、発信器の角周波数Ωを入力強度εのべき乗の展開式（［数１３］
(0)
(1)
(2)
式 ） で 表 し た と き の 入 力 強 度 ε の 係 数 を 順 に Ω (φ )、 Ω (φ )、 Ω (φ )、 ・ ・ ・(0)(1)と す る と 、 ０ 次 の 係 数 は Ω (φ )＝ ω に な り 、 １ 次 の 係 数 Ω (φ )は 、 ［ 数 １ ４ ］ 式 によって求めることができる（ステップＳ７８）。
【００８４】

30

【数１４】

40

【００８５】ま た 、 位 相 θ (τ )も ［ 数 １ ３ ］ 式 に 示 す よ う に 、 入 力 強 度 ε の べ き 乗 の 展 開 式 で 表 す こ(0)と が で き る が 、 そ の と き の ０ 次 の 係 数 は θ (τ )＝ τ と な る 。 な お 、 τ は 外 部 入 力 信 号の 角 周 波 数 Ω に 時 間 ｔ を 掛 け た Ω tを 置 き 換 え た も の 、 す な わ ち τ ＝ Ω tで あ る 。
【００８６】
位 相 θ (τ )を 入 力 強 度 ε の べ き 乗 で 展 開 し た ［ 数 １ ３ ］ 式 に お け る ε の 係 数 は 、 ［ 数 １５］式で求められる（ステップＳ７９）。
【００８７】50

(14)JP 7323911 B2 2023.8.9【数１５】【００８８】
２(2)
さ ら に 、 発 信 器 の 角 周 波 数 Ω (φ )の 入 力 強 度 ε の ２ 乗 （ ε ） の 係 数 Ω (φ )は 、 ［数１６］式で求めることができる（ステップＳ８０）。
【００８９】【数１６】【００９０】
以上のステップＳ７８〜Ｓ８０の計算がｉ＝０からｉ＝Ｎまで繰り返され、位相φが０〜２πの間でＮ分割されたそれぞれの値に対して行われる（ステップＳ８１）。すなわち、ステップＳ８１で、ｉ＜Ｎ（φ＜２π）であれば（ステップＳ８１のＹｅｓ）、ｉがＮ（φ＝２π）になるまで、ステップＳ７７〜Ｓ８０の処理が繰り返される。そして、ステップＳ８１でｉ≧Ｎになったと判断された場合に（ステップＳ８１のＮｏ）、次の入力強度限界εlinを求める処理に移行する。
【００９１】
こ こ で も 、 入 力 強 度 ε を 求 め る た め の 「 ｊ 」 が 「 ０ 」 に 設 定 さ れ る (ス テ ッ プ Ｓ ８ ２ )。
ここで「ｊ＝０」は、０から入力強度限界εlinまでの入力強度εをＭ段階で分割したときの初期値を示している。また、入力強度εjは、入力強度限界εlinを使って、εj＝ｊ／Ｍ＊εlinで表される（ステップＳ８３）。この式は、入力強度限界εlinをＭ分割した
ときのｊ（ｊ＝０〜Ｍ）番目の入力強度εjを示す式である。
【００９２】(1)(2)
そ し て 、 ス テ ッ プ Ｓ ７ ８ で 求 め た Ω (φ )と 、 ス テ ッ プ Ｓ ８ ０ で 求 め た Ω (φ )を 使って、外部入力信号Ω（φ）を近似式で求める（ステップＳ８４）。このステップＳ８４の式は、発信器と外部入力信号の角周波数差（離調）Δω（＝Ω（φ）−ω）が入力強度(1)(2)
ε の ２ 次 式 で 表 さ れ る こ と を 意 味 し て お り 、 １ 次 係 数 が Ω (φ )、 ２ 次 係 数 が Ω (φ )である。
【００９３】
次 に 、 φ （ ０ 〜 ２ π ） 間 の 最 小 値 の 外 部 信 号 の 角 周 波 数 Ω （ φ ） （ minΩ （ φ ） ） と 発振器の角周波数ωとの差を求め、ロックレンジの左端Δωleftを算出する。さらに、φ（
０ 〜 ２ π ） 間 の 最 大 値 の 外 部 信 号 の 角 周 波 数 Ω （ φ ） （ maxΩ （ φ ） ） と 発 振 器 の 角 周 波
数ωとの差を求め、ロックレンジの右端Δωrightを算出する（ステップＳ８５）。
【００９４】
以上のステップＳ８３〜Ｓ８５の計算をｊ＝０からｊ＝Ｍまで繰り返し、入力強度εの計算を０から入力強度限界εlinの間でＭ分割されたそれぞれの値に対して行う（ステップＳ８６）。すなわち、ステップＳ８６で、ｊ＜Ｍ（εj＜εlin）であれば（ステップＳ
８８のＹｅｓ）、ｊがＭ（εj＝εlin）になるまで、ステップＳ８３〜Ｓ８５の処理が繰
り 返 さ れ る 。 そ し て 、 ス テ ッ プ Ｓ ８ ６ で j≧ Ｍ に な っ た と 判 断 さ れ た 場 合 に は （ ス テ ッ プ
Ｓ８８のＮｏ）、（ε、Δω）を記憶部に格納して（ステップＳ８７）処理を終了する（ステップＳ８８）。そして、図２のステップＳ１８に移行する。
【００９５】＜発振器の構成例＞
図１０は、図１に示す発振器２３の具体的な構成例を示す図である。ここでは、発振器
２３として、リングオシレータとした例である。
信号源（パルスジェネレータ）２１からのパルスが入力信号生成器（バイアスサーキッ
ト）２２に供給され、入力信号生成器２２で生成された信号が、発振器（リングオシレー
タ）２３の入力端子２３ａに供給される。
発振器２３は、３段に直列接続されたインバータ２３ｂ、２３ｃ、２３ｄとスイッチン

20

30

40

50

(15)

JP 7323911 B2 2023.8.9

グ素子２３ｆとを備え、３段目のインバータ２３ｄの出力信号が、１段目のインバータ２３ｂの入力に戻されると共に、発振器２３の出力として出力端子２３ｅから出力される。
【００９６】
入力端子２３ａに得られる信号は、スイッチング素子２３ｆのゲートに供給される。ス
イッチング素子２３ｆのドレインは、２段目のインバータ２３ｃと３段目のインバータ２３ｄの間に接続され、スイッチング素子２３ｆのソースは、１段目のインバータ２３ｂの入力側に接続される。
本実施の形態例では、このように注入同期発振器として構成された発振器２３の同期引込周波数帯域を、高精度に算出することができる。
【００９７】＜同期引込周波数帯域演算精度の例＞
図１１は、本実施の形態例による同期引込周波数帯域演算装置１０で、図１０に示す発振器（リングオシレータ）２３の同期引込周波数帯域を演算した場合を示す。
この図１１において、従来方法と記載された特性は、位相方程式で求めた特性であり、本発明と記載された特性は、本実施の形態例による同期引込周波数帯域演算装置１０で求めた特性である。
【００９８】
図１１の（Ａ）は、ロックレンジの左端におけるリミットサイクル（ＬＣ）の乖離度ｄleft（ε）を示し、図１１の（Ｂ）は、ロックレンジの右端におけるリミットサイクル（
ＬＣ）の乖離度ｄright（ε）を示す。
図１１の（Ｃ）は、ロックレンジ（ＬＲ）の左端における精度Δleft（ε）を示し、図１１の（Ｄ）は、ロックレンジ（ＬＲ）の右端における精度Δright（ε）を示す。図１
１の（Ｃ），（Ｄ）に示す△印の特性は、位相方程式（従来方法）で求めた精度であり、○印の特性は、本実施の形態例により求めた精度である。これらの精度の数値はパーセン
トである。
図１１において、横軸は外部入力信号の入力強度（電圧振幅ε［Ｖ］）である。
【００９９】
なお、図１１の（Ａ）および（Ｃ）は、外部入力信号の角周波数を、入力強度毎に注入同期が成立する最小のものにした場合であり、図１１の（Ｂ）および（Ｄ）は、外部入力
信号の角周波数を、入力強度毎に注入同期が成立する最大のものにした場合である。
ここでは、図１１の（Ａ）および（Ｂ）に示すデータ点（□印）がリミットサイクル（
ＬＣ）の乖離度であり、図１１の（Ｃ）および（Ｄ）に示すデータ点（△印）が従来方法（位相方程式）により得られたロックレンジ（ＬＲ）の精度、データ点（○印）が本実施の形態例（本発明）により得られたロックレンジ（ＬＲ）の精度である。
【０１００】
図１１の（Ａ）および（Ｂ）から判るように、入力強度ε＝０．０９０［Ｖ］およびε＝０．１０５［Ｖ］が、ロックレンジ（ＬＲ）の左端と右端における、それぞれのリミッ
トサイクル（ＬＣ）の乖離度が線形的に振る舞う入力強度限界εlinである。本実施の形
態例の場合、図１１の（Ｃ）および（Ｄ）に示すように、ロックレンジ（ＬＲ）の精度が従来方法（位相方程式）よりも優れており、入力強度がε≦εlinであれば、予め設定した精度１％よりも良くなることがわかる。
【０１０１】
次に、同期引込周波数帯域の算出にかかる時間（ＣＰＵ時間）を、直接シミュレーショ
ン法と位相方程式（従来手法）と本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合で比較する。
直 接 シ ミ ュ レ ー シ ョ ン 法 で は 、 Synopsys社 の 回 路 シ ミ ュ レ ー タ HSPICEを イ ン テ ル 社 製 の2.6GHz× 12コ ア Xeonプ ロ セ ッ サ E5-2630v2（ た だ し 、 デ ル 社 の PRECISION T5610を 使 用 ） で実行した。
位相方程式（従来手法）および本実施の形態例における同期引込周波数帯域演算装置１０ （ 本 発 明 ） で は 、 C言 語 で 実 装 し た プ ロ グ ラ ム を イ ン テ ル 社 製 の 3.6GHz× 8コ ア Xeonプ ロ

10

20

30

40

50

(16)

JP 7323911 B2 2023.8.9

セ ッ サ Core i7-4790（ た だ し 、 デ ル 社 の OPTIPLEX9020を 使 用 ） で 実 行 し た 。
【０１０２】
同期引込周波数帯域の算出にかかる時間（ＣＰＵ時間）は、
直接シミュレーション法の場合は約６時間であり、位相方程式（従来手法）の場合は約１μ秒であり、本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合は約２〜３μ秒である。
【０１０３】
本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合の計算時間
は、直接シミュレーション法の場合よりも非常に短く、従来手法である位相方程式の場合と比べても同等程度であることがわかる。
【０１０４】
このように本実施の形態例の同期引込周波数帯域演算装置１０によると、従来よりも非常に少ない計算コストで、非常に精度の高い同期引込周波数帯域を得ることができるようになる。すなわち、位相方程式により同期引込周波数帯域を得る場合よりも高い精度で算出ができると共に、直接シミュレーション法よりも大幅に少ない計算で行える。
【０１０５】＜本発明の第２の実施の形態例＞
次に、本発明の第２の実施の形態例を、図１２〜図１４を参照して説明する。
本発明の第２の実施の形態例では、生体の電気活動を算出する同期引込周波数帯域演算装置としたものである。同期引込周波数帯域演算装置としての構成や演算処理は、第１の実施の形態例で説明した同期引込周波数帯域演算装置１０と同じである。
【０１０６】
図 １ ２ は 、 生 体 の 神 経 細 胞 の 電 気 活 動 を 表 す 代 表 的 モ デ ル で あ る 、 Hodgkin-Huxleyモ デルを示す。図１２の左側に示すように、細胞体１０１は、核１０２、樹状突起１０３、シナプス１０４、軸索丘１０５、ミエリン鞘１０６、および軸索抹消１０７を備える。
ここで、図１２の右側に示すように、神経細胞の電気活動が伝導する箇所である神経線
維１１０に、電位計１４０が接続された微細電極１２０および電極１３０を装着して、生体の神経細胞の電気活動による電位を計測する。つまり、生体が備える神経細胞を発振器と し て 、 こ の Hodgkin-Huxleyモ デ ル に よ る 同 期 引 込 周 波 数 帯 域 を 、 同 期 引 込 周 波 数 帯 域 演
算装置１０が演算する。
【０１０７】
こ の Hodgkin-Huxleyモ デ ル は 、 図 １ ２ に 示 す よ う な 神 経 細 胞 に 見 ら れ る 活 動 電 位 に つ いて 、 現 象 論 的 に 記 述 し た ｖ ， ｍ ， ｈ ， ｎ の 4つ の 変 数 か ら な る 連 立 方 程 式 （ ［ 数 １ ７ ］ 式）である。
【０１０８】【数１７】【０１０９】［数１７］式の連立方程式において、各変数の意味は以下の通りである。
ｖ：細胞膜の電位（膜電位）＋ｍ：Ｎａ チャネルの活性化ゲートが開いている割合＋
ｈ：Ｎａ チャネルの不活性化ゲートが開いている割合＋
ｎ：Ｋ チャネルの活性化ゲートが開いている割合【０１１０】＋
図１３は、活性化ゲートや不活性化ゲートに生じるナトリウムイオン（Ｎａ ）やカリ

10

20

30

50

(17)

JP 7323911 B2 2023.8.9＋

ウムイオン（Ｋ ）を示すものである。
【０１１１】
連立方程式でのｍ，ｈ，ｎの時間変化は、それぞれのゲートが開く速度定数αと閉じる
速度定数βからなる方程式に従う。速度定数α，βは膜電位ｖに依存し、その依存性は、
イカ巨大軸索を用いた実験値により便宜的に以下の［数１８］式で表される。
【０１１２】【数１８】【０１１３】
Hodgkin-Huxleyモ デ ル で の 、 そ の 他 の 定 数 は 以 下 の 通 り で あ る 。
Ｃｍ：細胞膜のコンデンサーとしての容量
Ｉｄ：細胞内に注入する電流の直流成分
Ｉ（ωｔ）：細胞内に注入する電流の交流成分
また、各チャネルを流れる電流は、以下の式で定義される。
【０１１４】【数１９】【０１１５】＋
ＩＮａ：Ｎａ チャネル全体を流れる電流（３つの活性化ゲートと１つの不活性ゲー
トによって開閉される）＋
ｇＮａ：Ｎａ のイオン電流の最大コンダクタンス＋
ｖＮａ：Ｎａ のイオン電流がゼロになる平衡電位＋
ＩＫ：Ｋ チャネル全体を流れる電流（４つの活性化ゲートによって開閉される）＋
ｇＫ：Ｋ のイオン電流の最大コンダクタンス＋
ｖＫ：Ｋ のイオン電流がゼロになる平衡電位
ＩＬ：その他のイオン電流（リーク電流と呼ばれる）
ｇＬ：リーク電流の最大コンダクタンス
ｖＬ：リーク電流がゼロになる平衡電位【０１１６】
こ の よ う な Hodgkin-Huxleyモ デ ル を 適 用 し て 、 第 １ の 実 施 の 形 態 例 で 説 明 し た 処 理 で 、同期引込周波数帯域演算装置１０は、生体の同期引込周波数帯域を算出することができる。
図１４は、本実施の形態例による同期引込周波数帯域演算装置１０で、図１２に示す神経細胞の同期引込周波数帯域を演算した場合を示す。
この図１４において、従来方法と記載された特性は、位相方程式で求めた特性であり、本発明と記載された特性は、本実施の形態例で求めた特性である。
【０１１７】
図１４の（Ａ）は、ロックレンジの左端におけるリミットサイクル（ＬＣ）の乖離度ｄleft（ε）を示し、図１４の（Ｂ）は、ロックレンジの右端におけるリミットサイクル（
ＬＣ）の乖離度ｄright（ε）を示す。
図１４の（Ｃ）は、ロックレンジ（ＬＲ）の左端における精度Δleft（ε）を示し、図１１の（Ｄ）は、ロックレンジ（ＬＲ）の右端における精度Δright（ε）を示す。図１
４の（Ｃ），（Ｄ）に示す△印の特性は、位相方程式（従来方法）で求めた精度であり、○印の特性は、本実施の形態例により求めた精度である。これらの精度の数値はパーセン

30

40

50

(18)

JP 7323911 B2 2023.8.9

トである。
図１４において、横軸は外部入力信号の入力強度（電流振幅ε［Ａ］）である。
【０１１８】
図１４の（Ａ）および（Ｃ）は、外部入力信号の角周波数を、入力強度毎に注入同期が成立する最小のものにした場合であり、図１４の（Ｂ）および（Ｄ）は、外部入力信号の角周波数を、入力強度毎に注入同期が成立する最大のものにした場合である。
ここでは、図１４の（Ａ）および（Ｂ）に示すデータ点（□印）がリミットサイクル（
ＬＣ）の乖離度であり、図１４の（Ｃ）および（Ｄ）に示すデータ点（△印）が従来方法（位相方程式）により得られたロックレンジ（ＬＲ）の精度、データ点（○印）が本実施の形態例（本発明）により得られたロックレンジ（ＬＲ）の精度である。
【０１１９】
図１４の（Ａ）および（Ｂ）から判るように、入力強度ε＝１．０［Ａ］が、ロックレンジ（ＬＲ）の左端と右端における、それぞれのリミットサイクル（ＬＣ）の乖離度が線形的に振る舞う入力強度限界εlinである。本実施の形態例の場合、図１４の（Ｃ）および（Ｄ）に示すように、ロックレンジ（ＬＲ）の精度が従来方法（位相方程式）よりも優れており、入力強度がε≦εlinであれば、予め設定した精度１％よりも良くなることがわかる。
【０１２０】
次に、同期引込周波数帯域の算出にかかる時間（ＣＰＵ時間）を、直接シミュレーショ
ン法と位相方程式（従来手法）と本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合で比較する。
直接シミュレーション法、位相方程式（従来手法）および本実施の形態例における同期
引 込 周 波 数 帯 域 演 算 装 置 １ ０ （ 本 発 明 ） で は 、 い ず れ も C言 語 で 実 装 し た プ ロ グ ラ ム を イ
ン テ ル 社 製 の 3.6GHz× 8コ ア Xeonプ ロ セ ッ サ Core i7-4790（ た だ し 、 デ ル 社 の OPTIPLEX902
0を 使 用 ） で 実 行 し た 。
【０１２１】
同期引込周波数帯域の算出にかかる時間（ＣＰＵ時間）は、
直接シミュレーション法の場合は約２時間であり、位相方程式（従来手法）の場合は約１μ秒であり、本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合は約２〜３μ秒である。
【０１２２】
本実施の形態例における同期引込周波数帯域演算装置１０（本発明）の場合の計算時間
は、直接シミュレーション法の場合よりも非常に短く、従来手法である位相方程式の場合と比べても同等程度であることがわかる。
【０１２３】
このように本実施の形態例の同期引込周波数帯域演算装置１０においても、従来よりも非常に少ない計算コストで、非常に精度の高い同期引込周波数帯域を得ることができるようになる。すなわち、位相方程式により同期引込周波数帯域を得る場合よりも高い精度で算出ができると共に、直接シミュレーション法よりも大幅に少ない計算で行える。
【０１２４】＜変形例＞
なお、上述した各実施の形態例は、本発明の好適な例を示すものであり、本発明の同期
引込周波数帯域演算装置は、様々な発振器を備える注入同期系の同期引込周波数帯域の演算に用いることができる。また、同期引込周波数帯域を演算する際に、上述した実施の形
態例では、求めた乖離度の平均二乗誤差の変化から乖離度の線形領域となる入力強度限界を取得する処理まで、演算で求めるようにした。これに対して、線形領域か否かの判断は、演算結果の表示出力結果からユーザが判断し、その判断した結果を同期引込周波数帯域
演算装置に入力させるようにしてもよい。
【０１２５】

10

20

30

40

50

(19)

JP 7323911 B2 2023.8.9

また、図１に示す同期引込周波数帯域演算装置１０は、専用の演算装置で構成する他に、本発明の同期引込周波数帯域演算処理の各手順を実行するプログラムを作成して、そのプログラムのコンピュータ装置での実行で、同期引込周波数帯域演算装置として機能するようにしてもよい。
【符号の説明】【０１２６】
１０…同期引込周波数帯域演算装置、１１…バス、１２…ＲＡＭ、１３…ＲＯＭ、１４…データベース部、１５…演算部、１６…制御部、１７ａ，１７ｂ…入出力インターフェース（Ｉ／Ｏ）、１８…入力部、１９…表示部、２０…注入同期系、２１…信号源、２２…入力信号生成器、２３…発振器、３０…高周波変換器、３１…アンテナ、５０…テーブ
ル、１０１…細胞体、１０２…核、１０３…樹状突起、１０４…シナプス、１０５…軸索丘、１０６…ミエリン鞘、１０７…軸索抹消、１１０…神経線維、１２０…微細電極、１３０…電極、１４０…電位計

【図１】

【図２】

10

(20)

【図３】

【図４】

【図５】

【図６】

【図７】

JP 7323911 B2 2023.8.9

(21)

【図８】

【図９】

【図１０】

【図１１】

JP 7323911 B2 2023.8.9

(22)

【図１２】

【図１４】

【図１３】

JP 7323911 B2 2023.8.9

(23)JP 7323911 B2 2023.8.9

̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶̶
フロントページの続き(72)発明者 中川正基

東京都調布市調布ケ丘一丁目５番地１(72)発明者 慶田

東京都調布市調布ケ丘一丁目５番地１
審査官

及川

国立大学法人電気通信大学内朗
国立大学法人電気通信大学内尚人(56)参考文献 特開２０１５−１４６５３４（ＪＰ，Ａ）
特開２０１０−１４７５９９（ＪＰ，Ａ）
特開２００９−１１７８９４（ＪＰ，Ａ）
特開平０８−１６１４３４（ＪＰ，Ａ）(58)調査した分野(Int.Cl.，ＤＢ名)
Ｈ０３Ｌ

１／００−９／００