Wavelet | Nuclear Physics 101

A little summery of what I am doing

August 31, 2017

GoLuckyRyan atomic, Basic Clebsch-Gordon, Fourier transform, hydrogen, Laguerre, rotation, Schrödinger equation, Wavelet Leave a comment

My current position is developing a computational program that can measure the system of polarized target automatically and repeatedly. The program needs to connect with the microwave generator, voltage supply, power meter, and oscilloscope. That is purely technical and nothing special.

Later, after the program started to collecting data, another program is needed to analysis the data. Although there are many program that can do analysis, but those program are not so easy to use, in both control and display. So, I make an analysis program. The program is simply read the 2-D time-domain signal, and then determine some parameters of a specific function that fit the signal. However, the signal could be very noisy, so FFTW and wavelet analysis were implemented. That’s why the wavelet analysis appeared in this blog.

After that, the sample of polarized target is made from various pentacene derivatives, that no known energy level exist. One way to find out the energy levels of the singlet excited states is to measure the absorption spectrum. However, the energy of the triplet state is difficult to measure. And the energy level of the triplet state is critical for Dynamic Nuclear Polarization to be happened. Thus, one way to find out is doing computation chemistry.

My last chemistry class was like 15 years ago, when I was junior high school. But the basic of computation chemistry is solving the Schrodinger equation. That is what theoretical nuclear physicists do! The variation method, the Hartree-fock method, I heard it and somehow know it, but never do it with my own hand and computer. That is why I revisit the Hartree-fock method, and found out my previous understanding is so naive.

In the course of studying Hartree-Fock method, and one problem is evaluating the overlap integral

$\displaystyle \int \psi_a(r_1) \psi_b(r_2) \frac{1}{r_{12}} \psi_c(r_1) \psi_d(r_2) dr_1dr_2$

For a 3-D system, the integral involves product of multiple spherical harmonics. That is really troublesome. Therefore, I move to study the spherical harmonics, and the related rotational invariant, Wigner D-matrix, Clebcsh-Gordon series and Fourier series. The spherical harmonics arises from solving the Laplace equation in spherical coordinate. A general theory of the solution of Lapalace equation involves Legendre polynomial, which is a special case for Hypergeometric function. And a very interesting connection is that the elliptical function of the 1st and 2nd kind are also two special cases of hypergeometric function, that, the solution of Laplace equation for elliptical boundary condition is elliptical functions. That connects spherical harmonics and elliptical function! Wow!

I am now a bit off-track, that I am very interesting on the function of all functions. We know that there are many elementary functions, such as sin, cos, and Log. And even more kind of special functions, such as

Hermite — solving 1-D harmonic oscillator
Laguerre — the radial function of hydrogen
Legendre — the solution of the “ $\theta$ ” of Laplace equation
Gamma — a continuation of factorial
Bessel — solution of 3-D infinite square well
Elliptic — magnetic field of a solenoid
Dirac delta
Gaussian or Error function

For discrete argument

Clebsch-Gordon
Factoral
Binomial

As far as I know, the Hypergemetric function is like “mother of functions”, although not all special functions can be expressed as it. I am driven by curiosity, so, I am not sure where I will go. For instance, the transformation of Hypergeometric function is very interesting.

So, for now, as an ending, I found one article is very interesting. The History and Future of Special Functions, by Stephen Wolfram.

A wavelet presentation

May 18, 2017

GoLuckyRyan Basic Wavelet Leave a comment

Wavelet Transformation

Recently, I gave a mini-lecture on wavelet analysis for my colleagues. This is a 30-min compact lecture for introduction and the application.

Enjoy! and feel free to leave comments if you have questions.

Visualization of wavelet

April 29, 2017

GoLuckyRyan Basic Daubechies, MRA, 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 scaling coefficient, i calculated and plot the wavelet for $m = 4$ .

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

April 28, 2017

GoLuckyRyan Basic MRA, scaling coefficient, Wavelet 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 scaling 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

April 13, 2017

GoLuckyRyan Basic denoising, FDWT, FFT, Fourier transform, threshold, tran, Wavelet 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}$

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

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

Algorithm of Wavelet Transform (with Qt class)

March 24, 2017

GoLuckyRyan computer, Math, theory FDWT, Qt, transform, Wavelet 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

$\displaystyle 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$ .

$\displaystyle 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:

$\displaystyle \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,

$\displaystyle \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)$

