図面 (/)

技術 定Q変換の成分演算装置および定Q変換の成分演算方法

出願人 株式会社明電舎
発明者 林孝則外山達斎
出願日 2016年4月28日 (4年7ヶ月経過) 出願番号 2016-091694
公開日 2017年11月2日 (3年1ヶ月経過) 公開番号 2017-198619
状態 特許登録済
技術分野 機械的振動・音波の測定 複合演算
主要キーワード 記号法 合計演算 計算時刻 解析位置 周波数解析手法 周波数ゼロ 窓形状 矩形窓
関連する未来課題
重要な関連分野

この項目の情報は公開日時点(2017年11月2日)のものです。
また、この項目は機械的に抽出しているため、正しく解析できていない場合があります

図面 (9)

課題

より精度の高い定Q変換を、少ない演算量で高速に行うことができるようにする。

解決手段

設定した周波数成分KM以下では時系列データの離散フーリエ変換周波数軸係数との積和計算により求めた演算量と、前記KMを超える帯域では時系列データと時間軸係数との積和計算により求めた演算量とを含む合計演算量を算出し、前記KMを全周波数成分kに各々設定して各KM毎の前記合計演算量を求め、該KM毎の合計演算量が最小となるKMを周波数成分閾値として決定する計算法境界決定部14と、前記kがKM以下の場合は低周波帯域計算部23を選択して周波数軸係数による計算を行い、kがKMを超える場合は高周波帯域計算部22を選択して時間軸係数による計算を行い、両方の計算結果を定Q変換合成部24で合成する。

概要

背景

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

尚、オクターブ解析は、音響や振動の周波数を、1オクターブ(周波数比2:1の周波数)をN分割する等比の周波数帯域(1/Nオクターブバンド)に分割して解析する手法である。前記Nはビン数ともいい、1,3,6,12,24などがよく用いられる。

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

このような状況に対応するために定Q変換(Constant Q transform)という手法が提案されており、例えば非特許文献1〜4などが存在する。

尚、非特許文献3、4はそれぞれ非特許文献1、2の手法の日本語による解説とPython言語による検証例を報告しているWebサイトである。

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

非特許文献1によると、定Q変換は次で定式化される(表記は少し変更している)。

0以外の値は各kに対して、n=(N−Nk)/2,…,(N+Nk)/2−1の範囲でのみ取り、この範囲の窓関数の形状として具体的にはハミング窓ハニング窓、矩形窓など、多くの種類が提案されている
Nk:窓幅,Nk=Q(fs/fk)となるように取る
N≧N1≧…≧Nk…≧NK≧2Qとなる
Q:窓の周期数,定Q変換ではこれを全ての周波数で共通にする
非特許文献1、2ではQ=(21/B−1)-1程度に取るようにしている
Bはビン数で、1オクターブを分割する数のこと
fs:サンプリング周波数
fk:第k成分の周波数,fk=(21/B)k-1fmin
fmin= f1は定Q変換で解析する最小の周波数で、解析する範囲から決める
fkは解析する最大の周波数で、fk≦fs/2となるようにする。

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

これに対して非特許文献2では、高速フーリエ変換FFT;Fast Fourier transform)を利用した演算の高速化方法が提案されている。

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

尚、明細書中の以下の文章では、前記数式(1)、(2)の各左辺の定Q変換の成分を、「Xkcq」と表記することもある。

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

さらに(1/N)Sn,kは、計算するデータ長と解析する周波数範囲が決まっていれば、入力データ列xnに寄らず事前に計算できるためこれを予め用意しておけば、入力データの離散フーリエ変換はFFTで実行するとして、定Q変換を入力データに対する1回のFFTと、FFT結果ベクトルと予め用意した疎行列との積演算により計算できることになる。

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

概要

より精度の高い定Q変換を、少ない演算量で高速に行うことができるようにする。設定した周波数成分KM以下では時系列データの離散フーリエ変換と周波数軸係数との積和計算により求めた演算量と、前記KMを超える帯域では時系列データと時間軸係数との積和計算により求めた演算量とを含む合計演算量を算出し、前記KMを全周波数成分kに各々設定して各KM毎の前記合計演算量を求め、該KM毎の合計演算量が最小となるKMを周波数成分閾値として決定する計算法境界決定部14と、前記kがKM以下の場合は低周波帯域計算部23を選択して周波数軸係数による計算を行い、kがKMを超える場合は高周波帯域計算部22を選択して時間軸係数による計算を行い、両方の計算結果を定Q変換合成部24で合成する。

