Juliaで離散ウェーブレット変換を試してみる(備忘録)

導入

最近ウェーブレット変換について勉強しているが、いまいちウェーブレット変換の操作について実感がわかない。
そこでjuliaを使いながら(連続ウェーブレット変換の離散化ではなく)離散ウェーブレット変換を実行し、雰囲気をつかんでいきたい。

離散ウェーブレット変換

何やら画像処理や信号処理に使用されるらしい。
この分野の教科書は難しいものが多く(ドブシーとか)、なかなか前に進んでいかない感じがする…

さて、まずはN=8(サンプル数をＮとした)の場合に手計算で離散ウェーブレット変換を実行してみる。
基底として、
$$
\ket{\phi_b;t,n}
$$
というものを用意する。これを基本単位ステップという。
ここで$b$は信号のシフト量、$t$は信号の幅、$n$は始点とする。

図で示すとこんな感じ。
基底で表す信号の大きさは1か0の2値をとることにする。

この表現を使うことで、例えば
$$
\ket{\phi_0;4,1} = \ket{\phi_0,2,1} +\ket{\phi_0,2,3} = \ket{\phi_0,2,1}+\ket{\phi_2,2,1}
$$
などと表すことができる。この信号はサイト1から始まる幅4の信号で、サイト1から始まる幅2の信号と
サイト3から始まる幅2の信号の和で表される。もちろん、$\sum_{n=1}^4\ket{\phi_0;1,n}$と表現することもできる。
さて任意の信号をこれらの基底を使って表したい。しかしながら、
$\ket{\phi_0;4,1} = \ket{\phi_0,2,1} +\ket{\phi_0,2,3}$
という表現では基底の数が足りず、任意の信号を表せない。
そこで
$$
\ket{\psi_0;4,1} = \ket{\phi_0,2,1} - \ket{\phi_0,2,3}
$$
という基底の表現を考える。
これによって2つの基底が2つの基底に変換されるので完全である。
この$\ket{\psi_b;t,n}$を基本ウェーブレットという。

具体例として、あるベクトル(データ)$\ket{s}=s_1\ket{\phi_0,2,1} +s_2\ket{\phi_0,2,3}$
を基本単位ステップ$\ket{\phi_0;4,1}$と基本ウェーブレット$\ket{\psi_0;4,1}$を使って表すことを考えてみる。容易に
$$
\ket{s} = \frac{s_1 + s_2}{2}\ket{\phi_0;4,1} +\frac{s_1 - s_2}{2}\ket{\psi_0;4,1}
$$
となることがわかる。
すなわち、ある信号を基本ウェーブレットを使って表現すると、
信号の平均と、差分に分解することができる。
このような変換を離散ウェーブレット変換という。
もともと2点の基本単位ステップで表されていたものを4点で表すと、"解像度"が落ちるが、
ざっくりとした信号の値がわかることになるとも解釈できる。

さて、基底の変換を一般化してみる。
$m<n$として結果だけ書くと、
$$
\ket{\phi_0;2(n-m+1),m} = \ket{\phi_0;n-m+1,m} + \ket{\phi_0;n-m+1,n+1},
\ket{\psi_0;2(n-m+1),m} = \ket{\phi_0;n-m+1,m} - \ket{\phi_0;n-m+1,n+1}
$$
となる。試しに$m=1,n=2$を代入すると、上記の基底の変換が再現される。

信号分解

基本ウェーブレットを使って信号を分解することを考えてみる。
具体例として8サンプルの信号を扱う。
$$
\ket{s} = [2\ 4\ 3\ 3\ 5\ 1\ 9\ 9]^t= \sum_{n=1}^8 s_n \ket{\phi_0;1,n}
$$
としてみよう。
まず初めに2倍の幅を持つ基底$\ket{\phi_0,2,n},\ket{\psi_0;2,n}(n=1,3,5,7)$を
使って表現すると、
$$
\ket{\phi_0;1,n} = \frac{1}{2}(\ket{\phi_0,2,n}+\ket{\psi_0;2,n}),
\ket{\phi_0;1,2n} = \frac{1}{2}(\ket{\phi_0,2,n}-\ket{\psi_0;2,n})
$$
なので、

 \ket{s} = [3\ 3\ 3\ 9\ -1\ 0\ 2\ 0]^t= \sum_{n=1}^4 (s_n^\phi \ket{\phi_0;2,n} + s_n^\psi \ket{\psi_0;2,n})

