JP3836847B2 - 周波数解析方法および装置 - Google Patents

Description

本発明は、非定常的な信号の時間波形の解析を可能とする周波数解析方法および装置に関する。

信号波形の解析には、フーリエ変換によって得られるパワースペクトルを用いることが一般的である。しかし、フーリエ変換は、実数周波数を用いて時系列を特徴つけようとするものであるため、時系列に周期性のあることが前提とされる。このため、過渡現象の解析や、例えばランダムに発生するパルス的な信号や、減衰・発散などを伴う信号といった非定常的な信号の解析には不適であるという問題がある。

過渡現象の解析は、ウェーブレット解析を用いることにより可能であるが、解析によって得られる情報の物理的意味が不明確であるという問題がある。また、ウェーブレット解析は計算量が多く、リアルタイム解析が困難であるという問題がある。さらには、時間分解能を高めようとすると周波数分解能が低下するといった、時間・周波数分解能の両立が困難であるという問題がある。

一方、人の音声を解析する場合には、ホルマント解析が用いられるが、ホルマント解析を基礎とする技術は、信号が区間的に定常的であることを前提としており、区間的にも定常性を持たない信号の解析には用いられていなかった。人の音声以外へのホルマント解析の適用としては、例えば特許文献１が知られているが、この文献においても信号の定常性が前提とされている。
特許第３４５３１３０号公報

このように、従来の技術では、非定常的な信号を有効に解析できないという問題がある。

本発明は、上記に鑑みてなされたものであり、その目的とするとするところは、非定常的な時系列信号を解析し得る周波数解析方法および装置を提供することにある。

第１の本発明に係る周波数解析方法は、時系列信号を記憶した記憶手段から前記時系列信号を読み出すステップと、前記時系列信号を有限長の区間に分割するステップと、前記時系列信号を用いて各区間ごとに時系列信号をｅ^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を係数とする次式について考える。

１-Σ_ｋａ_ｋｚ^-ｋ …(2)
ここで、ｋの範囲は１〜Ｋである。式(2)を因数分解して１次の因数の積で表すことを考えると、方程式の根をｅ ^{２πｉＦｋΔＴ} ｚ ^-１ とするとき、式(2)は次のように書き換えることができる。

１-Σ_ｋａ_ｋｚ^-ｋ
＝Π_ｋ(１-ｅ^{２πｉＦｋΔＴ}ｚ^-１) …(3)
複素周波数算出部５は、予め求めておいた応答関数ａ_ｋをメモリから読み出し、この応答関数ａ_ｋを用いて、式(3)を満たすＫ個の複素周波数Ｆ_ｋ＝ｆ_ｋ＋ｉｂ_ｋを算出し、表示部６へ出力する。ここで、ｆ _ｋ は周波数項、ｂ _ｋ は減衰項である。

表示部６は、図２に示すように、横軸を時間、縦軸を大きさとするグラフを示す画面上で、ｆ_kとｂ_kを表示する。

時系列信号ｓ_lの波形は、Σ_kΣ_t'kc_t'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を算出し、メモリに格納する。

Σ_ｌ‖ｓ_ｌ-ｙ(ｔ_ｌ)‖ …(5)
ここで、ｌの範囲は０〜Ｌ（Ｌは正の整数）である。次に、未定係数ａ _ｋ を根とする線形微分方程式で表される次式の微分演算子Ｄを考える。

Ｄ＝Π_ｋ(d/dt-ａ_ｋ) …(6)
微分係数算出部７では、この微分演算子Ｄに含まれる未定係数ａ_ｋについて、メモリから読み出した未定係数ｃ_ｍを用いて、例えば最小自乗法により次式が最小となる未定係数ａ_ｋを算出し、メモリに格納する。

