(57)【要約】【課題】回転パラメータを精度良く推定する。【解決手段】回転推定システム１は、回転を伴って飛翔する球体１０に対して三軸の方向から電磁波を照射し、球体１０からの反射電磁波を受信して、三軸の反射電磁波の受信強度に関する時系列情報を計測する測定装置５と、三軸の時系列情報を周波数情報にそれぞれ変換する計算装置６の周波数領域変換部６１と、三軸の周波数情報におけるピーク間隔Δｆｘ、Δｆｙ、Δｆｚをそれぞれ特定する閾値処理演算部６２と、特定された三軸のピーク間隔に基づいて、球体１０の三軸の角速度ωｘ，ωｙ，ωｚを推定する角速度分離部６４と、三軸の角速度を成分とする角速度ベクトルωの大きさを球体１０の回転速度として、また、角速度ベクトルωの方向を球体１０の回転軸位置として出力する角速度ベクトル導出部６６と、を備える。【選択図】図２
2014182032.tif 000002
【特許請求の範囲】
【請求項１】
回転を伴って飛翔する回転飛翔体の回転速度及び回転軸位置を推定する回転推定システムであって、
前記回転飛翔体に対して三軸の方向から電磁波を照射し、前記回転飛翔体からの反射電磁波を前記三軸で受信して、前記三軸の反射電磁波の受信強度に関する時系列情報を計測する計測手段と、
前記計測手段により計測された前記三軸の時系列情報を周波数情報にそれぞれ変換する変換手段と、
前記変換手段により変換された前記三軸の周波数情報におけるピーク間隔をそれぞれ特定するピーク間隔特定手段と、
前記ピーク間隔特定手段により特定された前記三軸のピーク間隔に基づいて、前記回転飛翔体の前記三軸の角速度を推定する角速度推定手段と、
前記角速度推定手段により推定された前記三軸の角速度を成分とする角速度ベクトルの大きさを前記回転飛翔体の回転速度として、また、前記角速度ベクトルの方向を前記回転飛翔体の回転軸位置として出力する出力手段と、
を備えることを特徴とする回転飛翔体の回転推定システム。
【請求項２】
前記角速度推定手段は、
前記三軸のピーク間隔のそれぞれが、各軸方向から前記回転飛翔体を視た場合の、視点となる軸を除く他の二軸まわりの回転の角速度を合成した回転速度に対応する情報であることに基づいて、前記回転飛翔体の前記三軸の角速度を推定することを特徴とする、請求項１に記載の回転飛翔体の回転推定システム。
【請求項３】
前記角速度推定手段は、前記三軸のピーク間隔をそれぞれΔｆｘ，Δｆｙ，Δｆｚ、前記三軸の角速度をそれぞれωｘ，ωｙ，ωｚと表すとき、下記の数式を用いて前記三軸の角速度を推定することを特徴とする、請求項２に記載の回転飛翔体の回転推定システム。
【数１】
2014182032.tif 000024
【請求項４】
前記回転飛翔体の軌道変化を前記三軸の方向から取得し、前記回転飛翔体の飛翔方向から予想される軌道変化と比較することによって、前記角速度推定手段により推定された前記角速度の正負方向を判定する判定手段を備える
ことを特徴とする、請求項１〜３のいずれか１項に記載の回転飛翔体の回転推定システム。
【請求項５】
前記判定手段は、
前記回転飛翔体からの反射電磁波の測定開始点をＰ０（ｘ０，ｙ０，ｚ０）、軌道変化前の任意の測定点をＰｈ（ｘｈ，ｙｈ，ｚｈ）、軌道変化後の測定終了点をｐｉ（ｘｉ，ｙｉ，ｚｉ）、前記測定開始点ｐ０から前記測定終了点ｐｉまでの前記三軸の軌道変化量をΔｘ，Δｙ，Δｚ、重力加速度をｇ、初速度をｖ０、前記三軸の風速をｗｘ，ｗｙ，ｗｚと表すとき、下記の判別式によって前記角速度の正負方向を判定する
ことを特徴とする、請求項４に記載の回転飛翔体の回転推定システム。
【数２】
2014182032.tif 000025
【請求項６】
前記ピーク間隔特定手段が、
前記変換手段により変換された前記三軸の周波数情報から所定の閾値以上のピークを抽出して、隣接するピーク間のピーク間隔を算出し、
前記算出したピーク間隔のうち頻度の高いものを抽出し、前記抽出したピーク間隔に基づいて前記ピーク間隔を特定する
ことを特徴とする、請求項１〜５のいずれか１項に記載の回転飛翔体の回転推定システム。
【請求項７】
表示手段を備え、
前記回転飛翔体が球技用ボールであり、
前記出力手段により出力された前記球技用ボールの回転速度及び回転軸位置に基づいて、前記球技用ボールの球種または軌道に関する情報を前記表示手段に表示する
ことを特徴とする、請求項１〜６のいずれか１項に記載の回転飛翔体の回転推定システム。
【請求項８】
回転を伴って飛翔する回転飛翔体の回転速度及び回転軸位置を推定する回転推定方法であって、
前記回転飛翔体に対して三軸の方向から電磁波を照射し、前記回転飛翔体からの反射電磁波を前記三軸で受信して、前記三軸の反射電磁波の受信強度に関する時系列情報を計測する計測ステップと、
前記計測ステップにより計測された前記三軸の時系列情報を周波数情報にそれぞれ変換する変換ステップと、
前記変換ステップにより変換された前記三軸の周波数情報におけるピーク間隔をそれぞれ特定するピーク間隔特定ステップと、
前記ピーク間隔特定ステップにより特定された前記三軸のピーク間隔に基づいて、前記回転飛翔体の前記三軸の角速度を推定する角速度推定ステップと、
前記角速度推定ステップにより推定された前記三軸の角速度を成分とする角速度ベクトルの大きさを前記回転飛翔体の回転速度として、また、前記角速度ベクトルの方向を前記回転飛翔体の回転軸位置として出力する出力ステップと、
を備えることを特徴とする回転飛翔体の回転推定方法。
【請求項９】
回転を伴って飛翔する回転飛翔体の回転速度及び回転軸位置を推定する回転推定プログラムであって、
前記回転飛翔体に対して三軸の方向から電磁波を照射し、前記回転飛翔体からの反射電磁波を前記三軸で受信して、前記三軸の反射電磁波の受信強度に関する時系列情報を計測する計測機能と、
前記計測機能により計測された前記三軸の時系列情報を周波数情報にそれぞれ変換する変換機能と、
前記変換機能により変換された前記三軸の周波数情報におけるピーク間隔をそれぞれ特定するピーク間隔特定機能と、
前記ピーク間隔特定機能により特定された前記三軸のピーク間隔に基づいて、前記回転飛翔体の前記三軸の角速度を推定する角速度推定機能と、
前記角速度推定機能により推定された前記三軸の角速度を成分とする角速度ベクトルの大きさを前記回転飛翔体の回転速度として、また、前記角速度ベクトルの方向を前記回転飛翔体の回転軸位置として出力する出力機能と、
をコンピュータに実現させるための回転飛翔体の回転推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、回転飛翔体の回転推定システム、回転推定方法、及び回転推定プログラムに関する。
【背景技術】
【０００２】
従来、例えば投球時の野球ボールや打ち出し後のゴルフボールの軌道、回転速度、回転方向など、回転を伴って飛翔する回転飛翔体の回転パラメータを計測する手法が知られている。例えば特許文献１には、ドップラーシフトから球体の回転速度及び軌道を求める技術が開示されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特許第４８６５７３５号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、従来の回転飛翔体の回転パラメータを計測する手法では、回転パラメータの計測精度にさらなる改善の余地があった。
【０００５】
本発明は、上記に鑑みてなされたものであって、回転飛翔体の回転パラメータを精度良く推定できる回転飛翔体の回転推定システム、回転推定方法、及び回転推定プログラムを提供することを目的とする。
【課題を解決するための手段】
【０００６】
上記課題を解決するために、本発明に係る回転飛翔体の回転推定システムは、回転を伴って飛翔する回転飛翔体の回転速度及び回転軸位置を推定する回転推定システムであって、前記回転飛翔体に対して三軸の方向から電磁波を照射し、前記回転飛翔体からの反射電磁波を前記三軸で受信して、前記三軸の反射電磁波の受信強度に関する時系列情報を計測する計測手段と、前記計測手段により計測された前記三軸の時系列情報を周波数情報にそれぞれ変換する変換手段と、前記変換手段により変換された前記三軸の周波数情報におけるピーク間隔をそれぞれ特定するピーク間隔特定手段と、前記ピーク間隔特定手段により特定された前記三軸のピーク間隔に基づいて、前記回転飛翔体の前記三軸の角速度を推定する角速度推定手段と、前記角速度推定手段により推定された前記三軸の角速度を成分とする角速度ベクトルの大きさを前記回転飛翔体の回転速度として、また、前記角速度ベクトルの方向を前記回転飛翔体の回転軸位置として出力する出力手段と、を備えることを特徴とする。
【０００７】
また、上記の回転推定システムにおいて、前記角速度推定手段は、前記三軸のピーク間隔のそれぞれが、各軸方向から前記回転飛翔体を視た場合の、視点となる軸を除く他の二軸まわりの回転の角速度を合成した回転速度に対応する情報であるとの特性に基づいて、前記回転飛翔体の前記三軸の角速度を推定することが好ましい。
【０００８】
また、上記の回転推定システムにおいて、前記角速度推定手段は、前記三軸のピーク間隔をそれぞれΔｆｘ、Δｆｙ、Δｆｚ、前記三軸の角速度をそれぞれωｘ，ωｙ，ωｚと表すとき、下記の数式を用いて前記三軸の角速度を推定することが好ましい。
【数１】
2014182032.tif 000003
【０００９】
また、上記の回転推定システムは、前記回転飛翔体の軌道変化を前記三軸の方向から取得し、前記回転飛翔体の飛翔方向から予想される軌道変化と比較することによって、前記角速度推定手段により推定された前記角速度の正負方向を判定する判定手段を備えることが好ましい。
【００１０】
また、上記の回転推定システムにおいて、前記判定手段は、前記回転飛翔体からの反射電磁波の測定開始点をＰ０（ｘ０，ｙ０，ｚ０）、軌道変化前の任意の測定点をＰｈ（ｘｈ，ｙｈ，ｚｈ）、軌道変化後の測定終了点をｐｉ（ｘｉ，ｙｉ，ｚｉ）、前記測定開始点ｐ０から前記測定終了点ｐｉまでの前記三軸の軌道変化量をΔｘ，Δｙ，Δｚ、重力加速度をｇ、初速度をｖ０、前記三軸の風速をｗｘ，ｗｙ，ｗｚと表すとき、下記の判別式によって前記角速度の正負方向を判定することが好ましい。
【数２】
2014182032.tif 000004
【００１１】
また、上記の回転推定システムにおいて、前記ピーク間隔特定手段が、前記変換手段により変換された前記三軸の周波数情報から所定の閾値以上のピークを抽出して、隣接するピーク間のピーク間隔を算出し、前記算出したピーク間隔のうち頻度の高いものを抽出し、前記抽出したピーク間隔に基づいて前記ピーク間隔を特定することが好ましい。
【００１２】
また、上記の回転推定システムは、表示手段を備え、前記回転飛翔体が球技用ボールであり、前記出力手段により出力された前記球技用ボールの回転速度及び回転軸位置に基づいて、前記球技用ボールの球種または軌道に関する情報を前記表示手段に表示することが好ましい。
【００１３】
同様に、上記課題を解決するために、本発明に係る回転飛翔体の回転推定方法は、回転を伴って飛翔する回転飛翔体の回転速度及び回転軸位置を推定する回転推定方法であって、前記回転飛翔体に対して三軸の方向から電磁波を照射し、前記回転飛翔体からの反射電磁波を前記三軸で受信して、前記三軸の反射電磁波の受信強度に関する時系列情報を計測する計測ステップと、前記計測ステップにより計測された前記三軸の時系列情報を周波数情報にそれぞれ変換する変換ステップと、前記変換ステップにより変換された前記三軸の周波数情報におけるピーク間隔をそれぞれ特定するピーク間隔特定ステップと、前記ピーク間隔特定ステップにより特定された前記三軸のピーク間隔に基づいて、前記回転飛翔体の前記三軸の角速度を推定する角速度推定ステップと、前記角速度推定ステップにより推定された前記三軸の角速度を成分とする角速度ベクトルの大きさを前記回転飛翔体の回転速度として、また、前記角速度ベクトルの方向を前記回転飛翔体の回転軸位置として出力する出力ステップと、を備えることを特徴とする。
【００１４】
同様に、上記課題を解決するために、本発明に係る回転飛翔体の回転推定プログラムは、回転を伴って飛翔する回転飛翔体の回転速度及び回転軸位置を推定する回転推定プログラムであって、前記回転飛翔体に対して三軸の方向から電磁波を照射し、前記回転飛翔体からの反射電磁波を前記三軸で受信して、前記三軸の反射電磁波の受信強度に関する時系列情報を計測する計測機能と、前記計測機能により計測された前記三軸の時系列情報を周波数情報にそれぞれ変換する変換機能と、前記変換機能により変換された前記三軸の周波数情報におけるピーク間隔をそれぞれ特定するピーク間隔特定機能と、前記ピーク間隔特定機能により特定された前記三軸のピーク間隔に基づいて、前記回転飛翔体の前記三軸の角速度を推定する角速度推定機能と、前記角速度推定機能により推定された前記三軸の角速度を成分とする角速度ベクトルの大きさを前記回転飛翔体の回転速度として、また、前記角速度ベクトルの方向を前記回転飛翔体の回転軸位置として出力する出力機能と、をコンピュータに実現させることを特徴とする。
【発明の効果】
【００１５】
本発明に係る回転飛翔体の回転推定システム、回転推定方法、及び回転推定プログラムは、回転飛翔体の回転パラメータを精度良く推定できるという効果を奏する。
【図面の簡単な説明】
【００１６】
【図１】図１は、本発明の一実施形態に係る回転飛翔体の回転推定システムの概略構成を示す模式図である。
【図２】図２は、図１中の制御装置の機能ブロック図である。
【図３】図３は、本実施形態の回転推定システムにより実施される球体の角速度ベクトルの導出処理を示すフローチャートである。
【図４】図４は、野球用軟式球の１回転当りのレーダ断面積の角度変化を示す図である。
【図５】図５は、コーナーリフレクターを内蔵した球体の１回転当りのレーダ断面積の角度変化を示す図である。
【図６】図６は、レーダの球体追跡手法を説明するための図である。
【図７】図７は、ｘ軸の受信強度の時間変化の一例を示す図である。
【図８】図８は、ｙ軸の受信強度の時間変化の一例を示す図である。
【図９】図９は、ｚ軸の受信強度の時間変化の一例を示す図である。
【図１０】図１０は、ｘ軸の受信強度の周波数情報の一例を示す図である。
【図１１】図１１は、ｙ軸の受信強度の周波数情報の一例を示す図である。
【図１２】図１２は、ｚ軸の受信強度の周波数情報の一例を示す図である。
【図１３】図１３は、受信強度の周波数情報において、本来発生するスペクトルのピークが得られない場合を示す図である。
【図１４】図１４は、受信強度の周波数情報において、球体の回転以外の成分のスペクトルが得られた場合を示す図である。
【図１５】図１５は、図３のフローチャートのステップＳ１０３における閾値処理及び最頻値処理のサブルーチンを示すフローチャートである。
【図１６】図１６は、閾値処理における閾値設置手法の一例を示す図である。
【図１７】図１７は、ｘｙｚ座標系における球体の軌道変化を示す模式図である。
【図１８】図１８は、図１７中の球体の軌道をｚｘ平面座標に投射した図である。
【図１９】図１９は、図１７中の球体の軌道をｘｙ平面座標に投射した図である。
【図２０】図２０は、角速度ベクトルの概略を示す模式図である。
【図２１】図２１は、野球のカーブボールを３方向から測定した実施例の構成を示す模式図である。
【図２２】図２２は、実施例で取得されたｘ軸の受信強度の時間変化を示す図である。
【図２３】図２３は、実施例で取得されたｙ軸の受信強度の時間変化を示す図である。
【図２４】図２４は、実施例で取得されたｚ軸の受信強度の時間変化を示す図である。
【図２５】図２５は、ｘ軸の受信強度の周波数スペクトルを示す図である。
【図２６】図２６は、ｙ軸の受信強度の周波数スペクトルを示す図である。
【図２７】図２７は、ｚ軸の受信強度の周波数スペクトルを示す図である。
【図２８】図２８は、ピーク間隔を再取得したｙ軸の受信強度の周波数スペクトルを示す図である。
【図２９】図２９は、ピーク間隔を再取得したｚ軸の受信強度の周波数スペクトルを示す図である。
【図３０】図３０は、実施例における球体のｘ軸の軌道変化を示す図である。
【図３１】図３１は、実施例における球体のｙ軸の軌道変化を示す図である。
【図３２】図３２は、実施例における球体のｚ軸の軌道変化を示す図である。
【発明を実施するための形態】
【００１７】
以下に本発明に係る回転飛翔体の回転推定システムの実施形態を図面に基づいて説明する。なお、以下の図面において、同一または相当する部分には同一の参照番号を付し、その説明は繰り返さない。
【００１８】
［実施形態］
まず、図１、２を参照して、本実施形態に係る回転推定システム１の構成について説明する。図１は、本発明の一実施形態に係る回転推定システムの概略構成を示す模式図であり、図２は、図１中の制御装置の機能ブロック図である。
【００１９】
（回転推定システムの構成）
図１に示すように、回転推定システム１は、回転して飛翔する回転飛翔体１０に対して三軸方向から電磁波を照射し、受信される反射電磁波を用いて、この回転飛翔体の回転速度及び回転軸方向を推定するものである。本実施形態では、回転飛翔体１０の一例として、野球用ボールなどの球体１０を例示する。また、本実施形態では、球体１０に電磁波を照射し、反射電磁波を受信する三軸として、球体１０を原点として相互に直交するｘ軸、ｙ軸、ｚ軸を設定する。ここで、ｚ軸は鉛直上方を正方向とし、ｘ軸及びｙ軸の方向は、このｚ軸に直交する水平方向とする。球体１０は、ｘ軸に沿った初速で、ｘ軸の正方向（図１に「球体移動方向」と示す方向）に飛翔を開始するものとする。
【００２０】
図１に示すように、回転推定システム１は、ｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４、測定装置５（計測手段）、計算装置６、記録装置７、出力装置８（表示手段）を備える。
【００２１】
ｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４は、それぞれｘ、ｙ、ｚ軸上に配置され、ｘ、ｙ、ｚ軸方向から球体１０に対して電磁波１１，１３，１５を照射して、球体１０によって反射された反射電磁波１２，１４，１６を受信する。電磁波は、連続波またはパルス波を用いる。
【００２２】
測定装置５は、ｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４によって得られた反射電磁波１２，１４，１６の受信強度の時間変化（時系列情報）、電波到来時間、受信周波数などを測定する。
【００２３】
計算装置６は、測定装置５によって得られた受信信号の周波数領域変換、閾値処理、最頻値処理、角速度分離、角速度ベクトル導出、ドップラーシフト導出、ドップラーシフトあるいは電波到来時間を用いた初速度導出、軌道変化導出、角速度の曖昧さ除去などの処理を行う。
【００２４】
計算装置６は、上記の機能に関して、図２に示すように、周波数領域変換部６１（変換手段）、閾値処理演算部６２（ピーク間隔特定手段）、最頻値処理演算部６３（ピーク間隔特定手段）、角速度分離部６４（角速度推定手段）、角速度方向判定部６５（判定手段）、角速度ベクトル導出部６６（出力手段）の各機能を実現するよう構成されている。
【００２５】
周波数領域変換部６１は、測定装置５により取得されたｘ，ｙ，ｚ軸の時系列情報を周波数情報にそれぞれ変換する。
【００２６】
閾値処理演算部６２は、周波数領域変換部６１により変換されたｘ，ｙ，ｚ軸の周波数情報から所定の閾値以上のピークを抽出して、隣接するピーク間のピーク間隔を算出する。
【００２７】
最頻値処理演算部６３は、閾値処理演算部６２により算出されたピーク間隔のうち頻度の高いものを抽出し、この抽出したピーク間隔に基づいてピーク間隔Δｆｘ，Δｆｙ，Δｆｚを特定する。
【００２８】
角速度分離部６４は、特定されたｘ，ｙ，ｚ軸のピーク間隔Δｆｘ，Δｆｙ，Δｆｚに基づいて、球体１０のｘ，ｙ，ｚ軸の角速度ωｘ，ωｙ，ωｚを推定する。
【００２９】
角速度方向判定部６５は、角速度分離部６４により推定された角速度ωｘ，ωｙ，ωｚの正負方向を判定する。
【００３０】
角速度ベクトル導出部６６は、角速度分離部６４により推定されたｘ，ｙ，ｚ軸の角速度ωｘ，ωｙ，ωｚを成分とする角速度ベクトルωを導出する。なお、計算装置６の各機能ブロックの処理の詳細については図３などを参照して後述する。
【００３１】
記録装置７は、計算装置６によって得られた値を保存する装置であり、より詳細には、角速度ベクトル導出部６６により導出された角速度ベクトルに関する情報や、この角速度ベクトルから導出できる球体１０の回転速度や回転軸方向などの情報を保存する。
【００３２】
出力装置８は、計算装置６によって得られた値の出力を行い、ディスプレイへの角速度ベクトル表示や、インターネット配信などを行う。
【００３３】
測定装置５、計算装置６、記録装置７、出力装置８は、物理的には、ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、ＥＥＰＲＯＭ（Ｅｌｅｃｔｒｉｃａｌｌｙ  Ｅｒａｓａｂｌｅ  ａｎｄ  Ｐｒｏｇｒａｍｍａｂｌｅ  Ｒｅａｄ  Ｏｎｌｙ  Ｍｏｍｏｒｙ）やＨＤＤ（Ｈａｒｄ  Ｄｉｓｋ  Ｄｒｉｖｅ）等の記憶装置、装置内外の各部との通信を行なうインタフェース、スイッチ、キーボート、マウスなどの入力装置、ディスプレイなどの表示装置、等のハードウェアを有するマイクロコンピュータである。本実施形態で説明する測定装置５、計算装置６、記録装置７、出力装置８の各機能の全部または一部は、ＣＰＵ、ＲＡＭ、ＲＯＭ等のハードウェア上に所定のアプリケーションプログラムを読み込ませることにより、ＣＰＵの制御の元でインタフェース、入力装置、表示装置等を動作させると共に、ＲＡＭ、ＲＯＭ、記憶装置におけるデータの読み出し及び書き込みを行うことで実現される。
【００３４】
また、上記の測定装置５、計算装置６、記録装置７、出力装置８のアプリケーションプログラム（特に測定装置５や、計算装置６の周波数領域変換部６１、閾値処理演算部６２、最頻値処理演算部６３、角速度分離部６４、角速度方向判定部６５、角速度ベクトル導出部６６に関する機能）を、コンピュータで読み取り可能な記録媒体に格納しても良く、また、プログラム製品として構成することもできる。ここで、この「記録媒体」とは、メモリーカード、ＵＳＢメモリ、ＳＤカード、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＥＰＲＯＭ、ＥＥＰＲＯＭ、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤ、及び、Ｂｌｕ−ｒａｙ  Ｄｉｓｃ等の任意の「可搬用の物理媒体」を含むものとする。また、アプリケーションプログラムは、回転推定システム１に対して任意のネットワークを介して接続されたアプリケーションプログラムサーバに記憶されていてもよく、必要に応じてその全部または一部をダウンロードすることも可能である。
【００３５】
（回転推定システムの動作）
次に、図３〜２０を参照して、本実施形態の回転推定システム１の動作を説明する。図３は、本実施形態の回転推定システムにより実施される球体の角速度ベクトルの導出処理を示すフローチャートである。図３のフローチャートの処理は、計測対象の球体１０が飛翔している間、例えば所定時間ごとに実施される。
【００３６】
（反射電磁波の計測処理）
ステップＳ１０１では、測定装置５により、ｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４を利用して、ｘ、ｙ、ｚ軸方向から球体１０に対して電磁波１１，１３，１５が照射され、球体１０によって反射された反射電磁波１２，１４，１６が所定時間分受信される。測定装置５は、ｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４により受信された所定時間分の反射電磁波１２，１４，１６を用いて、ｘ、ｙ、ｚ軸の反射電磁波１２，１４，１６に関して受信強度の時間変化（時系列情報）を計測する。
【００３７】
ここで図４，５を参照して、本実施形態において、反射電磁波１２，１４，１６の受信強度の時間変化を計測する理由について説明する。図４は、野球用軟式球の１回転当りのレーダ断面積の角度変化を示す図である。図４の横軸は、球体１０の一例としての野球用軟式球の回転角度（ｄｅｇ）を示し、縦軸はレーダ断面積（ｄＢ・ｍ２）を示す。「レーダ断面積（Ｒａｄａｒ  Ｃｒｏｓｓ  Ｓｅｃｔｉｏｎ  ：ＲＣＳ）」とは、物体が電磁波を散乱する度合いを表す量であり、受信電磁波の受信強度に応じて変動するパラメータである。図４に示すように、野球用軟式球の回転角度に応じて、レーダ断面積は不均一に変化する。したがって、球体１０の回転の位相変化と、反射電磁波１２，１４，１６の受信強度の時間変化との間には強い相関関係があると考えることができる。本実施形態の回転推定システム１は、最終的に球体１０の回転軸方向及び回転速度を推定するものである。これらの推定を精度良く行なうためには、推定のための入力情報として、球体１０の回転位相に関する情報を含むことが好ましい。そこで、本実施形態では、球体１０の回転位相と強い相関関係にある反射電磁波１２，１４，１６の受信強度の時間変化を利用する構成としている。
【００３８】
なお、計測対象のＲＣＳが小さく、他の物体と見分けがつきにくい場合には、例えば特許第５１４２３６６号公報に開示される「遊戯用ボール」を使用することもできる。この遊戯用ボールは、８個のコーナーリフレクターを互いに直交するように構成された球体リフレクタである。図５は、この遊戯用ボールの１回転当りのレーダ断面積の角度変化を示す図である。回転推定システム１の計測対象に遊戯用ボールの構成を適用することで、ＲＣＳの大幅な改善と角度変化の安定性を図ることができる。また、反射電磁波をレーダで検知しやすくできる。
【００３９】
なお、ｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４と球体１０との間の距離が短い場合には、球体１０の移動に応じてｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４も移動させ、球体１０が常に三軸の原点に配置されるよう球体１０を追跡する構成としてもよい。
【００４０】
球体の追跡手法として、例えば電磁波１１，１３，１５の照射から反射電磁波１２，１４，１６の受信までの到来時間の時間変化を利用することができる。図６は、レーダの球体追跡手法を説明するための図である。図６に示すように、ｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４における到来時間の時間変化量をそれぞれΔｔｘ，Δｔｙ，Δｔｚとすると、ｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４と球体１０との間の距離の時間変化量ΔＲｘ，ΔＲｙ，ΔＲｚは下記の（１）式で算出できる。
【数３】
2014182032.tif 000005
ここで、ｃは光速度である。
【００４１】
本実施形態では、球体１０は、ｘ，ｙ，ｚ軸から成るｘｙｚ座標系において、飛翔中に常に原点に位置するよう設定される。上記（１）式で算出されたΔＲｘ，ΔＲｙ，ΔＲｚは、ｘｙｚ座標系における球体１０の原点から各軸方向への移動量である。したがって、この移動量に応じて、ｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４の位置を移動する制御を行うことで、球体１０をｘｙｚ座標系の原点に維持しつつ追尾することができる。
【００４２】
例えば、ｘ軸レーダ２の場合、ｙ軸方向の移動量ΔＲｙが正の場合にはｙ軸の正方向（図６では「＋Ｒｙ」と示す方向）に、また、ΔＲｙが負の場合にはｙ軸の負方向（図６では「−Ｒｙ」と示す方向）に、ｘ軸レーダ２の位置を｜ΔＲｙ｜分だけ移動する。また、ｚ軸方向の移動量ΔＲｚが正の場合にはｚ軸の正方向（図６では「＋Ｒｚ」と示す方向）に、また、ΔＲｚが負の場合にはｚ軸の負方向（図６では「−Ｒｚ」と示す方向）に、ｘ軸レーダ２の位置を｜ΔＲｚ｜分だけ移動する。同様に、ｙ軸レーダ３は、＋Ｒｘ方向または−Ｒｘ方向に｜ΔＲｘ｜分だけ移動し、また、＋Ｒｚ方向または−Ｒｚ方向に｜ΔＲｚ｜分だけ移動する。ｚ軸レーダ４は、＋Ｒｘ方向または−Ｒｘ方向に｜ΔＲｘ｜分だけ移動し、また、＋Ｒｙ方向または−Ｒｙ方向に｜ΔＲｙ｜分だけ移動する。
【００４３】
また、ｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４により受信された反射電磁波１２，１４，１６に含まれる球体１０以外から発生した雑音は、例えば移動目標指示装置（Ｍｏｖｉｎｇ  Ｔａｒｇｅｔ  Ｉｎｄｉｃａｔｏｒ：ＭＴＩ）などを用いることによって取り除くことができる。
【００４４】
ステップＳ１０１において計測されるｘ、ｙ、ｚ軸の受信強度の時間変化の一例を図７〜９に示す。図７は、ｘ軸の受信強度の時間変化の一例を示す図であり、図８は、ｙ軸の受信強度の時間変化の一例を示す図であり、図９は、ｚ軸の受信強度の時間変化の一例を示す図である。図７〜９の縦軸は相対受信強度（ｄＢ）を示し、横軸は計測時間（ｓ（秒））を示す。図７〜９に示すグラフは、ｚ軸を中心に秒間１回転させた野球用軟式球のｘ、ｙ、ｚ軸の相対受信強度の時間変化である。ここで「相対受信強度」とは、受信強度の最小値を基準とする受信強度の相対値である。ステップＳ１０１の処理が完了するとステップＳ１０２に移行する。
【００４５】
（周波数変換処理）
ステップＳ１０２では、計算装置６の周波数領域変換部６１により、ｘ、ｙ、ｚ軸の受信強度の時間変化（時系列情報）が周波数領域の情報（周波数情報）に変換される。周波数領域変換部６１は、フーリエ変換、ウェーブレット変換、ウィグナー分布等の周知の周波数変換手法を用いることができる。
【００４６】
また、周波数領域に変換する測定値は測定開始点と終了点の値が同一であり、かつ測定範囲内の周期が整数でなければならないという制約がある。しかし、実測においてはこの条件を常に満たすことは不可能であり、これを無視して周波数領域に変換すると周波数分解能が低下する、ピーク周辺に余計な値が発生するといった問題が起こる。よって、この問題を解消するために窓関数を用いる。窓関数とはある有限区間以外で０となる関数を指し、これを信号にかけることによって測定開始点と終了点の値が同一となり、かつ周期が整数となるため、周波数分解能が向上する。この窓関数にはハミング窓、ハニング窓、矩形窓などといった種類があるが、周波数領域変換部６１は、これらの周知の窓関数を用いることができる。
【００４７】
ステップＳ１０２において変換されたｘ，ｙ，ｚ軸の受信強度の周波数情報の一例を図１０〜１２に示す。図１０は、ｘ軸の受信強度の周波数情報の一例を示す図であり、図１１は、ｙ軸の受信強度の周波数情報の一例を示す図であり、図１２は、ｚ軸の受信強度の周波数情報の一例を示す図である。図１０〜１２の縦軸はパワースペクトル（ｄＢ）を示し、横軸は周波数（Ｈｚ／ｒｐｓ）を示す。図１０〜１２に示す周波数情報は、それぞれ図７〜９に示すｘ，ｙ，ｚ軸の受信強度の時間変化を周波数変換したものである。この例では、周波数変換手法としてフーリエ変換を使用し、窓関数にはハニング窓を使用した。図１０〜１２に示すように、周波数領域では、ｘ，ｙ，ｚ軸の受信強度は周期的なピークを得ることができる。なお、この時の隣接するピーク同士の周波数の間隔は、各軸における回転速度と一致するものである。ステップＳ１０２の処理が完了するとステップＳ１０３に移行する。
【００４８】
（閾値処理及び最頻値処理）
ステップＳ１０３では、計算装置６の閾値処理演算部６２により閾値処理が行われ、最頻値処理演算部６３により最頻値処理が行われる。図１０〜１２に示したように、ｘ，ｙ，ｚ軸の受信強度の周波数スペクトルは、ピーク間隔が各軸における回転速度と一致するような、周期的なピークを有するものである。しかしながら、例えば、図１３中に点線で示す領域のように、本来発生する周波数のピークが得られない場合や、図１４中に点線で示す領域のように、球体１０の回転とは関係のないピークが発生する場合が生じ得る。この場合、周波数スペクトルから正しいピーク間隔が得られず、回転速度の推定に悪影響が出る虞がある。本ステップにおける閾値処理及び最頻値処理は、図１３，１４に点線で示すような、正しいピーク間隔を取得するのに不要となる要素を取り除くために実施する。
【００４９】
ステップＳ１０３では、図１５に示すサブルーチンが実行される。図１５は、ステップＳ１０３における閾値処理及び最頻値処理のサブルーチンを示すフローチャートである。なお、図１５と以下の説明では、便宜上ｘ、ｙ、ｚ軸の符号が省略されているが、図１５のフローチャートに示す一連の処理は、ステップＳ１０２で取得されたｘ、ｙ、ｚ軸の受信強度の周波数情報のそれぞれについて個別に行なわれる。
【００５０】
図１５のフローチャートに沿ってステップＳ１０３の処理の詳細を説明すると、まず、閾値処理演算部６２により任意の閾値が設定され、受信強度の周波数情報からこの閾値以上のピークを抽出する閾値処理が行われる（ステップＳ２０１）。
【００５１】
ここで、図１６を参照して閾値設定手法の一例を説明する。図１６は、閾値処理における閾値設置手法の一例を示す図である。図１６の縦軸はパワースペクトル（ｄＢ）を示し、横軸は周波数（Ｈｚ／ｒｐｓ）を示す。一般に、レーダの受信強度の周波数スペクトルには、周波数０の周辺の領域に直流成分が含まれる。レーダの受信電力の増大に応じて、この直流成分は増大する。このような直流成分の影響により、例えば図１６に示すように、周波数０におけるパワースペクトルのピーク値が、他の成分に対して極端に大きくなる場合がある。閾値処理においてピークを精度良く抽出するためには、このような直流成分の影響を除外して閾値を設定することが望ましい。そこで、例えば図１６に示すように、直流成分が含まれる０〜ｆａまでの範囲を除外して、周波数ｆａ以上の領域でのピークの最大値をとる周波数ｆｍａｘを特定する。そして、周波数ｆａからｆｍａｘの区間におけるパワースペクトルの平均値を算出して、この平均値を閾値として設定する。直流成分を含む領域の境界の周波数ｆａは、例えばパワースペクトルの微分値が負から正に切り替わる周波数、すなわちパワースペクトルが周波数０における最大値から減少し増加に転じる周波数である。なお、図１６を参照して説明した閾値設定手法は一例であり、これ以外の手法を適宜用いることができる。
【００５２】
図１５に戻り、閾値処理演算部６２により、各軸の受信強度の周波数スペクトルからピークｆ１，ｆ２，ｆ３，…ｆｉが取得される（ステップＳ２０２）。閾値処理演算部６２は、例えば図１０〜１２に示すように、ｘ，ｙ，ｚ軸の周波数スペクトルごとに所定の閾値を設定し、この閾値以上の周波数スペクトルからピークとなる周波数をｆ１，ｆ２，ｆ３，…ｆｉがとして取得する。取得された各ピークの符号は、例えば、周波数が小さい方から大きい方へ順番に付与される。そして、隣り合うピーク同士のピーク間隔Δｆ１，Δｆ２，Δｆ３、・・・Δｆｉが取得される（ステップＳ２０３）。例えば、１番目のピーク間隔Δｆ１は、Δｆ１＝ｆ２−ｆ１で算出できる。
【００５３】
続いて最頻値処理演算部６３により最頻値処理が行われる。まずステップＳ２０３で取得されたピーク間隔Δｆ１，Δｆ２，Δｆ３、・・・Δｆｉのうち最小値Δｆｍｉｎ１を取得し、Δｆｍｉｎ１≦Δｆｍｉｎ１＋２Δｆｒの範囲にあるピーク間隔を取得し、それらをグループ１とする。ここで、Δｆｒは、フーリエ変換後の周波数分解能である。
【００５４】
次に、ステップＳ２０３で取得されたピーク間隔Δｆ１，Δｆ２，Δｆ３、・・・Δｆｉからグループ１を除いた中から最小値Δｆｍｉｎ２を取得し、Δｆｍｉｎ２≦Δｆｍｉｎ２＋２Δｆｒの範囲にあるピーク間隔を取得し、それらをグループ２とする。
【００５５】
この演算を、ステップＳ２０３で取得されたピーク間隔Δｆ１，Δｆ２，Δｆ３、・・・Δｆｉのすべてがいずれかのグループに属するまで行なった後、最も要素数の多いグループが、最頻値グループΔｆｍ１として取得される（ステップＳ２０４）。
【００５６】
次に、最頻値グループΔｆｍ１の要素数が４個未満か否かが判定される（ステップＳ２０５）。最頻値グループΔｆｍ１の要素数が４個未満の場合には（ステップＳ２０５のＹｅｓ）、ピーク数が不十分であり周波数情報に周期的な回転の傾向が無いため、その軸における回転速度は発生していないものとみなし、該当する軸の周波数スペクトルのピーク間隔Δｆが０と設定されて（ステップＳ２０６）、メインフローに戻る。
【００５７】
一方、最頻値グループΔｆｍ１の要素数が４個以上の場合には（ステップＳ２０５のＮｏ）、続いて、この最頻値グループΔｆｍ１の要素数が、ステップＳ２０３で取得されたピーク間隔Δｆ１，Δｆ２，Δｆ３、・・・Δｆｉの全体数の半分以上か否かが確認される（ステップＳ２０７）。最頻値グループΔｆｍ１の要素数が全体の半分以上である場合（ステップＳ２０７のＹｅｓ）、最頻値グループΔｆｍ１の各要素のピーク間隔の値の平均値（図１５では、Δｆｍ１の上にバーを付けて表記）が、該当する軸の周波数スペクトルのピーク間隔Δｆとして設定されて（ステップＳ２０８）、メインフローに戻る。
【００５８】
最頻値グループΔｆｍ１の要素数が全体の半分より少ない場合（ステップＳ２０７のＮｏ）、最頻値グループΔｆｍ１の次に要素数の多いグループが、２番目の最頻値グループΔｆｍ２として取得される（ステップＳ２０９）。
【００５９】
最頻値グループΔｆｍ１と、２番目の最頻値グループΔｆｍ２とを残し、これ以外のグループに属するピーク間隔に関するピークが削除される（ステップＳ２１０）。例えば、ピーク間隔Δｆ３（＝ｆ４−ｆ３）が、グループΔｆｍ１、Δｆｍ２のいずれにも属していない場合には、ピークｆ３が削除される。
【００６０】
ステップＳ２１０の処理が完了すると、ステップＳ２０２に戻り、ピーク間隔が再度取得されて上記の処理が繰り返される。
【００６１】
（角速度の推定）
図３に戻り、ステップＳ１０４では、角速度分離部６４により、ステップＳ１０３にて算出されたｘ，ｙ，ｚ軸の周波数スペクトルのピーク間隔Δｆｘ，Δｆｙ，Δｆｚを用いて、ｘ，ｙ，ｚ軸の角速度ωｘ，ωｙ，ωｚが算出される。ここで、ｘ，ｙ，ｚ軸の周波数スペクトルのピーク間隔Δｆｘ，Δｆｙ，Δｆｚは、球体１０を各軸方向から視たときの回転速度と対応する情報と考えることができる。このとき回転速度は、視点となる軸を除く他の２軸まわりの回転の角速度を合成したものである。例えば、ｘ軸の周波数スペクトルのピーク間隔Δｆｘは、球体１０の回転をｙｚ平面上に射影したときの、ｙｚ平面上の回転による回転速度に対応する情報である。このときの回転速度は、ｙ軸まわりの回転による角速度ωｙと、ｚ軸まわりの回転による角速度ωｚとを合成したものである。したがって、ｘ，ｙ，ｚ軸の周波数スペクトルのピーク間隔Δｆｘ，Δｆｙ，Δｆｚと、ｘ，ｙ，ｚ軸の角速度ωｘ，ωｙ，ωｚとは、以下の（２）〜（４）式の関係にある。このパラメータを用いることによって角速度を分離することができる。
【数４】
2014182032.tif 000006
【００６２】
上記の（２）〜（４）式を展開すると、角速度ωｘ，ωｙ，ωｚに関する（５）〜（７）式を導出できる。
【数５】
2014182032.tif 000007
【００６３】
角速度分離部６４は、ステップＳ１０３にて算出されたｘ，ｙ，ｚ軸の周波数スペクトルのピーク間隔Δｆｘ，Δｆｙ，Δｆｚを上記の（５）〜（７）式に代入することにより、ｘ，ｙ，ｚ軸の角速度ωｘ，ωｙ，ωｚを算出する。ステップＳ１０４の処理が完了するとステップＳ１０５に移行する。
【００６４】
（角速度の方向判定）
ステップＳ１０５では、角速度方向判定部６５により、ステップＳ１０４で算出された角速度ωｙ，ωｚの正負方向が判定される。角速度を導出する上記（５）〜（７）式には正負の曖昧性が残っている。角速度方向判定部６５は、球体１０の変化距離を用いて角速度の方向を特定する。空気中を飛翔する球体１０にはマグヌス効果や重力が働くため、回転方向によって飛翔軌道の変化方向が決まる。このときの軌道変化をｘ，ｙ，ｚ軸の三方向から同時に取得し（以下に記載する「軌道変化量Δｘ，Δｙ，Δｚ」）、球体１０の飛翔方向から予想される軌道変化（以下に記載する「基準軌道変化量ｘｍ，ｙｍ，ｚｍ」）と比較することによって角速度の正負を決定し、曖昧さを除去する。
【００６５】
図１７〜１９を参照して、角速度方向判定部６５における角速度方向判定手法を詳細に説明する。図１７は、ｘｙｚ座標系における球体１０の軌道変化を示す模式図である。上述のように、本実施形態では、球体１０は、ｘ軸の正方向に沿って飛翔を開始するものとする。ｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダ４は、球体１０の受信強度の測定と同時に球体１０の空間座標を取得する。図１７に示すように、このときの球体１０からの反射電磁波の測定開始点をｐ０（ｘ０，ｙ０，ｚ０）、マグヌス効果などの影響により球体１０の軌道が変化する前の任意の測定点をｐｈ（ｘｈ，ｙｈ，ｚｈ）、軌道変化後の測定終了点をｐｉ（ｘｉ，ｙｉ，ｚｉ）とし、ｘ，ｙ，ｚ軸における測定開始点ｐ０から測定終了点ｐｉまでの距離（軌道変化量）をそれぞれΔｘ，Δｙ，Δｚとする。
【００６６】
図１８は、図１７中の球体１０の軌道をｚｘ平面座標に投射した図である。図１８に示されるｚｘ平面座標での球体１０の軌道変化に基づいて、ｙ軸の角速度ωｙの正負を決定する。球体１０の回転による軌道変化が発生しないと仮定した場合の球体１０のｚ軸方向の軌道変化量ｚｍ（以下「基準軌道変化量」とも記載する）は、飛翔方向による軌道変化量ｚｓ、重力落下による軌道変化量ｚｆ、風力による軌道変化量ｚｗの和となる。この基準軌道変化量ｚｍと、実測された軌道変化量Δｚとを比較することによって、ｙ軸の角速度ωｙの正負を決定する。
【００６７】
ｚ軸の飛翔方向による軌道変化量ｚｓは、次のように算出する。まず、測定開始点ｐ０と、軌道変化前の任意の測定点Ｐｈと取得し、それらの差より測定点Ｐｈにおけるｚ方向の変化距離ｚｓ´を得る。球体１０が測定開始点ｐ０と測定点ｐｈの２点を通る等速直線運動を行なうとすると、球体１０の飛翔方向による軌道変化量ｚｓは、下記の（８）式で算出することができる。
【数６】
2014182032.tif 000008
【００６８】
これに加え、球体１０が空気中を飛翔する際には重力による落下が発生する。この重力による軌道変化量ｚｆは、重力加速度をｇ、球体１０の飛翔時間をｔとして下記の（９）式で算出することができる。
【数７】
2014182032.tif 000009
【００６９】
また、ｘ軸方向の軌道変化量Δｘは、下記の（１０）式で算出することができる。
【数８】
2014182032.tif 000010
ここで、ｖ０は初速度、ｗｘはｘ軸方向の風速を表す。なお、初速度ｖ０の導出には、ドップラーレーダを用いてドップラーシフトから算出する手法、または球体１０の軌道変化とその所要時間から算出する手法を用いることができる。
【００７０】
上記の（９）式及び（１０）式より、重力による軌道変化量Ｚｆは、下記の（１１）式で算出できる。
【数９】
2014182032.tif 000011
【００７１】
次に風による軌道変化量ｚｗを求める。実環境、特に屋外での測定の場合、風による軌道変化が発生する。この時のｚ軸方向の風速をｗｚとすると、この風力による軌道変化量ｚｗは、下記の（１２）式で算出できる。
【数１０】
2014182032.tif 000012
【００７２】
また、球体１０の飛翔時間ｔを、上記（１０）式に基づいて初速度ｖ０、ｘ軸方向の風力ｗｘ、ｘ軸方向の軌道変化量Δｘを用いて表し、（１２）式に代入すると、風力による軌道変化量ｚｗは、次の（１３）式で表すことができる。
【数１１】
2014182032.tif 000013
【００７３】
球体１０のｚ軸方向の基準軌道変化量ｚｍは、上述のとおり、飛翔方向による軌道変化量ｚｓ、重力落下による軌道変化量ｚｆ、及び風力による軌道変化量ｚｗの和であるため、上記の（８），（１１），（１３）式に基づき、下記の（１４）式で表すことができる。
【数１２】
2014182032.tif 000014
【００７４】
このように算出される球体１０のｚ軸方向の基準軌道変化量ｚｍが、球体１０のｙ軸の角速度ωｙの正負を決定するための閾値となり、このｚｍと、実測されたｚ軸方向の軌道変化量Δｚとを比較することによって、ｙ軸の角速度ωｙの方向を決定することができる。
【００７５】
軌道変化量Δｚが基準軌道変化量ｚｍより大きい場合（Δｚ＞ｚｍ）、すなわち図１８では球体１０がｚｍの下端より上方にある場合、球体１０はマグヌス力によって揚力を受けている状態にある。この力を受ける条件はバックスピンであり、ｙ軸の正方向側（図１８の紙面奥側）から視た場合、球体１０は時計回りに回転している。したがって、Δｚ＞ｚｍの場合、ｙ軸の角速度ωｙは負方向となる。
【００７６】
一方、軌道変化量Δｚが基準軌道変化量ｚｍより小さい場合（Δｚ＜ｚｍ）、すなわち図１８に示すように球体１０がｚｍの下端より下方にある場合、Δｚはｚｍと回転による下方向の変化の和となる。下方向の変化が発生する条件はトップスピンであり、ｙ軸の正方向側から視た場合、球体１０は反時計回りの回転となる。したがって、Δｚ＜ｚｍの場合、ｙ軸の角速度ωｙは正方向となる。以上より、ｙ軸の角速度ωｙの正負は、下記の（１５）式で決定することができる。
【数１３】
2014182032.tif 000015
【００７７】
図１９は、図１７中の球体１０の軌道をｘｙ平面座標に投射した図である。図１９に示されるｘｙ平面座標における球体１０の軌道変化に基づいて、ｚ軸の角速度ωｚの正負を決定する。球体１０の回転による軌道変化が発生しないと仮定した場合の球体１０のｙ軸方向の軌道変化量ｙｍ（以下「基準軌道変化量」とも記載する）は、飛翔方向による軌道変化量ｙｓと、風力による軌道変化量ｙｗの和となる。この基準軌道変化量ｙｍと、実測された軌道変化量Δｙとを比較をすることによって、ｚ軸の角速度ωｚの正負を決定する。
【００７８】
ｙ軸の飛翔方向による軌道変化量ｙｓは、次のように算出する。まず、測定開始点ｐ０と、軌道変化前の任意の測定点ｐｈとを取得し、それらの差より測定点ｐｈにおけるｙ方向の変化距離ｙｓ´を得る。このとき、球体１０の軌道が測定開始点ｐ０と測定点ｐｈの２点を通る延長線上とすると、球体１０の飛翔方向による軌道変化量ｙｓは、下記の（１６）式で算出できる。
【数１４】
2014182032.tif 000016
【００７９】
次に、風による軌道変化量ｙｗを求める。ｙ軸方向の風速をｗｙとすると、この風力による軌道変化量ｙｗは、下記の（１７）式で算出できる。
【数１５】
2014182032.tif 000017
【００８０】
また、球体１０の飛翔時間ｔを、初速度ｖ０、ｘ軸方向の風力ｗｘ、ｘ軸方向の軌道変化量Δｘを用いて表し、（１７）式に代入すると、風力による軌道変化量ｙｗは、次の（１８）式で表すことができる。
【数１６】
2014182032.tif 000018
【００８１】
球体のｙ軸方向の基準軌道変化量ｙｍは、上述のとおり、飛翔方向による軌道変化量ｙｓ、及び風力による軌道変化量ｙｗの和であるため、上記の（１６），（１８）式に基づき、下記の（１９）式で表すことができる。
【数１７】
2014182032.tif 000019
【００８２】
このように算出される球体１０のｙ軸方向の基準軌道変化量ｙｍが、球体１０のｚ軸の角速度ωｚの正負を決定するための閾値となり、このｙｍと、実測されたｙ軸方向の軌道変化量Δｙとを比較することによって、ｚ軸の角速度ωｚの方向を決定することができる。
【００８３】
軌道変化量Δｙが基準軌道変化量ｙｍより大きい場合（Δｙ＞ｙｍ）、すなわち図１９に示すように球体１０がｙｍの左端より左側にある場合、ｚ軸の正方向側（図１９の紙面手前側）から視た場合、球体１０は反時計回りに回転している。したがって、Δｙ＞ｙｍの場合、ｚ軸の角速度ωｚは正方向となる。
【００８４】
一方、軌道変化量Δｙが基準軌道変化量ｙｍより小さい場合（Δｙ＜ｙｍ）、すなわち図１９では球体１０がｙｍの左端より右側にある場合、ｚ軸の正方向側から視た場合、球体１０は時計回りに回転している。したがって、Δｙ＜ｙｍの場合、ｚ軸の角速度ωｚは負方向となる。以上より、ｚ軸の角速度ωｚの正負は、下記の（２０）式で決定することができる。
【数１８】
2014182032.tif 000020
【００８５】
なお、進行方向軸であるｘ軸の角速度ωｘは、正負どちらの場合でも球体１０は同様の軌道変化をする傾向にある。このため、本実施形態では、角速度方向判定部６５は、ｘ軸の角速度ωｘについては任意の方向を選択する。ステップＳ１０５の処理が完了するとステップＳ１０６に移行する。
【００８６】
（角速度ベクトルの導出）
ステップＳ１０６では、角速度ベクトル導出部６６により、ステップＳ１０４、Ｓ１０５で算出された各軸の角速度ωｘ，ωｙ，ωｚを用いて、球体１０の角速度ベクトルωが出力される。角速度ベクトルωは、以下の（２１）式で表すことができる。
【数１９】
2014182032.tif 000021
ここで、ｉ，ｊ，ｋは、それぞれｘ，ｙ，ｚ軸の単位ベクトルである。角速度ベクトルの概略を図２０に示す。図２０に示す角速度ベクトルωのベクトル量（長さ）が球体１０の回転速度、角速度ベクトルωの向きが球体１０の回転軸方向となる。ステップＳ１０６の処理が完了すると、本制御フローを終了する。
【００８７】
（角速度ベクトルの利用）
なお、計算装置６は、角速度ベクトルωを導出した後、この導出した角速度ベクトルωの情報を記憶装置７や出力装置８に出力し、記憶装置７に記憶したり、出力装置８に出力できる。また、出力装置８に、球体１０の回転速度や回転軸方向を文字情報や図形情報で表示することもできる。
【００８８】
次に、本実施形態に係る回転推定システム１の効果を説明する。
【００８９】
本実施形態の回転推定システム１は、回転を伴って飛翔する球体１０の回転速度及び回転軸位置を推定する。回転推定システム１は、球体１０に対してｘ，ｙ，ｚ軸の方向から電磁波１１，１３，１５を照射し、球体１０からの反射電磁波１２，１４，１６をｘ，ｙ，ｚ軸で受信して、ｘ，ｙ，ｚ軸の反射電磁波の受信強度に関する時系列情報を計測する測定装置５と、測定装置５により計測されたｘ，ｙ，ｚ軸の時系列情報を周波数情報にそれぞれ変換する計測装置６の周波数領域変換部６１と、周波数領域変換部６１により変換されたｘ，ｙ，ｚ軸の周波数情報におけるピーク間隔Δｆｘ，Δｆｙ，Δｆｚをそれぞれ特定する閾値処理演算部６２と、特定されたｘ，ｙ，ｚ軸のピーク間隔Δｆｘ，Δｆｙ，Δｆｚに基づいて、球体１０のｘ，ｙ，ｚ軸の角速度ωｘ，ωｙ，ωｚを推定する角速度分離部６４と、角速度分離部６４により推定されたｘ，ｙ，ｚ軸の角速度ωｘ，ωｙ，ωｚを成分とする角速度ベクトルωの大きさを球体１０の回転速度として、また、角速度ベクトルωの方向を球体１０の回転軸位置として出力する角速度ベクトル導出部６６と、を備える。
【００９０】
この構成により、回転を伴って飛翔する球体１０の回転位相変化に関する情報（反射電磁波の受信強度に関する時系列情報）を三方向から多角的に取得できるので、実際の球体１０の回転挙動について詳細かつ高精度に把握することできる。また、これらの多角的な情報に基づき球体１０の回転速度及び回転軸位置を推定するので、球体１０の実際の回転挙動に好適に則した回転パラメータを推定することが可能となり、球体１０の回転パラメータを精度良く推定できる。
【００９１】
また、本実施形態の回転推定システム１において、角速度分離部６４は、三軸のピーク間隔Δｆｘ，Δｆｙ，Δｆｚのそれぞれが、各軸方向から球体１０を視た場合の、視点となる軸を除く他の二軸まわりの回転の角速度を合成した回転速度に対応する情報であるとの特性に基づいて、球体１０の三軸の角速度ωｘ，ωｙ，ωｚを推定する。より詳細には、角速度分離部６４は、ｘ，ｙ，ｚ軸のピーク間隔をそれぞれΔｆｘ，Δｆｙ，Δｆｚ，ｘ，ｙ，ｚ軸の角速度をそれぞれωｘ，ωｙ，ωｚと表すとき、上記の（５）〜（７）式を用いて三軸の角速度を推定する。
【００９２】
上述のように、回転飛翔する球体１０のｘ，ｙ，ｚ軸の受信強度の周波数スペクトルは、ピーク間隔が各軸における回転速度と一致するような、周期的なピークを有する特徴がある。したがって、ｘ，ｙ，ｚ軸の受信強度の周波数スペクトルのピーク間隔Δｆｘ，Δｆｙ，Δｆｚを利用した上記の（５）〜（７）式を用いることで、球体１０のｘ，ｙ，ｚ軸の角速度を精度良く推定することができる。
【００９３】
また、本実施形態の回転推定システム１は、球体１０の軌道変化量Δｘ，Δｙ，Δｚをｘ，ｙ，ｚ軸の方向から取得し、球体１０の飛翔方向から予想される基準軌道変化量ｘｍ，ｙｍ，ｚｍと比較することによって、角速度分離部６４により推定された角速度ωｘ，ωｙ，ωｚの正負方向を判定する角速度方向判定部６５を備える。
【００９４】
上述のように、角速度分離部６４により上記の（５）〜（７）式を用いて算出される角速度ωｘ，ωｙ，ωｚは、正負両方を含むものであって曖昧性が残るものである。そこで、上記構成により球体１０の軌道変化に基づいて角速度の方向を特定することで、角速度分離部６４により算出された角速度の曖昧性を解消できる。これにより、最終的に導出される球体１０の回転パラメータの推定精度を向上できる。
【００９５】
また、本実施形態の回転推定システム１において、角速度方向判定部６５は、球体１０からの反射電磁波の測定開始点をＰ０（ｘ０，ｙ０，ｚ０）、軌道変化前の任意の測定点をＰｈ（ｘｈ，ｙｈ，ｚｈ）、軌道変化後の測定終了点をｐｉ（ｘｉ，ｙｉ，ｚｉ）、測定開始点ｐ０から測定終了点ｐｉまでのｘ，ｙ，ｚ軸の軌道変化量をΔｘ，Δｙ，Δｚ、重力加速度をｇ、初速度をｖ０、ｘ，ｙ，ｚ軸の風速をｗｘ，ｗｙ，ｗｚと表すとき、上記の（１５），（２０）式によって角速度の正負方向を判定する。
【００９６】
球体１０の回転方向が正負逆の場合、この回転に起因する球体１０の軌道変化は明確に異なるものとなる。上記構成により、球体の飛翔方向、重力落下、風力による軌道変化の影響を相殺して、球体１０の回転の影響によって発生した軌道変化量を切り出すことができる。したがって、球体１０の角速度ωｘ，ωｙ，ωｚの正負方向を精度良く判定できる。
【００９７】
また、本実施形態の回転推定システム１において、閾値処理演算部６２が、周波数領域変換部６１により変換されたｘ，ｙ，ｚ軸の周波数情報から所定の閾値以上のピークを抽出して、隣接するピーク間のピーク間隔を算出し、最頻値処理演算部６３が、閾値処理演算部６２により算出されたピーク間隔のうち頻度の高いものを抽出し、この抽出したピーク間隔に基づいてピーク間隔Δｆｘ、Δｆｙ、Δｆｚを特定する。
【００９８】
回転飛翔する球体１０に関するｘ、ｙ、ｚ軸の受信強度の周波数スペクトルは、ピーク間隔が各軸における回転速度と一致するような、周期的なピークを有するものである。しかし、この周波数スペクトルには、例えば反射電磁波の計測時のノイズなどの影響により、球体１０の回転とは関係のないピークが発生したり、本来発生する周波数のピークが得られない場合など、何らかの不確定要素が含まれる場合がある。このような不確定要素は最終的な回転パラメータの推定に悪影響を与える虞がある。そこで上記構成により、最頻値処理を施してピーク間隔Δｆｘ、Δｆｙ、Δｆｚを特定することで、ピーク間隔の情報から不確定要素を除外することが可能となる。これにより、球体１０の回転に関係する正しいピーク間隔を特定することが可能となり、回転パラメータの推定精度を向上できる。
【００９９】
本実施形態の回転推定システム１は、球技に応用することが考えられる。例えば、計測対象の球体１０を球技用ボールとすれば、角速度ベクトル導出部６６により出力された球技用ボールの回転速度及び回転軸位置に基づいて、回転軸位置からボールの球種を推定したり、回転速度から変化球の変化度合いを推定したり、ボールの起動を予測するなど、回転飛翔する球技用ボールの球種や軌道に関する情報を生成することができる。そして、このように生成した情報を、出力装置８や球技会場の表示装置に表示して、選手や観客に提示することもできる。
【０１００】
本実施形態の回転推定システム１を球技に適用することの効果をさらに説明する。本実施形態を球技に適用した場合の効果を、選手側と観客側の両方で説明する。
【０１０１】
まず、観客側の効果を説明する。試合に用いるボールは規格が決まっており、特に表面の加工が難しいために、従来は試合中に回転パラメータをリアルタイムで測定することが困難であったが、本実施形態の回転推定システム１は電磁波を用いて、受信強度の変化から角速度ベクトルを測定するため、無加工の球体でも回転パラメータをリアルタイムで測定することが可能となる。
【０１０２】
このようにリアルタイムで測定したボールの回転速度や回転軸位置などの情報を、会場のスクリーンに映し出して観客に情報提供することができれば、選手がどのような球を投げたのか、あるいは打っていたのかが確認できるようになり、臨場感の向上や戦術の推測等といったサービスを提供することができる。
【０１０３】
選手側の効果としては、リアルタイムで試合中に投げた、あるいは打った球の回転パラメータを取得することができるため、その場での技術修正を行うことができる。また、試合中に投げた球のパラメータをデータとして残すことが可能になるため、「重い球」、「伸びがある球」といった曖昧な定義だった球の解析、データを用いた弱点の改善やトレーニングの提案、好調時と不調時のパラメータを比べることによる調子の確認等といった技術の向上を図ることが出来、更なる人材育成につなげることができる。
【０１０４】
なお、球体の回転によって発生するドップラーシフトを利用して、球体の回転速度を導出する従来手法（特許文献１参照）が知られている。ドップラーシフトは、波の発生源（電磁波送信機）と、測定物（球体）との相対速度によって、反射波の周波数に変化が発生する現象である。この従来手法において、ドップラーシフトの発生要因は、球体の回転速度である。一方、本実施形態の回転推定システム１は、球体１０のレーダ断面積の変化によって発生する受信強度の時間変化に基づいて、球体１０の回転速度および回転軸を算出するものである。本実施形態において入力情報として用いる球体１０のレーダ断面積は、ドップラーシフトに影響を及ぼすものではないので、上記の従来手法では取得することができない情報である。このように、本実施形態の回転推定システム１は、ドップラーシフトを利用する従来手法と比較して、球体パラメータの導出に用いる入力情報が明らかに異なるものであり、従って、この入力情報に基づき球体パラメータを導出するための解析手法も異なるものである。
【実施例】
【０１０５】
次に、実施例に基づいて本発明をより具体的に説明する。図２１に、野球のカーブボールを３方向から測定した実施例の構成を示す。計測対象の球体１０として、特許第５１４２３６６号公報に記載の「遊戯用ボール」を使用した。送信電力は０（ｄＢｍ）、送信周波数は１０（ＧＨｚ）、測定時間は０．４０９６（秒）、測定ポイントは８１９２（点）、周波数分解能は２．４４（Ｈｚ）とし、ｘ軸レーダ２はバックネット、ｙ軸レーダ３は１塁側スタンド、ｚ軸レーダ４は天井に配置した。球体１０の初期位置（ピッチャーマウンド）とｘ軸レーダ２、ｙ軸レーダ３、ｚ軸レーダとの間の距離は、それぞれ３６．４（ｍ）、４０（ｍ）、５０（ｍ）であった。球体１０の各軸の回転速度は、すべて２０（ｒｐｓ）とした。風は吹いていないものとする。
【０１０６】
以上の条件下において取得した球体１０から反射電磁波に関する受信強度の時間変化を図２２，２３，２４に示す。図２２，２３，２４は、それぞれ取得されたｘ，ｙ，ｚ軸の受信強度の時間変化を示す図である。図２２〜２４の縦軸は受信強度（ｐＷ）を示し、横軸は計測時間（ｓ（秒））を示す。
【０１０７】
このとき得られた各軸の受信強度の時間変化をフーリエ変換すると、図２５，２６，２７で示される周波数スペクトルを得ることができる。図２５，２６，２７は、それぞれｘ、ｙ、ｚ軸の受信強度の周波数スペクトルを示す図である。図２５〜２７の縦軸はパワースペクトル（ｄＢ）を示し、横軸は周波数（Ｈｚ／ｒｐｓ）を示す。これらの周波数スペクトルにおいて、直流成分ピークを除いた最大値は、各軸いずれもｆ＝８０．５７（Ｈｚ／ｒｐｓ）の時である。よって、この範囲における平均値から閾値を求めると、図２５に示すｘ軸の周波数スペクトルでは、閾値は−１２３．０３（ｄＢ）となった。図２６に示すｙ軸の周波数スペクトルの場合、閾値は−１２８．９８（ｄＢ）となった。図２７に示すｚ軸の周波数スペクトルの場合、閾値は−１３２．７２（ｄＢ）となった。このように導出した閾値を用いて閾値処理を行うことによってピークとその間隔が得られる。この時得られたピーク間隔数は、図２５に示すｘ軸の周波数スペクトルでは１８個、図２６に示すｙ軸の周波数スペクトルでは２１個、図２７に示すｚ軸の周波数スペクトルでは２０個となった。
【０１０８】
次に、ピーク間隔のグループ分けを行った。図２５，２６，２７に示す各軸の周波数スペクトルにおける最小ピーク間隔は、すべて１９．５３（Ｈｚ／ｒｐｓ）となった。この値から周波数分解能２．４４（Ｈｚ）を２倍した値を加えた２４．４１（Ｈｚ／ｒｐｓ）までの値がグループ１となる。次にグループ１を除く最小ピーク間隔を求め、同様の処理を行った。この処理をすべての値がいずれかのグループに属するまで行った。
【０１０９】
この処理の結果、図２５，２７に示すｘ軸及びｚ軸の周波数スペクトルの場合、グループは、１９．５３  （Ｈｚ／ｒｐｓ）≦２４．４１（Ｈｚ／ｒｐｓ）のグループ１、５８．５９（Ｈｚ／ｒｐｓ）≦６３．４７（Ｈｚ／ｒｐｓ）のグループ２、７８．１３（Ｈｚ／ｒｐｓ）≦８３．０１（Ｈｚ／ｒｐｓ）のグループ３となった。また、図２６に示すｙ軸の周波数スペクトルの場合、グループは、１９．５３（Ｈｚ／ｒｐｓ）≦２４．４１（Ｈｚ／ｒｐｓ）のグループ１、４１．５０（Ｈｚ／ｒｐｓ）≦４６．３８（Ｈｚ／ｒｐｓ）のグループ２、６１．０４（Ｈｚ／ｒｐｓ）≦６５．９２（Ｈｚ／ｒｐｓ）のグループ３、７８．１３（Ｈｚ／ｒｐｓ）≦８３．０１（Ｈｚ／ｒｐｓ）のグループ４となった。
【０１１０】
グループ分けが完了したので、グループ内の要素の数をカウントし、最も数が多いグループと全体における割合を確認する。図２５に示すｘ軸の周波数スペクトルの場合、最も要素を多く含むグループはグループ３の９個であった。全要素数は１８個であるため、全要素数に対するグループ内の要素数の割合は５０％である。グループ内の要素が全体の半分以上を占めているので、このグループ内の要素の平均値を取得し、それをｘ軸のピーク間隔Δｆｘとする。グループ３の平均値は８０．０２（Ｈｚ／ｒｐｓ９であるため、Δｆｘ＝８０．０２（Ｈｚ／ｒｐｓ）となった。
【０１１１】
図２６に示すｙ軸の周波数スペクトルの場合、最も要素を多く含むグループはグループ１と４の７個であり、全要素数に対するグループ内の要素数の割合は４７．６２％となった。また、図２７に示すｚ軸の周波数スペクトルの場合、最も要素を多く含むグループはグループ２と３の７個であり、全要素数に対するグループ内の要素数の割合は３５％となった。これらは共にグループ内の要素数の割合が５０％を超えていないため、これらのグループ以外のピーク間隔Δｆｉの番号ｉに対応するピークｆｉを削除する。図２６に示すｙ軸の周波数スペクトルの場合、グループ２、３に属するピーク間隔は、Δｆ２，Δｆ４，Δｆ１０，Δｆ１２，Δｆ１５，Δｆ１７，Δｆ１９であるため、これらの番号に対応するピークｆ２，ｆ４，ｆ１０，ｆ１２，ｆ１５，ｆ１７，ｆ１９を削除した。同様に、図２７に示すｚ軸の周波数スペクトルの場合、グループ１に属するピーク間隔は、Δｆ２，Δｆ１０，Δｆ１２，Δｆ１３，Δｆ１５，Δｆ１７であるため、これらの番号に対応するピークｆ２，ｆ１０，ｆ１２，ｆ１３，ｆ１５，ｆ１７を削除した。
【０１１２】
次に、ピークを削除したｙ軸及びｚ軸の周波数スペクトルにおいて、ピーク番号ｉ及びピーク間隔Δｆｉを再取得した。図２８及び図２９は、ピーク間隔を再取得したｙ軸及びｚ軸の受信強度の周波数スペクトルを示す図である。図２８，２９の縦軸はパワースペクトル（ｄＢ）を示し、横軸は周波数（Ｈｚ／ｒｐｓ）を示す。
【０１１３】
図２８，２９に示すｙ軸及びｚ軸の受信強度の周波数スペクトルにおいて、ピーク間隔のグループ分けを再度行なった。図２８に示すｙ軸の周波数スペクトルの場合、ピーク間隔の全要素数は１５個、グループは、１９．５３（Ｈｚ／ｒｐｓ）≦２４．４１（Ｈｚ／ｒｐｓ）のグループ１、６１．０４（Ｈｚ／ｒｐｓ）≦６５．９２（Ｈｚ／ｒｐｓ）のグループ２、７８．１３（Ｈｚ／ｒｐｓ）≦８３．０１（Ｈｚ／ｒｐｓ）のグループ３、１３９．１６（Ｈｚ／ｒｐｓ）≦１４４．０１（Ｈｚ／ｒｐｓ）のグループ４となった。図２９に示すｚ軸の周波数スペクトルの場合、ピーク間隔の全要素数は１５個、グループは、５８．５９（Ｈｚ／ｒｐｓ）≦６３．４７（Ｈｚ／ｒｐｓ）のグループ１、７８．１３（Ｈｚ／ｒｐｓ）≦８３．０１（Ｈｚ／ｒｐｓ）のグループ２、１００．１０（Ｈｚ／ｒｐｓ）≦１０４．９８（Ｈｚ／ｒｐｓ）のグループ３となった。
【０１１４】
グループ分けが完了したので、グループ内の要素の数をカウントし、最も数が多いグループと全体における割合を確認する。図２８に示すｙ軸の周波数スペクトルの場合、最も要素を多く含むグループはグループ３の１１個なのに対し、全要素数は１５個であるため、全要素数に対するグループ内の要素数の割合は７３％となった。図２９に示すｚ軸の周波数スペクトルの場合、最も要素を多く含むグループはグループ２の１１個なのに対し、全要素数は１５個であるため、全要素数に対するグループ内の要素数の割合は７３％となった。これらのグループは共に要素数が全体の半分以上を占めているので、このグループ内の要素の平均値を取得し、これらをｙ軸及びｚ軸のピーク間隔Δｆｙ，Δｆｚする。図２８におけるグループ３の平均値は７９．５９（Ｈｚ／ｒｐｓ）、図２９におけるグループ２の平均値は７９．５９（Ｈｚ／ｒｐｓ）となるため、Δｆｙ＝７９．５９（Ｈｚ／ｒｐｓ）、Δｆｚ＝７９．５９（Ｈｚ／ｒｐｓ）となった。
【０１１５】
なお、この球体１０は１回転する際に４つの周期を有しているので、Δｆｘ，Δｆｙ，Δｆｚを４で割る必要がある。これによって各軸の周波数スペクトルのピーク間隔はΔｆｘ＝１９．９０（Ｈｚ／ｒｐｓ）、Δｆｙ＝１９．９０（Ｈｚ／ｒｐｓ）、Δｆｚ＝１９．９０（Ｈｚ／ｒｐｓ）となった。
【０１１６】
本実施例の測定条件は、上述のように球体１０の各軸の回転速度をすべて２０（ｒｐｓ）と設定した。上記計算により算出された各軸のピーク間隔Δｆｘ，Δｆｙ，Δｆｚを設定値（＝２０ｒｐｓ）と比較すると、設定値に対する推定誤差が０．５％と非常に小さく抑えられている。
【０１１７】
以上によって得られたピーク間隔Δｆｘ，Δｆｙ，Δｆｚから、（５）〜（７）式を用いて各軸の角速度ωｘ，ωｙ，ωｚに分離すると、ωｘ＝３９．７９（ｒａｄ／ｓ）、ωｙ＝３９．７９π（ｒａｄ／ｓ）、ωｚ＝３９．７９π（ｒａｄ／ｓ）となり、角速度ベクトルωは、
【数２０】
2014182032.tif 000022
となった。
【０１１８】
次に、球体１０の軌道変化に基づき曖昧性の除去を行った。図３０，３１，３２は、実施例における球体１０のｘ，ｙ，ｚ軸の軌道変化をそれぞれ示す図である。図３０〜３２の縦軸は球体１０の移動距離（ｍ）を示し、横軸は測定時間（ｓ（秒））を示す。図３０，３１，３２に示すように、この測定における測定開始点の座標をｐ０（０，０，０）とした時、各軸の軌道変化量は、Δｘ＝１３．５６（ｍ）、Δｙ＝０．２５（ｍ）、Δｚ＝−０．５（ｍ）となる。また、この時の測定点はｐ０＝０（ｓ）、ｐｈ＝０．１２（ｓ）、ｐｉ＝０．４０９６（ｓ）とした。なお、ｐｈにおけるｘ軸の軌道変化は４（ｍ）であるため、球体の初速度は１２０（ｋｍ／ｈ）となる。
【０１１９】
ｚｘ平面座標におけるｚｓ，ｚｆ，ｚｗは、それぞれ（８）式、（１１）式、（１３）式より、０．４０６８（ｍ）、−０．８１２５（ｍ）、０（ｍ）となった。これらより、球体１０の回転による変化が発生しない場合のｚ軸の軌道変化（基準軌道変化量）は、（１４）式より−０．４０５７（ｍ）となった。これに対し、Δｚ＝−０．５（ｍ）であるため、（１５）式よりωｙ＞０となった。
【０１２０】
ｘｙ平面座標におけるｙｓ，ｙｗは、（１６）式、（１８）式より−０．３４６５（ｍ）、０（ｍ）となった。これらより、球体１０の回転による変化が発生しない場合のｙ軸の軌道変化（基準軌道変化量）は、（１９）式より−０．４０５７（ｍ）となった。これに対してΔｙ＝０．２５（ｍ）であるため、（２０）式よりωｚ＞０となった。これより、角速度ベクトルωは、
【数２１】
2014182032.tif 000023
となった。
【０１２１】
以上の計算を実際にコンピュータプログラムで１０回動作させ、計算時間の最大値、最小値、平均値を計測した。なお、この計算に用いたコンピュータのＣＰＵは、Ｉｎｔｅｌ（登録商標）社のＣｏｒｅ（商標）２  Ｄｕｏ  Ｅ７３００（２．６６（ＧＨｚ））であり、メモリは２（ＧＢ）とした。プログラム言語には、Ｍａｔｌａｂ（登録商標）を使用して、この計算のためのコンピュータプログラムを作成した。このコンピュータ及びプログラム言語を用いた計算時間は、最大値０．６８７秒、最小値０．５７８秒、平均値０．６３１秒であり、いずれにおいても１秒以下で角速度ベクトルを算出することが出来た。なお，これらの計算時間は、マシンスペック、プログラム言語、プログラムの書き方に依存するため、上記の条件を変更することによって計算時間の短縮をさらに図ることができる。
【０１２２】
以上、本発明の実施形態を説明したが、上記実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。上記実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。上記実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【０１２３】
上記実施形態では、本発明に係る回転推定システム１が回転パラメータを推定する対象である回転飛翔体として球体１０を例示して説明したが、回転飛翔体は回転して飛翔する物体であればよく、例えば楕円体などの球体以外の形状の物体も回転飛翔体に含めることができる。また、上記実施形態では、回転飛翔体として、主に野球用ボールなど球技に用いるボールを想定していたが、例えば砲弾など球技以外の用途の物体も対象とすることができる。
【符号の説明】
【０１２４】
１  回転推定システム
５  測定装置（計測手段）
６  計算装置
６１  周波数領域変換部（変換手段）
６２  閾値処理演算部（ピーク間隔特定手段）
６３  最頻値処理演算部（ピーク間隔特定手段）
６４  角速度分離部（角速度推定手段）
６５  角速度方向判定部（判定手段）
６６  角速度ベクトル導出部（出力手段）
８  出力装置（表示手段）
１０  球体（回転飛翔体）
【図１】
2014182032.tif 000026
【図２】
2014182032.tif 000027
【図３】
2014182032.tif 000028
【図４】
2014182032.tif 000029
【図５】
2014182032.tif 000030
【図６】
2014182032.tif 000031
【図７】
2014182032.tif 000032
【図８】
2014182032.tif 000033
【図９】
2014182032.tif 000034
【図１０】
2014182032.tif 000035
【図１１】
2014182032.tif 000036
【図１２】
2014182032.tif 000037
【図１３】
2014182032.tif 000038
【図１４】
2014182032.tif 000039
【図１５】
2014182032.tif 000040
【図１６】
2014182032.tif 000041
【図１７】
2014182032.tif 000042
【図１８】
2014182032.tif 000043
【図１９】
2014182032.tif 000044
【図２０】
2014182032.tif 000045
【図２１】
2014182032.tif 000046
【図２２】
2014182032.tif 000047
【図２３】
2014182032.tif 000048
【図２４】
2014182032.tif 000049
【図２５】
2014182032.tif 000050
【図２６】
2014182032.tif 000051
【図２７】
2014182032.tif 000052
【図２８】
2014182032.tif 000053
【図２９】
2014182032.tif 000054
【図３０】
2014182032.tif 000055
【図３１】
2014182032.tif 000056
【図３２】
2014182032.tif 000057