連続ウェーブレット変換

Pythonで連続ウェーブレット変換を試みたことのまとめ。

背景 フーリエ変換について

ある音声データについて、どういった周波数成分が含まれているのか?を調べる方法として、まず最初に思いつくのはフーリエ変換だと思う。 フーリエ変換をざっくり説明すると、音声データのある区間についていろんな周波数のサイン波との相関を計算していくことで、各周波数ごとの強弱を調べることができるというもの。ちなみに下図はサイン波、矩形波、ノコギリ波のスペクトルを求めたもの。左側が波形、右側が周波数スペクトルとよばれるもので、横軸を周波数、縦軸を振幅でプロットしたもの。

このようなフーリエ変換を応用したもので一般的なのは下図のようなグラフだと思う。横軸を時間、縦軸を周波数、強度を色(青→弱、赤→強)で表したもので、スペクトログラムと呼んでいる。

ウェーブレットについて

フーリエ変換によるスペクトログラムは頻繁に使われているが、時間精度に限界があるという弱点がある。 というのも、スペクトログラムの作るために音声データをある時間間隔ごとに細かく区切ることから始まる。区切ったあと、それぞれにフーリエ変換を行い、それぞれ周波数スペクトルを導出してから時間・周波数・強度のグラフに再構成していく(厳密には単純に区切るだけでなく、オーバーラップもしているが)。つまり、時間に対する精度は音声信号を区切る間隔に制約される。

フーリエ変換ばかり述べているが、ここで対抗馬として挙がってくるのがウェーブレット変換という手法。 ウェーブレットはwave(波)+let(小さい)の合成語(Wikipedia引用)で、下図のように始めと終わりでフェードイン・アウトしたような小さな波形のことを指していて、マザーウェーブレットと呼んでいる。様々なものがあるが下図はその一例で、モルレーのウェーブレット(Morlet wavelet)と呼ばれるもの。個人的に馴染みがあるのがコレだけなので、以後もコレを使用していく。 ちなみに、青色とオレンジ色の線があるのは実部と虚部があるため。虚部は実部に対して90°遅れていて、例えば相対的な位相差を見ることができる。

(下図 マザーウェーブレットの例 Morlet wavelet, m=6)

上図をプロットするためのpythonコード(Jupyter notebookを推奨) この儀式は以下のPythonコードすべてに共通するため、ここでまとめた書いた。

#儀式
%matplotlib notebook

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 儀式その2 plotの体裁を整える https://qiita.com/rnhd/items/325d21621cfe9e553c17
plt.rcParams['font.family'] ='sans-serif'#使用するフォント
plt.rcParams['xtick.direction'] = 'in'#x軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams['ytick.direction'] = 'in'#y軸の目盛線が内向き('in')か外向き('out')か双方向か('inout')
plt.rcParams['xtick.major.width'] = 1.0#x軸主目盛り線の線幅
plt.rcParams['ytick.major.width'] = 1.0#y軸主目盛り線の線幅
plt.rcParams['font.size'] = 8 #フォントの大きさ
plt.rcParams['axes.linewidth'] = 1.0# 軸の線幅edge linewidth。囲みの太さ

本チャンはこちら。

#マザーウェーブレットの作成
fs=1000 # サンプリング周波数を定義[Hz]
m=6 #マザーウェーブレット(Morlet)のパラメータ、慣例的に6を使用
freq=20 #マザーウェーブレット(Morlet)の基本周波数[Hz]

#mw: mother wavelet
mw=signal.morlet(int(fs/freq*m*2), m,1.0, complete=True)

#mw_x: 横軸(時間軸)
mw_x = np.arange(0, len(mw)/fs, 1/fs) 

#マザーウェーブレットのプロット
fig =plt.figure(figsize=(8,3)) #plot in new figure
fig.subplots_adjust(bottom=0.2) #xlabelがはみでてしまうので、subplotのサイズ調整

