特異スペクトル解析/部分空間同定法による周波数解析

が、しかし、周波数解析には、GHAよりよい方法が存在する。時系列データから、Hankel行列を作り、特徴抽出するというのは、色々な分野で使われている手法らしいけれども、出所は不明。この方法によって、周波数解析を行うことができる。特に、名前とかはないっぽい。詳細は
Total Least Square based algorithm for time-domain NMR data fitting
http://www.sciencedirect.com/science/article/pii/S1064185884712095

Perceptual audio modeling with exponentially damped sinusoids
http://www.sciencedirect.com/science/article/pii/S0165168404002361


#部分空間同定法というものが1980年後半あたりから存在するらしく、発想としては同じものである気がする。元を正せば、Ho-Kalmanアルゴリズムとかと同じなのかもしれないけど、謎。また、別の分野で、特異スペクトル解析(SSA)という名前が付いている方法も、同じ発想に基づいている気がする。部分空間同定法と特異スペクトル解析は同じアイデアが違う分野で独立に発展したものなのかは、分からない

(参考)システム同定理論の動向
https://www.jstage.jst.go.jp/article/ieejeiss1987/115/7/115_7_860/_article/-char/ja/

Singular spectrum analysis
http://en.wikipedia.org/wiki/Singular_spectrum_analysis



実装は

import numpy as np

#-- exponential sinusoidal model
def esm(dat,Fs,Cmax=5,method=1):
   N = len(dat)
   L =N/2
   M = N+1-L
   Kmax = Cmax*2
   H = np.zeros(L*M).reshape((L,M))
   for n in xrange(L):
      for m in xrange(M):
          H[n,m] = dat[n+m]
   U,S,V = np.linalg.svd(H)
   for k,(s1,s2) in enumerate(zip(S,S[1:])):
       if s2<s1*1.0e-5 or k==len(S)-2 or k>=Kmax-1:
          K=k+1
          U_K,V_K = U[: , :K],V[: , :K]
          break
   if method==1:
      U_K_lastdelete = np.zeros( (L-1)*K ).reshape(L-1,K)
      U_K_topdelete = np.zeros( (L-1)*K ).reshape(L-1,K)
      for k1 in xrange(L-1):
         for k2 in xrange(K):
            U_K_lastdelete[k1,k2] = U_K[k1,k2]
            U_K_topdelete[k1,k2] = U_K[k1+1,k2]
      U_K_lastdelete = np.matrix(U_K_lastdelete)
      U_K_topdelete = np.matrix(U_K_topdelete)
      Q = np.linalg.inv(U_K_lastdelete.T*U_K_lastdelete)*U_K_lastdelete.T*U_K_topdelete
   else:
      U2 = np.zeros( (L-1)*2*K ).reshape( (L-1,2*K) )
      for k1 in xrange(L-1):
         for k2 in xrange(K):
            U2[k1,k2] = U_K[k1,k2]
            U2[k1,k2+K] = U_K[k1+1,k2]
      _,_,V2 = np.linalg.svd(U2)
      V12 = np.matrix(V2[K:,:K])
      V22 = np.matrix(V2[K:,K:])
      Q = -V12*np.linalg.inv(V22)
   evals , _ = np.linalg.eig(Q)
   A = np.zeros(N*K,dtype=np.complex).reshape((N,K))
   for n in xrange(N):
      for k,v in enumerate(evals):
          A[n,k] = v**n
   A = np.matrix(A)
   C = np.matrix(np.linalg.pinv(A))*np.matrix(dat).T
   C = np.array([C[k,0] for k in xrange(K)])
   amps = np.abs(C)
   phases = [np.arctan2(x.imag,x.real) for x in C]
   freqs = [np.arctan2(x.imag,x.real)*Fs/(2*np.pi) for x in evals]
   damps = np.log(np.abs(evals))
   ret = []
   for k in xrange(K):
       f = freqs[k]
       if f<0.0:continue
       ret.append( (freqs[k] , 2*amps[k] , phases[k] , damps[k]) )
   ret.sort(key=lambda x:x[0])
   return ret


例として

Fs = 400.0
N=10
print esm(np.sin(2*np.pi*22.6*np.arange(N)/Fs)+2.2*np.exp(-0.2*np.arange(N))*np.sin(0.3*np.pi+2*np.pi*28.4*np.arange(N)/Fs),Fs)

は、正弦波と指数関数的に減衰する正弦波の重ね合わせであるけども、10点で、位相・振幅・周波数・減衰係数を極めて正確に計算できる(必要なデータ点の数は、FFTの時と違って、サンプリング周波数とは関係ない。16kHzだろうが、44.1kHzだろうが、10点で計算できる)。上に書いたGHAの実装では、これはできないし、GHAは本質的には非線形な最適化問題を解く必要があるので原理的に難しいものがあるのに対して、こっちの方法は、完全に線形代数だけで済んでいる。


また、LPC分析/ARモデルのような方法では、モデル次数は任意で赤池情報量基準などを持ちだしてきて、別途モデルの次数を決めるという手続きがあるけど、この方法では、単純な正弦波の重ね合わせ程度であれば、Hankel行列の特異値の大きさから、適切な"次数"を見積もることができる(まぁ、現実の信号では、なかなか正弦波の時のように綺麗にはいかないが)



