JP5286524B2

JP5286524B2 - 周波数測定装置及び周波数測定方法

Info

Publication number: JP5286524B2
Application number: JP2008510844A
Authority: JP
Inventors: 雄作藤井
Original assignee: Gunma University NUC
Current assignee: Gunma University NUC
Priority date: 2006-03-22
Filing date: 2007-03-22
Publication date: 2013-09-11
Anticipated expiration: 2027-03-22
Also published as: JPWO2007119488A1; WO2007119488A1

Description

本発明は、周波数測定装置及び周波数測定方法に係り、特に、入力波形の周波数又は周期を測定する周波数測定装置及び周波数測定方法に関する。

従来より、周波数カウンタ等による周波数測定においては、入力波形のゼロクロス点（時刻）を測定し、複数（例えば、Ｎ＋１個）のゼロクロスを含むＮ個の周期の時間Ｐを計算し、その時間Ｐにおける平均周期Ｔ（＝Ｐ／Ｎ）、及び平均周波数ｆ（＝Ｎ／Ｐ＝１／Ｔ）を計算している。

しかしながら、従来の方法では、サンプリング間隔内に存在する複数のゼロクロスは、単に、波の数を数える目的にのみ使用されおり、時刻情報など使用していない情報が多いため、周波数測定における精度や分解能が低下してしまう、という問題がある。

また、周波数測定の分解能を向上させるためには、サンプリング間隔を定めるＮを大きくするしかなく、サンプリング間隔が拡大してしまう、という問題がある。

また、周波数の測定において全波形を記録した後に信号処理を行う方法として、自己相関法などがあるが、信号処理に膨大な量の計算を要するため、実現困難な場合が多い、という問題がある。

本発明は、上記の問題点を解決するためになされたもので、サンプリング間隔を拡大させることなく、高精度かつ高分解能で周波数又は周期を測定することができる周波数測定装置及び周波数測定方法を提供することを目的とする。

上記の目的を達成するために第１の発明に係る周波数測定装置は、入力された波形情報の第１のゼロクロス位置を含む第１の所定範囲内の複数のゼロクロス位置各々の時刻の平均値を、前記第１のゼロクロス位置に対する第１の時刻として算出する第１の算出手段と、前記波形情報の前記第１のゼロクロス位置とは異なる第２のゼロクロス位置を含む第２の所定範囲内の複数のゼロクロス位置各々の時刻の平均値を、前記第２のゼロクロス位置に対する第２の時刻として算出する第２の算出手段と、前記第１の時刻、前記第２の時刻、及び前記第１のゼロクロス位置と前記第２のゼロクロス位置との間に存在する波の数に基づいて、前記波形情報の周期及び周波数の少なくとも一方を測定する測定手段とを含んで構成されている。

また、第２の発明に係る周波数測定方法は、入力された波形情報の第１のゼロクロス位置を含む第１の所定範囲内の複数のゼロクロス位置各々の時刻の平均値を、前記第１のゼロクロス位置に対する第１の時刻として算出し、前記波形情報の前記第１のゼロクロス位置とは異なる第２のゼロクロス位置を含む第２の所定範囲内の複数のゼロクロス位置各々の時刻の平均値を、前記第２のゼロクロス位置に対する第２の時刻として算出し、前記第１の時刻、前記第２の時刻、及び前記第１のゼロクロス位置と前記第２のゼロクロス位置との間に存在する波の数に基づいて、前記波形情報の周期及び周波数の少なくとも一方を測定することを特徴としている。

第１の発明及び第２の発明によれば、入力された波形情報の第１のゼロクロス位置を含む第１の所定範囲内の複数のゼロクロス位置各々の時刻の平均値を、第１のゼロクロス位置に対する第１の時刻として算出し、また、波形情報の第１のゼロクロス位置とは異なる第２のゼロクロス位置を含む第２の所定範囲内の複数のゼロクロス位置各々の時刻の平均値を、第２のゼロクロス位置に対する第２の時刻として算出する。

そして、第１の時刻、第２の時刻、及び第１のゼロクロス位置と第２のゼロクロス位置との間に存在する波の数に基づいて、波形情報の周期及び周波数の少なくとも一方を測定する。

従って、時刻計測対象となるゼロクロス位置を含む所定範囲内の複数のゼロクロス位置各々の時刻の平均値を用いて、波形情報の周期又は周波数を測定することにより、多くのゼロクロス位置の時刻を用いて測定しているため、サンプリング間隔を拡大させることなく、高精度かつ高分解能で周波数又は周期を測定することができる。

なお、ゼロクロス位置とは、波形情報が振動を示している場合に、振動方向が一方向から他方向へ反転する位置をいう。

また、第３の発明に係る周波数測定装置は、入力された波形情報の第１のゼロクロス位置を含む第１の所定範囲内の複数のゼロクロス位置の各々と、前記第１のゼロクロス位置とは異なる第２のゼロクロス位置を含む第２の所定範囲内の複数のゼロクロス位置の各々との時刻の差を各々算出する算出手段と、前記算出手段によって算出された時刻の差の平均値及び前記第１のゼロクロス位置と前記第２のゼロクロス位置との間に存在する波の数に基づいて、前記波形情報の周期及び周波数の少なくとも一方を測定する測定手段とを含んで構成されている。

また、第４の発明に係る周波数測定方法は、入力された波形情報の第１のゼロクロス位置を含む第１の所定範囲内の複数のゼロクロス位置の各々と、前記第１のゼロクロス位置とは異なる第２のゼロクロス位置を含む第２の所定範囲内の複数のゼロクロス位置の各々との時刻の差を各々算出し、前記算出された時刻の差の平均値及び前記第１のゼロクロス位置と前記第２のゼロクロス位置との間に存在する波の数に基づいて、前記波形情報の周期及び周波数の少なくとも一方を測定することを特徴としている。