∫_ｔ‖Ｄ・ｙ(t)‖dt …(7)
ここで、ｔの範囲はｔ_ａ〜ｔ_ｂである。ｔ_ｌは時系列信号ｓ_ｌのｌ番目の値を与える時刻であり、ｔ_ａとｔ_ｂは、係数ａ_ｋを算出するのに用いられる時系列信号ｓ_ｌの先頭値を与える時刻ｔ_０と最終値を与える時刻ｔ_Ｌとを用いて、不等式ｔ_０≦ｔ_ａ≦ｔ_ｂ≦ｔ_Ｌを満たす範囲で与えられる。
式(7)を最小にする未定係数ａ _ｋ の算出は、∫ _ｔ ‖Ｄ・ｙ(t)‖dtをａ _ｋ について偏微分したときの∂/∂ａ _ｋ ∫ _ｔ ‖Ｄ・ｙ(t)‖dt＝０を満たすａ _ｋ を求めるようにする。

複素周波数算出部５では、メモリから係数ａ_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}｜² …(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ａ_k ^mｚ^-m) …(17)
＝Σ_kｄ_k+Σ_m(Σ_kｄ_kａ_k ^m)ｚ^-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ｓｉｎ２πωｔという形式の波形をランダムに発生させた時系列信号を示す波形を示す図である。 図９の時系列信号をフーリエ変換した波形を示す図である。 図９の時系列信号に対して、予測係数を用いて複素周波数を算出し、この複素周波数を用いて算出したスペクトル包絡線を示す図である。 図１１におけるピークｐの時間変化をプロットしたグラフを示す図である。

符号の説明

１…測定部
２…変換部
３…応答関数算出部
４…応答関数最適化部
５…複素周波数算出部
６…表示部
７…微分係数算出部
８…予測係数算出部
９…スペクトル包絡算出部
１０…Ｚ変換算出部
１１…記憶部

Claims (11)

1. 分割手段が、各波形の発生時刻をｔ' _ｋ 、各波形の大きさをｃ _ｔ'ｋ とし、ｔ＜０の範囲で０をとりｔ≧０の範囲で１をとる関数ｈ(ｔ)、複素周波数Ｆ _ｋ を用いてΣ _ｋ Σ _ｔ'k c _ｔ'k ｈ(ｔ-ｔ' _ｋ )ｅ ^{２πｉＦｋ(ｔ-t'k)} という形で表すことができる非定常的な時系列信号ｓ _ｌ を記憶した記憶手段から前記時系列信号を読み出して有限長の区間に分割するステップと、
   応答関数算出手段が、各区間ごとに前記時系列信号を用いて次の連立方程式
   ｓ _ｌ ＝Σ _ｋ ａ _ｋ ｓ _ｌ-ｋ
   ｓ _ｌ-１ ＝Σ _ｋ ａ _ｋ ｓ _ｌ-１-ｋ
   …
   ｓ _ｌ-Ｋ ＝Σ _ｋ ａ _ｋ ｓ _ｌ-Ｋ-ｋ
   をａ _ｋ について解くことにより応答関数ａ_ｋを算出して記憶手段に記憶させるステップと、
   複素周波数算出手段が、前記算出された応答関数ａ_ｋを係数とする次式１-Σ_ｋａ_ｋｚ^-ｋ を、ｅ ^{２πｉＦｋΔＴ} ｚ ^-１ を根としたときの１次の因数の積に因数分解して得られる１-Σ_ｋａ_ｋｚ^-ｋ＝Π_ｋ(１-ｅ^{２πｉＦｋΔＴ}ｚ^-１)の関係を満たす複素周波数Ｆ_ｋを算出して記憶手段に記憶させるステップと、
   を有することを特徴とする周波数解析方法。