#Hankel行列を考えると、うまくいく理由は(わたしの理解では)以下のようなことだと思う。定数係数の線形離散方程式で記述される系、例として、
a x_{n+2} + b x_{n+1} + c x_{n} = 0
のような場合、ベクトル$(x_{n+2},x_{n+1},x_n)$はベクトル$(a,b,c)$と直交していると読める。なので、長さ3の時系列の断片$(x_{n+2},x_{n+1},x_n)$は、(a,b,c)と直交する二次元部分空間内に、いつでも乗っていることになる。現実には、誤差があるので、人工的なデータでない限り、完全に、二次元部分空間内に乗らないかもしれない。けども、そういう部分空間を見つけるのは、主成分分析などでよくやられることである。結局、部分空間同定法/特異スペクトル解析というのは、時系列断片を特徴ベクトルに用いた主成分分析と言えると思う。長さ4の断片を切り出してきても、$(a,b,c,0)$及び$(0,a,b,c)$と直交するから、依然として、時系列断片は、二次元空間を張る。というわけで、ある程度、自然にモデルの次数が決まる仕組みは、ここにある(主成分分析を使う場合、現実の問題に対しては、部分空間の次元が、それほど容易に決まらないのが普通で、それは、部分空間同定法/特異スペクトル解析でも、現実の信号に対しては、同様の問題がある)。後述するように、部分空間の次元の選び方は、"分解能"にも関係してくる。結局、減衰する正弦波(周波数0の場合は指数関数も含む)の一次結合を等間隔でサンプリングしたデータ点は、定数係数線形離散方程式を満たし、逆に、定数係数線形離散方程式の解は、減衰する正弦波の一次結合で書けるという、単純な事実が、この方法の基本にあり、線形代数で済む理由となってる(調波構造を持っていなくても成立するという点が、FFTに勝る重要な点と思う)(一方で、現実の信号のほとんどは、定数係数線形離散方程式なんて単純なもので記述できるとは、到底思えないから、そのへんが限界でもある)



非定常的な場合は、どうか

Fs = 400.0
dat = np.concatenate( (np.zeros(10) , np.cos(2*np.pi*120.0*np.arange(160)/Fs)) )
print esm(dat,Fs)

の結果は

[(57.166653803298772, 2.9669267589952173, -1.5265697631493844, -0.35942213305489079), 
(93.565331327607836, 13.474944195396525, -0.079124818422766544, -0.49751142234087681), 
(119.99992098672045, 0.99961049374219579, -0.00018568477795987266, 9.5007879491267817e-08),
 (128.56540649013473, 19.305963491468344, 2.0126327519301537, -0.54374022011442635), 
(168.14658392535586, 6.7648582012261427, -2.7731003650922683, -0.46578445738769148)]

であって、57Hzや168Hzのようなとこにも、強い振幅が出ている。このへんの振る舞いは、GHAより悪い。ただし、減衰係数を見ると、120Hzのみが安定していて、比べると、他は減衰が大きい。そして、120Hzにおける振幅・位相は正確である(esmで、Qの固有値を絶対値1に正規化すると、減衰係数を0に設定することになり、こうすると、57Hzや168Hzの振幅は小さくなる。代わりに、120Hzの振幅の正確さが失われる)


ノイズを含む場合も事情は似ている

Fs=400.0
dat = np.cos(2*np.pi*120.0*np.arange(160)/Fs) + np.array([0.1*np.random.random() for _ in xrange(160)])
print esm(dat,Fs)

120Hzの位相・振幅・減衰係数は

(119.99921873872711, 1.0091465576607852, 0.0030668774699237926, -0.00017205681001066048)

で、位相は、多少ずれているけども、許容できる範囲と思われる

Fs=400.0
dat = np.concatenate( (2.1*np.cos(2*np.pi*81.0*np.arange(130)/Fs) , np.cos(2*np.pi*120.0*np.arange(160)/Fs)) )
print esm(dat,Fs)

とかだと、位相・振幅は、もはや信用できない



最後に、振幅が極めて近い場合

Fs=16000.0
dat1 = np.cos(2*np.pi*120.0*np.arange(200)/Fs) +  np.cos(2*np.pi*121.0*np.arange(200)/Fs + 0.135*np.pi)
dat2 = np.cos(2*np.pi*120.0*np.arange(200)/Fs) +  np.cos(2*np.pi*120.5*np.arange(200)/Fs + 0.135*np.pi)
dat3 = np.cos(2*np.pi*121.0*np.arange(800)/Fs) +  np.cos(2*np.pi*120.5*np.arange(800)/Fs + 0.135*np.pi)
print esm(dat1,Fs)
print esm(dat2,Fs)
print esm(dat3,Fs)

dat1とdat3は正しく計算できるけれども、dat2では、正しい答えを返してこない(サンプリング周波数は16kHzなので、現実的な問題に近い)。この方法であっても、十分な分解能を得るには、FFTと同様、十分なサンプルデータが必要であるということは変わらない。とはいえ、FFTであれば、サンプリング周波数16000Hzで800点でも、16000/800=20Hzの分解能なので、それに比べると単純比較で、数十倍分解能がよいと言える(振幅比を偏らせると、もっと分解能は悪くなるが)。この精度は、特異値を、どこで切るかに依存する。実装の"s2<s1*1.0e-5"の部分を"s2<s1*1.0e-7"とかに変更すると、分解能はあがる(まぁしかし、一方で、ノイズまで取り込んでしまう可能性も増えるので、そのへんで、トレードオフがある)し、逆に、"s2<s1*1.0e-3"とかだと、分解能は下がる


この方法は、全然有名でないっぽいのだけど、試した限りでは、一番良い方法である。まとめとしては、周波数解析をするなら、FFTでやるのはやめて、SSAでやるべきじゃね的な話(SSAステマブログ