(57)【特許請求の範囲】
【請求項１】
２つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記２つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記２つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする位相差計測方法。
【請求項２】
２つの被計測信号をそれぞれＡ／Ｄ変換するデジタル変換ステップと、デジタル変換後の２つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記２つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記２つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする位相差計測方法。
【請求項３】
被計測信号をフィルタリングし、このフィルタリング前後での信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後の２つの信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記２つのデジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする周波数計測方法。
【請求項４】
被計測信号をＡ／Ｄ変換するデジタル変換ステップと、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記デジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする周波数計測方法。
【請求項５】
２つの被計測信号をそれぞれＡ／Ｄ変換するＡ／Ｄ変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の２つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出し、前記２つの被計測信号それぞれの振幅値から前記面積を正規化し、正規化された前記面積から前記２つの被計測信号の位相差を算出することを特徴とする位相差計測システム。
【請求項６】
被計測信号をＡ／Ｄ変換するＡ／Ｄ変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出し、前記被計測信号の振幅値から前記面積を正規化し、正規化された前記面積から前記デジタル信号の位相差を算出し、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定することを特徴とする周波数計測システム。
【請求項７】
前記フィルタリングに、IIRくし型フィルタを用いることを特徴とする請求項６記載の周波数計測システム。
【請求項８】
前記フィルタリングに、線形位相FIRフィルタを用いることを特徴とする請求項６記載の周波数計測システム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、リサージュ図となる閉曲線から被計測信号の位相差を求める位相差計測方法と位相差計測システム、並びに、リサージュ図となる閉曲線から被計測信号の周波数を求める周波数計測方法及び周波数計測システムに関する。
【背景技術】
【０００２】
信号処理の分野において基準周波数との位相差（位相角）を計測することは重要なコア技術である。例えば、電力・電力量の位相差（位相角）測定、設備や機器の振動分析による異常診断や、静音性の評価や騒音原因の対策方法を検討するための騒音分析など様々な分野、更には音響処理における音源推定や楽音処理における倍音の分析において位相差の計測は重要である。
【０００３】
従来、位相差計測方法としては、リサージュ法(Lissajous method)によりリサージュ図を描画する方法や、ＦＦＴ（Fast Fourier Transform）アナライザを用いた方法や、波形から周期を測定する方法（電力・電力量の測定において使用されている）や、波形の部分間の相関を利用する方法（音声処理系で使用されている）などが知られている。
【０００４】
既知の位相差計測方法のうち、周波数校正においては、微小な周波数偏差を検出するための簡便な方法としてリサージュ法が用いられることがある。
リサージュ図は、計測表示装置をＸ−Ｙ入力モードに設定し、二つの単振動波形の順序対として得られる点の軌跡が描く平面図形である。周波数比が１の場合、位相差に応じて、形状は楕円・円・直線となり、周波数比が整数の場合は複数の山が描かれる。
【０００５】
リサージュ図から２つの被計測信号の位相差を求める従来の方法としては、例えば、計測表示装置のＸ軸とＹ軸に、それぞれ周波数の等しい正弦波を入力すると、楕円形のリサージュ図が現れる。このときの位相差を正確に計算するためには図形がＸ軸（又はＹ軸）を横切る点をそれぞれ正確に求めることとなる。そして、リサージュ図を描く閉曲線のＹ軸方向の最大振幅と、前記閉曲線がＹ軸を横切る際のＹ軸方向の振幅とから位相差を計算する。
【０００６】
特許文献１には、多軸加振試験機の位相差制御方法に関し、複数の系の信号を用いて表示装置上に正規化したリサージュ図を重ねて描かせることで入力信号の位相差と出力信号の位相差が揃っているか否かを検出することが記載されている。
【０００７】
特許文献２には、光干渉計の周期誤差低減方法に関し、極座標変換した直交２相信号の位相角と半径の変動値を、誤差を求める各式に代入して前記直交２相信号の振幅、零点及び位相差の補正値を求め、極座標の半径が一定となるように補正することが記載されている。そして、光干渉計の動作原理を説明するために、リサージュ図が説明されている。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開２０００−２９８１４９号公報
【特許文献２】特許第４４６５４５１号公報
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、上記従来のリサージュ法では、位相差を正確に計算するためには図形がＸ軸（又はＹ軸）を横切る点の座標をそれぞれ正確に求めなければならない。デジタル信号処理によって補間処理を行うことでゼロクロス点を求めることができるが、軌跡が原点に近づくにつれ原点近傍を通過する相対的な速度が大きくなり補間誤差が低下してしまう。
そして、従来、周波数カウンタでは「レシプロカル方式」を採用しており、小数点以下の周波数の計測が可能となっているが、周期を精密に計測するためには非常に高いサンプリング周波数を使用し周期を計測する必要がある。また周期の基準となるための「ゼロクロス点」の検出精度はＤＣオフセットの影響を受けてしまい検出精度が低下する懸念がある。
また、上記従来のリサージュ法では、微小な周波数偏差、例えば０．１／１０００Ｈｚを認識するためには３６０／１００００（度／秒）で変化する位相差を短時間に認識し周波数の微調整を繰り返す必要がある。
【００１０】
そこで、本願発明の目的は、従来のリサージュ法による処理の仕方を改良することで、精度良く位相差計測を行なうことが可能な位相差計測方法及び位相差計測システムと、精度良く周波数計測を行なうことが可能な周波数計測方法及び周波数計測システムを提供することにある。
【課題を解決するための手段】
【００１１】
本願発明者は、リサージュ図が描く閉曲線の囲む面積に着目し、短時間で精密な位相差を検出できる方法を見出し、本願発明を完成するに至った。
【００１２】
本発明の位相差計測方法は、２つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記２つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記２つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする。
【００１３】
従来のリサージュ法では、リサージュ図のうちのゼロクロスする近傍のデータしか計算に使用されていないが、本発明によれば、リサージュ図の閉曲線の囲む面積を算出することで、リサージュ図を構成する全データが計算に使用されることとなり、精密な位相差計測ができる。
【００１４】
本発明の位相差計測方法は、２つの被計測信号をそれぞれＡ／Ｄ変換するデジタル変換ステップと、デジタル変換後の２つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記２つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記２つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする。
【００１５】
本発明によれば、ゼロクロスのような特定の条件を抽出する処理が存在せず、アナログ信号をデジタル信号に変換する変換手段（Ａ／Ｄ変換）における誤差は、全体に正負が統計的に等しい状態で拡散され、誤差が大幅に低減されることとなる。
【００１６】
本発明の位相差計測システムは、２つの被計測信号をそれぞれＡ／Ｄ変換するＡ／Ｄ変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の２つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出し、前記２つの被計測信号それぞれの振幅値から前記面積を正規化し、正規化された前記面積から前記２つの被計測信号の位相差を算出することを特徴とする。
【００１７】
本発明によれば、リサージュ図を描画可能な１周期の被計測信号波形があれば、コンピュータ上で動作する前記計測プログラムによって前記２つの被計測信号の位相差が計算できることから、短時間で位相差計測を行なうことができる。
【００１８】
本発明の周波数計測方法は、被計測信号をフィルタリングし、このフィルタリング前後での信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後の２つの信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記２つのデジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする。
【００１９】
本発明によれば、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、前記周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を正確に求めることができる。
【００２０】
本発明の周波数計測方法は、被計測信号をＡ／Ｄ変換するデジタル変換ステップと、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記デジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする。
【００２１】
本発明の周波数計測システムは、被計測信号をＡ／Ｄ変換するＡ／Ｄ変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出し、前記被計測信号の振幅値から前記面積を正規化し、正規化された前記面積から前記デジタル信号の位相差を算出し、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定することを特徴とする。
【００２２】
本発明によれば、ゼロクロスのような特定の条件を抽出する処理が存在せず、アナログ信号をデジタル信号に変換する変換手段（Ａ／Ｄ変換）における誤差は、全体に正負が統計的に等しい状態で拡散され、誤差が大幅に低減されることとなる。
【００２３】
本発明は、前記フィルタリングに、IIRくし型フィルタを用いることを特徴とする。
【００２４】
本発明によれば、IIRくし型フィルタ（フィードバック型くし型フィルタとも呼ばれる）を用いることで、設定されたノッチ周波数近傍で入出力間の位相を急峻に変化させるので、周波数の最小分解能が向上する。
【００２５】
本発明は、前記フィルタリングに、線形位相FIRフィルタを用いることを特徴とする。
【００２６】
本発明によれば、線形位相FIRフィルタを用いることで、直線性が良くなり、計測可能な周波数帯域が広帯域となる。
【発明の効果】
【００２７】
本発明の位相差計測方法と位相差計測システムによれば、リサージュ図の閉曲線の囲む面積を算出することで、リサージュ図を構成する全データが計算に使用されることとなり、精密な位相差計測ができる。そして、ゼロクロスのような特定の条件を抽出する処理が存在せず、アナログ信号をデジタル信号に変換する変換手段（Ａ／Ｄ変換）における誤差は、全体に正負が統計的に等しい状態で拡散され、誤差が大幅に低減されることとなる。また、本発明によれば、閉曲線積分の性質から原点は軌跡の外であってもかまわない。信号にオフセットが加わっても面積の計算には反映されない。ノイズが加わっても面積は平均化され影響は緩和される。その結果、しきい値周辺のデータに依存するゼロクロス法の欠点が解消される。そして、本発明によれば、リサージュ図を描画可能な１周期の被計測信号波形があれば、コンピュータ上で動作する前記計測プログラムによって前記２つの被計測信号の位相差が計算できることから、短時間で位相差計測を行なうことができる。
【００２８】
本発明の周波数計測方法と周波数計測システムによれば、前記位相差と前記フィルタリングにおける移相特性から周波数偏差を算出し、前記周波数偏差と前記フィルタリングにおけるノッチ周波数とから前記被計測信号の周波数を正確に求めることができる。そして、前記フィルタリングに、IIRくし型フィルタを用いることで、設定されたノッチ周波数近傍で入出力間の位相を急峻に変化させるので、周波数の最小分解能が向上する。または、線形位相FIRフィルタを用いることで、直線性が良くなり、計測可能な周波数帯域が広帯域となる。
【図面の簡単な説明】
【００２９】
【図１】本発明の実施の形態の位相差計測システムを説明するブロック図である。
【図２】上記位相差計測システムに係る計測プログラムの構成を例示するブロック図である。
【図３】リサージュ図から２つの信号の位相差を求める本発明の方法を示す図である。
【図４】本発明の実施の形態の位相差計測手順を示すフローチャート図である。
【図５】本発明の実施の形態の周波数計測システムを説明するブロック図である。
【図６】上記周波数計測システムに係る計測プログラム構成を例示するブロック図である。
【図７】上記周波数計測システムに係る計測プログラム構成を例示するブロック図である。
【図８】本発明の実施の形態の周波数計測手順を示すフローチャート図である。
【図９】ＩＩＲくし型フィルタの構造図である。
【図１０】ＩＩＲくし型フィルタの特性図であり、（ａ）は移相特性図であり、（ｂ）は振幅特性図である。
【図１１】上記周波数計測システムによる位相差分解能のばらつきを示す時系列特性図である。
【図１２】上記周波数計測システムによる位相差分解能のばらつきを示す度数分布図である。
【図１３】上記周波数計測システムにおける測定毎の計測値を平均化した際の周波数偏差散布図である。
【図１４】微細周波数差によって発生する位相変化を示す図である。
【発明を実施するための形態】
【００３０】
本発明に係る位相差計測原理について以下に説明する。
【００３１】
（位相差計測原理）
信号ｘと信号ｙの周波数比が１：１，位相差φである場合のリサージュ図は、位相差及び振幅に応じて楕円／円／直線となる。これを閉曲線と捉えて、その軌道の描く面積から位相差sin(φ)を求める方法を以下に説明する。
【００３２】
まず、閉曲線の囲む面積の一般式を数１に示す。
【００３３】
【数１】
0005674042.tif 000002
【００３４】
図３に示すリサージュ図となる閉曲線の囲む面積を符号Ｓとすると、軌跡中、隣接する２点と原点がつくる微小な三角形の部分面積ｄＳは、軌跡と軌跡の移動ベクトルとの外積を用いて、次の数式２または数式３で表される。この部分面積ｄＳは、軌跡のベクトル（＝Ｐ）と移動ベクトル（＝ｄＰ）の外積をとり１/２にすることで得られる。
【００３５】
【数２】
0005674042.tif 000003
【００３６】
【数３】
0005674042.tif 000004
【００３７】
部分面積ｄＳを外積の計算式にしたがって計算すると、軌跡が反時計回りの時はｄＳ＞０、時計回りの時はｄＳ＜０となる。また, 全体の面積ＳはｄＳを総計し次の数式４で表される。周期Ｔは被計測信号の周波数とサンプリング条件によって変化する。
【００３８】
【数４】
0005674042.tif 000005
【００３９】
ここでは「位相差φの信号が描くリサージュ面積」を導く。位相差φのある信号ｑ１，ｑ２を考える。符号ａ１，ａ２は振幅、符号ω０は角速度であるとすると、次の数式５で表される。
【００４０】
【数５】
0005674042.tif 000006
【００４１】
三角関数の加法定理を用いて、数式５およびその微分形を数式１に代入・変形すると面積Ｓは次の数式６で表される。ここで、積分区間は１周期とする。描画される面積は位相差と振幅で決まることが示される。
【００４２】
【数６】
0005674042.tif 000007
【００４３】
前述の振幅（ａ１，ａ２）が別の方法で計測されたと仮定し、得られた閉曲線面積Ｓを、次の数式７で表すように正規化し、正規化リサージュ面積（ＬｓＲｏｔ）を得る。
【００４４】
【数７】
0005674042.tif 000008
【００４５】
上記数式７で得られた面積は、位相差＝０度の時に面積＝０、位相差＝±９０度の時に面積＝πとなるよう正規化された面積である。−π≦ＬｓＲｏｔ≦πとなり無次元量となる。上記数式６と数式７より、次の数式８が得られる。
【００４６】
【数８】
0005674042.tif 000009
【００４７】
上記数式８から、位相差ｓｉｎ（φ）は正規化リサージュ面積（ＬｓＲｏｔ）から容易に得られることがわかる。
【００４８】
ここで、上述した、図３に示すリサージュ図となる閉曲線の囲む部分面積ｄＳは、軌跡と軌跡の移動ベクトルとの外積を用いて、前述の数式２または数式３で表される。しかし、外積の性質から、Ｐ×Ｐ≡０（恒等的にゼロ）であり、前記差分をとる計算は不要であることから、数式９を導くことができる。
【００４９】
【数９】
0005674042.tif 000010
【００５０】
本発明に係る位相差計測手順、位相差計測システム、及び位相差計測プログラムについて以下に説明する。
【００５１】
（位相差計測手順）
図４は、本発明の実施形態の位相差計測手順を示すフローチャート図である。本実施形態の位相差計測方法は、２つの被計測信号ｑ１，ｑ２をそれぞれＡ／Ｄ変換するデジタル変換ステップＳ１１と、デジタル変換後の２つの被計測信号ｑｄ１，ｑｄ２からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップＳ２１と、前記２つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップＳ３１と、正規化された前記面積から前記２つの被計測信号の位相差を算出する位相差算出ステップＳ４１からなる。  本実施形態によれば、リサージュ図を描画可能な１周期の被計測信号波形があれば、コンピュータ３上で動作する計測プログラム２によって前記２つの被計測信号の位相差が計算でき、短時間で位相差計測を行なうことができる。
【００５２】
（位相差計測システム）
図１は、本発明の実施の形態の位相差計測システム１を機能的に説明するブロック図である。本実施形態の位相差計測システム１は、２つの被計測信号ｑ１、ｑ２をそれぞれＡ／Ｄ変換するＡ／Ｄ変換器７と、計測プログラム２を備えたコンピュータ３から構成される（図１）。
【００５３】
図１に示す例では、コンピュータ３には、データ入力を受け付ける入力手段８１と、入力されたデータに演算処理を行うＣＰＵ４と、演算処理後のデータを蓄積保存するデータベース５とデータ処理前後のデータを画面表示するディスプレイ８が備わっており、コンピュータ３上で計測プログラム２が動作する。計測プログラム２は、デジタルフィルタプログラム２２と位相差計測プログラム２３からなる。
【００５４】
（位相差計測プログラム）
図２は、上記実施形態の位相差計測システム１に係る計測プログラム２の構成を例示するブロック図である。図２に示す例では、計測プログラム２は、デジタルフィルタプログラム２２と位相差計測プログラム２３からなる。デジタル変換ステップＳ１１にてデジタル変換された２つの被計測信号ｑｄ１，ｑｄ２がコンピュータ３に入力されると、それぞれデジタルフィルタ２２２を通して、リサージュ図となる閉曲線の囲む面積を算出する面積計算処理を行い（符号Ｓ２１）、２つの被計測信号ｑｄ１，ｑｄ２のそれぞれの振幅計算処理を行い（符号Ｓ２１１）、前記２つの被計測信号それぞれの振幅値から前記面積を正規化する面積正規化処理を行い（符号Ｓ３１）、正規化された前記面積から前記２つの被計測信号の位相差を算出する位相差算出処理（符号Ｓ４１）を行なった後、位相差信号データａｎ１２として出力する。なお、ここでは、リサージュ図を直接的に描画することはしていないが、例えば、面積計算処理Ｓ２１の前段階でリサージュ図を描画してディスプレイ表示する構成としても良い。
【００５５】
次に、本発明の周波数計測におけるフィルタについて説明する。
【００５６】
（ＩＩＲくし型フィルタ）
一般に、ＩＩＲくし型フィルタの構造は、図９のように示される。遅延段数をＫ、フィードバックゲインをａとすると、伝達関数は次の数式１０で表される。ｊは虚数単位である。
【００５７】
【数１０】
0005674042.tif 000011
【００５８】
周波数と角速度をサンプリング周波数ｆｓで正規化すると、正規化周波数Ｆ、正規化角速度Ωは次の数式１１で表される。
【００５９】
【数１１】
0005674042.tif 000012
【００６０】
前述の正規化角周波数と正規化角速度を用い、伝達関数は正規化された角速度を用いて次の数式１２で表される。以後、この正規化した表現を基準に説明する。
【００６１】
【数１２】
0005674042.tif 000013
【００６２】
ＩＩＲくし型フィルタの移相φは前述の伝達関数から、次の数式１３で表される。次に、移相特性φがΩの関数として表せるように数１３を変形し、周波数Ωによってどのくらいの位相変化φが発生するかを、フィルタ移相特性φ（Ω）で表す。フィルタ移相特性φ(Ω)は、次の数式１４で表される。
【００６３】
【数１３】
0005674042.tif 000014
【００６４】
【数１４】
0005674042.tif 000015
【００６５】
図１０は、ＩＩＲくし型フィルタの特性図である。図１０（ａ）は移相特性図であり、グラフの縦軸は入出力位相差を示しており、グラフの横軸は相対周波数を示している。図１０（ｂ）は振幅特性図であり、グラフの縦軸は振幅を示しており、グラフの横軸は相対周波数を示している。ここでは−１＜ａ＜０を使用しており、ノッチ周波数の奇数倍の周波数において極めて鋭敏な移相特性とノッチ周波数近傍のみを通過させる振幅特性となっている。
【００６６】
上記数式１４をテーラー展開して得られる３次項以上を切り捨て、１次近似する。ノッチ周波数近傍での移相φの近似式は、次の数式１５で表される。
【００６７】
【数１５】
0005674042.tif 000016
【００６８】
ここで、位相変化率を求める。ｆ０近傍における位相変化率は、φをΩで微分することで得られ次の数式１６で表される。
【００６９】
【数１６】
0005674042.tif 000017
【００７０】
次に、得られた位相変化率から周波数偏差を求める。Δｆはｆ０からの変位とする。ΔφとΔｆの関係は、次の数式１７で表される。
【００７１】
【数１７】
0005674042.tif 000018
【００７２】
上記数式１７から、計測される位相差は周波数偏差がπｇ０倍されて計測できることがわかる（感度の向上）。上記数式１７を変形すると、次の数式１８のようにも書ける。
【００７３】
【数１８】
0005674042.tif 000019
【００７４】
上記数式１８を移項し周波数偏差（Δｆ／ｆ０）は、次の数式１９で表される。
【００７５】
【数１９】
0005674042.tif 000020
【００７６】
くし形フィルタのノッチ周波数ｆ０はサンプリング周波数ｆｓと遅延段数のみで決定できる（数式２１を参照）。  前述の正規化リサージュ面積を精密に定量したとすると、数式１９と周波数偏差Δｆとｆ０の関係を合わせ、数式２０の式が得られる。数式２０，数式２１のように単純な式を用いて精密な周波数判別が可能である。
【００７７】
【数２０】
0005674042.tif 000021
【００７８】
【数２１】
0005674042.tif 000022
【００７９】
得られた数式２０，数式２１を見るとサンプリング周波数が基準の周波数となっている（国家標準を用いて校正可能である）。トレーサブルな周波数標準を計測することでサンプリング周波数のキャリブレーションが可能である。つまり、得られた数式２０，数式２１を見ると、Ｋ，ｇ０π２は純粋な定数である。  正規化リサージュ面積（ＬｓＲｏｔ）は計測値である。ｆｓはハードウェアに実装される水晶発振子の精度により実際の装置で動作している周波数の実測値が決まる。  よってｆｓを正しく校正することで（国家標準を用いて校正可能である）、実際に動作しているサンプリング周波数（ｆｓ）の確度を高めることができ、計算結果（理論値）と実際の値（計測値）を近づけることができる。
【００８０】
本発明に係る周波数計測手順、周波数計測システム、及び周波数計測プログラムについて以下に説明する。
【００８１】
（周波数計測手順）
図８は、本発明の実施形態の周波数計測手順を示すフローチャート図である。本実施形態の周波数計測方法は、被計測信号（例えばｑ１）をＡ／Ｄ変換するデジタル変換ステップＳ１１と、デジタル変換後の被計測信号（例えばｑｄ１）をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力する位相差生成ステップＳ１５と、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップＳ２１と、前記被計測信号の振幅値から前記面積を正規化する正規化ステップＳ３１と、正規化された前記面積から前記デジタル信号の位相差を算出する位相差算出ステップＳ４１と、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップＳ５１からなる。
【００８２】
（周波数計測システム）
図５は、本発明の実施形態の周波数計測システム１１を説明するブロック図である。本実施形態の周波数計測システム１１は、被計測信号ｑ１，ｑ２，・・・，ｑｋ（ｋは正の整数）をそれぞれＡ／Ｄ変換するＡ／Ｄ変換器７と、計測プログラム２を備えたコンピュータ３から構成される（図５）。ここで、同一の符号は同じ機能を表しており、その説明を適宜省略する。図５に示す例では、多チャンネルの被計測信号を入力する場合で説明しているが、一つの被計測信号（例えばｑ１）としても良いし、二つの被計測信号（例えばｑ１，ｑ２）としても良い。
【００８３】
（第１の周波数計測プログラム）
図６は、上記実施形態の周波数計測システム１１に係る計測プログラム２の構成を例示するブロック図である。図６は、一つの被計測信号ｑ１を入力する場合の例であり、計測プログラム２は、デジタルフィルタプログラム２２と周波数計測プログラム２４からなる。デジタル変換ステップＳ１１にてデジタル変換された被計測信号ｑｄ１がコンピュータ３に入力されると、デジタルフィルタ２２２を通し、このフィルタリング前後でのデジタル信号の位相差を変換出力して位相差を生成し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積計算処理を行い（符号Ｓ２１）、被計測信号の振幅計算処理を行い（符号Ｓ２１１）、前記被計測信号の振幅値から前記面積を正規化する面積正規化処理を行い（符号Ｓ３１）、正規化された前記面積から前記２つの被計測信号の位相差を算出する位相差算出処理（符号Ｓ４１）を行い、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定処理を行ない（符号Ｓ５１）、その結果を周波数信号データａｎ１として出力する。なお、ここでは、リサージュ図を直接的に描画することはしていないが、例えば、面積計算処理Ｓ２１の前段階でリサージュ図を描画してディスプレイ表示する構成としても良い。
【００８４】
本実施形態の周波数計測システム１１では、前記デジタルフィルタとしてＩＩＲくし型フィルタを用いている。これは、ＩＩＲくし型フィルタを用いることで、設定されたノッチ周波数近傍で入出力間の位相を急峻に変化させて、周波数の最小分解能を向上させるためである。
【００８５】
（第２の周波数計測プログラム）
図７は、上記実施形態の周波数計測システム１１に係る計測プログラム２の構成の他の例を示すブロック図である。図７に示す例では、計測プログラム２は、デジタルフィルタプログラム２２と周波数計測プログラム２４からなる。図７は、２チャンネルの被計測信号ｑ１，ｑ２を入力する場合の例である。デジタル変換ステップＳ１１にてデジタル変換された被計測信号ｑｄ１、ｑｄ２がそれぞれコンピュータ３に入力されると、前処理フィルタ２２１を通してデジタルフィルタ２２２（くし型フィルタ２２２）を通す。このくし型フィルタ２２２によるフィルタリング前後でのデジタル信号の位相差を変換出力して位相差を生成し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積計算処理を行い（符号Ｓ２１）、被計測信号の振幅計算処理を行い（符号Ｓ２１１）、前記被計測信号の振幅値から前記面積を正規化する面積正規化処理を行い（符号Ｓ３１）、正規化された前記面積から前記２つの被計測信号の位相差を算出する位相差算出処理（符号Ｓ４１）を行い、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定処理を行ない（符号Ｓ５１）、その結果をそれぞれ周波数信号データａｎ１，ａｎ２として出力する。そして、それぞれデジタルフィルタ２２２を通されたデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積計算処理を行い（符号Ｓ２１）、２つの被計測信号ｑｄ１，ｑｄ２のそれぞれの振幅計算処理を行い（符号Ｓ２１１）、前記２つの被計測信号それぞれの振幅値から前記面積を正規化する面積正規化処理を行い（符号Ｓ３１）、正規化された前記面積から前記２つの被計測信号の位相差を算出する位相差算出処理（符号Ｓ４１）を行なった後、位相差信号データａｎ１２として出力する。
【００８６】
（実施例）
図５と図７に示すブロック図の周波数計測システムを試作した。２チャンネルの入力をもち、それぞれの精密周波数とチャンネル間の精密位相差を計測する構成である。
【００８７】
前処理フィルタ２２１は目的外の周波数帯域をカットするために用いた。くし型フィルタ２２２は、特定のノッチ周波数が遅延段数Ｋによって設定されるＩＩＲくし型フィルタを用いた。周波数計測システム１１は信号周波数ｆ１，ｆ２を計測し、それと同時に２つの信号の位相差を得る。サンプリングは１６ｂｉｔ，４４１００Ｈｚとした。
【００８８】
次に、上記実施例の周波数計測システム１１を用いて最小位相差分解能・分解能ばらつき・精密周波数の精度について評価を行なった。  試験信号としてルビジウム（Rb）発振器（10MHz）を外部周波数基準としたファンクションジェネレータ（以下FGと略す）を使用した。信号には２チャンネル４４１００Ｈｚを用いた。Ａ／Ｄ変換器７を外部接続して使用した。
【００８９】
今回、意図的に位相差ゼロの状態の信号を作成し本ソフトウェアによる計測を行った。測定は１０分間行った。数値の平均と標準偏差から分解能を推定した。短時間平均（約０．６秒）で１．１μラジアン（±３σ）、長時間平均（約６秒）で０．３３μラジアン（±３σ）の分解能が確認できた。
【００９０】
図１１は、上記周波数計測システム１１による位相差分解能のばらつきを示す時系列特性図である。これは信号発生器で生成した２つの信号（ｑ１、ｑ２）の位相差を実際に観測したものである。短時間平均（約０．６秒）と長時間平均（約６秒）で挙動をみたところ、おおむね±１μラジアンの範囲に収まっていることが確認された。  図１２は、上記実施例の周波数計測システム１１による位相差分解能のばらつきを示す度数分布図である。図１１の値から長時間平均データをヒストグラムにプロットしたものである。標準偏差σは１．０８５Ｅ−０７であった。
【００９１】
次に、周波数精度を検証した。図１３は、上記周波数計測システム１１における測定毎の計測値を平均化した際の周波数偏差散布図である。２チャンネルで周波数測定を行い、周波数のばらつきを分析したところ、１．５μＨｚ（±３σ）であった。小数点以下の周波数は、μＨｚの桁まで一致することがわかった。周波数の平均値として４４０．９９８，９５８，１７Ｈｚ（σ=０．５２μＨｚ）が得られた。
設定した周波数４４１．００Ｈｚとの誤差はＡ／Ｄ変換器に搭載されるサンプリングクロック用水晶発振子の周波数偏差やドリフトが原因である。Ａ／Ｄ変換器の安定化が十分でない場合、基準となるサンプリング周波数がまずドリフトする。そのため、本実装例における両チャンネルの瞬時周波数も同一方向にドリフトし、相関係数は１に近い値（強い相関）となる。さまざまな安定化の注意を払い、最終的に安定化後の相関係数は０．４６であった。このことは、両チャンネルの誤差の一部が独立したプロセスによって発生していることを示唆している。
【００９２】
図１４は、微細周波数変動による位相変化を示す図である。  ある線路長おける伝送遅延による位相差の検出を試みる。光速は３．０ｅ＋８ｍ／ｓｅｃ．と有限であり, わずかに伝送遅延を発生させる。
試験に用いた周波数４４１．００Ｈｚの波長λは６８０２７２ｍである。もし1μラジアンの位相差が検出可能ならば、線路長の約１０ｃｍの差を検出できる。実験では長さが５７ｃｍ異なる２本のケーブルを用いてファンクションジェネレータとＡ／Ｄ変換器を接続し、位相差φ１を計測した。次にケーブルを入れ替えて接続し位相差φ２を計測した。φ１＝５２．５μラジアン、φ２＝８１０μラジアンの結果が得られた。位相差Δφ=７５７μラジアンとなり、相対的線路長差（＝１．１４ｍ）における理論値である１０．５μラジアンをはるかに上回る。この結果は光速の限界による伝送遅延に由来する位相差ではなく、おそらく線路の持つインピーダンスによる伝送遅延の位相差を捉えている。実験ではケーブル長の差や遅延量も十分に考慮しなければならないことを改めて示している。
【００９３】
以上、本発明が適用できる計測システムとしては、周波数カウンタに代わる周波数計測器、信号発生器に内蔵される位相差計測処理部、位相差計測を利用するインピーダンスメータ、音響機器の周波数特性のうち特に位相遅延特性の計測、アレーマイクロホンを使用する信号処理における位相差計測、レーザー光線を使用した距離計測器における位相差計測を用いる時間遅延測定、音響周波数帯における微小なドップラー効果を利用する計測器、歌声や楽器の精密音程判別、音声処理における基本周波数推定、ピアノの調律で使用される精密チューナー、ピアノ調律におけるインハーモニシティの精密計測など、多彩な用途がある。
特に、低周波では１波長〜数波長のデータでも小数点以下の周波数が短時間で得られることから、電力測定器における周波数測定・力率測定、電力系統周波数（５０Ｈｚ、６０Ｈｚ）の短時間精密計測、電力系統における特定地点の位相差を精密に測定することで得られる電力潮流推定などに適用可能である。
【符号の説明】
【００９４】
１          位相差計測システム、
２          計測プログラム、
２２        デジタルフィルタプログラム、
２３        位相差計測プログラム、
２４        周波数計測プログラム、
３          コンピュータ、
４          ＣＰＵ、
７          Ａ／Ｄ変換器、
１１        周波数計測システム、
ｑ１，ｑ２，・・・，ｑｋ      被計測信号（ｋ＝正の整数）
【図１】
0005674042.tif 000023
【図２】
0005674042.tif 000024
【図３】
0005674042.tif 000025
【図４】
0005674042.tif 000026
【図５】
0005674042.tif 000027
【図６】
0005674042.tif 000028
【図７】
0005674042.tif 000029
【図８】
0005674042.tif 000030
【図９】
0005674042.tif 000031
【図１０】
0005674042.tif 000032
【図１１】
0005674042.tif 000033
【図１２】
0005674042.tif 000034
【図１３】
0005674042.tif 000035
【図１４】
0005674042.tif 000036