位相差計測方法及び位相差計測システム、並びに周波数計測方法及び周波数計測システム
【課題】 リサージュ法による処理の仕方を改良することで、精度良く位相差計測を行なうことが可能な位相差計測システムを提供する。
【解決手段】 位相差計測システム1は、2つの被計測信号q1、q2をそれぞれA/D変換するA/D変換器7と、計測プログラム2を備えたコンピュータ3から構成され、前記計測プログラム2が、デジタル変換後の2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出し、前記2つの被計測信号それぞれの振幅値から前記面積を正規化し、正規化された前記面積から前記2つの被計測信号の位相差を算出する。
【解決手段】 位相差計測システム1は、2つの被計測信号q1、q2をそれぞれA/D変換するA/D変換器7と、計測プログラム2を備えたコンピュータ3から構成され、前記計測プログラム2が、デジタル変換後の2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出し、前記2つの被計測信号それぞれの振幅値から前記面積を正規化し、正規化された前記面積から前記2つの被計測信号の位相差を算出する。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、リサージュ図となる閉曲線から被計測信号の位相差を求める位相差計測方法と位相差計測システム、並びに、リサージュ図となる閉曲線から被計測信号の周波数を求める周波数計測方法及び周波数計測システムに関する。
【背景技術】
【0002】
信号処理の分野において基準周波数との位相差(位相角)を計測することは重要なコア技術である。例えば、電力・電力量の位相差(位相角)測定、設備や機器の振動分析による異常診断や、静音性の評価や騒音原因の対策方法を検討するための騒音分析など様々な分野、更には音響処理における音源推定や楽音処理における倍音の分析において位相差の計測は重要である。
【0003】
従来、位相差計測方法としては、リサージュ法(Lissajous method)によりリサージュ図を描画する方法や、FFT(Fast Fourier Transform)アナライザを用いた方法や、波形から周期を測定する方法(電力・電力量の測定において使用されている)や、波形の部分間の相関を利用する方法(音声処理系で使用されている)などが知られている。
【0004】
既知の位相差計測方法のうち、周波数校正においては、微小な周波数偏差を検出するための簡便な方法としてリサージュ法が用いられることがある。
リサージュ図は、計測表示装置をX−Y入力モードに設定し、二つの単振動波形の順序対として得られる点の軌跡が描く平面図形である。周波数比が1の場合、位相差に応じて、形状は楕円・円・直線となり、周波数比が整数の場合は複数の山が描かれる。
【0005】
リサージュ図から2つの被計測信号の位相差を求める従来の方法としては、例えば、計測表示装置のX軸とY軸に、それぞれ周波数の等しい正弦波を入力すると、楕円形のリサージュ図が現れる。このときの位相差を正確に計算するためには図形がX軸(又はY軸)を横切る点をそれぞれ正確に求めることとなる。そして、リサージュ図を描く閉曲線のY軸方向の最大振幅と、前記閉曲線がY軸を横切る際のY軸方向の振幅とから位相差を計算する。
【0006】
特許文献1には、多軸加振試験機の位相差制御方法に関し、複数の系の信号を用いて表示装置上に正規化したリサージュ図を重ねて描かせることで入力信号の位相差と出力信号の位相差が揃っているか否かを検出することが記載されている。
【0007】
特許文献2には、光干渉計の周期誤差低減方法に関し、極座標変換した直交2相信号の位相角と半径の変動値を、誤差を求める各式に代入して前記直交2相信号の振幅、零点及び位相差の補正値を求め、極座標の半径が一定となるように補正することが記載されている。そして、光干渉計の動作原理を説明するために、リサージュ図が説明されている。
【先行技術文献】
【特許文献】
【0008】
【特許文献1】特開2000−298149号公報
【特許文献2】特許第4465451号公報
【発明の概要】
【発明が解決しようとする課題】
【0009】
しかしながら、上記従来のリサージュ法では、位相差を正確に計算するためには図形がX軸(又はY軸)を横切る点の座標をそれぞれ正確に求めなければならない。デジタル信号処理によって補間処理を行うことでゼロクロス点を求めることができるが、軌跡が原点に近づくにつれ原点近傍を通過する相対的な速度が大きくなり補間誤差が低下してしまう。
そして、従来、周波数カウンタでは「レシプロカル方式」を採用しており、小数点以下の周波数の計測が可能となっているが、周期を精密に計測するためには非常に高いサンプリング周波数を使用し周期を計測する必要がある。また周期の基準となるための「ゼロクロス点」の検出精度はDCオフセットの影響を受けてしまい検出精度が低下する懸念がある。
また、上記従来のリサージュ法では、微小な周波数偏差、例えば0.1/1000Hzを認識するためには360/10000(度/秒)で変化する位相差を短時間に認識し周波数の微調整を繰り返す必要がある。
【0010】
そこで、本願発明の目的は、従来のリサージュ法による処理の仕方を改良することで、精度良く位相差計測を行なうことが可能な位相差計測方法及び位相差計測システムと、精度良く周波数計測を行なうことが可能な周波数計測方法及び周波数計測システムを提供することにある。
【課題を解決するための手段】
【0011】
本願発明者は、リサージュ図が描く閉曲線の囲む面積に着目し、短時間で精密な位相差を検出できる方法を見出し、本願発明を完成するに至った。
【0012】
本発明の位相差計測方法は、2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする。
【0013】
従来のリサージュ法では、リサージュ図のうちのゼロクロスする近傍のデータしか計算に使用されていないが、本発明によれば、リサージュ図の閉曲線の囲む面積を算出することで、リサージュ図を構成する全データが計算に使用されることとなり、精密な位相差計測ができる。
【0014】
本発明の位相差計測方法は、2つの被計測信号をそれぞれA/D変換するデジタル変換ステップと、デジタル変換後の2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする。
【0015】
本発明によれば、ゼロクロスのような特定の条件を抽出する処理が存在せず、アナログ信号をデジタル信号に変換する変換手段(A/D変換)における誤差は、全体に正負が統計的に等しい状態で拡散され、誤差が大幅に低減されることとなる。
【0016】
本発明の位相差計測システムは、2つの被計測信号をそれぞれA/D変換するA/D変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出し、前記2つの被計測信号それぞれの振幅値から前記面積を正規化し、正規化された前記面積から前記2つの被計測信号の位相差を算出することを特徴とする。
【0017】
本発明によれば、リサージュ図を描画可能な1周期の被計測信号波形があれば、コンピュータ上で動作する前記計測プログラムによって前記2つの被計測信号の位相差が計算できることから、短時間で位相差計測を行なうことができる。
【0018】
本発明の周波数計測方法は、被計測信号をフィルタリングし、このフィルタリング前後での信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後の2つの信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つのデジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする。
【0019】
本発明によれば、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、前記周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を正確に求めることができる。
【0020】
本発明の周波数計測方法は、被計測信号をA/D変換するデジタル変換ステップと、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記デジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする。
【0021】
本発明の周波数計測システムは、被計測信号をA/D変換するA/D変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出し、前記被計測信号の振幅値から前記面積を正規化し、正規化された前記面積から前記デジタル信号の位相差を算出し、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定することを特徴とする。
【0022】
本発明によれば、ゼロクロスのような特定の条件を抽出する処理が存在せず、アナログ信号をデジタル信号に変換する変換手段(A/D変換)における誤差は、全体に正負が統計的に等しい状態で拡散され、誤差が大幅に低減されることとなる。
【0023】
本発明は、前記フィルタリングに、IIRくし型フィルタを用いることを特徴とする。
【0024】
本発明によれば、IIRくし型フィルタ(フィードバック型くし型フィルタとも呼ばれる)を用いることで、設定されたノッチ周波数近傍で入出力間の位相を急峻に変化させるので、周波数の最小分解能が向上する。
【0025】
本発明は、前記フィルタリングに、線形位相FIRフィルタを用いることを特徴とする。
【0026】
本発明によれば、線形位相FIRフィルタを用いることで、直線性が良くなり、計測可能な周波数帯域が広帯域となる。
【発明の効果】
【0027】
本発明の位相差計測方法と位相差計測システムによれば、リサージュ図の閉曲線の囲む面積を算出することで、リサージュ図を構成する全データが計算に使用されることとなり、精密な位相差計測ができる。そして、ゼロクロスのような特定の条件を抽出する処理が存在せず、アナログ信号をデジタル信号に変換する変換手段(A/D変換)における誤差は、全体に正負が統計的に等しい状態で拡散され、誤差が大幅に低減されることとなる。また、本発明によれば、閉曲線積分の性質から原点は軌跡の外であってもかまわない。信号にオフセットが加わっても面積の計算には反映されない。ノイズが加わっても面積は平均化され影響は緩和される。その結果、しきい値周辺のデータに依存するゼロクロス法の欠点が解消される。そして、本発明によれば、リサージュ図を描画可能な1周期の被計測信号波形があれば、コンピュータ上で動作する前記計測プログラムによって前記2つの被計測信号の位相差が計算できることから、短時間で位相差計測を行なうことができる。
【0028】
本発明の周波数計測方法と周波数計測システムによれば、前記位相差と前記フィルタリングにおける移相特性から周波数偏差を算出し、前記周波数偏差と前記フィルタリングにおけるノッチ周波数とから前記被計測信号の周波数を正確に求めることができる。そして、前記フィルタリングに、IIRくし型フィルタを用いることで、設定されたノッチ周波数近傍で入出力間の位相を急峻に変化させるので、周波数の最小分解能が向上する。または、線形位相FIRフィルタを用いることで、直線性が良くなり、計測可能な周波数帯域が広帯域となる。
【図面の簡単な説明】
【0029】
【図1】本発明の実施の形態の位相差計測システムを説明するブロック図である。
【図2】上記位相差計測システムに係る計測プログラムの構成を例示するブロック図である。
【図3】リサージュ図から2つの信号の位相差を求める本発明の方法を示す図である。
【図4】本発明の実施の形態の位相差計測手順を示すフローチャート図である。
【図5】本発明の実施の形態の周波数計測システムを説明するブロック図である。
【図6】上記周波数計測システムに係る計測プログラム構成を例示するブロック図である。
【図7】上記周波数計測システムに係る計測プログラム構成を例示するブロック図である。
【図8】本発明の実施の形態の周波数計測手順を示すフローチャート図である。
【図9】IIRくし型フィルタの構造図である。
【図10】IIRくし型フィルタの特性図であり、(a)は移相特性図であり、(b)は振幅特性図である。
【図11】上記周波数計測システムによる位相差分解能のばらつきを示す時系列特性図である。
【図12】上記周波数計測システムによる位相差分解能のばらつきを示す度数分布図である。
【図13】上記周波数計測システムにおける測定毎の計測値を平均化した際の周波数偏差散布図である。
【図14】微細周波数差によって発生する位相変化を示す図である。
【発明を実施するための形態】
【0030】
本発明に係る位相差計測原理について以下に説明する。
【0031】
(位相差計測原理)
信号xと信号yの周波数比が1:1,位相差φである場合のリサージュ図は、位相差及び振幅に応じて楕円/円/直線となる。これを閉曲線と捉えて、その軌道の描く面積から位相差sin(φ)を求める方法を以下に説明する。
【0032】
まず、閉曲線の囲む面積の一般式を数1に示す。
【0033】
【数1】
【0034】
図3に示すリサージュ図となる閉曲線の囲む面積を符号Sとすると、軌跡中、隣接する2点と原点がつくる微小な三角形の部分面積dSは、軌跡と軌跡の移動ベクトルとの外積を用いて、次の数式2または数式3で表される。この部分面積dSは、軌跡のベクトル(=P)と移動ベクトル(=dP)の外積をとり1/2にすることで得られる。
【0035】
【数2】
【0036】
【数3】
【0037】
部分面積dSを外積の計算式にしたがって計算すると、軌跡が反時計回りの時はdS>0、時計回りの時はdS<0となる。また, 全体の面積SはdSを総計し次の数式4で表される。周期Tは被計測信号の周波数とサンプリング条件によって変化する。
【0038】
【数4】
【0039】
ここでは「位相差φの信号が描くリサージュ面積」を導く。位相差φのある信号q1,q2を考える。符号a1,a2は振幅、符号ω0は角速度であるとすると、次の数式5で表される。
【0040】
【数5】
【0041】
三角関数の加法定理を用いて、数式5およびその微分形を数式1に代入・変形すると面積Sは次の数式6で表される。ここで、積分区間は1周期とする。描画される面積は位相差と振幅で決まることが示される。
【0042】
【数6】
【0043】
前述の振幅(a1,a2)が別の方法で計測されたと仮定し、得られた閉曲線面積Sを、次の数式7で表すように正規化し、正規化リサージュ面積(LsRot)を得る。
【0044】
【数7】
【0045】
上記数式7で得られた面積は、位相差=0度の時に面積=0、位相差=±90度の時に面積=πとなるよう正規化された面積である。−π≦LsRot≦πとなり無次元量となる。上記数式6と数式7より、次の数式8が得られる。
【0046】
【数8】
【0047】
上記数式8から、位相差sin(φ)は正規化リサージュ面積(LsRot)から容易に得られることがわかる。
【0048】
ここで、上述した、図3に示すリサージュ図となる閉曲線の囲む部分面積dSは、軌跡と軌跡の移動ベクトルとの外積を用いて、前述の数式2または数式3で表される。しかし、外積の性質から、P×P≡0(恒等的にゼロ)であり、前記差分をとる計算は不要であることから、数式9を導くことができる。
【0049】
【数9】
【0050】
本発明に係る位相差計測手順、位相差計測システム、及び位相差計測プログラムについて以下に説明する。
【0051】
(位相差計測手順)
図4は、本発明の実施形態の位相差計測手順を示すフローチャート図である。本実施形態の位相差計測方法は、2つの被計測信号q1,q2をそれぞれA/D変換するデジタル変換ステップS11と、デジタル変換後の2つの被計測信号qd1,qd2からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップS21と、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップS31と、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出ステップS41からなる。 本実施形態によれば、リサージュ図を描画可能な1周期の被計測信号波形があれば、コンピュータ3上で動作する計測プログラム2によって前記2つの被計測信号の位相差が計算でき、短時間で位相差計測を行なうことができる。
【0052】
(位相差計測システム)
図1は、本発明の実施の形態の位相差計測システム1を機能的に説明するブロック図である。本実施形態の位相差計測システム1は、2つの被計測信号q1、q2をそれぞれA/D変換するA/D変換器7と、計測プログラム2を備えたコンピュータ3から構成される(図1)。
【0053】
図1に示す例では、コンピュータ3には、データ入力を受け付ける入力手段81と、入力されたデータに演算処理を行うCPU4と、演算処理後のデータを蓄積保存するデータベース5とデータ処理前後のデータを画面表示するディスプレイ8が備わっており、コンピュータ3上で計測プログラム2が動作する。計測プログラム2は、デジタルフィルタプログラム22と位相差計測プログラム23からなる。
【0054】
(位相差計測プログラム)
図2は、上記実施形態の位相差計測システム1に係る計測プログラム2の構成を例示するブロック図である。図2に示す例では、計測プログラム2は、デジタルフィルタプログラム22と位相差計測プログラム23からなる。デジタル変換ステップS11にてデジタル変換された2つの被計測信号qd1,qd2がコンピュータ3に入力されると、それぞれデジタルフィルタ222を通して、リサージュ図となる閉曲線の囲む面積を算出する面積計算処理を行い(符号S21)、2つの被計測信号qd1,qd2のそれぞれの振幅計算処理を行い(符号S211)、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する面積正規化処理を行い(符号S31)、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出処理(符号S41)を行なった後、位相差信号データan12として出力する。なお、ここでは、リサージュ図を直接的に描画することはしていないが、例えば、面積計算処理S21の前段階でリサージュ図を描画してディスプレイ表示する構成としても良い。
【0055】
次に、本発明の周波数計測におけるフィルタについて説明する。
【0056】
(IIRくし型フィルタ)
一般に、IIRくし型フィルタの構造は、図9のように示される。遅延段数をK、フィードバックゲインをaとすると、伝達関数は次の数式10で表される。jは虚数単位である。
【0057】
【数10】
【0058】
周波数と角速度をサンプリング周波数fsで正規化すると、正規化周波数F、正規化角速度Ωは次の数式11で表される。
【0059】
【数11】
【0060】
前述の正規化角周波数と正規化角速度を用い、伝達関数は正規化された角速度を用いて次の数式12で表される。以後、この正規化した表現を基準に説明する。
【0061】
【数12】
【0062】
IIRくし型フィルタの移相φは前述の伝達関数から、次の数式13で表される。次に、移相特性φがΩの関数として表せるように数13を変形し、周波数Ωによってどのくらいの位相変化φが発生するかを、フィルタ移相特性φ(Ω)で表す。フィルタ移相特性φ(Ω)は、次の数式14で表される。
【0063】
【数13】
【0064】
【数14】
【0065】
図10は、IIRくし型フィルタの特性図である。図10(a)は移相特性図であり、グラフの縦軸は入出力位相差を示しており、グラフの横軸は相対周波数を示している。図10(b)は振幅特性図であり、グラフの縦軸は振幅を示しており、グラフの横軸は相対周波数を示している。ここでは−1<a<0を使用しており、ノッチ周波数の奇数倍の周波数において極めて鋭敏な移相特性とノッチ周波数近傍のみを通過させる振幅特性となっている。
【0066】
上記数式14をテーラー展開して得られる3次項以上を切り捨て、1次近似する。ノッチ周波数近傍での移相φの近似式は、次の数式15で表される。
【0067】
【数15】
【0068】
ここで、位相変化率を求める。f0近傍における位相変化率は、φをΩで微分することで得られ次の数式16で表される。
【0069】
【数16】
【0070】
次に、得られた位相変化率から周波数偏差を求める。Δfはf0からの変位とする。ΔφとΔfの関係は、次の数式17で表される。
【0071】
【数17】
【0072】
上記数式17から、計測される位相差は周波数偏差がπg0倍されて計測できることがわかる(感度の向上)。上記数式17を変形すると、次の数式18のようにも書ける。
【0073】
【数18】
【0074】
上記数式18を移項し周波数偏差(Δf/f0)は、次の数式19で表される。
【0075】
【数19】
【0076】
くし形フィルタのノッチ周波数f0はサンプリング周波数fsと遅延段数のみで決定できる(数式21を参照)。 前述の正規化リサージュ面積を精密に定量したとすると、数式19と周波数偏差Δfとf0の関係を合わせ、数式20の式が得られる。数式20,数式21のように単純な式を用いて精密な周波数判別が可能である。
【0077】
【数20】
【0078】
【数21】
【0079】
得られた数式20,数式21を見るとサンプリング周波数が基準の周波数となっている(国家標準を用いて校正可能である)。トレーサブルな周波数標準を計測することでサンプリング周波数のキャリブレーションが可能である。つまり、得られた数式20,数式21を見ると、K,g0π2は純粋な定数である。 正規化リサージュ面積(LsRot)は計測値である。fsはハードウェアに実装される水晶発振子の精度により実際の装置で動作している周波数の実測値が決まる。 よってfsを正しく校正することで(国家標準を用いて校正可能である)、実際に動作しているサンプリング周波数(fs)の確度を高めることができ、計算結果(理論値)と実際の値(計測値)を近づけることができる。
【0080】
本発明に係る周波数計測手順、周波数計測システム、及び周波数計測プログラムについて以下に説明する。
【0081】
(周波数計測手順)
図8は、本発明の実施形態の周波数計測手順を示すフローチャート図である。本実施形態の周波数計測方法は、被計測信号(例えばq1)をA/D変換するデジタル変換ステップS11と、デジタル変換後の被計測信号(例えばqd1)をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力する位相差生成ステップS15と、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップS21と、前記被計測信号の振幅値から前記面積を正規化する正規化ステップS31と、正規化された前記面積から前記デジタル信号の位相差を算出する位相差算出ステップS41と、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップS51からなる。
【0082】
(周波数計測システム)
図5は、本発明の実施形態の周波数計測システム11を説明するブロック図である。本実施形態の周波数計測システム11は、被計測信号q1,q2,・・・,qk(kは正の整数)をそれぞれA/D変換するA/D変換器7と、計測プログラム2を備えたコンピュータ3から構成される(図5)。ここで、同一の符号は同じ機能を表しており、その説明を適宜省略する。図5に示す例では、多チャンネルの被計測信号を入力する場合で説明しているが、一つの被計測信号(例えばq1)としても良いし、二つの被計測信号(例えばq1,q2)としても良い。
【0083】
(第1の周波数計測プログラム)
図6は、上記実施形態の周波数計測システム11に係る計測プログラム2の構成を例示するブロック図である。図6は、一つの被計測信号q1を入力する場合の例であり、計測プログラム2は、デジタルフィルタプログラム22と周波数計測プログラム24からなる。デジタル変換ステップS11にてデジタル変換された被計測信号qd1がコンピュータ3に入力されると、デジタルフィルタ222を通し、このフィルタリング前後でのデジタル信号の位相差を変換出力して位相差を生成し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積計算処理を行い(符号S21)、被計測信号の振幅計算処理を行い(符号S211)、前記被計測信号の振幅値から前記面積を正規化する面積正規化処理を行い(符号S31)、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出処理(符号S41)を行い、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定処理を行ない(符号S51)、その結果を周波数信号データan1として出力する。なお、ここでは、リサージュ図を直接的に描画することはしていないが、例えば、面積計算処理S21の前段階でリサージュ図を描画してディスプレイ表示する構成としても良い。
【0084】
本実施形態の周波数計測システム11では、前記デジタルフィルタとしてIIRくし型フィルタを用いている。これは、IIRくし型フィルタを用いることで、設定されたノッチ周波数近傍で入出力間の位相を急峻に変化させて、周波数の最小分解能を向上させるためである。
【0085】
(第2の周波数計測プログラム)
図7は、上記実施形態の周波数計測システム11に係る計測プログラム2の構成の他の例を示すブロック図である。図7に示す例では、計測プログラム2は、デジタルフィルタプログラム22と周波数計測プログラム24からなる。図7は、2チャンネルの被計測信号q1,q2を入力する場合の例である。デジタル変換ステップS11にてデジタル変換された被計測信号qd1、qd2がそれぞれコンピュータ3に入力されると、前処理フィルタ221を通してデジタルフィルタ222(くし型フィルタ222)を通す。このくし型フィルタ222によるフィルタリング前後でのデジタル信号の位相差を変換出力して位相差を生成し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積計算処理を行い(符号S21)、被計測信号の振幅計算処理を行い(符号S211)、前記被計測信号の振幅値から前記面積を正規化する面積正規化処理を行い(符号S31)、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出処理(符号S41)を行い、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定処理を行ない(符号S51)、その結果をそれぞれ周波数信号データan1,an2として出力する。そして、それぞれデジタルフィルタ222を通されたデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積計算処理を行い(符号S21)、2つの被計測信号qd1,qd2のそれぞれの振幅計算処理を行い(符号S211)、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する面積正規化処理を行い(符号S31)、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出処理(符号S41)を行なった後、位相差信号データan12として出力する。
【0086】
(実施例)
図5と図7に示すブロック図の周波数計測システムを試作した。2チャンネルの入力をもち、それぞれの精密周波数とチャンネル間の精密位相差を計測する構成である。
【0087】
前処理フィルタ221は目的外の周波数帯域をカットするために用いた。くし型フィルタ222は、特定のノッチ周波数が遅延段数Kによって設定されるIIRくし型フィルタを用いた。周波数計測システム11は信号周波数f1,f2を計測し、それと同時に2つの信号の位相差を得る。サンプリングは16bit,44100Hzとした。
【0088】
次に、上記実施例の周波数計測システム11を用いて最小位相差分解能・分解能ばらつき・精密周波数の精度について評価を行なった。 試験信号としてルビジウム(Rb)発振器(10MHz)を外部周波数基準としたファンクションジェネレータ(以下FGと略す)を使用した。信号には2チャンネル44100Hzを用いた。A/D変換器7を外部接続して使用した。
【0089】
今回、意図的に位相差ゼロの状態の信号を作成し本ソフトウェアによる計測を行った。測定は10分間行った。数値の平均と標準偏差から分解能を推定した。短時間平均(約0.6秒)で1.1μラジアン(±3σ)、長時間平均(約6秒)で0.33μラジアン(±3σ)の分解能が確認できた。
【0090】
図11は、上記周波数計測システム11による位相差分解能のばらつきを示す時系列特性図である。これは信号発生器で生成した2つの信号(q1、q2)の位相差を実際に観測したものである。短時間平均(約0.6秒)と長時間平均(約6秒)で挙動をみたところ、おおむね±1μラジアンの範囲に収まっていることが確認された。 図12は、上記実施例の周波数計測システム11による位相差分解能のばらつきを示す度数分布図である。図11の値から長時間平均データをヒストグラムにプロットしたものである。標準偏差σは1.085E−07であった。
【0091】
次に、周波数精度を検証した。図13は、上記周波数計測システム11における測定毎の計測値を平均化した際の周波数偏差散布図である。2チャンネルで周波数測定を行い、周波数のばらつきを分析したところ、1.5μHz(±3σ)であった。小数点以下の周波数は、μHzの桁まで一致することがわかった。周波数の平均値として440.998,958,17Hz(σ=0.52μHz)が得られた。
設定した周波数441.00Hzとの誤差はA/D変換器に搭載されるサンプリングクロック用水晶発振子の周波数偏差やドリフトが原因である。A/D変換器の安定化が十分でない場合、基準となるサンプリング周波数がまずドリフトする。そのため、本実装例における両チャンネルの瞬時周波数も同一方向にドリフトし、相関係数は1に近い値(強い相関)となる。さまざまな安定化の注意を払い、最終的に安定化後の相関係数は0.46であった。このことは、両チャンネルの誤差の一部が独立したプロセスによって発生していることを示唆している。
【0092】
図14は、微細周波数変動による位相変化を示す図である。 ある線路長おける伝送遅延による位相差の検出を試みる。光速は3.0e+8m/sec.と有限であり, わずかに伝送遅延を発生させる。
試験に用いた周波数441.00Hzの波長λは680272mである。もし1μラジアンの位相差が検出可能ならば、線路長の約10cmの差を検出できる。実験では長さが57cm異なる2本のケーブルを用いてファンクションジェネレータとA/D変換器を接続し、位相差φ1を計測した。次にケーブルを入れ替えて接続し位相差φ2を計測した。φ1=52.5μラジアン、φ2=810μラジアンの結果が得られた。位相差Δφ=757μラジアンとなり、相対的線路長差(=1.14m)における理論値である10.5μラジアンをはるかに上回る。この結果は光速の限界による伝送遅延に由来する位相差ではなく、おそらく線路の持つインピーダンスによる伝送遅延の位相差を捉えている。実験ではケーブル長の差や遅延量も十分に考慮しなければならないことを改めて示している。
【0093】
以上、本発明が適用できる計測システムとしては、周波数カウンタに代わる周波数計測器、信号発生器に内蔵される位相差計測処理部、位相差計測を利用するインピーダンスメータ、音響機器の周波数特性のうち特に位相遅延特性の計測、アレーマイクロホンを使用する信号処理における位相差計測、レーザー光線を使用した距離計測器における位相差計測を用いる時間遅延測定、音響周波数帯における微小なドップラー効果を利用する計測器、歌声や楽器の精密音程判別、音声処理における基本周波数推定、ピアノの調律で使用される精密チューナー、ピアノ調律におけるインハーモニシティの精密計測など、多彩な用途がある。
特に、低周波では1波長〜数波長のデータでも小数点以下の周波数が短時間で得られることから、電力測定器における周波数測定・力率測定、電力系統周波数(50Hz、60Hz)の短時間精密計測、電力系統における特定地点の位相差を精密に測定することで得られる電力潮流推定などに適用可能である。
【符号の説明】
【0094】
1 位相差計測システム、
2 計測プログラム、
22 デジタルフィルタプログラム、
23 位相差計測プログラム、
24 周波数計測プログラム、
3 コンピュータ、
4 CPU、
7 A/D変換器、
11 周波数計測システム、
q1,q2,・・・,qk 被計測信号(k=正の整数)
【技術分野】
【0001】
本発明は、リサージュ図となる閉曲線から被計測信号の位相差を求める位相差計測方法と位相差計測システム、並びに、リサージュ図となる閉曲線から被計測信号の周波数を求める周波数計測方法及び周波数計測システムに関する。
【背景技術】
【0002】
信号処理の分野において基準周波数との位相差(位相角)を計測することは重要なコア技術である。例えば、電力・電力量の位相差(位相角)測定、設備や機器の振動分析による異常診断や、静音性の評価や騒音原因の対策方法を検討するための騒音分析など様々な分野、更には音響処理における音源推定や楽音処理における倍音の分析において位相差の計測は重要である。
【0003】
従来、位相差計測方法としては、リサージュ法(Lissajous method)によりリサージュ図を描画する方法や、FFT(Fast Fourier Transform)アナライザを用いた方法や、波形から周期を測定する方法(電力・電力量の測定において使用されている)や、波形の部分間の相関を利用する方法(音声処理系で使用されている)などが知られている。
【0004】
既知の位相差計測方法のうち、周波数校正においては、微小な周波数偏差を検出するための簡便な方法としてリサージュ法が用いられることがある。
リサージュ図は、計測表示装置をX−Y入力モードに設定し、二つの単振動波形の順序対として得られる点の軌跡が描く平面図形である。周波数比が1の場合、位相差に応じて、形状は楕円・円・直線となり、周波数比が整数の場合は複数の山が描かれる。
【0005】
リサージュ図から2つの被計測信号の位相差を求める従来の方法としては、例えば、計測表示装置のX軸とY軸に、それぞれ周波数の等しい正弦波を入力すると、楕円形のリサージュ図が現れる。このときの位相差を正確に計算するためには図形がX軸(又はY軸)を横切る点をそれぞれ正確に求めることとなる。そして、リサージュ図を描く閉曲線のY軸方向の最大振幅と、前記閉曲線がY軸を横切る際のY軸方向の振幅とから位相差を計算する。
【0006】
特許文献1には、多軸加振試験機の位相差制御方法に関し、複数の系の信号を用いて表示装置上に正規化したリサージュ図を重ねて描かせることで入力信号の位相差と出力信号の位相差が揃っているか否かを検出することが記載されている。
【0007】
特許文献2には、光干渉計の周期誤差低減方法に関し、極座標変換した直交2相信号の位相角と半径の変動値を、誤差を求める各式に代入して前記直交2相信号の振幅、零点及び位相差の補正値を求め、極座標の半径が一定となるように補正することが記載されている。そして、光干渉計の動作原理を説明するために、リサージュ図が説明されている。
【先行技術文献】
【特許文献】
【0008】
【特許文献1】特開2000−298149号公報
【特許文献2】特許第4465451号公報
【発明の概要】
【発明が解決しようとする課題】
【0009】
しかしながら、上記従来のリサージュ法では、位相差を正確に計算するためには図形がX軸(又はY軸)を横切る点の座標をそれぞれ正確に求めなければならない。デジタル信号処理によって補間処理を行うことでゼロクロス点を求めることができるが、軌跡が原点に近づくにつれ原点近傍を通過する相対的な速度が大きくなり補間誤差が低下してしまう。
そして、従来、周波数カウンタでは「レシプロカル方式」を採用しており、小数点以下の周波数の計測が可能となっているが、周期を精密に計測するためには非常に高いサンプリング周波数を使用し周期を計測する必要がある。また周期の基準となるための「ゼロクロス点」の検出精度はDCオフセットの影響を受けてしまい検出精度が低下する懸念がある。
また、上記従来のリサージュ法では、微小な周波数偏差、例えば0.1/1000Hzを認識するためには360/10000(度/秒)で変化する位相差を短時間に認識し周波数の微調整を繰り返す必要がある。
【0010】
そこで、本願発明の目的は、従来のリサージュ法による処理の仕方を改良することで、精度良く位相差計測を行なうことが可能な位相差計測方法及び位相差計測システムと、精度良く周波数計測を行なうことが可能な周波数計測方法及び周波数計測システムを提供することにある。
【課題を解決するための手段】
【0011】
本願発明者は、リサージュ図が描く閉曲線の囲む面積に着目し、短時間で精密な位相差を検出できる方法を見出し、本願発明を完成するに至った。
【0012】
本発明の位相差計測方法は、2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする。
【0013】
従来のリサージュ法では、リサージュ図のうちのゼロクロスする近傍のデータしか計算に使用されていないが、本発明によれば、リサージュ図の閉曲線の囲む面積を算出することで、リサージュ図を構成する全データが計算に使用されることとなり、精密な位相差計測ができる。
【0014】
本発明の位相差計測方法は、2つの被計測信号をそれぞれA/D変換するデジタル変換ステップと、デジタル変換後の2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする。
【0015】
本発明によれば、ゼロクロスのような特定の条件を抽出する処理が存在せず、アナログ信号をデジタル信号に変換する変換手段(A/D変換)における誤差は、全体に正負が統計的に等しい状態で拡散され、誤差が大幅に低減されることとなる。
【0016】
本発明の位相差計測システムは、2つの被計測信号をそれぞれA/D変換するA/D変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出し、前記2つの被計測信号それぞれの振幅値から前記面積を正規化し、正規化された前記面積から前記2つの被計測信号の位相差を算出することを特徴とする。
【0017】
本発明によれば、リサージュ図を描画可能な1周期の被計測信号波形があれば、コンピュータ上で動作する前記計測プログラムによって前記2つの被計測信号の位相差が計算できることから、短時間で位相差計測を行なうことができる。
【0018】
本発明の周波数計測方法は、被計測信号をフィルタリングし、このフィルタリング前後での信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後の2つの信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つのデジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする。
【0019】
本発明によれば、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、前記周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を正確に求めることができる。
【0020】
本発明の周波数計測方法は、被計測信号をA/D変換するデジタル変換ステップと、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記デジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする。
【0021】
本発明の周波数計測システムは、被計測信号をA/D変換するA/D変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出し、前記被計測信号の振幅値から前記面積を正規化し、正規化された前記面積から前記デジタル信号の位相差を算出し、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定することを特徴とする。
【0022】
本発明によれば、ゼロクロスのような特定の条件を抽出する処理が存在せず、アナログ信号をデジタル信号に変換する変換手段(A/D変換)における誤差は、全体に正負が統計的に等しい状態で拡散され、誤差が大幅に低減されることとなる。
【0023】
本発明は、前記フィルタリングに、IIRくし型フィルタを用いることを特徴とする。
【0024】
本発明によれば、IIRくし型フィルタ(フィードバック型くし型フィルタとも呼ばれる)を用いることで、設定されたノッチ周波数近傍で入出力間の位相を急峻に変化させるので、周波数の最小分解能が向上する。
【0025】
本発明は、前記フィルタリングに、線形位相FIRフィルタを用いることを特徴とする。
【0026】
本発明によれば、線形位相FIRフィルタを用いることで、直線性が良くなり、計測可能な周波数帯域が広帯域となる。
【発明の効果】
【0027】
本発明の位相差計測方法と位相差計測システムによれば、リサージュ図の閉曲線の囲む面積を算出することで、リサージュ図を構成する全データが計算に使用されることとなり、精密な位相差計測ができる。そして、ゼロクロスのような特定の条件を抽出する処理が存在せず、アナログ信号をデジタル信号に変換する変換手段(A/D変換)における誤差は、全体に正負が統計的に等しい状態で拡散され、誤差が大幅に低減されることとなる。また、本発明によれば、閉曲線積分の性質から原点は軌跡の外であってもかまわない。信号にオフセットが加わっても面積の計算には反映されない。ノイズが加わっても面積は平均化され影響は緩和される。その結果、しきい値周辺のデータに依存するゼロクロス法の欠点が解消される。そして、本発明によれば、リサージュ図を描画可能な1周期の被計測信号波形があれば、コンピュータ上で動作する前記計測プログラムによって前記2つの被計測信号の位相差が計算できることから、短時間で位相差計測を行なうことができる。
【0028】
本発明の周波数計測方法と周波数計測システムによれば、前記位相差と前記フィルタリングにおける移相特性から周波数偏差を算出し、前記周波数偏差と前記フィルタリングにおけるノッチ周波数とから前記被計測信号の周波数を正確に求めることができる。そして、前記フィルタリングに、IIRくし型フィルタを用いることで、設定されたノッチ周波数近傍で入出力間の位相を急峻に変化させるので、周波数の最小分解能が向上する。または、線形位相FIRフィルタを用いることで、直線性が良くなり、計測可能な周波数帯域が広帯域となる。
【図面の簡単な説明】
【0029】
【図1】本発明の実施の形態の位相差計測システムを説明するブロック図である。
【図2】上記位相差計測システムに係る計測プログラムの構成を例示するブロック図である。
【図3】リサージュ図から2つの信号の位相差を求める本発明の方法を示す図である。
【図4】本発明の実施の形態の位相差計測手順を示すフローチャート図である。
【図5】本発明の実施の形態の周波数計測システムを説明するブロック図である。
【図6】上記周波数計測システムに係る計測プログラム構成を例示するブロック図である。
【図7】上記周波数計測システムに係る計測プログラム構成を例示するブロック図である。
【図8】本発明の実施の形態の周波数計測手順を示すフローチャート図である。
【図9】IIRくし型フィルタの構造図である。
【図10】IIRくし型フィルタの特性図であり、(a)は移相特性図であり、(b)は振幅特性図である。
【図11】上記周波数計測システムによる位相差分解能のばらつきを示す時系列特性図である。
【図12】上記周波数計測システムによる位相差分解能のばらつきを示す度数分布図である。
【図13】上記周波数計測システムにおける測定毎の計測値を平均化した際の周波数偏差散布図である。
【図14】微細周波数差によって発生する位相変化を示す図である。
【発明を実施するための形態】
【0030】
本発明に係る位相差計測原理について以下に説明する。
【0031】
(位相差計測原理)
信号xと信号yの周波数比が1:1,位相差φである場合のリサージュ図は、位相差及び振幅に応じて楕円/円/直線となる。これを閉曲線と捉えて、その軌道の描く面積から位相差sin(φ)を求める方法を以下に説明する。
【0032】
まず、閉曲線の囲む面積の一般式を数1に示す。
【0033】
【数1】
【0034】
図3に示すリサージュ図となる閉曲線の囲む面積を符号Sとすると、軌跡中、隣接する2点と原点がつくる微小な三角形の部分面積dSは、軌跡と軌跡の移動ベクトルとの外積を用いて、次の数式2または数式3で表される。この部分面積dSは、軌跡のベクトル(=P)と移動ベクトル(=dP)の外積をとり1/2にすることで得られる。
【0035】
【数2】
【0036】
【数3】
【0037】
部分面積dSを外積の計算式にしたがって計算すると、軌跡が反時計回りの時はdS>0、時計回りの時はdS<0となる。また, 全体の面積SはdSを総計し次の数式4で表される。周期Tは被計測信号の周波数とサンプリング条件によって変化する。
【0038】
【数4】
【0039】
ここでは「位相差φの信号が描くリサージュ面積」を導く。位相差φのある信号q1,q2を考える。符号a1,a2は振幅、符号ω0は角速度であるとすると、次の数式5で表される。
【0040】
【数5】
【0041】
三角関数の加法定理を用いて、数式5およびその微分形を数式1に代入・変形すると面積Sは次の数式6で表される。ここで、積分区間は1周期とする。描画される面積は位相差と振幅で決まることが示される。
【0042】
【数6】
【0043】
前述の振幅(a1,a2)が別の方法で計測されたと仮定し、得られた閉曲線面積Sを、次の数式7で表すように正規化し、正規化リサージュ面積(LsRot)を得る。
【0044】
【数7】
【0045】
上記数式7で得られた面積は、位相差=0度の時に面積=0、位相差=±90度の時に面積=πとなるよう正規化された面積である。−π≦LsRot≦πとなり無次元量となる。上記数式6と数式7より、次の数式8が得られる。
【0046】
【数8】
【0047】
上記数式8から、位相差sin(φ)は正規化リサージュ面積(LsRot)から容易に得られることがわかる。
【0048】
ここで、上述した、図3に示すリサージュ図となる閉曲線の囲む部分面積dSは、軌跡と軌跡の移動ベクトルとの外積を用いて、前述の数式2または数式3で表される。しかし、外積の性質から、P×P≡0(恒等的にゼロ)であり、前記差分をとる計算は不要であることから、数式9を導くことができる。
【0049】
【数9】
【0050】
本発明に係る位相差計測手順、位相差計測システム、及び位相差計測プログラムについて以下に説明する。
【0051】
(位相差計測手順)
図4は、本発明の実施形態の位相差計測手順を示すフローチャート図である。本実施形態の位相差計測方法は、2つの被計測信号q1,q2をそれぞれA/D変換するデジタル変換ステップS11と、デジタル変換後の2つの被計測信号qd1,qd2からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップS21と、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップS31と、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出ステップS41からなる。 本実施形態によれば、リサージュ図を描画可能な1周期の被計測信号波形があれば、コンピュータ3上で動作する計測プログラム2によって前記2つの被計測信号の位相差が計算でき、短時間で位相差計測を行なうことができる。
【0052】
(位相差計測システム)
図1は、本発明の実施の形態の位相差計測システム1を機能的に説明するブロック図である。本実施形態の位相差計測システム1は、2つの被計測信号q1、q2をそれぞれA/D変換するA/D変換器7と、計測プログラム2を備えたコンピュータ3から構成される(図1)。
【0053】
図1に示す例では、コンピュータ3には、データ入力を受け付ける入力手段81と、入力されたデータに演算処理を行うCPU4と、演算処理後のデータを蓄積保存するデータベース5とデータ処理前後のデータを画面表示するディスプレイ8が備わっており、コンピュータ3上で計測プログラム2が動作する。計測プログラム2は、デジタルフィルタプログラム22と位相差計測プログラム23からなる。
【0054】
(位相差計測プログラム)
図2は、上記実施形態の位相差計測システム1に係る計測プログラム2の構成を例示するブロック図である。図2に示す例では、計測プログラム2は、デジタルフィルタプログラム22と位相差計測プログラム23からなる。デジタル変換ステップS11にてデジタル変換された2つの被計測信号qd1,qd2がコンピュータ3に入力されると、それぞれデジタルフィルタ222を通して、リサージュ図となる閉曲線の囲む面積を算出する面積計算処理を行い(符号S21)、2つの被計測信号qd1,qd2のそれぞれの振幅計算処理を行い(符号S211)、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する面積正規化処理を行い(符号S31)、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出処理(符号S41)を行なった後、位相差信号データan12として出力する。なお、ここでは、リサージュ図を直接的に描画することはしていないが、例えば、面積計算処理S21の前段階でリサージュ図を描画してディスプレイ表示する構成としても良い。
【0055】
次に、本発明の周波数計測におけるフィルタについて説明する。
【0056】
(IIRくし型フィルタ)
一般に、IIRくし型フィルタの構造は、図9のように示される。遅延段数をK、フィードバックゲインをaとすると、伝達関数は次の数式10で表される。jは虚数単位である。
【0057】
【数10】
【0058】
周波数と角速度をサンプリング周波数fsで正規化すると、正規化周波数F、正規化角速度Ωは次の数式11で表される。
【0059】
【数11】
【0060】
前述の正規化角周波数と正規化角速度を用い、伝達関数は正規化された角速度を用いて次の数式12で表される。以後、この正規化した表現を基準に説明する。
【0061】
【数12】
【0062】
IIRくし型フィルタの移相φは前述の伝達関数から、次の数式13で表される。次に、移相特性φがΩの関数として表せるように数13を変形し、周波数Ωによってどのくらいの位相変化φが発生するかを、フィルタ移相特性φ(Ω)で表す。フィルタ移相特性φ(Ω)は、次の数式14で表される。
【0063】
【数13】
【0064】
【数14】
【0065】
図10は、IIRくし型フィルタの特性図である。図10(a)は移相特性図であり、グラフの縦軸は入出力位相差を示しており、グラフの横軸は相対周波数を示している。図10(b)は振幅特性図であり、グラフの縦軸は振幅を示しており、グラフの横軸は相対周波数を示している。ここでは−1<a<0を使用しており、ノッチ周波数の奇数倍の周波数において極めて鋭敏な移相特性とノッチ周波数近傍のみを通過させる振幅特性となっている。
【0066】
上記数式14をテーラー展開して得られる3次項以上を切り捨て、1次近似する。ノッチ周波数近傍での移相φの近似式は、次の数式15で表される。
【0067】
【数15】
【0068】
ここで、位相変化率を求める。f0近傍における位相変化率は、φをΩで微分することで得られ次の数式16で表される。
【0069】
【数16】
【0070】
次に、得られた位相変化率から周波数偏差を求める。Δfはf0からの変位とする。ΔφとΔfの関係は、次の数式17で表される。
【0071】
【数17】
【0072】
上記数式17から、計測される位相差は周波数偏差がπg0倍されて計測できることがわかる(感度の向上)。上記数式17を変形すると、次の数式18のようにも書ける。
【0073】
【数18】
【0074】
上記数式18を移項し周波数偏差(Δf/f0)は、次の数式19で表される。
【0075】
【数19】
【0076】
くし形フィルタのノッチ周波数f0はサンプリング周波数fsと遅延段数のみで決定できる(数式21を参照)。 前述の正規化リサージュ面積を精密に定量したとすると、数式19と周波数偏差Δfとf0の関係を合わせ、数式20の式が得られる。数式20,数式21のように単純な式を用いて精密な周波数判別が可能である。
【0077】
【数20】
【0078】
【数21】
【0079】
得られた数式20,数式21を見るとサンプリング周波数が基準の周波数となっている(国家標準を用いて校正可能である)。トレーサブルな周波数標準を計測することでサンプリング周波数のキャリブレーションが可能である。つまり、得られた数式20,数式21を見ると、K,g0π2は純粋な定数である。 正規化リサージュ面積(LsRot)は計測値である。fsはハードウェアに実装される水晶発振子の精度により実際の装置で動作している周波数の実測値が決まる。 よってfsを正しく校正することで(国家標準を用いて校正可能である)、実際に動作しているサンプリング周波数(fs)の確度を高めることができ、計算結果(理論値)と実際の値(計測値)を近づけることができる。
【0080】
本発明に係る周波数計測手順、周波数計測システム、及び周波数計測プログラムについて以下に説明する。
【0081】
(周波数計測手順)
図8は、本発明の実施形態の周波数計測手順を示すフローチャート図である。本実施形態の周波数計測方法は、被計測信号(例えばq1)をA/D変換するデジタル変換ステップS11と、デジタル変換後の被計測信号(例えばqd1)をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力する位相差生成ステップS15と、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップS21と、前記被計測信号の振幅値から前記面積を正規化する正規化ステップS31と、正規化された前記面積から前記デジタル信号の位相差を算出する位相差算出ステップS41と、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップS51からなる。
【0082】
(周波数計測システム)
図5は、本発明の実施形態の周波数計測システム11を説明するブロック図である。本実施形態の周波数計測システム11は、被計測信号q1,q2,・・・,qk(kは正の整数)をそれぞれA/D変換するA/D変換器7と、計測プログラム2を備えたコンピュータ3から構成される(図5)。ここで、同一の符号は同じ機能を表しており、その説明を適宜省略する。図5に示す例では、多チャンネルの被計測信号を入力する場合で説明しているが、一つの被計測信号(例えばq1)としても良いし、二つの被計測信号(例えばq1,q2)としても良い。
【0083】
(第1の周波数計測プログラム)
図6は、上記実施形態の周波数計測システム11に係る計測プログラム2の構成を例示するブロック図である。図6は、一つの被計測信号q1を入力する場合の例であり、計測プログラム2は、デジタルフィルタプログラム22と周波数計測プログラム24からなる。デジタル変換ステップS11にてデジタル変換された被計測信号qd1がコンピュータ3に入力されると、デジタルフィルタ222を通し、このフィルタリング前後でのデジタル信号の位相差を変換出力して位相差を生成し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積計算処理を行い(符号S21)、被計測信号の振幅計算処理を行い(符号S211)、前記被計測信号の振幅値から前記面積を正規化する面積正規化処理を行い(符号S31)、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出処理(符号S41)を行い、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定処理を行ない(符号S51)、その結果を周波数信号データan1として出力する。なお、ここでは、リサージュ図を直接的に描画することはしていないが、例えば、面積計算処理S21の前段階でリサージュ図を描画してディスプレイ表示する構成としても良い。
【0084】
本実施形態の周波数計測システム11では、前記デジタルフィルタとしてIIRくし型フィルタを用いている。これは、IIRくし型フィルタを用いることで、設定されたノッチ周波数近傍で入出力間の位相を急峻に変化させて、周波数の最小分解能を向上させるためである。
【0085】
(第2の周波数計測プログラム)
図7は、上記実施形態の周波数計測システム11に係る計測プログラム2の構成の他の例を示すブロック図である。図7に示す例では、計測プログラム2は、デジタルフィルタプログラム22と周波数計測プログラム24からなる。図7は、2チャンネルの被計測信号q1,q2を入力する場合の例である。デジタル変換ステップS11にてデジタル変換された被計測信号qd1、qd2がそれぞれコンピュータ3に入力されると、前処理フィルタ221を通してデジタルフィルタ222(くし型フィルタ222)を通す。このくし型フィルタ222によるフィルタリング前後でのデジタル信号の位相差を変換出力して位相差を生成し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積計算処理を行い(符号S21)、被計測信号の振幅計算処理を行い(符号S211)、前記被計測信号の振幅値から前記面積を正規化する面積正規化処理を行い(符号S31)、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出処理(符号S41)を行い、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定処理を行ない(符号S51)、その結果をそれぞれ周波数信号データan1,an2として出力する。そして、それぞれデジタルフィルタ222を通されたデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積計算処理を行い(符号S21)、2つの被計測信号qd1,qd2のそれぞれの振幅計算処理を行い(符号S211)、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する面積正規化処理を行い(符号S31)、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出処理(符号S41)を行なった後、位相差信号データan12として出力する。
【0086】
(実施例)
図5と図7に示すブロック図の周波数計測システムを試作した。2チャンネルの入力をもち、それぞれの精密周波数とチャンネル間の精密位相差を計測する構成である。
【0087】
前処理フィルタ221は目的外の周波数帯域をカットするために用いた。くし型フィルタ222は、特定のノッチ周波数が遅延段数Kによって設定されるIIRくし型フィルタを用いた。周波数計測システム11は信号周波数f1,f2を計測し、それと同時に2つの信号の位相差を得る。サンプリングは16bit,44100Hzとした。
【0088】
次に、上記実施例の周波数計測システム11を用いて最小位相差分解能・分解能ばらつき・精密周波数の精度について評価を行なった。 試験信号としてルビジウム(Rb)発振器(10MHz)を外部周波数基準としたファンクションジェネレータ(以下FGと略す)を使用した。信号には2チャンネル44100Hzを用いた。A/D変換器7を外部接続して使用した。
【0089】
今回、意図的に位相差ゼロの状態の信号を作成し本ソフトウェアによる計測を行った。測定は10分間行った。数値の平均と標準偏差から分解能を推定した。短時間平均(約0.6秒)で1.1μラジアン(±3σ)、長時間平均(約6秒)で0.33μラジアン(±3σ)の分解能が確認できた。
【0090】
図11は、上記周波数計測システム11による位相差分解能のばらつきを示す時系列特性図である。これは信号発生器で生成した2つの信号(q1、q2)の位相差を実際に観測したものである。短時間平均(約0.6秒)と長時間平均(約6秒)で挙動をみたところ、おおむね±1μラジアンの範囲に収まっていることが確認された。 図12は、上記実施例の周波数計測システム11による位相差分解能のばらつきを示す度数分布図である。図11の値から長時間平均データをヒストグラムにプロットしたものである。標準偏差σは1.085E−07であった。
【0091】
次に、周波数精度を検証した。図13は、上記周波数計測システム11における測定毎の計測値を平均化した際の周波数偏差散布図である。2チャンネルで周波数測定を行い、周波数のばらつきを分析したところ、1.5μHz(±3σ)であった。小数点以下の周波数は、μHzの桁まで一致することがわかった。周波数の平均値として440.998,958,17Hz(σ=0.52μHz)が得られた。
設定した周波数441.00Hzとの誤差はA/D変換器に搭載されるサンプリングクロック用水晶発振子の周波数偏差やドリフトが原因である。A/D変換器の安定化が十分でない場合、基準となるサンプリング周波数がまずドリフトする。そのため、本実装例における両チャンネルの瞬時周波数も同一方向にドリフトし、相関係数は1に近い値(強い相関)となる。さまざまな安定化の注意を払い、最終的に安定化後の相関係数は0.46であった。このことは、両チャンネルの誤差の一部が独立したプロセスによって発生していることを示唆している。
【0092】
図14は、微細周波数変動による位相変化を示す図である。 ある線路長おける伝送遅延による位相差の検出を試みる。光速は3.0e+8m/sec.と有限であり, わずかに伝送遅延を発生させる。
試験に用いた周波数441.00Hzの波長λは680272mである。もし1μラジアンの位相差が検出可能ならば、線路長の約10cmの差を検出できる。実験では長さが57cm異なる2本のケーブルを用いてファンクションジェネレータとA/D変換器を接続し、位相差φ1を計測した。次にケーブルを入れ替えて接続し位相差φ2を計測した。φ1=52.5μラジアン、φ2=810μラジアンの結果が得られた。位相差Δφ=757μラジアンとなり、相対的線路長差(=1.14m)における理論値である10.5μラジアンをはるかに上回る。この結果は光速の限界による伝送遅延に由来する位相差ではなく、おそらく線路の持つインピーダンスによる伝送遅延の位相差を捉えている。実験ではケーブル長の差や遅延量も十分に考慮しなければならないことを改めて示している。
【0093】
以上、本発明が適用できる計測システムとしては、周波数カウンタに代わる周波数計測器、信号発生器に内蔵される位相差計測処理部、位相差計測を利用するインピーダンスメータ、音響機器の周波数特性のうち特に位相遅延特性の計測、アレーマイクロホンを使用する信号処理における位相差計測、レーザー光線を使用した距離計測器における位相差計測を用いる時間遅延測定、音響周波数帯における微小なドップラー効果を利用する計測器、歌声や楽器の精密音程判別、音声処理における基本周波数推定、ピアノの調律で使用される精密チューナー、ピアノ調律におけるインハーモニシティの精密計測など、多彩な用途がある。
特に、低周波では1波長〜数波長のデータでも小数点以下の周波数が短時間で得られることから、電力測定器における周波数測定・力率測定、電力系統周波数(50Hz、60Hz)の短時間精密計測、電力系統における特定地点の位相差を精密に測定することで得られる電力潮流推定などに適用可能である。
【符号の説明】
【0094】
1 位相差計測システム、
2 計測プログラム、
22 デジタルフィルタプログラム、
23 位相差計測プログラム、
24 周波数計測プログラム、
3 コンピュータ、
4 CPU、
7 A/D変換器、
11 周波数計測システム、
q1,q2,・・・,qk 被計測信号(k=正の整数)
【特許請求の範囲】
【請求項1】
2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする位相差計測方法。
【請求項2】
2つの被計測信号をそれぞれA/D変換するデジタル変換ステップと、デジタル変換後の2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする位相差計測方法。
【請求項3】
被計測信号をフィルタリングし、このフィルタリング前後での信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後の2つの信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つのデジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする周波数計測方法。
【請求項4】
被計測信号をA/D変換するデジタル変換ステップと、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記デジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする周波数計測方法。
【請求項5】
2つの被計測信号をそれぞれA/D変換するA/D変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出し、前記2つの被計測信号それぞれの振幅値から前記面積を正規化し、正規化された前記面積から前記2つの被計測信号の位相差を算出することを特徴とする位相差計測システム。
【請求項6】
被計測信号をA/D変換するA/D変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出し、前記被計測信号の振幅値から前記面積を正規化し、正規化された前記面積から前記デジタル信号の位相差を算出し、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定することを特徴とする周波数計測システム。
【請求項7】
前記フィルタリングに、IIRくし型フィルタを用いることを特徴とする請求項6記載の周波数計測システム。
【請求項8】
前記フィルタリングに、線形位相FIRフィルタを用いることを特徴とする請求項6記載の周波数計測システム。
【請求項1】
2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする位相差計測方法。
【請求項2】
2つの被計測信号をそれぞれA/D変換するデジタル変換ステップと、デジタル変換後の2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記2つの被計測信号それぞれの振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つの被計測信号の位相差を算出する位相差算出ステップを有することを特徴とする位相差計測方法。
【請求項3】
被計測信号をフィルタリングし、このフィルタリング前後での信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後の2つの信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記2つのデジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする周波数計測方法。
【請求項4】
被計測信号をA/D変換するデジタル変換ステップと、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力する位相差生成ステップと、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出する面積算出ステップと、前記被計測信号の振幅値から前記面積を正規化する正規化ステップと、正規化された前記面積から前記デジタル信号の位相差を算出する位相差算出ステップと、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定する周波数推定ステップを有することを特徴とする周波数計測方法。
【請求項5】
2つの被計測信号をそれぞれA/D変換するA/D変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の2つの被計測信号からリサージュ図となる閉曲線の囲む面積を算出し、前記2つの被計測信号それぞれの振幅値から前記面積を正規化し、正規化された前記面積から前記2つの被計測信号の位相差を算出することを特徴とする位相差計測システム。
【請求項6】
被計測信号をA/D変換するA/D変換器と、計測プログラムを備えたコンピュータから構成され、前記計測プログラムが、デジタル変換後の被計測信号をフィルタリングし、このフィルタリング前後でのデジタル信号の位相差を変換出力し、前記フィルタリング前後のデジタル信号からリサージュ図となる閉曲線の囲む面積を算出し、前記被計測信号の振幅値から前記面積を正規化し、正規化された前記面積から前記デジタル信号の位相差を算出し、前記位相差と前記フィルタリングにおけるフィルタ移相特性から周波数偏差を算出し、この周波数偏差と前記フィルタリングのノッチ周波数とから前記被計測信号の周波数を推定することを特徴とする周波数計測システム。
【請求項7】
前記フィルタリングに、IIRくし型フィルタを用いることを特徴とする請求項6記載の周波数計測システム。
【請求項8】
前記フィルタリングに、線形位相FIRフィルタを用いることを特徴とする請求項6記載の周波数計測システム。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【公開番号】特開2013−53929(P2013−53929A)
【公開日】平成25年3月21日(2013.3.21)
【国際特許分類】
【出願番号】特願2011−192346(P2011−192346)
【出願日】平成23年9月5日(2011.9.5)
【新規性喪失の例外の表示】特許法第30条第1項適用申請有り 平成23年3月3日 社団法人電子情報通信学会発行の「電子情報通信学会技術研究報告(信学技報,vol.110,no.465)に発表。 平成23年7月29日 社団法人電気学会発行の「電気学会研究会資料 計測研究会講演要旨集」をもって発表。
【出願人】(511216271)
【Fターム(参考)】
【公開日】平成25年3月21日(2013.3.21)
【国際特許分類】
【出願日】平成23年9月5日(2011.9.5)
【新規性喪失の例外の表示】特許法第30条第1項適用申請有り 平成23年3月3日 社団法人電子情報通信学会発行の「電子情報通信学会技術研究報告(信学技報,vol.110,no.465)に発表。 平成23年7月29日 社団法人電気学会発行の「電気学会研究会資料 計測研究会講演要旨集」をもって発表。
【出願人】(511216271)
【Fターム(参考)】
[ Back to top ]