(57)【要約】【課題】短時間フーリエ変換に基づく時系列信号の特徴分析において、不確定性が少ない形で、時系列信号特徴を算出することができる時系列信号特徴推定技術を提供する。【解決手段】離散時間tn、離散角周波数ωkにおける振幅An,kから、離散角周波数ωkについて離散時間tnにおける対数振幅logAn,kの時間差分である対数振幅時間差分ΔtlogAn,k/Δtnを算出する対数振幅時間差分算出部１１０と、対数振幅時間差分ΔtlogAn,k/Δtnから、所定の式を用いて離散時間tn、離散角周波数ωkにおける群遅延特性GDn,kを時系列信号特徴として算出する時系列信号特徴算出部１２０とを含む。【選択図】図１
2018097430.tif 000002
【特許請求の範囲】
【請求項１】
窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間tn、離散角周波数ωk(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
【数９６】
2018097430.tif 000098
【数９７】
2018097430.tif 000099
(ただし、fsは標本化周波数、Nは離散フーリエ変換の点数である）、
離散時間tn、離散角周波数ωkにおける振幅をAn,kとし、
前記振幅An,k(n1,n2はn2-n1≧1を満たす整数、n=n1, n1+1,…, n2、kは整数)から、離散角周波数ωkについて離散時間tnにおける対数振幅logAn,kの時間差分である対数振幅時間差分ΔtlogAn,k/Δtnを算出する対数振幅時間差分算出部と、
前記対数振幅時間差分ΔtlogAn,k/Δtnから、次式を用いて離散時間tn、離散角周波数ωkにおける群遅延特性GDn,kを時系列信号特徴として算出する時系列信号特徴算出部と
【数９８】
2018097430.tif 000100
(ただし、Tは窓関数の有効幅であり、T=N/fs)
を含む時系列信号特徴推定装置。
【請求項２】
請求項１に記載の時系列信号特徴推定装置であって、
前記対数振幅時間差分算出部は、前記対数振幅時間差分ΔtlogAn,k/Δtnを、n1<n<n2を満たすnについては中心差分近似により算出し、
δ1、δ2を非零の数とし、
前記時系列信号特徴算出部は、前記群遅延特性GDn,kを、n1<n<n2を満たすnについては次式により算出する
【数９９】
2018097430.tif 000101
ことを特徴とする時系列信号特徴推定装置。
【請求項３】
窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間tn、離散角周波数ωk(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
【数１００】
2018097430.tif 000102
【数１０１】
2018097430.tif 000103
(ただし、fsは標本化周波数、Nは離散フーリエ変換の点数である）、
離散時間tn、離散角周波数ωkにおける振幅をAn,kとし、
前記振幅An,k(nは整数、k1,k2はk2-k1≧1を満たす整数、k=k1, k1+1,…, k2)から、離散時間tnについて離散角周波数ωkにおける対数振幅logAn,kの角周波数差分である対数振幅角周波数差分ΔωlogAn,k/Δωkを算出する対数振幅角周波数差分算出部と、
前記対数振幅角周波数差分ΔωlogAn,k/Δωkから、次式を用いて離散時間tn、離散角周波数ωkにおける瞬時周波数IFn,kを時系列信号特徴として算出する時系列信号特徴算出部と
【数１０２】
2018097430.tif 000104
を含む時系列信号特徴推定装置。
【請求項４】
請求項３に記載の時系列信号特徴推定装置であって、
前記対数振幅角周波数差分算出部は、前記対数振幅角周波数差分ΔωlogAn,k/Δωkを、k1<k<k2を満たすkについては中心差分近似により算出し、
δ4、δ5を非零の数とし、
前記時系列信号特徴算出部は、前記瞬時周波数IFn,kを、k1<k<k2を満たすkについては次式により算出する
【数１０３】
2018097430.tif 000105
ことを特徴とする時系列信号特徴推定装置。
【請求項５】
窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間tn、離散角周波数ωk(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
【数１０４】
2018097430.tif 000106
【数１０５】
2018097430.tif 000107
(ただし、fsは標本化周波数、Nは離散フーリエ変換の点数である）、
離散時間tn、離散角周波数ωkにおける振幅をAn,kとし、
前記振幅An,k(n1,n2はn2-n1≧1を満たす整数、n=n1, n1+1,…, n2、k=0, 1,…, N-1)から、k=0,1,…, N-2の各kについて離散角周波数ωkからωk+1にかけての群遅延特性の符号反転値の積分値である群遅延特性積分値Σ(-GD)n,kを算出する群遅延特性積分値算出部と、
k0を0以上N-1以下のある整数、位相の初期値をφn,k_0とし、
前記群遅延特性積分値Σ(-GD)n,k(k=k0, k0+1,…, N-2)から、次式を用いて、k=k0+1から昇順に位相φn,k(k=k0+1,…, N-1)を時系列信号特徴として算出し、
【数１０６】
2018097430.tif 000108
前記群遅延特性積分値Σ(-GD)n,k(k=0, 1,…,k0-1)から、次式を用いて、k=k0-1から降順に位相φn,k(k=0, 1,…,k0-1)を時系列信号特徴として算出する時系列信号特徴算出部と
【数１０７】
2018097430.tif 000109
を含む時系列信号特徴推定装置。
【請求項６】
請求項５に記載の時系列信号特徴推定装置であって、
前記群遅延特性積分値算出部は、請求項１に記載の時系列信号特徴推定装置を用いて、前記振幅An,k(n=n1, n1+1,…, n2、k=0, 1,…, N-1)から、k=0,1,…, N-1の各kについて群遅延特性GDn,kを算出し、前記群遅延特性GDn,kから、k=0,1,…, N-2の各kについて前記群遅延特性積分値Σ(-GD)n,kを算出する
ことを特徴とする時系列信号特徴推定装置。
【請求項７】
窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間tn、離散角周波数ωk(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
【数１０８】
2018097430.tif 000110
【数１０９】
2018097430.tif 000111
(ただし、fsは標本化周波数、Nは離散フーリエ変換の点数である）、
離散時間tn、離散角周波数ωkにおける振幅をAn,kとし、
前記振幅An,k(n0を整数、mを正の整数、n=n0, n0+1,…, n0+m、k1,k2はk2-k1≧1を満たす整数、k=k1, k1+1,…, k2)から、n=n0, n0+1,…, n0+m-1の各nについて離散時間tnからtn+1にかけての瞬時周波数の2π倍の積分値である瞬時周波数積分値Σ(2πIF)n,kを算出する瞬時周波数算出部と、
位相の初期値をφn_0,kとし、
前記瞬時周波数積分値Σ(2πIF)n,k(n=n0, n0+1,…, n0+m-1)から、次式を用いて、n=n0+1から昇順に位相φn,k(n=n0+1,…, n0+m)を時系列信号特徴として算出する時系列信号特徴算出部と
【数１１０】
2018097430.tif 000112
を含む時系列信号特徴推定装置。
【請求項８】
請求項１ないし７のいずれか１項に記載の時系列信号特徴推定装置としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、時系列信号特徴推定技術に関し、特に短時間フーリエ変換に基づく時系列信号の特徴を推定する技術に関する。
【背景技術】
【０００２】
音声信号や音響信号などの時系列信号を分析する手法の１つとして、短時間フーリエ変換がある（非特許文献１参照）。分析対象となる時系列信号は、短時間フーリエ変換された結果、時間と角周波数（時間と周波数）のインデックスを持つ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次元の系列となってしまう。このため、群遅延特性や瞬時角周波数の値を算出するためには、位相成分の不連続を解消しアンラップされた位相成分の特性を知る必要があるが、一般に位相成分の特性のアンラップ処理は一意な解を持たず不確定な問題となる。
【０００７】
そこで本発明は、短時間フーリエ変換に基づく時系列信号の特徴分析において、不確定性が少ない形で、時系列信号特徴を算出することができる時系列信号特徴推定技術を提供することを目的とする。
【課題を解決するための手段】
【０００８】
本発明の一態様は、窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間tn、離散角周波数ωk(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
【０００９】
【数１】
2018097430.tif 000003
【００１０】
【数２】
2018097430.tif 000004
【００１１】
(ただし、fsは標本化周波数、Nは離散フーリエ変換の点数である）、離散時間tn、離散角周波数ωkにおける振幅をAn,kとし、前記振幅An,k(n1,n2はn2-n1≧1を満たす整数、n=n1, n1+1,…, n2、kは整数)から、離散角周波数ωkについて離散時間tnにおける対数振幅logAn,kの時間差分である対数振幅時間差分ΔtlogAn,k/Δtnを算出する対数振幅時間差分算出部と、前記対数振幅時間差分ΔtlogAn,k/Δtnから、次式を用いて離散時間tn、離散角周波数ωkにおける群遅延特性GDn,kを時系列信号特徴として算出する時系列信号特徴算出部と
【００１２】
【数３】
2018097430.tif 000005
【００１３】
(ただし、Tは窓関数の有効幅であり、T=N/fs)を含む。
【００１４】
本発明の一態様は、窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間tn、離散角周波数ωk(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
【００１５】
【数４】
2018097430.tif 000006
【００１６】
【数５】
2018097430.tif 000007
【００１７】
(ただし、fsは標本化周波数、Nは離散フーリエ変換の点数である）、離散時間tn、離散角周波数ωkにおける振幅をAn,kとし、前記振幅An,k(nは整数、k1,k2はk2-k1≧1を満たす整数、k=k1, k1+1,…, k2)から、離散時間tnについて離散角周波数ωkにおける対数振幅logAn,kの角周波数差分である対数振幅角周波数差分ΔωlogAn,k/Δωkを算出する対数振幅角周波数差分算出部と、前記対数振幅角周波数差分ΔωlogAn,k/Δωkから、次式を用いて離散時間tn、離散角周波数ωkにおける瞬時周波数IFn,kを時系列信号特徴として算出する時系列信号特徴算出部と
【００１８】
【数６】
2018097430.tif 000008
【００１９】
を含む。
【００２０】
本発明の一態様は、窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間tn、離散角周波数ωk(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
【００２１】
【数７】
2018097430.tif 000009
【００２２】
【数８】
2018097430.tif 000010
【００２３】
(ただし、fsは標本化周波数、Nは離散フーリエ変換の点数である）、離散時間tn、離散角周波数ωkにおける振幅をAn,kとし、前記振幅An,k(n1,n2はn2-n1≧1を満たす整数、n=n1, n1+1,…, n2、k=0, 1,…, N-1)から、k=0,1,…, N-2の各kについて離散角周波数ωkからωk+1にかけての群遅延特性の符号反転値の積分値である群遅延特性積分値Σ(-GD)n,kを算出する群遅延特性積分値算出部と、k0を0以上N-1以下のある整数、位相の初期値をφn,k_0とし、前記群遅延特性積分値Σ(-GD)n,k(k=k0, k0+1,…, N-2)から、次式を用いて、k=k0+1から昇順に位相φn,k(k=k0+1,…, N-1)を時系列信号特徴として算出し、
【００２４】
【数９】
2018097430.tif 000011
【００２５】
前記群遅延特性積分値Σ(-GD)n,k(k=0, 1,…,k0-1)から、次式を用いて、k=k0-1から降順に位相φn,k(k=0, 1,…,k0-1)を時系列信号特徴として算出する時系列信号特徴算出部と
【００２６】
【数１０】
2018097430.tif 000012
【００２７】
を含む。
【００２８】
本発明の一態様は、窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換の離散時間tn、離散角周波数ωk(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)を
【００２９】
【数１１】
2018097430.tif 000013
【００３０】
【数１２】
2018097430.tif 000014
【００３１】
(ただし、fsは標本化周波数、Nは離散フーリエ変換の点数である）、離散時間tn、離散角周波数ωkにおける振幅をAn,kとし、前記振幅An,k(n0を整数、mを正の整数、n=n0, n0+1,…, n0+m、k1,k2はk2-k1≧1を満たす整数、k=k1, k1+1,…, k2)から、n=n0, n0+1,…, n0+m-1の各nについて離散時間tnからtn+1にかけての瞬時周波数の2π倍の積分値である瞬時周波数積分値Σ(2πIF)n,kを算出する瞬時周波数算出部と、位相の初期値をφn_0,kとし、前記瞬時周波数積分値Σ(2πIF)n,k(n=n0, n0+1,…, n0+m-1)から、次式を用いて、n=n0+1から昇順に位相φn,k(n=n0+1,…, n0+m)を時系列信号特徴として算出する時系列信号特徴算出部と
【００３２】
【数１３】
2018097430.tif 000015
【００３３】
を含む。
【発明の効果】
【００３４】
本発明によれば、時系列信号の短時間フーリエ変換の振幅成分に基づき、位相成分、群遅延特性、瞬時周波数を算出することにより、これらの値を不確定性の少ない形で算出することが可能となる。
【図面の簡単な説明】
【００３５】
【図１】時系列信号特徴推定装置１００の構成の一例を示す図。
【図２】時系列信号特徴推定装置１００の動作の一例を示す図。
【図３】時系列信号特徴推定装置２００の構成の一例を示す図。
【図４】時系列信号特徴推定装置２００の動作の一例を示す図。
【図５】時系列信号特徴推定装置３００の構成の一例を示す図。
【図６】時系列信号特徴推定装置３００の動作の一例を示す図。
【図７】時系列信号特徴推定装置４００の構成の一例を示す図。
【図８】時系列信号特徴推定装置４００の動作の一例を示す図。
【図９】離散角周波数インデックスk=0における真の位相と群遅延特性の推定値を示す図。
【発明を実施するための形態】
【００３６】
以下、本発明の実施の形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【００３７】
＜表記方法＞
_（アンダースコア）は下付き添字を表す。例えば、xy_zはyzがxに対する上付き添字であり、xy_zはyzがxに対する下付き添字であることを表す。
【００３８】
＜短時間フーリエ変換に基づく時系列信号の特徴推定の原理＞
まず、短時間フーリエ変換に基づく時系列信号の特徴推定の原理となる群遅延特性、瞬時角周波数、位相が満たす方程式について説明する。これらの方程式を用いることにより、時系列信号の短時間フーリエ変換により得られる2次元の複素数信号の位相成分を、複素数信号の値そのものから直接計算する代わりに、複素数信号の振幅成分から推定することが可能となる。また、群遅延特性、瞬時角周波数も振幅成分から推定することが可能となる。
【００３９】
［群遅延特性、瞬時角周波数、位相が満たす方程式］
まず、時系列信号である連続時間信号x(τ)について、式(1)で与えられる時間tと角周波数ωにおける短時間フーリエ変換Xt,ωを考える。
【００４０】
【数１４】
2018097430.tif 000016
【００４１】
ただし、w(t)はt=0において最大値を取る窓関数、Tは窓関数の有効幅とする。なお、jは虚数単位である。
【００４２】
このとき、短時間フーリエ変換Xt,ωを時間tで偏微分すると式(2)が得られる。
【００４３】
【数１５】
2018097430.tif 000017
【００４４】
同様に、短時間フーリエ変換Xt,ωを角周波数ωで偏微分すると式(3)が得られる。
【００４５】
【数１６】
2018097430.tif 000018
【００４６】
ここで、窓関数として式(4)で与えられるガウス窓を選択する。
【００４７】
【数１７】
2018097430.tif 000019
【００４８】
ただし、ガウス窓の標準偏差σは0より大きい実数(正の実数)とする。
【００４９】
このとき、窓関数w(τ-t)を時間tで偏微分すると式(5)が得られる。
【００５０】
【数１８】
2018097430.tif 000020
【００５１】
このとき、式(4)を式(3)に、式(5)を式(2)に適用し、共通項を消去すると、式(6)が得られる。
【００５２】
【数１９】
2018097430.tif 000021
【００５３】
短時間フーリエ変換Xt,ωの振幅をAt,ω、位相をφt,ωとすると、式(7)が成り立つ。また、式(7)の対数をとると、式(8)が得られる。
【００５４】
【数２０】
2018097430.tif 000022
【００５５】
【数２１】
2018097430.tif 000023
【００５６】
ここで、logXt,ωを時間tで偏微分すると式(9)が得られる。
【００５７】
【数２２】
2018097430.tif 000024
【００５８】
したがって、式(8)、式(9)より、式(10)が成立する。
【００５９】
【数２３】
2018097430.tif 000025
【００６０】
同様に、logXt,ωを角周波数ωで偏微分すると式(11)が得られる。
【００６１】
【数２４】
2018097430.tif 000026
【００６２】
したがって、式(8)、式(11)より、式(12)が成立する。
【００６３】
【数２５】
2018097430.tif 000027
【００６４】
式(10)と式(12)を式(6)に代入すると、式(13)が得られる。
【００６５】
【数２６】
2018097430.tif 000028
【００６６】
ここで、式(13)の左辺の実数部と虚数部は独立して0となることから、式(14)と式(15)が得られる。
【００６７】
【数２７】
2018097430.tif 000029
【００６８】
【数２８】
2018097430.tif 000030
【００６９】
式(14)の左辺は群遅延特性の定義そのものであり、式(14)の右辺から群遅延特性が位相を用いることなく振幅情報である対数振幅のみから計算できることが分かる。また、式(15)の左辺は瞬時角周波数の定義そのものであり、式(15)の右辺から瞬時角周波数も位相を用いることなく対数振幅のみから計算できることが分かる。瞬時周波数は瞬時角周波数を2πで除することで求められることから、瞬時周波数も対数振幅のみから計算できることが分かる。
【００７０】
さらに、式(14)を角周波数について積分すると式(16)が得られる。
【００７１】
【数２９】
2018097430.tif 000031
【００７２】
ただし、C1は積分定数である。
【００７３】
同様に、式(15)を時間について積分すると、式(17)が得られる。
【００７４】
【数３０】
2018097430.tif 000032
【００７５】
ただし、C2は積分定数である。
【００７６】
なお、積分定数C1、C2は別途決定する必要がある。
【００７７】
式(16)、式(17)の左辺は位相φt,ωであり、位相φt,ωもまた対数振幅のみから計算できることが分かる。
【００７８】
以上より、短時間フーリエ変換の窓関数としてガウス窓を選択することで、群遅延特性は式(14)を、瞬時角周波数（瞬時周波数）は式(15)を、位相は式(16)または式(17)を用いてそれぞれ対数振幅から計算できることが分かる。
【００７９】
群遅延特性や瞬時角周波数（瞬時周波数）は、ラップされた位相成分の代わりに振幅成分を用いて計算することで、不確定性の問題の影響を受けにくくなる。また、位相そのものの計算についても、複素数信号から直接計算する場合にはラップされた状態でしか得られないが、式(16)または式(17)を用いることで、物理的特徴の理解しやすいアンラップされた状態で計算することができる。
【００８０】
以上の短時間フーリエ変換に基づく時系列信号の特徴推定の原理となる群遅延特性、瞬時角周波数、位相が満たす方程式は、連続時間信号についての短時間フーリエ変換について説明したものである。しかし、実際にコンピュータで計算する場合には、離散化されたデジタル信号についての短時間フーリエ変換の結果を利用する必要がある。
【００８１】
そこで、以下では、連続値を取る時間tと角周波数ωをそれぞれ離散化し、群遅延特性、瞬時周波数、位相を推定するための近似式について説明する。
【００８２】
［群遅延特性、瞬時周波数、位相を推定するための近似式］
時系列信号である連続時間信号x(τ)の、標準偏差がσ(σは正の実数)であるガウス窓を窓関数とする短時間フーリエ変換Xt,ωに対して、離散フーリエ変換を適用して得られる離散短時間フーリエ変換Xn,k(nは離散時間インデックスを表す整数、kは離散角周波数インデックスを表す整数)、つまり、窓関数として標準偏差σ(σは正の実数)のガウス窓を用いて分析される時系列信号の離散短時間フーリエ変換Xn,kを考える。
【００８３】
これにより、時間tを離散化した離散時間tn、角周波数ωを離散化した離散角周波数ωkは、それぞれ式(18)、式(19)で表される。
【００８４】
【数３１】
2018097430.tif 000033
【００８５】
【数３２】
2018097430.tif 000034
【００８６】
ただし、fsは標本化周波数、Nは離散フーリエ変換(DFT; Discrete Fourier Transform)の点数である。
【００８７】
また、このとき、窓関数の有効幅Tは、式(20)で表される。
【００８８】
【数３３】
2018097430.tif 000035
【００８９】
離散化されたデジタル信号x(τi)について、離散時間tn、離散角周波数ωkにおける離散短時間フーリエ変換Xn,kは、式(21)で表される。
【００９０】
【数３４】
2018097430.tif 000036
【００９１】
式(21)は、式(1)に対応するものである。
【００９２】
（群遅延特性を推定するための近似式）
式(14)を離散化することにより、群遅延特性を推定するための近似式を求める。
【００９３】
離散時間tn、離散角周波数ωkにおける振幅をAn,kとする。また、離散角周波数ωkについて離散時間tnにおける対数振幅logAn,kの時間差分である対数振幅時間差分をΔtlogAn,k/Δtnとする。
【００９４】
離散角周波数ωkについて離散時間tnにおける対数振幅時間差分ΔtlogAn,k/Δtnとして中心差分近似を用いると、式(22)を得る。
【００９５】
【数３５】
2018097430.tif 000037
【００９６】
ただし、δ1、δ2は数値安定化のために導入した非零の数とする。
【００９７】
式(14)と式(22)から、離散時間tn、離散角周波数ωkにおける群遅延特性GDn,kを推定するための近似式である式(23)を得る。
【００９８】
【数３６】
2018097430.tif 000038
【００９９】
また、式(24)の関係から、中心差分近似を用いた対数振幅時間差分ΔtlogAn,k/Δtnとして式(22)の代わりに式(25)を得る。
【０１００】
【数３７】
2018097430.tif 000039
【０１０１】
【数３８】
2018097430.tif 000040
【０１０２】
ただし、δ3は数値安定化のために導入した非零の数とする。なお、ΔtAn,k/Δtnは離散角周波数ωkについて離散時間tnにおける振幅An,kの時間差分である振幅時間差分である。
【０１０３】
式(14)と式(25)から、離散時間tn、離散角周波数ωkにおける群遅延特性GDn,kを推定するための近似式である式(26)を得る。
【０１０４】
【数３９】
2018097430.tif 000041
【０１０５】
なお、ここでは、中心差分近似を用いて対数振幅時間差分ΔtlogAn,k/Δtnを算出したが、中心差分近似の代わりに、前進差分近似や後進差分近似を用いてもよい。例えば、式(22)に対応する前進差分近似、後進差分近似はそれぞれ式(22a)、式(22b)のようになる。
【０１０６】
【数４０】
2018097430.tif 000042
【０１０７】
【数４１】
2018097430.tif 000043
【０１０８】
同様に、式(25)に対応する前進差分近似、後進差分近似はそれぞれ式(25a)、式(25b)のようになる。
【０１０９】
【数４２】
2018097430.tif 000044
【０１１０】
【数４３】
2018097430.tif 000045
【０１１１】
ここで、ガウス窓の標準偏差σについて、σのα倍(ただし、αは１以上の整数とする)が窓関数の有効幅の1/2に相当すると仮定すると、ガウス窓の標準偏差σは式(27)を満たす。
【０１１２】
【数４４】
2018097430.tif 000046
【０１１３】
例えば、α=4とすると、このαに対応した標準偏差σは、式(28)で求めることができる。
【０１１４】
【数４５】
2018097430.tif 000047
【０１１５】
上記仮定を導入すると、選択したσの妥当性を直観的に解釈しやすくなるというメリットがある。
【０１１６】
（瞬時周波数を推定するための近似式）
式(15)を離散化し2πで除することにより、瞬時周波数を推定するための近似式を求める。
【０１１７】
離散時間tnについて離散角周波数ωkにおける対数振幅logAn,kの角周波数差分である対数振幅角周波数差分をΔωlogAn,k/Δωkとする。
【０１１８】
離散時間tnについて離散角周波数ωkにおける対数振幅角周波数差分ΔωlogAn,k/Δωkとして中心差分近似を用いると、式(29)を得る。
【０１１９】
【数４６】
2018097430.tif 000048
【０１２０】
ただし、δ4、δ5は数値安定化のために導入した非零の数とする。
【０１２１】
式(15)と式(29)から、離散時間tn、離散角周波数ωkにおける瞬時周波数IFn,kを推定するための近似式である式(30)を得る。
【０１２２】
【数４７】
2018097430.tif 000049
【０１２３】
また、式(31)の関係から、中心差分近似を用いた対数振幅角周波数差分ΔωlogAn,k/Δωkとして式(29)の代わりに式(32)を得る。
【０１２４】
【数４８】
2018097430.tif 000050
【０１２５】
【数４９】
2018097430.tif 000051
【０１２６】
ただし、δ6は数値安定化のために導入した非零の数とする。なお、ΔωAn,k/Δωkは離散時間tnについて離散角周波数ωkにおける振幅An,kの角周波数差分である振幅角周波数差分である。
【０１２７】
式(15)と式(32)から、離散時間tn、離散角周波数ωkにおける瞬時周波数IFn,kを推定するための近似式である式(33)を得る。
【０１２８】
【数５０】
2018097430.tif 000052
【０１２９】
なお、ここでは、中心差分近似を用いて対数振幅角周波数差分ΔωlogAn,k/Δωkを算出したが、中心差分近似の代わりに、前進差分近似や後進差分近似を用いてもよい。例えば、式(29)に対応する前進差分近似、後進差分近似はそれぞれ式(29a)、式(29b)のようになる。
【０１３０】
【数５１】
2018097430.tif 000053
【０１３１】
【数５２】
2018097430.tif 000054
【０１３２】
同様に、式(32)に対応する前進差分近似、後進差分近似はそれぞれ式(32a)、式(32b)のようになる。
【０１３３】
【数５３】
2018097430.tif 000055
【０１３４】
【数５４】
2018097430.tif 000056
【０１３５】
（位相を推定するための近似式）
位相を推定するための近似式を求めるため、以下では、式(23)、式(26)、式(30)、式(33)について台形近似により積分計算をし、再帰式の形式にて求める。
【０１３６】
式(23)、式(26)については、周波数方向に積分近似する。周波数方向の台形近似による積分計算の式は式(34)のようになる。
【０１３７】
【数５５】
2018097430.tif 000057
【０１３８】
式(34)は、離散時間tn、離散角周波数ωkにおける位相φn,kが、離散時間tn、離散角周波数ωk-1における位相φn,k-1に、離散角周波数ωk-1からωkにかけての群遅延特性の符号反転値（式(14)参照）の積分値を加えたものとなることを示している。式(34)の右辺の第二項が台形の面積に相当する。
【０１３９】
離散時間tn、離散角周波数ωkにおける位相φn,kは、式(34)に式(23)を適用することで、式(35)の近似式で表される。
【０１４０】
【数５６】
2018097430.tif 000058
【０１４１】
ただし、δ7、δ8は数値安定化のために導入した非零の数とする。
【０１４２】
同様に、離散時間tn、離散角周波数ωkにおける位相φn,kは、式(34)に式(26)を適用することで、式(36)の近似式で表される。
【０１４３】
【数５７】
2018097430.tif 000059
【０１４４】
また、式(34)から式(37)を得る。
【０１４５】
【数５８】
2018097430.tif 000060
【０１４６】
式(37)は、離散時間tn、離散角周波数ωkにおける位相φn,kが、離散時間tn、離散角周波数ωk+1における位相φn,k+1から、離散角周波数ωkからωk+1にかけての群遅延特性の符号反転値（式(14)参照）の積分値を差し引いたものとなることを示している。式(37)の右辺の第二項が台形の面積に相当する。式(37)を用いると、周波数方向について負の方向に計算を進めていく形で位相を求めることができることが分かる。
【０１４７】
離散時間tn、離散角周波数ωkにおける位相φn,kはを求める近似式として、式(37)に対応する形で式(35)、式(36)からそれぞれ式(38)、式(39)を得る。
【０１４８】
【数５９】
2018097430.tif 000061
【０１４９】
【数６０】
2018097430.tif 000062
【０１５０】
次に、式(30)、式(33)については、時間方向に積分近似する。時間方向の台形近似による積分計算の式は式(40)のようになる。
【０１５１】
【数６１】
2018097430.tif 000063
【０１５２】
式(40)は、離散時間tn、離散角周波数ωkにおける位相φn,kが、離散時間tn-1、離散角周波数ωkにおける位相φn-1,kに、離散時間tn-1からtnにかけての瞬時周波数の2π倍（式(15)参照）の積分値を加えたものとなることを示している。式(40)の右辺の第二項が台形の面積に相当する。
【０１５３】
離散時間tn、離散角周波数ωkにおける位相φn,kは、式(40)に式(30)を適用することで、式(41)の近似式で表される。
【０１５４】
【数６２】
2018097430.tif 000064
【０１５５】
ただし、δ9、δ10は数値安定化のために導入した非零の数とする。
【０１５６】
同様に、離散時間tn、離散角周波数ωkにおける位相φn,kは、式(40)に式(33)を適用することで、式(42)の近似式で表される。
【０１５７】
【数６３】
2018097430.tif 000065
【０１５８】
式(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は、群遅延特性GDn,k-1と群遅延特性GDn,kを用いて計算される値である。
【０１６０】
【数６４】
2018097430.tif 000066
【０１６１】
なお、式(23)や式(26)（あるいはその導出に用いる式(22)、式(22a)、式(22b)や式(25)、式(25a)、式(25b)）を見ればわかるように、群遅延特性GDn,kは、一般に振幅An-1,k、振幅An,k、振幅An+1,kのうち、少なくとも2つの値を用いて計算することができる。したがって、群遅延特性積分値Σ(-GD)n,k-1は、振幅An-1,k-1、振幅An,k-1、振幅An+1,k-1のうち、少なくとも2つの値と、振幅An-1,k、振幅An,k、振幅An+1,kのうち、少なくとも2つの値を用いて計算される値でもある。
【０１６２】
同様に、離散角周波数ωkからωk+1にかけての群遅延特性の符号反転値の積分値（群遅延特性積分値）をΣ(-GD)n,kとすると、式(37)に相当する式として式(44)が得られる。ここで、群遅延特性積分値Σ(-GD)n,kは、群遅延特性GDn,kと群遅延特性GDn,k+1を用いて計算される値である。
【０１６３】
【数６５】
2018097430.tif 000067
【０１６４】
なお、群遅延特性積分値Σ(-GD)n,kは、振幅An-1,k、振幅An,k、振幅An+1,kのうち、少なくとも2つの値と、振幅An-1,k+1、振幅An,k+1、振幅An+1,k+1のうち、少なくとも2つの値を用いて計算される値でもある。
【０１６５】
また、離散時間tn-1からtnにかけての瞬時周波数の2π倍の積分値（以下、瞬時周波数の2π倍の積分値を単に瞬時周波数積分値という）をΣ(2πIF)n-1,kとすると、式(40)に相当する式として式(45)が得られる。ここで、瞬時周波数積分値Σ(2πIF)n-1,kは、瞬時周波数IFn-1,kと瞬時周波数IFn,kを用いて計算される値である。
【０１６６】
【数６６】
2018097430.tif 000068
【０１６７】
なお、式(30)や式(33)（あるいはその導出に用いる式(29)、式(29a)、式(29b)や式(32)、式(32a)、式(32b)）を見ればわかるように、瞬時周波数IFn,kは、一般に振幅An,k-1、振幅An,k、振幅An,k+1のうち、少なくとも2つの値を用いて計算することができる。したがって、瞬時周波数積分値Σ(2πIF)n-1,kは、振幅An-1,k-1、振幅An-1,k、振幅An-1,k+1のうち、少なくとも2つの値と、振幅An,k-1、振幅An,k、振幅An,k+1のうち、少なくとも2つの値を用いて計算される値でもある。
【０１６８】
＜第一実施形態＞
以下、図１〜図２を参照して時系列信号特徴推定装置１００について説明する。図１に示すように時系列信号特徴推定装置１００は、対数振幅時間差分算出部１１０、時系列信号特徴算出部１２０、記録部１９０を含む。記録部１９０は、時系列信号特徴推定装置１００の処理に必要な情報を適宜記録する構成部である。例えば、記録部１９０には、標本化周波数fs、離散フーリエ変換の点数N、ガウス窓の標準偏差σを事前に記録しておく。もちろん、標本化周波数fs、離散フーリエ変換の点数N、ガウス窓の標準偏差σを事前に記録部１９０に記録しておく代わりに、時系列信号特徴推定装置１００の入力として与えるのでもよい。
【０１６９】
時系列信号特徴推定装置１００の入力となる振幅An,k(n1,n2はn2-n1≧1を満たす整数、n=n1, n1+1,…, n2、kは整数)は、事前に得られているものとする。例えば、メモリや伝送容量を削減するため、事前に時系列信号を短時間フーリエ変換した結果のうち、振幅An,kのみを記録部１９０に記録しておいてもよい。また、時系列信号特徴推定装置１００による処理開始前に外部から振幅An,kが伝送されてくるのでもよい。
【０１７０】
時系列信号特徴推定装置１００は、振幅An,k(n=n1, n1+1,…,n2、kはある整数)から、離散時間tn、離散角周波数ωkにおける時系列信号特徴である群遅延特性GDn,kを推定し、出力する。
【０１７１】
図２に従い時系列信号特徴推定装置１００の動作について説明する。対数振幅時間差分算出部１１０は、振幅An,k(n=n1, n1+1,…,n2、kはある整数)から、離散角周波数ωkについて離散時間tnにおける対数振幅logAn,kの時間差分である対数振幅時間差分ΔtlogAn,k/Δtnを算出する（Ｓ１１０）。対数振幅時間差分ΔtlogAn,k/Δtnは、n1<n<n2を満たすnについては中心差分近似、端点であるn=n1またはn=n2については前進差分近似または後進差分近似により算出すればよい。もちろん、n1<n<n2を満たすnについて、中心差分近似の代わりに前進差分近似や後進差分近似を用いてもよい。対数振幅時間差分ΔtlogAn,k/Δtnを、n1<n<n2を満たすnについては中心差分近似、n=n1については前進差分近似、n=n2については後進差分近似により算出する場合、例えば、式(22)、式(22b)、式(22b)を用いて算出することができる。
【０１７２】
時系列信号特徴算出部１２０は、Ｓ１１０で算出した対数振幅時間差分ΔtlogAn,k/Δtnから、離散時間tn、離散角周波数ωkにおける群遅延特性GDn,kを時系列信号特徴として算出する（Ｓ１２０）。具体的には、ガウス窓の標準偏差σ、窓関数の有効幅T=N/fsを用いて
【０１７３】
【数６７】
2018097430.tif 000069
【０１７４】
と計算すればよい。例えば、式(23)を用いて、n1<n<n2を満たすnについては中心差分近似、n=n1については前進差分近似、n=n2については後進差分近似により対数振幅時間差分ΔtlogAn,k/Δtnを計算する場合は、それぞれ
【０１７５】
【数６８】
2018097430.tif 000070
【０１７６】
【数６９】
2018097430.tif 000071
【０１７７】
【数７０】
2018097430.tif 000072
【０１７８】
となる。また、式(26)を用いて、n1<n<n2を満たすnについては中心差分近似、n=n1については前進差分近似、n=n2については後進差分近似により対数振幅時間差分ΔtlogAn,k/Δtnを計算する場合は、それぞれ
【０１７９】
【数７１】
2018097430.tif 000073
【０１８０】
【数７２】
2018097430.tif 000074
【０１８１】
【数７３】
2018097430.tif 000075
【０１８２】
となる。
【０１８３】
本実施形態の発明によれば、短時間フーリエ変換による信号分析において、群遅延特性の推定値を本来の定義である位相を用いて計算する代わりに、安定的に計算可能な振幅情報（対数振幅や振幅）を用いて計算することで、位相のアンラップ処理の不確定性の問題を回避することができる。
【０１８４】
また、群遅延特性を振幅情報のみから推定できるため、短時間フーリエ変換の結果として得られる複素数信号の情報のうち、位相情報が欠損し、振幅情報のみが与えられた状況においても、完全な複素数信号が与えられた場合と同等の群遅延特性の推定が可能となる。
【０１８５】
＜第二実施形態＞
以下、図３〜図４を参照して時系列信号特徴推定装置２００について説明する。図３に示すように時系列信号特徴推定装置２００は、対数振幅角周波数差分算出部２１０、時系列信号特徴算出部２２０、記録部１９０を含む。記録部１９０は、時系列信号特徴推定装置２００の処理に必要な情報を適宜記録する構成部である。
【０１８６】
時系列信号特徴推定装置２００の入力となる振幅An,k(nは整数、k1,k2はk2-k1≧1を満たす整数、k=k1, k1+1,…, k2)は、事前に得られているものとする。
【０１８７】
時系列信号特徴推定装置２００は、振幅An,k(nはある整数、k=k1, k1+1,…, k2)から、離散時間tn、離散角周波数ωkにおける時系列信号特徴である瞬時周波数IFn,kを推定し、出力する。
【０１８８】
図４に従い時系列信号特徴推定装置２００の動作について説明する。対数振幅角周波数差分算出部２１０は、振幅An,k(nはある整数、k=k1, k1+1,…, k2)から、離散時間tnについて離散角周波数ωkにおける対数振幅logAn,kの角周波数差分である対数振幅角周波数差分ΔωlogAn,k/Δωkを算出する（Ｓ２１０）。対数振幅角周波数差分ΔωlogAn,k/Δωkは、k1<k<k2を満たすkについては中心差分近似、端点であるk=k1またはk=k2については前進差分近似または後進差分近似により算出すればよい。もちろん、k1<k<k2を満たすkについて、中心差分近似の代わりに前進差分近似や後進差分近似を用いてもよい。対数振幅角周波数差分ΔωlogAn,k/Δωkを、k1<k<k2を満たすkについては中心差分近似、k=k1については前進差分近似、k=k2については後進差分近似により算出する場合、例えば、式(29)、式(29b)、式(29b) を用いて算出することができる。
【０１８９】
時系列信号特徴算出部２２０は、Ｓ２１０で算出した対数振幅角周波数差分ΔωlogAn,k/Δωkから、離散時間tn、離散角周波数ωkにおける瞬時周波数IFn,kを時系列信号特徴として算出する（Ｓ２２０）。具体的には、ガウス窓の標準偏差σを用いて
【０１９０】
【数７４】
2018097430.tif 000076
【０１９１】
と計算すればよい。例えば、式(30)を用いて、k1<k<k2を満たすkについては中心差分近似、k=k1については前進差分近似、k=k2については後進差分近似により対数振幅角周波数差分ΔωlogAn,k/Δωkを計算する場合は、それぞれ
【０１９２】
【数７５】
2018097430.tif 000077
【０１９３】
【数７６】
2018097430.tif 000078
【０１９４】
【数７７】
2018097430.tif 000079
【０１９５】
となる。また、式(33)を用いて、k1<k<k2を満たすkについては中心差分近似、k=k1については前進差分近似、k=k2については後進差分近似により対数振幅角周波数差分ΔωlogAn,k/Δωkを計算する場合は、それぞれ
【０１９６】
【数７８】
2018097430.tif 000080
【０１９７】
【数７９】
2018097430.tif 000081
【０１９８】
【数８０】
2018097430.tif 000082
【０１９９】
となる。
【０２００】
本実施形態の発明によれば、短時間フーリエ変換による信号分析において、瞬時周波数の推定値を本来の定義である位相を用いて計算する代わりに、安定的に計算可能な振幅情報（対数振幅や振幅）を用いて計算することで、位相のアンラップ処理の不確定性の問題を回避することができる。
【０２０１】
また、瞬時周波数を振幅情報のみから推定できるため、短時間フーリエ変換の結果として得られる複素数信号の情報のうち、位相情報が欠損し、振幅情報のみが与えられた状況においても、完全な複素数信号が与えられた場合と同等の瞬時周波数の推定が可能となる。
【０２０２】
＜第三実施形態＞
以下、図５〜図６を参照して時系列信号特徴推定装置３００について説明する。図５に示すように時系列信号特徴推定装置３００は、群遅延特性積分値算出部３１０、時系列信号特徴算出部３２０、記録部１９０を含む。記録部１９０は、時系列信号特徴推定装置３００の処理に必要な情報を適宜記録する構成部である。
【０２０３】
時系列信号特徴推定装置３００の入力となる振幅An,k(n1,n2はn2-n1≧1を満たす整数、n=n1, n1+1,…, n2、k=0, 1,…, N-1)は、事前に得られているものとする。
【０２０４】
時系列信号特徴推定装置３００は、振幅An,k(n=n1, n1+1,…, n2、k=0, 1,…, N-1)から、離散時間tn、離散角周波数ωkにおける時系列信号特徴である位相φn,kを推定し、出力する。
【０２０５】
図６に従い時系列信号特徴推定装置３００の動作について説明する。群遅延特性積分値算出部３１０は、振幅An,k(n=n1, n1+1,…,n2、k=0, 1,…, N-1)から、k=0,1,…, N-2の各kについて離散角周波数ωkからωk+1にかけての群遅延特性の符号反転値の積分値である群遅延特性積分値Σ(-GD)n,kを算出する（Ｓ３１０）。例えば、群遅延特性積分値算出部３１０は、時系列信号特徴推定装置１００を用いて、振幅An,k(n=n1, n1+1,…,n2、k=0, 1,…, N-1)から、k=0,1,…, N-1の各kについて群遅延特性GDn,kを算出したうえで、群遅延特性積分値Σ(-GD)n,kを算出するのでよい。
【０２０６】
また、別の例として、群遅延特性積分値算出部３１０は、n1<n<n2を満たすnについては、以下の式を用いて算出し、
【０２０７】
【数８１】
2018097430.tif 000083
【０２０８】
n=n1については、以下の式を用いて算出し、
【０２０９】
【数８２】
2018097430.tif 000084
【０２１０】
n=n2については、以下の式を用いて算出してもよい。
【０２１１】
【数８３】
2018097430.tif 000085
【０２１２】
また、別の例として、群遅延特性積分値算出部３１０は、n1<n<n2を満たすnについては、以下の式を用いて算出し、
【０２１３】
【数８４】
2018097430.tif 000086
【０２１４】
n=n1については、以下の式を用いて算出し、
【０２１５】
【数８５】
2018097430.tif 000087
【０２１６】
n=n2については、以下の式を用いて算出してもよい。
【０２１７】
【数８６】
2018097430.tif 000088
【０２１８】
時系列信号特徴算出部３２０は、位相の初期値をφn,k_0(k0は0以上N-1以下のある整数) とＳ３１０で算出した群遅延特性積分値Σ(-GD)n,k(k=0, 1,…, N-2)から、位相φn,k(k=0, 1,…, N-1、ただし、k=k0を除く)を時系列信号特徴として算出する（Ｓ３２０）。整数k0と位相の初期値φn,k_0は記録部１９０に事前に記録しておいてもよいし、時系列信号特徴推定装置３００の入力として与えるのでもよい。具体的には、以下のような手順で位相φn,k(k=0, 1,…, N-1、ただし、k=k0を除く)を算出する。
【０２１９】
まず、k=k0+1とし、位相の初期値φn,k_0と群遅延特性積分値Σ(-GD)n,k(k=k0, k0+1,…, N-2)を用いて次式（式(43)）より、位相φn,k_0+1、位相φn,k_0+2、…、位相φn,N-1と昇順に算出していく。
【０２２０】
【数８７】
2018097430.tif 000089
【０２２１】
次に、k=k0-1とし、位相の初期値φn,k_0と群遅延特性積分値Σ(-GD)n,k(k=0, 1,…, k0-1)を用いて次式（式(44)）より、位相φn,k_0-1、位相φn,k_0-2、…、位相φn,0と降順に算出していく。
【０２２２】
【数８８】
2018097430.tif 000090
【０２２３】
時系列信号特徴算出部３２０による位相の算出には、位相の初期値φn,k_0を与える必要がある。各nについて、k0=0のときの位相φn,0を初期値とすることができる。分析対象の信号が実数である場合、位相φn,0の値は0かπのいずれかに限定される。離散時間インデックスnの変化に従い、位相φn,0の値が0からπまたはπから0に変化するポイントにおいて、式(23)または式(26)より推定された離散角周波数インデックスk=0あるいはその周辺の離散角周波数インデックスにおける群遅延特性の推定値には図９(b)に示すような、大きなスパイク性の変化が観測される。このような変化を検出することで、位相φn,0の値を0かπかで切り替え、適切な初期値を与えることができる。
【０２２４】
なお、時系列信号特徴算出部３２０で用いる再帰式の開始点は、k0=0に限定されるものではないから、任意のk0に対して初期値として適切な値を獲得できる場合には、それを初期値としてもよい。この場合、先ほど説明したように、kの値について昇順に再帰式を計算するだけでなく、降順に再帰式を計算することにより位相の推定値を得ることができる。
【０２２５】
本実施形態の発明によれば、短時間フーリエ変換の結果として得られる複素数信号から位相を求める従来の方法では位相の値が-πからπの間にラップされた状態でしか計算できないのに対し、振幅情報を用いることで、アンラップされた状態での位相の値を推定することができる。
【０２２６】
また、位相を振幅情報のみから推定できるため、短時間フーリエ変換の結果として得られる複素数信号の情報のうち、位相情報が欠損し、振幅情報のみが与えられた状況においても、完全な複素数信号が与えられた場合と同等の位相の推定が可能となる。
【０２２７】
さらに、振幅情報のみが与えられた条件において位相の推定値を復元することで、一部の情報が欠損した中で複素数信号を再構成することが可能となる。
【０２２８】
＜第四実施形態＞
以下、図７〜図８を参照して時系列信号特徴推定装置４００について説明する。図７に示すように時系列信号特徴推定装置４００は、瞬時周波数積分値算出部４１０、時系列信号特徴算出部４２０、記録部１９０を含む。記録部１９０は、時系列信号特徴推定装置４００の処理に必要な情報を適宜記録する構成部である。
【０２２９】
時系列信号特徴推定装置４００の入力となる振幅An,k(n0を整数、mを正の整数、n=n0, n0+1,…, n0+m、k1,k2はk2-k1≧1を満たす整数、k=k1, k1+1,…, k2)は、事前に得られているものとする。
【０２３０】
時系列信号特徴推定装置４００は、振幅An,k(n=n0, n0+1,…, n0+m、k=k1, k1+1,…, k2)から、離散時間tn、離散角周波数ωkにおける時系列信号特徴である位相φn,kを推定し、出力する。
【０２３１】
図８に従い時系列信号特徴推定装置４００の動作について説明する。瞬時周波数積分値算出部４１０は、振幅An,k(n=n0, n0+1,…, n0+m、k=k1, k1+1,…, k2)から、n=n0, n0+1,…, n0+m-1の各nについて離散時間tnからtn+1にかけての瞬時周波数の2π倍の積分値である瞬時周波数積分値Σ(2πIF)n,kを算出する（Ｓ４１０）。例えば、瞬時周波数積分値算出部４１０は、時系列信号特徴推定装置２００を用いて、振幅An,k(n=n0, n0+1,…, n0+m、k=k1, k1+1,…, k2)から、n=n0, n0+1,…, n0+mの各nについて瞬時周波数IFn,kを算出したうえで、周波数積分値Σ(2πIF)n,kを算出するのでよい。
【０２３２】
また、別の例として、瞬時周波数積分値算出部４１０は、k1<k<k2を満たすkについては、以下の式を用いて算出し、
【０２３３】
【数８９】
2018097430.tif 000091
【０２３４】
k=k1については、以下の式を用いて算出し、
【０２３５】
【数９０】
2018097430.tif 000092
【０２３６】
k=k2については、以下の式を用いて算出してもよい。
【０２３７】
【数９１】
2018097430.tif 000093
【０２３８】
また、別の例として、瞬時周波数積分値算出部４１０は、k1<k<k2を満たすkについては、以下の式を用いて算出し、
【０２３９】
【数９２】
2018097430.tif 000094
【０２４０】
k=k1については、以下の式を用いて算出し、
【０２４１】
【数９３】
2018097430.tif 000095
【０２４２】
k=k2については、以下の式を用いて算出してもよい。
【０２４３】
【数９４】
2018097430.tif 000096
【０２４４】
時系列信号特徴算出部４２０は、位相の初期値φn_0,kとＳ４１０で算出した瞬時周波数積分値Σ(2πIF)n,k(n=n0, n0+1,…, n0+m-1)から、位相φn,k(n=n0+1,…, n0+m)を時系列信号特徴として算出する（Ｓ４２０）。整数n0と位相の初期値φn_0,kは記録部１９０に事前に記録しておいてもよいし、時系列信号特徴推定装置４００の入力として与えるのでもよい。具体的には、以下のような手順で位相φn,k(n=n0+1,…, n0+m)を算出する。
【０２４５】
n=n0+1とし、位相の初期値φn_0,kと瞬時周波数積分値Σ(2πIF)n,k(n=n0, n0+1,…, n0+m-1)を用いて次式（式(45)）より、位相φn_0+1,k、位相φn_0+2,k、…、位相φn_0+m,kと昇順に算出していく。
【０２４６】
【数９５】
2018097430.tif 000097
【０２４７】
時系列信号特徴算出部４２０による位相の算出には、位相の初期値φn_0,kを与える必要がある。各kについて、n0=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）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０２６１】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ハードウェアエンティティを構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【図１】
2018097430.tif 000113
【図２】
2018097430.tif 000114
【図３】
2018097430.tif 000115
【図４】
2018097430.tif 000116
【図５】
2018097430.tif 000117
【図６】
2018097430.tif 000118
【図７】
2018097430.tif 000119
【図８】
2018097430.tif 000120
【図９】
2018097430.tif 000121