第３の発明及び第４の発明によれば、入力された波形情報の第１のゼロクロス位置を含む第１の所定範囲内の複数のゼロクロス位置の各々と、第１のゼロクロス位置とは異なる第２のゼロクロス位置を含む第２の所定範囲内の複数のゼロクロス位置の各々との時刻の差を各々算出する。

そして、算出された時刻の差の平均値及び第１のゼロクロス位置と第２のゼロクロス位置との間に存在する波の数に基づいて、波形情報の周期及び周波数の少なくとも一方を測定する。

従って、時刻計測対象となるゼロクロス位置の各々を含む各所定範囲内の複数のゼロクロス位置各々の時刻に基づいて算出される複数の時刻の差の平均値を用いて、波形情報の周期又は周波数を測定することにより、多くのゼロクロス位置の時刻を用いて測定しているため、サンプリング間隔を拡大させることなく、高精度かつ高分解能で周波数又は周期を測定することができる。

また、上記の第１の所定範囲を、第１のゼロクロス位置の前後同数のゼロクロス位置を含む範囲とし、第２の所定範囲を、第２のゼロクロス位置の前後同数のゼロクロス位置を含む範囲とすることができる。これにより、時刻計測対象となるゼロクロス位置の前後同数のゼロクロス位置の時刻情報を有効活用して、サンプリング間隔を拡大させることなく、高精度かつ高分解能で周波数又は周期を測定することができる。

以上説明したように、本発明の周波数測定装置及び周波数測定方法によれば、時刻計測対象となるゼロクロス位置を含む所定範囲内の複数のゼロクロス位置各々の時刻の平均値、又は時刻計測対象となるゼロクロス位置の各々を含む各所定範囲内の複数のゼロクロス位置各々の時刻に基づいて算出される複数の時刻の差の平均値を用いて、波形情報の周期又は周波数を測定することにより、多くのゼロクロス位置の時刻を用いて測定しているため、サンプリング間隔を拡大させることなく、高精度かつ高分解能で周波数又は周期を測定することができる、という効果が得られる。

本発明の第１の実施の形態に係る周波数測定装置を示すブロック図である。本発明の第１の実施の形態の周波数を計測する計測ルーチンの内容を示す流れ図である。本発明の第１の実施の形態の波形情報を示すグラフである。ビート周波数及び基準周波数の測定データを示すグラフである。ビート周波数及び基準周波数の測定データから算出した錘体の変位速度を示すグラフである。変位速度から算出した錘体の変位を示すグラフである。変位速度から算出した錘体の加速度を示すグラフである。加速度から算出した錘体の慣性力を示すグラフである。ビート周波数及び基準周波数の瞬時値に基づく慣性力とビート周波数の瞬時値及び基準周波数の平均値に基づく慣性力とを示すグラフである。デジタイザを用いて計測した慣性力と周波数カウンタを用いて計測した慣性力とを示すグラフである。慣性力と位置との関係を示すグラフである。通常の周波数測定アルゴリズムで測定した慣性力を示すグラフである。通常の周波数測定アルゴリズムにおいて、サンプリング数を大きくとって測定した慣性力を示すグラフである。本発明の第１の実施の形態に係る周波数測定アルゴリズムで測定した慣性力を示すグラフである。本発明の第３の実施の形態の波形情報を示すグラフである。

以下、図面を参照して本発明の実施の形態を詳細に説明する。本実施の形態では、浮上質量法に基づいて、衝撃力に対する材料の力学的応答を求めるために周波数を高精度に測定する周波数測定装置に、本発明を適用した場合を例に説明する。

ここで、浮上質量法は、浮上支持した質量（慣性質量）に作用する慣性力を光波干渉計を利用して高精度に取り出すことを特徴とする、変動する力の発生方法及び計測方法である。浮上質量法では、浮上支持した物体（慣性質量）に作用する慣性力を、光波干渉計により、超高精度に測定する。測定中に光波干渉計で測定するのは、物体にあてた信号光のドップラーシフト周波数ｆのみである。このドップラーシフト周波数ｆから、物体の速度ｖ、位置ｘ、加速度ａ、及び慣性力Ｆを数値微分や数値積分などによって求める。これにより、物体の速度ｖ、位置ｘ、加速度ａ、及び慣性力Ｆを、高精度に、かつ、完全に同期したデータとして、測定することができる。

また、本実施の形態に係る周波数測定装置は、ＰＣベースの計測システム（ＮａｔｉｏｎａｌＩｎｓｔｒｕｍｅｎｔｓのＬａｂＶｉｅｗなど）であり、ＰＣのＰＣＩバスなどに後付で取り付けたデジタイザ（高速ＡＤボード）で入力波形を丸ごと取り込み，その上で，様々な用途に向け，その波形を処理する。例えば，デジタルオシロスコープとして利用する場合や、スペクトラムアナライザとして利用することが想定される。

図１に示すように、本発明の第１の実施の形態に係る周波数測定装置１０には、浮上支持した物体を取り付けるための浮上物体支持部１２、及び光波干渉計１４を備えた計測ユニットが設けられている。浮上物体支持部１２は、摩擦抵抗の小さい静圧空気軸受を備えており、静圧空気軸受の可動部１６が、軸受支持部１８によって支持されている。

可動部１６は、水平方向を向くように軸受支持部１８によって支持されており、また、軸受支持部１８には、貫通孔が穿設されており、静圧空気軸受の可動部１６は、この貫通孔を貫通するように配置されている。