目的

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

効果

実績

技術文献被引用数
0件
牽制数
0件

この技術が所属する分野

ライセンス契約や譲渡などの可能性がある特許掲載中! 開放特許随時追加・更新中 詳しくはこちら

請求項1

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

請求項2

前記定Q変換における時間軸係数を、によって算出する時間軸係数算出部と、前記定Q変換における周波数軸係数を、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換であるから(1/N)Sn,kとして算出する周波数軸係数算出部と、を備え、前記定Q変換手段は、前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第1の計算法を実施する低周波帯域計算部と、前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第2の計算法を実施する高周波帯域計算部と、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求める合成部と、を備えていることを特徴とする請求項1に記載の定Q変換の成分演算装置。

請求項3

を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)毎に算出して時間軸係数Tn,k,mを求め、該Tn,k,mを前記定Q変換における時間軸係数とする時間軸係数算出部と、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換であるから得た(1/N)Sn,kを、前記解析位置pm毎に算出して周波数軸係数(1/N)Sn,k,mを求め、該(1/N)Sn,k,mを前記定Q変換における周波数軸係数とする周波数軸係数算出部と、を備え、前記定Q変換手段は、時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、前記解析位置pm毎に、前記周波数軸係数算出部で算出された周波数軸係数(1/N)Sn,k,mと前記フーリエ係数Xnとを積和計算して前記第1の計算法を実施する低周波帯域計算部と、前記解析位置pm毎に、前記時間軸係数算出部で算出された時間軸係数Tn,k,mと時系列データxnとを積和計算して前記第2の計算法を実施する高周波帯域計算部と、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求める合成部と、を備えていることを特徴とする請求項1に記載の定Q変換の成分演算装置。

請求項4

を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)のうち1つの解析位置p1において算出して時間軸係数Tn,k,1を求め、該Tn,k,1を前記定Q変換における時間軸係数とする時間軸係数算出部と、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換であるから得た(1/N)Sn,kを、前記1つの解析位置p1において算出して周波数軸係数(1/N)Sn,k,1を求め、該(1/N)Sn,k,1を前記定Q変換における周波数軸係数とする周波数軸係数算出部と、を備え、前記定Q変換手段は、時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、前記解析位置pm毎に、前記周波数軸係数算出部で求められた周波数軸係数(1/N)Sn,k,1と、前記解析位置pm毎の周波数軸係数(1/N)Sn,k,mは前記1つの解析位置p1における周波数軸係数(1/N)Sn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗と、前記フーリエ係数Xnとを積和計算して前記第1の計算法を実施する低周波帯域計算部と、前記解析位置pm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Tn,k,1と時系列データxn-p1+pmとを積和計算して前記第2の計算法を実施する高周波帯域計算部と、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求める合成部と、を備えていることを特徴とする請求項1に記載の定Q変換の成分演算装置。

請求項5

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

請求項6

時間軸係数算出部が、前記定Q変換における時間軸係数を、によって算出するステップと、周波数軸係数算出部が、前記定Q変換における周波数軸係数を、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換であるから(1/N)Sn,kとして算出するステップと、を備え、前記定Q変換ステップは、低周波帯域計算部が、前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第1の計算法を実施するステップと、高周波帯域計算部が、前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第2の計算法を実施するステップと、合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求めるステップと、を備えていることを特徴とする請求項5に記載の定Q変換の成分演算方法。

請求項7

時間軸係数算出部が、を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)毎に算出して時間軸係数Tn,k,mを求め、該Tn,k,mを前記定Q変換における時間軸係数とするステップと、周波数軸係数算出部が、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換であるから得た(1/N)Sn,kを、前記解析位置pm毎に算出して周波数軸係数(1/N)Sn,k,mを求め、該(1/N)Sn,k,mを前記定Q変換における周波数軸係数とするステップと、を備え、前記定Q変換ステップは、低周波帯域計算部が、時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、前記解析位置pm毎に、前記周波数軸係数算出部で算出された周波数軸係数(1/N)Sn,k,mと前記フーリエ係数Xnとを積和計算して前記第1の計算法を実施するステップと、高周波帯域計算部が、前記解析位置pm毎に、前記時間軸係数算出部で算出された時間軸係数Tn,k,mと時系列データxnとを積和計算して前記第2の計算法を実施するステップと、合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求めるステップと、を備えていることを特徴とする請求項5に記載の定Q変換の成分演算方法。