となる。ここで上付きの$\phi,\psi$はそれぞれ基本単位ステップ、基本ウェーブレットを表している。
ここで基本単位ステップはまだ基本ウェーブレットに分解できるので、分解できなくなるまで計算を実行すると、

 \ket{s} = \left[\frac{9}{2}\ -\frac{3}{2}\ 0\ -3\ -1\ 0\ 2\ 0
 \right]^t
 =\frac{9}{2} \ket{\phi_0;8,1}
  -\frac{3}{2} \ket{\psi_0;8,1}
  +0 \ket{\psi_0;4,1} -3 \ket{\psi_0;4,5}
  -1 \ket{\psi_0;2,1} + 0  \ket{\psi_0;2,3}
  +2 \ket{\psi_0;2,5} + 0  \ket{\psi_0;2,7}

となる。
右辺第1項は全データの平均の数を表している。
第2項、第3項＋第4項、第5項-第8項はそれぞれ組となっており、
順番に低周波数から高周波数を表している。

プログラム例

juliaで離散ウェーブレット変換してみる。
愚直に実装した冗長なプログラムなのでより最適なものはあると思う。
久々にjuliaを使ったらいろいろ忘れている…
ver1.7.2です。

実装の流れとしては以下の通り。
まずウェーブレット変換をかける信号を用意

次にウェーブレット変換を実装する。
このとき、前節で行った基底の変換を一般化し、隣通しの和と差の半分をとり、
和の半分(平均)は配列の番号が小さい側に、差は配列の番号が大きい側に代入している。
つまり基底として基本単位ステップは前側、基本ウェーブレットは後ろ側に置いている。

さらに逆変換も用意した。これはウェーブレット変換の反対をたどればよい。

最後に、ローパスフィルタを実装して高周波成分を除去した場合に原信号がどのように変化するか確認した。

dwt.jl

using Plots

fs = 8 #サンプリング周波数
T = 1 #周期
ω = 2π/T #角周波数
dt = 1/fs #サンプリング間隔
N = 64 #データ数

xn = [dt*i for i in 0:N-1] #データ列
f(x,ω) = sin(ω*x)
y(x) = f(x,ω) 
yn = y.(xn) #信号列、今回はsin波

"""
離散ウェーブレット変換の関数
配列を並べ替えて計算を実行
"""
function dwt(arr)
    function wt(arr)
        new_arr = zeros(length(arr))
        for i in eachindex(arr)
            if i%2==1
                new_arr[div(i+1,2)] = 0.5*(arr[i]+arr[i+1])
            else
                new_arr[div(i+length(arr),2)] = 0.5*(arr[i-1]-arr[i])
            end
        end
        return new_arr
    end

    w = []
    new_arr = zeros(length(arr))
    N = length(arr)
    while N>1
        wlist = []
        new_arr = wt(arr[1:N])
        for j in div(N,2)+1:N
            push!(wlist,new_arr[j])
        end
        N = div(N,2)
        arr = new_arr
        w = vcat(wlist,w)
        if N==1
            w = vcat(arr[1],w)
        end
    end
    #println(wlist)
    return w
end

"""
逆離散ウェーブレット変換
もとに戻るように実装
"""
function inv_dwt(arr)
    function inv1(arr1)     
        N = length(arr1)
        new_arr = zeros(N)
        for n in 1:div(N,2)
            new_arr[2n-1] = arr1[n] + arr1[n+div(N,2)]
            new_arr[2n] = arr1[n] - arr1[n+div(N,2)]
        end
        return new_arr
    end

    N=length(arr)
    nlist = []
    n = N
    while 0<n
        push!(nlist,n)
        n = div(n,2)
    end
    nlist = reverse(nlist)
    pop!(nlist)
    new_arr = zeros(N)
    c = copy(arr)
    for n in nlist
        c[1:2n] = inv1(c[1:2n])
    end
    new_arr = c

    return new_arr
end

"""
ローパスフィルタ
高周波成分を無理やり0にしてカットする
"""
function lpf(arr)
    N = length(arr)
    for j in div(N,2):N
        arr[j] = 0.0
    end
    return arr
end


plot(xn,yn)
plot!(xn,inv_dwt(lpf(dwt(yn))))

グラフを描画。y1は原信号、y2はウェーブレット変換後LPFにかけて逆変換したもの。
高周波成分がなくなることでなめらかではなくなり、角ばった波形が得られた。

ウェーブレット変換したデータのプロット。
64点データがあるので、
1,1,2,4,8,16,32個の点のグループに分かれていることがわかる。
上記のLPFは32個の点からなるグループを排除したことになる。
どうやら各グループは周波数の順番に並んでおり、
またそれぞれのグループで全時間N/fs秒すべてを覆いつくしているらしい(本当？)。
しかしこの波形の意味はよくわかっていない。課題である。

まとめ

プログラムに落とし込むことで何となく離散ウェーブレット変換の操作が理解できた。
しかしながらその結果の解釈はまだ理解が及んでいない。

参考文献

1. 三谷政昭 著、やり直しのための通信数学(CQ出版社)
2. 和田成夫 著、よくわかる信号処理(森北出版)
3. https://qiita.com/kaityo256/items/70dc20658ef98d229de9