(57)【特許請求の範囲】
【請求項１】
時系列データを定比周波数成分へデータ変換する定Ｑ変換の成分演算装置であって、
設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、
前記ＫMを全周波数成分ｋに各々設定して各ＫM毎の前記合計演算量を求め、
前記各ＫM毎の合計演算量を評価し、該合計演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段と、
定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段と、
時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段と、
を備えたことを特徴とする定Ｑ変換の成分演算装置。
【請求項２】
前記定Ｑ変換における時間軸係数を、
【数４】
0006677069.tif 000034
によって算出する時間軸係数算出部と、
前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【数５】
0006677069.tif 000035
から（１／Ｎ）Ｓn,kとして算出する周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする請求項１に記載の定Ｑ変換の成分演算装置。
【請求項３】
【数４】
0006677069.tif 000036
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）毎に算出して時間軸係数Ｔn,k,mを求め、該Ｔn,k,mを前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【数５】
0006677069.tif 000037
から得た（１／Ｎ）Ｓn,kを、前記解析位置ｐm毎に算出して周波数軸係数（１／Ｎ）Ｓn,k,mを求め、該（１／Ｎ）Ｓn,k,mを前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、前記解析位置ｐm毎に、前記周波数軸係数算出部で算出された周波数軸係数（１／Ｎ）Ｓn,k,mと前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,mと時系列データｘnとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする請求項１に記載の定Ｑ変換の成分演算装置。
【請求項４】
【数４】
0006677069.tif 000038
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）のうち１つの解析位置ｐ1において算出して時間軸係数Ｔn,k,1を求め、該Ｔn,k,1を前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【数５】
0006677069.tif 000039
から得た（１／Ｎ）Ｓn,kを、前記１つの解析位置ｐ1において算出して周波数軸係数（１／Ｎ）Ｓn,k,1を求め、該（１／Ｎ）Ｓn,k,1を前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、  前記解析位置ｐm毎に、
前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓn,k,1と、
前記解析位置ｐm毎の周波数軸係数（１／Ｎ）Ｓn,k,mは前記１つの解析位置ｐ1における周波数軸係数（１／Ｎ）Ｓn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗
【数６】
0006677069.tif 000040
と、
前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,1と時系列データｘn-p1+pmとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする請求項１に記載の定Ｑ変換の成分演算装置。
【請求項５】
時系列データを定比周波数成分へデータ変換する定Ｑ変換の成分演算方法であって、
計算法境界決定手段が、設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、前記ＫMを全周波数成分ｋに各々設定して各ＫM毎の前記合計演算量を求めるステップと、
計算法境界決定手段が、前記各ＫM毎の合計演算量を評価し、該合計演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定するステップと、
計算法選択手段が、定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択するステップと、
定Ｑ変換手段が、時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換ステップと、
を備えたことを特徴とする定Ｑ変換の成分演算方法。
【請求項６】
時間軸係数算出部が、前記定Ｑ変換における時間軸係数を、
【数４】
0006677069.tif 000041
によって算出するステップと、
周波数軸係数算出部が、前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【数５】
0006677069.tif 000042
から（１／Ｎ）Ｓn,kとして算出するステップと、を備え、
前記定Ｑ変換ステップは、
低周波帯域計算部が、前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施するステップと、
高周波帯域計算部が、前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施するステップと、
合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求めるステップと、
を備えていることを特徴とする請求項５に記載の定Ｑ変換の成分演算方法。
【請求項７】
時間軸係数算出部が、
【数４】
0006677069.tif 000043
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）毎に算出して時間軸係数Ｔn,k,mを求め、該Ｔn,k,mを前記定Ｑ変換における時間軸係数とするステップと、
周波数軸係数算出部が、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【数５】
0006677069.tif 000044
から得た（１／Ｎ）Ｓn,kを、前記解析位置ｐm毎に算出して周波数軸係数（１／Ｎ）Ｓn,k,mを求め、該（１／Ｎ）Ｓn,k,mを前記定Ｑ変換における周波数軸係数とするステップと、を備え、
前記定Ｑ変換ステップは、
低周波帯域計算部が、時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、前記解析位置ｐm毎に、前記周波数軸係数算出部で算出された周波数軸係数（１／Ｎ）Ｓn,k,mと前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施するステップと、
高周波帯域計算部が、前記解析位置ｐm毎に、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,mと時系列データｘnとを積和計算して前記第２の計算法を実施するステップと、
合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求めるステップと、
を備えていることを特徴とする請求項５に記載の定Ｑ変換の成分演算方法。
【請求項８】
時間軸係数算出部が、
【数４】
0006677069.tif 000045
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）のうち１つの解析位置ｐ1において算出して時間軸係数Ｔn,k,1を求め、該Ｔn,k,1を前記定Ｑ変換における時間軸係数とするステップと、
周波数軸係数算出部が、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【数５】
0006677069.tif 000046
から得た（１／Ｎ）Ｓn,kを、前記１つの解析位置ｐ1において算出して周波数軸係数（１／Ｎ）Ｓn,k,1を求め、該（１／Ｎ）Ｓn,k,1を前記定Ｑ変換における周波数軸係数とするステップと、を備え、
前記定Ｑ変換ステップは、
低周波帯域計算部が、時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、
前記解析位置ｐm毎に、
前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓn,k,1と、
前記解析位置ｐm毎の周波数軸係数（１／Ｎ）Ｓn,k,mは前記１つの解析位置ｐ1における周波数軸係数（１／Ｎ）Ｓn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗
【数６】
0006677069.tif 000047
と、
前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施するステップと、
高周波帯域計算部が、前記解析位置ｐm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,1と時系列データｘn-p1+pmとを積和計算して前記第２の計算法を実施するステップと、
合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求めるステップと、
を備えていることを特徴とする請求項５に記載の定Ｑ変換の成分演算方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、振動や音響などの時系列データの周波数解析手法の一つである、定比周波数成分へのデータ変換に関し、特に定Ｑ変換の成分演算装置および方法に関する。
【背景技術】
【０００２】
時系列データの周波数解析には不確定性原理があり、周波数解像度を上げようとすればより長時間のデータを解析する必要がある。フーリエ変換は全ての周波数帯域で同一の周波数解像度にて周波数成分に分割して解析するが、音楽解析や騒音解析、振動解析などではオクターブ解析など、周波数を等比の周波数成分に分割して解析することがある。
【０００３】
尚、オクターブ解析は、音響や振動の周波数を、１オクターブ（周波数比２：１の周波数）をＮ分割する等比の周波数帯域（１／Ｎオクターブバンド）に分割して解析する手法である。前記Ｎはビン数ともいい、１，３，６，１２，２４などがよく用いられる。
【０００４】
前記のように周波数を分割して解析を行う場合、高周波帯域の成分は低周波数帯域の成分と同じ周波数解像度は必要ない上に、長時間のデータを解析すると成分がその期間で平均的な値に均されるので、高周波帯域での早い変化を見逃す恐れもある。
【０００５】
このような状況に対応するために定Ｑ変換（Ｃｏｎｓｔａｎｔ  Ｑ  ｔｒａｎｓｆｏｒｍ）という手法が提案されており、例えば非特許文献１〜４などが存在する。
【０００６】
尚、非特許文献３、４はそれぞれ非特許文献１、２の手法の日本語による解説とＰｙｔｈｏｎ言語による検証例を報告しているＷｅｂサイトである。
【０００７】
この定Ｑ変換は、全ての周波数帯域で同じ期間のデータを解析するのではなく、全ての周波数帯域で同じ周期数になるように、周波数毎に参照するデータ数を変えて解析する。この方法では、高周波帯域ほど短い期間のデータで解析するため、低周波帯域で細かい周波数解像度を得ると同時に高周波帯域での早い変化を捉えることができる。
【０００８】
非特許文献１によると、定Ｑ変換は次で定式化される（表記は少し変更している）。
【０００９】
【数１】
0006677069.tif 000002
【００１０】
０以外の値は各ｋに対して、ｎ＝（Ｎ−Ｎk）／２，…，（Ｎ＋Ｎk）／２−１の範囲でのみ取り、この範囲の窓関数の形状として具体的にはハミング窓、ハニング窓、矩形窓など、多くの種類が提案されている
Ｎk：窓幅，Ｎk＝Ｑ（ｆs／ｆk）となるように取る
Ｎ≧Ｎ1≧…≧Ｎk…≧ＮK≧２Ｑとなる
Ｑ：窓の周期数，定Ｑ変換ではこれを全ての周波数で共通にする
非特許文献１、２ではＱ＝（２1/B−１）-1程度に取るようにしている
Ｂはビン数で、１オクターブを分割する数のこと
ｆs：サンプリング周波数
ｆk：第ｋ成分の周波数，ｆk＝（２1/B）k-1ｆmin
ｆmin＝  ｆ1は定Ｑ変換で解析する最小の周波数で、解析する範囲から決める
ｆkは解析する最大の周波数で、ｆk≦ｆs／２となるようにする。
【００１１】
定Ｑ変換は特に低周波帯域の成分で窓幅Ｎkが非常に大きくなることにより演算量が多い。
【００１２】
これに対して非特許文献２では、高速フーリエ変換（ＦＦＴ；Ｆａｓｔ  Ｆｏｕｒｉｅｒ  ｔｒａｎｓｆｏｒｍ）を利用した演算の高速化方法が提案されている。
【００１３】
この方法ではパーセバルの公式により、定Ｑ変換を次の数式（２）のように時間軸での計算から周波数軸での計算に置き換える。
【００１４】
【数２】
0006677069.tif 000003
【００１５】
尚、明細書中の以下の文章では、前記数式（１）、（２）の各左辺の定Ｑ変換の成分を、「Ｘkcq」と表記することもある。
【００１６】
ここでＳn,kは（小さなｋに対して特に）ほとんどの要素が０に近いので、一定値以下の項を０として疎行列表現すると（実質非ゼロ項の積和計算だけになるため）非常に高速に計算できるようになる。
【００１７】
さらに（１／Ｎ）Ｓn,kは、計算するデータ長と解析する周波数範囲が決まっていれば、入力データ列ｘnに寄らず事前に計算できるためこれを予め用意しておけば、入力データの離散フーリエ変換はＦＦＴで実行するとして、定Ｑ変換を入力データに対する１回のＦＦＴと、ＦＦＴ結果ベクトルと予め用意した疎行列との積演算により計算できることになる。
【００１８】
【数３】
0006677069.tif 000004
【００１９】
離散でないフーリエ変換でも同等の公式が成り立ち、これらはパーセバルの定理、プランシュレルの定理など様々な呼び方がされる。
【先行技術文献】
【非特許文献】
【００２０】
【非特許文献１】Ｊｕｄｉｔｈ  Ｃ．Ｂｒｏｗｎ：Ｃａｌｃｕｌａｔｉｏｎ  ｏｆ  ａ  ｃｏｎｓｔａｎｔ  Ｑ  ｓｐｅｃｔｒａｌ  ｔｒａｎｓｆｏｒｍ，Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ．８９（１）：４２５−４３４，１９９１
【非特許文献２】Ｊｕｄｉｔｈ  Ｃ．Ｂｒｏｗｎ  ａｎｄ  Ｍｉｌｌｅｒ  Ｓ．Ｐｕｃｋｅｔｔｅ：Ａｎ  ｅｆｆｉｃｉｅｎｔ  ａｌｇｏｒｉｔｈｍ  ｆｏｒ  ｔｈｅ  ｃａｌｃｕｌａｔｉｏｎ  ｏｆ  ａ  ｃｏｎｓｔａｎｔ  Ｑ  ｔｒａｎｓｆｏｒｍ，Ｊ．Ａｃｏｕｓｔ．Ｓｏｃ．Ａｍ．９２（５）：２６９８−２７０１，１９９２
【非特許文献３】ｙｕｋａｒａ＿１３：[Ｐｙｔｈｏｎ]Ｃｏｎｓｔａｎｔ−Ｑ変換（対数周波数スペクトログラム）、音楽プログラミングの超入門（仮）  ｉｎ  Ｈａｔｅｎａ  Ｂｌｏｇ，２０１３−１２−０１，２０１３、インターネット＜ＵＲＬ：ｈｔｔｐ：／／ｙｕｋａｒａ−１３．ｈａｔｅｎａｂｌｏｇ．ｃｏｍ／ｅｎｔｒｙ／２０１３／１２／０１／２２２７４２＞．[２０１６−０２−２２  アクセス]
【非特許文献４】ｙｕｋａｒａ＿１３：[Ｐｙｔｈｏｎ]  高速なＣｏｎｓｔａｎｔ−Ｑ変換（ｗｉｔｈ  ＦＦＴ）、音楽プログラミングの超入門（仮）  ｉｎ  Ｈａｔｅｎａ  Ｂｌｏｇ，２０１４−０１−０５，２０１４、インターネット＜ＵＲＬ：ｈｔｔｐ：／／ｙｕｋａｒａ−１３．ｈａｔｅｎａｂｌｏｇ．ｃｏｍ／ｅｎｔｒｙ／２０１４／０１／０５／０６２４１４＞．[２０１６−０２−２２  アクセス]
【発明の概要】
【発明が解決しようとする課題】
【００２１】
非特許文献２による高速化のアイデアは周波数軸での係数Ｓn,kが疎行列とみなせることに依拠しているが、実際にはＳn,kが疎行列というのは粗い評価である。確かにその値は一部を除いて小さいが、精度を高めると必ずしも疎にならない。特に高周波帯域でこの傾向が顕著である。実際、非特許文献２、４で係数Ｓn,kが十分小さいとしている閾値はそれぞれ０．１５、０．００５４であり、通常の評価には問題ないが有効数字７桁以上の精度を得るには十分とは言えない。このため、要求精度によってはこの高速化手法がかえって遅くなる可能性がある。
【００２２】
実際に周波数軸での計算でどの程度の演算量を必要とするかを調べるため、非特許文献４に準じた設定により時間軸での計算と周波数軸での計算でそれぞれ計算対象となる非ゼロの項数を評価する。
【００２３】
設定パラメータは、本明細書の記号法で表記してＢ＝２４、Ｑ＝２８、ｆs＝１６０００、ｆmin＝６０、Ｋ＝１６０とし、定Ｑ変換をそれぞれ矩形窓、ハミング窓、ハニング窓で計算する場合において、時間軸での計算に必要な非ゼロ項数と、周波数軸での計算で所定の精度の係数まで計算する際に必要な項数を表１に示す。
【００２４】
【表１】
0006677069.tif 000005
【００２５】
なお、非特許文献１、２に合わせるとＢ＝２４のときＱ＝３４となるが、非特許文献３、４では追加の係数ｆｒａｔｉｏを導入してＱの値を変えているので、その設定値Ｑ＝２８を使っている。また解析する最大周波数６０００ＨｚからＫ＝１６０を導いた。ｆsの値は非特許文献３，４に明記されていなかったので、最大周波数６０００Ｈｚまで解析できるサンプリング周波数１２０００Ｈｚ以上のＷＡＶ音声データとしてｆs＝１６０００Ｈｚを選択した。
【００２６】
表１中の周波数軸右の数値は係数をゼロと評価する係数の大きさの閾値であり、各数値右のパーセンテージは矩形窓の時間軸計算の項数を基準としている。
【００２７】
また、非ゼロの項数と周波数の関係を、時間軸での計算の場合のグラフと周波数軸（閾値ごと）での計算の場合のグラフとして図５〜図８に示す。図５は矩形窓、図６はハミング窓、図７はハニング窓、図８は複素Ｍｏｒｌｅｔウェーブレット変換を用いて計算した場合を各々示している。図５〜図８において、各曲線の下部分の面積が演算量に相当する（ＦＦＴの演算は別である）。
【００２８】
表１から、窓形状により程度が違うが、精度が高いほど（係数をゼロと評価する係数Ｓn,kの大きさの閾値が小さいほど）周波数軸での計算に必要な項数が増えることが分かる。特に矩形窓やハミング窓では精度が高くなると時間軸での計算より周波数軸での計算の方が積和すべき項数が大きくなり、周波数軸での計算で高速にならなくなる。
【００２９】
尚、周波数軸での計算にはＦＦＴ計算と疎行列の非ゼロ係数抽出のオーバーヘッドもあるので、必要項数に差がなければ時間軸計算が有利である。
【００３０】
より精度の高い定Ｑ変換を高速に行うには工夫が必要になる。  
【００３１】
本発明は上記課題を解決するものであり、その目的は、より精度の高い定Ｑ変換を、少ない演算量で高速に行うことができる定Ｑ変換の成分演算装置および方法を提供することにある。
【課題を解決するための手段】
【００３２】
上記課題を解決するための請求項１に記載の定Ｑ変換の成分演算装置は、時系列データを定比周波数成分へデータ変換する定Ｑ変換の成分演算装置であって、
設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、
前記ＫMを全周波数成分ｋに各々設定して各ＫM毎の前記合計演算量を求め、
前記各ＫM毎の合計演算量を評価し、該合計演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段と、
定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段と、
時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段と、
を備えたことを特徴とする。
【００３３】
また、請求項５に記載の定Ｑ変換の成分演算方法は、時系列データを定比周波数成分へデータ変換する定Ｑ変換の成分演算方法であって、
計算法境界決定手段が、設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、前記ＫMを全周波数成分ｋに各々設定して各ＫM毎の前記合計演算量を求めるステップと、
計算法境界決定手段が、前記各ＫM毎の合計演算量を評価し、該合計演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定するステップと、
計算法選択手段が、定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択するステップと、
定Ｑ変換手段が、時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換ステップと、
を備えたことを特徴とする。
【００３４】
また、請求項２に記載の定Ｑ変換の成分演算装置は、請求項１において、前記定Ｑ変換における時間軸係数を、
【００３５】
【数４】
0006677069.tif 000006
【００３６】
によって算出する時間軸係数算出部と、
前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【００３７】
【数５】
0006677069.tif 000007
【００３８】
から（１／Ｎ）Ｓn,kとして算出する周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。
【００３９】
また、請求項６に記載の定Ｑ変換の成分演算方法は、請求項５において、時間軸係数算出部が、前記定Ｑ変換における時間軸係数を、
【００４０】
【数４】
0006677069.tif 000008
【００４１】
によって算出するステップと、
周波数軸係数算出部が、前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【００４２】
【数５】
0006677069.tif 000009
【００４３】
から（１／Ｎ）Ｓn,kとして算出するステップと、を備え、
前記定Ｑ変換ステップは、
低周波帯域計算部が、前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施するステップと、
高周波帯域計算部が、前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施するステップと、
合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求めるステップと、
を備えていることを特徴とする。
【００４４】
定Ｑ変換の成分を演算する場合、低周波帯域では時系列データの離散フーリエ変換と周波数軸係数との積和計算の方が演算量が少なく、高周波帯域では時系列データと時間軸係数との積和計算の方が演算量が少ない。
【００４５】
このため上記構成によれば、低周波帯域、高周波帯域のどちらの場合でも演算量の少ない計算法で定Ｑ変換の成分を演算することができ、演算の高速化が実現される。
【００４６】
また、請求項３に記載の定Ｑ変換の成分演算装置は、請求項１において、
【００４７】
【数４】
0006677069.tif 000010
【００４８】
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）毎に算出して時間軸係数Ｔn,k,mを求め、該Ｔn,k,mを前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【００４９】
【数５】
0006677069.tif 000011
【００５０】
から得た（１／Ｎ）Ｓn,kを、前記解析位置ｐm毎に算出して周波数軸係数（１／Ｎ）Ｓn,k,mを求め、該（１／Ｎ）Ｓn,k,mを前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、前記解析位置ｐm毎に、前記周波数軸係数算出部で算出された周波数軸係数（１／Ｎ）Ｓn,k,mと前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,mと時系列データｘnとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。
【００５１】
また、請求項７に記載の定Ｑ変換の成分演算方法は、請求項５において、時間軸係数算出部が、
【００５２】
【数４】
0006677069.tif 000012
【００５３】
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）毎に算出して時間軸係数Ｔn,k,mを求め、該Ｔn,k,mを前記定Ｑ変換における時間軸係数とするステップと、
周波数軸係数算出部が、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【００５４】
【数５】
0006677069.tif 000013
【００５５】
から得た（１／Ｎ）Ｓ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と各解析位置での周波数軸係数との積和計算でよいため、全体としての計算が高速化される。
【００５７】
また、請求項４に記載の定Ｑ変換の成分演算装置は、請求項１において、
【００５８】
【数４】
0006677069.tif 000014
【００５９】
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）のうち１つの解析位置ｐ1において算出して時間軸係数Ｔn,k,1を求め、該Ｔn,k,1を前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【００６０】
【数５】
0006677069.tif 000015
【００６１】
から得た（１／Ｎ）Ｓn,kを、前記１つの解析位置ｐ1において算出して周波数軸係数（１／Ｎ）Ｓn,k,1を求め、該（１／Ｎ）Ｓn,k,1を前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、  前記解析位置ｐm毎に、
前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓn,k,1と、
前記解析位置ｐm毎の周波数軸係数（１／Ｎ）Ｓn,k,mは前記１つの解析位置ｐ1における周波数軸係数（１／Ｎ）Ｓn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗
【００６２】
【数６】
0006677069.tif 000016
【００６３】
と、
前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,1と時系列データｘn-p1+pmとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。
【００６４】
また、請求項８に記載の定Ｑ変換の成分演算方法は、請求項５において、時間軸係数算出部が、
【００６５】
【数４】
0006677069.tif 000017
【００６６】
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）のうち１つの解析位置ｐ1において算出して時間軸係数Ｔn,k,1を求め、該Ｔn,k,1を前記定Ｑ変換における時間軸係数とするステップと、
周波数軸係数算出部が、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
【００６７】
【数５】
0006677069.tif 000018
【００６８】
から得た（１／Ｎ）Ｓn,kを、前記１つの解析位置ｐ1において算出して周波数軸係数（１／Ｎ）Ｓn,k,1を求め、該（１／Ｎ）Ｓn,k,1を前記定Ｑ変換における周波数軸係数とするステップと、を備え、
前記定Ｑ変換ステップは、
低周波帯域計算部が、時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、
前記解析位置ｐm毎に、
前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓn,k,1と、
前記解析位置ｐm毎の周波数軸係数（１／Ｎ）Ｓn,k,mは前記１つの解析位置ｐ1における周波数軸係数（１／Ｎ）Ｓn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗
【００６９】
【数６】
0006677069.tif 000019
【００７０】
と、
前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施するステップと、
高周波帯域計算部が、前記解析位置ｐm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,1と時系列データｘn-p1+pmとを積和計算して前記第２の計算法を実施するステップと、
合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求めるステップと、
を備えていることを特徴とする。
【００７１】
上記構成によれば、複数の時刻（複数の解析位置ｐm（ｍ＝１，…，Ｍ））で各々定Ｑ変換を行う際に、１つの解析位置ｐ1における時間軸係数および周波数軸係数を利用して各時刻での計算が可能となるので、複数の時刻毎（複数の解析位置毎）に時間軸係数および周波数軸係数を用意する必要はなく、データ量が少なくなってメモリを節約することができる。
【発明の効果】
【００７２】
（１）請求項１〜８に記載の発明によれば、低周波帯域、高周波帯域のどちらの場合でも演算量の少ない計算法で定Ｑ変換の成分を演算することができ、演算の高速化が実現される。
（２）請求項３、７に記載の発明によれば、複数の時刻で各々定Ｑ変換を行う際に、低周波帯域における周波数軸での計算は、時系列データの離散フーリエ変換を、複数の時刻毎に行わなくても、ｘn全体に対して一括して離散フーリエ変換したフーリエ係数Ｘnと各解析位置での周波数軸係数との積和計算でよいため、全体としての計算が高速化される。
（３）請求項４、８に記載の発明によれば、複数の時刻で各々定Ｑ変換を行う際に、１つの解析位置ｐ1における時間軸係数および周波数軸係数を利用して各時刻での計算が可能となるので、複数の時刻毎（複数の解析位置毎）に時間軸係数および周波数軸係数を用意する必要はなく、データ量が少なくなってメモリを節約することができる。
【図面の簡単な説明】
【００７３】
【図１】本発明の実施例による定Ｑ変換の成分演算装置の構成図。
【図２】本発明の実施例２における時間軸・周波数軸の非ゼロ項数を示すグラフ。
【図３】本発明の実施例に１おける時間軸・周波数軸の非ゼロ項数を示すグラフ。
【図４】従来手法と本発明の手法での定Ｑ変換の複数時刻計算の演算量を示すグラフ。
【図５】従来手法の矩形窓における時間軸・周波数軸の項数を示すグラフ。
【図６】従来手法のハミング窓における時間軸・周波数軸の項数を示すグラフ。
【図７】従来手法のハニング窓における時間軸・周波数軸の項数を示すグラフ。
【図８】複素Ｍｏｒｌｅｔウェーブレットの計算項数を示すグラフ。
【発明を実施するための形態】
【００７４】
以下、図面を参照しながら本発明の実施の形態を説明するが、本発明は下記の実施形態例に限定されるものではない。
【００７５】
本実施形態例では、定Ｑ変換を高精度でも演算量を少なく計算するために、定Ｑ変換の成分Ｘkcqを数式（２）の時間軸での計算と周波数軸での計算のどちらの方法で計算するか、予め周波数成分ｋ毎に演算量の小さい方を選択しておくことで実現する。
【００７６】
定Ｑ変換の成分Ｘkcqの計算はおおよそ低周波帯域では周波数軸での計算の演算量が小さく、高周波帯域では時間軸での計算の演算量が小さいから、ｋの境界値ＫMを決めて、ｋ≦ＫMでは周波数軸での計算、ｋ＞ＫMでは時間軸での計算、というようにすれば良い。
【実施例１】
【００７７】
図１は本発明の実施例１による定Ｑ変換の成分演算装置の構成を示している。図１において、１０は予め定Ｑ変換を行うためのパラメータを決定して定Ｑ変換計算の準備をする準備部であり、２０は個別の時系列データに対する定Ｑ変換を計算する個別計算部である。
【００７８】
準備部１０は、
設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、前記ＫMを全周波数成分ｋに各々設定して各ＫM毎の前記合計演算量を求め、前記各ＫM毎の合計演算量を評価し、該演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段の機能と、
定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段の機能とを具備している。
【００７９】
前記計算法境界決定手段および計算法選択手段の各機能を実行する準備部１０は、具体的には、定Ｑ変換のパラメータを決定するパラメータ決定部１１、時間軸係数Ｔn,kを計算する時間軸係数算出部１２、周波数軸係数（１／Ｎ）Ｓn,kを計算する周波数軸係数算出部１３、時間軸での計算と周波数軸での計算との境界のｋ値ＫMを決定する計算法境界決定部１４からなる。
【００８０】
個別計算部２０は、時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段の機能を具備している。
【００８１】
前記定Ｑ変換手段の機能を実行する個別計算部２０は、具体的には、時系列データの長さを整える前処理部２１、高周波帯域の定Ｑ変換を計算する高周波帯域計算部２２、低周波帯域の定Ｑ変換を計算する低周波帯域計算部２３、高周波帯域での結果と低周波帯域での結果をまとめる定Ｑ変換合成部２４からなる。
【００８２】
図１の定Ｑ変換の成分演算装置は、例えばコンピュータにより構成され、通常のコンピュータのハードウェアリソース、例えばＲＯＭ、ＲＡＭ、ＣＰＵ、入力装置、出力装置、通信インターフェース、ハードディスク、記録媒体およびその駆動装置を備えている。
【００８３】
このハードウェアリソースとソフトウェアリソース（ＯＳ、アプリケーションなど）との協働の結果、定Ｑ変換の成分演算装置は、図１に示すように、パラメータ決定部１１、時間軸係数算出部１２、周波数軸係数算出部１３、計算法境界決定部１４、前処理部２１、高周波帯域計算部２２、低周波帯域計算部２３、定Ｑ変換合成部２４を実装する。
【００８４】
また各部での処理結果は、ハードディスク或いはＲＡＭなどの保存手段・記憶手段に格納され、随時利用されるものとする。
【００８５】
次に、図１の各部で実施される内容を詳述する。  
【００８６】
まず、準備部１０のパラメータ決定部１１では、予め定Ｑ変換で解析する時系列データのサンプリング周波数ｆs、解析する周波数範囲（ｆmin，ｆmax）と１オクターブの周波数を区切るビン数Ｂ、解析対象データの周期数Ｑ、窓関数を決定する。このときｆmax≦ｆs／２が必要で、Ｋ＝ｆｌｏｏｒ（Ｂ＊ｌｏｇ2（ｆmax／ｆmin））＋１となる。これを基に計算する入力データの長さをＮ＝Ｎ1＝Ｑ（ｆs／ｆmin）とする。
【００８７】
また周波数軸係数をゼロとみなすべき係数の大きさの閾値Ｔｈもここで決定しておく。
【００８８】
時間軸係数算出部１２では、パラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔn,kを次の定義式に従って計算する。
【００８９】
【数４】
0006677069.tif 000020
【００９０】
周波数軸係数算出部１３では、パラメータ決定部１１で決定したパラメータを基に、周波数軸係数（１／Ｎ）Ｓn,kを、
【００９１】
【数５】
0006677069.tif 000021
【００９２】
の定義式に従って計算する。この際、｜Ｓn,k｜≦Ｔｈとなる項は周波数軸係数（１／Ｎ）Ｓn,k＝０とする。
【００９３】
計算法境界決定部１４では、周波数軸での計算を行う最大のｋ値ＫM（周波数成分閾値）を求める。ＫMは０，…，Ｋの何れかの値を持ち、ＫM＝０は周波数軸での計算を行わずに非特許文献１、２と同様に全て時間軸での計算で求める場合を意味する。
【００９４】
周波数成分閾値ＫMの決定方法としては、ＫM＝０，…，Ｋの各場合について、定Ｑ変換の成分Ｘkcqを求めるのに必要な演算量ＣX（ＫM）を評価してそれが最小となるようにＫMを決定する。
【００９５】
定Ｑ変換の成分Ｘkcqの演算量ＣX（ＫM）（合計演算量）は、長さＮのＦＦＴの演算量ＣFFT（Ｎ）、非ゼロの時間軸係数毎の演算量ＣT、非ゼロの時間軸係数の数ＮT（ｋ）、非ゼロの周波数軸係数毎の演算量ＣF、ｋで非ゼロとみなす周波数軸係数の数ＮF（ｋ）を使っておおよそ
【００９６】
【数７】
0006677069.tif 000022
【００９７】
と表せる。
【００９８】
ここでＣFFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、ＣFFT（Ｎ）＝（１／４）Ｎｌｏｇ2（Ｎ）、Ｃ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の中央部分に一致するようにする。
【０１０５】
次に定Ｑ変換の成分Ｘkcqを、前記計算法境界決定部１４で決定した周波数成分閾値ＫMにより、ｋ≦ＫMについては低周波帯域計算部２３で計算（第１の計算法）し、ｋ＞ＫMについては高周波帯域計算部２２で計算（第２の計算法）する。
【０１０６】
高周波帯域計算部２２は、前記時間軸係数算出部１２で算出された時間軸係数Ｔn,kと時系列データｘnを積和計算して（第２の計算法を実施して）定Ｑ変換の成分Ｘkcqを求める。各ｋについて非ゼロの時間軸係数Ｔn,kは高々Ｎk個なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。
【０１０７】
低周波帯域計算部２３は、まず時系列データｘnをＦＦＴしてフーリエ係数Ｘnを求める。次に前記周波数軸係数算出部１３で算出された周波数軸係数（１／Ｎ）Ｓn,kとフーリエ係数Ｘnを積和計算して（第１の計算法を実施して）定Ｑ変換の成分Ｘkcqを求める。各ｋについて非ゼロの周波数軸係数（１／Ｎ）Ｓn,kは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。
【０１０８】
尚、周波数軸係数は多くの場合、両端（周波数ゼロ近辺とサンプリング周波数近辺）部分では非ゼロ項とゼロ項が入り混じるが中間部分ではゼロ項が続くので、非ゼロ項の範囲を少し広げて、両端部分は全て非ゼロ項として計算しても良い。この場合、図３の細かい点線で示す周波数軸簡易版の特性カーブのように、非ゼロ項が少し増えるがその数はそれほど多くはなく、また非ゼロ判定ロジックが簡略化できるので演算時間は大差なくなる。
【０１０９】
最後に定Ｑ変換合成部２４で、高周波帯域計算部２２と低周波帯域計算部２３とで計算した定Ｑ変換の成分Ｘkcqをまとめて、定Ｑ変換を完了する。
【０１１０】
上記のように本実施例１によれば、従来の高速化技術の問題点である、高い精度の定Ｑ変換を計算しようとすると演算量が膨れ上り通常の時間軸での計算よりもむしろ遅くなるという問題が改善される。
【０１１１】
その効果を示すため、定Ｑ変換の計算方法（時間軸計算、周波数軸計算、提案手法）毎に、必要な演算量を表２のように試算した。試算は、時系列データのパラメータを表１と同じにして、演算量パラメータをＣFFT（Ｎ）＝（１／４）Ｎｌｏｇ2（Ｎ）、ＣT＝１、ＣF＝１として計算した定Ｑ変換の成分Ｘkcqの演算量ＣX（ＫM）を表２に示す。
【０１１２】
【表２】
0006677069.tif 000023
【０１１３】
表２中の周波数軸右の数値は係数をゼロと評価する係数の大きさの閾値であり、各数値右のパーセンテージは矩形窓の時間軸計算の演算量を基準としている。
【０１１４】
表２では、ＦＦＴ演算の演算量が原因で高速化の効果が小さめに出ているが、ハニング窓では閾値１．０ｅ−５でもおおよそ１／５（２２．５％）の演算量になっている。またＦＦＴ演算の演算量が重く周波数軸での計算で効果が出ない場合でも、自動的に時間軸での計算になるので反って遅くなるという問題は発生しない。
【実施例２】
【０１１５】
実施例１は、ある一時刻における定Ｑ変換の各成分値を求める際には効率的であるが、実際の定Ｑ変換では、時系列データの周波数成分の時間変化を捉えるために少しずつ時刻をずらしながら繰り返し定Ｑ成分を計算する。このときずらす時刻は、最短ではサンプリング周波数で、時系列データの情報を余すことなく計算するには最長でも最高周波数での窓幅にしたい。もちろん、時系列データの特徴サンプル抽出などでは、これよりも長い間隔で時刻をずらして定Ｑ変換する場合もあるが、いずれにせよ時系列データ全体に対しては多くの繰り返し計算をすることになり、まだまだ負荷が高い。
【０１１６】
多数の時刻で同時に定Ｑ変換を計算する際には更なる効率化が必要である。  
【０１１７】
そこで本実施例２では、数式（２）の周波数軸での計算で定Ｑ変換を行う場合に、定Ｑ変換を計算する時刻ごとに計算に使用する時系列データの範囲でＤＦＴ（あるいはＦＦＴ）する代わりに、複数時刻での定Ｑ変換に必要な範囲全体で一括してＤＦＴ（あるいはＦＦＴ）を行うことで高速化を図った。
【０１１８】
数式（１）で定義する定Ｑ変換はデータ列ｘn，ｎ＝０，…，Ｎ−１の中央ｘN/2での成分Ｘkcqを計算するように定式化しているが、窓関数の非ゼロ範囲を時刻方向に平行移動することで他の時刻での成分を計算するように調整できる。すなわちＴn,kとｘnとの積和のインデックスをずらすことで実現できる。
【０１１９】
数式（２）による周波数軸での計算に関しては、初めに定Ｑ変換する時刻全体で必要な全データに対して一括してＤＦＴ（あるいはＦＦＴ）を行い、予め各時刻用に用意した（１／Ｎ）Ｓn,kと順次積和することで実現する。
【０１２０】
本実施例２による定Ｑ変換の成分演算装置の具体的な計算の構成は、図１と同様であるが、各部で実行される機能が以下に示すように異なる。
【０１２１】
まず、準備部１０のパラメータ決定部１１では、予め定Ｑ変換で解析する時系列データのサンプリング周波数ｆs、解析する周波数範囲（ｆmin，ｆmax）と１オクターブの周波数を区切るビン数Ｂ、解析対象データの周期数Ｑ、窓関数を決定する。このときｆmax≦ｆs／２が必要で、Ｋ＝ｆｌｏｏｒ（Ｂ＊ｌｏｇ2（ｆmax／ｆmin））＋１となる。
【０１２２】
さらに入力データｘnの全長Ｎと、解析位置ｐm，ｍ＝１，…，Ｍ（時間方向の解析位置；時刻）を決定しておく。このとき、各位置ｐmでの定Ｑ変換はｘpmを中心とする長さＮ1＝Ｑ（ｆs／ｆmin）のデータ列をもとに計算し、（Ｎ1／２）≦ｐ1≦…≦ｐM≦Ｎ−（Ｎ1／２）であるものとする。すなわち、全ての解析位置ｐm，ｍ＝１，…，Ｍで定Ｑ変換の計算に必要なデータがｘ0，…，ｘN-1に含まれるものとする。尚、入力データの全点で定Ｑ変換を計算したいような場合は予め入力データを前後にゼロ拡張して上記条件を満たすように変換しておけばよい。
【０１２３】
また、周波数軸係数をゼロとみなすべき係数の大きさの閾値Ｔｈもここで決定しておく。
【０１２４】
時間軸係数算出部１２ではパラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔn,kを各解析位置ｐmについて次の定義式に従って計算する。
【０１２５】
【数４】
0006677069.tif 000024
【０１２６】
このときｐ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＝０，…，Ｋの各場合について、定Ｑ変換の成分Ｘkcqを求めるのに必要な演算量ＣX（ＫM）を評価してそれが最小となるようにＫMを決定する。尚、この条件は複数のｐmに寄らないので、ｐ1（１つの解析位置）について求めておけば十分である。
【０１３０】
定Ｑ変換の成分Ｘkcqの演算量ＣX（ＫM）（合計演算量）は、長さＮのＦＦＴの演算量ＣFFT（Ｎ）、非ゼロの時間軸係数毎の演算量ＣT、非ゼロの時間軸係数の数ＮT（ｋ）、非ゼロの周波数軸係数毎の演算量ＣF、ｋで非ゼロとみなす周波数軸係数の数ＮF（ｋ）を使っておおよそ、前記実施例１で述べた数式（３）
【０１３１】
【数７】
0006677069.tif 000025
【０１３２】
と表せる。
【０１３３】
ここでＣFFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、ＣFFT（Ｎ）＝（１／４）Ｎｌｏｇ2（Ｎ）、ＣT＝１、ＣF＝１としても良い。
【０１３４】
この演算量ＣX（ＫM）は、時間軸係数および周波数軸係数の、ハニング窓での項数と周波数（周波数成分ｋ；対数周波数）の関係を示す図２のように、非ゼロの時間軸係数の数ＮT（ｋ）がｋに関して指数的に減少し、非ゼロの周波数軸係数の数ＮF（ｋ）がｋに関しておおよそ指数的に増加するため、ＫM＝１，…，Ｋの範囲ではおおよそ下に凸の特性カーブとなる。ただし、ＫM＝０ではＣFFTがなくなるためＣX（０）が小さくなる。結果として演算量ＣX（ＫM）はＫM＝０あるいはＫM＝１，…，Ｋでの最小点の何れかで最小値を取ることになる。
【０１３５】
個別計算部２０では入力された時系列の定Ｑ変換を次のようにして求める。  
【０１３６】
前処理部２１では、入力された時系列データの長さがＮより短ければ前後にゼロを補完して、Ｎより長ければ前後の値を除去して長さがＮの時系列データｘnになるようにする。ゼロを補完する場合には、入力された時系列データの中央部分が時系列データｘnの中央部分に一致するようにする。データを除去する場合は、解析したい範囲に寄るが、単純に後ろのデータを除去するのでも良い。
【０１３７】
次に定Ｑ変換の成分Ｘkcqを、計算法境界決定部１４で決定したＫMにより、ｋ≦ＫMについては低周波帯域計算部２３で、ｋ＞ＫMについては高周波帯域計算部２２で計算する。
【０１３８】
【数８】
0006677069.tif 000026
【０１３９】
尚、前記定Ｑ変換の成分は以下、「Ｘk,mcq」と表記することもある。
【０１４０】
各ｋについて、前記時間軸係数算出部１２で算出された非ゼロの時間軸係数Ｔn,k,mは高々Ｎk個なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。
【０１４１】
低周波帯域計算部２３は、まず時系列データｘnを全体に対して一度ＦＦＴしてフーリエ係数Ｘnを求める。次にｐm毎に、前記周波数軸係数算出部１３で算出された周波数軸係数（１／Ｎ）Ｓn,k,mとフーリエ係数Ｘnを積和計算して定Ｑ変換の成分Ｘk,mcqを求める。
【０１４２】
各ｋについて非ゼロの周波数軸係数（１／Ｎ）Ｓn,k,mは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。
【０１４３】
尚、周波数軸係数は多くの場合、両端（周波数ゼロ近辺とサンプリング周波数近辺）部分では非ゼロ項とゼロ項が入り混じるが中間部分ではゼロ項が続くので、非ゼロ項の範囲を少し広げて、両端部分は全て非ゼロ項として計算しても良い。この場合、図２の点線で示す周波数軸簡易版の特性カーブのように、非ゼロ項が少し増えるがその数はそれほど多くはなく、また非ゼロ判定ロジックが簡略化できるので演算時間は大差なくなる。
【０１４４】
最後に定Ｑ変換合成部２４で、高周波帯域計算部２２と低周波帯域計算部２３とで計算した定Ｑ変換の成分Ｘk,mcqをまとめて、定Ｑ変換を完了する。
【実施例３】
【０１４５】
前記実施例２では、時間軸係数Ｔn,k,m、周波数軸係数（１／Ｎ）Ｓn,k,mは解析位置ｐm毎に用意するので、事前に用意して保管しておくデータ量が大きくなる。本実施例３では、係数データのメモリ使用量を大幅に減らすように以下のような対策を講じた。
【０１４６】
【数９】
0006677069.tif 000027
【０１４７】
尚、数式（４）中の
【０１４８】
【数６】
0006677069.tif 000028
【０１４９】
は、複素指数関数の時刻および周波数毎に異なるべき乗を表している。
【０１５０】
本実施例３による定Ｑ変換の成分演算装置の具体的な構成は図１と同様であり、以下に図１と異なる部分のみを説明する。
【０１５１】
時間軸係数算出部１２では、パラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔn,k,1を定義式に従って求める。すなわち、
【０１５２】
【数４】
0006677069.tif 000029
【０１５３】
を、予め設定した解析位置ｐmのうち１つの解析位置（例えば１番目の解析位置）ｐ1において算出し、時間軸係数Ｔn,k,1を求める。
【０１５４】
周波数軸係数算出部１３では、パラメータ決定部１１で決定したパラメータを基に、周波数軸係数（１／Ｎ）Ｓn,k,1を定義式に従って計算する。
【０１５５】
すなわち、
【０１５６】
【数５】
0006677069.tif 000030
【０１５７】
から得た（１／Ｎ）Ｓn,kを、前記１つの解析位置ｐ1において算出して周波数軸係数（１／Ｎ）Ｓn,k,1を求める。
【０１５８】
この際、｜Ｓn,k,1｜≦Ｔｈとなる項は周波数軸係数（１／Ｎ）Ｓn,k,1＝０とする。
【０１５９】
計算法境界決定部１４では、周波数軸計算の際に複素べき乗の積和が追加されるために非ゼロの周波数軸係数毎の演算量ＣFを実施例２の場合より大きく、例えば５（複素べき乗の演算に＋３、複素べき乗の積和に＋１）とする。
【０１６０】
高周波帯域計算部２２は、複数の解析位置ｐm毎に範囲をずらして、前記時間軸係数算出部１２で算出された時間軸係数Ｔn,k,1と時系列データｘn-p1+pmを積和計算して定Ｑ変換の成分Ｘk,mcqを求める。
【０１６１】
低周波帯域計算部２３は、まず時系列データｘnを全体に対して一度ＦＦＴしてフーリエ係数Ｘnを求める。
【０１６２】
【数１０】
0006677069.tif 000031
【０１６３】
以上、定Ｑ変換の高速計算方式を説明したが、他の定比周波数成分へのデータ変換方法、例えばウェーブレット変換も計算の定式化はほとんど同じ形式にすることが可能であり、本発明と同様の計算方法を適用することは可能である。しかし周波数軸係数が疎行列に近くなるか否かは使用するウェーブレットに依存するため、必ずしも高速計算方式になるものではない。ただし例えば複素Ｍｏｒｌｅｔ  Ｗａｖｅｌｅｔに関しては、少なくとも一つの時刻について計算する場合は高速化が期待できる。
【０１６４】
前記実施例２、３の手法によれば、定Ｑ変換を高速に実行するにあたって、先行の高速化技術（実施例１）が時間軸での計算でも周波数軸での計算でもそれらを併用する計算方法でも、一時刻毎に定Ｑ変換を計算するのに対して、複数時刻での定Ｑ変換を同時に計算する際に、ＤＦＴ（あるいはＦＦＴ）を一時刻の計算範囲毎ではなく全体に対して一度だけ実施することで、全体としての計算が高速化される。
【０１６５】
実施例３では、時間軸と周波数軸との係数データを計算時刻毎に持つ代わりに、それらを複素数のべき乗計算に置き換えて（数式（４））、演算量は少し増える代わりに係数データのメモリ使用量を大幅に減らしている。
【０１６６】
効果を示すために、定Ｑ変換を時間軸計算で時刻毎に逐次計算した場合、実施例１の手法で逐次計算した場合、実施例２の手法で一括計算した場合、実施例３の方法でメモリを節約する場合の各方法にて演算量を評価した。尚、評価における設定パラメータは以下の通りである。
【０１６７】
【表３】
0006677069.tif 000032
【０１６８】
また、演算量のパラメータは実機での計測から以下の通りである。
【０１６９】
【表４】
0006677069.tif 000033
【０１７０】
これらを基に、１から１０００個の時刻で定Ｑ変換を計算した際に掛かる演算量を試算してグラフにすると図４のようになる。
【０１７１】
図４より、実施例１は時間軸計算より高速であり、少なくとも１０時刻以上を一括で計算する場合には実施例２、３がさらに高速になることが分かる。
【０１７２】
おおむね１００時刻以上の一括計算では、実施例２で実施例１の１０倍、実施例３で実施例１の５倍の高速化が見込める。
【０１７３】
尚、時間軸係数Ｔn,k,1の非ゼロ項数は実施例１でも実施例２、３でも変わりないが、周波数軸係数（１／Ｎ）Ｓn,k,1の非ゼロ項数については、実施例１と実施例２、３でＦＦＴするデータ数（＝行列（１／Ｎ）Ｓn,k,1の次元）が異なるために、実施例１の場合より実施例２、３の場合の方が倍以上多くなる。また一回のＦＦＴに掛かる時間も同様に大幅に大きくなる。これらはもちろん図４に反映している。
【０１７４】
それでも実施例２、３の手法が高速になるのはＦＦＴの実行回数が大きく異なるためである。
【符号の説明】
【０１７５】
１０…準備部
１１…パラメータ決定部
１２…時間軸係数算出部
１３…周波数軸係数算出部
１４…計算法境界決定部
２０…個別計算部
２１…前処理部
２２…高周波帯域計算部
２３…低周波帯域計算部
２４…定Ｑ変換合成部
【図１】
0006677069.tif 000048
【図２】
0006677069.tif 000049
【図３】
0006677069.tif 000050
【図４】
0006677069.tif 000051
【図５】
0006677069.tif 000052
【図６】
0006677069.tif 000053
【図７】
0006677069.tif 000054
【図８】
0006677069.tif 000055