請求項8

時間軸係数算出部が、を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)のうち1つの解析位置p1において算出して時間軸係数Tn,k,1を求め、該Tn,k,1を前記定Q変換における時間軸係数とするステップと、周波数軸係数算出部が、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換であるから得た(1/N)Sn,kを、前記1つの解析位置p1において算出して周波数軸係数(1/N)Sn,k,1を求め、該(1/N)Sn,k,1を前記定Q変換における周波数軸係数とするステップと、を備え、前記定Q変換ステップは、低周波帯域計算部が、時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、前記解析位置pm毎に、前記周波数軸係数算出部で求められた周波数軸係数(1/N)Sn,k,1と、前記解析位置pm毎の周波数軸係数(1/N)Sn,k,mは前記1つの解析位置p1における周波数軸係数(1/N)Sn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗と、前記フーリエ係数Xnとを積和計算して前記第1の計算法を実施するステップと、高周波帯域計算部が、前記解析位置pm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Tn,k,1と時系列データxn-p1+pmとを積和計算して前記第2の計算法を実施するステップと、合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求めるステップと、を備えていることを特徴とする請求項5に記載の定Q変換の成分演算方法。

技術分野

0001

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

背景技術

0002

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

0003

尚、オクターブ解析は、音響や振動の周波数を、1オクターブ(周波数比2:1の周波数)をN分割する等比の周波数帯域(1/Nオクターブバンド)に分割して解析する手法である。前記Nはビン数ともいい、1,3,6,12,24などがよく用いられる。

0004

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

0005

このような状況に対応するために定Q変換(Constant Q transform)という手法が提案されており、例えば非特許文献1〜4などが存在する。

0006

尚、非特許文献3、4はそれぞれ非特許文献1、2の手法の日本語による解説とPython言語による検証例を報告しているWebサイトである。

0007

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

0008

非特許文献1によると、定Q変換は次で定式化される(表記は少し変更している)。

0009

0010

0以外の値は各kに対して、n=(N−Nk)/2,…,(N+Nk)/2−1の範囲でのみ取り、この範囲の窓関数の形状として具体的にはハミング窓ハニング窓、矩形窓など、多くの種類が提案されている
Nk:窓幅,Nk=Q(fs/fk)となるように取る
N≧N1≧…≧Nk…≧NK≧2Qとなる
Q:窓の周期数,定Q変換ではこれを全ての周波数で共通にする
非特許文献1、2ではQ=(21/B−1)-1程度に取るようにしている
Bはビン数で、1オクターブを分割する数のこと
fs:サンプリング周波数
fk:第k成分の周波数,fk=(21/B)k-1fmin
fmin= f1は定Q変換で解析する最小の周波数で、解析する範囲から決める
fkは解析する最大の周波数で、fk≦fs/2となるようにする。

0011

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

0012

これに対して非特許文献2では、高速フーリエ変換FFT;Fast Fourier transform)を利用した演算の高速化方法が提案されている。

0013

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

0014

0015

尚、明細書中の以下の文章では、前記数式(1)、(2)の各左辺の定Q変換の成分を、「Xkcq」と表記することもある。

0016

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

0017

さらに(1/N)Sn,kは、計算するデータ長と解析する周波数範囲が決まっていれば、入力データ列xnに寄らず事前に計算できるためこれを予め用意しておけば、入力データの離散フーリエ変換はFFTで実行するとして、定Q変換を入力データに対する1回のFFTと、FFT結果ベクトルと予め用意した疎行列との積演算により計算できることになる。

0018

0019

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

先行技術

0020

Judith C.Brown:Calculation of a constant Q spectral transform,J.Acoust.Soc.Am.89(1):425−434,1991
Judith C.Brown and Miller S.Puckette:An efficient algorithm for the calculation of a constant Q transform,J.Acoust.Soc.Am.92(5):2698−2701,1992
yukara_13:[Python]Constant−Q変換(対数周波数スペクトログラム)、音楽プログラミングの超入門(仮) in Hatena Blog,2013−12−01,2013、インターネット<URL:http://yukara−13.hatenablog.com/entry/2013/12/01/222742>.[2016−02−22アクセス]
yukara_13:[Python]高速なConstant−Q変換(withFFT)、音楽プログラミングの超入門(仮) in Hatena Blog,2014−01−05,2014、インターネット<URL:http://yukara−13.hatenablog.com/entry/2014/01/05/062414>.[2016−02−22 アクセス]

