Visualization of wavelet

Leave a comment

Many wavelet does not have functional form, but defined by the MRA coefficient.

The visualization of wavelet can be done by using wavelet construction.

\displaystyle v_{j+1,k} = \sum_{n} g_0(k-2n)v_{j,n} + g_1(k-2n) w_{j,n}

For scaling function, we can define v_0 = {1} and w_0 = {0} .

\displaystyle v_{1,k} = \sum_{n} g_0(k-2n)v_{0,n} = g_0(k)

Similarly, the wavelet can be started with v_0 = {0} and w_0 = {1} .

\displaystyle v_{1,k} = \sum_{n} g_1(k-2n)w_{0,n} = g_1(k)

Then build by iteration,

\displaystyle v_{j+1,k} = \sum_{n} g_0(k-2n) v_{j,n}

From last post on the MRA coefficient, i calculated and plot the wavelet for m = 4.

Capture.PNG

we can see the wavelet becomes the Haar wavelet as the free parameter goes to 1. In fact, it becomes a shifted Haar wavelet when the free parameter goes to 0, as we can imagine.

When the free parameter is 0.683013, it is the Daubechies-2 wavelet. Notes that some people will absorbed a factor $latex 1/ \sqrt{2} $ into the coefficient, so that their free parameter is 0.683013/\sqrt{2} = 0.482963.

Orthonormal MRA Coefficient

Leave a comment

A multi-resolution analysis is defined by scaling function and the corresponding wavelet. From the scaling relations

\displaystyle \phi_{j,k}(x) = \frac{1}{\sqrt{2}} \sum_{l} g_0(l) \phi_{j+1,2k+l}(x)

\displaystyle \psi_{j,k}(x) = \frac{1}{\sqrt{2}} \sum_{l} g_1(l) \phi_{j+1,2k+l}(x)

the scaling function and wavelet can be defined from the MRA coefficient g_0, g_1

The coefficients are constrained due to the properties of wavelet and scaling function.

\displaystyle \int \phi(x) dx = 1

\displaystyle \int \phi_{j,k}(x) \phi_{j,k'}(x) dx = \delta_{kk'}

\displaystyle \int \psi(x) dx = 0

\displaystyle \int \psi_{j,k}(x) \psi_{j,k'}(x) dx = \delta_{kk'}

