JP6640702B2 - 時系列信号特徴推定装置、プログラム - Google Patents

Description

本発明は、時系列信号特徴推定技術に関し、特に短時間フーリエ変換に基づく時系列信号の特徴を推定する技術に関する。

音声信号や音響信号などの時系列信号を分析する手法の１つとして、短時間フーリエ変換がある（非特許文献１参照）。分析対象となる時系列信号は、短時間フーリエ変換された結果、時間と角周波数（時間と周波数）のインデックスを持つ2次元の複素数信号に変換される。

この複素数信号の実数部と虚数部をそのまま分析に用いる代わりに、振幅成分と位相成分に変換し、これらの両方の成分の特徴またはいずれか一方の成分のみの特徴に注目して、分析を行うことがある。また、位相成分の角周波数についての微分の正負を反転させた値である群遅延特性や、位相成分の時間についての微分の値である瞬時角周波数（瞬時角周波数を2πで除した値である瞬時周波数）に基づき、分析を行うこともある。

M. Parchami, W. P. Zhu, B. Champagne, E. Plourde, "Recent Developments in Speech Enhancement in the Short-Time Fourier Transform Domain", IEEE Circuits and Systems Magazine, Vol.16, No.3, pp.45-77, 2016.

２つの成分のうち、振幅成分は、得られる特徴の物理的な解釈のしやすさや数値的な計算の簡潔さから、採用されることが多い。

一方、位相成分や、その角周波数微分や時間微分から得られる群遅延特性や瞬時角周波数（瞬時周波数）についても、有益な特徴が得られることが期待されるものの、複素数信号から直接的に計算される位相成分は、その値が-πからπの間にラップされた状態となるために、時間軸方向にも周波数軸方向にも不連続点を有する2次元の系列となってしまう。このため、群遅延特性や瞬時角周波数の値を算出するためには、位相成分の不連続を解消しアンラップされた位相成分の特性を知る必要があるが、一般に位相成分の特性のアンラップ処理は一意な解を持たず不確定な問題となる。

そこで本発明は、短時間フーリエ変換に基づく時系列信号の特徴分析において、不確定性が少ない形で、時系列信号特徴を算出することができる時系列信号特徴推定技術を提供することを目的とする。

