JP3873721B2 - 周波数解析装置及び打検装置 - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明は、打検の反響音等の振動のモード周波数を、高速フーリエ変換法（Fast Fourier Transform：ＦＦＴ）等の離散フーリエ級数展開を用いて同定する周波数解析技術に関し、特に、ＦＦＴ等の離散フーリエ級数展開を実行する際に用いる窓関数として複素窓関数を用いた周波数解析技術に関する。
【０００２】
【従来の技術】
内容物が正常に充填・密封された缶詰においては、缶内圧力が陰圧等の所定の圧力に保たれている。これに対し、充填・密封が不完全な場合には缶内圧は大気圧と等しくなり、また内容物が腐敗・発酵してガスが発生すると缶内圧力が上昇する。このように、異常缶と正常缶とでは缶内圧力が一般に異なっている。そして、缶内圧力が異なれば、缶体の振動の固有振動数が異なる。そこで、従来から打検により缶内圧力を検査し、異常缶を排除している。
【０００３】
打検としては、一般に、インライン打検とケース打検とが実施されている。
インライン打検は、主に充填・密封の不良による漏洩を検知するため、充填・巻締・レトルト工程を経てカートン詰めされる直前の缶詰に対して実施される。
一方、ケース打検は、主に変敗による缶内圧力変化を検知するため、カートン詰め後、数日間の保存期間を経過した缶詰に対し、カートンに詰めたままで実施される。ケース打検は、カートンを開けずに検査できるので、異常が発見されて回収された製品を再選別する際に特に有用である。
【０００４】
打検の際には、一般にＦＦＴを用いて缶体の振動の固有振動数を求める。ＦＦＴとは、離散フーリエ級数展開を高速に実行するための計算アルゴリズムのことであり、その計算原理は、離散フーリエ級数展開とまったく同じである。
ＦＦＴでは、観測した振動波形が無限に繰り返されるとみなして観測データを解析する。ところが、データ列の先頭のデータと最後のデータとでは、通常、位相不整合が生じている。このため、ＦＦＴにより求めるスペクトルには誤差が生じる。
【０００５】
そこで、ＦＦＴによる周波数解析では、データ列に実数関数を乗じる前処理を行うことにより、位相不整合を解消し、スペクトルの誤差を低減させる。この前処理において、位相不整合を解消する目的で用いられる関数を窓関数と称する。窓関数としては、例えばハニング窓関数やハミング窓関数が知られているが、このほかにも解析対象とする振動の種類に応じて種々の実数関数が用いられる。
【０００６】
ところで、ＦＦＴをはじめとする離散フーリエ級数展開におけるスペクトルの周波数分解能は、直交関数系の基底モードの周波数、すなわち、全観測時間の逆数により制約される。この制約を相反則という。相反則における周波数分解能とは、近接する二つのモードの振動周波数を弁別する能力をいう。
【０００７】
全観測時間がＴの場合、ＦＦＴ等の離散フーリエ級数展開により求められるスペクトルの周波数間隔は１／Ｔである。二つのモードを弁別するためには、スペクトルが＜山−谷−山＞のパターンを有していることが必要である。このため、この場合の周波数分解能は２／Ｔとなる。したがって、振動時間が短い場合、高い周波数分解能を得ることは原理的に困難である。例えば、打検による振動のように振動時間がわずか数ミリ秒（ｍｓ）間しか続かない場合、周波数分解能は高々数百ヘルツ（Ｈｚ）にすぎない。
【０００８】
これに対し、打検による振動は、通常、振動モードが離散的であり、多くの場合、振動モードは単一モードである。このため、打検では、一つのモード周波数を決定する能力、言い換えれば、モード周波数を同定する能力が要求される。すなわち、打検の際の周波数解析において要求される周波数精度は、上述した周波数弁別能力とは異なるものである。このため、周波数同定能力については、上述した相反則は適用されない。したがって、観測時間の短い波形データからであっても、周波数弁別能力以上の精度でモード周波数を同定することは可能である。
【０００９】
そこで、従来の打検システムにおいては、高い周波数同定能力を得るため、サンプリングデータの後ろに「０」値を付加してサンプリング点を拡張し、見かけ上の観測時間を長くしてＦＦＴを行っている。例えば、サンプリング周波数２０ｋＨｚで１２８点をサンプリングした観測データに「０」値データを付加し、サンプリング点を２０４８点に拡張したうえでＦＦＴを行っている。この方法（以下、「０点付加法」とも称する。）により、１モード振動のモード周波数を十分な精度で同定することができる。その結果、例えば、１０Ｈｚ程度の精度で振動のモード周波数を求め、缶内圧力変化を検知することができる。
【００１０】
【発明が解決しようとする課題】
上述した０点付加法による周波数解析方法は、観測時間が短い場合であっても高い精度でモード周波数を同定することができる点で優れているが、サンプリング点を拡張した分だけ余分な計算を行う必要があり、演算負荷が増大してしまう点で改良の余地がある。
【００１１】
なお、高速演算速度を有するＣＰＵを採用したり、ＦＦＴ専用素子を利用することにより、演算処理速度の向上を図ることも考えられる。しかしながら、高速演算ＣＰＵは、発熱量の増大や、それに伴う動作の安定性限界等の問題を有している。また、ＦＦＴ専用素子を採用すると、打検装置の性能が特定の素子によって決まってしまう。その結果、打検装置の寿命は一般的な素子の生産期間よりも長いため、素子の生産中止により打検装置の仕様変更を余儀なくされることになる。
【００１２】
本発明は、上記の事情にかんがみてなされたものであり、演算負荷の増加を抑制しつつ、高い精度でモード周波数を同定することができる周波数解析技術の提供を目的とする。
【００１３】
【課題を解決するための手段】
上記第一の目的の達成を図るため、本出願に係る発明者は、種々の検討及び実験を重ねた結果、ＦＦＴ等の離散フーリエ級数展開により求めたフーリエスペクトルを複素三角関数で表示した場合、周波数原点が位相項として他の表示周波数と分離して表せることに着目した。そして、この位相項を前処理において窓関数として用いれば、スペクトルの表示周波数の原点を周波数軸上で任意にずらすことができることに想到した。すなわち、表示周波数原点をずらすことにより、任意の周波数のスペクトルを求めることができることに想到した。
【００１４】
そこで、本発明の請求項１記載の周波数解析装置は、振動現象から一定のサンプリング時間間隔でサンプリングされた波形データ列に、窓関数を乗じて切出波形データ列を生成する前処理手段と、前記切出波形データ列から高速フーリエ変換法等の離散フーリエ級数展開によりスペクトルデータ列を生成するフーリエ変換手段と、前記スペクトルデータ列に含まれるスペクトルデータのうち、最大値となるスペクトルデータの周波数を前記振動音のモード周波数として同定する同定手段とを備える周波数解析装置であって、前記前処理手段は、前記窓関数として、実数窓関数（実数関数の窓関数）に複素位相項を乗じた複素窓関数を使用し、前記複素窓関数の前記複素位相項の値を一定間隔で変化させた各複素窓関数を前記波形データ列にそれぞれ乗じることにより、互いに位相をずらした複数の切出波形データ列を生成し、前記フーリエ変換手段は、前記各切出波形データ列から、周波数原点を一定周波数ずつずらした複数のスペクトルデータ列を追加生成する構成としてある。
【００１５】
ＦＦＴ等の離散フーリエ級数展開により求めたフーリエスペクトルを複素三角関数で表示した場合、周波数原点を位相項として他の表示周波数と分離して表示することができる。このため、ＦＦＴ等の離散フーリエ級数展開を行う前に位相項だけを計算することができる。すなわち、位相項を前処理において複素数の窓関数として用いることができる。その結果、複素窓関数は、単に、観測データ列の端部どうしの位相不整合を低減する目的にとどまらず、表示周波数の原点を周波数軸上で移動させる目的でも利用することができる。
【００１６】
そして、本発明の周波数解析装置によれば、複素窓関数の複素位相項の値を変化させ、スペクトルの周波数原点を周波数軸上で移動させる。その結果、任意の周波数を原点とするスペクトルデータ列を生成することができる。このため、スペクトルの表示周波数間隔は原理的に無限に細かくすることができる。すなわち、周波数同定能力は、原理的に制約を受けない。これにより、短時間の波形データであっても高い精度で打検のモード周波数を同定することができる。
【００１７】
また、本発明によれば、通常のＦＦＴに加え、周波数原点をシフトさせた切出波形データに対してもＦＦＴ等の離散フーリエ級数展開を行うため、通常のＦＦＴに比べて演算量は増加する。しかしながら、上述した０点付加法と比較すると、より少ない演算量で同程度の精度でモード周波数を同定することができる。このため、０点付加法を用いた場合に比べて、本発明では演算負荷の増加を抑制しつつ、高精度でモード周波数を同定することができる。
【００１８】
また、前処理手段は、複素窓関数の複素位相項の値を一定間隔で変化させた各複素窓関数を波形データ列にそれぞれ乗じることにより、互いに位相をずらした複数の切出波形データ列を生成し、フーリエ変換手段は、各切出波形データ列から、周波数原点を一定周波数ずつずらした複数のスペクトルデータ列を追加生成する構成としてあるので、上述したように、周波数原点をずらしたスペクトルデータ列を生成すれば、通常のＦＦＴによりスペクトルデータ列を生成する場合に比べ、演算量が増加する。しかしながら、上述した０点付加法と比較すると、より少ない演算量で同程度の精度でモード周波数を同定することができる。このため、０点付加法を用いた場合に比べて、本発明では演算負荷の軽減を図ることができる。その結果、処理時間の短縮を図ることができる。
【００１９】
また、請求項２記載の発明は、振動現象から一定のサンプリング時間間隔でサンプリングされた波形データ列に、窓関数を乗じて切出波形データ列を生成する前処理手段と、前記切出波形データ列から高速フーリエ変換法等の離散フーリエ級数展開によりスペクトルデータ列を生成するフーリエ変換手段と、前記スペクトルデータ列に含まれるスペクトルデータのうち、最大値となるスペクトルデータの周波数を前記振動音のモード周波数として同定する同定手段とを備える周波数解析装置であって、前記同定手段は、周波数原点をずらしていない基礎スペクトルデータ列のうち、最大値となるスペクトルデータの周波数を中央周波数とし、前記前処理手段は、前記中央周波数を前記基礎スペクトルデータ列の周波数間隔の半分ずつ前後にずらした候補区間内の周波数へ前記基礎スペクトルデータ列の周波数原点を移動させる複素窓関数を、前記波形データ列に乗じて切出波形データ列を生成し、前記フーリエ変換手段は、前記切出波形データ列から生成されるスペクトルデータ列のうち、周波数原点のスペクトルデータを選択的に追加生成し、前記同定手段は、中央周波数のスペクトルデータ及び追加生成されたスペクトルデータのうちで最大値となるスペクトルデータの周波数をモード周波数として同定する構成としてある。
【００２０】
このようにすれば、通常のＦＦＴによるスペクトルデータ列を計算した後、そのうちの最大値のスペクトルの周波数付近のみで周波数間隔を分割してモード周波数を求めることができる。さらに、候補区間内の周波数へ周波数原点を移動させるので、追加生成されるスペクトルデータ列のうち、周波数原点に対応するスペクトルデータのみを計算してモード周波数を求めることができる。これにより、演算負荷のより一層の軽減を図ることができる。その結果、処理時間の一層の短縮を図ることができる。
【００２１】
また、請求項３記載の発明は、振動現象から一定のサンプリング時間間隔でサンプリングされた波形データ列に、窓関数を乗じて切出波形データ列を生成する前処理手段と、前記切出波形データ列から高速フーリエ変換法等の離散フーリエ級数展開によりスペクトルデータ列を生成するフーリエ変換手段と、前記スペクトルデータ列に含まれるスペクトルデータのうち、最大値となるスペクトルデータの周波数を前記振動音のモード周波数として同定する同定手段とを備える周波数解析装置であって、（ａ）前記同定手段は、周波数原点をずらしていない基礎スペクトルデータ列のうち、最大値となるスペクトルデータの周波数を基準周波数とし、（ｂ）前記前処理手段は、前記基準周波数を前記基礎スペクトルデータ列の周波数間隔の半分である分割周波数間隔ずつ前後にずらした二つの周波数へ周波数原点をそれぞれ移動させる複素窓関数を、前記波形データ列にそれぞれ乗じて一対の切出波形データ列をそれぞれ生成し、（ｃ）前記フーリエ変換手段は、前記一対の切出波形データ列からそれぞれ生成されるスペクトルデータ列のうち、二つの周波数原点のスペクトルデータを選択的に追加生成し、（ｄ）前記同定手段は、追加生成された二つのスペクトルデータ及び前記基準周波数のスペクトルデータのうち、最大値をとるスペクトルデータの周波数を新たな基準周波数とし、（ｅ）前記前処理手段は、最新の基準周波数を直前の分割周波数間隔のさらに半分の新たな分割周波数ずつ前後にずらした二つの周波数へ周波数原点をそれぞれ移動させる複素窓関数を、前記波形データ列にそれぞれ乗じて一対の新たな切出波形データ列を生成し、（ｆ）前記フーリエ変換手段は、最新の一対の切出波形データ列からそれぞれ生成されるスペクトルデータ列のうち、二つの周波数原点のスペクトルデータを選択的に追加生成し、（ｇ）前記同定手段は、追加生成された最新の二つのスペクトルデータ及び最新の基準周波数のスペクトルデータのうちから、最大値をとるスペクトルデータの周波数を新たな基準周波数とし、前記前処理手段、フーリエ変換手段及び同定手段は、前記（ｅ）〜（ｇ）の処理を繰り返した後、最新の基準周波数を前記モード周波数として同定する構成としてある。
【００２２】
このようにすれば、通常のＦＦＴ等の離散フーリエ級数展開によるスペクトルデータ列を計算した後、最大値のスペクトルの周波数付近のスペクトルのみを、その範囲を徐々に狭めながら選択的に計算してモード周波数を求めることができる。さらに、いったん基準周波数を求めた後は、最新の基準周波数を直前の分割周波数間隔のさらに半分の新たな分割周波数ずつ前後にずらした二つの周波数へ周波数原点をそれぞれ移動させるので、追加生成されるスペクトルデータ列のうち、周波数原点に対応するスペクトルデータのみを計算してモード周波数を求めることができる。これにより、演算負荷のより一層の軽減を図ることができる。その結果、処理時間の一層の短縮を図ることができる。
【００２３】
また、請求項４記載の発明によれば、前処理手段は、複素窓関数として、下記の（１）式に示すＣ（ｎ，ｍ）を含む複素関数を使用する構成としてある。
Ｃ（ｎ，ｍ）＝Ｈ（ｎ）・ｅｘｐ（−ｉ２π（ｎ・ｍ）／（ＮＭ））…（１）
ただし、ｎ＝０，１，…，Ｎ−１を満たす。Ｎは、波形データのサンプリング数を表し、１以上の整数である。また、ｍ＝０，１，…，Ｍ−１を満たす。Ｍは、スペクトルデータ列の周波数間隔の分割数を表し、１以上の整数である。また、Ｈ（ｎ）は、実数関数の窓関数を表す。
【００２４】
このように、従来の実数窓関数Ｈ（ｎ）に、複素位相項ｅｘｐ（−ｉ２π（ｎ・ｍ）／（Ｎ／Ｍ））を乗じた複素窓関数Ｃ（ｎ，ｍ）を用いることが望ましい。波形データ列｛ｘ（ｎ）｝に、この複素窓関数Ｃ（ｎ，ｍ）を乗じると、切出波形データ列｛Ｃ（ｎ，ｍ）・ｘ（ｎ）｝が得られる。この切出波形データ列｛Ｃ（ｎ，ｍ）・ｘ（ｎ）｝に対して離散フーリエ級数展開を行って得たスペクトルデータ列｛Ｘ（ｋ，ｍ）｝の周波数ｆ_k,mは、下記の（３）式で与えられる。
ｆ_k,m＝（ｆ_s／Ｎ）（ｋ＋ｍ／Ｍ）…（３）
【００２５】
上記の（３）式に示すように、Ｎ点のスペクトルデータ列｛Ｘ（ｋ，ｍ）｝の周波数間隔は、（ｆ_s／Ｎ）である。さらに、ｍ＝０，１，…，Ｍ−１と変化させることにより、スペクトルデータ列の各周波数を（ｆ_s／（Ｍ・Ｎ））ずつずらすことができる。「Ｍ」の値は無制限に大きくすることができるので、（ｆ_s／（Ｍ・Ｎ））の値も無制限に小さくすることができる。すなわち、複素窓関数を用いることにより、単一モード、又は、他と十分に独立したモードの振動波形のモード周波数は、原理的にいくらでも精度よく求めることができる。
【００２６】
また、請求項５記載の発明によれば、前処理手段は、複素窓関数として、Ｃ（ｎ，ｍ）に、第二の複素位相項ｅｘｐ（−ｉ２πｎｋ／Ｎ）を乗じた下記の（２）式に示すＣ_k（ｎ，ｍ）を使用する構成としてある。
Ｃ_k（ｎ，ｍ）＝Ｃ（ｎ，ｍ）・ｅｘｐ（−ｉ２πｎｋ／Ｎ）…（２）
ただし、ｋは、スペクトルデータ列に含まれるスペクトルデータの、周波数原点のスペクトルデータから数えた番号を表し、０以上の整数である。
【００２７】
上記の（１）式に示すＣ（ｎ，ｍ）により、周波数間隔内で周波数原点を移動させることができる。これに加えて、上記の（２）式に示すＣ_k（ｎ，ｍ）によれば、第二複素位相項により、周波数間隔を越えて、任意の周波数は周波数原点を移動させることができる。
【００２８】
また、本発明の請求項６記載の打検装置は、被検査体に衝撃を与える衝撃手段と、打撃を与えられた前記被検査体の振動音を集音して電気的な反響信号を生成する集音手段と、前記反響信号をアナログ／デジタル変換し、一定のサンプリング時間間隔でサンプリングされた波形データ列を生成するデジタル変換手段と、同定されたモード周波数に基づいて被検査体の良否を判定する判定手段とを備える打検装置であって、前記請求項１〜５のいずれかに記載の周波数解析装置を備え、当該周波数解析装置によって前記モード周波数を同定する構成としてある。
【００２９】
このように、本発明の打検装置は、上述の周波数解析装置を利用したものであり、前処理手段が、複素窓関数の複素位相項の値を変化させ、スペクトルの周波数原点を周波数軸上で移動させる。その結果、フーリエ変換手段は、任意の周波数を原点とするスペクトルデータ列を生成することができる。このため、フーリエ変換手段により生成されるスペクトルの表示周波数間隔は原理的に無限に細かくすることができる。すなわち、周波数同定能力は、原理的に制約を受けない。これにより、高い精度でモード周波数を同定することができる。
【００３０】
また、本発明によれば、通常のＦＦＴに加え、周波数原点をシフトさせた切出波形データに対してもＦＦＴ等の離散フーリエ級数展開を行うため、通常のＦＦＴに比べて演算量は増加する。しかしながら、上述した０点付加法と比較すると、より少ない演算量で同程度の精度でモード周波数を同定することができる。このため、０点付加法を用いた場合に比べて、本発明では演算負荷の増加を抑制しつつ、高精度でモード周波数を同定することができる。
【００３１】
また、請求項７記載の発明は、フーリエ変換手段が、スペクトルデータ列に含まれるスペクトルデータのうち、１０００〜３０００Ｈｚの周波数範囲の少なくとも一部分を含む特定の周波数範囲内のスペクトルデータのみを選択的に生成する構成としてある。
【００３２】
通常のＦＦＴにおいては、サンプリング周波数が２０ｋＨｚの場合、０〜１００００Ｈｚの周波数範囲にわたって周波数解析を行っている。これに対して、打検による振動のモード周波数は、１０００Ｈｚ〜３０００Ｈｚの範囲内に含まれることが多い。このため、この周波数範囲の少なくとも一部分を含む特定の周波数範囲内のスペクトルデータを選択的に計算すれば、演算負荷を一層軽減し、処理速度の向上を図ることができる。
【００３３】
【発明の実施の形態】
以下、本発明の周波数解析装置及び打検装置の実施の形態について、図面を参照して説明する。本発明の各実施形態の説明に先立ち、本発明の理解を容易にするため、離散的に観測されたデータ列に対する、ＦＦＴ等の離散フーリエ級数展開をはじめとする離散フーリエ級数展開による離散的スペクトルの推定方法について説明する。
【００３４】
サンプリング間隔τ＝１／ｆ_s（ｆ_sはサンプリング周波数）でＮ回観測された波形データ列｛ｘ（ｔ_ｎ）：ｎ＝０，１，２，…，Ｎ―１｝を、関数系列｛φ_k（ｔ）：ｋ＝０，１，２，…｝を用いて最小二乗近似する。すなわち、二乗ノルムは、下記の（４）式で与えられる。
【００３５】
【数１】