\displaystyle \int \psi_{j,k}(x) \phi_{j,k'}(x) dx = 0

These properties lead to

\displaystyle \sum g_0(l) = 2

\displaystyle \sum g_1(l) = 0

\displaystyle \sum_{l,n} g_0(l) g_0(l+2n) = \begin{matrix} 2, & n=0 \\ 0, & else  \end{matrix}

\displaystyle \sum_{l,n} g_1(l) g_1(l+2n) = \begin{matrix} 2, & n=0 \\ 0, & else  \end{matrix}

\displaystyle \sum_{l,n} g_0(l) g_1(l+2n) = 0

The 3rd and 4th constrains requires the numbers of non-zero element in g_0, g_1 are even.

One of the solution is setting

g_1(k) = (-1)^k g_0 (1-k)

so that we don’t need to worry g_1 and the 4th constrain becomes the 3rd constrain, and the 5th constrain is always satisfied. Now, only the 1st, 2nd, and 3rd constrains are needed. This is equivalent to 1+m/2 equations with number of non-zero elements in g_0 is m.

m 1 + \frac{m}{2} Degree of Freedom
2 2 0
4 3 1
6 4 2
8 5 3

For size of 4, the solution is

\displaystyle g_0 = \left(a,  \frac{1-\sqrt{1+4a-4a^2}}{2},  1-a,  \frac{1+\sqrt{1+4a-4a^2}}{2}\right)

In fact, the coefficient for g_0 can be grouped as even and odd, so that

\displaystyle \sum g_0(2k) = \sum g_0(2k+1) = 1

and the constrain 3rd can lead to,

\displaystyle (\sum g_0(2k))^2 + (\sum g_0(2k+1)^2 = 2 ,

which is automatically fulfill.

 

Hard & soft thresholding

Leave a comment

In usual Fourier transform (FT), the filter is cut-off certain frequency.

This trick is also suitable for wavelet transform (WT). However, there could be some “features” located in high frequency scale (or octave) , a simply cut-off would remove these features.

If the signal to noise level is large, that means the noise has smaller amplitude than that the signal, we can use hard or soft thresholding, which zero any coefficient, which is after the FT or WT,  less then a threshold.

Lets X be the coefficient. The hard thresholding is

Y=\begin{cases} 0, & |X| <\sigma \\ X, & \mbox{else} \end{cases}

ht

The soft thresholding is

Y = \begin{cases} 0, & |X| < \sigma \\ sign(X) f(|X|, \sigma), & \mbox{else} \end{cases}

A popular function

\displaystyle f(x, \sigma) = \frac{x - \sigma}{ X_{max} - \sigma } X_{max}

st.PNG

or

\displaystyle f(x, \sigma) = x - \sigma

st2.PNG

Coulomb Energy in Nuclear Unit

Leave a comment

Coulomb energy between 2 charge q_1, q_2 in SI unit is

\displaystyle U_c = \frac{1}{4\pi\epsilon_0} \frac{q_1 q_2}{r}  [J]

\displaystyle U_c = \frac{e^2}{4\pi\epsilon_0} \frac{Z_1 Z_2}{r [m]} [J \cdot m] = \frac{Z_1 Z_2}{r [m]} (\frac{e^2}{4\pi \epsilon_0}) [J \cdot m] = \frac{Z_1 Z_2}{r [m]} 2.30708 \times 10^{-28} [J \cdot m]

We need to convert the SI unit into nuclear unit:

1 [J] = \frac{10^{-6}}{e} [MeV]

1 [m] = 10^{15} [fm]

Since the unit of the r depends on the nominator, the change of the unit does not need to be compensated.

\displaystyle U_c = \frac{Z_1 Z_2}{r [fm]} 1.43996 [MeV \cdot fm]

Therefore, a simple expression

\displaystyle U_c = \frac{e^2}{r} Z_1 Z_2

where e^2 = 1.44 [MeV\cdot fm]

Other useful quantities are:

  • \hbar c = 197.327 [MeV\cdot fm]
  • e^2/\hbar c = 1/137.036
  • \hbar = 6.58212 [MeV\cdot s]
  • c = 2.99792458\times 10^{23} [fm/s]

Qt with FFTW

Leave a comment

To use Qt with FFTW (Fast Fourier Transfrom in the West) in windows.

  1. download the FFTW pre-compiled *.dll file
    http://www.fftw.org/install/windows.html
  2. extract the zip. Copy the fftw3.h and libfftw3-3.dll and put them into your project folder, the folder same as *.pro file.
  3. in *.pro, add
    LIBS += "$$PWD/libfftw3-3.dll"
  4. in QT creator, on the “projects” column, right click on the project name and select “Add Existing Files…”, choose the fftw.h
  5. in main.cpp, add
    #include "fftw.h"
  6. got to the build directory, usually named as “Build-…..”, you will see “debug” or “release” folders. copy the libfftw3-3.dll and paste it in the “debug” or “release” folder depends on you are debug or release.

The step 6 is very important, without so, the program will not even started, while the compilation is completely normal.

Electron spin resonance in power sample

Leave a comment

The topic was discussed in a paper “Dynamic nuclear polarization by photoexcited-triplet electron spins in polycrystalline samples” by K. Takeda, K. Takegoshi, and T. Terao. I am going to give some supplementary on how the calculation can be done.


In a spin-1 system, the spin operator S_x, S_y, S_z and the rotation operator are easy to calculate.

The zero-field splitting Hamiltonian in the ordinary basis of S_z is

\displaystyle H_{zfs} = \begin{pmatrix} D/3 & 0 & E \\ 0 & -2D/3 & 0 \\ E & 0 & D/3 \end{pmatrix}

The eigen vectors in the paper are

\displaystyle \left|X\right> = \frac{1}{\sqrt{2}}( \left|-\right> - \left|+\right> )

\displaystyle \left|Y\right> = \frac{i}{\sqrt{2}}( \left|-\right> + \left|+\right> )

\displaystyle \left|Z\right> = \left|0\right>

with eigen values X = D/3 -E, Y = D/3 +E, and Z = -2D/3 . The convention that the \left| Y \right> has an imaginary number i is not necessary. In the paper, there are typos on the eigen values.

In other way, the zero-field splitting Hamiltonian can be diagonalized into

\displaystyle H_{zfs} =  P\cdot H_{XYZ} \cdot P^{-1}

\displaystyle H_{XYZ} = \begin{pmatrix} X & 0 & 0 \\ 0 & Y & 0 \\ 0 & 0 & Z \end{pmatrix}

\displaystyle P = \frac{1}{\sqrt{2}} \begin{pmatrix} -1 & i & 0 \\ 0 & 0 & \sqrt{2} \\ 1 & i & 0 \end{pmatrix}

where P is the transformation changing from ordinary basis to eigen basis.


Given that the population of the eigen basis is w = ( w_X, w_Y, w_Z) , the density matrix is

\displaystyle \rho_{XYZ} = \sum_{i = X, Y, Z} \left| i \right>\left< i \right| w_i

The corresponding population in ordinary basis is the diagonal elements of

\rho = P \cdot \rho_{XYZ} \cdot P^{-1}

In the paper, it gives w = (0.76, 0.16, 0.08). In order to get the population when the molecule X-axis is aligned with the Z-axis Lab frame, the \rho has to be rotated by 90 degree, and we have the population w' = ( 0.12, 0.76, 0.12) .


The Zeeman Hamiltonian is H_B = \omega S_Z , where \omega = \gamma_e  B  \hbar . In the paper, the author fix the crystal and rotating the external magnetic field. To express the H_B into eigen basis. We simply apply a transformation as the density matrix. But first we need to rotate the magnetic field. The rotation is simple, as the direction of the magnetic field is (\theta, \phi) . Thus, the Zeeman Hamiltonian is H_B = \omega \hat{n} \cdot \vec{S} .

The transformation is

H_{B_XYZ} = P^{-1} \cdot H_B \cdot P

again, the total Hamiltonian is difference from the author.

\displaystyle H = \begin{pmatrix} X & -i \omega cos(\theta) & i \omega cos(\phi) sin(\theta) \\ i \omega cos(\theta) & Y & -i\omega sin(\phi) sin(\theta) \\ -i\omega cos(\phi) sin(\theta) & i \omega sin(\phi) sin(\theta) & Z \end{pmatrix}

The difference is because the normal vector \hat{n} = ( sin(\phi) sin(\theta), cos(\phi) sin(\theta), cos(\theta) ) is used in the paper.

After that, the calculation for the power sample is the same.


Instead of using integration, we can generated the power using isotropic distribution. I generated 10k molecule, the result of 9.65 GHz microwave is

power_ESR.PNG

where the y axis is the magnetic field in Tesla. The result is similar to the paper. I think that their calculation somehow changed the eigen value, so that the eigen vectors were matching with the eigen energy. And the notation of the normal vector and the constant factor of the eigen vectors does not matter.

Phase shift of elastics scattering

Leave a comment

I found that the derivation of most “google result” is not clear enough. So here is my derivation. Before process, people may need to review the pervious post.

Most people start from,

u_l = A( \hat{h}_l^- - S_l \hat{h}_l^+ )

that annoying me because it is somehow not “natural”. Why there is a “minus” sign? Why the \hat{h}_l^- is the first term? For my self, a more natural way is,

u_l = a \hat{h}_l^+ + b \hat{h}_l^-

where a, b are complex numbers, but that is still not so natural, because in numerical calculation, for simplicity, there is no complex number, we only have,

u_l = \alpha \hat{j}_l + \beta \hat{n}_l

The first term is alway there as it is the free solution and bounded at r = 0. the second term is caused by the potential.


The goal is to find a solution take the form

\displaystyle \psi = A \left( e^{i \vec{k} \cdot \vec{r}} + f(\theta) \frac{e^{ikr}}{r} \right)

where the first term is free wave and the second term is scattered wave. The solution for elastics scattering is

\displaystyle \psi = \sum C_l P_l (\cos\theta) \frac{u_l}{kr} \rightarrow \sum C_l P_l(\cos\theta) (\alpha \hat{j}_l + \beta \hat{n}_l)

we used the substitution,

\displaystyle R_l(r) = \frac{u_l(\rho)}{\rho},  \rho = kr.

The radial function can be solved using Rungu-Kutta method on the equation,

\displaystyle \frac{d^2}{d\rho^2} u_l = \frac{2 m_\mu}{\hbar^2} (V-E) u_l + \frac{l(l+1)}{\rho^2}

and the solution of u_l at far away is,

u_l \rightarrow  \alpha \hat{j}_l + \beta \hat{n}_l.

the arrow means r \rightarrow \infty. So, the problem is how to rewrite the solution. In the way, we will see how the phase shift or the S-matrix was found.


The free solution is the spherical wave,

\displaystyle e^{i \vec{k} \cdot \vec{r}} = \sum_l (2l+1) i^l P_l(\cos\theta) j_l(kr)

The spherical Bessel function j_l(kr) cna be express as Heankel function

h_l^{\pm} = n_l \pm i j_l \rightarrow e^{\pm i (kr - l \frac{\pi}{2})}

The + sign is outgoing wave.


\displaystyle u_l \rightarrow (\alpha \hat{j}_l + \beta \hat{n}_l)

\displaystyle = \frac{\alpha}{2i} (\hat{h}_l^{+} - \hat{h}_l^{-}) + \frac{\beta}{2}(\hat{h}_l^{+} + \hat{h}_l^{-})

\displaystyle = \frac{\alpha + i \beta}{2i} (\hat{h}_l^{+} - \hat{h}_l^{-}) + \beta \hat{h}_l^{+}

\displaystyle = (\alpha - i \beta ) \left( \frac{\hat{h}_l^+ - \hat{h}_l^-}{2i} + \frac{\beta}{\alpha - i \beta} \hat{h}_l^+\right)

\displaystyle = (\alpha - i \beta ) \left( \hat{j}_l + \frac{\beta}{\alpha - i \beta} \hat{h}_l^+\right)

Since the u_l should be normalized, we can se \alpha = \cos \delta and \beta = \sin\delta.

\displaystyle \frac{\beta}{\alpha - i \beta } = \sin(\delta) e^{i\delta}

We put u_l back

\displaystyle \psi \rightarrow \sum_l C_l P_l (cos\theta)(\alpha - i \beta ) \left( \hat{j}_l + \sin(\delta) e^{i\delta} \hat{h}_l^+\right)

By setting

\displaystyle C_l = A i^l \frac{2l+1}{\alpha - i \beta} ,

we have the first term is the free wave function. In the second term, \hat{h}_l^+ \rightarrow e^{i(kr - l \frac{\pi}{2}}) / kr . Notice that

e^{i l \frac{\pi}{2}} = i^{-l}

That cancel the i^l term in C_l . And we have

 \displaystyle f(\theta) = \sum (2l+1) P_l (\cos\theta) \frac{\sin(\delta) e^{i\delta}}{k}


some people will write the u_l as \hat{h}_l^{\pm} and the S-matrix,

\displaystyle u_l = \frac{\alpha + i \beta} {2i} \hat{h}_l^+ - \frac{\alpha - i \beta}{2i} \hat{h}_l^-

\displaystyle = -\frac{\alpha - i \beta}{2i} \left( \hat{h}_l^- - \frac{\alpha + i \beta}{\alpha - i \beta} \hat{h}_l^+ \right)

\displaystyle = A' (\hat{h}_l^- - S_l \hat{h}_l^+)

where

\displaystyle S_l =\frac{\alpha + i \beta}{\alpha - i \beta} = e^{2i\delta} .

Remember that this is the S-matrix for elastics scattering.

 

Older Entries