pythonで時間周波数解析~ウェーブレット変換~

以前、STFTを用いて時間周波数解析する記事を書きましたが、今回はウェーブレット変換を用いてやってみます。

連続ウェーブレット変換 (CWT) の絶対値を取ったものはスカログラムと呼ばれます。
スカログラムは、非周期信号に対して良い分解能を示します。
ウェーブレット変換の詳細は割愛しますので、Wikipedia等で見ていただければと思います。

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

PyWaveletsでウェーブレット変換

PyWaveletsはpythonでウェーブレット変換をするためのライブラリです。
pipコマンドでインストールできます。

pip install PyWavelets

1	pip install PyWavelets

ライブラリをインポートします。
PyWaveletsのライブラリ名はpywtです。

import matplotlib.pyplot as plt
import numpy as np
import pywt

import matplotlib.pyplot as plt

import numpy as np

import pywt

ここにあるようにしてwavファイルを読み込みます。
ウェーブレット変換してプロットします。

if __name__ == "__main__":
    wave, fs = wav_read(path_to_wavfile)
    wavelet_type = 'cmor1.5-1.0'
    scale = np.arange(1,257)
    coef, _ = pywt.cwt(wave, scales=scale, wavelet=wavelet_type)
    t = np.arange(len(wave))/fs
    frq = pywt.scale2frequency(scale=scale, wavelet=wavelet_type)*fs
    plt.pcolormesh(t, frq, 10*np.log(np.abs(coef)), cmap='jet')
    plt.colorbar()
    plt.show()
    plt.close()

if __name__ == "__main__":

wave, fs = wav_read(path_to_wavfile)

wavelet_type = 'cmor1.5-1.0'

scale = np.arange(1,257)

coef, _ = pywt.cwt(wave, scales=scale, wavelet=wavelet_type)

t = np.arange(len(wave))/fs

frq = pywt.scale2frequency(scale=scale, wavelet=wavelet_type)*fs

plt.pcolormesh(t, frq, 10*np.log(np.abs(coef)), cmap='jet')

plt.colorbar()

plt.show()

plt.close()

実行結果がこちらになります。
※見やすくするために少し間引いたりしています。

STFTで解析した場合と雰囲気は似ていますね。

TypeError: alpha must be a float or None

自分が実行した時に上記のエラーが吐かれました。
https://github.com/matplotlib/matplotlib/issues/17805
ウェーブレット変換を行うライブラリのファイルである cws.py を開いて、
COI_DEFAULTS = { 'alpha': '0.5', 'hatch':'/', }の'0.5'の部分を0.5に置き換えて保存すると正しく動作しました。