離散ウェーブレット変換の勉強の備忘録としてまとめておきます。
尚、間違っている箇所を含む可能性があります。ご了承ください。

ウェーブレット変換とは

ウェーブレット変換とは周波数解析手法の一つで、信号をマザーウェーブレットと呼ばれる基底関数の重ね合わせ（拡大縮小・平行移動含む）で表現する手法である。

マザーウェーブレットは条件を満たせばどんな関数でもよい。PyWaveletsには以下が用意されている。

['haar', 'db', 'sym', 'coif', 'bior', 'rbio', 'dmey', 'gaus', 'mexh', 'morl', 'cgau', 'shan', 'fbsp', 'cmor']

例えば以下のような関数がある(db2)。

sin波やcos波とは異なり有限の長さの波であることに注目したい。

ウェーブレット変換の種類

ウェーブレット変換は（他の周波数変換と同様に）以下の2種類が存在する。

連続ウェーブレット変換(Continuous Wavelet Transform; CWT)

連続ウェーブレット変換では、信号$f(t)$に対して以下のように定義される。

W(a,b)=\frac{1}{\sqrt{a}}\int_{-\infty}^{\infty}f(t)\psi\bigg(\frac{t-b}{a}\bigg)dt

$a,b$はそれぞれ拡大縮小、平行移動を表すパラメータであり、$\psi$はウェーブレット関数である。

離散ウェーブレット変換(Discrete Wavelet Transform; DWT)

DWTは多重解像度解析として実装されていることが多い。厳密には別のものであるが、多くのソフトウエアでは同じものとして扱われるようである。

多重解像度解析では、ウェーブレット関数$\psi$に加えてスケーリング関数$\phi$が存在し、これらの関数の重ね合わせで信号を表現する。$\psi$をマザーウェーブレットと呼ぶのに対し、$\phi$をファーザーウェーブレットやスケーリング関数と呼び、この2つの関数にはトゥースケールの関係が成立している。

ウェーブレット変換の結果として得られる出力は、詳細係数(detail coefficients) と 近似係数(approximation coefficients) であり、これらは繰り返されるフィルタの畳み込みによって得られる。(ここで用いられるフィルタの値がウェーブレット関数の種類に依存している。)したがって、離散ウェーブレット変換は言うなればフィルタリング処理である。

画像に対する2次元離散ウェーブレット変換

ウェーブレット変換は画像の圧縮に用いられる。
画像に対して2次元ウェーブレット変換を行うと以下のようになる。

LLは画像の低周波成分、LH,HLは低周波と高周波、HHは高周波成分である。
LLに対して繰り返し2次元ウェーブレット変換を行うとさらに低周波成分と高周波成分に分離されていく。

高周波成分のデータ量を圧縮することで画像全体の品質ある程度保ちながらデータ量を削減することが可能になっており、JPEG2000という規格で採用されている。

PyWaveletsを見ていく

PyWaveletsとは

