JP6677069B2 - 定ｑ変換の成分演算装置および定ｑ変換の成分演算方法 - Google Patents

Description

本発明は、振動や音響などの時系列データの周波数解析手法の一つである、定比周波数成分へのデータ変換に関し、特に定Ｑ変換の成分演算装置および方法に関する。

時系列データの周波数解析には不確定性原理があり、周波数解像度を上げようとすればより長時間のデータを解析する必要がある。フーリエ変換は全ての周波数帯域で同一の周波数解像度にて周波数成分に分割して解析するが、音楽解析や騒音解析、振動解析などではオクターブ解析など、周波数を等比の周波数成分に分割して解析することがある。

尚、オクターブ解析は、音響や振動の周波数を、１オクターブ（周波数比２：１の周波数）をＮ分割する等比の周波数帯域（１／Ｎオクターブバンド）に分割して解析する手法である。前記Ｎはビン数ともいい、１，３，６，１２，２４などがよく用いられる。

前記のように周波数を分割して解析を行う場合、高周波帯域の成分は低周波数帯域の成分と同じ周波数解像度は必要ない上に、長時間のデータを解析すると成分がその期間で平均的な値に均されるので、高周波帯域での早い変化を見逃す恐れもある。

このような状況に対応するために定Ｑ変換（Ｃｏｎｓｔａｎｔ Ｑ ｔｒａｎｓｆｏｒｍ）という手法が提案されており、例えば非特許文献１〜４などが存在する。

尚、非特許文献３、４はそれぞれ非特許文献１、２の手法の日本語による解説とＰｙｔｈｏｎ言語による検証例を報告しているＷｅｂサイトである。

この定Ｑ変換は、全ての周波数帯域で同じ期間のデータを解析するのではなく、全ての周波数帯域で同じ周期数になるように、周波数毎に参照するデータ数を変えて解析する。この方法では、高周波帯域ほど短い期間のデータで解析するため、低周波帯域で細かい周波数解像度を得ると同時に高周波帯域での早い変化を捉えることができる。

非特許文献１によると、定Ｑ変換は次で定式化される（表記は少し変更している）。

０以外の値は各ｋに対して、ｎ＝（Ｎ−Ｎ_k）／２，…，（Ｎ＋Ｎ_k）／２−１の範囲でのみ取り、この範囲の窓関数の形状として具体的にはハミング窓、ハニング窓、矩形窓など、多くの種類が提案されている
Ｎ_k：窓幅，Ｎ_k＝Ｑ（ｆ_s／ｆ_k）となるように取る
Ｎ≧Ｎ₁≧…≧Ｎ_k…≧Ｎ_K≧２Ｑとなる
Ｑ：窓の周期数，定Ｑ変換ではこれを全ての周波数で共通にする
非特許文献１、２ではＱ＝（２^1/B−１）^-1程度に取るようにしている
Ｂはビン数で、１オクターブを分割する数のこと
ｆ_s：サンプリング周波数
ｆ_k：第ｋ成分の周波数，ｆ_k＝（２^1/B）^k-1ｆ_min
ｆ_min＝ ｆ₁は定Ｑ変換で解析する最小の周波数で、解析する範囲から決める
ｆ_kは解析する最大の周波数で、ｆ_k≦ｆ_s／２となるようにする。

定Ｑ変換は特に低周波帯域の成分で窓幅Ｎ_kが非常に大きくなることにより演算量が多い。

これに対して非特許文献２では、高速フーリエ変換（ＦＦＴ；Ｆａｓｔ Ｆｏｕｒｉｅｒ ｔｒａｎｓｆｏｒｍ）を利用した演算の高速化方法が提案されている。

この方法ではパーセバルの公式により、定Ｑ変換を次の数式（２）のように時間軸での計算から周波数軸での計算に置き換える。

尚、明細書中の以下の文章では、前記数式（１）、（２）の各左辺の定Ｑ変換の成分を、「Ｘ_k ^cq」と表記することもある。

ここでＳ_n,kは（小さなｋに対して特に）ほとんどの要素が０に近いので、一定値以下の項を０として疎行列表現すると（実質非ゼロ項の積和計算だけになるため）非常に高速に計算できるようになる。

さらに（１／Ｎ）Ｓ_n,kは、計算するデータ長と解析する周波数範囲が決まっていれば、入力データ列ｘ_nに寄らず事前に計算できるためこれを予め用意しておけば、入力データの離散フーリエ変換はＦＦＴで実行するとして、定Ｑ変換を入力データに対する１回のＦＦＴと、ＦＦＴ結果ベクトルと予め用意した疎行列との積演算により計算できることになる。

離散でないフーリエ変換でも同等の公式が成り立ち、これらはパーセバルの定理、プランシュレルの定理など様々な呼び方がされる。

Ｊｕｄｉｔｈ Ｃ．Ｂｒｏｗｎ：Ｃａｌｃｕｌａｔｉｏｎ ｏｆ ａ ｃｏｎｓｔａｎｔ Ｑ ｓｐｅｃｔｒａｌ ｔｒａｎｓｆｏｒｍ，Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ．８９（１）：４２５−４３４，１９９１ Ｊｕｄｉｔｈ Ｃ．Ｂｒｏｗｎ ａｎｄ Ｍｉｌｌｅｒ Ｓ．Ｐｕｃｋｅｔｔｅ：Ａｎ ｅｆｆｉｃｉｅｎｔ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｔｈｅ ｃａｌｃｕｌａｔｉｏｎ ｏｆ ａ ｃｏｎｓｔａｎｔ Ｑ ｔｒａｎｓｆｏｒｍ，Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ．９２（５）：２６９８−２７０１，１９９２ ｙｕｋａｒａ＿１３：[Ｐｙｔｈｏｎ]Ｃｏｎｓｔａｎｔ−Ｑ変換（対数周波数スペクトログラム）、音楽プログラミングの超入門（仮） ｉｎ Ｈａｔｅｎａ Ｂｌｏｇ，２０１３−１２−０１，２０１３、インターネット＜ＵＲＬ：ｈｔｔｐ：／／ｙｕｋａｒａ−１３．ｈａｔｅｎａｂｌｏｇ．ｃｏｍ／ｅｎｔｒｙ／２０１３／１２／０１／２２２７４２＞．[２０１６−０２−２２ アクセス] ｙｕｋａｒａ＿１３：[Ｐｙｔｈｏｎ] 高速なＣｏｎｓｔａｎｔ−Ｑ変換（ｗｉｔｈ ＦＦＴ）、音楽プログラミングの超入門（仮） ｉｎ Ｈａｔｅｎａ Ｂｌｏｇ，２０１４−０１−０５，２０１４、インターネット＜ＵＲＬ：ｈｔｔｐ：／／ｙｕｋａｒａ−１３．ｈａｔｅｎａｂｌｏｇ．ｃｏｍ／ｅｎｔｒｙ／２０１４／０１／０５／０６２４１４＞．[２０１６−０２−２２ アクセス]