2. 分割手段が、各波形の発生時刻をｔ' _ｋ 、各波形の大きさをｃ _ｔ'ｋ とし、ｔ＜０の範囲で０をとりｔ≧０の範囲で１をとる関数ｈ(ｔ)、複素周波数Ｆ _ｋ を用いてΣ _ｋ Σ _ｔ'k c _ｔ'k ｈ(ｔ-ｔ' _ｋ )ｅ ^{２πｉＦｋ(ｔ-t'k)} という形で表すことができる非定常的な時系列信号ｓ _ｌ を記憶した記憶手段から前記時系列信号を読み出して有限長の区間に分割するステップと、
   微分係数算出手段が、各区間ごとに多項式ｙ(ｔ)＝Σ_ｍｃ_ｍ(ｔ-ｔ_０)^ｍを用いたときのΣ_ｌ‖ｓ_ｌ-ｙ(ｔ)‖が最小となるように未定係数ｃ_ｍを算出し、この未定係数ｃ_ｍとａ _ｋ を根とする微分演算子Ｄ＝Π_ｋ(d/dt-ａ_ｋ)とを用いた∫_ｔ‖Ｄ・ｙ(t)‖dtについて、ａ _ｋ で偏微分したときの∂/∂ａ _ｋ ∫ _ｔ ‖Ｄ・ｙ(t)‖dt＝０を満たす未定係数ａ_ｋを算出して記憶手段に記憶させるステップと、
   複素周波数算出手段が、算出されたａ_ｋを用いてａ_ｋ＝２πｉＦ_ｋの関係を満たす複素周波数Ｆ_ｋを算出して記憶手段に記憶させるステップと、
   を有することを特徴とする周波数解析方法。
3. 分割手段が、各波形の発生時刻をｔ' _ｋ 、各波形の大きさをｃ _ｔ'ｋ とし、ｔ＜０の範囲で０をとりｔ≧０の範囲で１をとる関数ｈ(ｔ)、複素周波数Ｆ _ｋ を用いてΣ _ｋ Σ _ｔ'k c _ｔ'k ｈ(ｔ-ｔ' _ｋ )ｅ ^{２πｉＦｋ(ｔ-t'k)} という形で表すことができる非定常的な時系列信号ｓ _ｌ を記憶した記憶手段から前記時系列信号を読み出して有限長の区間に分割するステップと、
   予測係数算出手段が、各区間ごとに前記時系列信号ｓ _ｌ を近似する時系列信号ｓ'_ｌをｓ'_ｌ＝Σ_ｋａ_ｋｓ_ｌ-ｋとしたときのΣ_ｌ‖ｓ_ｌ-ｓ'_ｌ‖が最小となる予測係数ａ_ｋを算出して記憶手段に記憶させるステップと、
   複素周波数算出手段が、前記算出された予測係数ａ_ｋを係数とする１-Σ_ｋａ_ｋｚ^-ｋ を、ｅ ^{２πｉＦｋΔＴ} ｚ ^-１ を根としたときの１次の因数の積に因数分解して得られる１-Σ_ｋａ_ｋｚ^-ｋ＝Π_ｋ（１-ｅ^{２πＦｋΔＴ}ｚ^-１）の関係を満たす複素周波数Ｆ_ｋを算出して記憶手段に記憶させるステップと、
   を有することを特徴とする周波数解析方法。