For discrete data $x_i$ , it can be decomposed into the MRA space. We start by the largest $V_0$ space, where the wavelet is most shrunken.

$\displaystyle x_i = \sum_{k} v_{0,k} \phi(i-k)$

to decompose to the $V_{-1}$ and $W_{-1}$ space. We can use the nested property of the MRA space, $\phi(2t)$ can be decomposed into $\phi(t-k)$ and $\psi(t-k)$ ,

$\displaystyle \psi(2t-l) = \sum_{k} h_0(2k-l) \phi(t-k) + h_1(2k-l) \psi(t-k)$

where (given that $\phi(t)$ and $\psi(t)$ are orthonormal ),

$h_0(2k-l) = \left< \phi(2t-l) | \phi(t-k) \right>$

$h_1(2k-l) = \left< \phi(2t-l) | \psi(t-k) \right>$

Therefore, using the coefficient of $h_0$ and $h_1$ , the wavelet coefficient $v_{0,k}$ can be decomposed to

$\displaystyle v_{s-1,k} = \sum_{l} h_0(2k-l) v_{s,l}$

$\displaystyle w_{s-1,k} = \sum_{l} h_1(2k-l) v_{s,l}$

in graphic representation

This is a fast discrete wavelet transform.

Due to the nested space of MRA, we also expect that the coefficient $h_0$ and $h_1$ are related to $g_0$ . For orthonormal wavelet,

$\displaystyle h_0(k) = \frac{1}{2} g_0(-k)$

$\displaystyle h_1(k) = \frac{1}{2} (-1)^{k} g_0 (k+1)$

Since the $g_0$ is finite, the $g_1, h_0, h_1$ are all finite. That greatly reduce the computation cost of the discrete wavelet transform.

To reconstruct the discrete data $x_i$ , we don’t need to use

$\displaystyle v_{s+1,l} = \sum_{k} v_{s,k} \phi(l - k) + w_{s,k} \psi(l-k)$

using the nested space of MRA, $\psi(t) = \sum_{k} g_1(k) \psi(2t-k)$ ,

$\displaystyle v_{s+1,l} = \sum_{k} g_0(l-2k) v_{s,k} + g_1(l-2k) w_{s,k}$

in graphical representation,

I attached the wavelet transfrom class for Qt, feel free to modify.

https://drive.google.com/file/d/0BycN9tiDv9kmMVRfcVFMbDFWS0k/view?usp=sharing

https://drive.google.com/file/d/0BycN9tiDv9kmUmlmNk1kaVJCbEU/view?usp=sharing

in the code, the data did not transform to MRA space. The code treats the data already in the MRA space. Some people said this is a “crime”. But for the seek of “speed”, it is no need to map the original discrete data into MRA space. But i agree, for continuous function, we must map to MRA space.

Wavelet Analysis or MRA

March 21, 2017

GoLuckyRyan Math, theory FFTW, Fourier transform, MRA, Short time Fourier Transfrom, STFT, uncertainly, Wavelet 2 Comments

Although the Fourier transform is a very powerful tool for data analysis, it has some limit due to lack of time information. From physics point of view, any time-data should live in time-frequency space. Since the Fourier transform has very narrow frequency resolution, according to uncertainty principle, the time resolution will be very large, therefore, no time information can be given by Fourier transform.

Usually, such limitation would not be a problem. However, when analysis musics, long term performance of a device, or seismic survey, time information is very crucial.

To over come this difficulty, there a short-time Fourier transform (STFT) was developed. The idea is the applied a time-window (a piecewise uniform function, or Gaussian) on the data first, then FT. By applying the time-window on difference time of the data (or shifting the window), we can get the time information. However, since the frequency range of the time-window always covers the low frequency, this means the high frequency signal is hard to extract.

To improve the STFT, the time-window can be scaled (usually by 2). When the time window is shrink by factor of 2, the frequency range is expanded by factor of 2. If we can subtract the frequency ranges for the time-window and the shrink-time-window, the high frequency range is isolated.

To be more clear, let say the time-window function be

$\phi_{[0,1)}(t) = 1 , 0 \leq t < 1$

its FT is

$\hat{\phi}(\omega) = sinc(\pi \omega)$

Lets also define a dilation operator

$Df(t) = \sqrt{2} f(2t)$

the factor $\sqrt{2}$ is for normalization.

The FT of $D\phi(t)$ has smaller frequency range, like the following graph.