本発明の一態様は、窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間t_n、離散角周波数ω_k(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を

(ただし、f_sは標本化周波数、Nは離散フーリエ変換の点数である）、離散時間t_n、離散角周波数ω_kにおける振幅をA_n,kとし、前記振幅A_n,k(n₁,n₂はn₂-n₁≧1を満たす整数、n=n₁, n₁+1,…, n₂、kは整数)から、離散角周波数ω_kについて離散時間t_nにおける対数振幅logA_n,kの時間差分である対数振幅時間差分Δ_tlogA_n,k/Δt_nを算出する対数振幅時間差分算出部と、前記対数振幅時間差分Δ_tlogA_n,k/Δt_nから、次式を用いて離散時間t_n、離散角周波数ω_kにおける群遅延特性GD_n,kを時系列信号特徴として算出する時系列信号特徴算出部と

(ただし、Tは窓関数の有効幅であり、T=N/f_s)を含む。

本発明の一態様は、窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間t_n、離散角周波数ω_k(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を

(ただし、f_sは標本化周波数、Nは離散フーリエ変換の点数である）、離散時間t_n、離散角周波数ω_kにおける振幅をA_n,kとし、前記振幅A_n,k(nは整数、k₁,k₂はk₂-k₁≧1を満たす整数、k=k₁, k₁+1,…, k₂)から、離散時間t_nについて離散角周波数ω_kにおける対数振幅logA_n,kの角周波数差分である対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kを算出する対数振幅角周波数差分算出部と、前記対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kから、次式を用いて離散時間t_n、離散角周波数ω_kにおける瞬時周波数IF_n,kを時系列信号特徴として算出する時系列信号特徴算出部と

を含む。

本発明の一態様は、窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間t_n、離散角周波数ω_k(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を

(ただし、f_sは標本化周波数、Nは離散フーリエ変換の点数である）、離散時間t_n、離散角周波数ω_kにおける振幅をA_n,kとし、前記振幅A_n,k(n₁,n₂はn₂-n₁≧1を満たす整数、n=n₁, n₁+1,…, n₂、k=0, 1,…, N-1)から、k=0,1,…, N-2の各kについて離散角周波数ω_kからω_k+1にかけての群遅延特性の符号反転値の積分値である群遅延特性積分値Σ(-GD)_n,kを算出する群遅延特性積分値算出部と、k₀を0以上N-1以下のある整数、位相の初期値をφ_{n,k_0}とし、前記群遅延特性積分値Σ(-GD)_n,k(k=k₀, k₀+1,…, N-2)から、次式を用いて、k=k₀+1から昇順に位相φ_n,k(k=k₀+1,…, N-1)を時系列信号特徴として算出し、

前記群遅延特性積分値Σ(-GD)_n,k(k=0, 1,…,k₀-1)から、次式を用いて、k=k₀-1から降順に位相φ_n,k(k=0, 1,…,k₀-1)を時系列信号特徴として算出する時系列信号特徴算出部と

を含む。

本発明の一態様は、窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間t_n、離散角周波数ω_k(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を

(ただし、f_sは標本化周波数、Nは離散フーリエ変換の点数である）、離散時間t_n、離散角周波数ω_kにおける振幅をA_n,kとし、前記振幅A_n,k(n₀を整数、mを正の整数、n=n₀, n₀+1,…, n₀+m、k₁,k₂はk₂-k₁≧1を満たす整数、k=k₁, k₁+1,…, k₂)から、n=n₀, n₀+1,…, n₀+m-1の各nについて離散時間t_nからt_n+1にかけての瞬時周波数の2π倍の積分値である瞬時周波数積分値Σ(2πIF)_n,kを算出する瞬時周波数算出部と、位相の初期値をφ_{n_0,k}とし、前記瞬時周波数積分値Σ(2πIF)_n,k(n=n₀, n₀+1,…, n₀+m-1)から、次式を用いて、n=n₀+1から昇順に位相φ_n,k(n=n₀+1,…, n₀+m)を時系列信号特徴として算出する時系列信号特徴算出部と

を含む。

本発明によれば、時系列信号の短時間フーリエ変換の振幅成分に基づき、位相成分、群遅延特性、瞬時周波数を算出することにより、これらの値を不確定性の少ない形で算出することが可能となる。

時系列信号特徴推定装置１００の構成の一例を示す図。 時系列信号特徴推定装置１００の動作の一例を示す図。 時系列信号特徴推定装置２００の構成の一例を示す図。 時系列信号特徴推定装置２００の動作の一例を示す図。 時系列信号特徴推定装置３００の構成の一例を示す図。 時系列信号特徴推定装置３００の動作の一例を示す図。 時系列信号特徴推定装置４００の構成の一例を示す図。 時系列信号特徴推定装置４００の動作の一例を示す図。 離散角周波数インデックスk=0における真の位相と群遅延特性の推定値を示す図。

以下、本発明の実施の形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。

＜表記方法＞
_（アンダースコア）は下付き添字を表す。例えば、x^y_zはy_zがxに対する上付き添字であり、x_{y_z}はy_zがxに対する下付き添字であることを表す。

＜短時間フーリエ変換に基づく時系列信号の特徴推定の原理＞
まず、短時間フーリエ変換に基づく時系列信号の特徴推定の原理となる群遅延特性、瞬時角周波数、位相が満たす方程式について説明する。これらの方程式を用いることにより、時系列信号の短時間フーリエ変換により得られる2次元の複素数信号の位相成分を、複素数信号の値そのものから直接計算する代わりに、複素数信号の振幅成分から推定することが可能となる。また、群遅延特性、瞬時角周波数も振幅成分から推定することが可能となる。

［群遅延特性、瞬時角周波数、位相が満たす方程式］
まず、時系列信号である連続時間信号x(τ)について、式(1)で与えられる時間tと角周波数ωにおける短時間フーリエ変換X_t,ωを考える。

ただし、w(t)はt=0において最大値を取る窓関数、Tは窓関数の有効幅とする。なお、jは虚数単位である。

このとき、短時間フーリエ変換X_t,ωを時間tで偏微分すると式(2)が得られる。

同様に、短時間フーリエ変換X_t,ωを角周波数ωで偏微分すると式(3)が得られる。

ここで、窓関数として式(4)で与えられるガウス窓を選択する。

ただし、ガウス窓の標準偏差σは0より大きい実数(正の実数)とする。

このとき、窓関数w(τ-t)を時間tで偏微分すると式(5)が得られる。

このとき、式(4)を式(3)に、式(5)を式(2)に適用し、共通項を消去すると、式(6)が得られる。

短時間フーリエ変換X_t,ωの振幅をA_t,ω、位相をφ_t,ωとすると、式(7)が成り立つ。また、式(7)の対数をとると、式(8)が得られる。

ここで、logX_t,ωを時間tで偏微分すると式(9)が得られる。

したがって、式(8)、式(9)より、式(10)が成立する。

同様に、logX_t,ωを角周波数ωで偏微分すると式(11)が得られる。

したがって、式(8)、式(11)より、式(12)が成立する。

式(10)と式(12)を式(6)に代入すると、式(13)が得られる。

ここで、式(13)の左辺の実数部と虚数部は独立して0となることから、式(14)と式(15)が得られる。

式(14)の左辺は群遅延特性の定義そのものであり、式(14)の右辺から群遅延特性が位相を用いることなく振幅情報である対数振幅のみから計算できることが分かる。また、式(15)の左辺は瞬時角周波数の定義そのものであり、式(15)の右辺から瞬時角周波数も位相を用いることなく対数振幅のみから計算できることが分かる。瞬時周波数は瞬時角周波数を2πで除することで求められることから、瞬時周波数も対数振幅のみから計算できることが分かる。

さらに、式(14)を角周波数について積分すると式(16)が得られる。

ただし、C₁は積分定数である。

同様に、式(15)を時間について積分すると、式(17)が得られる。

ただし、C₂は積分定数である。

なお、積分定数C₁、C₂は別途決定する必要がある。

式(16)、式(17)の左辺は位相φ_t,ωであり、位相φ_t,ωもまた対数振幅のみから計算できることが分かる。

以上より、短時間フーリエ変換の窓関数としてガウス窓を選択することで、群遅延特性は式(14)を、瞬時角周波数（瞬時周波数）は式(15)を、位相は式(16)または式(17)を用いてそれぞれ対数振幅から計算できることが分かる。

群遅延特性や瞬時角周波数（瞬時周波数）は、ラップされた位相成分の代わりに振幅成分を用いて計算することで、不確定性の問題の影響を受けにくくなる。また、位相そのものの計算についても、複素数信号から直接計算する場合にはラップされた状態でしか得られないが、式(16)または式(17)を用いることで、物理的特徴の理解しやすいアンラップされた状態で計算することができる。

以上の短時間フーリエ変換に基づく時系列信号の特徴推定の原理となる群遅延特性、瞬時角周波数、位相が満たす方程式は、連続時間信号についての短時間フーリエ変換について説明したものである。しかし、実際にコンピュータで計算する場合には、離散化されたデジタル信号についての短時間フーリエ変換の結果を利用する必要がある。

そこで、以下では、連続値を取る時間tと角周波数ωをそれぞれ離散化し、群遅延特性、瞬時周波数、位相を推定するための近似式について説明する。

［群遅延特性、瞬時周波数、位相を推定するための近似式］
時系列信号である連続時間信号x(τ)の、標準偏差がσ(σは正の実数)であるガウス窓を窓関数とする短時間フーリエ変換X_t,ωに対して、離散フーリエ変換を適用して得られる離散短時間フーリエ変換X_n,k(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)、つまり、窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換X_n,kを考える。

これにより、時間tを離散化した離散時間t_n、角周波数ωを離散化した離散角周波数ω_kは、それぞれ式(18)、式(19)で表される。

ただし、f_sは標本化周波数、Nは離散フーリエ変換(DFT; Discrete Fourier Transform)の点数である。

また、このとき、窓関数の有効幅Tは、式(20)で表される。

離散化されたデジタル信号x(τ_i)について、離散時間t_n、離散角周波数ω_kにおける離散短時間フーリエ変換X_n,kは、式(21)で表される。

式(21)は、式(1)に対応するものである。

（群遅延特性を推定するための近似式）
式(14)を離散化することにより、群遅延特性を推定するための近似式を求める。

離散時間t_n、離散角周波数ω_kにおける振幅をA_n,kとする。また、離散角周波数ω_kについて離散時間t_nにおける対数振幅logA_n,kの時間差分である対数振幅時間差分をΔ_tlogA_n,k/Δt_nとする。

離散角周波数ω_kについて離散時間t_nにおける対数振幅時間差分Δ_tlogA_n,k/Δt_nとして中心差分近似を用いると、式(22)を得る。

ただし、δ₁、δ₂は数値安定化のために導入した非零の数とする。

式(14)と式(22)から、離散時間t_n、離散角周波数ω_kにおける群遅延特性GD_n,kを推定するための近似式である式(23)を得る。

また、式(24)の関係から、中心差分近似を用いた対数振幅時間差分Δ_tlogA_n,k/Δt_nとして式(22)の代わりに式(25)を得る。

ただし、δ₃は数値安定化のために導入した非零の数とする。なお、Δ_tA_n,k/Δt_nは離散角周波数ω_kについて離散時間t_nにおける振幅A_n,kの時間差分である振幅時間差分である。

式(14)と式(25)から、離散時間t_n、離散角周波数ω_kにおける群遅延特性GD_n,kを推定するための近似式である式(26)を得る。

なお、ここでは、中心差分近似を用いて対数振幅時間差分Δ_tlogA_n,k/Δt_nを算出したが、中心差分近似の代わりに、前進差分近似や後進差分近似を用いてもよい。例えば、式(22)に対応する前進差分近似、後進差分近似はそれぞれ式(22a)、式(22b)のようになる。

同様に、式(25)に対応する前進差分近似、後進差分近似はそれぞれ式(25a)、式(25b)のようになる。

ここで、ガウス窓の標準偏差σについて、σのα倍(ただし、αは１以上の整数とする)が窓関数の有効幅の1/2に相当すると仮定すると、ガウス窓の標準偏差σは式(27)を満たす。

例えば、α=4とすると、このαに対応した標準偏差σは、式(28)で求めることができる。

上記仮定を導入すると、選択したσの妥当性を直観的に解釈しやすくなるというメリットがある。

（瞬時周波数を推定するための近似式）
式(15)を離散化し2πで除することにより、瞬時周波数を推定するための近似式を求める。

離散時間t_nについて離散角周波数ω_kにおける対数振幅logA_n,kの角周波数差分である対数振幅角周波数差分をΔ_ωlogA_n,k/Δω_kとする。

離散時間t_nについて離散角周波数ω_kにおける対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kとして中心差分近似を用いると、式(29)を得る。

ただし、δ₄、δ₅は数値安定化のために導入した非零の数とする。

式(15)と式(29)から、離散時間t_n、離散角周波数ω_kにおける瞬時周波数IF_n,kを推定するための近似式である式(30)を得る。

また、式(31)の関係から、中心差分近似を用いた対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kとして式(29)の代わりに式(32)を得る。

ただし、δ₆は数値安定化のために導入した非零の数とする。なお、Δ_ωA_n,k/Δω_kは離散時間t_nについて離散角周波数ω_kにおける振幅A_n,kの角周波数差分である振幅角周波数差分である。

式(15)と式(32)から、離散時間t_n、離散角周波数ω_kにおける瞬時周波数IF_n,kを推定するための近似式である式(33)を得る。

なお、ここでは、中心差分近似を用いて対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kを算出したが、中心差分近似の代わりに、前進差分近似や後進差分近似を用いてもよい。例えば、式(29)に対応する前進差分近似、後進差分近似はそれぞれ式(29a)、式(29b)のようになる。

同様に、式(32)に対応する前進差分近似、後進差分近似はそれぞれ式(32a)、式(32b)のようになる。

（位相を推定するための近似式）
位相を推定するための近似式を求めるため、以下では、式(23)、式(26)、式(30)、式(33)について台形近似により積分計算をし、再帰式の形式にて求める。

式(23)、式(26)については、周波数方向に積分近似する。周波数方向の台形近似による積分計算の式は式(34)のようになる。

式(34)は、離散時間t_n、離散角周波数ω_kにおける位相φ_n,kが、離散時間t_n、離散角周波数ω_k-1における位相φ_n,k-1に、離散角周波数ω_k-1からω_kにかけての群遅延特性の符号反転値（式(14)参照）の積分値を加えたものとなることを示している。式(34)の右辺の第二項が台形の面積に相当する。

離散時間t_n、離散角周波数ω_kにおける位相φ_n,kは、式(34)に式(23)を適用することで、式(35)の近似式で表される。

ただし、δ₇、δ₈は数値安定化のために導入した非零の数とする。

同様に、離散時間t_n、離散角周波数ω_kにおける位相φ_n,kは、式(34)に式(26)を適用することで、式(36)の近似式で表される。

また、式(34)から式(37)を得る。

式(37)は、離散時間t_n、離散角周波数ω_kにおける位相φ_n,kが、離散時間t_n、離散角周波数ω_k+1における位相φ_n,k+1から、離散角周波数ω_kからω_k+1にかけての群遅延特性の符号反転値（式(14)参照）の積分値を差し引いたものとなることを示している。式(37)の右辺の第二項が台形の面積に相当する。式(37)を用いると、周波数方向について負の方向に計算を進めていく形で位相を求めることができることが分かる。

離散時間t_n、離散角周波数ω_kにおける位相φ_n,kはを求める近似式として、式(37)に対応する形で式(35)、式(36)からそれぞれ式(38)、式(39)を得る。

次に、式(30)、式(33)については、時間方向に積分近似する。時間方向の台形近似による積分計算の式は式(40)のようになる。

式(40)は、離散時間t_n、離散角周波数ω_kにおける位相φ_n,kが、離散時間t_n-1、離散角周波数ω_kにおける位相φ_n-1,kに、離散時間t_n-1からt_nにかけての瞬時周波数の2π倍（式(15)参照）の積分値を加えたものとなることを示している。式(40)の右辺の第二項が台形の面積に相当する。

離散時間t_n、離散角周波数ω_kにおける位相φ_n,kは、式(40)に式(30)を適用することで、式(41)の近似式で表される。

ただし、δ₉、δ₁₀は数値安定化のために導入した非零の数とする。

同様に、離散時間t_n、離散角周波数ω_kにおける位相φ_n,kは、式(40)に式(33)を適用することで、式(42)の近似式で表される。

式(35)と式(36)（式(38)と式(39)）、式(41)と式(42)の位相を推定するための再帰的近似式は、台形近似による積分計算により得られたものであるが、台形近似以外の積分計算方法を用いてもよい。例えば、シンプソンの公式など、一般にニュートン・コーツの公式として知られる、等間隔に離散化された被積分関数の数値積分公式を用いて、式(35)と式(36)（式(38)と式(39)）、式(41)と式(42)のような位相を推定するための再帰的近似式を求めてもよい。

一般に、離散角周波数ω_k-1からω_kにかけての群遅延特性の符号反転値の積分値（以下、群遅延特性の符号反転値の積分値を単に群遅延特性積分値という）をΣ(-GD)_n,k-1とすると、式(34)に相当する式として式(43)が得られる。ここで、群遅延特性積分値Σ(-GD)_n,k-1は、群遅延特性GD_n,k-1と群遅延特性GD_n,kを用いて計算される値である。

なお、式(23)や式(26)（あるいはその導出に用いる式(22)、式(22a)、式(22b)や式(25)、式(25a)、式(25b)）を見ればわかるように、群遅延特性GD_n,kは、一般に振幅A_n-1,k、振幅A_n,k、振幅A_n+1,kのうち、少なくとも2つの値を用いて計算することができる。したがって、群遅延特性積分値Σ(-GD)_n,k-1は、振幅A_n-1,k-1、振幅A_n,k-1、振幅A_n+1,k-1のうち、少なくとも2つの値と、振幅A_n-1,k、振幅A_n,k、振幅A_n+1,kのうち、少なくとも2つの値を用いて計算される値でもある。

同様に、離散角周波数ω_kからω_k+1にかけての群遅延特性の符号反転値の積分値（群遅延特性積分値）をΣ(-GD)_n,kとすると、式(37)に相当する式として式(44)が得られる。ここで、群遅延特性積分値Σ(-GD)_n,kは、群遅延特性GD_n,kと群遅延特性GD_n,k+1を用いて計算される値である。

なお、群遅延特性積分値Σ(-GD)_n,kは、振幅A_n-1,k、振幅A_n,k、振幅A_n+1,kのうち、少なくとも2つの値と、振幅A_n-1,k+1、振幅A_n,k+1、振幅A_n+1,k+1のうち、少なくとも2つの値を用いて計算される値でもある。

また、離散時間t_n-1からt_nにかけての瞬時周波数の2π倍の積分値（以下、瞬時周波数の2π倍の積分値を単に瞬時周波数積分値という）をΣ(2πIF)_n-1,kとすると、式(40)に相当する式として式(45)が得られる。ここで、瞬時周波数積分値Σ(2πIF)_n-1,kは、瞬時周波数IF_n-1,kと瞬時周波数IF_n,kを用いて計算される値である。

なお、式(30)や式(33)（あるいはその導出に用いる式(29)、式(29a)、式(29b)や式(32)、式(32a)、式(32b)）を見ればわかるように、瞬時周波数IF_n,kは、一般に振幅A_n,k-1、振幅A_n,k、振幅A_n,k+1のうち、少なくとも2つの値を用いて計算することができる。したがって、瞬時周波数積分値Σ(2πIF)_n-1,kは、振幅A_n-1,k-1、振幅A_n-1,k、振幅A_n-1,k+1のうち、少なくとも2つの値と、振幅A_n,k-1、振幅A_n,k、振幅A_n,k+1のうち、少なくとも2つの値を用いて計算される値でもある。

＜第一実施形態＞
以下、図１〜図２を参照して時系列信号特徴推定装置１００について説明する。図１に示すように時系列信号特徴推定装置１００は、対数振幅時間差分算出部１１０、時系列信号特徴算出部１２０、記録部１９０を含む。記録部１９０は、時系列信号特徴推定装置１００の処理に必要な情報を適宜記録する構成部である。例えば、記録部１９０には、標本化周波数f_s、離散フーリエ変換の点数N、ガウス窓の標準偏差σを事前に記録しておく。もちろん、標本化周波数f_s、離散フーリエ変換の点数N、ガウス窓の標準偏差σを事前に記録部１９０に記録しておく代わりに、時系列信号特徴推定装置１００の入力として与えるのでもよい。

時系列信号特徴推定装置１００の入力となる振幅A_n,k(n₁,n₂はn₂-n₁≧1を満たす整数、n=n₁, n₁+1,…, n₂、kは整数)は、事前に得られているものとする。例えば、メモリや伝送容量を削減するため、事前に時系列信号を短時間フーリエ変換した結果のうち、振幅A_n,kのみを記録部１９０に記録しておいてもよい。また、時系列信号特徴推定装置１００による処理開始前に外部から振幅A_n,kが伝送されてくるのでもよい。

時系列信号特徴推定装置１００は、振幅A_n,k(n=n₁, n₁+1,…,n₂、kはある整数)から、離散時間t_n、離散角周波数ω_kにおける時系列信号特徴である群遅延特性GD_n,kを推定し、出力する。

図２に従い時系列信号特徴推定装置１００の動作について説明する。対数振幅時間差分算出部１１０は、振幅A_n,k(n=n₁, n₁+1,…,n₂、kはある整数)から、離散角周波数ω_kについて離散時間t_nにおける対数振幅logA_n,kの時間差分である対数振幅時間差分Δ_tlogA_n,k/Δt_nを算出する（Ｓ１１０）。対数振幅時間差分Δ_tlogA_n,k/Δt_nは、n₁<n<n₂を満たすnについては中心差分近似、端点であるn=n₁またはn=n₂については前進差分近似または後進差分近似により算出すればよい。もちろん、n₁<n<n₂を満たすnについて、中心差分近似の代わりに前進差分近似や後進差分近似を用いてもよい。対数振幅時間差分Δ_tlogA_n,k/Δt_nを、n₁<n<n₂を満たすnについては中心差分近似、n=n₁については前進差分近似、n=n₂については後進差分近似により算出する場合、例えば、式(22)、式(22b)、式(22b)を用いて算出することができる。

時系列信号特徴算出部１２０は、Ｓ１１０で算出した対数振幅時間差分Δ_tlogA_n,k/Δt_nから、離散時間t_n、離散角周波数ω_kにおける群遅延特性GD_n,kを時系列信号特徴として算出する（Ｓ１２０）。具体的には、ガウス窓の標準偏差σ、窓関数の有効幅T=N/f_sを用いて

と計算すればよい。例えば、式(23)を用いて、n₁<n<n₂を満たすnについては中心差分近似、n=n₁については前進差分近似、n=n₂については後進差分近似により対数振幅時間差分Δ_tlogA_n,k/Δt_nを計算する場合は、それぞれ

となる。また、式(26)を用いて、n₁<n<n₂を満たすnについては中心差分近似、n=n₁については前進差分近似、n=n₂については後進差分近似により対数振幅時間差分Δ_tlogA_n,k/Δt_nを計算する場合は、それぞれ

となる。

本実施形態の発明によれば、短時間フーリエ変換による信号分析において、群遅延特性の推定値を本来の定義である位相を用いて計算する代わりに、安定的に計算可能な振幅情報（対数振幅や振幅）を用いて計算することで、位相のアンラップ処理の不確定性の問題を回避することができる。

また、群遅延特性を振幅情報のみから推定できるため、短時間フーリエ変換の結果として得られる複素数信号の情報のうち、位相情報が欠損し、振幅情報のみが与えられた状況においても、完全な複素数信号が与えられた場合と同等の群遅延特性の推定が可能となる。

＜第二実施形態＞
以下、図３〜図４を参照して時系列信号特徴推定装置２００について説明する。図３に示すように時系列信号特徴推定装置２００は、対数振幅角周波数差分算出部２１０、時系列信号特徴算出部２２０、記録部１９０を含む。記録部１９０は、時系列信号特徴推定装置２００の処理に必要な情報を適宜記録する構成部である。

時系列信号特徴推定装置２００の入力となる振幅A_n,k(nは整数、k₁,k₂はk₂-k₁≧1を満たす整数、k=k₁, k₁+1,…, k₂)は、事前に得られているものとする。

時系列信号特徴推定装置２００は、振幅A_n,k(nはある整数、k=k₁, k₁+1,…, k₂)から、離散時間t_n、離散角周波数ω_kにおける時系列信号特徴である瞬時周波数IF_n,kを推定し、出力する。

図４に従い時系列信号特徴推定装置２００の動作について説明する。対数振幅角周波数差分算出部２１０は、振幅A_n,k(nはある整数、k=k₁, k₁+1,…, k₂)から、離散時間t_nについて離散角周波数ω_kにおける対数振幅logA_n,kの角周波数差分である対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kを算出する（Ｓ２１０）。対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kは、k₁<k<k₂を満たすkについては中心差分近似、端点であるk=k₁またはk=k₂については前進差分近似または後進差分近似により算出すればよい。もちろん、k₁<k<k₂を満たすkについて、中心差分近似の代わりに前進差分近似や後進差分近似を用いてもよい。対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kを、k₁<k<k₂を満たすkについては中心差分近似、k=k₁については前進差分近似、k=k₂については後進差分近似により算出する場合、例えば、式(29)、式(29b)、式(29b) を用いて算出することができる。

時系列信号特徴算出部２２０は、Ｓ２１０で算出した対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kから、離散時間t_n、離散角周波数ω_kにおける瞬時周波数IF_n,kを時系列信号特徴として算出する（Ｓ２２０）。具体的には、ガウス窓の標準偏差σを用いて

と計算すればよい。例えば、式(30)を用いて、k₁<k<k₂を満たすkについては中心差分近似、k=k₁については前進差分近似、k=k₂については後進差分近似により対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kを計算する場合は、それぞれ

となる。また、式(33)を用いて、k₁<k<k₂を満たすkについては中心差分近似、k=k₁については前進差分近似、k=k₂については後進差分近似により対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kを計算する場合は、それぞれ

となる。

本実施形態の発明によれば、短時間フーリエ変換による信号分析において、瞬時周波数の推定値を本来の定義である位相を用いて計算する代わりに、安定的に計算可能な振幅情報（対数振幅や振幅）を用いて計算することで、位相のアンラップ処理の不確定性の問題を回避することができる。

また、瞬時周波数を振幅情報のみから推定できるため、短時間フーリエ変換の結果として得られる複素数信号の情報のうち、位相情報が欠損し、振幅情報のみが与えられた状況においても、完全な複素数信号が与えられた場合と同等の瞬時周波数の推定が可能となる。

＜第三実施形態＞
以下、図５〜図６を参照して時系列信号特徴推定装置３００について説明する。図５に示すように時系列信号特徴推定装置３００は、群遅延特性積分値算出部３１０、時系列信号特徴算出部３２０、記録部１９０を含む。記録部１９０は、時系列信号特徴推定装置３００の処理に必要な情報を適宜記録する構成部である。

時系列信号特徴推定装置３００の入力となる振幅A_n,k(n₁,n₂はn₂-n₁≧1を満たす整数、n=n₁, n₁+1,…, n₂、k=0, 1,…, N-1)は、事前に得られているものとする。

時系列信号特徴推定装置３００は、振幅A_n,k(n=n₁, n₁+1,…, n₂、k=0, 1,…, N-1)から、離散時間t_n、離散角周波数ω_kにおける時系列信号特徴である位相φ_n,kを推定し、出力する。

図６に従い時系列信号特徴推定装置３００の動作について説明する。群遅延特性積分値算出部３１０は、振幅A_n,k(n=n₁, n₁+1,…,n₂、k=0, 1,…, N-1)から、k=0,1,…, N-2の各kについて離散角周波数ω_kからω_k+1にかけての群遅延特性の符号反転値の積分値である群遅延特性積分値Σ(-GD)_n,kを算出する（Ｓ３１０）。例えば、群遅延特性積分値算出部３１０は、時系列信号特徴推定装置１００を用いて、振幅A_n,k(n=n₁, n₁+1,…,n₂、k=0, 1,…, N-1)から、k=0,1,…, N-1の各kについて群遅延特性GD_n,kを算出したうえで、群遅延特性積分値Σ(-GD)_n,kを算出するのでよい。

また、別の例として、群遅延特性積分値算出部３１０は、n₁<n<n₂を満たすnについては、以下の式を用いて算出し、

n=n₁については、以下の式を用いて算出し、

n=n₂については、以下の式を用いて算出してもよい。

また、別の例として、群遅延特性積分値算出部３１０は、n₁<n<n₂を満たすnについては、以下の式を用いて算出し、

n=n₁については、以下の式を用いて算出し、

n=n₂については、以下の式を用いて算出してもよい。

時系列信号特徴算出部３２０は、位相の初期値をφ_{n,k_0}(k₀は0以上N-1以下のある整数) とＳ３１０で算出した群遅延特性積分値Σ(-GD)_n,k(k=0, 1,…, N-2)から、位相φ_n,k(k=0, 1,…, N-1、ただし、k=k₀を除く)を時系列信号特徴として算出する（Ｓ３２０）。整数k₀と位相の初期値φ_{n,k_0}は記録部１９０に事前に記録しておいてもよいし、時系列信号特徴推定装置３００の入力として与えるのでもよい。具体的には、以下のような手順で位相φ_n,k(k=0, 1,…, N-1、ただし、k=k₀を除く)を算出する。

まず、k=k₀+1とし、位相の初期値φ_{n,k_0}と群遅延特性積分値Σ(-GD)_n,k(k=k₀, k₀+1,…, N-2)を用いて次式（式(43)）より、位相φ_{n,k_0+1}、位相φ_{n,k_0+2}、…、位相φ_n,N-1と昇順に算出していく。

次に、k=k₀-1とし、位相の初期値φ_{n,k_0}と群遅延特性積分値Σ(-GD)_n,k(k=0, 1,…, k₀-1)を用いて次式（式(44)）より、位相φ_{n,k_0-1}、位相φ_{n,k_0-2}、…、位相φ_n,0と降順に算出していく。

時系列信号特徴算出部３２０による位相の算出には、位相の初期値φ_{n,k_0}を与える必要がある。各nについて、k₀=0のときの位相φ_n,0を初期値とすることができる。分析対象の信号が実数である場合、位相φ_n,0の値は0かπのいずれかに限定される。離散時間インデックスnの変化に従い、位相φ_n,0の値が0からπまたはπから0に変化するポイントにおいて、式(23)または式(26)より推定された離散角周波数インデックスk=0あるいはその周辺の離散角周波数インデックスにおける群遅延特性の推定値には図９(b)に示すような、大きなスパイク性の変化が観測される。このような変化を検出することで、位相φ_n,0の値を0かπかで切り替え、適切な初期値を与えることができる。

なお、時系列信号特徴算出部３２０で用いる再帰式の開始点は、k₀=0に限定されるものではないから、任意のk₀に対して初期値として適切な値を獲得できる場合には、それを初期値としてもよい。この場合、先ほど説明したように、kの値について昇順に再帰式を計算するだけでなく、降順に再帰式を計算することにより位相の推定値を得ることができる。

本実施形態の発明によれば、短時間フーリエ変換の結果として得られる複素数信号から位相を求める従来の方法では位相の値が-πからπの間にラップされた状態でしか計算できないのに対し、振幅情報を用いることで、アンラップされた状態での位相の値を推定することができる。

また、位相を振幅情報のみから推定できるため、短時間フーリエ変換の結果として得られる複素数信号の情報のうち、位相情報が欠損し、振幅情報のみが与えられた状況においても、完全な複素数信号が与えられた場合と同等の位相の推定が可能となる。

さらに、振幅情報のみが与えられた条件において位相の推定値を復元することで、一部の情報が欠損した中で複素数信号を再構成することが可能となる。

＜第四実施形態＞
以下、図７〜図８を参照して時系列信号特徴推定装置４００について説明する。図７に示すように時系列信号特徴推定装置４００は、瞬時周波数積分値算出部４１０、時系列信号特徴算出部４２０、記録部１９０を含む。記録部１９０は、時系列信号特徴推定装置４００の処理に必要な情報を適宜記録する構成部である。

時系列信号特徴推定装置４００の入力となる振幅A_n,k(n₀を整数、mを正の整数、n=n₀, n₀+1,…, n₀+m、k₁,k₂はk₂-k₁≧1を満たす整数、k=k₁, k₁+1,…, k₂)は、事前に得られているものとする。

時系列信号特徴推定装置４００は、振幅A_n,k(n=n₀, n₀+1,…, n₀+m、k=k₁, k₁+1,…, k₂)から、離散時間t_n、離散角周波数ω_kにおける時系列信号特徴である位相φ_n,kを推定し、出力する。

図８に従い時系列信号特徴推定装置４００の動作について説明する。瞬時周波数積分値算出部４１０は、振幅A_n,k(n=n₀, n₀+1,…, n₀+m、k=k₁, k₁+1,…, k₂)から、n=n₀, n₀+1,…, n₀+m-1の各nについて離散時間t_nからt_n+1にかけての瞬時周波数の2π倍の積分値である瞬時周波数積分値Σ(2πIF)_n,kを算出する（Ｓ４１０）。例えば、瞬時周波数積分値算出部４１０は、時系列信号特徴推定装置２００を用いて、振幅A_n,k(n=n₀, n₀+1,…, n₀+m、k=k₁, k₁+1,…, k₂)から、n=n₀, n₀+1,…, n₀+mの各nについて瞬時周波数IF_n,kを算出したうえで、周波数積分値Σ(2πIF)_n,kを算出するのでよい。

また、別の例として、瞬時周波数積分値算出部４１０は、k₁<k<k₂を満たすkについては、以下の式を用いて算出し、

k=k₁については、以下の式を用いて算出し、

k=k₂については、以下の式を用いて算出してもよい。

また、別の例として、瞬時周波数積分値算出部４１０は、k₁<k<k₂を満たすkについては、以下の式を用いて算出し、

k=k₁については、以下の式を用いて算出し、

k=k₂については、以下の式を用いて算出してもよい。

時系列信号特徴算出部４２０は、位相の初期値φ_{n_0,k}とＳ４１０で算出した瞬時周波数積分値Σ(2πIF)_n,k(n=n₀, n₀+1,…, n₀+m-1)から、位相φ_n,k(n=n₀+1,…, n₀+m)を時系列信号特徴として算出する（Ｓ４２０）。整数n₀と位相の初期値φ_{n_0,k}は記録部１９０に事前に記録しておいてもよいし、時系列信号特徴推定装置４００の入力として与えるのでもよい。具体的には、以下のような手順で位相φ_n,k(n=n₀+1,…, n₀+m)を算出する。

n=n₀+1とし、位相の初期値φ_{n_0,k}と瞬時周波数積分値Σ(2πIF)_n,k(n=n₀, n₀+1,…, n₀+m-1)を用いて次式（式(45)）より、位相φ_{n_0+1,k}、位相φ_{n_0+2,k}、…、位相φ_{n_0+m,k}と昇順に算出していく。

時系列信号特徴算出部４２０による位相の算出には、位相の初期値φ_{n_0,k}を与える必要がある。各kについて、n₀=0のときの位相φ_0,kを初期値とすることができる。具体的には、離散時間インデックスn=0の時点において、信号が無音の状態であることなどを分析の前提とすることで、適切な初期値を与えることができる。

本実施形態の発明によれば、短時間フーリエ変換の結果として得られる複素数信号から位相を求める従来の方法では位相の値が-πからπの間にラップされた状態でしか計算できないのに対し、振幅情報を用いることで、アンラップされた状態での位相の値を推定することができる。

また、位相を振幅情報のみから推定できるため、短時間フーリエ変換の結果として得られる複素数信号の情報のうち、位相情報が欠損し、振幅情報のみが与えられた状況においても、完全な複素数信号が与えられた場合と同等の位相の推定が可能となる。

さらに、振幅情報のみが与えられた条件において位相の推定値を復元することで、一部の情報が欠損した中で複素数信号を再構成することが可能となる。

＜変形例＞
この発明は上述の実施形態に限定されるものではなく、この発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。上記実施形態において説明した各種の処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。

例えば、隣接するインデックスの位相の推定値の差分を、第三実施形態における式(43)や式(44)のような周波数方向の再帰式で得られる推定値と第四実施形態における式(45)のような時間方向の再帰式で得られる推定値との平均、重み付き平均により計算してもよい。

＜補記＞
本発明の装置は、例えば単一のハードウェアエンティティとして、キーボードなどが接続可能な入力部、液晶ディスプレイなどが接続可能な出力部、ハードウェアエンティティの外部に通信可能な通信装置（例えば通信ケーブル）が接続可能な通信部、ＣＰＵ（Central Processing Unit、キャッシュメモリやレジスタなどを備えていてもよい）、メモリであるＲＡＭやＲＯＭ、ハードディスクである外部記憶装置並びにこれらの入力部、出力部、通信部、ＣＰＵ、ＲＡＭ、ＲＯＭ、外部記憶装置の間のデータのやり取りが可能なように接続するバスを有している。また必要に応じて、ハードウェアエンティティに、ＣＤ−ＲＯＭなどの記録媒体を読み書きできる装置（ドライブ）などを設けることとしてもよい。このようなハードウェア資源を備えた物理的実体としては、汎用コンピュータなどがある。

ハードウェアエンティティの外部記憶装置には、上述の機能を実現するために必要となるプログラムおよびこのプログラムの処理において必要となるデータなどが記憶されている（外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくこととしてもよい）。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶される。

ハードウェアエンティティでは、外部記憶装置（あるいはＲＯＭなど）に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてメモリに読み込まれて、適宜にＣＰＵで解釈実行・処理される。その結果、ＣＰＵが所定の機能（上記、…部、…手段などと表した各構成要件）を実現する。

本発明は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、上記実施形態において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。

既述のように、上記実施形態において説明したハードウェアエンティティ（本発明の装置）における処理機能をコンピュータによって実現する場合、ハードウェアエンティティが有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記ハードウェアエンティティにおける処理機能がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。

また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。

また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ハードウェアエンティティを構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

Claims (8)

1. 窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間t_n、離散角周波数ω_k(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
   

   

   
   

   

   
   (ただし、f_sは標本化周波数、Nは離散フーリエ変換の点数である）とし、
   前記時系列信号を離散短時間フーリエ変換することにより得られる離散時間t _n 、離散角周波数ω _k における振幅A _n,k (n=n ₁ , n ₁ +1,…, n ₂ （ただし、n ₁ ,n ₂ はn ₂ -n ₁ ≧1を満たす整数）、kはある整数)、前記標準偏差σ、前記標本化周波数f _s 、前記離散フーリエ変換の点数Nは、いずれも事前に記録部に記録されているか、処理開始時に入力として与えられるものであり、
   前記振幅A_n,k (n=n ₁ , n ₁ +1,…, n ₂ )、前記標本化周波数f _s から、離散角周波数ω_kについて離散時間t_nにおける対数振幅logA_n,kの時間差分である対数振幅時間差分Δ_tlogA_n,k/Δt_nを算出する対数振幅時間差分算出部と、
   前記対数振幅時間差分Δ_tlogA_n,k/Δt_n (n=n ₁ , n ₁ +1,…, n ₂ )、前記標準偏差σ、前記標本化周波数f _s 、前記離散フーリエ変換の点数Nから、次式を用いて離散時間t_n、離散角周波数ω_kにおける群遅延特性GD_n,kを時系列信号特徴として算出する時系列信号特徴算出部と
   

   

   
   (ただし、Tは窓関数の有効幅であり、T=N/f_s)
   を含む時系列信号特徴推定装置。
2. 請求項１に記載の時系列信号特徴推定装置であって、
   前記対数振幅時間差分算出部は、前記対数振幅時間差分Δ_tlogA_n,k/Δt_nを、n₁<n<n₂を満たすnについては中心差分近似により算出し、
   δ₁、δ₂を非零の数とし、
   前記時系列信号特徴算出部は、前記群遅延特性GD_n,kを、n₁<n<n₂を満たすnについては次式により算出する
   

   

   
   ことを特徴とする時系列信号特徴推定装置。
3. 窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間t_n、離散角周波数ω_k(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
   

   

   
   

   

   
   (ただし、f_sは標本化周波数、Nは離散フーリエ変換の点数である）とし、
   前記時系列信号を離散短時間フーリエ変換することにより得られる離散時間t _n 、離散角周波数ω _k における振幅A _n,k (nはある整数、k=k ₁ , k ₁ +1,…, k ₂ （ただし、k ₁ ,k ₂ はk ₂ -k ₁ ≧1を満たす整数）)、前記標準偏差σ、前記標本化周波数f _s 、前記離散フーリエ変換の点数Nは、いずれも事前に記録部に記録されているか、処理開始時に入力として与えられるものであり、
   前記振幅A_n,k (k=k ₁ , k ₁ +1,…, k ₂ )、前記標本化周波数f _s 、前記離散フーリエ変換の点数Nから、離散時間t_nについて離散角周波数ω_kにおける対数振幅logA_n,kの角周波数差分である対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kを算出する対数振幅角周波数差分算出部と、
   前記対数振幅角周波数差分Δ_ωlogA_n,k/Δω_k (k=k ₁ , k ₁ +1,…, k ₂ )、前記標準偏差σ、前記標本化周波数f _s 、前記離散フーリエ変換の点数Nから、次式を用いて離散時間t_n、離散角周波数ω_kにおける瞬時周波数IF_n,kを時系列信号特徴として算出する時系列信号特徴算出部と
   

   

   
   を含む時系列信号特徴推定装置。
4. 請求項３に記載の時系列信号特徴推定装置であって、
   前記対数振幅角周波数差分算出部は、前記対数振幅角周波数差分Δ_ωlogA_n,k/Δω_kを、k₁<k<k₂を満たすkについては中心差分近似により算出し、
   δ₄、δ₅を非零の数とし、
   前記時系列信号特徴算出部は、前記瞬時周波数IF_n,kを、k₁<k<k₂を満たすkについては次式により算出する
   

   

   
   ことを特徴とする時系列信号特徴推定装置。
5. 窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間t_n、離散角周波数ω_k(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
   

   

   
   

   

   
   (ただし、f_sは標本化周波数、Nは離散フーリエ変換の点数である）とし、
   前記時系列信号を離散短時間フーリエ変換することにより得られる離散時間t _n 、離散角周波数ω _k における振幅A _n,k (n=n ₁ , n ₁ +1,…, n ₂ （ただし、n ₁ ,n ₂ はn ₂ -n ₁ ≧1を満たす整数）、k=0, 1,…, N-1)、前記標準偏差σ、前記標本化周波数f _s 、前記離散フーリエ変換の点数Nは、いずれも事前に記録部に記録されているか、処理開始時に入力として与えられるものであり、
   前記振幅A_n,k (n=n ₁ , n ₁ +1,…, n ₂ 、k=0, 1,…, N-1)、前記標準偏差σ、前記標本化周波数f _s 、前記離散フーリエ変換の点数Nから、k=0,1,…, N-2の各kについて離散角周波数ω_kからω_k+1にかけての群遅延特性の符号反転値の積分値である群遅延特性積分値Σ(-GD)_n,kを算出する群遅延特性積分値算出部と、
   k₀を0以上N-1以下のある整数、位相の初期値をφ_{n,k_0} (n=n ₁ , n ₁ +1,…, n ₂ )とし、
   前記整数k ₀ 、前記位相の初期値φ _{n,k_0} (n=n ₁ , n ₁ +1,…, n ₂ )は、いずれも事前に記録部に記録されているか、処理開始時に入力として与えられるものであり、
   n=n ₁ , n ₁ +1,…, n ₂ の各nについて、前記群遅延特性積分値Σ(-GD)_n,k(k=k₀, k₀+1,…, N-2)から、次式を用いて、k=k₀+1から昇順に位相φ_n,k(k=k₀+1,…, N-1)を時系列信号特徴として算出し、
   

   

   
   前記群遅延特性積分値Σ(-GD)_n,k(k=0, 1,…,k₀-1)から、次式を用いて、k=k₀-1から降順に位相φ_n,k(k=0, 1,…,k₀-1)を時系列信号特徴として算出する時系列信号特徴算出部と
   

   

   
   を含む時系列信号特徴推定装置。
6. 請求項５に記載の時系列信号特徴推定装置であって、
   前記群遅延特性積分値算出部は、請求項１に記載の時系列信号特徴推定装置を用いて、前記振幅A_n,k(n=n₁, n₁+1,…, n₂、k=0, 1,…, N-1)から、k=0,1,…, N-1の各kについて群遅延特性GD_n,kを算出し、前記群遅延特性GD_n,kから、k=0,1,…, N-2の各kについて前記群遅延特性積分値Σ(-GD)_n,kを算出する
   ことを特徴とする時系列信号特徴推定装置。
7. 窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間t_n、離散角周波数ω_k(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
   

   

   
   

   

   
   (ただし、f_sは標本化周波数、Nは離散フーリエ変換の点数である）とし、
   前記時系列信号を離散短時間フーリエ変換することにより得られる離散時間t _n 、離散角周波数ω _k における振幅A _n,k (n=n ₀ , n ₀ +1,…, n ₀ +m（ただし、n ₀ は整数、mは正の整数）、k=k ₁ , k ₁ +1,…, k ₂ （ただし、k ₁ ,k ₂ はk ₂ -k ₁ ≧1を満たす整数）)、前記標準偏差σ、前記標本化周波数f _s 、前記離散フーリエ変換の点数Nは、いずれも事前に記録部に記録されているか、処理開始時に入力として与えられるものであり、
   前記振幅A_n,k (n=n ₀ , n ₀ +1,…, n ₀ +m、k=k ₁ , k ₁ +1,…, k ₂ )、前記標準偏差σ、前記標本化周波数f _s 、前記離散フーリエ変換の点数Nから、n=n₀, n₀+1,…, n₀+m-1の各nについて離散時間t_nからt_n+1にかけての瞬時周波数の2π倍の積分値である瞬時周波数積分値Σ(2πIF)_n,kを算出する瞬時周波数算出部と、
   位相の初期値をφ_{n_0,k} (k=k ₁ , k ₁ +1,…, k ₂ )とし、
   前記位相の初期値φ _{n_0,k} (k=k ₁ , k ₁ +1,…, k ₂ )は、事前に記録部に記録されているか、処理開始時に入力として与えられるものであり、
   k=k ₁ , k ₁ +1,…, k ₂ の各kについて、前記瞬時周波数積分値Σ(2πIF)_n,k(n=n₀, n₀+1,…, n₀+m-1)から、次式を用いて、n=n₀+1から昇順に位相φ_n,k(n=n₀+1,…, n₀+m)を時系列信号特徴として算出する時系列信号特徴算出部と
   

   

   
   を含む時系列信号特徴推定装置。
8. 請求項１ないし７のいずれか１項に記載の時系列信号特徴推定装置としてコンピュータを機能させるためのプログラム。

Images