フーリエが最良とは限らないかもしれない正弦波への分解
$${ s(t) = \displaystyle \sum_{k=1}^{K} a_k \sin(2\pi f_k t + \theta_k) }$$
という正弦波の和で書ける信号があったとして、$${s(t)}$$のデータから振幅、周波数、位相$${a_k , f_k , \theta_{k}}$$を計算したいとする
多分、現代の大学レベルの数学を習得した人の頭に最初に浮かぶのはフーリエ...だと思う
FFTは高速な実装も研究されてるし数値計算ライブラリに標準で備わってることも多いので使いやすい
$${N=2}$$で、$${f_1,f_2 = 20 \space \mathrm{Hz},30 \space \mathrm{Hz}}$$の場合、この分解を行うには、10Hz単位でスペクトルが分かればいいので、100msのデータがあればいいだろう
簡単なpythonプログラムを書いて実験してみる
import numpy as np
import matplotlib.pyplot as plt
freqs = [20,30]
Fs = 44000. #-- サンプリング周波数
times = np.arange(Fs * 0.1)
wavs = np.sin(2 * np.pi * freqs[0] * times / Fs) + 2.0 * np.sin(2*np.pi*freqs[1] * times / Fs)
wavs += 0.1 * ( np.random.rand( wavs.size ) - 0.5) #--add noise
plt.plot(times/Fs, wavs)
plt.title("raw signal")
plt.show()
max_freq = 150.0
specs = 2 * np.abs(np.fft.fft(wavs)) / wavs.size
freqs = np.fft.fftfreq(wavs.size, d = 1/Fs)
sz = np.sum(np.logical_and(freqs>0 , freqs<max_freq)) #-- 見やすくするために高い周波数をグラフから間引く
plt.bar(freqs[:sz],specs[:sz] )
plt.title("fourier amplitude")
plt.show()実行すると、以下の図が表示される
特に何も問題はない
$${N=2}$$で$${f_1,f_2 = 23 \space \mathrm{Hz} , 33 \space \mathrm{Hz}}$$の時はどうか
二つの周波数の差はさっきと同じ10Hzだけど、100msのデータではスペクトルを分離できず、1secのデータが必要になると思われる
また簡単なプログラムで実験してみる
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
freqs = [23,33]
Fs = 44000. #-- サンプリング周波数
times = np.arange(Fs * 1)
wavs = np.sin(2 * np.pi * freqs[0] * times / Fs) + 2.0 * np.sin(2*np.pi*freqs[1] * times / Fs)
wavs += 0.1 * ( np.random.rand( wavs.size ) - 0.5)
plt.plot(times/Fs, wavs)
plt.title("raw signal for duration=1sec")
plt.show()
#-- 1秒の信号からFFT
max_freq = 150.0
specs = 2 * np.abs(np.fft.fft(wavs)) / wavs.size
freqs = np.fft.fftfreq(wavs.size, d = 1/Fs)
sz = np.sum(np.logical_and(freqs>0 , freqs<max_freq))
plt.bar(freqs[:sz],specs[:sz] )
plt.title("fourier amplitude for duration=1sec")
plt.show()
#-- 100msの信号からFFT
wavs = wavs[:wavs.size//10]
max_freq = 150.0
specs = 2 * np.abs(np.fft.fft(wavs)) / wavs.size
freqs = np.fft.fftfreq(wavs.size, d = 1/Fs)
sz = np.sum(np.logical_and(freqs>0 , freqs<max_freq))
plt.bar(freqs[:sz],specs[:sz] )
plt.title("fourier amplitude for duration=100ms")
plt.show()
#-- 100msの信号からHann窓＋FFT
max_freq = 150.0
specs = 2 * np.abs(np.fft.fft(signal.hann(wavs.size)*wavs)) / wavs.size
specs *= (wavs.size / np.sum(signal.hann(wavs.size)))
freqs = np.fft.fftfreq(wavs.size, d = 1/Fs)
sz = np.sum(np.logical_and(freqs>0 , freqs<max_freq))
plt.bar(freqs[:sz],specs[:sz] )
plt.title("fourier amplitude for duration=100ms with hanning window")
plt.show()
100msのデータでは何となく、20Hzから40Hzのあたりにピークがありそうな雰囲気はあるが、23Hzと33Hzという数字は取ってこれない
(工学的用途で使うような)人工的に生成される信号なら、周波数の決定権はこっちにあることも多く、便利な周波数を選べばいいかもしれない
一方で、自然現象を解析したい場合は、周波数成分が未知のまま、信号を分析しないといけないかもしれないし、信号自体に綺麗な周期性が見えるとは限らない
原理的には十分長い時間の信号データを持ってきてFFTすればいいが、一般的には、信号成分がいつまでも定常的とは限らない(つまり、振幅や周波数が定数と見なして問題ないのは、比較的短時間である可能性がある)ので、できるだけ短時間の信号から周波数や振幅を抽出できた方がいい
そんな問題意識を持ちつつ、改めて冒頭の式
$${ s(t) = \displaystyle \sum_{k=1}^{K} a_k \sin(2\pi f_k t + \theta_k) }$$
を眺めていると、こいつが定数係数常微分方程式の解になってることに気付く
幸運なことに定数係数常微分方程式は(等間隔で)差分化しても、定数係数差分方程式になる。つまり
$${s( n \cdot \delta t) = s_n }$$
として
$${\displaystyle \sum_{j=0}^{J} c_{j} s_{n-j} = 0}$$
を満たす定数$${c_j}$$が存在する。$${c_0 \neq 0}$$であるが、この差分方程式は
$${ s_{n} = -\displaystyle \sum_{j=1}^{J} \dfrac{c_{j}}{c_{0}} s_{n-j} }$$
と変形すれば過去のデータから現在の値を予測する式とも解釈できる。が、予測式としての解釈は今は重要じゃない
一般に定数係数常微分方程式の解は$${\exp(\alpha t)}$$という形だが、$${|\alpha|=1}$$は保証されないので、絶対値が1でなければ減衰もしくは指数関数的成長を表す項が出る
差分方程式の階数$${J-1}$$が既知なら、Hankel行列
$${H_{i,j} = s_{i+j-1}}$$ $${\space}$$ $${(j=1,\cdots,J)}$$
を作って(正方行列とは限らないが)、そのカーネルを計算すればいい
現実のデータでは、諸要因によって、必ず信号に誤差が入ってくるので、厳密なカーネルは存在しないが、SVDなどによって、最小二乗の意味で解を見つけることはできるだろう
大きな問題は、階数$${J-1}$$が一般に全然既知じゃない。それどころか、例えば矩形波なら原理的には無限個の正弦波を含む
サンプリング周波数の半分より上は見なくていいが、データ量の半分より大きな階数を取るとHankel行列の行数が列数より小さくなるので、データ量の半分に設定しておけばそれ以上は望めない
実験してみると分かるけど、ノイズが乗ってるような普通のデータに対しては、小さい階数では正しい答えが出ないことが殆どなので、大きな階数を取ることは必要でもある
もう一つの問題は計算量だけど、できる限り短時間のデータから正弦波成分を抽出したろうというのが目標なので、短時間のデータなら、問題ではないかもしれない
そんなにうまくいくのか先ほどと同じデータで実験
import numpy as np
import matplotlib.pyplot as plt
import math
freqs = [23,33]
Fs = 44000. #-- サンプリング周波数
times = np.arange(Fs * 0.1)
wavs = np.sin(2 * np.pi * freqs[0] * times / Fs) + 2. * np.sin(2 * np.pi * freqs[1] * times / Fs)
wavs += 0.1 * ( np.random.rand( wavs.size ) - 0.5) #--add noise
rk = wavs.size//2 - 1
hank = np.zeros((rk, wavs.size-rk))
for n in range(wavs.size-rk):
hank[: , n] = wavs[n:n+rk]
u,s,vh = np.linalg.svd(hank)
ans = u[:,-1]
freqlist = np.array([Fs * math.atan2(c.imag, c.real) / (2*math.pi) for c in np.roots(ans)])
print( freqlist.flat[np.abs(freqlist-freqs[0]).argmin()] )
print( freqlist.flat[np.abs(freqlist-freqs[1]).argmin()] )実行すると22.983563155912886と33.00889580388116が出るので、23Hzと33Hzが含まれてそうなことは確認できる
が、候補が沢山ありすぎて、どれが主要な周波数か調べるのが面倒くさい
特異値の方を大きい順に見ると、2072.02049189, 1658.40226164, 789.37354323, 728.30902977, 3.98187102 ,...となって最小値は4.65e-3程度で、小さい特異値が大量に並ぶ
大きい4つの特異値は23Hzと33Hzの成分に対応していると考えられ(負の周波数もあるので2倍の4つになる)、小さい特異値成分からの寄与を無視することで、denoisingができると思われる。が、今は別にノイズを除去したいわけじゃない
Hankel行列には、特殊な行列の分解が存在する。$${s_n}$$をsin関数の線形結合で書くのは理論的には不便なので指数関数を代わりに使って
$${s_{n} = \displaystyle \sum_{k=1}^{K} r_{k} z_{k}^n }$$
と書けるとすると
$${ H = \left( \begin{matrix} 1 & \cdots & 1 \\ z_1 & \cdots & z_K \\ & \cdots & \\ z_1^{N-J} & \cdots & z_K^{N-J} \end{matrix} \right) \left( \begin{matrix} z_1 r_1 & \cdots & z_1^{J} r_1 \\ & \cdots & \\ z_K r_K & \cdots & z_K^{J} r_K \end{matrix} \right) }$$
はすぐ分かる
$${\mathbf{r} = (r_{1} , \cdots , r_{K})^{T}}$$
$${Z = \mathrm{diag}(z_1 , \cdots , z_{K})}$$
と置くと
$${ H = \left( \begin{matrix} 1 & \cdots & 1 \\ z_1 & \cdots & z_K \\ & \cdots & \\ z_1^{N-J} & \cdots & z_K^{N-J} \end{matrix} \right) \left( \begin{matrix} Z \mathbf{r} & \cdots Z^{J} \mathbf{r} \end{matrix} \right) }$$
のように書ける。これを
$${H = H_{L} H_{R}}$$
と表すことにする
$${K \times K}$$可逆行列$${Q}$$を使って
$${\left( \begin{matrix} 1 & \cdots & 1 \end{matrix} \right) \mapsto \left( \begin{matrix} 1 & \cdots & 1 \end{matrix} \right)Q }$$
$${Z \mapsto Q^{-1} Z Q }$$
$${\mathbf{r} \mapsto Q^{-1} \mathbf{r}}$$
のようにすると、
$${ H_L H_R \mapsto (H_{L} Q)(Q^{-1} H_{R}) = H }$$
で、別の分解を得られる
truncateしたSVD
$${H \approx U_{K} \Sigma_{K} V_{K}^{H} }$$
と比べる($${\Sigma_K}$$は特異値を大きい順に$${K}$$個並べた対角行列で、$${U_K}$$と$${V_K}$$は残った特異値に対応する成分をとりだしている)と、"genericには"
$${U_{K} \approx H_{L} Q}$$
となる可逆行列$${Q}$$の存在が分かる
$${K}$$をいくつにするかは特異値を見てアタリを付けることができる
$${H_{L}}$$から一行目と最後の$${N-J}$$行目を除いた行列を$${H_{L,bot}}$$と$${H_{L,top}}$$とすると
$${H_{L,top} Z = H_{L,bot}}$$
だから、$${Z' = Q^{-1} Z Q}$$に対して
$${U_{K,top} Z' \approx U_{K,bot}}$$
となる。$${U_{K,top}}$$と$${U_{K,bot}}$$は行列$${U_{K}}$$から最後の行と一行目を除いた行列
この方法は1987年の論文Improved algorithm for noniterative time-domain model fitting to exponentially damped magnetic resonance signalsで提案されHSVDと呼ばれている。HはHankelの頭文字らしい
今までと同じデータでうまくいくか実験する
import numpy as np
import matplotlib.pyplot as plt
import math
freqs = [23,33]
Fs = 44000. #-- サンプリング周波数
times = np.arange(Fs * 0.1)
wavs = np.sin(2 * np.pi * freqs[0] * times / Fs) + 2. * np.sin(2 * np.pi * freqs[1] * times / Fs)
wavs += 0.1 * ( np.random.rand( wavs.size ) - 0.5) #--add noise
rk = wavs.size//2 - 1
hank = np.zeros((rk, wavs.size-rk))
for n in range(wavs.size-rk):
hank[: , n] = wavs[n:n+rk]
u,s,vh = np.linalg.svd(hank)
K=4
U_K_top = u[:-1,:K]
U_K_bot = u[1: ,:K]
Zp=np.dot(np.linalg.pinv(U_K_top) , U_K_bot)
for z in np.linalg.eig(Zp)[0]:
print(Fs*math.atan2(z.imag,z.real)/(2*math.pi))出力は
32.99490887189571
-32.99490887189571
23.019053467921363
-23.019053467921363
のようになる
無事(?)、±33Hzと±23Hzが得られる
$${Z'}$$の固有値は絶対値が1であるという保証はないので、そのへんの確認は本来は飛ばすべきではない
$${Z'}$$の固有値が1でない場合も含めて考えると、この方法は、一般には減衰項付きの
$${ s(t) = \displaystyle \sum_{k=1}^{K} e^{-b_k t} a_k \sin(2\pi f_k t + \theta_k) }$$
という形の表示を探索していることになる
まだ周波数しか得てないが、振幅と位相は、線形結合の結合係数を計算するだけの問題なので、これも線形代数でできる
効率的とはいえないが、全部まとめると、以下のような実装になる(簡単のため、$${K}$$は決め打ちしているけど、本来は特異値を見て決める処理を挟むべき。但し、これを決めるための正解というのはないので、こういう実装にしている)
import numpy as np
import matplotlib.pyplot as plt
import math
def hsvd(dat , Fs , K):
rk = dat.size // 2 -1
hank = np.zeros((rk, dat.size-rk))
for n in range(dat.size-rk):
hank[: , n] = dat[n:n+rk]
u,s,vh = np.linalg.svd(hank)
U_K_top = u[:-1,:K]
U_K_bot = u[1: ,:K]
Q = np.dot(np.linalg.pinv(U_K_top) , U_K_bot)
freqlist = []
for z in np.linalg.eig(Q)[0]:
#-- 本来はzの絶対値がほぼ1であるかどうかチェックしたほうがいい
f0 = Fs*math.atan2(z.imag,z.real)/(2*math.pi)
if f0>0:freqlist.append(f0)
return freqlist
#-- 周波数が分かったら、振幅と位相を計算する
def sosd(dat , Fs, freqs):
times = np.arange(dat.size) / Fs
components = [np.ones(times.size)]
for f in freqs:
components.append( np.cos(2*np.pi*f*times) )
components.append( np.sin(2*np.pi*f*times) )
c = np.stack(components).T
coeff = np.dot( np.linalg.pinv(c) , dat)
return coeff
if __name__=="__main__":
freqs = [23,33]
Fs = 44000. #-- サンプリング周波数
times = np.arange(Fs * 0.1)
wavs = np.sin(2 * np.pi * freqs[0] * times / Fs) + 2. * np.sin(2 * np.pi * freqs[1] * times / Fs)
wavs += 0.1 * ( np.random.rand( wavs.size ) - 0.5) #--add noise
freqlist = hsvd(wavs, Fs , K=4)
print( freqlist )
print( sosd(wavs , Fs, freqlist) )というわけで、100msのデータで23Hzと33Hzの正弦波が発見できるようになった
もっと短い信号でいけるかというと、直感的に、10Hzの差を区別するには、ある程度の長さのデータは必要だと思われる
例えば、30Hzの正弦波と31Hzの正弦波を一周期で判別できるかというと、ノイズが一切ないなら可能かもだが、ノイズがあったら、どこが周期点なのか正確に特定することは難しい
一方で、もし一秒分のデータがあるなら、波の数が30個か31個かという違いはかなり明瞭に区別できると思われる
実験してみると、50msのデータでも10Hzの差を検出できないことはなさそう
それを下回ると信頼度は落ちていく。FFTの時と違って、明確な閾値はないが
FFTよりも計算量が多いので、FFTを過去にするほど優れてるってわけじゃないが、使い道はあるだろう