非特許文献２による高速化のアイデアは周波数軸での係数Ｓ_n,kが疎行列とみなせることに依拠しているが、実際にはＳ_n,kが疎行列というのは粗い評価である。確かにその値は一部を除いて小さいが、精度を高めると必ずしも疎にならない。特に高周波帯域でこの傾向が顕著である。実際、非特許文献２、４で係数Ｓ_n,kが十分小さいとしている閾値はそれぞれ０．１５、０．００５４であり、通常の評価には問題ないが有効数字７桁以上の精度を得るには十分とは言えない。このため、要求精度によってはこの高速化手法がかえって遅くなる可能性がある。

実際に周波数軸での計算でどの程度の演算量を必要とするかを調べるため、非特許文献４に準じた設定により時間軸での計算と周波数軸での計算でそれぞれ計算対象となる非ゼロの項数を評価する。

設定パラメータは、本明細書の記号法で表記してＢ＝２４、Ｑ＝２８、ｆ_s＝１６０００、ｆ_min＝６０、Ｋ＝１６０とし、定Ｑ変換をそれぞれ矩形窓、ハミング窓、ハニング窓で計算する場合において、時間軸での計算に必要な非ゼロ項数と、周波数軸での計算で所定の精度の係数まで計算する際に必要な項数を表１に示す。

なお、非特許文献１、２に合わせるとＢ＝２４のときＱ＝３４となるが、非特許文献３、４では追加の係数ｆｒａｔｉｏを導入してＱの値を変えているので、その設定値Ｑ＝２８を使っている。また解析する最大周波数６０００ＨｚからＫ＝１６０を導いた。ｆ_sの値は非特許文献３，４に明記されていなかったので、最大周波数６０００Ｈｚまで解析できるサンプリング周波数１２０００Ｈｚ以上のＷＡＶ音声データとしてｆ_s＝１６０００Ｈｚを選択した。

表１中の周波数軸右の数値は係数をゼロと評価する係数の大きさの閾値であり、各数値右のパーセンテージは矩形窓の時間軸計算の項数を基準としている。

また、非ゼロの項数と周波数の関係を、時間軸での計算の場合のグラフと周波数軸（閾値ごと）での計算の場合のグラフとして図５〜図８に示す。図５は矩形窓、図６はハミング窓、図７はハニング窓、図８は複素Ｍｏｒｌｅｔウェーブレット変換を用いて計算した場合を各々示している。図５〜図８において、各曲線の下部分の面積が演算量に相当する（ＦＦＴの演算は別である）。

表１から、窓形状により程度が違うが、精度が高いほど（係数をゼロと評価する係数Ｓ_n,kの大きさの閾値が小さいほど）周波数軸での計算に必要な項数が増えることが分かる。特に矩形窓やハミング窓では精度が高くなると時間軸での計算より周波数軸での計算の方が積和すべき項数が大きくなり、周波数軸での計算で高速にならなくなる。

尚、周波数軸での計算にはＦＦＴ計算と疎行列の非ゼロ係数抽出のオーバーヘッドもあるので、必要項数に差がなければ時間軸計算が有利である。

より精度の高い定Ｑ変換を高速に行うには工夫が必要になる。

本発明は上記課題を解決するものであり、その目的は、より精度の高い定Ｑ変換を、少ない演算量で高速に行うことができる定Ｑ変換の成分演算装置および方法を提供することにある。

上記課題を解決するための請求項１に記載の定Ｑ変換の成分演算装置は、時系列データを定比周波数成分へデータ変換する定Ｑ変換の成分演算装置であって、
設定した周波数成分Ｋ_M以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記Ｋ_Mを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、
前記Ｋ_Mを全周波数成分ｋに各々設定して各Ｋ_M毎の前記合計演算量を求め、
前記各Ｋ_M毎の合計演算量を評価し、該合計演算量が最小となるＫ_Mを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段と、
定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値Ｋ_M以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値Ｋ_Mを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段と、
時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段と、
を備えたことを特徴とする。

また、請求項５に記載の定Ｑ変換の成分演算方法は、時系列データを定比周波数成分へデータ変換する定Ｑ変換の成分演算方法であって、
計算法境界決定手段が、設定した周波数成分Ｋ_M以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記Ｋ_Mを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、前記Ｋ_Mを全周波数成分ｋに各々設定して各Ｋ_M毎の前記合計演算量を求めるステップと、
計算法境界決定手段が、前記各Ｋ_M毎の合計演算量を評価し、該合計演算量が最小となるＫ_Mを、計算法境界を示す周波数成分閾値として決定するステップと、
計算法選択手段が、定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値Ｋ_M以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値Ｋ_Mを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択するステップと、
定Ｑ変換手段が、時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換ステップと、
を備えたことを特徴とする。

また、請求項２に記載の定Ｑ変換の成分演算装置は、請求項１において、前記定Ｑ変換における時間軸係数を、

によって算出する時間軸係数算出部と、
前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である

から（１／Ｎ）Ｓ_n,kとして算出する周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。

また、請求項６に記載の定Ｑ変換の成分演算方法は、請求項５において、時間軸係数算出部が、前記定Ｑ変換における時間軸係数を、

によって算出するステップと、
周波数軸係数算出部が、前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である

から（１／Ｎ）Ｓ_n,kとして算出するステップと、を備え、
前記定Ｑ変換ステップは、
低周波帯域計算部が、前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施するステップと、
高周波帯域計算部が、前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施するステップと、
合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求めるステップと、
を備えていることを特徴とする。

定Ｑ変換の成分を演算する場合、低周波帯域では時系列データの離散フーリエ変換と周波数軸係数との積和計算の方が演算量が少なく、高周波帯域では時系列データと時間軸係数との積和計算の方が演算量が少ない。

このため上記構成によれば、低周波帯域、高周波帯域のどちらの場合でも演算量の少ない計算法で定Ｑ変換の成分を演算することができ、演算の高速化が実現される。

また、請求項３に記載の定Ｑ変換の成分演算装置は、請求項１において、

を、予め設定した時間方向の複数の解析位置ｐ_m（ｍ＝１，…，Ｍ）毎に算出して時間軸係数Ｔ_n,k,mを求め、該Ｔ_n,k,mを前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である