【００３６】
上記二乗ノルムの値が最小となるように係数ｃ_kを決定することにより、関数系列｛φ_k（ｔ）：ｋ＝０，１，２，…｝の一次結合として、データ列｛ｘ（ｔ_ｎ）：ｎ＝０，１，２，…，Ｎ―１｝を最良近似する。関数系列を適切に選べば、最良近似は一通りに決まる。
【００３７】
関数系列を下記の（５）式で与えられる波動方程式の解である複素三角関数系｛ｅｘｐ｛ｉ２πｆ_kｔ｝：ｋ＝０，１，…｝で構成し、その係数をＸ（ｋ）とすると、Ｘ（ｋ）は振動周波数ｆ_kに対する振動振幅を表す。これをスペクトルと称する。
（ｄ²ｘ／ｄｔ²）＋ω²ｘ＝０ …（５）
【００３８】
関数系列｛φ_k（ｔ）｝を下記の（６）式で与えられる正規化直交条件を満たすように定めると、この関数系列は離散フーリエ級数展開となる。この正規直交系列に対する係数ｃ_kをフーリエ係数と称する。フーリエ係数ｃ_kは、上記の（４）式を各フーリエ係数で偏微分してその値を０とおくことにより、下記の（７）式で与えられる。
【００３９】
【数２】

【００４０】
さらに、関数系列｛φ_k（ｔ）｝が複素関数系｛ｅｘｐ（ｉ２πｆ_kｔ）：ｋ＝０，１，２，…｝であり、波形データ列｛ｘ（ｔ_ｎ）｝が、ｔ_ｎ＝ｎ・τの等間隔サンプリングである場合、上記の（６）式に示した正規直交系列条件式は、下記の（８）式で与えられる。
【００４１】
【数３】

【００４２】
これにより、周波数に関する直交条件は、下記の（９）式で与えられる。
ｆ_k＝ｆ₀＋（ｋ／Ｎτ）＝ｆ₀＋（ｋ／Ｔ）＝ｆ₀＋（ｋ／Ｎ）ｆ_s …（９）
また、フーリエスペクトル（スペクトルデータ）Ｘ（ｋ）は、下記の（１０）式で与えられる。
【００４３】
【数４】

