Algorithm of Wavelet Transform

Leave a comment

There are many kind of wavelet transform, and I think the names are quite confusing.

For instance, there are continuous and discrete wavelet transforms, in which, the “continuous” and “discrete” are for the wavelet parameters, not for the “data” itself. Therefore, for discrete data, there are “continuous” and “discrete” wavelet transforms, and for function, there are also “continuous” and “discrete” wavelet transforms.

In here, we will focus on discrete wavelet transform for function first. This discrete wavelet transform is also called as wavelet series, which express a compact support function into series of wavelet.

For simplicity, we also focus on orthonormal wavelet.

As the wavelet span the entire space, any compact function can be expressed as

f(t) = \sum_{j,k} \left<f(t)|\psi_{j,k}(t)\right> \psi_{j,k}(t)

\psi_{j,k}(t) = 2^{j/2} \psi(2^j t - k)

where j, k are integer.

Now, we move to discrete data discrete wavelet transform. The data is discrete, we can imagine only t_n = t_0 + n \Delta points are known with finite n .

f_n = f(t_n) = \sum_{j,k} \left<f_n|\psi_{j,k}(t_n) \right> \psi_{j,k}(t_n)

the integration becomes a finite sum.

Without loss of generality, we can set t_0 = 0, \Delta = 1, and then the time axis becomes an integer number axis. We found that j  \leq 0 as the wavelet can only be expand, not shrink. Because there are finite number of data point, i.e. n < \infty, -Log_2(n) < j \leq 0 .

However, this double summation for each f_n is very time consuming. There is a Fast Discrete Wavelet Transform. Before we continuous, we must study the wavelet.

From the last post, we know that the scaling function that generate a MRA must be:

\phi(t) = \sum_{k} g_0(k) \phi(2t-k)

\left<\phi(t-k) | \phi(t-k') \right> = \delta_{kk'}

, where k are integer. The set of shifted scaling function span a space V_0 . For the wavelet,

\psi(t) = \sum_{k} g_1(k) \psi(2t-k)

\left<\psi(t-k) | \psi(t-k') \right> = \delta_{kk'}

The set of shifted wavelet span a space W_0, so that W_0 \perp V_0, so that

\left<\phi(t-k)|\psi(t-k') \right> = 0

Since the wavelet is generated from the scaling function, we expect the coefficient of g_0(k) and g_1(k) are related. In fact, the relationship for orthonormal scaling function and wavelet is

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

I will continuous later….. very tired….


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]

\displaystyle U_c = \frac{Z_1 Z_2}{r [m]} 1.43996 \times 10^{-09} [MeV \cdot m]

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


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

but that is not so natural at all, 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^+)


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

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


Cross section IV

Leave a comment

A more microscopic view of the cross section is

\displaystyle  \frac{dN}{dS} \sigma \rho dt dS= dn,

where N is number of incident particle, \sigma is the total cross section, \rho is the particle density in the target, S is the area, t is the thickness of the target, and n is particle detected.

The cross section is a constant, so, after integration,

\displaystyle \sigma = \frac{n}{N \rho t} .

If the beam is not uniform, dN/dS = f(S) is a function of S. The integration has to be careful.


A reminder of Exponential function

Leave a comment

Suddenly found that exponential function is very wonderful. It is the eigen-function of differential operator. because of that, it appears in almost everywhere as differential operator constantly appears.

The definition of the exponential e is,

\displaystyle e = \lim_{n\to\infty} \left(1+\frac{1}{n} \right)^n

this is the 1% growth continually in a unit of time. It is easy to see that

\displaystyle e^x =  \lim_{n\to\infty} \left(1+\frac{x}{n} \right)^n

e^x e^y = e^{x+y}


\displaystyle \left(1+\frac{x}{n} \right)^n = \sum_{r=0}^{n} C_r^n \left(\frac{x}{n} \right)^r  = \sum_{r=0}^{n} \left(\prod_{k=0}^{r-1} \frac{n-k}{n} \right)\frac{x^r}{r!},

Take n\to\infty,

\displaystyle e^x = \sum_{r=0}^{\infty} \frac{x^r}{r!}

It is easy to see that

\displaystyle \frac{d}{dx} e^x = e^x


Older Entries