Pywavelets (GitHub:https://github.com/PyWavelets/pywt) はPythonでウェーブレット変換を行うためのオープンソースライブラリである。下記ではこのOSSを見ていく。

DWTで使われる変数

主な登場人物となりそうな変数は以下に並んでいる。

pywt/_extensions/c/wavelets.h

typedef struct {
    BaseWavelet base;
    double* dec_hi_double;  /* highpass decomposition */
    double* dec_lo_double;  /* lowpass decomposition */
    double* rec_hi_double;  /* highpass reconstruction */
    double* rec_lo_double;  /* lowpass reconstruction */
    float* dec_hi_float;
    float* dec_lo_float;
    float* rec_hi_float;
    float* rec_lo_float;
    size_t dec_len;   /* length of decomposition filter */
    size_t rec_len;   /* length of reconstruction filter */

    int vanishing_moments_psi;
    int vanishing_moments_phi;

} DiscreteWavelet;

ほとんどはフィルタに関するものである。hi,loはそれぞれhighpass,lowpassフィルタを表し、decomposition(分解）用とreconstruction(再構成)用がある。

vanishing_moments(消失モーメント）はウェーブレット関数の性質を表す指標であるらしい。bior2.2などの数字部分を表しているのではないかと推測した。

フィルタ係数の値

フィルタ係数の値はwaveletの種類によって決定される。
(こちらを参照：https://wavelets.pybytes.com/ )

これらの値はPyWaveletsではヘッダーファイルに定数として記述されている。
(https://github.com/PyWavelets/pywt/blob/master/pywt/_extensions/c/wavelets_coeffs.template.h)

この定数はどこから？と思い辿ってみると、計算によって算出されていることがわかる。（理解できていないので、計算方法については論文参照）

・https://repository.kulib.kyoto-u.ac.jp/dspace/handle/2433/60181
・https://jqichina.files.wordpress.com/2012/02/ten-lectures-of-waveletsefbc88e5b08fe6b3a2e58d81e8aeb2efbc891.pdf

このヘッダファイルから値を読み出し、dec_hi_doubleなどの変数に値を代入している箇所はpywt/_extensions/c/wavelets.cである。

pywt.dwtを追う

1次元のDWTを行う場合、使う関数はpywt.dwtとなりそうである。そこで、今回はpywt.dwtを追ってイメージを掴みたい。
（画像等で2次元のDWTを行う場合はpywt.dwt2になるだろうが、基本的にはdwtが2軸で行われるだけだと思う）

実装箇所はこちら(GitHub: https://github.com/PyWavelets/pywt/blob/master/pywt/_dwt.py)

詳細はさておき、まずは関数の説明を引用します。

"""
dwt(data, wavelet, mode='symmetric', axis=-1)
Single level Discrete Wavelet Transform.
Parameters
----------
data : array_like
Input signal
wavelet : Wavelet object or name
Wavelet to use
mode : str, optional
Signal extension mode, see :ref:Modes <ref-modes>.
axis: int, optional
Axis over which to compute the DWT. If not given, the
last axis is used.
Returns
-------
(cA, cD) : tuple
Approximation and detail coefficients.
Notes
-----
Length of coefficients arrays depends on the selected mode.
For all modes except periodization:
len(cA) == len(cD) == floor((len(data) + wavelet.dec_len - 1) / 2)
For periodization mode ("per"):
len(cA) == len(cD) == ceil(len(data) / 2)
Examples
--------
>>> import pywt
>>> (cA, cD) = pywt.dwt([1, 2, 3, 4, 5, 6], 'db1')
>>> cA
array([ 2.12132034, 4.94974747, 7.77817459])
>>> cD
array([-0.70710678, -0.70710678, -0.70710678])
"""

引数：入力データやウェーブレットの種類
返り値：近似係数と近似係数（ただし、サイズは入力データの半分）

であることがわかる。

内部では、

pywt/_dwt.py

 cA, cD = dwt_single(data, wavelet, mode)

のようにしてdwt_single関数が呼びだされている。

さらにdwt_single内では、

pywt/_extensions/_dwt.pyx

retval_a = c_wt.double_dec_a(&data[0], data_size, wavelet.w,
                            <double *>cA.data, output_len, mode)
retval_d = c_wt.double_dec_d(&data[0], data_size, wavelet.w,
                            <double *>cD.data, output_len, mode)

のような関数が呼びだされている。

c_の部分は実装c言語であることを表し、double_dec_aは倍精度浮動小数点型でdecomposition(分解)し、apporoximation coefficients(近似係数)やdetail coefficients(詳細係数）を返すという意味であろう。

ここまで何も進んでいないように見えるが、ウェーブレットに関して言えば実際その通りで、ここまでは分岐処理やエラー処理が主である。

上記をまとめると
dwt -> dwt_single -> double_dec_a
の順に関数が呼ばれている。

肝心の関数double_dec_aが探しても見つからず、どこにも定義されていないと思ったのだが、CATマクロで名前が変化していた。(CATはconcatinateのことで、引数2つの値を連結するっぽい）

pywt/_extensions/c/wt.template.c

int CAT(TYPE, _dec_a)(const TYPE * const restrict input, const size_t input_len,
                      const DiscreteWavelet * const restrict wavelet,
                      TYPE * const restrict output, const size_t output_len,
                      const MODE mode){

    /* check output length */
    if(output_len != dwt_buffer_length(input_len, wavelet->dec_len, mode)){
        return -1;
    }

    return CAT(TYPE, _downsampling_convolution)(input, input_len,
                                                wavelet->CAT(dec_lo_, REAL_TYPE),
                                                wavelet->dec_len, output,
                                                2, mode);
}


/* Decomposition of input with highpass filter */

int CAT(TYPE, _dec_d)(const TYPE * const restrict input, const size_t input_len,
                      const DiscreteWavelet * const restrict wavelet,
                      TYPE * const restrict output, const size_t output_len,
                      const MODE mode){

    /* check output length */
    if(output_len != dwt_buffer_length(input_len, wavelet->dec_len, mode))
        return -1;

    return CAT(TYPE, _downsampling_convolution)(input, input_len,
                                                wavelet->CAT(dec_hi_, REAL_TYPE),
                                                wavelet->dec_len, output,
                                                2, mode);
}

近似係数を返す関数も近似係数を返す関数もTYPE_downsampling_convolution(TYPEはdouble等）という名前である。
3つ目の引数に注目すると、dec_lo_とdec_hiと異なっており、これらはローパスフィルタとハイパスフィルタを表すと思われる。

さて、TYPE_downsampling_convolutionは何をしているか見る。
下記はその一部である。

pywt/_extensions/c/convolution.template.c

int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, const size_t N,
                                         const REAL_TYPE * const restrict filter, const size_t F,
                                         TYPE * const restrict output,
                                         const size_t step, MODE mode)
{
    /* This convolution performs efficient downsampling by computing every
     * step'th element of normal convolution (currently tested only for step=1
     * and step=2).
     */

    size_t i = step - 1, o = 0;

    if(mode == MODE_PERIODIZATION)
        return CAT(TYPE, _downsampling_convolution_periodization)(input, N, filter, F, output, step, 1);

    if (mode == MODE_SMOOTH && N < 2)
        mode = MODE_CONSTANT_EDGE;

    // left boundary overhang
    for(; i < F && i < N; i+=step, ++o){
        TYPE sum = 0;
        size_t j;
        for(j = 0; j <= i; ++j)
            sum += filter[j]*input[i-j];

        switch(mode) {
        case MODE_SYMMETRIC:
            while (j < F){
                size_t k;
                for(k = 0; k < N && j < F; ++j, ++k)
                    sum += filter[j]*input[k];
                for(k = 0; k < N && j < F; ++k, ++j)
                    sum += filter[j]*input[N-1-k];
            }
            break;

sum += filter[j]*input[k]の行から、やはり最終的にはフィルタ処理が行われていることがわかる。

まとめ

この記事ではPyWaveletsのソースコードから、離散ウェーブレット変換がどのようなものなのか調べた。

PyWaveletsの実装ではハイパスフィルタとローパスフィルタを、異なるサイズの画像に対して繰り返し適用するものであることがわかった。

参考

・やり直しのための通信数学