【００４４】
ただし、上記の（９）及び（１０）式において、Ｔ＝Ｎτは全観測時間を表す。また、複素三角関数の周期性からｋは、ｋ＝０，１，…，Ｎ−１に限定される。
【００４５】
そして、上記の（９）式に示した周波数に関する直交条件から相反則が導き出される。すなわち、フーリエスペクトルの周波数間隔Δｆ_kは、全観測時間Ｔの逆数で定められることが分かる。
【００４６】
また、サンプリング周波数ｆ_sで標本化した波形データからは、最大ｆ_s／２までの周波数のスペクトルデータしか得られない。その理由は、ｆ_s／２よりも高い周波数成分を含む振動をサンプリングすると、エリアジングによる誤差を招くためである。したがって、スペクトルデータの有効な周波数は、下記の（１１）式で与えられる。
ｆ₀＋（ｋ／Ｎ）ｆ_s＜ｆ_s／２ …（１１）
【００４７】
なお、ＦＦＴでは、形式上、ｆ_s／２からｆ_sまでの周波数に対するスペクトルデータも得られるが、これらデータは、下記の（１２）式に示すように、０からｆ_s／２に対する負の周波数のスペクトルを求めることになり意味を持たない。
ｅｘｐ｛ｉ２π（ｆ_s−ｆ_ｋ）ｎτ｝＝ｅｘｐ｛ｉ２π（−ｆ_ｋ）ｎτ｝…（１２）
【００４８】
このようにＦＦＴ等の離散フーリエ級数展開による周波数解析においては、隣接するモード周波数の弁別能力が全観測時間Ｔにより制約されるとともに、スペクトルデータの周波数の上限がサンプリング周波数ｆ_sにより決まる。
【００４９】
これに対し、打検では、周波数弁別能力とは異なり、モード周波数を同定する能力が要求される。そして、周波数同定能力については、上述した相反則は適用されない。このため、観測時間の短い波形データからであっても、周波数弁別能力以上の精度で、ＦＦＴ等の離散フーリエ級数展開を用いてモード周波数を同定することは可能である。
【００５０】
また、上記の（１０）式に示したフーリエスペクトルは、下記の（１３）式に変形することができる。
【００５１】
【数５】

【００５２】
上記の（１３）式中の（ｘ（ｎτ）ｅ^{-i2 πｆ 0 ・ｎτ}）は、ｋに無関係である。したがって、この括弧中はＦＦＴ等の離散フーリエ級数展開を行う前に計算することができる。これは、この括弧中の複素位相項を窓関数として扱えることを意味する。この複素位相項を実数窓関数と掛け合わせて生成した窓関数を、複素窓関数と称する。複素窓関数の実数窓関数部には、解析対象とする振動の種類に応じて、ハニング窓関数やハミング窓関数などの、種々の実数関数を用いることができる。
【００５３】
そして、複素窓関数の複素位相項には任意の値を与えることができる。すなわち、複素位相項中のｆ₀の値を任意に設定することができる。その結果、通常の窓関数が波形データ列の端部における位相不整合を解消する目的で用いられるのに対し、複素窓関数は、スペクトルデータの表示周波数の原点を周波数軸上で移動させるために用いることができる。これにより、後述の各実施形態に示すように、高精度でモード周波数を決定することが可能となる。
【００５４】
つぎに、ＦＦＴの演算処理速度について説明する。ＦＦＴの演算処理速度は、最も処理時間のかかる複素乗算の回数で表される。複素乗算は、前処理において波形データ列に窓関数をかける処理、及び、ＦＦＴアルゴリズムの中核であるバタフライ演算で実行される。なお、複素乗算の回数のカウントにあたっては、１回の（実数）×（実数）乗算は複素数どうしの乗算の１／４回としてカウントされる。また、１回の（実数）×（複素数）乗算は複素数どうしの乗算の１／２回としてカウントされる。
【００５５】
先ず、通常のＦＦＴにおける複素乗算回数について説明する。上記の（１３）式において、ｅ^{-i2 π}＝Ｗ：回転子（Ｒｏｔａｔｏｒ）と置き換えると、上記（１３）式は下記の（１４）式で表される。ただし、下記の（１４）式中のｘ₀（ｎ）は、波形データｘ（ｎ）に実数窓関数Ｈ（ｎ）を乗じて生成された切出波形データを表す。
【００５６】
【数６】

【００５７】
ここで、Ｎを２^Pに限定し、ｎ、ｋを下記の（１５）及び（１６）式のように２進数表記する。
ｎ＝［ｎ_p-1，ｎ_p-2，…，ｎ_r，…，ｎ₁，ｎ₀］
＝２^P-1ｎ_P-1＋２^P-2ｎ_P-2＋…＋２^rｎ_r＋…＋２¹ｎ₁＋２⁰ｎ₀ …（１５）
ｋ＝［ｋ_p-1，ｋ_p-2，…，ｋ_s，…，ｋ₁，ｋ₀］
＝２^P-1ｋ_P-1＋２^P-2ｋ_P-2＋…＋２^rｋ_r＋…＋２^sｋ_s＋２⁰ｋ₀ …（１６）
これにより、上記の（１４）式は、下記の（１７）式のように表される。
【００５８】
【数７】

【００５９】
そして、上記の（１７）式中の回転子Ｗの指数項を展開すると、下記の（１８）式のように表される。
【００６０】
【数８】

【００６１】
また、ｎ_r、ｋ_sの値は、０又は１であるので、上記の（１８）式は、下記の（１９）式で表される。
【００６２】
【数９】

【００６３】
上記の（１９）式より、Ｘ（ｋ）がｎ_rについて帰納的に演算できることが分かる。さらに、下記の（２０）式に示す、｛ｘ（ｎ）｝のビット反転（Bit Reversion）数列を導入すると、上記の（１９）式は、下記の（２１）式のように表される。
ｙ₀（ｍ）＝ｙ₀（［ｍ_p-1，…，ｍ_r，…，ｍ₀］）
＝ｘ（［ｎ₀，…，ｎ_p-r-1，…，ｎ_p-1］）
…（２０）
【００６４】
【数１０】

【００６５】
上記（２１）式をバタフライ演算の漸化式に書き表すと、下記の（２２）及び（２３）式のように表される。
【００６６】
【数１１】

【００６７】
データ点数２^P点に対して、上記のバタフライ演算の各段では、ｒ次に対応するｍ_p-r＝１の項に回転子を乗じる複素乗算が行われる。その回数は２^P／２回である。これが（Ｐ−１）段繰り返されるので、バタフライ演算全体としては、２^P-1（Ｐ−１）回の複素乗算が行われる。
【００６８】
また、通常のＦＦＴでは、窓関数処理には２^P回の（実数）×（実数）乗算が行われるので、窓関数処理における複素乗算回数は、２^P／４回となる。したがって、通常のＦＦＴにおける複素乗算回数は、これらを合わせて、下記の表１に示すように、２^P-1（Ｐ−１／２）回となる。
【００６９】
例えば、サンプリング周波数ｆ_s＝２０ｋＨｚの場合において、データ点数１２８点（Ｐ＝７）のときの０Ｈｚから１００００Ｈｚの範囲で複素乗算回数は、下記の表２に示すように、４１６回である。
【００７０】
また、０点付加法により、２^P点のデータ点数を２^P+Q点に拡張したＦＦＴの場合、バタフライ演算部の複素乗算回数は、２^P+Q-1（Ｐ＋Ｑ−１）である。また、窓関数処理における複素乗算回数は、通常のＦＦＴと同じ２^P／４回となる。したがって、０点付加法によるＦＦＴにおける複素乗算回数は、これらを合わせて、下記の表１に示すように、２^P+Q-1（Ｐ＋Ｑ−１＋（１／２^Q+1））回となる。例えば、サンプリング周波数ｆ_s＝２０ｋＨｚの場合において、データ点数１２８点（Ｐ＝７）、拡張後データ点数２０４８点（Ｑ＝４）のときの０Ｈｚから１００００Ｈｚの範囲で複素乗算回数は、下記の表２に示すように、１０２７２回である。
【００７１】
【表１】

【００７２】
【表２】

