双直交ウェーブレットのスケーリング係数の求め方

双直交ウェーブレットも言っても色々ありますが、今回は **Cohen-Daubechies-Feauveau wavelet** のスケーリング係数を求めます。

+ http://reference.wolfram.com/language/ref/CDFWavelet.html
+ http://en.wikipedia.org/wiki/Cohen-Daubechies-Feauveau_wavelet
+ http://wavelets.pybytes.com/wavelet/bior4.4/

# 使用するライブラリ
sympyを使用してPythonの対話モードでスケーリング係数を算出していきます。

# 精度の設定
sympyは任意精度浮動小数点数演算に[mpmath](http://mpmath.org)を使用しています。デフォルトの精度（仮数部のビット数）は倍精度浮動小数点数と同じ53になっていますので適当に256くらいに設定しておきます。

```py:端末
>>> import sympy
>>> sympy.mpmath.mp.prec
53
>>> sympy.mpmath.mp.prec = 256
```

# DKA法で因数分解

$N=4$と置いて次の方程式の解を求めます。

```math
q(y) =
\sum_{k=0}^{N-1}
\begin{pmatrix}
N-1+k \\
k
\end{pmatrix}
y^k = 0
```

具体的に書くと次のようになります。

```math
q(y) = 20y^3 + 10y^2 + 4y + 1 = 0
```

```py:端末
>>> qy = [sympy.mpmath.binomial(4-1+k,k) for k in [3,2,1,0]]
>>> qy
[mpf('20.0'), mpf('10.0'), mpf('4.0'), mpf('1.0')]
```

3次方程式なのでカルダノの公式で厳密な解を求めることも可能ですが、スケーリング係数は浮動小数点数で得られれば十分なので[DKA法](http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method)を使用します。

```py:端末
>>> y = sympy.mpmath.polyroots(qy)
>>> y
[mpf('-0.3423840948583691316993036540027816871936619136844427977504078911201017562570965'), mpc(real='-0.07880795257081543415034817299860915640316904315777860112479605443994912187145069', imag='0.3739306454336101226089715992265723406258901629829692932753904369966748441011043'), mpc(real='-0.07880795257081543415034817299860915640316904315777860112479605443994912187145069', imag='-0.3739306454336101226089715992265723406258901629829692932753904369966748441011043')]
```

１つの実数解と２つの共役な複素数解が求まるので $q(y)$ は実数解をから得られる1次式と共役な複素数解から得られる2次式の積に分解出来ます。

# スケーリング係数
実数解を使用して $h_0(z)$ を作成します。$y=-(1/4)z+(1/2)-(1/4)z^{-1}$ と置いて更に $h_0(1)=\sqrt{2}$ となるように調整しておきます。

```py:端末
>>> h0z = sympy.sympify('-sqrt(2.0)*(y-y0)/y0') \
...            .subs({'y':'-1/4*z+1/2-1/4/z', 'y0':y[0]})
>>> h0z
-1.032622122063015088779015687939062566223023943573534440451244148644675318677*z + 3.479457806499125323032653234616953582887408360779881380902488297289350637353 - 1.032622122063015088779015687939062566223023943573534440451244148644675318677/z
```

vanising momentを適用して $h(z)$ を作成します。

```math
h(z) = z^{-2} \left( \frac{1+z}{2} \right)^4 h_0(z)
```

```py:端末
>>> hz = (sympy.sympify('z**(-2)*((z+1)/2)**4')*h0z).expand()
>>> hz
-0.06453888262893844304868848049619141038893899647334590252820275929029220741729*z**3 - 0.04068941760955843950521309482120604262529296334464102380640551858058441483457*z**2 + 0.4180922732222122294173439451808985229992791148815490275282027592902922074173*z + 0.7884856164056644517477371190118263104712661635056882976128110371611688296691 + 0.4180922732222122294173439451808985229992791148815490275282027592902922074173/z - 0.04068941760955843950521309482120604262529296334464102380640551858058441483457/z**2 - 0.06453888262893844304868848049619141038893899647334590252820275929029220741729/z**3
```

係数を取り出します。

```py:端末
>>> scaling_coeff = [hz.coeff('z',k) for k in [-3,-2,-1,0,1,2,3]]
>>> scaling_coeff
[-0.06453888262893844304868848049619141038893899647334590252820275929029220741729, -0.04068941760955843950521309482120604262529296334464102380640551858058441483457, 0.4180922732222122294173439451808985229992791148815490275282027592902922074173, 0.7884856164056644517477371190118263104712661635056882976128110371611688296691, 0.4180922732222122294173439451808985229992791148815490275282027592902922074173, -0.04068941760955843950521309482120604262529296334464102380640551858058441483457, -0.06453888262893844304868848049619141038893899647334590252820275929029220741729]
```

|n    |h[n]|
|:---:|---:|
|0    |0.7884856164056644517477371190118263104713|
|-1,+1|0.4180922732222122294173439451808985229993|
|-2,+2|-0.04068941760955843950521309482120604262529|
|-3,+3|-0.06453888262893844304868848049619141038894|
|-4,+4|0.0|

# 双対スケーリング係数

共役な複素数解のペアを使用して $\tilde{h_0}(z)$ を作成します。$y=-(1/4)z+(1/2)-(1/4)z^{-1}$ と置いて更に $\tilde{h_0}(1)=\sqrt{2}$ となるように調整しておきます。

```py:端末
>>> d_h0z = sympy.sympify('sqrt(2.0)*(y-y1)/y1*(y-y2)/y2') \
...              .subs({'y':'-1/4*z+1/2-1/4/z', 'y1':y[1], 'y2':y[2]})
>>> d_h0z
(-0.353553390593274*z + 0.818558056535050389640616746751093460160688904372839526463833807943026504668 - 0.5288177901591365145695624730424192033800671915517564410974157359390278433361*I - 0.353553390593274/z)*(-z/4 + 0.5788079525708154341503481729986091564031690431577786011247960544399491218714 + 0.3739306454336101226089715992265723406258901629829692932753904369966748441011*I - 1/(4*z))/((-0.07880795257081543415034817299860915640316904315777860112479605443994912187145 - 0.3739306454336101226089715992265723406258901629829692932753904369966748441011*I)*(-0.07880795257081543415034817299860915640316904315777860112479605443994912187145 + 0.3739306454336101226089715992265723406258901629829692932753904369966748441011*I))
```

vanising momentを適用して $\tilde{h}(z)$ を作成します。

```math
\tilde{h}(z) = z^{-2} \left( \frac{1+z}{2} \right)^4 \tilde{h_0}(z)
```

```py:端末
>>> d_hz = (sympy.sympify('z**(-2)*((z+1)/2)**4')*d_h0z).expand()
>>> d_hz
0.03782845550699546397896088239109810588598928919558903883377596750890210864563*z**4 - 0.02384946501938000354347538567498536776364603312870487872179724070970779258271*z**3 - 0.1106244044184234045317958917055274640064847218822291494893245217456879810963*z**2 + 1.848054221112700632088327353557561911770062193154051589356698864673908481848e-78*I*z**2 + 0.3774028556126537624098752472272712265473030007883118118671540106333915390923*z - 0.00000000000000001776744107070853984299687632133069137868739252186023273178181181386681671503*I*z + 0.852698679009403501358319120148908609110388988411630399468592427310915899508 + 7.392216884450802528353309414230247647080248772616206357426795458695633927391e-78*I + 0.3774028556126537624098752472272712265473030007883118118671540106333915390923/z - 0.00000000000000001776744107070853984299687632133069137868739252186023273178180811775837448963*I/z - 0.1106244044184234045317958917055274640064847218822291494893245217456879810963/z**2 + 1.848054221112700632088327353557561911770062193154051589356698864673908481848e-78*I/z**2 - 0.02384946501938000354347538567498536776364603312870487872179724070970779258271/z**3 + 0.03782845550699546397896088239109810588598928919558903883377596750890210864563/z**4
```

係数を取り出します。虚数部は実質0なので実数部のみ使用します。

```py:端末
>>> d_scaling_coeff = [sympy.re(d_hz.coeff('z',k)) for k in [-4,-3,-2,-1,0,1,2,3,4]]
>>> d_scaling_coeff
[0.03782845550699546397896088239109810588598928919558903883377596750890210864563, -0.02384946501938000354347538567498536776364603312870487872179724070970779258271, -0.1106244044184234045317958917055274640064847218822291494893245217456879810963, 0.3774028556126537624098752472272712265473030007883118118671540106333915390923, 0.8526986790094035013583191201489086091103889884116303994685924273109158995080, 0.3774028556126537624098752472272712265473030007883118118671540106333915390923, -0.1106244044184234045317958917055274640064847218822291494893245217456879810963, -0.02384946501938000354347538567498536776364603312870487872179724070970779258271, 0.03782845550699546397896088239109810588598928919558903883377596750890210864563]
```

|n    |h~[n]|
|:---:|----:|
|0    |0.8526986790094035013583191201489086091104|
|-1,+1|0.3774028556126537624098752472272712265473|
|-2,+2|-0.1106244044184234045317958917055274640065|
|-3,+3|-0.02384946501938000354347538567498536776365|
|-4,+4|0.03782845550699546397896088239109810588599|

# ウェーブレット係数

ウェーブレット係数 $g_n$ は双対スケーリング係数 $\tilde{h}_n$ を1つおきに符号を反転したものになります。

```py:端末
>>> wavelet_coeff = [s*(-1)**k for k,s in zip([-4,-3,-2,-1,0,1,2,3,4], d_scaling_coeff)]
>>> wavelet_coeff
[0.03782845550699546397896088239109810588598928919558903883377596750890210864563, 0.02384946501938000354347538567498536776364603312870487872179724070970779258271, -0.1106244044184234045317958917055274640064847218822291494893245217456879810963, -0.3774028556126537624098752472272712265473030007883118118671540106333915390923, 0.8526986790094035013583191201489086091103889884116303994685924273109158995080, -0.3774028556126537624098752472272712265473030007883118118671540106333915390923, -0.1106244044184234045317958917055274640064847218822291494893245217456879810963, 0.02384946501938000354347538567498536776364603312870487872179724070970779258271, 0.03782845550699546397896088239109810588598928919558903883377596750890210864563]
```

|n    | g[n]|
|:---:|----:|
|0    |0.8526986790094035013583191201489086091104|
|-1,+1|-0.3774028556126537624098752472272712265473|
|-2,+2|-0.1106244044184234045317958917055274640065|
|-3,+3|0.02384946501938000354347538567498536776365|
|-4,+4|0.03782845550699546397896088239109810588599|

# 双対ウェーブレット係数

双対ウェーブレット係数 $\tilde{g}_n$ はスケーリング係数 $h_n$ を1つおきに符号を反転したものになります。

```py:端末
>>> d_wavelet_coeff = [s*(-1)**k for k,s in zip([-3,-2,-1,0,1,2,3], scaling_coeff)]
>>> d_wavelet_coeff
[0.06453888262893844304868848049619141038893899647334590252820275929029220741729, -0.04068941760955843950521309482120604262529296334464102380640551858058441483457, -0.4180922732222122294173439451808985229992791148815490275282027592902922074173, 0.7884856164056644517477371190118263104712661635056882976128110371611688296691, -0.4180922732222122294173439451808985229992791148815490275282027592902922074173, -0.04068941760955843950521309482120604262529296334464102380640551858058441483457, 0.06453888262893844304868848049619141038893899647334590252820275929029220741729]
```

|n    |g~[n]|
|:---:|----:|
|0    |0.7884856164056644517477371190118263104713|
|-1,+1|-0.4180922732222122294173439451808985229993|
|-2,+2|-0.04068941760955843950521309482120604262529|
|-3,+3|0.06453888262893844304868848049619141038894|
|-4,+4|0.0|

# カスケードアルゴリズムで描画
カスケードアルゴリズムでスケーリング関数の近似値を求めることが出来ます。
![cdf_9_7_scaling.png](https://qiita-image-store.s3.amazonaws.com/0/42791/35cbdda5-b79d-7028-2576-21029e72133c.png)
![cdf_9_7_wavelet.png](https://qiita-image-store.s3.amazonaws.com/0/42791/bace67c1-beb3-6a53-4878-9413088fbbbe.png)

+ [ソースコード](https://github.com/fukuroder/cdf_9_7_wavelet_coefficients)

# 参考文献

+  I.ドブシー、ウェーブレット10講(シュプリンガー・フェアラーク東京)