発明が解決しようとする課題

0021

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

0022

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

0023

設定パラメータは、本明細書の記号法で表記してB=24、Q=28、fs=16000、fmin=60、K=160とし、定Q変換をそれぞれ矩形窓、ハミング窓、ハニング窓で計算する場合において、時間軸での計算に必要な非ゼロ項数と、周波数軸での計算で所定の精度の係数まで計算する際に必要な項数を表1に示す。

0024

0025

なお、非特許文献1、2に合わせるとB=24のときQ=34となるが、非特許文献3、4では追加の係数fratioを導入してQの値を変えているので、その設定値Q=28を使っている。また解析する最大周波数6000HzからK=160を導いた。fsの値は非特許文献3,4に明記されていなかったので、最大周波数6000Hzまで解析できるサンプリング周波数12000Hz以上のWAV音声データとしてfs=16000Hzを選択した。

0026

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

0027

また、非ゼロの項数と周波数の関係を、時間軸での計算の場合のグラフと周波数軸(閾値ごと)での計算の場合のグラフとして図5図8に示す。図5は矩形窓、図6はハミング窓、図7はハニング窓、図8複素Morletウェーブレット変換を用いて計算した場合を各々示している。図5図8において、各曲線の下部分の面積が演算量に相当する(FFTの演算は別である)。

0028

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

0029

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

0030

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

0031

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

課題を解決するための手段

0032

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

0033

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

0034

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

0035

0036

によって算出する時間軸係数算出部と、
前記定Q変換における周波数軸係数を、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換である

0037

0038

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

0039

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

0040

0041

によって算出するステップと、
周波数軸係数算出部が、前記定Q変換における周波数軸係数を、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換である

0042

0043

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

0044

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

0045

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

0046

また、請求項3に記載の定Q変換の成分演算装置は、請求項1において、

0047

0048

を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)毎に算出して時間軸係数Tn,k,mを求め、該Tn,k,mを前記定Q変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換である

0049

0050

から得た(1/N)Sn,kを、前記解析位置pm毎に算出して周波数軸係数(1/N)Sn,k,mを求め、該(1/N)Sn,k,mを前記定Q変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Q変換手段は、
時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、前記解析位置pm毎に、前記周波数軸係数算出部で算出された周波数軸係数(1/N)Sn,k,mと前記フーリエ係数Xnとを積和計算して前記第1の計算法を実施する低周波帯域計算部と、
前記解析位置pm毎に、前記時間軸係数算出部で算出された時間軸係数Tn,k,mと時系列データxnとを積和計算して前記第2の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求める合成部と、
を備えていることを特徴とする。

0051

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

0052

0053

を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)毎に算出して時間軸係数Tn,k,mを求め、該Tn,k,mを前記定Q変換における時間軸係数とするステップと、
周波数軸係数算出部が、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換である

0054

0055

から得た(1/N)Sn,kを、前記解析位置pm毎に算出して周波数軸係数(1/N)Sn,k,mを求め、該(1/N)Sn,k,mを前記定Q変換における周波数軸係数とするステップと、を備え、
前記定Q変換ステップは、
低周波帯域計算部が、時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、前記解析位置pm毎に、前記周波数軸係数算出部で算出された周波数軸係数(1/N)Sn,k,mと前記フーリエ係数Xnとを積和計算して前記第1の計算法を実施するステップと、
高周波帯域計算部が、前記解析位置pm毎に、前記時間軸係数算出部で算出された時間軸係数Tn,k,mと時系列データxnとを積和計算して前記第2の計算法を実施するステップと、
合成部が、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求めるステップと、
を備えていることを特徴とする。

0056

上記構成によれば、複数の時刻(複数の解析位置pm(m=1,…,M))で各々定Q変換を行う際に、時間軸係数および周波数軸係数は複数の時刻毎(解析位置pm毎)に用意するが、低周波帯域における周波数軸での計算は、時系列データの離散フーリエ変換を、複数の時刻毎に行わなくても、xn全体に対して一括して離散フーリエ変換したフーリエ係数Xnと各解析位置での周波数軸係数との積和計算でよいため、全体としての計算が高速化される。

0057

また、請求項4に記載の定Q変換の成分演算装置は、請求項1において、

0058