から得た（１／Ｎ）Ｓ_n,kを、前記解析位置ｐ_m毎に算出して周波数軸係数（１／Ｎ）Ｓ_n,k,mを求め、該（１／Ｎ）Ｓ_n,k,mを前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘ_n全体に対して離散フーリエ変換を行ってフーリエ係数Ｘ_nを求め、前記解析位置ｐ_m毎に、前記周波数軸係数算出部で算出された周波数軸係数（１／Ｎ）Ｓ_n,k,mと前記フーリエ係数Ｘ_nとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐ_m毎に、前記時間軸係数算出部で算出された時間軸係数Ｔ_n,k,mと時系列データｘ_nとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。

また、請求項７に記載の定Ｑ変換の成分演算方法は、請求項５において、時間軸係数算出部が、

を、予め設定した時間方向の複数の解析位置ｐ_m（ｍ＝１，…，Ｍ）毎に算出して時間軸係数Ｔ_n,k,mを求め、該Ｔ_n,k,mを前記定Ｑ変換における時間軸係数とするステップと、
周波数軸係数算出部が、時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である

から得た（１／Ｎ）Ｓ_n,kを、前記解析位置ｐ_m毎に算出して周波数軸係数（１／Ｎ）Ｓ_n,k,mを求め、該（１／Ｎ）Ｓ_n,k,mを前記定Ｑ変換における周波数軸係数とするステップと、を備え、
前記定Ｑ変換ステップは、
低周波帯域計算部が、時系列データｘ_n全体に対して離散フーリエ変換を行ってフーリエ係数Ｘ_nを求め、前記解析位置ｐ_m毎に、前記周波数軸係数算出部で算出された周波数軸係数（１／Ｎ）Ｓ_n,k,mと前記フーリエ係数Ｘ_nとを積和計算して前記第１の計算法を実施するステップと、
高周波帯域計算部が、前記解析位置ｐ_m毎に、前記時間軸係数算出部で算出された時間軸係数Ｔ_n,k,mと時系列データｘ_nとを積和計算して前記第２の計算法を実施するステップと、
合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求めるステップと、
を備えていることを特徴とする。

上記構成によれば、複数の時刻（複数の解析位置ｐ_m（ｍ＝１，…，Ｍ））で各々定Ｑ変換を行う際に、時間軸係数および周波数軸係数は複数の時刻毎（解析位置ｐ_m毎）に用意するが、低周波帯域における周波数軸での計算は、時系列データの離散フーリエ変換を、複数の時刻毎に行わなくても、ｘ_n全体に対して一括して離散フーリエ変換したフーリエ係数Ｘ_nと各解析位置での周波数軸係数との積和計算でよいため、全体としての計算が高速化される。

また、請求項４に記載の定Ｑ変換の成分演算装置は、請求項１において、

を、予め設定した時間方向の複数の解析位置ｐ_m（ｍ＝１，…，Ｍ）のうち１つの解析位置ｐ₁において算出して時間軸係数Ｔ_n,k,1を求め、該Ｔ_n,k,1を前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である

から得た（１／Ｎ）Ｓ_n,kを、前記１つの解析位置ｐ₁において算出して周波数軸係数（１／Ｎ）Ｓ_n,k,1を求め、該（１／Ｎ）Ｓ_n,k,1を前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘ_n全体に対して離散フーリエ変換を行ってフーリエ係数Ｘ_nを求め、 前記解析位置ｐ_m毎に、
前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓ_n,k,1と、
前記解析位置ｐ_m毎の周波数軸係数（１／Ｎ）Ｓ_n,k,mは前記１つの解析位置ｐ₁における周波数軸係数（１／Ｎ）Ｓ_n,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗

と、
前記フーリエ係数Ｘ_nとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐ_m毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔ_n,k,1と時系列データｘ_n-p1+pmとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。

また、請求項８に記載の定Ｑ変換の成分演算方法は、請求項５において、時間軸係数算出部が、

を、予め設定した時間方向の複数の解析位置ｐ_m（ｍ＝１，…，Ｍ）のうち１つの解析位置ｐ₁において算出して時間軸係数Ｔ_n,k,1を求め、該Ｔ_n,k,1を前記定Ｑ変換における時間軸係数とするステップと、
周波数軸係数算出部が、時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である

から得た（１／Ｎ）Ｓ_n,kを、前記１つの解析位置ｐ₁において算出して周波数軸係数（１／Ｎ）Ｓ_n,k,1を求め、該（１／Ｎ）Ｓ_n,k,1を前記定Ｑ変換における周波数軸係数とするステップと、を備え、
前記定Ｑ変換ステップは、
低周波帯域計算部が、時系列データｘ_n全体に対して離散フーリエ変換を行ってフーリエ係数Ｘ_nを求め、
前記解析位置ｐ_m毎に、
前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓ_n,k,1と、
前記解析位置ｐ_m毎の周波数軸係数（１／Ｎ）Ｓ_n,k,mは前記１つの解析位置ｐ₁における周波数軸係数（１／Ｎ）Ｓ_n,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗

と、
前記フーリエ係数Ｘ_nとを積和計算して前記第１の計算法を実施するステップと、
高周波帯域計算部が、前記解析位置ｐ_m毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔ_n,k,1と時系列データｘ_n-p1+pmとを積和計算して前記第２の計算法を実施するステップと、
合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求めるステップと、
を備えていることを特徴とする。

上記構成によれば、複数の時刻（複数の解析位置ｐ_m（ｍ＝１，…，Ｍ））で各々定Ｑ変換を行う際に、１つの解析位置ｐ₁における時間軸係数および周波数軸係数を利用して各時刻での計算が可能となるので、複数の時刻毎（複数の解析位置毎）に時間軸係数および周波数軸係数を用意する必要はなく、データ量が少なくなってメモリを節約することができる。

（１）請求項１〜８に記載の発明によれば、低周波帯域、高周波帯域のどちらの場合でも演算量の少ない計算法で定Ｑ変換の成分を演算することができ、演算の高速化が実現される。
（２）請求項３、７に記載の発明によれば、複数の時刻で各々定Ｑ変換を行う際に、低周波帯域における周波数軸での計算は、時系列データの離散フーリエ変換を、複数の時刻毎に行わなくても、ｘ_n全体に対して一括して離散フーリエ変換したフーリエ係数Ｘ_nと各解析位置での周波数軸係数との積和計算でよいため、全体としての計算が高速化される。
（３）請求項４、８に記載の発明によれば、複数の時刻で各々定Ｑ変換を行う際に、１つの解析位置ｐ₁における時間軸係数および周波数軸係数を利用して各時刻での計算が可能となるので、複数の時刻毎（複数の解析位置毎）に時間軸係数および周波数軸係数を用意する必要はなく、データ量が少なくなってメモリを節約することができる。