plt.plot(mw_x.T,np.real(mw),label="mw-real")
plt.plot(mw_x.T,np.imag(mw),label="mw-imag")
plt.xlabel("Time[s]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.title("morlet (m=6)")

plt.legend()

マザーウェーブレットは一定の周期のため、横幅(スケール)を伸ばしたり縮めたりすることで任意の周波数成分をもつ波形に変形することができる。 下図は20Hzの周波数成分をもつマザーウェーブレットと40Hzのものを並べたもの。一見同じ波形に見えるが、時間軸が0〜0.6秒と0〜0.3秒なので2倍異なってるよ!といっている。

上図をプロットするためのpythonコード(Jupyter notebookを推奨)

fs=1000 #サンプリング周波数[Hz]
m=6 #マザーウェーブレット(Morlet)のパラメータ、慣例的に6を使用

#波形1(基本周波数=20Hz)
freq=20 #マザーウェーブレット(Morlet)の基本周波数[Hz]
mw=signal.morlet(int(fs/freq*m*2), m,1.0, complete=True)
mw_x = np.arange(0, len(mw)/fs, 1/fs) 

fig =plt.figure(figsize=(8,3)) #plot in new figure
fig.subplots_adjust(bottom=0.2) #xlabelがはみでてしまうので、subplotのサイズ調整

plt.subplot(1,2,1)
plt.plot(mw_x.T,np.real(mw),label="mw-real")
plt.plot(mw_x.T,np.imag(mw),label="mw-imag")
plt.xlabel("Time[s]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.title("morlet (20Hz)")
plt.legend()


#波形2(基本周波数=40Hz)
freq=40
mw=signal.morlet(int(fs/freq*m*2), m,1.0, complete=True)
mw_x = np.arange(0, len(mw)/fs, 1/fs) 

plt.subplot(1,2,2)
plt.plot(mw_x.T,np.real(mw),label="mw-real")
plt.plot(mw_x.T,np.imag(mw),label="mw-imag")
plt.xlabel("Time[s]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.title("morlet (40Hz)")
plt.legend()

ウェーブレット変換(単一の周波数解析)

ようやくここから本題……。 スペクトログラムの弱点である時間精度の低さに対して、ウェーブレット変換が優位に立つことができるのは、このマザーウェーブレットを使用する点にある。 例えばある音声データについて、特定の周波数(たとえば20Hz)の強度がどのように変化していくかを調べるとする。ウェーブレット変換では、まずマザーウェーブレットについて周波数が20Hzとなるようにスケール(マザーウェーブレットの幅)を調整しておく。そして、それを音声データに対して畳み込み演算を行っていくことで、その周波数の強度の経時変化を調べることができる。

畳み込み演算は、関数A(長い方)に関数B(短い方)を、少しずつズラしながら積分する手法のこと。内部的には2つの関数が重なった部分の(正負の)面積を計算していて、2つの関数が似ていれば畳み込みの値は大きくなり、相関が強いと解釈することができる。一方で2つの関数が全く似通っていなければ、畳込みの値は正負が入り混じって総和が小さくなるため、相関が弱いと解釈することができる。 下図は、あるサンプル波形にウェーブレット変換を試みたもの。 上段はサンプルとして作成した20Hzの周波数成分をもつ波がピョっとでる波形、下段はサンプル波形と20Hzのマザーウェーブレット(morlet)とのたたみこみ演算を行い、相関の経時変化を示したもの。 下図からわかるように2.5秒あたりでピークをもっていて、それは解析対象の波形とマザーウェーブレットがピッタリ重なるのがこのタイミングだからに他ならない。

上図をプロットするためのpythonコード(Jupyter notebookを推奨)

#解析対象の波形を作成

fs=1000 #サンプリング周波数[Hz]
s_length=5 # 解析対象波形の長さ[秒]

#otamesi: 解析対象の波形 まずゼロ配列を作成、サンプリング周波数 1kHz、5秒の信号を想定
otamesi = np.zeros(s_length*fs) 
#otamesi_x: 解析対象波形のサンプル点と時間[秒]を対照するための配列
otamesi_x = np.arange(0, s_length, 1/fs) 

#マザーウェーブレットの作成(解析対象の波形用)
m=6 #マザーウェーブレット(Morlet)のパラメータ、慣例的に6を使用
freq=20 #マザーウェーブレット(Morlet)の基本周波数[Hz]
mw_otamesi=np.real(signal.morlet(int(fs/freq*m*2), m,1.0, complete=True))

#先に作ったゼロ配列に、マザーウェーブレットを合成(2.5秒をピークにすることを想定)
otamesi[int(fs*s_length/2-len(mw_otamesi)/2):int(fs*s_length/2+len(mw_otamesi)/2)]=mw_otamesi
otamesi+=np.random.randn(s_length*fs)*0.01 #最大振幅±0.01のホワイトノイズを乗せる
#解析対象の波形を作成 おわり

#プロット
fig =plt.figure(figsize=(8,4)) #plot in new figure
fig.subplots_adjust(bottom=0.2) #xlabelがはみでてしまうので、subplotのサイズ調整

plt.subplot(2,1,1)
plt.plot(otamesi_x.T,otamesi,linewidth=1)
plt.title("analyzing_waveform")
plt.xlabel("Time[s]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.xlim([0,5])

#マザーウェーブレット
m=6 #マザーウェーブレット(Morlet)のパラメータ、慣例的に6を使用
freq=20 #マザーウェーブレット(Morlet)の基本周波数[Hz]


mw=signal.morlet(int(fs/freq*m*2), m,1.0, complete=True)
corr_tmp=np.correlate(otamesi,mw)/sum(abs(mw))
sup_first=np.zeros(int(len(mw)/2))
sup_end=np.zeros(fs*s_length-len(corr_tmp)-int(len(mw)/2))
corr_tmp=np.append(sup_first,corr_tmp)
corr_tmp=np.append(corr_tmp,sup_end)

plt.subplot(2,1,2)
#plt.plot(otamesi_x.T,abs(corr_tmp),linewidth=1)
plt.plot(otamesi_x.T,(abs(corr_tmp))**2,linewidth=1)
plt.title("wavelet_coherence@20Hz")
plt.xlabel("Time[s]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.xlim([0,5])

fig.tight_layout() 

ウェーブレット変換(スペクトログラム表示)

そして、今度はさきほどのサンプル波形について複数の周波数でウェーブレット変換を行い、スペクトログラム風に表示してみる。 下図は、マザーウェーブレットの基本周波数を10〜40Hzまで1Hz刻みで変化させて、相関の経時変化を計算していった結果。20Hzかつ2.5秒のところでピークを持っていることに変わりないことが示されている。

上図をプロットするためのpythonコード(Jupyter notebookを推奨)

m=6
#解析する周波数帯[Hz]
list1 = range(10, 50, 1)

print("start")
for freq in list1:
    print(freq,end="")
    print('-',end="")
    mw4=signal.morlet(int(fs/freq*m*2), m,1.0, complete=True)
    corr_tmp=np.correlate(otamesi,mw4)/sum(abs(mw4))
    
    sup_first=np.zeros(int(len(mw4)/2))
    sup_end=np.zeros(fs*s_length-len(corr_tmp)-int(len(mw4)/2))
    corr_tmp=np.append(sup_first,corr_tmp)
    corr_tmp=np.append(corr_tmp,sup_end)

    if freq==list1[0]:
        corr_stack=abs(corr_tmp[::10])**2
    else:
        corr_stack=np.vstack((corr_stack, abs(corr_tmp[::10])**2))
print()
print("end")

#イメージマップ出力の段

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec


xsl = otamesi_x.T[::10]
xmin=xsl.min()
xmax=xsl.max()
ymin=list1.start
ymax=list1.stop

fig =plt.figure(figsize=(8,4))

gs = gridspec.GridSpec(2,1)
ax1 = fig.add_subplot(gs[0,:])

plt.plot(otamesi_x.T,otamesi,linewidth=1)

plt.title("sample_waveform")
plt.xlabel("Time[s]")
plt.ylabel("Amplitude")


plt.grid(True)
ax1.axis('tight')

plt.xlim([0,5])

ax2 = fig.add_subplot(gs[1,:])

im = ax2.imshow(corr_stack,extent=[xmin,xmax,ymin,ymax],interpolation='none',cmap='jet',origin='lower')

plt.title("wavelet_coherence")
plt.xlabel("Time[s]")
plt.ylabel("Frequency[Hz]")
plt.grid(linestyle='--', linewidth=0.5)


# カラーバーの追加をしたい ↓参考
# https://matplotlib.org/2.0.2/examples/pylab_examples/colorbar_tick_labelling_demo.html

cmlb=[0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1]
cbar = fig.colorbar(im, ticks=cmlb, orientation='horizontal')
cbar.ax.set_xticklabels(cmlb)  # horizontal colorbar


ax2.axis('tight')
plt.xlim([0,5])

フーリエ変換とウェーブレット変換の比較

最後に、フーリエ変換とウェーブレット変換それぞれによるスペクトログラムと比較してみる。 下図がその結果で、上段はサンプル波形、下段の左側がウェーブレット変換、右側がフーリエ変換によるスペクトログラム。 サンプル波形はべつのものを用意した。0〜5秒までで周波数が0〜150Hzまでスウィープするサイン波に、振幅1、2Hzのサイン波をかけ合わせて、1秒毎に包絡線がゼロになるようにしたもの。

この比較でウェーブレット変換とフーリエ変換それぞれの特徴がでていると思う。ウェーブレット変換ではサンプル波形の包絡線がゼロになっているのがくっきり見えているのに対して、フーリエ変換ではその境界が曖昧になっている。ただし、解析する時間の区切りかたやオーバーラップさせる範囲の組み合わせで、用途によって調整していくことはできると思う。

一方、フーリエ変換については周波数の分解能が一定であるのに対して、ウェーブレット変換では高周波になるほど分解能が低くなっている。 ただし、これも使いかたによると思う。というのも、人間の知覚に関する法則で「ウェーバー・フェヒナーの法則」というものがあり、人間が知覚する刺激の強さは実際の刺激の強さの対数に比例するというもの。 たとえば人間が1,2,3,4……と段階的に刺激を知覚している場合、実際の刺激は1,2,4,8…と指数関数的に強度を高めている場合が多い。音階が1オクターブ上がるごとに周波数が2倍になるのが最たる例で、440±1Hzの誤差がギリギリ弁別できたとしても、その2オクターブ上の1760±1Hzの誤差を知覚するのは神業に等しい。だから、高周波になるほど周波数分解能が低くなるといっても、人間が知覚できる分には問題ないので問題ないとも言えるので、結局これも用途による。

さて、これはなんの図だったか…あ、フーリエ変換とウェーブレット変換のスペクトログラムだった。

上図をプロットするためのpythonコード(Jupyter notebookを推奨)

fs=1000 # Hz
s_length=5 # second

x = np.arange(0, s_length, 1/fs) #サンプリング周波数 1kHz、5秒の信号を想定
y=np.sin(np.pi*x**2*30)*np.sin(np.pi*x*1)+np.random.randn(s_length*fs)*0.01
#0~5秒で周波数が0〜150Hzにスウィープする関数。ホワイトノイズも付与。

#マザーウェーブレット
m=6 #6がオススメらしいので

#解析する周波数帯[Hz]
list1 = range(2, 200, 1)

print("start")
for freq in list1:
    print(freq,end="")
    print('-',end="")
    mw4=signal.morlet(int(fs/freq*m*2), m,1.0, complete=True)
    corr_tmp=(np.correlate(y,mw4)/sum(abs(mw4)))
    
    sup_first=np.zeros(int(len(mw4)/2))
    sup_end=np.zeros(fs*s_length-len(corr_tmp)-int(len(mw4)/2))
    corr_tmp=np.append(sup_first,corr_tmp)
    corr_tmp=np.append(corr_tmp,sup_end)

    if freq==list1[0]:
        corr_stack=abs(corr_tmp[::10])
    else:
        corr_stack=np.vstack((corr_stack, abs(corr_tmp[::10])))
print()
print("end")


fig =plt.figure(figsize=(8,5))

gs = gridspec.GridSpec(3,2)
ax1 = fig.add_subplot(gs[0,:])

plt.plot(x,y,linewidth=0.5)

plt.title("sample_waveform_2")
plt.xlabel("Time[s]")
plt.ylabel("Amplitude")
plt.grid(True)

ax1.axis('tight')
plt.xlim([0,5])


ax2 = fig.add_subplot(gs[1:,0])

im = ax2.imshow(corr_stack,extent=[xmin,xmax,ymin,ymax],interpolation='none',cmap='jet',origin='lower')

plt.title("Spectrogram(Wavelet)")
plt.xlabel("Time[s]")
plt.ylabel("Frequency[Hz]")
plt.grid(linestyle='--', linewidth=0.5)

cmlb=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1]
cbar = fig.colorbar(im, ticks=cmlb,orientation='horizontal')
cbar.ax.set_xticklabels(cmlb) 
ax2.axis('tight')
plt.xlim([0,5])
plt.ylim([0,200])
###################

ax3 = fig.add_subplot(gs[1:,1])
NFFT = 1024 # the length of the windowing segments
Fs = fs # the sampling frequency
plt.specgram(y,NFFT=NFFT, Fs=Fs, noverlap=1000,cmap='jet',scale='linear')
plt.grid(linestyle='--', linewidth=0.5)

plt.title("spectrogram(FFT)")
plt.xlabel("Time[s]")
plt.ylabel("Frequency[Hz]")

plt.colorbar(orientation='horizontal')

ax3.axis('tight')
plt.xlim([0,5])
plt.ylim([0,200])
fig.tight_layout() 

結論

フーリエ変換もウェーブレット変換も、用途によって使い分けるよという結論になってしまった。とくにウェーブレットを布教したいという意図もない(メモランダムなので。。。)

以前はMATLABでWavelet変換を使っていて、そのときはすべてツールボックスがお膳立てをしてくれていたので細かいことを気にしていなかった。 それから数年、その当時の記憶が消失。個人でMATLABなどもってないし、今はなんかPythonが主流なので、思い出しも兼ねて当時の解析方法をPythonで再現することを試みてみた次第。 このやり方で良かったのか曖昧なところがあるので、自信がないところは露骨にボカしている(露骨にボカすとは)。 そういうところがたくさんあるため、もし参考にしようとする奇特な方におかれましては、あくまで参考の一つにされることを強く推奨します。