We can subtract the orange and blue curve to get the green curve. Then FT back the green curve to get the high-frequency time-window.

We can see that, we can repeat the dilation, or anti-dilation infinite time. Because of this, we can drop the FT basis $Exp(-2\pi i t \omega)$ , only use the low-pass time-window to see the low-frequency behaviour of the data, and use the high-pass time-window to see the high-frequency behaviour of the data. Now, we stepped into the Multi-resolution analysis (MRA).

In MRA, the low-pass time-window is called scaling function $\phi(t)$ , and the high-pass time-window is called wavelet $\psi(t)$ .

Since the scaling function is craetd by dilation, it has the property

$\phi(t) = \sum_{k} g_{0}(k) \phi(2t-k)$

where $k$ is integer. This means the vector space span by ${\phi(t-k)}_{k}=V_0$ is a subspace of the dilated space $DV_0 =V_1$ . The dilation can be go one forever, so that the whole frequency domain will be covered by $V_{\infty}$ .

Also, the space span by the wavelet, ${\psi(t-k)}=W_0$ , is also a subspace of $V_1$ . Thus, we can imagine the structure of MRA is:

Therefore, any function $f(t)$ can also be expressed into the wavelet spaces. i.e.

$f(t) = \sum_{j,k} w_{j,k} 2^{j/2}\psi(2^j t - k)$

where $j, k$ are integers.

I know this introduction is very rough, but it gives a relatively smooth transition from FT to WT (wavelet transform), when compare to the available material on the web.

[Pol. p target] Meeting report (June 8th)

June 8, 2011

GoLuckyRyan Lab Log, others denoising, Pol. p target, Wavelet Leave a comment

Done

Calibrating NMR system with water

by changing NMR level
after found the peak on level = 150, measure the FID area in successively.

Result:

the pulse may not be optimized. the 1st pulse gave FID area = 25, and 2nd pulse gave FID area = 7, 3rd pulse gave 5.

Observed the Coil Relaxation signal depends on the coil impedance. there is a characteristic peak to indicate the change in impedance. it can be used to measure the impedance.

the Coil was wrapped by Teflon tape and fastened the copper wire. Adhesive was used to fix the join of the coil and cable. the insertion mechanism was fixed by optics mounting.

Result:

the characteristic peak does not change so much. thus, we have confidence that the variation on NMR signal is NOT by the coil.

Optimization

the microwave delay time was measured. since it is not easy to have trigger on -10 us after the laser pulse end. we use 0us instead.
the microwave power is optimized at 1.0W
Laser pulse duty and chopper frequency.

Laser polarization angle

the change is smaller then signal fluctuation.

Wakui San comments

Laser mode

the laser is running at multiline mode, but the power detector is at 514.5nm

Crystal expired

the crystal we use is about 5 years ago

ToDo

have to measure the statistics of the data, due to a improvement of the coil.
Crystal orientation
Optimization

laser pulse duty

NMR pulse calibration by water
T1 and T2 measurement
Q-factor of the coil
more detail measurement on each parameters
thermal polarization

Discussion

the coil is being fixed, we are able to have a more reliable data. we have to find out the statistics. we can compare with a previous data @ Puw = 1.0 on June 3rd, we had collected 15 data for same setting and the s.d. is 30 unit.
To have the absolute polarization measurement

we have to lower the variation of signal
we have to lower the noise level
we have to get same setting for the NMR system

impedance
level
pulse time
gain
Forward and Backward power

the change of FID area due to change of external H-field

the data shown the angular frequency is pi per 30us, about 0.1 rad per 1us.
the angular frequency for proton is 267.5 rad per Tesla
if the field change for 1%
the change of the angular frequency will be 2.67 rad per us
or to say, the fluctuation of the field should be less then 0.05%

In order to preform Fourier transform, or the Finite time Fourier Transform, we can use wavelet analysis.
Before polarization transfer to 13C, we have to optimize the system.
the sample NMR signal is not flowing sine and cosine function

it it due to crystal field

Nuclear Physics 101

A little summery of what I am doing

A wavelet presentation

Visualization of wavelet

Orthonormal Scaling Coefficient

Hard & soft thresholding

Algorithm of Wavelet Transform (with Qt class)

Wavelet Analysis or MRA

[Pol. p target] Meeting report (June 8th)

Done

Wakui San comments

ToDo

Discussion

Recent Posts

Calendar

Archives

Keywords

# of Visitors

Login