本当に一般調和解析(GHA)は離散フーリエ変換(DFT)より精度がよいのか?
タイトルは
一般化調和解析の周波数領域での振舞い : 本当にGHAはFFTより精度がよいのか?
http://ci.nii.ac.jp/naid/10016612781
から。中身は、オープンアクセスでないので、読むことができなかったから、自分で考えて、試してみた
一般調和解析(GHA)は、Wienerに始まり〜というのは、よく書いてあるけど、Wienerは、あんまり関係なくて、調波構造を持たない(基本周波数の整数倍の周波数の波の重ね合わせで表現できない)信号$x(t)$を
という形で近似しようという考え。離散フーリエ変換の場合と違って、各成分の周波数は、基本周波数の整数倍には固定されてない。近似のよさとしては、残差平方和
が最小になることを目指す($t_i$はサンプリングされた時間)
とはいえ、各周波数には、何の関係もないので、一発で、これを解くのは難しい。代わりに
が、K=1,...,Nに関して、それぞれ最小となるように順番に決めていくという方法を採用する。K=1の時
を、最小になるような$f,A_1,B_1$を決定するのが、目標となる。明らかに$f$が固定されていれば、$A_1,B_1$を決定するのは容易い。$f$の決め方は、とりあえず、適当な周波数を一定間隔ごとにテストしていくという方法が考えられる。こうしても、$A_1,B_1$の計算は、かなり高速にできるので、実用的に使えなくはない
Fast and Accurate Generalized Harmonic Analysis Using Newton's Method
http://eprints.lib.hokudai.ac.jp/dspace/handle/2115/39825
では、FFTのパワースペクトルで、大雑把な位置を決めた後、より精密化していくという手順を踏んでいる(また、"標準的な手順"で、$A\n,B\n,f_n$を決めた後、これを初期値としてNewton法で精密化を行っている)
一般調和解析の強みは、調波構造を持たない波にも適用できること、周波数分解能が(少なくとも見かけ上は)サンプル点数によって制限されないことがある。ところで、世の中には、調波構造を持つと分かっている場合もある。それで、基本周波数$f_0$が何らかの方法で見つかったと仮定する。その場合には、冒頭の近似は
という形になる。この場合は、最小二乗の意味で最適な解を見つけることは一発でできる。これは、サンプル点数x(2N)行列の疑似逆行列を計算する問題になる。特に、2N=サンプル点数という条件が成立する時には、逆行列の計算なので(可逆であれば)上の近似は厳密になる。離散フーリエ変換は、特に、$f_0$が、サンプリング周波数/サンプル点数である時に、相当する。この場合には、GHAとDFTは、同じ答えを与えるはずで、これを元にして、GHAとDFTを比較することは、一応妥当と言える(単純な比較では、明らかにGHAの方が計算量は大分不利である)
#一般調和解析は、非定常な信号(ここでは、周波数・振幅・位相が時間的に変動する波を指すこととする)にも適用できると言われることもあるが、元の定義を見れば、そのような場合には、さほど、よい結果を与えるという保証はないと思われる。
以下のコードにあるsinfitは、基本周波数を指定して、最小二乗推定を行っている。
test1で行っているFFTの計算と、1つ目のsinfitの計算は等価であるはずで、120Hzにピーク以外の点では、実質0と思える小さな値をとる。2つ目のsinfitの計算でも結果は同様であるけど、こっちは20Hzおきに200Hzまで見ているので、計算量は減っている
test2では、正しい周波数は16Hzであるけども、FFTでは、400/80=5Hzおきの値しか出ないので、15Hzにピークが出るものの、周辺の周波数ででも、無視できない値をとる。一方のsinfitでは、当然ではあるが、16Hzに明瞭なピークをとる。2つ目のsinfitでは、15Hzに明瞭なピークが出るものの、FFTと同様、10Hzや20Hzにも、見せかけの振幅が現れる(更に、位相の計算は、大きく間違っている)
test3は、正弦波ではなく、2つの正弦波の重ね合わせで、周波数は、整数比をとらない。結果は予想される通りで、3つ目のsinfitがよい結果を出す
test4は、サンプル点数がものすごく少ない場合で、FFTは解像度が40Hzおきなので、16Hzを特定することはできない。1つ目、2つ目のsinfitは正確に計算ができる。3つ目は、一応15Hzにピークが出る
# -*- coding:utf-8 -*-
import numpy as np
#-- これはGHAではない
#-- (周波数、振幅、位相)の計算を行う
def sinfit(dat , Fs , f0 , Nmax=10):
n,m = Nmax*2 , len(dat)
M = np.zeros(n*m).reshape(m,n)
for i in xrange(0,Nmax):
for t in xrange(m):
M[t,i] = np.cos(2*np.pi*i*f0*t/Fs)
M[t,i+Nmax] = np.sin(2*np.pi*i*f0*t/Fs)
A = np.dot( np.linalg.pinv(M) , np.asarray(dat) )
amps = []
for i in xrange(0,Nmax):
amps.append( (f0*i , np.sqrt(A[i]**2 + A[i+Nmax]**2) , np.arctan2(A[i+Nmax] , A[i])) )
return amps
def test1():
Fs=400.0 #-- サンプリング周波数的なもの
dat = np.cos(2*np.pi*120.0*np.arange(80)/Fs)
print np.abs(np.fft.fft(dat))[:40]
print sinfit(dat , Fs , 5.0 , Nmax=40)
print sinfit(dat , Fs , 20.0 , Nmax=10)
def test2():
Fs=400.0 #-- サンプリング周波数的なもの
dat = np.cos(2*np.pi*16.0*np.arange(80)/Fs)
print np.abs(np.fft.fft(dat))[:40]
print sinfit(dat , Fs , 16.0 , Nmax=10)
print sinfit(dat , Fs , 5.0 , Nmax=10)
print sinfit(dat , Fs , 4.0 , Nmax=10)
def test3():
Fs=400.0 #-- サンプリング周波数的なもの
dat = np.cos(2*np.pi*16.0*np.arange(80)/Fs) + 2.5*np.cos(2*np.pi*44.0*np.arange(80)/Fs)
print np.abs(np.fft.fft(dat))[:40]
print sinfit(dat , Fs , 8.0 , Nmax=10)
print sinfit(dat , Fs , 5.0 , Nmax=20)
print sinfit(dat , Fs , 4.0 , Nmax=20)
def test4():
Fs=400.0 #-- サンプリング周波数的なもの
dat = np.cos(2*np.pi*16.0*np.arange(10)/Fs)
print np.abs(np.fft.fft(dat))[:5]
print sinfit(dat , Fs , 16.0 , Nmax=5)
print sinfit(dat , Fs , 8.0 , Nmax=5)
print sinfit(dat , Fs , 5.0 , Nmax=5)
以上から、基本周波数が分かれば、調波構造を持つ波に対して、FFTより精度よく、周波数・振幅・位相を計算できる可能性がある。とはいえ、通常は、基本周波数そのものがわからないことが多い。それも、一般調和解析は、克服できる可能性がある。実際の一般調和解析のナイーブな実装は、例えば以下のようなものになる
import numpy as np
def sinfit(dat , Fs , freqs):
n,m = len(freqs)*2 , len(dat)
M = np.zeros(n*m).reshape(m,n)
for i,f in enumerate(freqs):
for t in xrange(m):
M[t,i] = np.cos(2*np.pi*f*t/Fs)
M[t,i+len(freqs)] = np.sin(2*np.pi*f*t/Fs)
A = np.dot( np.linalg.pinv(M) , np.asarray(dat) )
amps = []
for i,f in enumerate(freqs):
amps.append( (f , np.sqrt(A[i]**2 + A[i+len(freqs)]**2) , np.arctan2(A[i+len(freqs)] , A[i])) )
return amps
def gha(dat, Fs , fmin=0 , fmax=1000 , fstep=1 , maxiter=10 , eps=1.0e-3):
N = len(dat)
if len(dat)==0:return []
freqs = range(fmin,fmax,fstep)
rests = np.copy( dat )
maxP = sum(rests*rests)/N
timeline = np.arange( len(rests) )/Fs
peaks = set([])
for _ in range(maxiter):
Dmin = None
for n,f in enumerate(freqs):
if f<=0:continue
cosvals = np.cos(2 * np.pi * f * timeline)
sinvals = np.sin(2 * np.pi * f * timeline)
XX = sum(cosvals**2)
XY = sum(cosvals*sinvals)
YY = sum(sinvals**2)
XZ = sum(cosvals*rests)
YZ = sum(sinvals*rests)
D = XX*YY-XY*XY
if abs(D)<1.0e-10:continue
Cf = (YY*XZ - XY*YZ)/D
Sf = (-XY*XZ + XX*YZ)/D
R = sum( (rests - Cf*cosvals - Sf*sinvals)**2 )
if Dmin==None or R < Dmin:
Dmin = R
bestf = f
besta = Cf
bestb = Sf
assert(Dmin!=None),"GHA failed"
#-- binary search(閾値は適当)
cstep = fstep*0.5
cfreq = bestf
for _ in xrange(30): #-- 30回繰り返すと、精度は、fstep*10^(-9)Hz
cosvals = np.cos(2 * np.pi * cfreq * timeline)
sinvals = np.sin(2 * np.pi * cfreq * timeline)
C = sum(dat * cosvals)
S = sum(dat * sinvals)
dC = -sum(np.arange(N) * dat * sinvals)
dS = sum(np.arange(N) * dat * cosvals)
V = C*dC + S*dS
if abs((cfreq-bestf) - fstep)<1.0e-2:
#-- ここに来るのは、通常binary searchの失敗
break
elif abs(V/np.sqrt(C*C+S*S))<1.0:
XX = sum(cosvals**2)
XY = sum(cosvals*sinvals)
YY = sum(sinvals**2)
XZ = sum(cosvals*rests)
YZ = sum(sinvals*rests)
D = XX*YY-XY*XY
besta = (YY*XZ - XY*YZ)/D
bestb = (-XY*XZ + XX*YZ)/D
bestf = cfreq
break
elif V<0.0:
cfreq = cfreq - cstep
else:
cfreq = cfreq + cstep
cstep = cstep*0.5
rests = rests - besta*np.cos(2*math.pi*bestf*timeline) - bestb*np.sin(2*math.pi*bestf*timeline)
peaks.add( bestf )
if sum(rests*rests)<N*maxP*eps:break
if len(peaks)==0:peaks = [bestf]
ret = sinfit(dat , Fs , list(peaks))
ret.sort(key=lambda x:x[1],reverse=True)
return ret
if __name__=="__main__":
Fs = 400.0
dat = np.cos(2*np.pi*120.0*np.arange(160)/Fs)
print gha(dat,Fs,fmax=int(Fs/2))
dat = np.cos(2*np.pi*16.2*np.arange(160)/Fs) + 2.24*np.sin(2*np.pi*44.4*np.arange(160)/Fs)
print gha(dat,Fs,fmax=int(Fs/2))
こういう方法でも、例えば、有声音の基本周波数を計算するのには、役立つ(人間の声帯振動の周波数は、100~250Hz程度と言われるので、そのへんで振幅の大きい周波数を見つければいい)。FFTでは、サンプル点数で解像度が決まってしまうけれども、GHAだと、ある程度は、限界を超えられる。楽器音の場合は、周波数は等間隔より、等比間隔で見る方がよいかもしれない
最後に、信号が非定常的だと、何が起きるのか見る
Fs = 400.0 dat = np.concatenate( (np.zeros(160) , np.cos(2*np.pi*120.0*np.arange(160)/Fs)) ) print gha(dat,Fs,fmax=int(Fs/2))
これの出力は
[(119.9951171875, 0.50001094971516724, -0.012898725654374629), (119, 0.28709673969837646, 2.2072307129126307), (121, 0.28627734160698998, -2.2224076743624175), (116.4228515625, 0.086264024722059796, 2.0584357299209537), (123.57421875, 0.084800532509307178, -2.071810440145756), (118, 0.062597699071082008, -0.27172429125308256), (122, 0.061484157892231654, 0.29403956332817544), (113.8505859375, 0.04668031043667914, 1.9160936081986761), (126.146484375, 0.045291270390806194, -1.931859981372358)]
で、周波数120Hzは抽出できてるものの、何度も周辺の値が出現している。最初の120Hzが求まった後、残差は、周波数120Hzの波を2つ結合したようなもの(位相は、ずれている)になっているはず。こういうことが繰り返し起こって、(定常的でないので)残差は0にはならない。また、振幅の値も、1.0とはなっていない。人間の声も、たとえ母音であっても、完全に定常的ではないので、これに似たことが起こるように見える。
以下は、周波数が、より緩やかに変動していく場合
Fs = 400.0 dat = np.cos(2*np.pi*120.0*np.arange(160)/Fs + 0.0*2*np.pi*10.0*np.cos(np.arange(160)/Fs)) print gha(dat,Fs,fmax=int(Fs/2))
この場合、瞬時周波数は、110~130Hzにあると考えられ、GHAも、この範囲の複数の周波数を返してくる。振幅や位相は、まちまちであって、振幅や位相について、意味のある情報を引き出すのは難しいように思う
まとめると、
・GHAは、定常信号に対してはDFTより精度よく周波数解析を行うことができる
・非定常な波でも、ピーク周波数を計算できるものの、振幅・位相は、もはや信用できない