本発明の実施例による定Ｑ変換の成分演算装置の構成図。 本発明の実施例２における時間軸・周波数軸の非ゼロ項数を示すグラフ。 本発明の実施例に１おける時間軸・周波数軸の非ゼロ項数を示すグラフ。 従来手法と本発明の手法での定Ｑ変換の複数時刻計算の演算量を示すグラフ。 従来手法の矩形窓における時間軸・周波数軸の項数を示すグラフ。 従来手法のハミング窓における時間軸・周波数軸の項数を示すグラフ。 従来手法のハニング窓における時間軸・周波数軸の項数を示すグラフ。 複素Ｍｏｒｌｅｔウェーブレットの計算項数を示すグラフ。

以下、図面を参照しながら本発明の実施の形態を説明するが、本発明は下記の実施形態例に限定されるものではない。

本実施形態例では、定Ｑ変換を高精度でも演算量を少なく計算するために、定Ｑ変換の成分Ｘ_k ^cqを数式（２）の時間軸での計算と周波数軸での計算のどちらの方法で計算するか、予め周波数成分ｋ毎に演算量の小さい方を選択しておくことで実現する。

定Ｑ変換の成分Ｘ_k ^cqの計算はおおよそ低周波帯域では周波数軸での計算の演算量が小さく、高周波帯域では時間軸での計算の演算量が小さいから、ｋの境界値Ｋ_Mを決めて、ｋ≦Ｋ_Mでは周波数軸での計算、ｋ＞Ｋ_Mでは時間軸での計算、というようにすれば良い。

図１は本発明の実施例１による定Ｑ変換の成分演算装置の構成を示している。図１において、１０は予め定Ｑ変換を行うためのパラメータを決定して定Ｑ変換計算の準備をする準備部であり、２０は個別の時系列データに対する定Ｑ変換を計算する個別計算部である。

準備部１０は、
設定した周波数成分Ｋ_M以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記Ｋ_Mを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、前記Ｋ_Mを全周波数成分ｋに各々設定して各Ｋ_M毎の前記合計演算量を求め、前記各Ｋ_M毎の合計演算量を評価し、該演算量が最小となるＫ_Mを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段の機能と、
定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値Ｋ_M以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値Ｋ_Mを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段の機能とを具備している。

前記計算法境界決定手段および計算法選択手段の各機能を実行する準備部１０は、具体的には、定Ｑ変換のパラメータを決定するパラメータ決定部１１、時間軸係数Ｔ_n,kを計算する時間軸係数算出部１２、周波数軸係数（１／Ｎ）Ｓ_n,kを計算する周波数軸係数算出部１３、時間軸での計算と周波数軸での計算との境界のｋ値Ｋ_Mを決定する計算法境界決定部１４からなる。

個別計算部２０は、時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段の機能を具備している。

前記定Ｑ変換手段の機能を実行する個別計算部２０は、具体的には、時系列データの長さを整える前処理部２１、高周波帯域の定Ｑ変換を計算する高周波帯域計算部２２、低周波帯域の定Ｑ変換を計算する低周波帯域計算部２３、高周波帯域での結果と低周波帯域での結果をまとめる定Ｑ変換合成部２４からなる。

図１の定Ｑ変換の成分演算装置は、例えばコンピュータにより構成され、通常のコンピュータのハードウェアリソース、例えばＲＯＭ、ＲＡＭ、ＣＰＵ、入力装置、出力装置、通信インターフェース、ハードディスク、記録媒体およびその駆動装置を備えている。

このハードウェアリソースとソフトウェアリソース（ＯＳ、アプリケーションなど）との協働の結果、定Ｑ変換の成分演算装置は、図１に示すように、パラメータ決定部１１、時間軸係数算出部１２、周波数軸係数算出部１３、計算法境界決定部１４、前処理部２１、高周波帯域計算部２２、低周波帯域計算部２３、定Ｑ変換合成部２４を実装する。

また各部での処理結果は、ハードディスク或いはＲＡＭなどの保存手段・記憶手段に格納され、随時利用されるものとする。

次に、図１の各部で実施される内容を詳述する。

まず、準備部１０のパラメータ決定部１１では、予め定Ｑ変換で解析する時系列データのサンプリング周波数ｆ_s、解析する周波数範囲（ｆ_min，ｆ_max）と１オクターブの周波数を区切るビン数Ｂ、解析対象データの周期数Ｑ、窓関数を決定する。このときｆ_max≦ｆ_s／２が必要で、Ｋ＝ｆｌｏｏｒ（Ｂ＊ｌｏｇ₂（ｆ_max／ｆ_min））＋１となる。これを基に計算する入力データの長さをＮ＝Ｎ₁＝Ｑ（ｆ_s／ｆ_min）とする。

また周波数軸係数をゼロとみなすべき係数の大きさの閾値Ｔｈもここで決定しておく。

時間軸係数算出部１２では、パラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔ_n,kを次の定義式に従って計算する。

周波数軸係数算出部１３では、パラメータ決定部１１で決定したパラメータを基に、周波数軸係数（１／Ｎ）Ｓ_n,kを、

の定義式に従って計算する。この際、｜Ｓ_n,k｜≦Ｔｈとなる項は周波数軸係数（１／Ｎ）Ｓ_n,k＝０とする。

計算法境界決定部１４では、周波数軸での計算を行う最大のｋ値Ｋ_M（周波数成分閾値）を求める。Ｋ_Mは０，…，Ｋの何れかの値を持ち、Ｋ_M＝０は周波数軸での計算を行わずに非特許文献１、２と同様に全て時間軸での計算で求める場合を意味する。

周波数成分閾値Ｋ_Mの決定方法としては、Ｋ_M＝０，…，Ｋの各場合について、定Ｑ変換の成分Ｘ_k ^cqを求めるのに必要な演算量Ｃ_X（Ｋ_M）を評価してそれが最小となるようにＫ_Mを決定する。

定Ｑ変換の成分Ｘ_k ^cqの演算量Ｃ_X（Ｋ_M）（合計演算量）は、長さＮのＦＦＴの演算量Ｃ_FFT（Ｎ）、非ゼロの時間軸係数毎の演算量Ｃ_T、非ゼロの時間軸係数の数Ｎ_T（ｋ）、非ゼロの周波数軸係数毎の演算量Ｃ_F、ｋで非ゼロとみなす周波数軸係数の数Ｎ_F（ｋ）を使っておおよそ

と表せる。

ここでＣ_FFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、Ｃ_FFT（Ｎ）＝（１／４）Ｎｌｏｇ₂（Ｎ）、Ｃ_T＝１、Ｃ_F＝１としても良い。

上記演算量Ｃ_X（Ｋ_M）を示す数式（３）において、右辺第１項は、Ｋ_Mが非ゼロのときＣ_FFTとし、そうでないときは０とすることを表し、右辺第２項は、周波数成分ｋが周波数成分閾値Ｋ_M以下の低周波帯域における周波数軸係数毎の演算量を表し、右辺第３項は、周波数成分ｋが周波数成分閾値Ｋ_Mを超える高周波帯域における時間軸係数毎の演算量を表している。