この貫通孔の内周面と可動部１６の外周面との間には、可動部１６の外周面全周に亘って圧縮空気層（約５〜１０μｍ）が形成されている。これにより、可動部１６は、極めて小さい摩擦抵抗で可動部１６の長さ方向に直進運動することが可能である。なお、図示を省略したが、圧縮空気は可動部１６の一端から導入され、圧縮空気は、可動部１６の中央部に形成された出口から、軸受支持部１８の内周に形成された溝を介して、軸受支持部１８の中央部に形成された入口に導かれる。そして、この圧縮空気は、軸受支持部１８の内部配管を介して、軸受支持部１８の貫通孔の内周面と可動部１６の外周面との間に供給される。また、給気チューブを介して軸受支持部１８に圧縮空気を直接導入するようにしてもよい。

また、可動部１６の側面後方には、光波干渉計を構成するコーナ・キューブ・プリズム２０が光入出射側に後方を向けて固定されている。

また、浮上物体支持部１２には、軸受支持部１８を載せるためのベース２２が設けられ、ベース２２はＬ字型となっており、ベース２２の一部が、可動部１６の端部に対向するように立設されている。ベース２２の立設部の可動部１６側の側面上方には、可動部１６が当接するためのゴムブロック等で構成された緩衝材２４が固定されている。

ベース２２の立設部と軸受支持部１８との間に設定された基準位置には、基準位置で鉛直方向に沿った光路を形成するように、レーザダイオード２６Ａとフォトダイオード２６Ｂとで構成された基準位置センサ２６が配置されている。基準位置センサ２６のフォトダイオード２６Ｂは、アナログ信号をデジタル信号に変換するＡＤコンバータ及び逆変換するＤＡコンバータからなるコンバータ２８とＰＣＩバス３０とを介してパーソナルコンピュータ３２に接続されている。パーソナルコンピュータ３２には、可動部１６に作用する慣性力を計測するためのプログラムを記憶したＲＯＭ、計測された時系列データを記憶するハードディスク等で構成された記憶媒体、及び記憶した時系列データを表示するＬＣＤやＣＲＴ等で構成された表示装置が設けられている。

光波干渉計１４は、ゼーマンタイプのヘリウムネオンレーザ等から構成された光源３４を備えており、この光源３４から出射された光は、無偏光ビームスプリッタ３６で第１の測定光と第１の参照光とに分光される。また、第１の測定光は、第２の偏光ビームスプリッタ３８に入射され、第２の測定光と第２の参照光とに分光される。

第２の測定光は、可動部１６に固定された第１のコーナ・キューブ・プリズム２０に入射され、第１のコーナ・キューブ・プリズム２０で反転されて光路が１８０°変更され、再度第２の偏光ビームスプリッタ３８に入射される。ここで、第２の偏光ビームスプリッタ３８と第１のコーナ・キューブ・プリズム２０との間の光路の方向は、可動部１６の移動方向と平行になるように設定されている。また、第２の測定光は、第１のコーナ・キューブ・プリズム２０で反転する際に、可動部１６の変位速度に応じたドップラーシフト(ドップラ効果による周波数変化または位相変化)を生じ、反射光の状態が変化する。

一方、第２の参照光は、第２の偏光ビームスプリッタ３８から第２のコーナ・キューブ・プリズム４０方向に反射され、第２のコーナ・キューブ・プリズム４０で反転されて再度第２の偏光ビームスプリッタ３８に入射され、第２の偏光ビームスプリッタ３８で第１のコーナ・キューブ・プリズム２０からの反射光（信号光）と干渉し、信号光と第２の参照光との周波数の差分に相当するビート周波数を生ずる干渉光が生じる。すなわち、光波干渉計１４によって、信号光に生じた状態変化（信号光の周波数のドップラーシフト、または信号光の位相の変化）が、干渉光により検出される。

干渉光は、第２の偏光ビームスプリッタ３８から反射ミラー４２と第１の偏光板（例えば、グラン・トンプソン・プリズム)４４とを介して、フォトダイオードで構成された第１の光検出器４６に入射される。そして、第１の光検出器４６によって、干渉光のビート周波数ｆ_beatに応じた電気信号に変換されて、この電気信号が第１の周波数カウンタ（ｍｏｄｅｌ：Ｒ５３６３；ｍａｎｕｆａｃｔｕｒｅｄｂｙＡｄｖａｎｔｅｓｔＣｏｒｐ．，Ｊａｐａｎ）４８に入力されると共に、デジタイザ（ｍｏｄｅｌ：５１０２；ｍａｎｕｆａｃｔｕｒｅｄｂｙＮａｔｉｏｎａｌＩｎｓｔｒｕｍｅｎｔｓＣｏｒｐ．，ＵＳＡ）５０に入力される。そして、第１の周波数カウンタ４８によりビート周波数ｆ_beatの測定値を表すデジタルデータが生成され、このデジタルデータが、ＧＰＩＢカード５２とＰＣＩバス３０とを介して、パーソナルコンピュータ３２に入力される。また、デジタイザ５０では、第１の光検出器４６から出力される干渉光のビート周波数ｆ_beatに応じた電気信号を、サンプリングレート２０ＭＳ／ｓで記録し、ビート周波数ｆ_beatに応じた電気信号が示す波形情報を記憶する。ここで、ビート周波数ｆ_beatに応じた電気信号がサンプリングデータとして５０００個記録されるようになっており、また、計測時間は０．２５秒となっている。

また、デジタイザ５０に記録された波形情報は、ＰＣＩバス６０を介して、パーソナルコンピュータ６２に入力される。

