(57)【要約】【課題】非定常的な時系列信号を解析する。【解決手段】応答関数算出部３により、時系列信号ｓlをＬ（Ｌは正の整数）個ごとの区間に分割し、各区間ごとに連続する時系列信号ｓlを用いて応答関数ａkをＫ（Ｋは正の整数）次まで算出する。次に、複素周波数算出部５により、応答関数ａkを用いて、１-Σkａkｚ-k＝Πk(１-ｅ2πiFkΔTｚ-1)を満たすＫ個の複素周波数Ｆk＝ｆk＋ｉｂkを算出する。【選択図】図１
2005249967.tif 000002
【特許請求の範囲】
【請求項１】
時系列信号を記憶した記憶手段から前記時系列信号を読み出すステップと、
前記時系列信号を有限長の区間に分割するステップと、
前記時系列信号を用いて各区間ごとに時系列信号をｅ2πiFktを含む式で表される波形の集合体として表したときの複素周波数Ｆkを算出するステップと、
を有することを特徴とする周波数解析方法。
【請求項２】
前記時系列信号を用いて応答関数ａkを算出するステップと、
前記応答関数ａkを係数とする１-Σkａkｚ-kとこれを因数分解して得られるΠk(１-ｅ2πiFkΔTｚ-1)とを用いた１-Σkａkｚ-k＝Πk(１-ｅ2πiFkΔTｚ-1)の関係を満たす複素周波数Ｆｋを算出するステップと、
を有することを特徴とする請求項１記載の周波数解析方法。
【請求項３】
前記時系列信号をｓlとし、多項式ｙ(ｔ)＝Σmｃm(ｔ-ｔ0)mを用いたときのΣl‖ｓl-ｙ(ｔ)‖が最小となるように未定係数ｃmを算出するステップと、
この未定係数ｃmと演算子Ｄ＝Πk(d/dt-ａk)とを用いた∫t‖Ｄ・ｙ(t)‖dtが最小となる未定係数ａkを算出するステップと、
算出されたａkを用いてａk＝２πｉＦkの関係を満たす複素周波数Ｆkを算出するステップと、
を有することを特徴とする請求項１記載の周波数解析方法。
【請求項４】
前記時系列信号をｓlとし、これを近似する時系列信号ｓ'lをｓ'l＝Σkａkｓl-kとしたときのΣl‖ｓl-ｓ'l‖が最小となる予測係数ａkを算出するステップと、
前記予測係数ａkを係数とする１-Σkａkｚ-kとこれを因数分解して得られるΠk（１-ｅ2πFkΔTｚ-1）とを用いた１-Σkａkｚ-k＝Πk（１-ｅ2πFkΔTｚ-1）の関係を満たす複素周波数Ｆkを算出するステップと、
を有することを特徴とする請求項１記載の周波数解析方法。
【請求項５】
前記時系列信号をｓlとしたときのＺ変換Σkｓkｚ-kを算出するステップと、
Σkｓkｚ-kを規格化したｃkを係数とする１+Σkｃkｚ-kを用いてａk及びｄkを係数とするΣkｄk/(1-ａkｚ-1)を算出し、Σkｄk/(1-ａkｚ-1)＝Πk(1-ｅ2πiF"kΔTｚ-1)/Πk(1-ｅ2πiF'kΔTｚ-1)の関係を満たすＫ'個の複素周波数Ｆ'k=ｆ'k＋ｉｂ'kおよびＫ"個の複素周波数Ｆ"k=ｆ"k+ｉｂ"kを算出するステップと、
を有することを特徴とする請求項１記載の周波数解析方法。
【請求項６】
前記複素周波数を用いてスペクトル包絡線を算出するステップを有することを特徴とする請求項１乃至５のいずれかに記載の周波数解析方法。
【請求項７】
時系列信号を記憶する記憶手段と、
前記記憶手段から時系列信号を読み出して有限長の区間に分割する分割手段と、
前記時系列信号を用いて各区間ごとに時系列信号をｅ2πiFktを含む式で表される波形の集合体として表したときの複素周波数Ｆkを算出する複素周波数算出手段と、
を有することを特徴とする周波数解析装置。
【請求項８】
前記時系列信号を用いて応答関数ａkを算出する応答関数算出手段を有し、
前記複素周波数算出手段は、前記応答関数ａkを係数とする１-Σkａkｚ-kとこれを因数分解して得られるΠk(１-ｅ2πiFkΔTｚ-1)とを用いた１-Σkａkｚ-k＝Πk(１-ｅ2πiFkΔTｚ-1)の関係を満たす複素周波数Ｆkを算出することを特徴とする請求項７記載の周波数解析装置。
【請求項９】
前記時系列信号をｓlとし、多項式ｙ(ｔ)＝Σmｃm(ｔ-ｔO)mを用いたときのΣl‖ｓl-ｙ(ｔ)‖が最小となるように未定係数ｃmを算出し、この未定係数ｃmと演算子Ｄ＝Πk(d/dt-ａk)とを用いた∫t‖Ｄ・ｙ(t)‖dtが最小となる未定係数ａkを算出する微分係数算出手段を備え、
前記複素周波数算出手段は、算出されたａkを用いてａk＝２πｉＦkの関係を満たす複素周波数Ｆkを算出することを特徴とする請求項７記載の周波数解析装置。
【請求項１０】
前記時系列信号をｓlとし、これを近似する時系列信号ｓ'lをｓ'l＝Σkａkｓl-kとしたときのΣl‖ｓl−ｓ'l‖が最小となる予測係数ａkを算出する予測係数算出手段を備え、
前記複素周波数算出手段は、前記予測係数ａkを係数とする１-Σkａkｚ-kとこれを因数分解して得られるΠk（１-ｅ2πiFkΔTｚ-1）とを用いた１-Σkａkｚ-k＝Πk（１-ｅ2πiFkΔTｚ-1）の関係を満たす複素周波数Ｆkを算出することを特徴とする請求項７記載の周波数解析装置。
【請求項１１】
前記時系列信号をｓlとしたときのＺ変換Σkｓkｚ-kを算出するＺ変換算出手段を備え、
前記複素周波数算出手段は、Σkｓkｚ-kを規格化したｃkを係数とする１+Σkｃkｚ-kを用いてａk及びｄkを係数とするΣkｄk/(1-ａkｚ-1)を算出し、Σkｄk/(1-ａkｚ-1)＝Πk(1-ｅ2πiF"kΔTｚ-1)/Πk(1-ｅ2πiF'kΔTｚ-1)の関係を満たすＫ'個の複素周波数Ｆ'k=ｆ'k＋ｉｂ'kおよびＫ"個の複素周波数Ｆ"k=ｆ"k+ｉｂ"kを算出することを特徴とする請求項７記載の周波数解析装置。
【請求項１２】
前記複素周波数を用いてスペクトル包絡線を算出するスペクトル包絡算出手段を備えることを特徴とする請求項７乃至１１のいずれかに記載の周波数解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、非定常的な信号の時間波形の解析を可能とする周波数解析方法および装置に関する。
【背景技術】
【０００２】
信号波形の解析には、フーリエ変換によって得られるパワースペクトルを用いることが一般的である。しかし、フーリエ変換は、実数周波数を用いて時系列を特徴つけようとするものであるため、時系列に周期性のあることが前提とされる。このため、過渡現象の解析や、例えばランダムに発生するパルス的な信号や、減衰・発散などを伴う信号といった非定常的な信号の解析には不適であるという問題がある。
【０００３】
過渡現象の解析は、ウェーブレット解析を用いることにより可能であるが、解析によって得られる情報の物理的意味が不明確であるという問題がある。また、ウェーブレット解析は計算量が多く、リアルタイム解析が困難であるという問題がある。さらには、時間分解能を高めようとすると周波数分解能が低下するといった、時間・周波数分解能の両立が困難であるという問題がある。
【０００４】
一方、人の音声を解析する場合には、ホルマント解析が用いられるが、ホルマント解析を基礎とする技術は、信号が区間的に定常的であることを前提としており、区間的にも定常性を持たない信号の解析には用いられていなかった。人の音声以外へのホルマント解析の適用としては、例えば特許文献１が知られているが、この文献においても信号の定常性が前提とされている。
【特許文献１】特許第３４５３１３０号公報
【発明の開示】
【発明が解決しようとする課題】
【０００５】
このように、従来の技術では、非定常的な信号を有効に解析できないという問題がある。
【０００６】
本発明は、上記に鑑みてなされたものであり、その目的とするとするところは、非定常的な時系列信号を解析し得る周波数解析方法および装置を提供することにある。
【課題を解決するための手段】
【０００７】
第１の本発明に係る周波数解析方法は、時系列信号を記憶した記憶手段から前記時系列信号を読み出すステップと、前記時系列信号を有限長の区間に分割するステップと、前記時系列信号を用いて各区間ごとに時系列信号をｅ2πiFktの項を含む式で表される波形の集合体として表したときの複素周波数Ｆkを算出するステップと、を有することを特徴とする。
【０００８】
本発明にあっては、時系列信号を有限長の区間に分割し、各区間ごとに時系列信号をｅ2πiFktを含む式で表される波形の集合体として特徴づける１個以上の複素周波数Ｆk＝ｆk＋ｉｂkを算出することで、フーリエ変換では得ることのできない虚数部ｂkが得られるようにして、非定常的な信号の時間波形の解析を可能としている。
【０００９】
第２の本発明に係る周波数解析方法は、前記時系列信号を用いて応答関数ａkを算出するステップと、前記応答関数ａkを係数とする１-Σkａkｚ-kとこれを因数分解して得られるΠk(１-ｅ2πiFkΔTｚ-1)とを用いた１-Σkａkｚ-k＝Πk(１-ｅ2πiFkΔTｚ-1)の関係を満たす複素周波数Ｆkを算出するステップと、を有することを特徴とする。
【００１０】
本発明にあっては、時系列信号Ｓlを用いて応答関数ａkを算出し、１-Σkａkｚ-k＝Πk(１-ｅ2πFkΔTz-1)の関係を満たす複素周波数Ｆkを算出することで、時系列信号を特徴づける１個以上の複素周波数を算出するようにしている。
【００１１】
第３の本発明に係る周波数解析方法は、前記時系列信号をｓlとし、多項式ｙ(ｔ)＝Σmｃm(ｔ-ｔO)mを用いたときのΣl‖ｓl-ｙ(ｔ)‖が最小となるように未定係数ｃmを算出するステップと、この未定係数ｃmと演算子Ｄ＝Πk(d/dt-ａk)とを用いた∫t‖Ｄ・ｙ(t)‖dtが最小となる未定係数ａkを算出するステップと、算出されたａkを用いてａk＝２πｉＦkの関係を満たす複素周波数Ｆkを算出するステップと、を有することを特徴とする。
【００１２】
本発明にあっては、時系列信号ｓlと多項式ｙ(ｔ)＝Σmｃm(ｔ-ｔO)mとを用いたΣl‖ｓl-ｙ(ｔ)‖が最小となるように未定係数ｃmを算出するとともに、演算子Ｄ＝Πk(d/dt-ａk)を用いた∫t‖Ｄ・ｙ(t)‖dtが最小となる未定係数ａkを算出した後、算出された係数ａkを用いてａk＝２πｉＦkなる関係を満たす複素周波数Ｆkを算出することで、時系列信号を特徴づける１個以上の複素周波数を算出するようにしている。
【００１３】
第４の本発明に係る周波数解析方法は、前記時系列信号をｓlとし、これを近似する時系列信号ｓ'lをｓ'l＝Σkａkｓl-kとしたときのΣl‖ｓl-ｓ'l‖が最小となる予測係数ａkを算出するステップと、前記予測係数ａkを係数とする１-Σkａkｚ-kとこれを因数分解して得られるΠk（１-ｅ2πiFkΔTｚ-1）とを用いた１-Σkａkｚ-k＝Πk（１-ｅ2πiFkΔTｚ-1）の関係を満たす複素周波数Ｆkを算出するステップと、を有することを特徴とする。
【００１４】
本発明にあっては、時系列信号Ｓlを用いて予測係数ａkを算出し、１-Σkａkｚ-k＝Πk(１-ｅ2πFkΔTｚ-1)の関係を満たす複素周波数Ｆkを算出することで、時系列信号を特徴づける１個以上の複素周波数を算出するようにしている。
【００１５】
第５の本発明に係る周波数解析方法は、前記時系列信号をｓlとしたときのＺ変換Σkｓkｚ-kを算出するステップと、Σkｓkｚ-kを規格化したｃkを係数とする１+Σkｃkｚ-kを用いてａk及びｄkを係数とするΣkｄk/(1-ａkｚ-1)を算出し、Σkｄk/(1-ａkｚ-1)＝Πk(1-ｅ2πiF"kΔTｚ-1)/Πk(1-ｅ2πiF'kΔTｚ-1)の関係を満たすＫ'個の複素周波数Ｆ'k=ｆ'k＋ｉｂ'kおよびＫ"個の複素周波数Ｆ"k=ｆ"k+ｉｂ"kを算出するステップと、を有することを特徴とする。
【００１６】
本発明にあっては、時系列信号ｓlをＺ変換してΣkｓkｚ-kを算出し、このΣkｓkｚ-kを規格化したｃkを係数とする１+Σkｃkｚ-kを用いてａk及びｄkを係数とするΣkｄk/(1-ａkｚ-1)を算出し、Σkｄk/(1-ａkｚ-1)＝Πk(1-ｅ2πiF"kΔTｚ-1)/Πk(1-ｅ2πiF'kΔTｚ-1)の関係を満たすＫ'個の複素周波数Ｆ'k=ｆ'k＋ｉｂ'kおよびＫ"個の複素周波数Ｆ"k=ｆ"k+ｉｂ"kを算出することで、時系列信号を特徴づける１個以上の複素周波数を算出するようにしている。
【００１７】
第６の本発明に係る周波数解析方法は、前記複素周波数を用いてスペクトル包絡線を算出するステップを有することを特徴とする。
【００１８】
本発明にあっては、複素周波数を用いてスペクトル包絡線を算出することで、特徴的な周波数をスペクトル包絡線のピークに位置する周波数として見出すことができる。
【００１９】
第７の本発明に係る周波数解析装置は、時系列信号を記憶する記憶手段と、前記記憶手段から時系列信号を読み出して有限長の区間に分割する分割手段と、前記時系列信号を用いて各区間ごとに時系列信号をｅ2πiFktを含む式で表される波形の集合体として表したときの複素周波数Ｆkを算出する複素周波数算出手段と、を有することを特徴とする。
【００２０】
第８の本発明に係る周波数解析装置は、前記時系列信号を用いて応答関数ａkを算出する応答関数算出手段を有し、前記複素周波数算出手段は、前記応答関数ａkを係数とする１-Σkａkｚ-kとこれを因数分解して得られるΠk(１-ｅ2πiFkΔTｚ-1)とを用いた１-Σkａkｚ-k＝Πk(１-ｅ2πiFkΔTｚ-1)の関係を満たす複素周波数Ｆkを算出することを特徴とする。
【００２１】
第９の本発明に係る周波数解析装置は、前記時系列信号をｓlとし、多項式ｙ(ｔ)＝Σmｃm(ｔ-ｔO)mを用いたときのΣl‖ｓl-ｙ(ｔ)‖が最小となるように未定係数ｃmを算出し、この未定係数ｃmと演算子Ｄ＝Πk(d/dt-ａk)とを用いた∫t‖Ｄ・ｙ(t)‖dtが最小となる未定係数ａkを算出する微分係数算出手段を備え、前記複素周波数算出手段は、算出されたａkを用いてａk＝２πｉＦkの関係を満たす複素周波数Ｆkを算出することを特徴とする。
【００２２】
第１０の本発明に係る周波数解析装置は、前記時系列信号をｓlとし、これを近似する時系列信号ｓ'lをｓ'l＝Σkａkｓl-kとしたときのΣl‖ｓl−ｓ'l‖が最小となる予測係数ａkを算出する予測係数算出手段を備え、前記複素周波数算出手段は、前記予測係数ａkを係数とする１-Σkａkｚ-kとこれを因数分解して得られるΠk（１-ｅ2πiFkΔTｚ-1）とを用いた１-Σkａkｚ-k＝Πk（１-ｅ2πiFkΔTｚ-1）の関係を満たす複素周波数Ｆkを算出することを特徴とする。
【００２３】
第１１の本発明に係る周波数解析装置は、前記時系列信号をｓlとしたときのＺ変換Σkｓkｚ-kを算出するＺ変換算出手段を備え、前記複素周波数算出手段は、Σkｓkｚ-kを規格化したｃkを係数とする１+Σkｃkｚ-kを用いてａk及びｄkを係数とするΣkｄk/(1-ａkｚ-1)を算出し、Σkｄk/(1-ａkｚ-1)＝Πk(1-ｅ2πiF"kΔTｚ-1)/Πk(1-ｅ2πiF'kΔTｚ-1)の関係を満たすＫ'個の複素周波数Ｆ'k=ｆ'k＋ｉｂ'kおよびＫ"個の複素周波数Ｆ"k=ｆ"k+ｉｂ"kを算出することを特徴とする。
【００２４】
第１２の本発明に係る周波数解析装置は、前記複素周波数を用いてスペクトル包絡線を算出するスペクトル包絡算出手段を備えることを特徴とする。
【発明の効果】
【００２５】
本発明に係る周波数解析方法および装置によれば、非定常的な時系列信号を解析することができる。
【発明を実施するための最良の形態】
【００２６】
以下、本発明の実施の形態について図面を用いて説明する。
【００２７】
［第１の実施の形態］
図１は、本実施の形態における周波数解析装置の構成を示す機能ブロック図である。同図の周波数解析装置は、測定部１と、アナログ／デジタル（Ａ／Ｄ）変換部２と、応答関数算出部３と、応答関数最適化部４と、複素周波数算出部５と、表示部６と、記憶部１１とを有する構成である。
【００２８】
本周波数解析装置は、専用のハードウェアによって構成するものとしてもよいし、コンピュータにより構成し、各部における処理を、コンピュータにインストールされたプログラムによって実行するようにしてもよい。記憶部１１は、磁気記憶装置、光ディスク装置、メモリ、といったデータを記憶可能な装置により構成される。
【００２９】
測定部１は、非定常的なアナログ信号ｓ(ｔ)を測定してＡ／Ｄ変換部２へ出力する。Ａ／Ｄ変換部２では、例えば２０ｋＨｚまでの周波数解析に対してはサンプリング周波数を４０ｋＨｚ以上とし、５００ＭＨｚまでの周波数解析に対してはサンプリング周波数を１ＧＨｚ以上に設定するなど、解析を行う周波数の上限値の２倍以上の周波数に設定し、アナログ信号ｓ(ｔ)をデジタル化して時系列信号ｓlを得て、この時系列信号ｓlを記憶部１１に一時的に記憶させる。ここでは、サンプリング周波数から定まるサンプリング間隔をΔＴとする。Ａ／Ｄ変換部２は、記憶部１１から時系列信号ｓlを読み出し、応答関数算出部３へ出力する。
【００３０】
応答関数算出部３は、時系列信号ｓlをＬ（Ｌは正の整数）個ごとの区間に分割し、各区間ごとに、連続する時系列信号ｓlを用いて次の連立方程式をａkについて解くことにより応答関数ａkをＫ（Ｋは正の整数）次まで算出し、メモリに格納する。
【００３１】
ｓl＝Σkａkｓl-k
ｓl-1＝Σkａkｓl-1-k
…
ｓl-K＝Σkａkｓl-K-k  …(1)
ここで、ｋの範囲は１〜Ｋである。Ｋは、算出しようとする複素周波数の個数であり、２以上の値とする。Ｌは、必要に応じてＫの２倍以上の値とする。Ｌの値をＫよりも大きくとる場合には、応答関数最適化部４により最小自乗法などを用いて最適な応答関数ａｋを算出する。
【００３２】
次に、応答関数ａkを係数とする次式について考える。
【００３３】
１-Σkａkｚ-k        …(2)
ここで、ｋの範囲は１〜Ｋである。式(2)を因数分解することを考えると次式が成立する。
【００３４】
１-Σkａkｚ-k
＝Πk(１-ｅ2πiFkΔTｚ-1) …(3)
複素周波数算出部５は、予め求めておいた応答関数ａkをメモリから読み出し、この応答関数ａkを用いて、式(3)を満たすＫ個の複素周波数Ｆk＝ｆk＋ｉｂkを算出し、表示部６へ出力する。
【００３５】
表示部６は、図２に示すように、横軸を時間、縦軸を大きさとするグラフを示す画面上で、ｆkとｂkを表示する。
【００３６】
時系列信号ｓlの波形は、ΣkΣt'kct'kｈ(ｔ-ｔ'k)ｅ2πiFk(t-t'k)という形で表すことができる。ここで、ｔ'kは各波形の発生時刻であり、Ｆkごとに複数回あってもよい。ｃt'kは、各波形の大きさ、ｈ(ｔ)はｔ＜０の範囲で０、ｔ≧０の範囲で１をとる関数である。これは、以降の実施の形態においても同様である。
【００３７】
したがって、本実施の形態によれば、時系列信号ｓlを有限長の区間に分割し、各区間ごとに時系列信号をｅ2πiFktを含む式で表される波形の集合体として特徴づける１個以上の複素周波数Ｆk＝ｆk＋ｉｂkを算出することで、フーリエ変換では得ることのできない虚数部ｂkが得られ、非定常的な信号の時間波形を解析することができる。
【００３８】
本実施の形態によれば、応答関数算出部３により時系列信号ｓlを用いて応答関数ａkを算出し、複素周波数算出部５により１-Σkａkｚ-k＝Πk(１-ｅ2πiFkΔTｚ-1)の関係を満たす複素周波数Ｆkを算出することで、時系列信号を特徴づける１個以上の複素周波数が算出され、非定常的な信号の時間波形を解析することができる。
【００３９】
［第２の実施の形態］
図３は、本実施の形態における周波数解析装置の構成を示す機能ブロック図である。同図の周波数解析装置は、図１の応答関数算出部３と応答関数最適化部４に代えて、微分係数算出部７を備えた構成である。その他、図１と同一物あるいは近似物については同一の符号を付すものとし、ここでは重複した説明は省略する。
【００４０】
時系列信号ｓlを、ｃmを未定係数とする次式の多項式でモデル化することを考える。
【００４１】
ｙ(t)＝Σmｃm(ｔ-ｔO)m …(4)
ここで、ｍの範囲は０〜Ｍ（Ｍは正の整数）である。微分係数算出部７は、記憶部１１から読み出された時系列信号ｓlを用いて、例えば最小自乗法により次式が最小となる未定係数ｃmを算出し、メモリに格納する。
【００４２】
Σl‖ｓl-ｙ(ｔl)‖     …(5)
ここで、ｌの範囲は０〜Ｌ（Ｌは正の整数）である。次に、ａkを未定係数とする線形微分方程式で表される次式の演算子Ｄを考える。
【００４３】
Ｄ＝Πk(d/dt-ａk)      …(6)
微分係数算出部７では、この演算子Ｄに含まれる未定係数ａkについて、メモリから読み出した未定係数ｃmを用いて、例えば最小自乗法により次式が最小となる未定係数ａkを算出し、メモリに格納する。
【００４４】
∫t‖Ｄ・ｙ(t)‖dt    …(7)
ここで、ｔの範囲はｔa〜ｔbである。ｔlは時系列信号ｓlのｌ番目の値を与える時刻であり、ｔaとｔbは、係数ａkを算出するのに用いられる時系列信号ｓlの先頭値を与える時刻ｔOと最終値を与える時刻ｔＬとを用いて、不等式ｔO≦ｔa≦ｔb≦ｔＬを満たす範囲で与えられる。
【００４５】
複素周波数算出部５では、メモリから係数ａkを読み出し、この係数ａkを用いて次の関係式を満たす複素周波数Ｆk＝ｆk＋ｉｂkを算出する。
【００４６】
ａk＝２πｉＦk     …(8)
表示部６は、このようにして各区間ごとに得られたＫ個の複素周波数Ｆk＝ｆk＋ｉｂkを表示する。ここでは、図４に示すように、複素周波数の実部を横軸、虚部を縦軸とするグラフを示す画面上で複素周波数Ｆkを表示する。もちろん、図２に示したように、時間軸を横軸、大きさを縦軸にして表示するようにしてもよい。
【００４７】
したがって、本実施の形態によれば、微分係数算出部７により、時系列信号ｓlと多項式ｙ(ｔ)＝Σmｃm(ｔ-ｔO)mとを用いたΣl‖ｓl-ｙ(ｔ)‖が最小となるように未定係数ｃmを算出するとともに、演算子Ｄ＝Πk(d/dt-ａk)を用いた∫t‖Ｄ・ｙ(t)‖dtが最小となる未定係数ａkを算出した後、複素周波数算出部５により、係数ａkを用いてａk＝２πｉＦkなる関係を満たす複素周波数Ｆkを算出することで、時系列信号を特徴づける１個以上の複素周波数が算出され、非定常的な信号の時間波形を解析することができる。
【００４８】
［第３の実施の形態］
図５は、本実施の形態における周波数解析装置の構成を示す機能ブロック図である。同図の周波数解析装置は、図１の応答関数算出部３と応答関数最適化部４に代えて、予測係数算出部８を備えるとともに、複素周波数算出部５の出力段にスペクトル包絡算出部９を備えた構成である。その他、図１と同一物あるいは近似物については同一の符号を付すものとし、ここでは重複した説明は省略する。
【００４９】
予測係数算出部８では、線形予測法（Linear Predictive Coding, LPC）による予測係数ａkを算出する。具体的には、次の処理を行う。
【００５０】
まず、Ａ／Ｄ変換部２により記憶部１１から読み出された時系列信号ｓlを近似する次の時系列信号ｓ'lを考える。
【００５１】
ｓ'l＝Σkａkｓl-k   …(9)
ここで、ｋの範囲は１〜Ｋである。予測係数算出部８では、この時系列信号ｓ'lを用いて、次式が最小となる予測係数ａkを算出し、メモリに格納する。
【００５２】
Σl‖ｓl-ｓ'l‖    …(10)
ここで、ｌの範囲はｌO−Ｌ＋１〜ｌOである。次に、予測係数ａkを用いた次式を考える。
【００５３】
１-Σkａkｚ-k        …(11)
ここで、ｋの範囲は１〜Ｋである。式(11)を因数分解することを考えると次式が成立する。
【００５４】
１-Σkａkｚ-k
＝Πk(１-ｅ2πiFkΔTｚ-1)  …(12)
複素周波数算出部５は、メモリから予測係数ａkを読み出し、この予測係数ａkを用いて式(12)を満たすＫ個の複素周波数Ｆk＝ｆk＋ｉｂkを算出し、スペクトル包絡算出部９へ出力する。
【００５５】
スペクトル包絡算出部９は、次式によってスペクトル包絡線Ｇ(ｆ)を算出する。
【００５６】
Ｇ(f)=1/｜1-ｅ2πi(Fk-f)ΔT｜2  …(13)
表示部６は、図６に示すように、横軸を周波数ｆとするグラフを示す画面上で、スペクトル包絡線Ｇ(ｆ)を表示する。
【００５７】
したがって、本実施の形態によれば、予測係数算出部８により時系列信号ｓlを用いて予測係数ａkを算出し、複素周波数算出部５により１-Σkａkｚ-k＝Πk(１-ｅ2πiFkΔTｚ-1)の関係を満たす複素周波数Ｆkを算出することで、時系列信号を特徴づける１個以上の複素周波数が算出され、非定常的な信号の時間波形を解析することができる。
【００５８】
本実施の形態によれば、スペクトル包絡算出部９により、Ｋ個の複素周波数Ｆk＝ｆk＋ｉｂkを用いてスペクトル包絡線Ｇ(ｆ)を算出するようにしたことで、特徴的な周波数ｆをスペクトル包絡線のピークに位置する周波数として見出すことができる。
【００５９】
なお、表示部６では、上記各実施の形態で説明したような表示を行ってもよい。また、上記各実施の形態において、表示部６で、スペクトル包絡線を表示するようにしてもよい。
【００６０】
［第４の実施の形態］
図７は、本実施の形態における周波数解析装置の構成を示す機能ブロック図である。同図の周波数解析装置は、図１の応答関数算出部３と応答関数最適化部４に代えて、Ｚ変換算出部１０を備えた構成である。その他、図１と同一物あるいは近似物については同一の符号を付すものとし、ここでは重複した説明は省略する。
【００６１】
Ｚ変換算出部１０では、Ａ／Ｄ変換部２により記憶部１１から読み出された連続するＫ＋１個の時系列信号ｓlを次式によりＺ変換する。
【００６２】
Σkｓkｚ-k  …(14)
ここで、ｋの範囲は０〜Ｋである。Ｚ変換算出部１０は、これを規格化してｃkを係数とする次式を与える。
【００６３】
１＋Σkｃkｚ-k   …(15)
ここで、ｋの範囲は１〜Ｋである。Ｋは、求めようとする複素周波数の個数であり、２以上の値とする。
【００６４】
複素周波数算出部５では、式(15)を用いてａk及びｄkを係数とする次式を算出する。
【００６５】
Σkｄk/(1-ａkｚ-1)  …(16)
ここで、ｋの範囲は１〜Ｋ'(Ｋ'は正の整数)である。式(16)の算出は、次式の対応関係に従って行う。
【００６６】
Σkｄk/(1-ａkｚ-1)
＝Σkｄk(1+Σmａkmｚ-m)       …(17)
＝Σkｄk+Σm(Σkｄkａkm)ｚ-m    …(18)
＝１+Σkｃkｚ-k+Ｏ(ｚ-K-1)     …(19)
ここで、式(17)，式(18)におけるｋの範囲は１〜Ｋ'であり、ｍの範囲は１〜∞である。式(19)におけるｋの範囲は１〜Ｋである。また、Ｏ(ｚ-K-1)は、複素周波数の算出にあたっての残差項である。
【００６７】
複素周波数算出部５では、式(16)を変形して次式を得る。
【００６８】
Σkｄk/(1-ａkｚ-1)
＝ΣkｄkΠk'≠k(1-ａk'ｚ-1)/Πk(1-ａkｚ-1)      …(20)
＝Πk(1-ａ'kｚ-1)/Πk(1-ａkｚ-1)             …(21)
＝Πk(1-ｅ2πiF"kΔTｚ-1)/Πk(1-ｅ2πiF'kΔTｚ-1)  …(22)
ここで、式(20)におけるｋ，ｋ'の範囲は、１〜Ｋである。式(21)，(22)の分子におけるｋの範囲は１〜Ｋ"、分母におけるｋの範囲は１〜Ｋ'である。
【００６９】
続いて、複素周波数算出部５は、式(22)の関係を満たすＫ'個の複素周波数Ｆ'k=ｆ'k＋ｉｂ'kおよびＫ"個の複素周波数Ｆ"k=ｆ"k+ｉｂ"kを算出する。
【００７０】
表示部６では、横軸を実部、縦軸を虚部とするグラフ、および横軸を時間、縦軸を虚部とするグラフを示す画面上でＫ'個の複素周波数Ｆ'kと、Ｋ"個の複素周波数Ｆ"kを表示する。
【００７１】
したがって、本実施の形態によれば、Ｚ変換算出部１０により時系列信号ｓlをＺ変換してΣkｓkｚ-kを算出し、複素周波数算出部５により、Σkｓkｚ-kを規格化したｃkを係数とする１+Σkｃkｚ-kを用いてａk及びｄkを係数とするΣkｄk/(1-ａkｚ-1)を算出し、Σkｄk/(1-ａkｚ-1)＝Πk(1-ｅ2πiF"kΔTｚ-1)/Πk(1-ｅ2πiF'kΔTｚ-1)の関係を満たすＫ'個の複素周波数Ｆ'k=ｆ'k＋ｉｂ'kおよびＫ"個の複素周波数Ｆ"k=ｆ"k+ｉｂ"kを算出することで、時系列信号を特徴づける１個以上の複素周波数が算出され、非定常的な信号の時間波形を解析することができる。
【００７２】
［実施例］
以下、実際の時系列信号を用いて複素周波数を算出したときの結果について説明する。
【００７３】
図９は、ｙ(ｔ)＝ｅ‐λtｓｉｎ２πωｔで表される波形をランダムに発生させた時系列信号の波形を示す図である。ここでは、周波数ωをゆっくり上下させている。また、波形の減衰速度λは定数で与えている。
【００７４】
図１０は、図９の時系列信号をフーリエ変換した波形を示す図である。図１０に示すように、フーリエ解析では、非定常的な信号に含まれる特徴的な周波数ωを見出すことができない。
【００７５】
図１１は、図９の時系列信号に対して、第３の実施の形態で説明したように、予測係数を用いて複素周波数を算出し、その複素周波数を用いて算出したスペクトル包絡線を示す図である。図１１に示すように、本実施例では、非定常的な信号に含まれる特徴的な周波数ωを、スペクトル包絡線のピークｐに位置する周波数として見出すことができる。
【００７６】
図１２は、図１１におけるピークｐの時間変化をプロットしたグラフである。図１１に周波数軸と示すラインは、図１１における周波数軸に対応するものである。このように、本実施例では、非定常的な信号に含まれる特徴的な周波数ωの時間変化も見出すことができる。
【００７７】
ホルマント解析は、区間ごとに定常的な可聴周波数帯域の音声信号の解析に用いるものと認識されているが、第３の実施の形態および本実施例によれば、定常性を持たないＧＨｚ帯の電波や電気信号などの解析に適用できることが明らかとなった。また、それらの単発のパルス信号に対しても適用することができる。
【００７８】
なお、上記各実施の形態における技術は、時間軸上の時系列信号の解析だけではなく、模様やパターンを示す空間的な時系列信号の解析にも適用することができる。
【図面の簡単な説明】
【００７９】
【図１】第１の実施の形態における周波数解析装置の構成を示す機能ブロック図である。
【図２】第１の実施の形態で算出した複素周波数を表示した様子を示す図である。
【図３】第２の実施の形態における周波数解析装置の構成を示す機能ブロック図である。
【図４】第２の実施の形態で算出した複素周波数を表示した様子を示す図である。
【図５】第３の実施の形態における周波数解析装置の構成を示す機能ブロック図である。
【図６】第３の実施の形態で算出した複素周波数を表示した様子を示す図である。
【図７】第４の実施の形態における周波数解析装置の構成を示す機能ブロック図である。
【図８】第４の実施の形態で算出した複素周波数を表示した様子を示す図である。
【図９】ｙ(ｔ)＝ｅ‐λtｓｉｎ２πωｔという形式の波形をランダムに発生させた時系列信号を示す波形を示す図である。
【図１０】図９の時系列信号をフーリエ変換した波形を示す図である。
【図１１】図９の時系列信号に対して、予測係数を用いて複素周波数を算出し、この複素周波数を用いて算出したスペクトル包絡線を示す図である。
【図１２】図１１におけるピークｐの時間変化をプロットしたグラフを示す図である。
【符号の説明】
【００８０】
１…測定部
２…変換部
３…応答関数算出部
４…応答関数最適化部
５…複素周波数算出部
６…表示部
７…微分係数算出部
８…予測係数算出部
９…スペクトル包絡算出部
１０…Ｚ変換算出部
１１…記憶部
【図１】
2005249967.tif 000003
【図２】
2005249967.tif 000004
【図３】
2005249967.tif 000005
【図４】
2005249967.tif 000006
【図５】
2005249967.tif 000007
【図６】
2005249967.tif 000008
【図７】
2005249967.tif 000009
【図８】
2005249967.tif 000010
【図９】
2005249967.tif 000011
【図１０】
2005249967.tif 000012
【図１１】
2005249967.tif 000013
【図１２】
2005249967.tif 000014