【００７３】
このように、０点付加法では、分解能は高くなるが、複素乗算回数が大幅に増加するため演算負荷が重たくなる。そこで、短時間の波形データであっても高い精度でモード周波数を同定するとともに、演算負荷の増加を抑制できる周波数解析技術の例について以下に説明する。
【００７４】
［第一実施形態］
まず、図１を参照して、第一実施形態の打検装置の構成について説明する。なお、本実施形態の打検装置は、本発明の周波数解析装置の一例でもある。本実施形態では、ベルトコンベア１に載せられて搬送中の缶２を被検査体とし、打検プローブ４により打検を行う。
【００７５】
打検プローブ４は、図２に示すように、缶２の缶底に衝撃を与えるエキサイタコイル４１と、缶２の振動音を集音して電気的な反響信号を生成するマイク４２とにより構成されている。エキサイタコイル４１に駆動部１２からインパルスが印加されると、磁力により缶２の缶底面が強制的に変位させられる。その結果、缶２は振動音を発する。
【００７６】
なお、図２では、側壁部と底部とが一体成形された缶胴と、上蓋都とから構成される２ピース缶の断面図を缶２として示したが、本発明では、側壁部のみからなる缶胴、底蓋、及び上蓋とから構成される３ピース缶の場合に対しても、同様に打検を行うことができる。
【００７７】
駆動部１２は、存在検知器３により缶２が検知されてから一定時間経過後、検知された缶２が打検プローブ４の真下を通過するタイミングで、エキサイタコイル４１にインパルスを印加する。存在検知器３は、ベルトコンベア１を挟んで設置された光源３１と受光素子３２とから構成されている。
【００７８】
マイク２により生成された反響信号は、フィルタ５を介してコンピュータ１０に送られる。フィルタ５は、反響信号に含まれる固有振動予測領域を大きく外れる高調波を含む高周波領域や低周波ノイズ領域を除去する。
【００７９】
コンピュータ１０は、メモリ（図示せず）に格納されているプログラムに従って周波数解析装置として動作する。ここで、図３に、コンピュータ１０の機能ブロックを示す。図３に示すように、コンピュータ１０は、デジタル変換部１１、前処理部１２、フーリエ変換部１３、同定部１４及び判定部１５により構成されている。
【００８０】
デジタル変換部１１は、反響信号をアナログ／デジタル変換し、一定のサンプリング時間間隔でサンプリングされた波形データ列｛ｘ（ｎ）｝を生成する。波形データ列は、いったんコンピュータ１０内のメモリ（図示せず）に格納される。
【００８１】
ここで、図４に、波形データの一例のグラフを示す。グラフ中のプロットは、波形データ列を構成する個々の波形データである。そして、グラフ中、各波形データのプロットを曲線Ｉで結んで示す。この波形データは、全観測時間６．４ｍｓの反響信号をサンプリング周波数２０ｋＨｚでサンプリングした１２８点のデータからなる。
【００８２】
また、前処理部１２は、メモリに格納されていた波形データ列を読み出し、波形データ列に窓関数を乗じて切出波形データ列を生成する。ここでは、窓関数として、実数窓関数Ｈ（ｎ）に複素位相項ｅｘｐ｛−ｉ２π（ｎ・ｍ）／（Ｎ・Ｍ）｝を乗じた、下記の（１）式に示す複素窓関数Ｃ（ｎ，ｍ）を使用する。
Ｃ（ｎ，ｍ）＝Ｈ（ｎ）・ｅｘｐ（−ｉ２π（ｎ・ｍ）／（Ｎ・Ｍ））…（１）
ただし、ｎ＝０，１，…，Ｎ−１を満たす。Ｎは波形データのサンプリング数を表し、１以上の整数を表す。また、ｍ＝０，１，…，Ｍ−１を満たす。Ｍは、スペクトルデータ列の周波数間隔の分割数を表し、１以上の整数である。また、Ｈ（ｎ）は、例えばハニング窓関数等の実数関数の窓関数を表す。
【００８３】
波形データ列｛ｘ（ｎ）｝に、この複素窓関数Ｃ（ｎ，ｍ）を乗じると、切出波形データ列｛Ｃ（ｎ，ｍ）・ｘ（ｎ）｝が得られる。生成された切出波形データ列は、いったんコンピュータ１０内のメモリ（図示せず）に格納される。そして、複素位相項中のｍ値を変化させた複素窓関数を波形データ列｛ｘ（ｎ）｝に乗じることにより、位相をずらした切出波形データ列｛Ｃ（ｎ，ｍ）・ｘ（ｎ）｝を生成することができる。
【００８４】
ここで、図５の（Ａ）に、Ｍ＝２^Q＝１６、ｍ＝０の場合の複素窓関数Ｃ（ｎ，０）を曲線IIで示し、切出波形データ列｛Ｃ（ｎ，０）・ｘ（ｎ）｝のグラフを曲線Iaに示す。ｍ＝０の場合、複素窓関数Ｃ（ｎ，０）は、複素位相項の値が１となるので、窓関数は実数関数Ｈ（ｎ）と同じとなる。
【００８５】
また、図５の（Ｂ）に、ｍ＝５の場合の複素窓関数Ｃ（ｎ，５）の実部関数を曲線IIIで示し、虚部関数を破線IVで示す。ただし、虚部関数は、見やすくするため、符号を反転して表している。また、切出波形データ列｛Ｃ（ｎ，５）・ｘ（ｎ）｝の実部のグラフを曲線Ibに示す。なお、図５の（Ｂ）には図示していないが、虚部関数で処理した虚部波形データも、実部波形データとともに、ＦＦＴ演算に用いられる。
【００８６】
また、フーリエ変換部１３は、メモリから切出波形データ列を読み出し、切出波形データ列｛Ｃ（ｎ，ｍ）・ｘ（ｎ）｝からＦＦＴによりスペクトルデータ列｛Ｘ（ｋ）｝を生成する。その際、フーリエ変換手段１３は、複素窓関数の位相項の値を変化させて位相をずらした切出波形データ列から、周波数原点ｆ₀をずらしたスペクトルデータ列を追加生成する。生成されたスペクトルデータ列は、コンピュータ１０内のメモリ（図示せず）に格納される。
【００８７】
ここで、図６の（Ａ）のグラフに、図５の（Ａ）に示したｍ＝０の場合の切出波形データ列に対してＦＦＴを行って求めたスペクトルデータを示す。なお、ここでは、ピーク周波数の近傍のみを示す。また、図６の（Ｂ）のグラフに、図５の（Ｂ）に示したｍ＝５の場合の切出波形データ列に対してＦＦＴを行って求めたスペクトルデータを示す。図６の（Ａ）及び（Ｂ）のグラフから、ｍ＝５の場合のスペクトルデータの周波数はｍ＝０の場合のスペクトルデータの周波数よりも高周波数側にシフトしていることが分かる。
【００８８】
そして、通常のＦＦＴにより求めたＮ点のスペクトルデータ列の周波数間隔ｆ_s／ＮをＭ等分したスペクトルデータ列を求めるためには、周波数原点ｆ₀を（１／Ｍ）（ｆ_s／Ｎ）ずつずらせばよい。その場合の周波数原点は、下記の（２４）式で表される。
ｆ₀（ｍ）＝（ｍ／Ｍ）（ｆ_s／Ｎ） …（２４）
ただし、ｍ＝０，１，２，…，Ｍ−１
【００８９】
さらに、切出波形データ列｛Ｃ（ｎ，ｍ）・ｘ（ｎ）｝に対して通常のＦＦＴ解析を行って得たスペクトルデータ列｛Ｘ（ｋ，ｍ）｝の周波数は、下記の（３）式で与えられる。
ｆ_k,m＝（ｆ_s／Ｎ）（ｋ＋ｍ／Ｍ）…（３）
ただし、ｍ＝０，１，…，Ｍ−１
【００９０】
上記（３）式において、ｍの値を変化させることにより、スペクトルデータ列の各周波数ｆ_k,mを（ｆ_s／Ｎ）・（１／Ｍ）ずつずらすことができる。「Ｍ」の値は無制限に大きくすることができるので、（１／Ｍ）の値も無制限に小さくすることができる。
【００９１】
ここで、図７のグラフに、通常のＦＦＴにより求められるスペクトルデータの周波数間隔をＭ＝１６分割したスペクトルデータを示す。グラフ中の各ドットはスペクトルデータを示す。また、太い縦線を伴うドットは、図６の（Ａ）に示したｍ＝０の場合のスペクトルデータである。また、細い縦線を伴うドットは、図６の（Ｂ）に示したｍ＝５の場合のスペクトルデータである。
【００９２】
また、同定部１４は、メモリに格納されていたスペクトルデータを読み出し、最大値となるスペクトルデータの周波数（ピーク周波数）を振動音のモード周波数として同定する。このように、複素窓関数を用いることにより、周波数間隔を分割したスペクトルデータを生成できるので、単一モード又は他と十分に独立したモードの振動波形のモード周波数を精度よく求めることができる。
【００９３】
また、判定部１５は、同定されたモード周波数に基づいて被検査体の良否を判定する。判定は、被検査体の缶２の種類（材質、形式、大きさ、内容物等）に応じて適正周波数の範囲を示すテーブルを設定しておき、モード周波数が適正周波数の範囲に含まれるか否かで行う。モード周波数が適正基準の範囲外である場合、判定部１５により異常缶と判定される。異常缶と判定された缶２は、リジェクタ７によってベルトコンベア１上から排除される。
【００９４】
なお、モード周波数は、缶内圧力を反映したものであり、缶内圧力は温度の影響を受ける。このため、温度センサ６の検出温度により、同定されたモード周波数の補正を行う。
【００９５】
また、本実施形態の打検装置においては、コンピュータ１０に、設定入力部８とディスプレイ９が接続されている。設定入力部８からは、例えば、波形データのサンプリング間隔や、被検査体である缶２の種類を入力することができる。また、ディスプレイ９には、例えば、設定入力画面や、周波数分析されたスペクトルや適正周波数範囲を画像表示することができる。
【００９６】
つぎに、図８を参照して、第一実施形態におけるコンピュータ１０による周波数解析方法（周波数分割法）について説明する。周波数分割法においては、周波数原点ｆ₀（ｍ）をｆ_s／（ＮＭ）ずつずらしたスペクトルデータ列を生成することにより、通常のＦＦＴで得られるスペクトルデータ列の表示周波数間隔ｆ_s／ＮをＭ等分する。そのために、先ず、引数ｍ＝０と設定する（Ｓ１０１）。
【００９７】
続いて、引数ｍ＝０の場合の前処理を行う（Ｓ１０２）。すなわち、上記の（１）式に示した複素窓関数Ｃ（ｎ，ｍ）を波形データ列｛ｘ（ｎ）｝に乗じ、切出波形データ列｛Ｃ（ｎ，ｍ）・ｘ（ｎ）｝を生成する。
【００９８】
引数ｍ＝０の場合は、上記の（１）式中の複素位相項の値が１となるので、実質的に通常の実数の窓関数Ｈ（ｎ）を乗じることになる。続いて、切出波形データ列｛Ｃ（ｎ，０）・ｘ（ｎ）｝に対してＦＦＴを実行し、スペクトルデータ列を生成する（Ｓ１０３）。
【００９９】
ここで、図９の（Ａ）に、スペクトルデータ列の一例をグラフで示す。グラフ中、１回目のＦＦＴで生成された周波数ｆ₀，ｆ₁，ｆ₂，ｆ₃，…のスペクトルデータを黒丸印で示す。これらの周波数間隔はｆ_s／Ｎである。なお、スペクトルデータ列の各周波数のうち、周波数原点ｆ₀の周波数は、必ずしも０Ｈｚに限定されない。また、図９においては、周波数原点ｆ₀を縦軸からずらして示す。
【０１００】
つぎに、引数ｍの値をインクリメントし、ｍ＝１とする（Ｓ１０４）。さらに、インクリメント後の引数ｍの値が、分割数Ｍに達するまで、上述したＳ１０２〜Ｓ１０４の処理を繰り返す（Ｓ１０５）。
【０１０１】
その結果、引数ｍの値をインクリメントするごとに、複素窓関数の複素位相項の値が一定間隔ｅｘｐ（−ｉ２πｎ／（ＮＭ））で変化する。そして、前処理では、複素窓関数Ｃ（ｎ，ｍ）を波形データ列｛ｘ（ｎ）｝にそれぞれ乗じることにより、位相をずらした複数の切出波形データ列｛Ｃ（ｎ，ｍ）・ｘ（ｎ）｝が生成される（Ｓ１０２）。
【０１０２】
また、ＦＦＴでは、切出波形データ列から、周波数原点ｆ₀（ｍ）を一定周波数ｆ_s／（ＮＭ）ずつずらした複数のスペクトルデータ列｛Ｘ（ｋ，ｍ）｝が追加生成される（Ｓ１０３）。
【０１０３】
例えば、Ｍ＝４として、表示周波数間隔ｆ_s／Ｎを４分割した結果を図９の（Ｂ）に示す。２回目のＦＦＴで、周波数ｆ₀＋（１／４）（ｆ_s／Ｎ），ｆ₁＋（１／４）（ｆ_s／Ｎ），ｆ₂＋（１／４）（ｆ_s／Ｎ），…のスペクトルデータが生成される。グラフ中、２回目のＦＦＴにより追加生成されたスペクトルデータを白丸印で示す。
【０１０４】
また、３回目のＦＦＴで、周波数ｆ₀＋（２／４）（ｆ_s／Ｎ），ｆ₁＋（２／４）（ｆ_s／Ｎ），ｆ₂＋（２／４）（ｆ_s／Ｎ），…のスペクトルデータが生成される。グラフ中、３回目のＦＦＴにより追加生成されたスペクトルデータを白三角印で示す。
【０１０５】
また、４回目のＦＦＴで、周波数ｆ₀＋（３／４）（ｆ_s／Ｎ），ｆ₁＋（３／４）（ｆ_s／Ｎ），ｆ₂＋（３／４）（ｆ_s／Ｎ），…のスペクトルデータが生成される。グラフ中、４回目のＦＦＴにより追加生成されたスペクトルデータを白四角印で示す。
【０１０６】
そして、インクリメント後の引数ｍの値がＭに達するまで処理を繰り返すことにより、通常のＦＦＴで得られるスペクトルデータ列の表示周波数間隔ｆ_s／ＮをＭ等分した表示周波数間隔ｆ_s／（ＮＭ）でスペクトルデータが生成される。図９の（Ｂ）に示す例では、表示周波数間隔ｆ_s／（４Ｎ）でスペクトルデータが生成されている。
【０１０７】
つぎに、このようにして狭い表示周波数間隔で生成されたスペクトルデータのうちから、スペクトルが最大値を示すピーク周波数が、打検の振動のモード周波数として同定される（Ｓ１０６）。図９の（Ｂ）に示す例では、周波数ｆ₁＋（３／４）（ｆ_s／Ｎ）のスペクトルが最大値となり、この周波数がモード周波数として同定される。
【０１０８】
ところで、表示周波数間隔をＭ分割する方法は、Ｎ点ＦＦＴをＭ回繰り返すことに相当する。このため、通常のＮ点ＦＦＴよりも計算回数は増加する。しかし、本実施形態の方法（周波数分割法）は、０点付加法により計算した場合に比べると、より少ない計算回数で同一の周波数精度を得ることができる。
【０１０９】
波形データ点数Ｎ＝２^P点、分割数Ｍ＝２^Qに対して、周波数分割法におけるバタフライ演算部の複素乗算回数は、２^P-1・（Ｐ−１）・２^Q回である。また、窓関数処理における複素乗算回数は、２^P-1（２^Q−１）＋２^P／４回となる。したがって、周波数分割法における複素乗算回数は、これらを合わせて、上記の表１に示すように、２^P+Q-1（Ｐ−（１／２^Q+1））回となる。
【０１１０】
例えば、サンプリング周波数ｆ_s＝２０ｋＨｚの場合において、データ点数１２８点（Ｐ＝７）、Ｑ＝４のときの０Ｈｚから１００００Ｈｚの範囲で複素乗算回数は、上記の表２に示すように、７１３６回であり、０点付加法の１０２７２回よりも少ない。したがって、第一実施形態の周波数分割法によれば、演算負荷の増加を抑制しつつ、短時間の波形データであっても高い精度でモード周波数を同定することができることが分かる。
【０１１１】
［第二実施形態］
つぎに、図１０を参照して、第二実施形態について説明する。第二実施形態においても、上述した第一実施形態の打検装置と同一の構成で打検を行う。ただし、第二実施形態においては、第一実施形態よりも演算負荷を軽減するため、フォーカス法として、以下のようにして周波数解析を行う。
【０１１２】
フォーカス法においては、周波数原点ｆ₀（ｍ）を周波数軸上でｆ_s／（ＮＭ）ずつずらしたスペクトルデータ列を追加生成するにあたり、通常のＦＦＴで得られた最大値を有するスペクトルデータの周波数の近傍に含まれるスペクトルデータのみを選択的に求める。
【０１１３】
そこで、まず、周波数原点をずらしていない基礎スペクトルデータ列のうち、最大値となるスペクトルデータの周波数を中央周波数として求める。、そのために、先ず、引数ｍ＝０と設定する（Ｓ２０１）。続いて、引数ｍ＝０の場合の前処理を行う（Ｓ２０２）。すなわち、上述の第一実施形態において（１）式に示した複素窓関数Ｃ（ｎ，ｍ）を波形データ列｛ｘ（ｎ）｝に乗じ、切出波形データ列｛Ｃ（ｎ，ｍ）・ｘ（ｎ）｝を生成する。引数ｍ＝０の場合、上記（１）式中の複素三角関数の位相項の値が０となるので、実質的に通常の実数の窓関数Ｈ（ｎ）を乗じることになる。
【０１１４】
続いて、切出波形データ列｛Ｃ（ｎ，０）・ｘ（ｎ）｝に対してＦＦＴを実行する。最初のＦＦＴは、通常のＦＦＴと同じとなる。このＦＦＴにより、周波数原点ｆ₀をずらしていない基礎スペクトルデータ列が生成される（Ｓ２０３）。ここで、図１１の（Ａ）に、基礎スペクトルデータ列の一例をグラフで示す。グラフ中、最初のＦＦＴで生成された周波数ｆ₀，ｆ₁，ｆ₂，ｆ₃，…のスペクトルデータを黒丸印で示す。これらの周波数間隔はｆ_s／Ｎである。なお、スペクトルデータ列の各周波数のうち、周波数原点ｆ₀の周波数は、必ずしも０Ｈｚに限定されない。また、図１１においては、周波数原点ｆ₀を縦軸からずらして示す。
【０１１５】
そして、基礎スペクトルデータ列のうち、最大値となるスペクトルデータの周波数と、それに対応するｋの値を求め、その周波数を中央周波数ｆ_cとする。ｋ値は、基礎スペクトルデータ列に含まれるスペクトルデータの、周波数原点のスペクトルデータから数えた番号を表し、０以上の整数である。なお、周波数原点のスペクトルデータはｋ＝０に相当する。また、例えば、図１１の（Ａ）に示したグラフにおいては、周波数ｆ₂のスペクトルが最大値であるので、周波数ｆ₂が中央周波数となり、ｋ＝２となる。
【０１１６】
つぎに、中央周波数ｆ_cを基礎スペクトルデータ列の周波数間隔（ｆ_s／Ｎ）の半分ずつ前後にずらした候補区間［ｆ_c−ｆ_s／２Ｎ，ｆ_c＋ｆ_s／２Ｎ］内の周波数へ基礎スペクトルデータ列の周波数原点を移動させる。そのために、下記の（２ａ）式の複素窓関数を、波形データ列｛ｘ（ｎ）｝に乗じて切出波形データ列を生成する。なお、下記の（２ａ）式は、（１）式に示した複素窓関数Ｃ（ｎ，ｍ）に、第二の複素位相項ｅｘｐ（−ｉ２πｎｋ／Ｎ）を乗じたものである。
【０１１７】
Ｃ_k（ｎ，ｍ）＝Ｃ（ｎ，ｍ）・ｅｘｐ（−ｉ２πｎｋ／Ｎ）＝Ｈ（ｎ）・ｅｘｐ（−ｉ２π（ｎ・ｍ）／（ＮＭ））・ｅｘｐ（−ｉ２πｎｋ／Ｎ）＝Ｈ（ｎ）・ｅｘｐ（−ｉ２πｎ（ｋＭ＋ｍ）／（ＮＭ））…（２ａ）
【０１１８】
ただし、候補区間［ｆ_c−ｆ_s／２Ｎ，ｆ_c＋ｆ_s／２Ｎ］内へ周波数原点を移動させるため、ｍ≦Ｍ／２の場合は、上記（２ａ）式を用い、ｍ≧Ｍ／２の場合は、下記の（２ｂ）式を用いる。
Ｃ_k-1（ｎ，ｍ）＝Ｃ（ｎ，ｍ）・ｅｘｐ（−ｉ２πｎ（ｋ−１）／Ｎ）＝Ｈ（ｎ）・（−ｉ２πｎ（（ｋ−１）Ｍ＋ｍ）／（ＮＭ））…（２ｂ）
なお、ｍ＝Ｍ／２のときは、上記の（２ａ）式と（２ｂ）式の両方でそれぞれ用いてスペクトルデータを求める。
【０１１９】
そして、まず、上記の（２ａ）式において、引数ｍの値をインクリメントし、ｍ＝１とする（Ｓ２０５）。続いて、引数ｍ＝１の場合の前処理を行う（Ｓ２０６）。すなわち、上記（２ａ）式に示した複素窓関数Ｃ_k（ｎ，１）を波形データ列｛ｘ（ｎ）｝に乗じ、切出波形データ列｛Ｃ_k（ｎ，１）・ｘ（ｎ）｝を生成する。引数ｍ＝１の場合、上記（２ａ）式中の複素三角関数の位相項の値は、ｅｘｐ（−ｉ２πｎ（ｋＭ＋１）／（ＮＭ））となる。
【０１２０】
続いて、切出波形データ列｛Ｃ_k（ｎ，１）・ｘ（ｎ）｝と、中央周波数に対応するｋを用いて、上記（１４）式によりフーリエスペクトルを計算する。この場合、離散フーリエ級数展開により、周波数原点ｆ₀をｆ_c＋（１／Ｍ）（ｆ_s／Ｎ）へ移動させたスペクトルデータ列｛Ｘ（ｋ，１）｝を生成することが可能である。しかし、第二実施形態においては、このスペクトルデータ列｛Ｘ（ｎ，１）｝のうち、周波数原点に対応するスペクトルのみを選択的に追加生成する（Ｓ２０７）。なぜならば、スペクトルデータ列｛Ｘ（ｎ，１）｝のうち、周波数原点のみが、中央周波数ｆ_cを基礎スペクトルデータ列の周波数間隔（ｆ_s／Ｎ）の半分（１／２）（ｆ_s／Ｎ）ずつ前後にずらした候補区間［ｆ_c−ｆ_s／２Ｎ，ｆ_c＋ｆ_s／２Ｎ］に含まれるからである。
【０１２１】
例えば、図１１の（Ａ）に示した例では、中央周波数ｆ₂の前後、［ｆ₂−（１／２）（ｆ_s／Ｎ），ｆ₂＋（１／２）（ｆ_s／Ｎ）］の候補区間内に含まれる周波数原点のスペクトルデータのみが生成される。Ｍ＝４の場合、図１１の（Ｂ）に示すように、新たな周波数原点ｆ₂＋（１／４）（ｆ_s／Ｎ）に対するスペクトルデータが生成される。グラフ中、引数ｍ＝１に対応する離散フーリエ級数展開により追加生成されたスペクトルデータを白丸印で示す。
【０１２２】
つぎに、引数ｍの値をインクリメントし、ｍ＝２とする（Ｓ２０８）。さらに、インクリメント後の引数ｍの値が、分割数Ｍに達するまで、上述したＳ２０６〜Ｓ２０７の処理を繰り返す（Ｓ２０９）。ただし、引数ｍの値がＭ／２を越えたときは、前記のｋの値を、中央周波数に対応する値から１減じた（ｋ−１）として、上述したＳ２０６〜Ｓ２０７のフォーカス解析処理を繰り返す。また、引数ｍの値がＭ／２に等しいときは、前記のｋの値が中央周波数に対応する値と、中央周波数に対応する値から１減じた値とについて、各々上述したＳ２０６〜Ｓ２０７のフォーカス解析処理を繰り返す。
【０１２３】
その結果、周波数原点ｆ₀（ｍ）を一定周波数ｆ_s／（ＮＭ）ずつずらした各スペクトルデータ列｛Ｘ（ｋ，ｍ）｝のうち、中央周波数ｆ_cを周波数間隔（ｆ_s／Ｎ）の半分（１／２）（ｆ_s／Ｎ）ずつ前後にずらした区間に含まれる周波数原点に対応するスペクトルデータだけが選択的に追加生成される。
【０１２４】
ここで、Ｍ＝４の場合のフォーカス解析結果を図１１の（Ｂ）に示す。２回目（ｍ＝２）のフォーカス解析処理では、引数ｍ＝２がＭ／２に等しい。このため、上記の（２ａ）式により、ｆ₂＋（１／２）（ｆ_s／Ｎ）に移動した周波数原点に対応するスペクトルデータを追加生成するとともに、上記の（２ｂ）式により、ｆ₁＋（１／２）（ｆ_s／Ｎ）に移動した周波数原点に対応するスペクトルデータを生成する。グラフ中２回目のフォーカス解析処理により追加生成されたスペクトルデータを白三角印で示す。なお、周波数ｆ₁＋（１／２）（ｆ_s／Ｎ）は、ｆ_２−（１／２）（ｆ_s／Ｎ）と同一である。
【０１２５】
また、３回目（ｍ＝３）のフォーカス解析処理では、ｍ≧Ｍ／２となる。このため、上記の（２ｂ）式により、ｆ₁＋（３／４）（ｆ_s／Ｎ）に移動した周波数原点に対応するスペクトルデータを生成する。グラフ中、３回目のフォーカス解析処理により追加生成されたスペクトルデータを白四角印で示す。
【０１２６】
このようにして、インクリメント後の引数ｍの値がＭ−１に達するまで処理を繰り返すことにより、中央周波数の前後、表示周波数間隔の半分の区間内でのみ、基礎スペクトルデータ列の表示周波数間隔ｆ_s／ＮをＭ等分した表示周波数間隔ｆ_s／（ＮＭ）でスペクトルデータが生成される。例えば、図１１の（Ｂ）に示す例では、区間［ｆ₂−（１／２）（ｆ_s／Ｎ），ｆ₂＋（１／２）（ｆ_s／Ｎ）］の範囲内でのみ、表示周波数間隔ｆ_s／（４Ｎ）でスペクトルデータが生成されている。
【０１２７】
つぎに、このようにして狭い表示周波数間隔で生成されたスペクトルデータのうちから、スペクトルが最大値を示すピーク周波数を、打検の振動のモード周波数として同定する（Ｓ２１０）。図１１の（Ｂ）に示す例では、区間［ｆ₂−（１／２）（ｆ_s／Ｎ），ｆ₂＋（１／２）（ｆ_s／Ｎ）］のうち、周波数ｆ₁＋（３／４）（ｆ_s／Ｎ）のスペクトルが最大値となり、この周波数がモード周波数として同定される。
【０１２８】
このように、フォーカス解析処理では、ｋの値として、中央周波数ｆ_cに対応する値（ｋ）と、この値から１を減じた値（ｋ−１）とのみが用いられる。このことを利用すると、前処理（Ｓ２０６）と、離散フーリエ級数展開（Ｓ２０７）の一部を同時に行うことができ、これによって複素乗算回数を大幅に軽減することができる。
【０１２９】
第二実施形態における上述した周波数解析方法をまとめれば、つぎのようになる。すなわち、引数ｍ≦（Ｍ／２）の場合、引数ｍに対する複素窓関数Ｃ（ｎ，ｍ）と、中央周波数に対応したｋに対する複素位相項ｅｘｐ（−ｉ２πｎｋ／Ｎ）とをかけた、新たな複素窓関数Ｃ_k（ｎ，ｍ）を生成する。この複素窓関数Ｃ_k（ｎ，ｍ）の複素位相項はｅｘｐ（−ｉ２πｎ（ｋＭ＋ｍ）／（ＮＭ））となる。この新たな複素窓関数を波形データ列｛ｘ（ｎ）｝に乗じて切出波形データ列｛Ｃ_k（ｎ，ｍ）・ｘ（ｎ）｝を生成する。この切出波形データ列から求められるフーリエスペクトルデータ列の周波数原点は、中央周波数ｆ_cとすると、ｆ_c＋（ｍ／Ｍ）（ｆ_s／Ｎ）となる。
【０１３０】
同様に、引数ｍ≧（Ｍ／２）の場合は、引数ｍに対する複素窓関数Ｃ（ｎ，ｍ）と、中央周波数に対応したｋから１を減じた値に対する複素位相項とをかけた、新たな複素窓関数Ｃ_k-1（ｎ，ｍ）を波形データ列｛ｘ（ｎ）｝に乗じて切出波形データ列｛Ｃ_k-1（ｎ，ｍ）・ｘ（ｎ）｝を生成する。この切出波形データ列から求められるフーリエスペクトルデータ列の周波数原点は、中央周波数ｆ_cとすると、ｆ_c−（１−（ｍ／Ｍ））（ｆ_s／Ｎ）となる。
【０１３１】
これらの周波数原点は中央周波数ｆ_cの前後、［ｆ_c−（１／２）（ｆ_s／Ｎ），ｆ_c＋（１／２）（ｆ_ｓ／Ｎ）］の範囲内にある。このため、生成した切出波形データ列からスペクトルデータを生成するにあたっては、各々のスペクトルデータ列のうち、周波数原点のスペクトルデータを計算することとなり、この場合は、複素乗算を行うことなく、単に波形データ列の各項を加算すればよいので、複素乗算回数を軽減することができる。
【０１３２】
ところで、本実施形態のフォーカス法によれば、通常のＮ点ＦＦＴよりも計算回数は増加する。しかし、フォーカス法は、０点付加法により計算した場合に比べると、より少ない計算回数で同一の周波数精度を得ることができる。
【０１３３】
すなわち、データ点数Ｎ＝２^P点、分割数Ｍ＝２^Qに対して、フォーカス法においては、最初に行う通常のＦＦＴにおける複素乗算回数が、２^P-1（Ｐ−１／２）回である。また、その後に行う離散フーリエ級数展開では、中央周波数付近の区間だけで、２^Ｑ個の２^P回複素乗算を行う。したがって、フォーカス法における複素乗算回数は、これらを合わせて、上記の表１に示したように、２^P-1（Ｐ−１／２＋２^{Ｑ +1}）回となる。
【０１３４】
例えば、Ｐ＝７、Ｑ＝４の場合、合計の複素乗算回数は、上記の表２に示すように、２４６４回となるしたがって、フォーカス法によれば、演算負荷の増加を一層抑制しつつ、短時間の波形データであっても高い精度でモード周波数を同定することができることが分かる。
【０１３５】
［第三実施形態］
つぎに、図１２を参照して、第三実施形態について説明する。第三実施形態においても、上述した第一実施形態の打検装置と同一の構成で打検を行う。ただし、第三実施形態においては、第一及び第二実施形態よりも演算負荷を軽減するため、分割比較法として、以下のようにして周波数解析を行う。
【０１３６】
分割比較法によりモード周波数を求めるにあたり、まず、上述の第一実施形態において（１）式に示した複素窓関数Ｃ（ｎ，ｍ）を波形データ列｛ｘ（ｎ）｝に乗じ、切出波形データ列｛Ｃ（ｎ，ｍ）・ｘ（ｎ）｝を生成する。その際、ｍ＝０とし、実質的に通常の実数の窓関数Ｈ（ｎ）を乗じる（Ｓ３０１）。
【０１３７】
続いて、切出波形データ列｛Ｃ（ｎ，０）・ｘ（ｎ）｝に対してＦＦＴを実行する。このＦＦＴは通常のＦＦＴと同じことになり、周波数原点ｆ₀をずらしていない基礎スペクトルデータ列｛Ｘ（ｋ，０）｝が生成される（Ｓ３０２）ここで、図１３の（Ａ）に、基礎スペクトルデータ列の一例のグラフで示す。グラフ中、上記のＦＦＴで生成された周波数ｆ₀，ｆ₁，ｆ₂，ｆ₃，…のスペクトルデータを黒丸印で示す。これら周波数間隔はｆ_s／Ｎである。なお、スペクトルデータ列の各周波数のうち、周波数原点ｆ₀の周波数は、必ずしも０Ｈｚに限定されない。また、図１３においては、周波数原点ｆ₀を縦軸からずらして示す。
【０１３８】
そして、周波数原点ｆ₀をずらしていない基礎スペクトルデータ列のうち、最大値となるスペクトルデータを求め、その周波数を基準周波数とする（Ｓ３０３）。例えば、図１３の（Ａ）に示した例においては、周波数ｆ₂のスペクトルが最大値であるので、周波数ｆ₂が基準周波数ｆ_cとなる。この場合、ｋ＝２となる。
【０１３９】
つぎに、基準周波数ｆ_cを基礎スペクトルデータ列の周波数間隔の半分である分割周波数間隔（１／２）（ｆ_ｓ／Ｎ）ずつ前後にずらした二つの周波数ｆ_c±（１／２）（ｆ_ｓ／Ｎ）へ周波数原点をそれぞれ移動させる。そのために、周波数原点ｆ₀を分割周波数間隔（１／２）（ｆ_ｓ／Ｎ）ずつ前後にずら複素窓関数Ｃ（ｎ，ｍ）＝Ｈ（ｎ）・ｅｘｐ（−ｉ２πｎ／２Ｎ）に、周波数間隔のｋ倍ずらす複素位相項ｅｘｐ（−ｉ２πｎｋ／Ｎ）を乗じた下記の（２ｃ）式に示す複素関数Ｃ_k（ｎ，ｍ）を使用して前処理を行う。
Ｃ_k（ｎ，ｍ）＝Ｈ（ｎ）・（−ｉ２πｎ（ｋ±（１／２））／Ｎ）…（２ｃ）
【０１４０】
そして、上記の（２ｃ）式に示した複素窓関数Ｃ_k（ｎ，ｍ）を、波形データ列｛ｘ（ｎ）｝に乗じ、一対の切出波形データ列｛Ｃ_k（ｎ，ｍ）・ｘ（ｎ）｝をそれぞれ生成する（Ｓ３０４）。なお、上記（２ｃ）式中では、ｍ／Ｍ＝±１／２に相当する。
【０１４１】
続いて、一対の切出波形データ列｛Ｃ_k（ｎ，ｍ）・ｘ（ｎ）｝に対して、上記（１４）式により、離散フーリエ級数展開を行う。この離散フーリエ級数展開によって、周波数原点ｆ₀を（１／２）（ｆ_ｓ／Ｎ）ずつ前後にシフトさせた二つのスペクトルデータ列を生成することが可能である。しかし、本実施形態では、二つスペクトルデータ列のうち、それぞれ周波数原点に対応する二つのスペクトルデータだけを選択的に追加生成する（Ｓ３０５）。
【０１４２】
例えば、図１３の（Ｂ）に示す例では、新たな周波数原点ｆ₂±（１／２）（ｆ_s／Ｎ）に対応する二つのスペクトルデータのみが生成される。グラフ中、この処理により追加生成された二つのスペクトルデータをそれぞれ白丸印で示す。
【０１４３】
さらに、追加生成された二つのスペクトルデータ及び基準周波数のスペクトルデータのうち、最大値をとるスペクトルデータを求め、その周波数を新たな基準周波数とする（Ｓ３０３）。例えば、図１３の（Ｂ）に示す例では、基準周波数ｆ₂、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）、及び、周波数ｆ₂＋（１／２）（ｆ_s／Ｎ）の三つのスペクトルデータのうち、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）のスペクトルデータの値が最も高い。このため、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）を新たな基準周波数とする。
【０１４４】
つぎに、最新の基準周波数を最新の分割後の周波数間隔のさらに半分である分割周波数間隔（１／４）（ｆ_ｓ／Ｎ）ずつ前後にずらした二つの周波数ｆ_c±（１／２）（ｆ_ｓ／Ｎ）へ、周波数原点をそれぞれ移動させる。そのために、図１３の（Ｂ）に示す例の場合、下記の（２ｄ）式に示す複素関数Ｃ_k（ｎ，ｍ）を使用して前処理を行う。
Ｃ_k（ｎ，ｍ）＝Ｈ（ｎ）・（−ｉ２πｎ（ｋ−（１／２）±（１／４））／Ｎ）…（２ｄ）
上記の（２ｄ）式では、ｍ／Ｍ＝−１／２±１／４に相当する。
【０１４５】
そして、上記の（２ｄ）式に示した複素窓関数Ｃ_k（ｎ，ｍ）を、波形データ列｛ｘ（ｎ）｝に乗じ、一対の切出波形データ列｛Ｃ_k（ｎ，ｍ）・ｘ（ｎ）｝をそれぞれ生成する（Ｓ３０４）。
【０１４６】
続いて、最新の一対の切出波形データ列｛Ｃ_k（ｎ，ｍ）・ｘ（ｎ）｝からそれぞれ生成されるスペクトルデータ列のうち、二つの周波数原点のスペクトルデータを選択的に追加生成する（Ｓ３０５）。ここでは、最新の基準周波数を（１／４）（ｆ_ｓ／Ｎ）ずつ前後にシフトさせた二つの周波数原点に対応するスペクトルデータだけを、上記（１４）式により、選択的に追加生成する。
【０１４７】
例えば、図１４の（Ａ）に示す例では、周波数原点ｆ_０が（（−１／２）±（１／４））（ｆ_s／Ｎ）シフトしたスペクトルデータ列のうち、最新の基準周波数ｆ₂−（１／２）（ｆ_s／Ｎ）を挟む、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）−（１／４）（ｆ_s／Ｎ）、及び、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）の二つの周波数原点に対するスペクトルデータのみが追加生成される。グラフ中、この処理により追加生成された二つのスペクトルデータをそれぞれ白三角印で示す。
【０１４８】
さらに、追加生成された二つのスペクトルデータ及び基準周波数のスペクトルデータのうち、最大値をとるスペクトルデータを求め、その周波数を新たな基準周波数とする（Ｓ３０３）。例えば、図１４の（Ａ）に示す例では、最新の基準周波数ｆ₂−（１／２）（ｆ_s／Ｎ）、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）−（１／４）（ｆ_s／Ｎ）、及び、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）の三つのスペクトルデータのうち、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）のスペクトルデータの値が最も高い。このため、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）を新たな基準周波数とする。
【０１４９】
つぎに、最新の基準周波数のスペクトルデータを含むスペクトルデータ列の周波数原点を、最新の分割後の周波数間隔（１／４）（ｆ_ｓ／Ｎ）のさらに半分にあたる（１／８）（ｆ_ｓ／Ｎ）ずつ前後にシフトさせる。そのため、図１４の（Ａ）に示す例では、下記の（２ｅ）式に示した複素窓関数Ｃ_k（ｎ，ｍ）を、波形データ列｛ｘ（ｎ）｝に乗じ、一対の切出波形データ列｛Ｃ_k（ｎ，ｍ）・ｘ（ｎ）｝を生成する前処理を行う（Ｓ３０４）。
Ｃ_k（ｎ，ｍ）＝Ｈ（ｎ）・（−ｉ２πｎ（ｋ−（１／２）＋（１／４）±（１／８））／Ｎ）…（２ｅ）
なお、上記（２ｅ）式では、ｍ／Ｍ＝−（１／２）＋（１／４）±（１／８）に相当する。
【０１５０】
続いて、最新の一対の切出波形データ列｛Ｃ_k（ｎ，ｍ）・ｘ（ｎ）｝からそれぞれ生成されるスペクトルデータ列のうち、二つの周波数原点のスペクトルデータを選択的に追加生成する（Ｓ３０５）。ここでは、最新の基準周波数を（１／８）（ｆ_ｓ／Ｎ）ずつ前後にシフトさせた二つの周波数原点に対応するスペクトルデータだけを、上記（１４）式により、選択的に追加生成する。
【０１５１】
例えば、図１４の（Ｂ）に示す例では、周波数原点を（−１／２＋１／４±１／８）（ｆ_s／Ｎ）シフトしたスペクトルデータ列のうち、最新の基準周波数ｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）を挟む、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）＋（１／８）（ｆ_s／Ｎ）、及び、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）−（１／８）（ｆ_s／Ｎ）の二つの周波数に対するスペクトルデータのみが生成される。グラフ中、この処理により追加生成された二つのスペクトルデータをそれぞれ白四角印で示す。
【０１５２】
そして、本実施形態では、分割後の周波数間隔が基準値以下となるまで、上述したＳ３０３〜Ｓ３０５の処理を繰り返す。そして、分割後の周波数間隔が基準値以下となった場合、その段階で、追加生成されている最新のスペクトルデータ及び最新の基準周波数のスペクトルデータのうちから、最大値をとるスペクトルデータの周波数を新たな基準周波数とし、最新の基準周波数をモード周波数として同定する。
【０１５３】
このようにして同定されたモード周波数は、例えば、下記の（２５）式で与えられる。ただし、下記の（２５）式中、ｆ_cは、最初のＦＦＴで生成した基礎スペクトルデータ列から得られた最初の基準周波数を表し、ｊ_mは、−１，０又は１を表す。また、Ｑは、基礎スペクトルデータ列の周波数間隔の分割数に相当する。
【０１５４】
【数１２】