また、第１の参照光は、第２の偏光板５４を介して、フォトダイオードで構成された第２の光検出器５６に入射され、第２の光検出器５６で、可動部１６の静止状態に対応する基準周波数f_restの電気信号に変換される。そして、この電気信号が入力される第２の周波数カウンタ５８によって基準周波数ｆ_restの値を表すデジタルデータが生成され、このデジタルデータが、ＧＰＩＢカード５２とＰＣＩバス３０とを介して、パーソナルコンピュータ３２に入力される。また、第２の光検出器５６から出力される基準周波数ｆ_restの電気信号が、デジタイザ５０でサンプル記録され、また、ＰＣＩバス６０を介して、基準周波数ｆ_restの電気信号が示す波形情報が、パーソナルコンピュータ６２に入力される。

本実施の形態に係る周波数測定装置１０の光源３４では、周波数が３ＭＨｚ程度異なり、かつ、一方が縦偏光の偏光成分を有し、他方が横偏光の偏光成分を有する２つの周波数レーザを用いている。無偏光ビームスプリッタ３６は、両方の周波数レーザの偏光成分を均等に２分する。一方、第２の偏光ビームスプリッタ３８は縦偏光成分を反射し、横偏光成分を透過する。

無偏光ビームスプリッタ３６によって２分された光の一方を、第２の偏光板５４により４５度の角度で切り出すことにより、光源３４が発振する２つの周波数レーザの差周波数が、ビートとして観測される。

一方、第２の偏光ビームスプリッタ３８によって２分された一方の偏光（この場合、横偏光）の成分だけが、可動部１６において反射されて、可動部１６の速度に比例したドップラーシフトを受ける。第２の偏光ビームスプリッタ３８で２分された２つの偏光成分は、再び、第２の偏光ビームスプリッタ３８で合波され、この合波された光を第１の偏光板４４によって４５度の角度で切出すことにより、２つの周波数レーザの差周波数がビートとして観測される。

ここで、横偏光成分は、可動部１６の速度に比例したドップラーシフトを受けているため、結果として、第１の光検出器４６で観測されるビート周波数ｆ_beatは、第２の光検出器５６で観測される周波数ｆ_restと比べると、ドップラーシフト分だけ異なっている。

なお、２台の周波数カウンタ４８、５８の各々は、第１の光検出器４６及び第２の光検出器５６が出力する電気信号を、サンプリングレート（Ｔ＝４００／ｆ）で、サンプリングデータとして、それぞれ２０００個記録する。本実施の形態の設定では、２台の周波数カウンタ４８、５８は、休止時間なしで、入力信号４００周期ごとに平均周波数を計算し、記録する。

また、パーソナルコンピュータ６２には、可動部１６に作用する慣性力を計測する後述する計測ルーチンを実行するためのプログラムを記憶したＲＯＭ、計測された時系列データを記憶するハードディスク等で構成された記憶媒体、及び記憶した時系列データを表示するＬＣＤやＣＲＴ等で構成された表示装置が設けられている。

次に、本実施の形態の計測ルーチンについて図２を用いて説明する。まず、パーソナルコンピュータ６２を立ち上げ、図２に示す計測ルーチンを起動させると共に、フォトダイオード２６Ｂ、第１の光検出器４６、及び第２の検出部５６にレーザ光が入射されるように、レーザダイオード２６Ａ及び光源３４からレーザ光を発振させる。また、パーソナルコンピュータ３２も同様に立ち上げる。このような状況で、初期段階として可動部１６を長さ方向に手で押して、可動部１６を慣性で移動させる。これにより、第１のコーナ・キューブ・プリズム２０及び可動部１６を一体とした物体（以下、錘体という）が慣性で移動する。なお、アクチュエータを設けて、このアクチュエータで錘体に初速（または、初期運動量、または初期運動エネルギー）を与えても良い。

ステップ１００において、基準位置センサ２６から基準位置信号が入力されたか否かを判断し、レーザダイオード２６Ａからフォトダイオード２６Ｂにいたるレーザ光を、可動部１６が遮り、基準位置センサ２６によって基準位置信号が生成されて、基準位置信号がコンピュータ６２に入力されると、ステップ１０２において、トリガ信号をデジタイザ５０に入力し、デジタイザ５０を起動する。また、トリガ信号がデジタイザ５０に入力されるときに、第１の周波数カウンタ４８及び第２の周波数カウンタ５８にもトリガ信号を入力し、これらを同時に起動する。これにより、デジタイザ５０、第１の周波数カウンタ４８、及び第２の周波数カウンタ５８が、第１の光検出器４６及び第２の検出部５６からの電気信号を同期して検出するようになる。

ここで、サンプリング数（例えば、５０００個）を示すデータを、ＰＣＩバス６０を介してパーソナルコンピュータ６２からデジタイザ５０に予め送信しておく。デジタイザ５０は、第１の検出器４６及び第２の光検出器５６の各々に対して、予め定められたサンプリング数の電気信号を記録し、図３に示すような、サンプリング数の電気信号に基づく波形情報が、ビート周波数ｆ_beat及び基準周波数ｆ_restの各々について、デジタイザ５０に記録される。

次のステップ１０４では、デジタイザ５０から出力されるビート周波数ｆ_beatの波形情報及び基準周波数ｆ_restの波形情報を全波形分取込む。

ここで、図３に示すように、波形情報におけるサンプリング測定間隔Ｐ₀、Ｐ₁、・・・Ｐ_nの各々は、波の数で定められている。