このためＫ_Mが非ゼロのときは、数式（３）の右辺第１項のＦＦＴの演算量および右辺第２項の演算量（周波数軸係数による演算量）からなる第１の演算量と、数式（３）の右辺第３項の時間軸係数による第２の演算量との合計が演算量Ｃ_X（Ｋ_M）となる。

またＫ_Mがゼロのときは、数式（３）の右辺第１項および右辺第２項がともに０となって右辺第３項の時間軸係数による第２の演算量のみが演算量Ｃ_X（Ｋ_M）となる。

この演算量Ｃ_X（Ｋ_M）は、時間軸係数および周波数軸係数の、ハニング窓での項数と周波数（周波数成分ｋ；対数周波数）の関係を示す図３のように、非ゼロの時間軸係数の数Ｎ_T（ｋ）がｋに関して指数的に減少し、非ゼロの周波数軸係数の数Ｎ_F（ｋ）がｋに関しておおよそ指数的に増加するため、Ｋ_M＝１，…，Ｋの範囲ではおおよそ下に凸の特性カーブとなる。ただし、Ｋ_M＝０ではＣ_FFTがなくなるためＣ_X（０）が小さくなる。結果として演算量Ｃ_X（Ｋ_M）はＫ_M＝０あるいはＫ_M＝１，…，Ｋでの最小点の何れかで最小値を取ることになる。

個別計算部２０では、入力された時系列データの定Ｑ変換を次のようにして求める。

前処理部２１では、入力された時系列データの長さがＮより短ければ前後に０を補完して、Ｎより長ければ前後の値を除去して長さがＮの時系列データｘ_nになるようにする。このとき入力された時系列データの中央部分が時系列データｘ_nの中央部分に一致するようにする。

次に定Ｑ変換の成分Ｘ_k ^cqを、前記計算法境界決定部１４で決定した周波数成分閾値Ｋ_Mにより、ｋ≦Ｋ_Mについては低周波帯域計算部２３で計算（第１の計算法）し、ｋ＞Ｋ_Mについては高周波帯域計算部２２で計算（第２の計算法）する。

高周波帯域計算部２２は、前記時間軸係数算出部１２で算出された時間軸係数Ｔ_n,kと時系列データｘ_nを積和計算して（第２の計算法を実施して）定Ｑ変換の成分Ｘ_k ^cqを求める。各ｋについて非ゼロの時間軸係数Ｔ_n,kは高々Ｎ_k個なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。

低周波帯域計算部２３は、まず時系列データｘ_nをＦＦＴしてフーリエ係数Ｘ_nを求める。次に前記周波数軸係数算出部１３で算出された周波数軸係数（１／Ｎ）Ｓ_n,kとフーリエ係数Ｘ_nを積和計算して（第１の計算法を実施して）定Ｑ変換の成分Ｘ_k ^cqを求める。各ｋについて非ゼロの周波数軸係数（１／Ｎ）Ｓ_n,kは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。

尚、周波数軸係数は多くの場合、両端（周波数ゼロ近辺とサンプリング周波数近辺）部分では非ゼロ項とゼロ項が入り混じるが中間部分ではゼロ項が続くので、非ゼロ項の範囲を少し広げて、両端部分は全て非ゼロ項として計算しても良い。この場合、図３の細かい点線で示す周波数軸簡易版の特性カーブのように、非ゼロ項が少し増えるがその数はそれほど多くはなく、また非ゼロ判定ロジックが簡略化できるので演算時間は大差なくなる。

最後に定Ｑ変換合成部２４で、高周波帯域計算部２２と低周波帯域計算部２３とで計算した定Ｑ変換の成分Ｘ_k ^cqをまとめて、定Ｑ変換を完了する。

上記のように本実施例１によれば、従来の高速化技術の問題点である、高い精度の定Ｑ変換を計算しようとすると演算量が膨れ上り通常の時間軸での計算よりもむしろ遅くなるという問題が改善される。

その効果を示すため、定Ｑ変換の計算方法（時間軸計算、周波数軸計算、提案手法）毎に、必要な演算量を表２のように試算した。試算は、時系列データのパラメータを表１と同じにして、演算量パラメータをＣ_FFT（Ｎ）＝（１／４）Ｎｌｏｇ₂（Ｎ）、Ｃ_T＝１、Ｃ_F＝１として計算した定Ｑ変換の成分Ｘ_k ^cqの演算量Ｃ_X（Ｋ_M）を表２に示す。

表２中の周波数軸右の数値は係数をゼロと評価する係数の大きさの閾値であり、各数値右のパーセンテージは矩形窓の時間軸計算の演算量を基準としている。

表２では、ＦＦＴ演算の演算量が原因で高速化の効果が小さめに出ているが、ハニング窓では閾値１．０ｅ−５でもおおよそ１／５（２２．５％）の演算量になっている。またＦＦＴ演算の演算量が重く周波数軸での計算で効果が出ない場合でも、自動的に時間軸での計算になるので反って遅くなるという問題は発生しない。

実施例１は、ある一時刻における定Ｑ変換の各成分値を求める際には効率的であるが、実際の定Ｑ変換では、時系列データの周波数成分の時間変化を捉えるために少しずつ時刻をずらしながら繰り返し定Ｑ成分を計算する。このときずらす時刻は、最短ではサンプリング周波数で、時系列データの情報を余すことなく計算するには最長でも最高周波数での窓幅にしたい。もちろん、時系列データの特徴サンプル抽出などでは、これよりも長い間隔で時刻をずらして定Ｑ変換する場合もあるが、いずれにせよ時系列データ全体に対しては多くの繰り返し計算をすることになり、まだまだ負荷が高い。

多数の時刻で同時に定Ｑ変換を計算する際には更なる効率化が必要である。

そこで本実施例２では、数式（２）の周波数軸での計算で定Ｑ変換を行う場合に、定Ｑ変換を計算する時刻ごとに計算に使用する時系列データの範囲でＤＦＴ（あるいはＦＦＴ）する代わりに、複数時刻での定Ｑ変換に必要な範囲全体で一括してＤＦＴ（あるいはＦＦＴ）を行うことで高速化を図った。

数式（１）で定義する定Ｑ変換はデータ列ｘ_n，ｎ＝０，…，Ｎ−１の中央ｘ_N/2での成分Ｘ_k ^cqを計算するように定式化しているが、窓関数の非ゼロ範囲を時刻方向に平行移動することで他の時刻での成分を計算するように調整できる。すなわちＴ_n,kとｘ_nとの積和のインデックスをずらすことで実現できる。