4. 分割手段が、各波形の発生時刻をｔ' _ｋ 、各波形の大きさをｃ _ｔ'ｋ とし、ｔ＜０の範囲で０をとりｔ≧０の範囲で１をとる関数ｈ(ｔ)、複素周波数Ｆ _ｋ を用いてΣ _ｋ Σ _ｔ'k c _ｔ'k ｈ(ｔ-ｔ' _ｋ )ｅ ^{２πｉＦｋ(ｔ-t'k)} という形で表すことができる非定常的な時系列信号ｓ _ｌ を記憶した記憶手段から前記時系列信号を読み出して有限長の区間に分割するステップと、
   Ｚ変換算出手段が、各区間ごとに前記時系列信号ｓ _ｌ についてのＺ変換Σ_ｋｓ_ｋｚ^-ｋを算出して記憶手段に記憶させるステップと、
   複素周波数算出手段が、前記算出されたＺ変換Σ_ｋｓ_ｋｚ^-ｋを規格化したｃ_ｋを係数とする１+Σ_ｋｃ_ｋｚ^-ｋを用いてａ_ｋ及びｄ_ｋを係数とするΣ_ｋｄ_ｋ/(1-ａ_ｋｚ^-１)を算出し、Σ_ｋｄ_ｋ/(1-ａ_ｋｚ^-１)＝Π_ｋ(1-ｅ^{２πｉＦ"ｋΔＴ}ｚ^-１)/Π_ｋ(1-ｅ^{２πｉＦ'ｋΔＴ}ｚ^-１)の関係を満たすＫ'個の複素周波数Ｆ'_ｋ=ｆ'_ｋ＋ｉｂ'_ｋおよびＫ"個の複素周波数Ｆ"_ｋ=ｆ"_ｋ+ｉｂ"_ｋを算出して記憶手段に記憶させるステップと、
   を有することを特徴とする周波数解析方法。
5. スペクトル包絡算出手段が、前記算出された複素周波数Ｆ _ｋ を、次式Ｇ(f)=1/｜1-ｅ ^{２πｉ（Ｆｋ-ｆ）ΔＴ} ｜ ^２ に適用することによりスペクトル包絡線Ｇ(f)を算出して記憶手段に記憶させるステップを有することを特徴とする請求項１乃至４のいずれかに記載の周波数解析方法。
6. 各波形の発生時刻をｔ' _ｋ 、各波形の大きさをｃ _ｔ'ｋ とし、ｔ＜０の範囲で０をとりｔ≧０の範囲で１をとる関数ｈ(ｔ)、複素周波数Ｆ _ｋ を用いてΣ _ｋ Σ _ｔ'k c _ｔ'k ｈ(ｔ-ｔ' _ｋ )ｅ ^{２πｉＦｋ(ｔ-t'k)} という形で表すことができる非定常的な時系列信号ｓ _ｌ を記憶する記憶手段と、
   前記記憶手段から前記時系列信号を読み出して有限長の区間に分割する分割手段と、
   各区間ごとに前記時系列信号を用いて次の連立方程式
   ｓ _ｌ ＝Σ _ｋ ａ _ｋ ｓ _ｌ-ｋ
   ｓ _ｌ-１ ＝Σ _ｋ ａ _ｋ ｓ _ｌ-１-ｋ
   …
   ｓ _ｌ-Ｋ ＝Σ _ｋ ａ _ｋ ｓ _ｌ-Ｋ-ｋ
   をａ _ｋ について解くことにより応答関数ａ_ｋを算出して記憶手段に記憶させる応答関数算出手段と、
   前記算出された応答関数ａ_ｋを係数とする次式１-Σ_ｋａ_ｋｚ^-ｋ を、ｅ ^{２πｉＦｋΔＴ} ｚ ^-１ を根としたときの１次の因数の積に因数分解して得られる１-Σ_ｋａ_ｋｚ^-ｋ＝Π_ｋ(１-ｅ^{２πｉＦｋΔＴ}ｚ^-１)の関係を満たす複素周波数Ｆ_ｋを算出して記憶手段に記憶させる複素周波数算出手段と、
   を有することを特徴とする周波数解析装置。
7. 各波形の発生時刻をｔ' _ｋ 、各波形の大きさをｃ _ｔ'ｋ とし、ｔ＜０の範囲で０をとりｔ≧０の範囲で１をとる関数ｈ(ｔ)、複素周波数Ｆ _ｋ を用いてΣ _ｋ Σ _ｔ'k c _ｔ'k ｈ(ｔ-ｔ' _ｋ )ｅ ^{２πｉＦｋ(ｔ-t'k)} という形で表すことができる非定常的な時系列信号ｓ _ｌ を記憶する記憶手段と、
   前記記憶手段から時系列信号を読み出して有限長の区間に分割する分割手段と、
   各区間ごとに多項式ｙ(ｔ)＝Σ_ｍｃ_ｍ(ｔ-ｔ_０)^ｍを用いたときのΣ_ｌ‖ｓ_ｌ-ｙ(ｔ)‖が最小となるように未定係数ｃ_ｍを算出し、この未定係数ｃ_ｍとａ _ｋ を根とする微分演算子Ｄ＝Π_ｋ(d/dt-ａ_ｋ)とを用いた∫_ｔ‖Ｄ・ｙ(t)‖dtについて、ａ _ｋ で偏微分したときの∂/∂ａ _ｋ ∫ _ｔ ‖Ｄ・ｙ(t)‖dt＝０を満たす未定係数ａ_ｋを算出して記憶手段に記憶させる微分係数算出手段と、
   算出されたａ_ｋを用いてａ_ｋ＝２πｉＦ_ｋの関係を満たす複素周波数Ｆ_ｋを算出して記憶手段に記憶させる複素周波数算出手段と、
   を有することを特徴とする周波数解析装置。
