[PYTHON] How to find the scaling factor of a biorthogonal wavelet

There are various types of biorthogonal wavelets, but this time we will find the scaling factor of ** 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/

Library to use

We will use sympy to calculate the scaling factor in Python's interactive mode.

Precision setting

sympy uses mpmath for arbitrary precision floating point arithmetic. The default precision (the number of bits in the mantissa) is 53, which is the same as the double precision floating point number, so set it to about 256 as appropriate.

Terminal

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

Factorization by DKA method

Set $ N = 4 $ to find the solution of the following equation.

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

The specific description is as follows.

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

Terminal

>>> 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')]

Since it is a cubic equation, it is possible to find an exact solution with Cardano's formula, but since it is sufficient to obtain the scaling coefficient with a floating point number, DKA method Use E2% 80% 93Kerner_method).

Terminal

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

Since one real solution and two conjugate complex solutions are obtained, $ q (y) $ can be decomposed into the product of the linear solution obtained from the real solution and the quadratic expression obtained from the conjugate complex solution.

Scaling coefficient

Create $ h_0 (z) $ using a real solution. Set $ y =-(1/4) z + (1/2)-(1/4) z ^ {-1} $ and adjust so that $ h_0 (1) = \ sqrt {2} $ I will leave it.

Terminal

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

Apply vanishing moments to create $ h (z) $.

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

Terminal

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

Take out the coefficient.

Terminal

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

Dual scaling factor

Create $ \ tilde {h_0} (z) $ using a pair of conjugate complex solutions. Put $ y =-(1/4) z + (1/2)-(1/4) z ^ {-1} $ and then $ \ tilde {h_0} (1) = \ sqrt {2} $ I will adjust it so that.

Terminal

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

Apply vanising moment to create $ \ tilde {h} (z) $.

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

Terminal

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

Take out the coefficient. Since the imaginary part is practically 0, only the real part is used.

Terminal

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

Wavelet coefficient

The wavelet coefficient $ g_n $ is the inverted sign of every other dual scaling coefficient $ \ tilde {h} _n $.

Terminal

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

Dual wavelet coefficient

The dual wavelet coefficient $ \ tilde {g} _n $ is the inverted sign of every other scaling coefficient $ h_n $.

Terminal

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

Draw with cascade algorithm

You can find the approximate value of the scaling function with the Cascade algorithm.

• Source code

References

• I. Dobsey, Wavelet 10 Lectures (Springer Fairlark Tokyo)