数式（２）による周波数軸での計算に関しては、初めに定Ｑ変換する時刻全体で必要な全データに対して一括してＤＦＴ（あるいはＦＦＴ）を行い、予め各時刻用に用意した（１／Ｎ）Ｓ_n,kと順次積和することで実現する。

本実施例２による定Ｑ変換の成分演算装置の具体的な計算の構成は、図１と同様であるが、各部で実行される機能が以下に示すように異なる。

まず、準備部１０のパラメータ決定部１１では、予め定Ｑ変換で解析する時系列データのサンプリング周波数ｆ_s、解析する周波数範囲（ｆ_min，ｆ_max）と１オクターブの周波数を区切るビン数Ｂ、解析対象データの周期数Ｑ、窓関数を決定する。このときｆ_max≦ｆ_s／２が必要で、Ｋ＝ｆｌｏｏｒ（Ｂ＊ｌｏｇ₂（ｆ_max／ｆ_min））＋１となる。

さらに入力データｘ_nの全長Ｎと、解析位置ｐ_m，ｍ＝１，…，Ｍ（時間方向の解析位置；時刻）を決定しておく。このとき、各位置ｐ_mでの定Ｑ変換はｘ_pmを中心とする長さＮ₁＝Ｑ（ｆ_s／ｆ_min）のデータ列をもとに計算し、（Ｎ₁／２）≦ｐ₁≦…≦ｐ_M≦Ｎ−（Ｎ₁／２）であるものとする。すなわち、全ての解析位置ｐ_m，ｍ＝１，…，Ｍで定Ｑ変換の計算に必要なデータがｘ₀，…，ｘ_N-1に含まれるものとする。尚、入力データの全点で定Ｑ変換を計算したいような場合は予め入力データを前後にゼロ拡張して上記条件を満たすように変換しておけばよい。

また、周波数軸係数をゼロとみなすべき係数の大きさの閾値Ｔｈもここで決定しておく。

時間軸係数算出部１２ではパラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔ_n,kを各解析位置ｐ_mについて次の定義式に従って計算する。

このときｐ_mに対するＴ_n,kをＴ_n,k,mと書くと、非ゼロ項に関してＴ_n,k,m＝Ｔ_n-pm+p1,k,1である。

周波数軸係数算出部１３ではパラメータ決定部１１で決定したパラメータを基に、各解析位置ｐ_mについて周波数軸係数（１／Ｎ）Ｓ_n,k,mを定義式に従って計算する。この際、｜Ｓ_n,k,m｜≦Ｔｈとなる項は周波数軸係数（１／Ｎ）Ｓ_n,k,m＝０とする。

計算法境界決定部１４では、周波数軸での計算を行う最大のｋ値Ｋ_M（周波数成分閾値）を求める。Ｋ_Mは０，…，Ｋの何れかの値を持ち、Ｋ_M＝０は周波数軸での計算を行わずに非特許文献１、２と同様に全て時間軸での計算で求める場合を意味する。

周波数成分閾値Ｋ_Mの決定方法としては、Ｋ_M＝０，…，Ｋの各場合について、定Ｑ変換の成分Ｘ_k ^cqを求めるのに必要な演算量Ｃ_X（Ｋ_M）を評価してそれが最小となるようにＫ_Mを決定する。尚、この条件は複数のｐ_mに寄らないので、ｐ₁（１つの解析位置）について求めておけば十分である。

定Ｑ変換の成分Ｘ_k ^cqの演算量Ｃ_X（Ｋ_M）（合計演算量）は、長さＮのＦＦＴの演算量Ｃ_FFT（Ｎ）、非ゼロの時間軸係数毎の演算量Ｃ_T、非ゼロの時間軸係数の数Ｎ_T（ｋ）、非ゼロの周波数軸係数毎の演算量Ｃ_F、ｋで非ゼロとみなす周波数軸係数の数Ｎ_F（ｋ）を使っておおよそ、前記実施例１で述べた数式（３）

と表せる。

ここでＣ_FFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、Ｃ_FFT（Ｎ）＝（１／４）Ｎｌｏｇ₂（Ｎ）、Ｃ_T＝１、Ｃ_F＝１としても良い。

この演算量Ｃ_X（Ｋ_M）は、時間軸係数および周波数軸係数の、ハニング窓での項数と周波数（周波数成分ｋ；対数周波数）の関係を示す図２のように、非ゼロの時間軸係数の数Ｎ_T（ｋ）がｋに関して指数的に減少し、非ゼロの周波数軸係数の数Ｎ_F（ｋ）がｋに関しておおよそ指数的に増加するため、Ｋ_M＝１，…，Ｋの範囲ではおおよそ下に凸の特性カーブとなる。ただし、Ｋ_M＝０ではＣ_FFTがなくなるためＣ_X（０）が小さくなる。結果として演算量Ｃ_X（Ｋ_M）はＫ_M＝０あるいはＫ_M＝１，…，Ｋでの最小点の何れかで最小値を取ることになる。

個別計算部２０では入力された時系列の定Ｑ変換を次のようにして求める。

前処理部２１では、入力された時系列データの長さがＮより短ければ前後にゼロを補完して、Ｎより長ければ前後の値を除去して長さがＮの時系列データｘ_nになるようにする。ゼロを補完する場合には、入力された時系列データの中央部分が時系列データｘ_nの中央部分に一致するようにする。データを除去する場合は、解析したい範囲に寄るが、単純に後ろのデータを除去するのでも良い。

次に定Ｑ変換の成分Ｘ_k ^cqを、計算法境界決定部１４で決定したＫ_Mにより、ｋ≦Ｋ_Mについては低周波帯域計算部２３で、ｋ＞Ｋ_Mについては高周波帯域計算部２２で計算する。

尚、前記定Ｑ変換の成分は以下、「Ｘ_k,m ^cq」と表記することもある。

各ｋについて、前記時間軸係数算出部１２で算出された非ゼロの時間軸係数Ｔ_n,k,mは高々Ｎ_k個なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。

低周波帯域計算部２３は、まず時系列データｘ_nを全体に対して一度ＦＦＴしてフーリエ係数Ｘ_nを求める。次にｐ_m毎に、前記周波数軸係数算出部１３で算出された周波数軸係数（１／Ｎ）Ｓ_n,k,mとフーリエ係数Ｘ_nを積和計算して定Ｑ変換の成分Ｘ_k,m ^cqを求める。

各ｋについて非ゼロの周波数軸係数（１／Ｎ）Ｓ_n,k,mは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。

