メインコンテンツへスキップ
見出し画像

書記が数学やるだけ#820 ウェーブレット解析

    ウェーブレット解析に必要な計算などをまとめておく。


    問題


    画像


    説明

    フーリエ変換では専ら周波数特性が解析され,時間領域の情報は残らなかった。これに窓関数をかける短時間フーリエ変換では,時間領域の解析が可能になるが,窓幅を周波数に合わせて固定する必要があるため,広い周波数領域の解析には向かなかった。これを解決するのがウェーブレット変換で,基底関数の拡大縮小を行うので,広い周波数領域の解析が可能となる。

    画像


    ウェーブレットとは波を拡大縮小・平行移動して足し合わせたものであり,これを用いて変換する。

    画像


    離散ウェーブレット変換について,特にハールウェーブレットが用いられる。

    画像


    解答

    ウェーブレット級数展開を以下に示しておく。

    画像


    ハールウェーブレットの概形を以下に示す,元の形に対し拡大縮小・平行移動を行うことができる。

    画像


    ハールウェーブレットの正規直交性について,内積の計算から示す。

    画像


    まずp>=2,つまりm≠kの場合でS=0であることを示す。

    画像


    次にqの値によりSの値を決める。

    画像


    ハールウェーブレットによるウェーブレット級数展開の係数を計算する。

    画像


    スカログラムとは信号の連続ウェーブレット変換 (CWT) の絶対値であり,以下のコードで可視化できる。

    import numpy as np
    import matplotlib.pyplot as plt
    import pywt
    import scipy
    
    # 与えられた信号
    t_start = 0
    t_end = 2
    num_points = 1000
    tm = np.linspace(t_start, t_end, num_points)
    amp = np.abs(np.sin(2.0 * np.pi * 2.0 * tm)) + np.tanh(tm)
    freq = 10.0 + 10.0 * tm * tm + np.cos(2.0 * np.pi * 3.0 * tm)
    sig = amp * np.sin(2.0 * np.pi * freq * tm) + amp * 5.0e-2 * np.random.randn(len(tm))
    
    # スケール毎の連続ウェーブレット変換
    wavelet = 'cmor'  # 連続モロー ウェーブレットを使用
    scales = np.arange(1, 128)  # 使用するスケールの範囲
    coefficients, frequencies = pywt.cwt(sig, scales, wavelet, sampling_period=1/(t_end - t_start))
    
    # スカログラムの表示
    plt.imshow(np.abs(coefficients), aspect='auto', extent=[t_start, t_end, scales[-1], scales[0]], cmap='jet', interpolation='bilinear')
    plt.colorbar(label='Magnitude')
    plt.title('Wavelet Transform Spectrogram')
    plt.xlabel('Time')
    plt.ylabel('Scale')
    plt.show()
    画像


    本記事のもくじはこちら:


    あなたへのおすすめ