There are various types of biorthogonal wavelets, but this time we will find the scaling factor of ** Cohen-Daubechies-Feauveau wavelet **.
We will use sympy to calculate the scaling factor in Python's interactive mode.
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
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.
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]
n | h[n] |
---|---|
0 | 0.7884856164056644517477371190118263104713 |
-1,+1 | 0.4180922732222122294173439451808985229993 |
-2,+2 | -0.04068941760955843950521309482120604262529 |
-3,+3 | -0.06453888262893844304868848049619141038894 |
-4,+4 | 0.0 |
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]
n | h~[n] |
---|---|
0 | 0.8526986790094035013583191201489086091104 |
-1,+1 | 0.3774028556126537624098752472272712265473 |
-2,+2 | -0.1106244044184234045317958917055274640065 |
-3,+3 | -0.02384946501938000354347538567498536776365 |
-4,+4 | 0.03782845550699546397896088239109810588599 |
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]
n | g[n] |
---|---|
0 | 0.8526986790094035013583191201489086091104 |
-1,+1 | -0.3774028556126537624098752472272712265473 |
-2,+2 | -0.1106244044184234045317958917055274640065 |
-3,+3 | 0.02384946501938000354347538567498536776365 |
-4,+4 | 0.03782845550699546397896088239109810588599 |
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]
n | g~[n] |
---|---|
0 | 0.7884856164056644517477371190118263104713 |
-1,+1 | -0.4180922732222122294173439451808985229993 |
-2,+2 | -0.04068941760955843950521309482120604262529 |
-3,+3 | 0.06453888262893844304868848049619141038894 |
-4,+4 | 0.0 |
You can find the approximate value of the scaling function with the Cascade algorithm.
Recommended Posts