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

使用するライブラリ

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

精度の設定

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

端末

>>> import sympy
>>> sympy.mpmath.mp.prec
53
>>> sympy.mpmath.mp.prec = 256

DKA法で因数分解

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

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

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

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

端末

>>> 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法を使用します。

端末

>>> 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}$ となるように調整しておきます。

端末

>>> 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)$ を作成します。

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

端末

>>> 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

係数を取り出します。

端末

>>> 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}$ となるように調整しておきます。

端末

>>> 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)$ を作成します。

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

端末

>>> 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なので実数部のみ使用します。

端末

>>> 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つおきに符号を反転したものになります。

端末

>>> 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つおきに符号を反転したものになります。

端末

>>> 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

カスケードアルゴリズムで描画

カスケードアルゴリズムでスケーリング関数の近似値を求めることが出来ます。

ソースコード

参考文献

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