pythonで瞬時周波数推定~ヒルベルト変換~

今回は信号の瞬時周波数の推定(と、位相推定)をします。
用いるのはヒルベルト変換です。

前回と同様に、音源はESC-50というデータセットからお借りしました。
データセットのリポジトリはこちら➡https://github.com/karolpiczak/ESC-50
siren(サイレン)のタグがついたファイルを借りました➡1-76831-A-42.wav

scipy.signal.hilbertでヒルベルト変換

scipyの中にヒルベルト変換用のライブラリがあるのでそれを用いて変換します。
numpyのunwrap関数を用いて瞬時位相を求め、それを用いて瞬時周波数を求めます。

if __name__ == "__main__":
    wave, fs = wav_read(path_to_wavfile)
    z = signal.hilbert(wave)
    phase = np.unwrap(np.angle(z))
    frq = np.diff(phase) / (2 * np.pi) * fs
    plt.subplot(2,1,1)
    plt.title("Instantaneous Frequency")
    plt.plot(t[1:],frq)
    plt.subplot(2,1,2)
    plt.title("Instantaneous Phase")
    plt.plot(t,phase)
    plt.show()
    plt.close()

if __name__ == "__main__":

wave, fs = wav_read(path_to_wavfile)

z = signal.hilbert(wave)

phase = np.unwrap(np.angle(z))

frq = np.diff(phase) / (2 * np.pi) * fs

plt.subplot(2,1,1)

plt.title("Instantaneous Frequency")

plt.plot(t[1:],frq)

plt.subplot(2,1,2)

plt.title("Instantaneous Phase")

plt.plot(t,phase)

plt.show()

plt.close()

出力結果はこんな感じです。

上側が瞬時周波数、下側が瞬時位相になります。
少しガタガタしてますね。

STFT、ウェーブレット変換との比較

同じ信号について、それぞれの手法で時間周波数解析を行っていたのでそれに重ねてプロットします。
ガタガタし過ぎていたので適当なフィルタをかけています。
STFTの場合。

ウェーブレット変換の場合。

どちらの手法でも一番濃ゆいところに大体合うように推定できていますね。
今回のサイレンの音のように音源が少ない時には、時間周波数解析をしなくても情報が十分得られるかもしれないと思いました。