尚、周波数軸係数は多くの場合、両端（周波数ゼロ近辺とサンプリング周波数近辺）部分では非ゼロ項とゼロ項が入り混じるが中間部分ではゼロ項が続くので、非ゼロ項の範囲を少し広げて、両端部分は全て非ゼロ項として計算しても良い。この場合、図２の点線で示す周波数軸簡易版の特性カーブのように、非ゼロ項が少し増えるがその数はそれほど多くはなく、また非ゼロ判定ロジックが簡略化できるので演算時間は大差なくなる。

最後に定Ｑ変換合成部２４で、高周波帯域計算部２２と低周波帯域計算部２３とで計算した定Ｑ変換の成分Ｘ_k,m ^cqをまとめて、定Ｑ変換を完了する。

前記実施例２では、時間軸係数Ｔ_n,k,m、周波数軸係数（１／Ｎ）Ｓ_n,k,mは解析位置ｐ_m毎に用意するので、事前に用意して保管しておくデータ量が大きくなる。本実施例３では、係数データのメモリ使用量を大幅に減らすように以下のような対策を講じた。

尚、数式（４）中の

は、複素指数関数の時刻および周波数毎に異なるべき乗を表している。

本実施例３による定Ｑ変換の成分演算装置の具体的な構成は図１と同様であり、以下に図１と異なる部分のみを説明する。

時間軸係数算出部１２では、パラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔ_n,k,1を定義式に従って求める。すなわち、

を、予め設定した解析位置ｐ_mのうち１つの解析位置（例えば１番目の解析位置）ｐ₁において算出し、時間軸係数Ｔ_n,k,1を求める。

周波数軸係数算出部１３では、パラメータ決定部１１で決定したパラメータを基に、周波数軸係数（１／Ｎ）Ｓ_n,k,1を定義式に従って計算する。

すなわち、

から得た（１／Ｎ）Ｓ_n,kを、前記１つの解析位置ｐ₁において算出して周波数軸係数（１／Ｎ）Ｓ_n,k,1を求める。

この際、｜Ｓ_n,k,1｜≦Ｔｈとなる項は周波数軸係数（１／Ｎ）Ｓ_n,k,1＝０とする。

計算法境界決定部１４では、周波数軸計算の際に複素べき乗の積和が追加されるために非ゼロの周波数軸係数毎の演算量Ｃ_Fを実施例２の場合より大きく、例えば５（複素べき乗の演算に＋３、複素べき乗の積和に＋１）とする。

高周波帯域計算部２２は、複数の解析位置ｐ_m毎に範囲をずらして、前記時間軸係数算出部１２で算出された時間軸係数Ｔ_n,k,1と時系列データｘ_n-p1+pmを積和計算して定Ｑ変換の成分Ｘ_k,m ^cqを求める。

低周波帯域計算部２３は、まず時系列データｘ_nを全体に対して一度ＦＦＴしてフーリエ係数Ｘ_nを求める。

以上、定Ｑ変換の高速計算方式を説明したが、他の定比周波数成分へのデータ変換方法、例えばウェーブレット変換も計算の定式化はほとんど同じ形式にすることが可能であり、本発明と同様の計算方法を適用することは可能である。しかし周波数軸係数が疎行列に近くなるか否かは使用するウェーブレットに依存するため、必ずしも高速計算方式になるものではない。ただし例えば複素Ｍｏｒｌｅｔ Ｗａｖｅｌｅｔに関しては、少なくとも一つの時刻について計算する場合は高速化が期待できる。

前記実施例２、３の手法によれば、定Ｑ変換を高速に実行するにあたって、先行の高速化技術（実施例１）が時間軸での計算でも周波数軸での計算でもそれらを併用する計算方法でも、一時刻毎に定Ｑ変換を計算するのに対して、複数時刻での定Ｑ変換を同時に計算する際に、ＤＦＴ（あるいはＦＦＴ）を一時刻の計算範囲毎ではなく全体に対して一度だけ実施することで、全体としての計算が高速化される。

実施例３では、時間軸と周波数軸との係数データを計算時刻毎に持つ代わりに、それらを複素数のべき乗計算に置き換えて（数式（４））、演算量は少し増える代わりに係数データのメモリ使用量を大幅に減らしている。

効果を示すために、定Ｑ変換を時間軸計算で時刻毎に逐次計算した場合、実施例１の手法で逐次計算した場合、実施例２の手法で一括計算した場合、実施例３の方法でメモリを節約する場合の各方法にて演算量を評価した。尚、評価における設定パラメータは以下の通りである。

また、演算量のパラメータは実機での計測から以下の通りである。

これらを基に、１から１０００個の時刻で定Ｑ変換を計算した際に掛かる演算量を試算してグラフにすると図４のようになる。

図４より、実施例１は時間軸計算より高速であり、少なくとも１０時刻以上を一括で計算する場合には実施例２、３がさらに高速になることが分かる。

おおむね１００時刻以上の一括計算では、実施例２で実施例１の１０倍、実施例３で実施例１の５倍の高速化が見込める。

尚、時間軸係数Ｔ_n,k,1の非ゼロ項数は実施例１でも実施例２、３でも変わりないが、周波数軸係数（１／Ｎ）Ｓ_n,k,1の非ゼロ項数については、実施例１と実施例２、３でＦＦＴするデータ数（＝行列（１／Ｎ）Ｓ_n,k,1の次元）が異なるために、実施例１の場合より実施例２、３の場合の方が倍以上多くなる。また一回のＦＦＴに掛かる時間も同様に大幅に大きくなる。これらはもちろん図４に反映している。

それでも実施例２、３の手法が高速になるのはＦＦＴの実行回数が大きく異なるためである。

１０…準備部
１１…パラメータ決定部
１２…時間軸係数算出部
１３…周波数軸係数算出部
１４…計算法境界決定部
２０…個別計算部
２１…前処理部
２２…高周波帯域計算部
２３…低周波帯域計算部
２４…定Ｑ変換合成部

Claims (8)

1. 時系列データを定比周波数成分へデータ変換する定Ｑ変換の成分演算装置であって、
   設定した周波数成分Ｋ_M以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記Ｋ_Mを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、
   前記Ｋ_Mを全周波数成分ｋに各々設定して各Ｋ_M毎の前記合計演算量を求め、
   前記各Ｋ_M毎の合計演算量を評価し、該合計演算量が最小となるＫ_Mを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段と、
   定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値Ｋ_M以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値Ｋ_Mを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段と、
   時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段と、
   を備えたことを特徴とする定Ｑ変換の成分演算装置。
2. 前記定Ｑ変換における時間軸係数を、
   

   

   によって算出する時間軸係数算出部と、
   前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
   

   

   から（１／Ｎ）Ｓ_n,kとして算出する周波数軸係数算出部と、を備え、
   前記定Ｑ変換手段は、
   前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
   前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
   前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
   を備えていることを特徴とする請求項１に記載の定Ｑ変換の成分演算装置。