0059

を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)のうち1つの解析位置p1において算出して時間軸係数Tn,k,1を求め、該Tn,k,1を前記定Q変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換である

0060

0061

から得た(1/N)Sn,kを、前記1つの解析位置p1において算出して周波数軸係数(1/N)Sn,k,1を求め、該(1/N)Sn,k,1を前記定Q変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Q変換手段は、
時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、 前記解析位置pm毎に、
前記周波数軸係数算出部で求められた周波数軸係数(1/N)Sn,k,1と、
前記解析位置pm毎の周波数軸係数(1/N)Sn,k,mは前記1つの解析位置p1における周波数軸係数(1/N)Sn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗

0062

0063

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

0064

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

0065

0066

を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)のうち1つの解析位置p1において算出して時間軸係数Tn,k,1を求め、該Tn,k,1を前記定Q変換における時間軸係数とするステップと、
周波数軸係数算出部が、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換である

0067

0068

から得た(1/N)Sn,kを、前記1つの解析位置p1において算出して周波数軸係数(1/N)Sn,k,1を求め、該(1/N)Sn,k,1を前記定Q変換における周波数軸係数とするステップと、を備え、
前記定Q変換ステップは、
低周波帯域計算部が、時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、
前記解析位置pm毎に、
前記周波数軸係数算出部で求められた周波数軸係数(1/N)Sn,k,1と、
前記解析位置pm毎の周波数軸係数(1/N)Sn,k,mは前記1つの解析位置p1における周波数軸係数(1/N)Sn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗

0069

0070

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

0071

上記構成によれば、複数の時刻(複数の解析位置pm(m=1,…,M))で各々定Q変換を行う際に、1つの解析位置p1における時間軸係数および周波数軸係数を利用して各時刻での計算が可能となるので、複数の時刻毎(複数の解析位置毎)に時間軸係数および周波数軸係数を用意する必要はなく、データ量が少なくなってメモリを節約することができる。

発明の効果

0072

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

図面の簡単な説明

0073

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

0074

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

0075

本実施形態例では、定Q変換を高精度でも演算量を少なく計算するために、定Q変換の成分Xkcqを数式(2)の時間軸での計算と周波数軸での計算のどちらの方法で計算するか、予め周波数成分k毎に演算量の小さい方を選択しておくことで実現する。

0076

定Q変換の成分Xkcqの計算はおおよそ低周波帯域では周波数軸での計算の演算量が小さく、高周波帯域では時間軸での計算の演算量が小さいから、kの境界値KMを決めて、k≦KMでは周波数軸での計算、k>KMでは時間軸での計算、というようにすれば良い。

0077

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

0078

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

0079

前記計算法境界決定手段および計算法選択手段の各機能を実行する準備部10は、具体的には、定Q変換のパラメータを決定するパラメータ決定部11、時間軸係数Tn,kを計算する時間軸係数算出部12、周波数軸係数(1/N)Sn,kを計算する周波数軸係数算出部13、時間軸での計算と周波数軸での計算との境界のk値KMを決定する計算法境界決定部14からなる。

0080

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

0081

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

0082

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

0083

このハードウェアリソースとソフトウェアリソース(OS、アプリケーションなど)との協働の結果、定Q変換の成分演算装置は、図1に示すように、パラメータ決定部11、時間軸係数算出部12、周波数軸係数算出部13、計算法境界決定部14、前処理部21、高周波帯域計算部22、低周波帯域計算部23、定Q変換合成部24を実装する。

0084

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

0085

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

0086

まず、準備部10のパラメータ決定部11では、予め定Q変換で解析する時系列データのサンプリング周波数fs、解析する周波数範囲(fmin,fmax)と1オクターブの周波数を区切るビン数B、解析対象データの周期数Q、窓関数を決定する。このときfmax≦fs/2が必要で、K=floor(B*log2(fmax/fmin))+1となる。これを基に計算する入力データの長さをN=N1=Q(fs/fmin)とする。

0087

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

0088

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

0089

0090

周波数軸係数算出部13では、パラメータ決定部11で決定したパラメータを基に、周波数軸係数(1/N)Sn,kを、

0091

0092

の定義式に従って計算する。この際、|Sn,k|≦Thとなる項は周波数軸係数(1/N)Sn,k=0とする。

0093