【０１５５】
例えば、基準周波数がｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）の場合、図１４の（Ｂ）に示す例では、基準周波数ｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）＋（１／８）（ｆ_s／Ｎ）、及び、周波数ｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）−（１／８）（ｆ_s／Ｎ）のうちの周波数ｆ₂−（１／２）（ｆ_s／Ｎ）＋（１／４）（ｆ_s／Ｎ）、を再び新たな基準周波数とし、この周波数をモード周波数として同定する。
【０１５６】
このようにすれば、通常のＦＦＴによるスペクトルデータ列を計算した後、モード周波数近傍の周波数のスペクトルデータのみを選択的に計算してモード周波数を求めることができる。これにより、演算負荷のより一層の軽減を図ることができる。その結果、処理時間の一層の短縮を図ることができる。
【０１５７】
このように、分割比較法では、ｋの値は最初のＦＦＴで求めたスペクトルデータ列から得られる、最初の基準周波数ｆ_cに対応する値のみが用いられる。このことを利用すると、複素窓関数を用いて切出波形データ列を生成する前処理と、フーリエ変換処理を同時に行うことができ、これによって複素乗算回数をさらに軽減することができる。
【０１５８】
例えば、最初の基準周波数ｆ_cに対して、周波数原点ｆ₀を（１／２）（ｆ_s／Ｎ）ずつ前後にシフトさせた二つのスペクトルデータ列のうち、基準周波数を（１／２）（ｆ_s／Ｎ）ずつ前後にシフトさせた二つの周波数に対するスペクトルデータだけを、上記（１４）式により選択的に追加生成する処理においては、周波数原点ｆ₀を（１／２）（ｆ_ｓ／Ｎ）ずつ前後にシフトさせた複素窓関数Ｃ（ｎ，ｍ）と、基準周波数ｆ_cに対応したｋに対する複素位相項ｅｘｐ（−ｉ２πｎｋ／Ｎ）とをかけた、新たな複素窓関数Ｃ_k（ｎ，ｍ）を生成する。この複素窓関数Ｃ_k（ｎ，ｍ）の複素位相項はｅｘｐ（−ｉ２πｎ（ｋ±（１／２））／Ｎ）となる。この新たな複素窓関数を波形データ列｛ｘ（ｎ）｝に乗じて切出波形データ列｛Ｃ_k（ｎ，ｍ）・ｘ（ｎ）｝を生成する。この切出波形データ列から求められるフーリエスペクトルデータ列の周波数原点は、ｆ_c±（１／２）（ｆ_s／Ｎ）となる。
【０１５９】
これらの周波数原点は、選択的に追加生成するスペクトルデータの周波数と一致している。このため、生成した切出波形データ列からスペクトルデータを生成するにあたっては、周波数原点のスペクトルデータを計算することとなり、この場合は、複素乗算を行うことなく、単に波形データ列の各項を加算すればよいので、複素乗算回数を軽減することができる。
【０１６０】
ところで、本実施形態の分割比較法によれば、通常のＮ点ＦＦＴよりも計算回数は増加する。しかし、分割比較法は、０点付加法やフォーカス法により計算した場合に比べると、より少ない計算回数で同一の周波数精度を得ることができる。
【０１６１】
すなわち、データ点数Ｎ＝２^P点、分割数Ｍ＝２^Qに対して、分割比較法においては、最初の通常のＦＦＴにおける複素乗算回数が、２^P-1（Ｐ−１／２）回である。また、以降の離散フーリエ級数展開では、２^Q個の２Ｐ回複素乗算を行う。したがって、分割比較法における複素乗算回数は、これらを合わせて、上記の表１に示したように、２^P-1（Ｐ−１／２＋４Ｑ）回となる。例えば、Ｐ＝７、Ｑ＝４の場合、合計の複素乗算回数は、上記の表２に示すように、１４４０回となる。したがって、分割比較法によれば、演算負荷の増加を一層抑制しつつ、短時間の波形データであっても高い精度でモード周波数を同定することができることが分かる。
【０１６２】
［第四実施形態］
つぎに、第四実施形態について説明する。第四実施形態においても、上述した第一実施形態の打検装置と同一の構成で打検を行う。ただし、第四実施形態においては、第一、第二及び第三実施形態よりも演算負荷を軽減するため、以下のようにして周波数解析を行う。
【０１６３】
すなわち、サンプリング周波数が２０ｋＨｚの場合、通常のＦＦＴにおいては、０〜１００００Ｈｚの周波数範囲にわたって周波数解析を行っている。これに対して、打検による振動のモード周波数は、１０００Ｈｚ〜３０００Ｈｚの範囲内に含まれることが多い。そこで、本実施形態では、この周波数範囲の少なくとも一部分を含む特定の周波数範囲内のスペクトルデータのみを選択的に生成する。これにより、演算負荷を一層軽減し、処理速度の向上を図ることができる。
【０１６４】
なお、特定周波数範囲は、缶種ごとに設定するのが好ましい。例えば、ある種類の缶では１２００〜２０００Ｈｚ、他のある種類の缶では１５００〜３０００Ｈｚといった範囲を特定周波数範囲として設定するとよい。また、１０００〜３０００Ｈｚの全範囲を含むより広い範囲を特定周波数範囲として設定してもよい。例えば、１０００〜３５００Ｈｚを特定周波数範囲として設定してもよい。
【０１６５】
以下、ＦＦＴによる周波数解析範囲を２^S点に限定する場合の複素乗算回数について検討する。これは、上述したバタフライ演算の漸化式である（２２）及び（２３）式において、ｋ_p-1＝ｋ_p-2＝…＝ｋ_S＝０とした場合に対応する。そして、この条件を適用すると、下記の（２６）、（２７）及び（２８）式が得られる。
【０１６６】
【数１３】