3. 

   を、予め設定した時間方向の複数の解析位置ｐ_m（ｍ＝１，…，Ｍ）毎に算出して時間軸係数Ｔ_n,k,mを求め、該Ｔ_n,k,mを前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
   時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
   

   

   から得た（１／Ｎ）Ｓ_n,kを、前記解析位置ｐ_m毎に算出して周波数軸係数（１／Ｎ）Ｓ_n,k,mを求め、該（１／Ｎ）Ｓ_n,k,mを前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
   前記定Ｑ変換手段は、
   時系列データｘ_n全体に対して離散フーリエ変換を行ってフーリエ係数Ｘ_nを求め、前記解析位置ｐ_m毎に、前記周波数軸係数算出部で算出された周波数軸係数（１／Ｎ）Ｓ_n,k,mと前記フーリエ係数Ｘ_nとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
   前記解析位置ｐ_m毎に、前記時間軸係数算出部で算出された時間軸係数Ｔ_n,k,mと時系列データｘ_nとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
   前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
   を備えていることを特徴とする請求項１に記載の定Ｑ変換の成分演算装置。
4. 

   を、予め設定した時間方向の複数の解析位置ｐ_m（ｍ＝１，…，Ｍ）のうち１つの解析位置ｐ₁において算出して時間軸係数Ｔ_n,k,1を求め、該Ｔ_n,k,1を前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
   時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
   

   

   から得た（１／Ｎ）Ｓ_n,kを、前記１つの解析位置ｐ₁において算出して周波数軸係数（１／Ｎ）Ｓ_n,k,1を求め、該（１／Ｎ）Ｓ_n,k,1を前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
   前記定Ｑ変換手段は、
   時系列データｘ_n全体に対して離散フーリエ変換を行ってフーリエ係数Ｘ_nを求め、 前記解析位置ｐ_m毎に、
   前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓ_n,k,1と、
   前記解析位置ｐ_m毎の周波数軸係数（１／Ｎ）Ｓ_n,k,mは前記１つの解析位置ｐ₁における周波数軸係数（１／Ｎ）Ｓ_n,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗
   

   

   と、
   前記フーリエ係数Ｘ_nとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
   前記解析位置ｐ_m毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔ_n,k,1と時系列データｘ_n-p1+pmとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
   前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
   を備えていることを特徴とする請求項１に記載の定Ｑ変換の成分演算装置。
5. 時系列データを定比周波数成分へデータ変換する定Ｑ変換の成分演算方法であって、
   計算法境界決定手段が、設定した周波数成分Ｋ_M以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記Ｋ_Mを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、前記Ｋ_Mを全周波数成分ｋに各々設定して各Ｋ_M毎の前記合計演算量を求めるステップと、
   計算法境界決定手段が、前記各Ｋ_M毎の合計演算量を評価し、該合計演算量が最小となるＫ_Mを、計算法境界を示す周波数成分閾値として決定するステップと、
   計算法選択手段が、定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値Ｋ_M以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値Ｋ_Mを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択するステップと、
   定Ｑ変換手段が、時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換ステップと、
   を備えたことを特徴とする定Ｑ変換の成分演算方法。
6. 時間軸係数算出部が、前記定Ｑ変換における時間軸係数を、
   

   

   によって算出するステップと、
   周波数軸係数算出部が、前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
   

   

   から（１／Ｎ）Ｓ_n,kとして算出するステップと、を備え、
   前記定Ｑ変換ステップは、
   低周波帯域計算部が、前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施するステップと、
   高周波帯域計算部が、前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施するステップと、
   合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求めるステップと、
   を備えていることを特徴とする請求項５に記載の定Ｑ変換の成分演算方法。
7. 時間軸係数算出部が、
   

   

   を、予め設定した時間方向の複数の解析位置ｐ_m（ｍ＝１，…，Ｍ）毎に算出して時間軸係数Ｔ_n,k,mを求め、該Ｔ_n,k,mを前記定Ｑ変換における時間軸係数とするステップと、
   周波数軸係数算出部が、時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
   

   

   から得た（１／Ｎ）Ｓ_n,kを、前記解析位置ｐ_m毎に算出して周波数軸係数（１／Ｎ）Ｓ_n,k,mを求め、該（１／Ｎ）Ｓ_n,k,mを前記定Ｑ変換における周波数軸係数とするステップと、を備え、
   前記定Ｑ変換ステップは、
   低周波帯域計算部が、時系列データｘ_n全体に対して離散フーリエ変換を行ってフーリエ係数Ｘ_nを求め、前記解析位置ｐ_m毎に、前記周波数軸係数算出部で算出された周波数軸係数（１／Ｎ）Ｓ_n,k,mと前記フーリエ係数Ｘ_nとを積和計算して前記第１の計算法を実施するステップと、
   高周波帯域計算部が、前記解析位置ｐ_m毎に、前記時間軸係数算出部で算出された時間軸係数Ｔ_n,k,mと時系列データｘ_nとを積和計算して前記第２の計算法を実施するステップと、
   合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求めるステップと、
   を備えていることを特徴とする請求項５に記載の定Ｑ変換の成分演算方法。
8. 時間軸係数算出部が、
   

   

   を、予め設定した時間方向の複数の解析位置ｐ_m（ｍ＝１，…，Ｍ）のうち１つの解析位置ｐ₁において算出して時間軸係数Ｔ_n,k,1を求め、該Ｔ_n,k,1を前記定Ｑ変換における時間軸係数とするステップと、
   周波数軸係数算出部が、時間軸係数Ｔ_n,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
   

   

   から得た（１／Ｎ）Ｓ_n,kを、前記１つの解析位置ｐ₁において算出して周波数軸係数（１／Ｎ）Ｓ_n,k,1を求め、該（１／Ｎ）Ｓ_n,k,1を前記定Ｑ変換における周波数軸係数とするステップと、を備え、
   前記定Ｑ変換ステップは、
   低周波帯域計算部が、時系列データｘ_n全体に対して離散フーリエ変換を行ってフーリエ係数Ｘ_nを求め、
   前記解析位置ｐ_m毎に、
   前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓ_n,k,1と、
   前記解析位置ｐ_m毎の周波数軸係数（１／Ｎ）Ｓ_n,k,mは前記１つの解析位置ｐ₁における周波数軸係数（１／Ｎ）Ｓ_n,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗
   

   

   と、
   前記フーリエ係数Ｘ_nとを積和計算して前記第１の計算法を実施するステップと、
   高周波帯域計算部が、前記解析位置ｐ_m毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔ_n,k,1と時系列データｘ_n-p1+pmとを積和計算して前記第２の計算法を実施するステップと、
   合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求めるステップと、
   を備えていることを特徴とする請求項５に記載の定Ｑ変換の成分演算方法。

Images