そして、ステップ１０６において、ビート周波数ｆ_beatの波形情報及び基準周波数ｆ_restの波形情報の各々に対して、各サンプリング測定間隔Ｐ_j（サンプリング測定間隔は、Ｎ_j個の波のゼロクロス位置の時刻ｔによって定められる時間範囲（ｔ_jNj〜ｔ_(j+1)Nj）とする）について、第１の時刻としてのサンプリング開始時刻Ｔ_j及び第２の時刻としてのサンプリング終了時刻Ｔ_j+1を算出する。ここで、サンプリング開始時刻Ｔ_jは、サンプリング開始位置のゼロクロス位置（例えば、時刻ｔ２の位置）及び前後同数の２ｎ個のゼロクロス位置を含む（２ｎ＋１）個のゼロクロス位置（ゼロクロス平均個数の半幅ｎ＝２とすると、ｔ₀、ｔ₁、ｔ₂、ｔ₃、ｔ₄、ｔ₅の位置）の平均時刻を、以下の（１）式によって算出される。

また、サンプリング終了時刻Ｔ_j+1も同様に、サンプリング終了位置のゼロクロス位置（例えば、時刻ｔ₁₂の位置）及び前後同数の２ｎ個のゼロクロス位置を含む（２ｎ＋１）個のゼロクロス位置（ｔ₁₀、ｔ₁₁、ｔ₁₂、ｔ₁₃、ｔ₁₄の位置）の平均時刻を算出する。

なお、本実施の形態では、ゼロクロス位置を、マイナス側からプラス側へ移行する隣接する２つのサンプリング測定点の線形補間により求めている。また、ノイズの影響を考慮し、ゼロに非常に近いサンプリング測定点をゼロクロス位置と見なし、そのサンプリング測定点の時刻をそのまま用いている。また、サンプリング測定間隔内の波の個数Ｎ_jは、サンプリング開始時刻Ｔ_jに対応するゼロクロス位置とサンプリング終了時刻Ｔ_j+1のゼロクロス位置との間に存在する波の数となっている。

そして、ステップ１０８で、ビート周波数ｆ_beatの波形情報及び基準周波数ｆ_restの波形情報の各々に対して、各サンプリング測定間隔Ｐ_jを、ステップ１０６で算出されたサンプリング開始時刻Ｔ_j及びサンプリング終了時刻Ｔ_j+1に基づいて、以下の（２）式によって算出する。
Ｐ_j＝Ｔ_j+1−Ｔ_j (２)
なお、本実施の形態では、Ｎ_jがサンプリング測定間隔を定める周期を示しており、例えば、Ｎ_j＝１０（ｊ＝０、１、・・・）となっている。

そして、ステップ１１０で、ビート周波数ｆ_beatの波形情報及び基準周波数ｆ_restの波形情報の各々に対して、サンプリング測定間隔Ｐ_j毎に、各サンプリング測定間隔Ｐ_j及び各サンプリング測定間隔Ｐ_j内の波の個数Ｎ_jに基づいて、以下の（３）式によって、周波数ｆ_jを算出する。
ｆ_j ＝Ｎ_j ／Ｐ_j (３)
そして、ステップ１１２において、サンプリング測定間隔毎に、錘体（第１のコーナ・キューブ・プリズム２０及び可動部１６を一体とした物体）の速度ｖを以下の（４）式及び（５）式に従って演算し、可動部１６の長さ方向の錘体の変位速度ｖを検出する。
ｖ＝λ_air・ｆ_Doppler／２（４）
ｆ_Doppler＝−（ｆ_beat−ｆ_rest）（５）
ここで、λ_airは、空気の屈折率、ｆ_Dopplerは、可動部１６で反射したビーム光のドップラーシフト周波数である。

次のステップ１１４では、この変位速度ｖを、サンプリング測定間隔で積分することにより、可動部１６の長さ方向の錘体の変位（すなわち、基準位置を基準とした位置)ｘをサンプリング測定間隔毎に検出する。また、ステップ１１６において、可動部１６の長さ方向の変位速度ｖを、サンプリング測定間隔で微分することにより、錘体の加速度αをサンプリング測定間隔毎に検出する。

ここで、錘体が軸受支持部１８で支持されていて、錘体と軸受支持部１８との間に作用する摩擦力を無視できる場合には、錘体の質量をＭ、錘体の加速度をαとすると、錘体の慣性力Ｆは、Ｆ＝Ｍ・αで表すことができる。

したがって、ステップ１１８では、検出された加速度αと予め記憶されている錘体の質量Ｍとを乗算することにより、錘体の慣性力Ｆをサンプリング測定間隔毎に検出する。

次のステップ１２０において、上記のステップ１１２〜１１８によって算出された変位速度ｖ、加速度α、位置ｘ、及び慣性力Ｆの各データを時系列的にハードディスク等の記憶媒体に記憶保持する。

そして、ステップ１２２において、キーボードの操作等により時系列データの表示指示があったか否かを判断し、表示指示があった場合にはステップ１２４において記憶媒体から時系列データを読み出して表示装置に表示する。また、時系列データをプリンタ等から適宜出力するようにしてもよい。

また、上記では、デジタイザ５０に記録されている波形情報に基づいて、パーソナルコンピュータ６２で周波数等が計測される様子を説明した。パーソナルコンピュータ３２においては、上記の処理と同時に、計測データの比較のために、周波数カウンタ４８、５８から取り込まれる周波数ｆ_beat、ｆ_restに基づいて、各種データが計測される。