計算法境界決定部14では、周波数軸での計算を行う最大のk値KM(周波数成分閾値)を求める。KMは0,…,Kの何れかの値を持ち、KM=0は周波数軸での計算を行わずに非特許文献1、2と同様に全て時間軸での計算で求める場合を意味する。

0094

周波数成分閾値KMの決定方法としては、KM=0,…,Kの各場合について、定Q変換の成分Xkcqを求めるのに必要な演算量CX(KM)を評価してそれが最小となるようにKMを決定する。

0095

定Q変換の成分Xkcqの演算量CX(KM)(合計演算量)は、長さNのFFTの演算量CFFT(N)、非ゼロの時間軸係数毎の演算量CT、非ゼロの時間軸係数の数NT(k)、非ゼロの周波数軸係数毎の演算量CF、kで非ゼロとみなす周波数軸係数の数NF(k)を使っておおよそ

0096

0097

と表せる。

0098

ここでCFFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、CFFT(N)=(1/4)Nlog2(N)、CT=1、CF=1としても良い。

0099

上記演算量CX(KM)を示す数式(3)において、右辺第1項は、KMが非ゼロのときCFFTとし、そうでないときは0とすることを表し、右辺第2項は、周波数成分kが周波数成分閾値KM以下の低周波帯域における周波数軸係数毎の演算量を表し、右辺第3項は、周波数成分kが周波数成分閾値KMを超える高周波帯域における時間軸係数毎の演算量を表している。

0100

このためKMが非ゼロのときは、数式(3)の右辺第1項のFFTの演算量および右辺第2項の演算量(周波数軸係数による演算量)からなる第1の演算量と、数式(3)の右辺第3項の時間軸係数による第2の演算量との合計が演算量CX(KM)となる。

0101

またKMがゼロのときは、数式(3)の右辺第1項および右辺第2項がともに0となって右辺第3項の時間軸係数による第2の演算量のみが演算量CX(KM)となる。

0102

この演算量CX(KM)は、時間軸係数および周波数軸係数の、ハニング窓での項数と周波数(周波数成分k;対数周波数)の関係を示す図3のように、非ゼロの時間軸係数の数NT(k)がkに関して指数的に減少し、非ゼロの周波数軸係数の数NF(k)がkに関しておおよそ指数的に増加するため、KM=1,…,Kの範囲ではおおよそ下に凸の特性カーブとなる。ただし、KM=0ではCFFTがなくなるためCX(0)が小さくなる。結果として演算量CX(KM)はKM=0あるいはKM=1,…,Kでの最小点の何れかで最小値を取ることになる。

0103

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

0104

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

0105

次に定Q変換の成分Xkcqを、前記計算法境界決定部14で決定した周波数成分閾値KMにより、k≦KMについては低周波帯域計算部23で計算(第1の計算法)し、k>KMについては高周波帯域計算部22で計算(第2の計算法)する。

0106

高周波帯域計算部22は、前記時間軸係数算出部12で算出された時間軸係数Tn,kと時系列データxnを積和計算して(第2の計算法を実施して)定Q変換の成分Xkcqを求める。各kについて非ゼロの時間軸係数Tn,kは高々Nk個なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。

0107

低周波帯域計算部23は、まず時系列データxnをFFTしてフーリエ係数Xnを求める。次に前記周波数軸係数算出部13で算出された周波数軸係数(1/N)Sn,kとフーリエ係数Xnを積和計算して(第1の計算法を実施して)定Q変換の成分Xkcqを求める。各kについて非ゼロの周波数軸係数(1/N)Sn,kは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。

0108

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

0109

最後に定Q変換合成部24で、高周波帯域計算部22と低周波帯域計算部23とで計算した定Q変換の成分Xkcqをまとめて、定Q変換を完了する。

0110

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

0111

その効果を示すため、定Q変換の計算方法(時間軸計算、周波数軸計算、提案手法)毎に、必要な演算量を表2のように試算した。試算は、時系列データのパラメータを表1と同じにして、演算量パラメータをCFFT(N)=(1/4)Nlog2(N)、CT=1、CF=1として計算した定Q変換の成分Xkcqの演算量CX(KM)を表2に示す。

0112

0113

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

0114

表2では、FFT演算の演算量が原因で高速化の効果が小さめに出ているが、ハニング窓では閾値1.0e−5でもおおよそ1/5(22.5%)の演算量になっている。またFFT演算の演算量が重く周波数軸での計算で効果が出ない場合でも、自動的に時間軸での計算になるので反って遅くなるという問題は発生しない。