8. 各波形の発生時刻をｔ' _ｋ 、各波形の大きさをｃ _ｔ'ｋ とし、ｔ＜０の範囲で０をとりｔ≧０の範囲で１をとる関数ｈ(ｔ)、複素周波数Ｆ _ｋ を用いてΣ _ｋ Σ _ｔ'k c _ｔ'k ｈ(ｔ-ｔ' _ｋ )ｅ ^{２πｉＦｋ(ｔ-t'k)} という形で表すことができる非定常的な時系列信号ｓ _ｌ を記憶する記憶手段と、
   前記記憶手段から時系列信号を読み出して有限長の区間に分割する分割手段と、
   各区間ごとに前記時系列信号ｓ _ｌ を近似する時系列信号ｓ'_ｌをｓ'_ｌ＝Σ_ｋａ_ｋｓ_ｌ−ｋとしたときのΣ_ｌ‖ｓ_ｌ−ｓ'_ｌ‖が最小となる予測係数ａ_ｋを算出して記憶手段に記憶させる予測係数算出手段と、
   前記算出された予測係数ａ_ｋを係数とする１-Σ_ｋａ_ｋｚ^-ｋ を、ｅ ^{２πｉＦｋΔＴ} ｚ ^-１ を根としたときの１次の因数の積に因数分解して得られる１-Σ_ｋａ_ｋｚ^-ｋ＝Π_ｋ（１-ｅ^{２πｉＦｋΔＴ}ｚ^-１）の関係を満たす複素周波数Ｆ_ｋを算出して記憶手段に記憶させる複素周波数算出手段と、
   を有することを特徴とする周波数解析装置。
9. 各波形の発生時刻をｔ' _ｋ 、各波形の大きさをｃ _ｔ'ｋ とし、ｔ＜０の範囲で０をとりｔ≧０の範囲で１をとる関数ｈ(ｔ)、複素周波数Ｆ _ｋ を用いてΣ _ｋ Σ _ｔ'k c _ｔ'k ｈ(ｔ-ｔ' _ｋ )ｅ ^{２πｉＦｋ(ｔ-t'k)} という形で表すことができる非定常的な時系列信号ｓ _ｌ を記憶する記憶手段と、
   前記記憶手段から時系列信号を読み出して有限長の区間に分割する分割手段と、
   各区間ごとに前記時系列信号ｓ _ｌ についてのＺ変換Σ_ｋｓ_ｋｚ^-ｋを算出して記憶手段に記憶させるＺ変換算出手段と、
   前記算出されたＺ変換Σ_ｋｓ_ｋｚ^-ｋを規格化したｃ_ｋを係数とする１+Σ_ｋｃ_ｋｚ^-ｋを用いてａ_ｋ及びｄ_ｋを係数とするΣ_ｋｄ_ｋ/(1-ａ_ｋｚ^-１)を算出し、Σ_ｋｄ_ｋ/(1-ａ_ｋｚ^-１)＝Π_ｋ(1-ｅ^{２πｉＦ"ｋΔＴ}ｚ^-１)/Π_ｋ(1-ｅ^{２πｉＦ'ｋΔＴ}ｚ^-１)の関係を満たすＫ'個の複素周波数Ｆ'_ｋ=ｆ'_ｋ＋ｉｂ'_ｋおよびＫ"個の複素周波数Ｆ"_ｋ=ｆ"_ｋ+ｉｂ"_ｋを算出して記憶手段に記憶させる複素周波数算出手段と、
   を有することを特徴とする周波数解析装置。
10. 前記算出された複素周波数Ｆ _ｋ を、次式Ｇ(f)=1/｜1-ｅ ^{２πｉ（Ｆｋ-ｆ）ΔＴ} ｜ ^２ に適用することによりスペクトル包絡線Ｇ(f)を算出して記憶手段に記憶させるスペクトル包絡算出手段を備えることを特徴とする請求項６乃至９のいずれかに記載の周波数解析装置。
11. 請求項６乃至１０のいずれかに記載の周波数解析装置における各手段をコンピュータに機能させるためのコンピュータプログラム。

Images