次に、上記の周波数測定装置１０によって、周波数ｆ_beat、ｆ_rest、変位速度ｖ、加速度α、位置ｘ、及び慣性力Ｆを計測した結果を図４Ａ〜図４Ｅに示す。この計測では、ゼロクロス平均個数の半幅ｎ＝１００、サンプリング測定間隔内の波の数Ｎ_j＝４００（ｊ＝０、１、・・・）とした。図４Ａ〜図４Ｅにビート周波数ｆ_beat、基準周波数ｆ_rest、これらの周波数から算出した錘体の変位速度ｖ、変位速度ｖを積分して算出した錘体の変位ｘ、変位速度ｖを微分して算出した錘体の加速度α、錘体の質量Ｍに加速度αを乗算して算出した錘体の慣性力Ｆの測定結果を示す。

また、図５に、測定された慣性力Ｆと時間との関係を示す。Ｆ_digitizerは、ｆ_beatの瞬時値とf_restの瞬時値とを用いて計算したものである。Ｆ_digitizer、_mは、f_beatの瞬時値とf_restの平均値（全測定時間における平均値）とを用いて計算したものである。両者の差（Ｆ_digitizer−Ｆ_digitizer、_m）が、ほぼゼロであることから、測定期間を通してf_restの変動が無視できる程度であることが分かる。このことより、図５に示したパルス状の波形の最後部付近で観測される波形の乱れが、f_restの変動に起因するものでなく、f_beatの変動に起因するものであることが分かる。

次に、図６に、デジタイザ５０を用いた周波数測定から算出した慣性力Ｆ_digitizerと、第１の周波数カウンタ４８及び第２の周波数カウンタ５８を用いた周波数測定から算出した慣性力Ｆ_counterとの比較を示す。測定結果を比較すると、上述したデジタイザ５０を用いた周波数測定が高精度かつ低ノイズに行われていることが分かる。なお、デジタイザ５０における周波数のサンプリングレートは、Ｎ_j＝４００であるので、周波数カウンタのサンプリングレートと全く同一である。すなわち、図６中で、Ｆ_counterの測定点である白丸がある時刻に、必ず、Ｆ_digitizerの測定点がある。

また、図７に、上記の処理により測定した慣性力Ｆ_digitizerと、位置（試料の押し込み深さ）ｘとの関係を示す。高精度かつ高分解能な測定が実現していることが分かる。

次に、図８Ａ〜図８Ｃに、図４Ａ〜図７と同一の測定データを用いた場合のＮ及びｎの効果を示す。図８Ａ〜図８Ｃは、図１の構成の周波数測定装置１０で、デジタイザ５０を用いて測定した周波数から算出した慣性力Ｆ_digitizerと時刻ｔとの関係を示し、図８Ａは、通常の周波数測定アルゴリズムで、Ｎ＝４００、ｎ＝０とした場合であり、図８Ｂは、通常の周波数測定アルゴリズムで、Ｎ＝４０００、ｎ＝０とした場合である。図８Ａと図８Ｂとを比較すると、図８Ｂでは、サンプリング測定間隔内の波の数Ｎを大きく取ったことにより、ノイズは大幅に低減されているが、その分、測定間隔が粗くなっている。また、図８Ｃは、本実施の形態のように、ゼロクロス平均個数の半幅ｎを導入し、Ｎ＝４００、ｎ＝１００とした場合であり、図８Ａと同様に測定間隔が細かく、かつ、図８Ｂと同程度のノイズ低減を達成している。また、この実験例では、図８Ｃにおいて、パルスピーク付近の力の振動を捉えることに成功したことを示している。

ここで、本発明の原理をより一般的に表現すると、以下のように表現される。すなわち、被測定信号（入力信号）について、所定の時刻帯である時刻・位相決定時刻帯に含まれる波形情報から、入力信号の位相θ₁（ｒａｄ／ｓ）と時刻ｔ₁（ｓ）とを求める。次に、別の時刻帯に含まれる情報から、入力信号の位相θ₂（ｒａｄ／ｓ)と時刻ｔ₂（ｓ）とを求める。その結果から、２つの時刻ｔ₁、ｔ₂の間に存在する波の数は、（θ₂−θ₁)／２π と計算される。したがって、周期Ｔ（ｓ）及び周波数ｆ（Ｈｚ）は、以下のように計算できる。
Ｔ＝（２π）（ｔ₂−ｔ₁）／（θ₂−θ₁）
ｆ＝１／Ｔ
また、時刻・位相決定時刻帯を入力信号のゼロクロス位置の個数で規定し、例えば、あるゼロクロス位置の前後ｎ個、合計（２ｎ＋１）個のゼロクロス位置が存在する範囲を時刻・位相決定時刻帯とする。そして、位相θ＝ｍπ（ｍは整数）のゼロクロス位置の時刻ｔを、（２ｎ＋１）個のゼロクロス位置の時刻の平均値を用いて求めている。

以上説明したように、本発明の第１の実施の形態に係る周波数測定装置によれば、サンプリング開始位置に対応するゼロクロス位置を含む所定範囲内の複数のゼロクロス位置の時刻の平均値と、サンプリング終了位置に対応するゼロクロス位置を含む所定範囲内の複数のゼロクロス位置の時刻の平均値とを用いて、波形情報の周波数を測定することにより、多くのゼロクロス位置の時刻を用いているため、サンプリング間隔を拡大させることなく、高精度かつ高分解能で周波数を測定することができる。

また、サンプリング開始位置及びサンプリング終了位置の各々のゼロクロス位置の前後同数のゼロクロス位置の時刻情報を有効活用して、サンプリング間隔を拡大させることなく、高精度かつ高分解能で周波数を測定することができる。