【０１６７】
ここで、１＜ｒ≦Ｐ−Ｓの各段においては、ｋ_r＝０であることから、演算に関与するｙ_r（ｍ）の個数を順次１／２ずつ減らしていくことになり、Ｐ−Ｓ＜ｒ≦Ｐ−１の各段においては、２^S個のｙ_r（ｍ）を用いて演算を行っている。したがって、この演算における複素乗算回数は、下記の（２９）式により与えられる。ただし、（２９）式中、Ｒ＝Ｐ−Ｓである。
【０１６８】
【数１４】

【０１６９】
そして、周波数解析範囲を限定する手法は、上述した各実施形態と組み合わせることができる。上述した第一実施形態において説明した周波数分割法と併用した場合、複素乗算回数は、上記の表１に示すように、２^P+Q（（３／２）−（１／２^R））＋２^P+Q-R-1（Ｐ−Ｒ−１）回となる。例えば、Ｐ＝７、Ｑ＝４の場合、合計の複素乗算回数は、上記の表２に示すように、３５８４回となるしたがって、周波数分割法において周波数解析範囲を限定すれば、より一層演算負荷の増加を抑制することができる。
【０１７０】
また、上述した第二実施形態において説明したフォーカス法と併用した場合、複素乗算回数は、上記の表１に示すように、２^P（３／２−１／２^R＋２^Q）＋２^P-R-1（Ｐ−Ｒ−１）回となる。例えば、Ｐ＝７、Ｑ＝４の場合、合計の複素乗算回数は、上記の表２に示すように、２２７２回となる。したがって、フォーカス法において周波数解析範囲を限定すれば、より一層演算負荷の増加を抑制することができる。
【０１７１】
また、上述した第三実施形態において説明した分割比較法と併用した場合、複素乗算回数は、上記の表１に示すように、２^P（（３／２）−（１／２^R）＋２Ｑ）＋２^P-R-1（Ｐ−Ｒ−１）回となる。例えば、Ｐ＝７、Ｑ＝４の場合、合計の複素乗算回数は、上記の表２に示すように、１２４８回となるしたがって、分割比較法において周波数解析範囲を限定すれば、より一層演算負荷の増加を抑制することができる。
【０１７２】
上述した実施の形態においては、本発明を特定の条件で構成した例について説明したが、本発明は、種々の変更を行うことができる。例えば、上述した実施形態においては、打検による短時間の波形データからモード周波数を求める例について説明したが、本発明の周波数解析方法及び装置では、波形データは打検によるものに限定されない。また、波形データは短時間のものに限定されない。例えば、本発明をスペクトルアナライザに適用し、機械的振動等の解析にも利用することができる。そして、モード周波数に基づいて、例えば機械の故障を発見することも期待できる。
【０１７３】
なお、上記の実施形態における周波数解析処理は、プログラムに制御されたコンピュータにより実行することができる。このプログラムは、例えば、記録媒体により提供される。記録媒体としては、例えば、磁気ディスク、半導体メモリ等の磁気記録媒体、ＣＤ−ＲＯＭ等の光学式記録媒体、又は、光磁気記録媒体その他の任意の、コンピュータで読み取り可能なものを使用することができる。
【０１７４】
【発明の効果】
以上、詳細に説明したように、本発明によれば、複素窓関数の複素位相項の値を変化させ、スペクトルの周波数原点を周波数軸上で移動させる。その結果、任意の周波数を原点とするスペクトルデータ列を生成することができる。このため、スペクトルの表示周波数間隔は原理的に無限に細かくすることができる。すなわち、周波数同定能力は、原理的に制約を受けない。これにより、短時間の波形データであっても高い精度で打検のモード周波数を同定することができる。
【０１７５】
また、本発明によれば、通常のＦＦＴに加え、周波数原点をシフトさせた切出波形データ列に対してもＦＦＴ等の離散フーリエ級数展開を行うため、通常のＦＦＴに比べて演算量は増加する。しかしながら、上述した０点付加法と比較すると、より少ない演算量で同程度の精度でモード周波数を同定することができる。このため、０点付加法を用いた場合に比べて、本発明では演算負荷の増加を抑制しつつ、高精度でモード周波数を同定することができる。また、本発明は１モード振動のピーク周波数の検出に用いて特に好適である。また、本発明によれば、専用ＬＳＩチップを用いなくとも、汎用コンピュータによって、ＦＦＴを高精度かつ高速に実施することが可能となる。
【図面の簡単な説明】
【図１】第一実施形態の打検装置の構成図である。
【図２】打検プローブの模式図である。
【図３】打検装置を構成するコンピュータの機能ブロック図である。
【図４】第一実施形態における波形データ例を示すグラフである。
【図５】（Ａ）及び（Ｂ）は、窓関数及び切出波形データ例を示すグラフである。
【図６】（Ａ）及び（Ｂ）は、ＦＦＴにより得られたスペクトルデータ例を示すグラフである。
【図７】ＦＦＴにより得られたスペクトルデータ例を示すグラフである。
【図８】第一実施形態における周波数解析方法（周波数分割法）を説明するためのフローチャートである。
【図９】（Ａ）及び（Ｂ）は、第一実施形態における周波数解析方法を説明するためのグラフである。
【図１０】第二実施形態における周波数解析方法（フォーカス法）を説明するためのフローチャートである。
【図１１】（Ａ）及び（Ｂ）は、第二実施形態における周波数解析方法を説明するためのグラフである。
【図１２】第三実施形態における周波数解析方法（分割比較法）を説明するためのフローチャートである。
【図１３】（Ａ）及び（Ｂ）は、第三実施形態における周波数解析方法を説明するためのグラフである。
【図１４】（Ａ）及び（Ｂ）は、第三実施形態における周波数解析方法を説明するための、図１３に続くグラフである。
【符号の説明】
１ ベルトコンベア
２ 缶
３ 存在検知器
４ 打検プローブ
５ フィルタ
６ 温度センサ
７ リジェクタ
８ 設定入力部
９ ディスプレイ
１０ コンピュータ
１１ デジタル変換部
１２ 前処理部
１３ フーリエ変換部
１４ 同定部
１５ 判定部
４１ エキサイタコイル
４２ マイク