0115

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

0116

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

0117

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

0118

数式(1)で定義する定Q変換はデータ列xn,n=0,…,N−1の中央xN/2での成分Xkcqを計算するように定式化しているが、窓関数の非ゼロ範囲を時刻方向に平行移動することで他の時刻での成分を計算するように調整できる。すなわちTn,kとxnとの積和のインデックスをずらすことで実現できる。

0119

数式(2)による周波数軸での計算に関しては、初めに定Q変換する時刻全体で必要な全データに対して一括してDFT(あるいはFFT)を行い、予め各時刻用に用意した(1/N)Sn,kと順次積和することで実現する。

0120

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

0121

まず、準備部10のパラメータ決定部11では、予め定Q変換で解析する時系列データのサンプリング周波数fs、解析する周波数範囲(fmin,fmax)と1オクターブの周波数を区切るビン数B、解析対象データの周期数Q、窓関数を決定する。このときfmax≦fs/2が必要で、K=floor(B*log2(fmax/fmin))+1となる。

0122

さらに入力データxnの全長Nと、解析位置pm,m=1,…,M(時間方向の解析位置;時刻)を決定しておく。このとき、各位置pmでの定Q変換はxpmを中心とする長さN1=Q(fs/fmin)のデータ列をもとに計算し、(N1/2)≦p1≦…≦pM≦N−(N1/2)であるものとする。すなわち、全ての解析位置pm,m=1,…,Mで定Q変換の計算に必要なデータがx0,…,xN-1に含まれるものとする。尚、入力データの全点で定Q変換を計算したいような場合は予め入力データを前後にゼロ拡張して上記条件を満たすように変換しておけばよい。

0123

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

0124

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

0125

0126

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

0127

周波数軸係数算出部13ではパラメータ決定部11で決定したパラメータを基に、各解析位置pmについて周波数軸係数(1/N)Sn,k,mを定義式に従って計算する。この際、|Sn,k,m|≦Thとなる項は周波数軸係数(1/N)Sn,k,m=0とする。

0128

計算法境界決定部14では、周波数軸での計算を行う最大のk値KM(周波数成分閾値)を求める。KMは0,…,Kの何れかの値を持ち、KM=0は周波数軸での計算を行わずに非特許文献1、2と同様に全て時間軸での計算で求める場合を意味する。

0129

周波数成分閾値KMの決定方法としては、KM=0,…,Kの各場合について、定Q変換の成分Xkcqを求めるのに必要な演算量CX(KM)を評価してそれが最小となるようにKMを決定する。尚、この条件は複数のpmに寄らないので、p1(1つの解析位置)について求めておけば十分である。

0130

定Q変換の成分Xkcqの演算量CX(KM)(合計演算量)は、長さNのFFTの演算量CFFT(N)、非ゼロの時間軸係数毎の演算量CT、非ゼロの時間軸係数の数NT(k)、非ゼロの周波数軸係数毎の演算量CF、kで非ゼロとみなす周波数軸係数の数NF(k)を使っておおよそ、前記実施例1で述べた数式(3)

0131

0132

と表せる。

0133

ここでCFFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、CFFT(N)=(1/4)Nlog2(N)、CT=1、CF=1としても良い。

0134

この演算量CX(KM)は、時間軸係数および周波数軸係数の、ハニング窓での項数と周波数(周波数成分k;対数周波数)の関係を示す図2のように、非ゼロの時間軸係数の数NT(k)がkに関して指数的に減少し、非ゼロの周波数軸係数の数NF(k)がkに関しておおよそ指数的に増加するため、KM=1,…,Kの範囲ではおおよそ下に凸の特性カーブとなる。ただし、KM=0ではCFFTがなくなるためCX(0)が小さくなる。結果として演算量CX(KM)はKM=0あるいはKM=1,…,Kでの最小点の何れかで最小値を取ることになる。

0135

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

0136

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

0137

次に定Q変換の成分Xkcqを、計算法境界決定部14で決定したKMにより、k≦KMについては低周波帯域計算部23で、k>KMについては高周波帯域計算部22で計算する。

0138

0139

尚、前記定Q変換の成分は以下、「Xk,mcq」と表記することもある。

0140

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

0141

低周波帯域計算部23は、まず時系列データxnを全体に対して一度FFTしてフーリエ係数Xnを求める。次にpm毎に、前記周波数軸係数算出部13で算出された周波数軸係数(1/N)Sn,k,mとフーリエ係数Xnを積和計算して定Q変換の成分Xk,mcqを求める。

