多重解像度解析 リフティングスキーム

ウィキペディア

多重解像度解析

出典: フリー百科事典『ウィキペディア（Wikipedia）』 (2021/02/16 06:10 UTC 版)

リフティングスキーム

1996年に Wim Sweldens が新しい離散ウェーブレット変換の計算方法であるリフティングスキーム（英: lifting scheme）を発表した。信号を偶数番目の要素と奇数番目の要素に分け、偶数番目の要素で奇数番目の要素を予測し、予測とのズレを高周波成分とし、残差を低周波成分とする手法である。論文の中で、自由に双直交ウェーブレット関数を定義できること、計算が高速化する事を長所に挙げている。ハールウェーブレットと双直交ウェーブレットを扱える。

定義

閉部分空間の系列 $\left\{V_{j}:j\in \mathbf {Z} \right\}\subset L^{2}(\mathbf {R} )$ が次の条件1〜5を満たすとき、$\left\{V_{j}\right\}_{j\in \mathbf {Z} }$ は多重解像度解析を成すという。$A^{c}$は$A$ の閉包を表す。

1. $V_{j}\subset V_{j+1},\ j\in \mathbf {Z}$
2. $\cap _{j\in \mathbf {Z} }V_{j}=\{0\},\ (\cup _{j\in \mathbf {Z} }V_{j})^{c}=L^{2}(\mathbf {R} )$
3. $f(x)\in V_{j}\Longleftrightarrow f(2x)\in V_{j+1}$
4. $f(x)\in V_{0}$ ならば $f(x-k)\in V_{0},\ \forall k\in \mathbf {Z}$
5. $\phi (x)\in V_{0}$ が存在して $\left\{\phi (x-k):k\in \mathbf {Z} \right\}$ が $V_{0}$ においてRiesz基底を成す。すなわち、任意の$f\in V_{0}$ に対して数列 $\left\{a_{k}:k\in \mathbf {Z} \right\}\in l^{2}$ がただ一つ存在して、$f=\sum _{k}a_{k}\phi \left(x-k\right)$ と表される。さらに定数 $C_{2}\geq C_{1}>0$ が存在して、任意の $\left\{a_{k}:k\in \mathbf {Z} \right\}\in l^{2}$ に対して不等式 $C_{1}\|a\|_{2}\leq \|\sum _{j}a_{j}\phi (x-j)\|\leq C_{2}\|a\|_{2}$ が成立する。

最後の条件はより厳しい条件として、

5'. $\phi (x)\in V_{0}$ が存在して、$\left\{\phi (x-k):k\in \mathbf {Z} \right\}$ が $V_{0}$ の正規直交基底となる。

が用いられる事も多い。ここではこの条件5'を用いる。

条件1は空間が入れ子になっていることを意味する。条件2は $P_{j}$ を空間 $V_{j}$ への直交射影作用素とすると、全ての $f\in L^{2}(\mathbf {R} )$ に対して $\lim _{j\rightarrow \infty }P_{j}f=f$ であることを意味する。条件3は条件1の全ての空間がスケール変換で得られる事を意味する。条件4は平行移動に対して空間が普遍であることを意味する。

閉部分空間の集まりが条件1〜5'を満たすときには、いつでも正規直交ウェーブレット $\psi _{j,k}(x)=2^{j/2}\psi (2^{j}x-k)$ による $L^{2}(\mathbf {R} )$ の基底 $\left\{\psi _{j,k}:j,k\in \mathbf {Z} \right\}$ が存在し、$P_{j+1}f=P_{j}f+\sum _{k\in \mathbf {Z} }\left\langle f,\psi _{j,k}\right\rangle \psi _{j,k}$ が成立する。

プログラム例

Python（ハールウェーブレット）

ハールウェーブレットの離散ウェーブレット変換は下記のリフティングスキームの数式にて行える。入力を x として、$x_{e}(z)$ は入力 x の先頭を1番目として偶数番目の要素、つまり、x[1::2]、$x_{o}(z)$ は奇数番目の要素、つまり、x[0::2]。c が低周波成分（基底関数がスケーリング関数）、d が高周波成分（基底関数がウェーブレット関数）。

