Python -swanライブラリを用いたウェーブレット変換のコードの解説をお願いします。

現在、加速度の時系列データをpythonを用いてウェーブレット変換することを試みています。

前提・実現したいこと

現在、加速度の時系列データをpythonを用いてウェーブレット変換することを試みています。
swanライブラリを用いたウェーブレット変換を参考にしてまずはウェーブレット変換をするにあたっての大まかな流れをつかもうとしています。
コードの解読をしていますが、解説が少ない、ウェーブレット変換の理論をしっかりと理解していないこともあり何をしているのかわかりません。
ウェーブレット変換の理論を本やいろいろなサイトを見たりして勉強していますが、内容が難しく理解できていないので、コードの解説と同時に理論の説明もしていただけたら助かります。

pycwt.cwt_fの使い方がわかりません。ネットで調べたのですが参考になりそうなサイトを見つけることができませんでした。生成した信号に対して連続ウェーブレット変換をしているとは思うのですが、freqsというものがどういうものなのかがわからず困っています。pycwt.Morlet(omega0)はマザーウェーブレットを指定していると認識しています。
pycwt.cwt_fをr変数に格納していますが、このrはある時間に対しての周波数の値が格納されているということでしょうか。
ウェーブレット変換の理論、時系列データに対するウェーブレット変換を用いた実装について参考になりそうなサイト等ありましたら教えていただきたいです。

初学者のため、内容が初歩的な質問であると思いますが回答よろしくお願いいたします。

ソースコード

Python
1import matplotlib.pyplot as plt
2import numpy as np
3from swan import pycwt
4
5x = np.arange(0, 20, 0.01)
6y = np.sin(2 * np.pi * 1 * x) * 2 + np.sin(2 * np.pi * 2 * x) * 2  + np.cos(2 * np.pi * 10 * x)   
7
8plt.plot(x, y)
9plt.show()
10
11Fs = 1/0.01
12omega0 = 8
13# (1)　Freqを指定してcwt
14freqs=np.arange(0.1,20,0.025)
15r=pycwt.cwt_f(y,freqs,Fs,pycwt.Morlet(omega0))
16rr=np.abs(r)
17
18plt.rcParams['figure.figsize'] = (20, 6)
19fig = plt.figure()
20ax1 = fig.add_axes([0.1, 0.75, 0.7, 0.2])
21ax2 = fig.add_axes([0.1, 0.1, 0.7, 0.60], sharex=ax1)
22ax3 = fig.add_axes([0.83, 0.1, 0.03, 0.6])
23
24ax1.plot(x, y, 'k')
25
26img = ax2.imshow(np.flipud(rr), extent=[0, 20, freqs[0], freqs[-1]],
27                 aspect='auto', interpolation='nearest')
28
29fig.colorbar(img, cax=ax3)
30
31plt.show()

jbpb0

2020/12/03 09:21

> ウェーブレット変換の理論、時系列データに対するウェーブレット変換を用いた実装について参考になりそうなサイト等ありましたら教えていただきたいです。理論というか、変換した結果の見方は、このあたりが参考になるかもしれません http://www.kobe-kosen.ac.jp/~michi/akamatsu/1/2-2_the_comparison_between_wavelet_and_fourie_transform_using_test_waveform.html https://friedrice-mushroom.hatenablog.com/entry/2019/08/31/113915

AI_engineer

2020/12/07 06:45

jbpb0さん回答ありがとうございます。 URL見させていただきました。マザーウェーブレットの作成をする際に、 #マザーウェーブレットの作成(解析対象の波形用) m=6 #マザーウェーブレット(Morlet)のパラメータ、慣例的に6を使用 freq=20 #マザーウェーブレット(Morlet)の基本周波数[Hz] mw_otamesi=np.real(signal.morlet(int(fs/freq*m*2), m,1.0, complete=True)) という部分のコードなのですが、3行目の mw_otamesi=np.real(signal.morlet(int(fs/freq*m*2), m,1.0, complete=True)) が何をしているかわかりません。どっちかというとプログラミングというより物理よりの質問になってしまいますがもしよろしければ教えていただきたいです。よろしくお願いします。

jbpb0

2020/12/10 05:56

https://friedrice-mushroom.hatenablog.com/entry/2019/08/31/113915 でのそこの流れですが、otamesi が解析対象の1次元時系列データ(0～5秒)、mw が解析に使うウェーブレットです解析に使うウェーブレット mw は、scipy.signal.morlet で作られた Complex Morlet wavelet です https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.morlet.html mw=... の式を見たら、そのままですよね解析対象のデータ otamesi は、「たまたま※」データの真ん中付近(全5秒の2.5秒のところ)の形が Morlet wavelet の実部と同じで、それにノイズが加わっています mw_otamesi が、otamesi の「たまたま※」Morlet wavelet の実部と同じ形をしているところで、それを otamesi の真ん中にはめ込んでいますその後にノイズを加えて、otamesi の完成です以上が、「#解析対象の波形を作成」から「#解析対象の波形を作成おわり」までの間でやっていることです (※：もちろん、本当にたまたまなのではなく、意図して狙っているのですが) corr_tmp=np.correlate(... が、一つの周波数だけですが、ウェーブレット変換ですその下の数行は、otamesi と長さを揃えるために前後に0を足してるだけです numpy.correlate https://numpy.org/doc/stable/reference/generated/numpy.correlate.html によって、解析対象データ otamesi の0～5秒の間の各部分と、ウェーブレット mw の相互相関が計算されます otamesi の2.5秒のところは「たまたま」mw (の実部)と同じ形をしているので、とても相関が高くなりますそれ以外のところは単なるノイズなので、相関は低くなりますそれによって、otamesi の2.5秒付近には、mw の周波数と近い周波数の何かがある、ということが分かりますそのWebページの、そこのちょっと下の「ウェーブレット変換(スペクトログラム表示)」では、周波数をいろいろ変えながら上記と同じことをやってます list1 =... が解析するいろいろな周波数で、for freq in list1: のループが list1 から一つの周波数だけ取り出してウェーブレット変換をしてますループの中でやってることは、上記で書いたのと基本的には同じことですループ内でウェーブレット mw4 を作っている式の freq が毎回変わるので、作られるウェーブレットの周波数が変わりますそれとの相互相関を計算しているわけだから、いろいろな周波数の成分がそれぞれどれくらいあるのか、が分かりますある周波数のウェーブレットとの相関が、データのある位置で高ければ、データのその位置には、ウェーブレットと同じ周波数の成分がたくさんあるはず、ということです