また、従来の手法（すなわち、サンプリング開始位置及び終了位置のゼロクロス位置の時刻のみを用いて周波数を測定する方法）では、デジタイザの性能が一定という条件下において、周波数の測定分解能をあげるには、Ｎを大きくするしかなく、Ｎを大きくすることは、サンプリング間隔の拡大（悪化）を意味していた。しかし、本実施の形態によれば、ゼロクロス平均個数の半幅ｎの導入により、サンプリング間隔の拡大を招くことなく、周波数の測定分解能を向上させることが可能となる。

また、本発明の実施の形態において、サンプリング開始位置又はサンプリング終了位置のゼロクロス位置及び前後同数のゼロクロス位置を含む「平均化を実施する区間」の導入により、「測定区間」の定義が、通常の周波数カウンタのアルゴリズムによるものとは、異なったものになる。このことは、図８Ａ〜図８Ｃ、および、その説明において、端的に示されている。すなわち、本発明のアルゴリズムを用いることにより、ごく僅かな「測定区間」の拡大（平均化区間程度の拡大）により、Ａ／Ｄ変換され取得された波形情報に対して、時間軸上、及び電圧軸上における大きな平均化効果を得ることができる。このことは、Ａ／Ｄ変換における時間軸上のノイズ、及び電圧軸上のノイズの大幅な低減が達成されると解釈することもできる。

本発明の実施の形態に係る周波数測定装置によれば、Ａ／Ｄ変換における時間軸上のノイズ、あるいは、電圧軸上のノイズが問題になるような極限域における、周波数測定（又は周期測定）における時間分解能（又はサンプリング間隔）、及び周波数分解能（又は周期分解能）の両方を向上させる。一方、通常の周波数測定や周期測定器においては、通常、Ａ／Ｄ変換における時間軸上のノイズ、あるいは、電圧軸上のノイズが問題になっており、周波数測定（又は周期測定）における時間分解能（又はサンプリング間隔）、及び周波数分解能（又は周期分解能）がトレードオフの関係にある。

なお、上記の実施の形態では、デジタイザによって記録された全波形測定値に基づいて波形情報の周波数を測定する方法を示したが、上記の（１）式〜（３）式によるアルゴリズムを、通常の周波数カウンタに内蔵するアルゴリズムとしても適用しても良い。この場合には、ｎ＝２、Ｎ＝１０とすると、ハードウェア（電気回路）として、以下のように実現することができる。まず、全部で（２ｎ＋１）×２＝１０個の時刻情報を記録する１０個のレジスタを用意し、５個のレジスタの記憶値を平均する平均化回路を用意する。そして、平均化回路の出力が、通常のカウンタで取り扱う測定時間の両端の時刻に相当するように構成すればよい。

また、サンプリング開始位置における（２ｎ＋１）個のゼロクロス位置とサンプリング終了位置における（２ｎ＋１）個のゼロクロス位置とが重複していない場合を例に説明したが、サンプリング開始位置に対応するゼロクロス位置とサンプリング終了位置に対応するゼロクロス位置とが異なっていれば、各々の（２ｎ＋１）個のゼロクロス位置が一部重複していてもよい。

また、周波数を算出する場合を例に説明したが、波形情報の周期を算出するようにしてもよい。その場合には、サンプリング測定間隔を、サンプリング測定間隔内に存在する波の数で除算して、波形情報の周期を算出すればよい。これにより、サンプリング間隔を拡大させることなく、高精度かつ高分解能で波形情報の周期を測定することができる。

また、タイムインターバル測定を行うようにしてもよい。この場合には、タイムインターバルアナライザの内部動作アルゴリズムを、測定時間の開始位置に対応するゼロクロス位置を含む所定範囲内の複数のゼロクロス位置の時刻の平均値と、測定時間の終了位置に対応するゼロクロス位置を含む所定範囲内の複数のゼロクロス位置の時刻の平均値とを用いて、タイムインターバル測定を行うアルゴリズムにすればよい。

次に、第２の実施の形態について説明する。第２の実施の形態では、サインフィッティングを適用して、サンプリング測定間隔内の波形の時刻及び位相を取得している点が第１の実施の形態と異なっている。なお、周波数測定装置の構成は、第１の実施の形態の構成と同様であるため、同一符号を付して説明を省略する。

周波数測定装置１０のパーソナルコンピュータ６２は、デジタイザ５０から波形情報を取り込むと、時刻及び位相決定時刻帯であるサンプリング測定間隔内の波形を、最小二乗法などにより適当な近似曲線（sin(ωt+θ₀)など）でフィッティングする。これにより、サンプリング測定間隔内における時刻及び位相が既知となるため、第１の実施の形態と同様に、所定のゼロクロス位置の時刻を平均化して、サンプリング開始時刻及びサンプリング終了時刻の各々を算出し、周波数等のデータを計測する。

また、２つのサンプリング測定間隔の位相接続（２πの何倍だけ離れているかの判定）については、測定間隔間のゼロクロスの数を数える方法を用いることが好ましい。

なお、サインフィッティングを用いる方法として、次の２つの場合が考えられる。サインフィッティングを用いる方法の一方は、個々のゼロクロス位置を求める目的で、第１の実施の形態における線形補間の代わりに、サインフィッティングを用いる方法であり、この場合、回帰に用いる範囲は、１／４周期前後となっている。なお、フィッティングする目的関数は正弦波でなく、多項式等でもよい。３次程度の多項式であれば、フィッティングに要する計算時間を節約できるメリットがある。

また、サインフィッティングを用いる方法のもう一方は、数周期程度以上の多くの区間でｙ＝Ａｓｉｎ（ωｔ＋θ₀）でフィッティングする方法であり、Ａ、ω、θ₀を最小二乗法などにより求める。この場合、１回のサインフィッティングにおいて、多周期に渡る平均化操作を行う形になる。この操作を、少しずつ重複を持たせながらシフトした区間に対して行うことで、より高い効果が得られる。

