多重解像度解析のソースを表示

'''多重解像度解析'''（たじゅうかいぞうどかいせき、{{lang-en-short|multiresolution analysis, MRA}}）とは、2倍毎の解像度の[[ウェーブレット]]を用いて[[離散ウェーブレット変換]]により解析する手法。スケーリング関数で[[基底展開]]された信号列を、半分の解像度のスケーリング関数とウェーブレット関数による基底展開の和に分解する。1989年に Stephane G. Mallat が発表した<ref>{{cite journal
|title=A theory for multiresolution signal decomposition: the wavelet representation
|author=Stephane G. Mallat
|year=1989
|journal=Pattern Analysis and Machine Intelligence, IEEE Transactions on
|volume=11
|issue=7
|pages=674-693
|doi=10.1109/34.192463
|url=http://www.cmap.polytechnique.fr/~mallat/papiers/MallatTheory89.pdf
}}</ref>。

本来は異なる物だが、[[Mathematica]]<ref>[https://reference.wolfram.com/language/ref/DiscreteWaveletTransform.html DiscreteWaveletTransform—Wolfram言語ドキュメント]</ref> や [[MATLAB]]<ref>[http://jp.mathworks.com/help/wavelet/ref/dwt.html Single-level discrete 1-D wavelet transform - MATLAB dwt - MathWorks 日本]</ref> をはじめとして、多くのソフトウェアでは多重解像度解析の事を離散ウェーブレット変換と呼んでいる。離散ウェーブレット変換の本来の定義は、[[離散ウェーブレット変換]]の項目を参照。

== 概要 ==
関数 <math>f_j(x)</math> をスケーリング関数 <math>\phi</math> で展開した上で、
: <math>f_j(x) = \sum_k c_k^{(j)} 2^{j/2} \phi(2^j x - k) = \sum_k c_k^{(j)} \phi_{j, k}(x)</math>
下記のウェーブレット関数 <math>\psi</math> への展開を用いて、
: <math>g_j(x) = \sum_k d_k^{(j)} 2^{j/2} \psi(2^j x - k) = \sum_k d_k^{(j)} \psi_{j, k}(x)</math>
関数 <math>f_j(x)</math> を異なる解像度（レベル）のウェーブレット関数に展開していく手法を多重解像度解析という。
: <math>f_j(x) = g_{j-1}(x) + g_{j-2}(x) + \cdots + g_{j-m}(x) + f_{j-m}(x)</math>
下記関数集合が[[基底関数]]として使われる。
: <math>\{ \phi_{j-m, k}(t), \psi_{j-m, k}(t), \cdots, \psi_{j-1, k}(t) \mid k \in \mathbf{Z} \}</math>
スケーリング関数とウェーブレット関数の間にはトゥースケール関係が成立する事が必要である。

== トゥースケール関係 ==
以下の関係が成立する場合、トゥースケール関係と呼ぶ。
: <math>\phi(x) = \sum_{k \in \mathbf{Z}} p_k \phi(2x - k)</math>
: <math>\psi(x) = \sum_{k \in \mathbf{Z}} q_k \phi(2x - k)</math>

<math>p_k, q_k</math> の値はウェーブレット関数ごとに異なる。例えば[[ハールウェーブレット]]の場合は、<math>p_0 = p_1 = q_0 = 1,\ q_1 = -1</math> 。

トゥースケール関係が成立していると、下記の式が成立する。
: <math>\phi(2x - l) = \frac{1}{2} \sum_{k \in \mathbf{Z}} ( g_{2k - l} \phi(x - k) + h_{2k - l} \psi(x- k) ),\quad l \in \mathbf{Z}</math>

<math>\{g_k\}, \{h_k\}</math> を分解数列と呼ぶ。例えば、直交ウェーブレットの場合、<math>g_k = \overline{p}_{-k},\ h_k = \overline{q}_{-k}</math> 。

分解数列が分かると、<math>f_j(x) = g_{j-1}(x) + f_{j-1}(x)</math> という変形が可能になり、これを再帰的に繰り返すと多重解像度解析になる。

== 2次元 ==
2次元の場合は、まず、関数 <math>f_j(x, y)</math> をスケーリング関数 <math>\phi</math> で基底展開する。
: <math>f_j(x, y) = \sum_{k_1} \sum_{k_2} c_{k_1, k_2}^{(j)} 2^j \phi(2^j x - k_1) \phi(2^j y - k_2) = \sum_{k_1} \sum_{k_2} c_{k_1, k_2}^{(j)} \phi_{j, k_1}(x) \phi_{j, k2}(y)</math>
<math>\phi_{j, k_1}(x) \phi_{j, k_2}(y)</math> に対して、
: <math>\phi_{j, k_1}(x)</math> は <math>\phi_{j - 1, k_1}(x)</math> と <math>\psi_{j - 1, k_1}(x)</math> に分解し、
: <math>\phi_{j, k_2}(y)</math> は <math>\phi_{j - 1, k_2}(y)</math> と <math>\psi_{j - 1, k_2}(y)</math> に分解し、
合わせて、
: <math>\{ \phi_{j - 1, k_1}(x) \phi_{j - 1, k_2}(y),\ \phi_{j - 1, k_1}(x) \psi_{j - 1, k_2}(y),\ \psi_{j - 1, k_1}(x) \phi_{j - 1, k_2}(y),\ \psi_{j - 1, k_1}(x) \psi_{j - 1, k_2}(y) \}</math> 
の4つに分解する。そして、<math>\phi_{j - 1, k_1}(x) \phi_{j - 1, k_2}(y)</math> を同じように再帰的に分解していく。

結果として、m 回繰り返すと、下記の関数集合が[[基底関数]]となる。
: <math>\{ \phi_{j-m, k_1}(x) \phi_{j-m, k_2}(y) \mid k_1 \in \mathbf{Z}, k_2 \in \mathbf{Z} \}\ \cup\ \{ \phi_{j_1, k_1}(x) \psi_{j_1, k_2}(y),\ \psi_{j_1, k_1}(x) \phi_{j_1, k_2}(y),\ \psi_{j_1, k_1}(x) \psi_{j_1, k_2}(y) \mid k_1 \in \mathbf{Z}, k_2 \in \mathbf{Z}, j - 1 \le j_1 \le j - m \}</math>

3次元以上も同じ。

== フィルタバンクによる表現 ==
=== 1レベルの変換 ===
信号<math>x</math>の離散ウェーブレット変換は、一組のフィルターを通すことによって計算される。ここではリフティングスキームではなく1989年に発表された手法を説明する。

下記の関係性を満たす <math>y_{\mathrm{low}}[k]</math> と <math>y_{\mathrm{high}}[k]</math> を <math>x[k]</math> から計算するのが離散ウェーブレット変換である。j は入力依存の何らかの整数。<math>\phi(t)</math> はスケーリング関数、<math>\psi(t)</math> はウェーブレット関数。
: <math>\sum_k x[k] \phi(2^j t - k) = \sum_k y_{\mathrm{low}}[k] \phi(2^{j-1} t - k) + \sum_k y_{\mathrm{high}}[k] \psi(2^{j-1} t - k)</math>

最初に入力信号列から <math>x[k]</math> を計算しないといけないが、[[ハールウェーブレット]]や Bior2.2 双直交ウェーブレットの場合は、入力信号の解像度と同じ解像度のスケーリング関数を使えば、<math>x[k]</math> = 入力信号列になり、特に計算は不要となる。

信号は、[[インパルス応答]]が<math>g</math>である[[ローパスフィルタ]]と、同じく<math>h</math>である[[ハイパスフィルタ]]を利用して分解される。得られた出力は、ハイパスフィルタから得たものを詳細係数(detail coefficients)、ローパスフィルタから得たものを近似係数(approximation coefficients)とよぶ。これら2つのフィルタは互いに密接な関係があり、[[直交ミラーフィルタ]]として知られたものである。

:<math>y_{\mathrm{low}} [n] = \sum\limits_{k =  - \infty }^\infty  {x[k] g[2n - k]} </math>
:<math>y_{\mathrm{high}} [n] = \sum\limits_{k =  - \infty }^\infty  {x[k] h[2n - k]} </math>

次に、半分にダウンサンプリングを行う。この分解によって、時間解像度は元の信号の半分になり、周波数解像度は2倍となる。これは、ハイゼンベルクの[[不確定性原理]]を満たしている。

[[Image:Wavelets_-_DWT.png|frame|none|フィルタ解析のブロックダイアグラム]]

ダウンサンプリングのオペレータとして<math>\downarrow</math>を使う。

:<math>(y \downarrow k)[n] = y[k\cdot n] </math>

以上の式をまとめると以下のようになる。

:<math>y_{\mathrm{low}} = (x*g)\downarrow 2 </math>
:<math>y_{\mathrm{high}} = (x*h)\downarrow 2 </math>

しかしながら、<math>x*g</math>の計算の後にダウンサンプリングがあるため、計算に無駄がある。リフティングスキームは、この点を改善している。

=== 多重解像度解析 ===
この分解は、近似係数を再び分解することによって繰り返され、これを多重解像度解析と呼ぶ。これは、異なる時間周波数の局在性を持った部分空間を枝とする二分木によって表すことが出来る。これは[[フィルタバンク]]として知られているものである。

[[Image:Wavelets_-_Filter_Bank.png|frame|none|3段階のフィルタバンク]] 

各々のレベルにおいて、低周波と高周波に分解される。<math>2^n</math> の長さの入力信号は、[[ハールウェーブレット]]の場合は <math>n</math> のレベルに分解される。他のウェーブレットの場合も、おおよそ <math>n</math> になるが、多少ずれる。

[[ハールウェーブレット]]の場合での、32サンプルの信号を例にとる。周波数レンジは0から <math>f_n</math> であり、3レベルまで分解するとすると、以下の表のようになる。

{| class="wikitable"
! レベル
! 周波数
! 係数の数
! 基底関数
|-
| rowspan="2" | 3
| <math>0</math> 〜 <math>{{f_n}}/8</math>
| 4
| <math>\phi_{j-3,k}</math>
|-
| <math>{{f_n}}/8</math> 〜 <math>{{f_n}}/4</math>
| 4
| <math>\psi_{j-3,k}</math>
|-
| 2
| <math>{{f_n}}/4</math> 〜 <math>{{f_n}}/2</math>
| 8
| <math>\psi_{j-2,k}</math>
|-
| 1
| <math>{{f_n}}/2</math> 〜 <math>f_n</math>
| 16
| <math>\psi_{j-1,k}</math>
|}

[[Image:Wavelets_-_DWT Freq.png|frame|none|多重解像度解析の周波数領域]]

== リフティングスキーム ==
1996年に Wim Sweldens が新しい離散ウェーブレット変換の計算方法である'''{{仮リンク|リフティングスキーム|en|Lifting scheme}}'''（{{lang-en-short|lifting scheme}}）を発表した<ref>{{cite journal
|title=The Lifting Scheme: A Custom-Design Construction of Biorthogonal Wavelets
|author=Wim Sweldens
|journal=Applied and Computational Harmonic Analysis
|volume=3
|issue=2
|pages=186–200
|year=1996
|month=April
|doi=10.1006/acha.1996.0015
}}</ref>。信号を偶数番目の要素と奇数番目の要素に分け、偶数番目の要素で奇数番目の要素を予測し、予測とのズレを高周波成分とし、残差を低周波成分とする手法である。論文の中で、自由に双直交ウェーブレット関数を定義できること、計算が高速化する事を長所に挙げている。[[ハールウェーブレット]]と双直交ウェーブレットを扱える。

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

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

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

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

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

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

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

== プログラム例 ==
=== Python（ハールウェーブレット） ===
[[ハールウェーブレット]]の離散ウェーブレット変換は下記の{{仮リンク|リフティングスキーム|en|Lifting scheme}}の数式にて行える<ref name="lifting_matlab">[http://jp.mathworks.com/help/wavelet/ug/lifting-method-for-constructing-wavelets.html Lifting Method for Constructing Wavelets - MATLAB & Simulink - MathWorks 日本]</ref>。入力を x として、<math>x_e(z)</math> は入力 x の先頭を1番目として偶数番目の要素、つまり、<code>x[1::2]</code>、<math>x_o(z)</math> は奇数番目の要素、つまり、<code>x[0::2]</code>。c が低周波成分（基底関数がスケーリング関数）、d が高周波成分（基底関数がウェーブレット関数）。
: <math>\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}</math>

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

ハールウェーブレットの逆離散ウェーブレット変換のソースコード。
<syntaxhighlight lang="python">
x[1::2] = c - d / 2
x[0::2] = d + x[1::2]
</syntaxhighlight>

=== Python（Bior2.2双直交ウェーブレット） ===
Bior2.2双直交ウェーブレット（2階Bスプライン、線形スプライン）の離散ウェーブレット変換は下記のリフティングスキームの数式にて行える<ref name="lifting_matlab"/>。<math>z</math> は一つ次の要素、<math>z^{-1}</math> は一つ前の要素を表す。
: <math>\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}</math>