Claims (7)

1. 振動現象から一定のサンプリング時間間隔でサンプリングされた波形データ列に、窓関数を乗じて切出波形データ列を生成する前処理手段と、前記切出波形データ列から高速フーリエ変換法等の離散フーリエ級数展開によりスペクトルデータ列を生成するフーリエ変換手段と、前記スペクトルデータ列に含まれるスペクトルデータのうち、最大値となるスペクトルデータの周波数を前記振動音のモード周波数として同定する同定手段とを備える周波数解析装置であって、
   前記前処理手段は、前記窓関数として、実数窓関数に複素位相項を乗じた複素窓関数を使用し、前記複素窓関数の前記複素位相項の値を一定間隔で変化させた各複素窓関数を前記波形データ列にそれぞれ乗じることにより、互いに位相をずらした複数の切出波形データ列を生成し、
   前記フーリエ変換手段は、前記各切出波形データ列から、周波数原点を一定周波数ずつずらした複数のスペクトルデータ列を追加生成することを特徴とする周波数解析装置。
2. 振動現象から一定のサンプリング時間間隔でサンプリングされた波形データ列に、窓関数を乗じて切出波形データ列を生成する前処理手段と、前記切出波形データ列から高速フーリエ変換法等の離散フーリエ級数展開によりスペクトルデータ列を生成するフーリエ変換手段と、前記スペクトルデータ列に含まれるスペクトルデータのうち、最大値となるスペクトルデータの周波数を前記振動音のモード周波数として同定する同定手段とを備える周波数解析装置であって、
   前記同定手段は、周波数原点をずらしていない基礎スペクトルデータ列のうち、最大値となるスペクトルデータの周波数を中央周波数とし、
   前記前処理手段は、前記中央周波数を前記基礎スペクトルデータ列の周波数間隔の半分ずつ前後にずらした候補区間内の周波数へ前記基礎スペクトルデータ列の周波数原点を移動させる複素窓関数を、前記波形データ列に乗じて切出波形データ列を生成し、
   前記フーリエ変換手段は、前記切出波形データ列から生成されるスペクトルデータ列のうち、周波数原点のスペクトルデータを選択的に追加生成し、
   前記同定手段は、中央周波数のスペクトルデータ及び追加生成されたスペクトルデータのうちで最大値となるスペクトルデータの周波数をモード周波数として同定することを特徴とする周波数解析装置。
