[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 **.

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.


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

Factorization by DKA method

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

q(y) =
N-1+k \\
y^k = 0

The specific description is as follows.

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

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


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


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


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


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

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.


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


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


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

Wavelet coefficient

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


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

Dual wavelet coefficient

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


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

Draw with cascade algorithm

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


Recommended Posts

How to find the scaling factor of a biorthogonal wavelet
How to find the memory address of a Pandas dataframe value
How to calculate the volatility of a brand
How to find the area of the Voronoi diagram
[Circuit x Python] How to find the transfer function of a circuit using Lcapy
[Ubuntu] How to delete the entire contents of a directory
How to find the optimal number of clusters in k-means
How to connect the contents of a list into a string
How to find the average amount of information (entropy) of the original probability distribution from a sample
How to determine the existence of a selenium element in Python
How to check the memory size of a variable in Python
How to check the memory size of a dictionary in Python
How to output the output result of the Linux man command to a file
How to get the vertex coordinates of a feature in ArcPy
[NNabla] How to remove the middle tier of a pre-built network
[Python] A simple function to find the center coordinates of a circle
How to check the version of Django
[Introduction to Python] How to sort the contents of a list efficiently with list sort
[NNabla] How to add a quantization layer to the middle layer of a trained model
How to put a line number at the beginning of a CSV file
How to create a wrapper that preserves the signature of the function to wrap
How to play a video while watching the number of frames (Mac)
A simple example of how to use ArgumentParser
Combinatorial optimization to find the hand of "Millijan"
How to find the correlation for categorical variables
How to pass the execution result of a shell command in a list in Python
How to mention a user group in slack notification, how to check the id of the user group
[NNabla] How to get the output (variable) of the middle layer of a pre-built network
How to count the number of elements in Django and output to a template
How to access the contents of a Linux disk on a Mac (but read-only)
[python] How to sort by the Nth Mth element of a multidimensional array
[Numpy, scipy] How to calculate the square root of a semi-fixed definite matrix
How to find the coefficient of the trendline that passes through the vertices in Python
How to make a Raspberry Pi that speaks the tweets of the specified user
How to get a list of files in the same directory with python
[Introduction to Python] How to get the index of data with a for statement
Read the Python-Markdown source: How to create a parser
How to know the port number of the xinetd service
How to write a GUI using the maya command
How to get the number of digits in Python
A memo of how to use AIST supercomputer ABCI
A memo to visually understand the axis of pandas.Panel
How to create a submenu with the [Blender] plugin
How to visualize the decision tree model of scikit-learn
How to write a list / dictionary type of Python3
Steps to calculate the likelihood of a normal distribution
[Blender] How to dynamically set the selection of EnumProperty
Basics of PyTorch (2) -How to make a neural network-
[Python] Summary of how to specify the color of the figure
How to hit the document of Magic Function (Line Magic)
How to access the global variable of the imported module
How to post a ticket from the Shogun API
[Selenium] How to specify the relative path of chromedriver?
Python Note: The mystery of assigning a variable to a variable
How to display the modification date of a file in C language up to nanoseconds
How to identify the element with the smallest number of characters in a Python list?
I can't find the clocksource tsc! ?? The story of trying to write a kernel patch
[Ruby] How to replace only a part of the string matched by the regular expression?
How to check in Python if one of the elements of a list is in another list
Find a guideline for the number of processes / threads to set in the application server
The world's most easy-to-understand explanation of how to make a LINE BOT (1) [Account preparation]