AI_engineer

2020/12/11 05:48

詳しい回答ありがとうございます。 https://friedrice-mushroom.hatenablog.com/entry/2019/08/31/113915 でコードを読んでいる際にほかにもわからなかったところが出てきたのでよろしければ教えていただきたいです。 (1)https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.morlet.html を見させていただきました。 mw=signal.morlet(int(fs/freq*m*2), m,1.0, complete=True) はscipy.signal.morletを用いて実際に作られたマザーウェーブレットであることは理解できました。 signal.morlet関数の第一引数はウェーブレットの長さと書いてありますが、ウェーブレットの長さとは何を指しているのでしょうか。 https://friedrice-mushroom.hatenablog.com/entry/2019/08/31/113915 では20Hzのmorletがグラフで図示されていますが、この場合0の値をとっている部分を含めた0.6がウェーブレットの長さになっているのでしょうか。 int(fs/freq*m*2) の計算によってなぜウェーブレットの長さを計算することができるのかがわかりません。 (本などを読みましたがサンプリング周波数、基本周波数、omega0(角周波数)がどのようなものかいまいち理解が及んでいないためそこと関連付けて教えていただけたら幸いです) (2) (1)で述べたmwには実際にどのような値が格納されているのでしょうか。 morlet waveletをmwに格納していると理解していますが、具体的にはどのような数値が格納されているのかいまいちイメージがつかめません。なので次の行の #mw_x: 横軸(時間軸) mw_x = np.arange(0, len(mw)/fs, 1/fs) の第2引数もこの計算式で何をしているのかが理解できません。 (3) mw=signal.morlet(int(fs/freq*m*2), m,1.0, complete=True) corr_tmp=np.correlate(otamesi,mw)/sum(abs(mw)) の部分でnp.correlateでotamesi(解析する信号)とmw(マザーウェーブレット)の相関係数を計算していると思いますが、sum(abs(mw))で割っているのはどのような理由があるのでしょうか。初歩的な質問ばかりで申し訳ありませんが回答していただけたら幸いです。よろしくお願いいたします。

AI_engineer

2020/12/11 07:16

申し訳ありません。追記させてください。 (4) https://friedrice-mushroom.hatenablog.com/entry/2019/08/31/113915 のウェーブレット変換(スペクトグラム表示)のソースコードの一部に if freq==list1[0]: corr_stack=abs(corr_tmp[::10])**2 else: corr_stack=np.vstack((corr_stack, abs(corr_tmp[::10])**2)) とありますが、これも何を意味しているのかわかりません。corr_tmpのデータ形式はnumpy.ndarrayとなっていますがcorr_tmpにはotamesi(解析する信号)とmw4(マザーウェーブレット)の相関係数が1次元配列で格納されているのではないのでしょうか。 (5) イメージマップを出力する際の、 xsl = otamesi_x.T[::10] は何を意味していますか。私の理解が及ばず多くの質問をしてしまい申し訳ありません。回答できるものだけでもよいのでお力をお借りしたいです。

jbpb0

2020/12/13 09:22

(1) > この場合0の値をとっている部分を含めた0.6がウェーブレットの長さになっているのでしょうか。そうです 0.6秒間のデータ(mw)の個数が int(fs/freq*m*2) です print(len(mw)) print(int(fs/freq*m*2)) (2) > mwには実際にどのような値が格納されているのでしょうか。 print(mw) で表示されます (3) > sum(abs(mw))で割っているのはどのような理由があるのでしょうか。そのWebページの一番下の「フーリエ変換とウェーブレット変換の比較」のコードを、一旦そのまま実行してみてくださいただし、下記をグラフのプロットの前に挿入してください xmin=x.min() xmax=x.max() ymin=list1.start ymax=list1.stop Webページのようなグラフ・図が表示されたなら、次に、コード内の sum(abs(mw4)) で割っているのを止めるようにコードを修正してから、再度実行してみてください「Spectrogram(Wavelet)」の図がだいぶ変わりますので、割り算有り/無しの図を見比べて、割り算の意味をご自分で考えてください (4) > if freq==list1[0]: はループの初回なので、まだ corr_stack が存在してませんから、代入 > corr_stack=abs(corr_tmp[::10])**2 して corr_stack を作ります > else: はループの2回目以降なので、既に存在している corr_stack に追加 > corr_stack=np.vstack((corr_stack, abs(corr_tmp[::10])**2)) しますなお、corr_tmpに付いてる[::10]は、データを10個飛ばしにするという意味ですが、無くても表示される図の見た目は変わりません (5) > xsl = otamesi_x.T[::10] は何を意味していますか。 corr_tmpに合わせて、otamesi_xも10個飛ばしにしてます xminとxmaxを計算するのだけに使ってますが、corr_tmpと同様に10個飛ばしにせず、下記のようにしても、やはり図の見た目は変わりません xmin=otamesi_x.min() xmax=otamesi_x.max()