0142

各kについて非ゼロの周波数軸係数(1/N)Sn,k,mは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。

0143

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

0144

最後に定Q変換合成部24で、高周波帯域計算部22と低周波帯域計算部23とで計算した定Q変換の成分Xk,mcqをまとめて、定Q変換を完了する。

0145

前記実施例2では、時間軸係数Tn,k,m、周波数軸係数(1/N)Sn,k,mは解析位置pm毎に用意するので、事前に用意して保管しておくデータ量が大きくなる。本実施例3では、係数データのメモリ使用量を大幅に減らすように以下のような対策を講じた。

0146

0147

尚、数式(4)中の

0148

0149

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

0150

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

0151

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

0152

0153

を、予め設定した解析位置pmのうち1つの解析位置(例えば1番目の解析位置)p1において算出し、時間軸係数Tn,k,1を求める。

0154

周波数軸係数算出部13では、パラメータ決定部11で決定したパラメータを基に、周波数軸係数(1/N)Sn,k,1を定義式に従って計算する。

0155

すなわち、

0156

0157

から得た(1/N)Sn,kを、前記1つの解析位置p1において算出して周波数軸係数(1/N)Sn,k,1を求める。

0158

この際、|Sn,k,1|≦Thとなる項は周波数軸係数(1/N)Sn,k,1=0とする。

0159

計算法境界決定部14では、周波数軸計算の際に複素べき乗の積和が追加されるために非ゼロの周波数軸係数毎の演算量CFを実施例2の場合より大きく、例えば5(複素べき乗の演算に+3、複素べき乗の積和に+1)とする。

0160

高周波帯域計算部22は、複数の解析位置pm毎に範囲をずらして、前記時間軸係数算出部12で算出された時間軸係数Tn,k,1と時系列データxn-p1+pmを積和計算して定Q変換の成分Xk,mcqを求める。

0161

低周波帯域計算部23は、まず時系列データxnを全体に対して一度FFTしてフーリエ係数Xnを求める。

0162

0163

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

0164

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

0165

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

0166

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

0167

0168

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

0169

0170

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

0171

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

0172

おおむね100時刻以上の一括計算では、実施例2で実施例1の10倍、実施例3で実施例1の5倍の高速化が見込める。

0173

尚、時間軸係数Tn,k,1の非ゼロ項数は実施例1でも実施例2、3でも変わりないが、周波数軸係数(1/N)Sn,k,1の非ゼロ項数については、実施例1と実施例2、3でFFTするデータ数(=行列(1/N)Sn,k,1の次元)が異なるために、実施例1の場合より実施例2、3の場合の方が倍以上多くなる。また一回のFFTに掛かる時間も同様に大幅に大きくなる。これらはもちろん図4に反映している。

実施例

0174

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

0175

10…準備部
11…パラメータ決定部
12…時間軸係数算出部
13…周波数軸係数算出部
14…計算法境界決定部
20…個別計算部
21…前処理部
22…高周波帯域計算部
23…低周波帯域計算部
24…定Q変換合成部

ページトップへ

この技術を出願した法人

この技術を発明した人物

ページトップへ

関連する挑戦したい社会課題

関連する公募課題

ページトップへ

技術視点だけで見ていませんか?

この技術の活用可能性がある分野

分野別動向を把握したい方- 事業化視点で見る -

ページトップへ

おススメ サービス

おススメ astavisionコンテンツ

新着 最近 公開された関連が強い技術

この 技術と関連性が強い人物

関連性が強い人物一覧

この 技術と関連する社会課題

関連する挑戦したい社会課題一覧

この 技術と関連する公募課題

関連する公募課題一覧

astavision 新着記事

サイト情報について

本サービスは、国が公開している情報(公開特許公報、特許整理標準化データ等)を元に構成されています。出典元のデータには一部間違いやノイズがあり、情報の正確さについては保証致しかねます。また一時的に、各データの収録範囲や更新周期によって、一部の情報が正しく表示されないことがございます。当サイトの情報を元にした諸問題、不利益等について当方は何ら責任を負いかねることを予めご承知おきのほど宜しくお願い申し上げます。

主たる情報の出典

特許情報…特許整理標準化データ(XML編)、公開特許公報、特許公報、審決公報、Patent Map Guidance System データ