${\begin{pmatrix}c\\d\end{pmatrix}}={\begin{pmatrix}1&1/2\\0&1\end{pmatrix}}{\begin{pmatrix}1&0\\-1&1\end{pmatrix}}{\begin{pmatrix}x_{e}(z)\\x_{o}(z)\end{pmatrix}}$

ハールウェーブレットの離散ウェーブレット変換のソースコードは下記のようになる。入力は x で NumPy の配列で与える。多重解像度解析をしたい場合は、x = c して長さが1になるまで繰り返す。$\phi _{n,k}(t)={\sqrt {2^{n}}}\phi (2^{n}t-k)$ の形の基底関数の係数にしたい場合は、c *= sqrt(2) と d /= sqrt(2) をする。

d = x[0::2] - x[1::2]
c = x[1::2] + d / 2

ハールウェーブレットの逆離散ウェーブレット変換のソースコード。

x[1::2] = c - d / 2
x[0::2] = d + x[1::2]

Python（Bior2.2双直交ウェーブレット）

Bior2.2双直交ウェーブレット（2階Bスプライン、線形スプライン）の離散ウェーブレット変換は下記のリフティングスキームの数式にて行える。$z$ は一つ次の要素、$z^{-1}$ は一つ前の要素を表す。

${\begin{pmatrix}c\\d\end{pmatrix}}={\begin{pmatrix}1&(z^{-1}+1)/4\\0&1\end{pmatrix}}{\begin{pmatrix}1&0\\-(1+z)/2&1\end{pmatrix}}{\begin{pmatrix}x_{e}(z)\\x_{o}(z)\end{pmatrix}}$

Bior2.2の離散ウェーブレット変換のソースコード。配列の境界で足りない分は対称パディングで埋めている。

y = np.pad(x, (4, 2), "symmetric")
d = y[2::2] - (y[1:-2:2] + y[3::2]) / 2
c = y[3:-2:2] + (d[:-1] + d[1:]) / 4
d = d[1:]

Bior2.2の逆離散ウェーブレット変換のソースコード。

x[1::2] = c[1:] - (d[:-1] + d[1:]) / 4
x[0] = 2 * d[0] + x[1]
x[2::2] = d[1:-1] + (x[1:-2:2] + x[3::2]) / 2

Java

もっとも単純な例である。ハールウェーブレットによる多重解像度解析をJavaで記述した。

/**
 * 注意：このメソッドは input 配列を破壊する。
 * input.length は2以上の2のべき乗でないといけない。
 */
public static int[] calc(int[] input) {
    int[] output = new int[input.length];

    // length は 2^n で n は1つずつ減っていく
    for (int length = input.length >> 1; ; length >>= 1) {
        for (int i = 0; i < length; i++) {
            int a = input[i * 2];
            int b = input[i * 2 + 1];
            output[i] = a + b;
            output[i + length] = a - b;
        }

        if (length == 1)
            return output;

        // 次の反復のために配列を入れ替える
        System.arraycopy(output, 0, input, 0, length << 1);
    }
}

下記コードは上記を逆ウェーブレット変換する。

/**
 * 注意：このメソッドは output 配列を破壊する。
 * output.length は2以上の2のべき乗でないといけない。
 */
public static int[] inverse(int[] output) {
    int[] input = new int[output.length];

    for (int length = 1; ; length *= 2) {
        for (int i = 0; i < length; i++) {
            int a = output[i];
            int b = output[i + length];
            input[i * 2] = (a + b) / 2;
            input[i * 2 + 1] = (a - b) / 2;
        }

        if (length == output.length >> 1)
            return input;

        // 次の反復のために配列を入れ替える
        System.arraycopy(input, 0, output, 0, length << 1);
    }
}

C言語

C言語による、CDF 9/7 ウェーブレット変換（JPEG-2000で利用）のリフティングを用いた高速な実装は、dwt97.c を参照。