次に、第３の実施の形態について説明する。なお、第１の実施の形態と周波数測定装置の構成は同様であるため、同一符号を付して説明を省略する。

図９に示すように、サンプリング測定間隔は、複数の測定間隔に基づいて算出され、以下の（６）式及び（７）式に基づいて、サンプリング測定間隔Ｐ_jは、２ｎ＋１個（例えば５個）の測定間隔の平均として算出される。
Ｐ_j,k ＝ｔ_Nj+N+k − ｔ_Nj+k (6)

サンプリング開始位置のゼロクロス位置を含む５個のゼロクロス位置の時刻の各々と、サンプリング終了位置のゼロクロス位置を含む５個のゼロクロス位置の時刻の各々とを用いて、ゼロクロス位置を一つずつずらしながら、各々測定間隔を算出し、算出された測定間隔を平均化することにより、サンプリング測定間隔Ｐ_jを求める。

そして、サンプリング測定間隔Ｐ_j及びゼロクロス位置の個数Ｎ_jに基づいて、以下の（8）式によって、周波数ｆ_jが算出される。
ｆ_j = Ｎ_j ／Ｐ_j (8)
このように、本実施の形態では、サンプリング測定間隔を求める段階で平均操作を行う。これに対し、第１の実施の形態では、ゼロクロス位置の時間を求める段階で平均操作を行っているが、第１の実施の形態と第３の実施の形態とにおける周波数を算出するアルゴリズムは数学的に等価である。

なお、Ｎの大きさの設定方法については、大きく分けて、通常の周波数カウンタと同様に２通りある。すなわち、Ｎを設定値として一定とする方法と、目標サンプリング時間を設定値としてＮを可変とする方法とがある。Ｎを一定として設定する方法では、入力周波数に応じて、サンプリング時間は変化する。（サンプリング時間はＮと入力周波数の周期との積である。）一方、目標サンプリング時間を設定値とする方法では、サンプリング時間が目標時間に最も近くなるように、Ｎをその都度変化させる。

以上説明したように、第３の実施の形態に係る周波数測定装置によれば、サンプリング開始位置となるゼロクロス位置を含む所定範囲内の複数のゼロクロス位置の各々と、サンプリング終了位置となるゼロクロス位置を含む所定範囲内の複数のゼロクロス位置の各々とを用いて、ゼロクロス位置をずらしながらサンプリング測定間隔を各々算出し、サンプリング測定間隔の平均値を用いて、波形情報の周波数を測定することにより、多くのゼロクロス位置の時刻を用いているため、サンプリング間隔を拡大させることなく、高精度かつ高分解能で周波数を測定することができる。

なお、（６）式〜（８）式に基づく上記のアルゴリズムを、通常の周波数カウンタに内蔵するアルゴリズムとして適用してもよい。しかし、その場合には、周波数カウンタの必要数が増えてしまうため、周波数カウンタに内蔵するアルゴリズムとして適用する場合には、第１の実施の形態に係る（１）式〜（３）式に基づくアルゴリズムを適用することが好ましい。

衝撃力に対する材料の力学的応答を求めるために周波数を計測する計測装置に適用することによって、材料に関して高精度な測定データを得ることができる。また、本発明を、周波数カウンタや、ＰＣベースの周波数測定システムの、動作アルゴリズムとして適用することにより、計算時間を過大にすることなく、単一周波数測定の精度を大幅に向上することができる。

符号の説明

１０周波数測定装置
３２、６２パーソナルコンピュータ
４６、５６光検出器
４８、５８周波数カウンタ
５０デジタイザ
ｆ_j 周波数
ｆ_beat ビート周波数
ｆ_rest 基準周波数
Ｐ_j サンプリング測定間隔
Ｔ_j サンプリング開始時刻
Ｔ_j サンプリング終了時刻

Claims

入力された波形情報の第１のゼロクロス位置を含む第１の所定範囲内の複数のゼロクロス位置各々の時刻の平均値を、前記第１のゼロクロス位置に対する第１の時刻として算出する第１の算出手段と、
前記波形情報の前記第１のゼロクロス位置とは異なる第２のゼロクロス位置を含む第２の所定範囲内の複数のゼロクロス位置各々の時刻の平均値を、前記第２のゼロクロス位置に対する第２の時刻として算出する第２の算出手段と、
前記第１の時刻、前記第２の時刻、及び前記第１のゼロクロス位置と前記第２のゼロクロス位置との間に存在する波の数に基づいて、前記波形情報の周期及び周波数の少なくとも一方を測定する測定手段と、
を含む周波数測定装置。
前記第１の所定範囲は、第１のゼロクロス位置の前後同数のゼロクロス位置を含む範囲であり、
前記第２の所定範囲は、第２のゼロクロス位置の前後同数のゼロクロス位置を含む範囲である請求項１記載の周波数測定装置。
入力された波形情報の第１のゼロクロス位置を含む第１の所定範囲内の複数のゼロクロス位置各々の時刻の平均値を、前記第１のゼロクロス位置に対する第１の時刻として算出し、
前記波形情報の前記第１のゼロクロス位置とは異なる第２のゼロクロス位置を含む第２の所定範囲内の複数のゼロクロス位置各々の時刻の平均値を、前記第２のゼロクロス位置に対する第２の時刻として算出し、
前記第１の時刻、前記第２の時刻、及び前記第１のゼロクロス位置と前記第２のゼロクロス位置との間に存在する波の数に基づいて、前記波形情報の周期及び周波数の少なくとも一方を測定する
ことを特徴とする周波数測定方法。