3. 振動現象から一定のサンプリング時間間隔でサンプリングされた波形データ列に、窓関数を乗じて切出波形データ列を生成する前処理手段と、前記切出波形データ列から高速フーリエ変換法等の離散フーリエ級数展開によりスペクトルデータ列を生成するフーリエ変換手段と、前記スペクトルデータ列に含まれるスペクトルデータのうち、最大値となるスペクトルデータの周波数を前記振動音のモード周波数として同定する同定手段とを備える周波数解析装置であって、
   （ａ）前記同定手段は、周波数原点をずらしていない基礎スペクトルデータ列のうち、最大値となるスペクトルデータの周波数を基準周波数とし、
   （ｂ）前記前処理手段は、前記基準周波数を前記基礎スペクトルデータ列の周波数間隔の半分である分割周波数間隔ずつ前後にずらした二つの周波数へ周波数原点をそれぞれ移動させる複素窓関数を、前記波形データ列にそれぞれ乗じて一対の切出波形データ列をそれぞれ生成し、
   （ｃ）前記フーリエ変換手段は、前記一対の切出波形データ列からそれぞれ生成されるスペクトルデータ列のうち、二つの周波数原点のスペクトルデータを選択的に追加生成し、
   （ｄ）前記同定手段は、追加生成された二つのスペクトルデータ及び前記基準周波数のスペクトルデータのうち、最大値をとるスペクトルデータの周波数を新たな基準周波数とし 、
   （ｅ）前記前処理手段は、最新の基準周波数を直前の分割周波数間隔のさらに半分の新たな分割周波数ずつ前後にずらした二つの周波数へ周波数原点をそれぞれ移動させる複素窓関数を、前記波形データ列にそれぞれ乗じて一対の新たな切出波形データ列を生成し、
   （ｆ）前記フーリエ変換手段は、最新の一対の切出波形データ列からそれぞれ生成されるスペクトルデータ列のうち、二つの周波数原点のスペクトルデータを選択的に追加生成し、
   （ｇ）前記同定手段は、追加生成された最新の二つのスペクトルデータ及び最新の基準周波数のスペクトルデータのうちから、最大値をとるスペクトルデータの周波数を新たな基準周波数とし、
   前記前処理手段、フーリエ変換手段及び同定手段は、前記（ｅ）〜（ｇ）の処理を繰り返した後、最新の基準周波数を前記モード周波数として同定することを特徴とする周波数解析装置。
4. 前記前処理手段は、前記複素窓関数として、実数窓関数Ｈ（ｎ）に第一の複素位相項を乗じた下記の（１）式に示すＣ（ｎ，ｍ）を含む複素関数を使用することを特徴とする請求項１〜３のいずれかに記載の周波数解析装置。
   Ｃ（ｎ，ｍ）＝Ｈ（ｎ）・ｅｘｐ（−ｉ２π（ｎ・ｍ）／（ＮＭ））…（１）
   ただし、ｎ＝０，１，…，Ｎ−１を満たす。Ｎは、前記波形データのサンプリング数を表し、１以上の整数である。また、ｍ＝０，１，…，Ｍ−１を満たす。Ｍは、前記スペクトルデータ列の周波数間隔の分割数を表し、１以上の整数である。
5. 前記前処理手段は、前記複素窓関数として、前記Ｃ（ｎ，ｍ）に、第二の複素位相項ｅｘｐ（−ｉ２πｎｋ／Ｎ）を乗じた下記の（２）式に示すＣ_k（ｎ，ｍ）を使用することを特徴とする請求項４記載の周波数解析装置。
   Ｃ_k（ｎ，ｍ）＝Ｃ（ｎ，ｍ）・ｅｘｐ（−ｉ２πｎｋ／Ｎ）…（２）
   ただし、ｋは、前記スペクトルデータ列に含まれるスペクトルデータの、周波数原点のスペクトルデータから数えた番号を表し、０以上の整数である。
6. 被検査体に衝撃を与える衝撃手段と、
   打撃を与えられた前記被検査体の振動音を集音して電気的な反響信号を生成する集音手段と、
   前記反響信号をアナログ／デジタル変換し、一定のサンプリング時間間隔でサンプリングされた波形データ列を生成するデジタル変換手段と、
   同定されたモード周波数に基づいて被検査体の良否を判定する判定手段とを備える打検装置であって、
   前記請求項１〜５のいずれかに記載の周波数解析装置を備え、当該周波数解析装置によって前記モード周波数を同定することを特徴とする打検装置。
7. 前記フーリエ変換手段は、前記スペクトルデータ列に含まれるスペクトルデータのうち、１０００〜３０００Ｈｚの周波数範囲の少なくとも一部分を含む特定の周波数範囲内のスペクトルデータのみを選択的に生成することを特徴とする請求項６に記載の打検装置。

Images