Bior2.2の離散ウェーブレット変換のソースコード。配列の境界で足りない分は対称パディングで埋めている。
<syntaxhighlight lang="python">
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:]
</syntaxhighlight>

Bior2.2の逆離散ウェーブレット変換のソースコード。
<syntaxhighlight lang="python">
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
</syntaxhighlight>

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

<syntaxhighlight lang="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);
    }
}
</syntaxhighlight>

下記コードは上記を逆ウェーブレット変換する。
<syntaxhighlight lang="java">
/**
 * 注意：このメソッドは 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);
    }
}
</syntaxhighlight>

=== C言語 ===
C言語による、[[:en:Cohen-Daubechies-Feauveau_wavelet|CDF]] 9/7 ウェーブレット変換（[[JPEG-2000]]で利用）のリフティングを用いた高速な実装は、[https://code.google.com/p/axonlib/source/browse/trunk/extern/dwt97.c?spec=svn19&r=19 dwt97.c] を参照。

== 関連項目 ==
* [[ウェーブレット]]
* [[ウェーブレット変換]]
* [[離散ウェーブレット変換]]

== 参照 ==
{{reflist}}

{{DEFAULTSORT:たしゆうかいそうとかいせき}}
[[Category:ウェーブレット]]
[[Category:数学に関する記事]]