## Orthonormal Scaling Coefficient

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$

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

## Absolute polarization measurement by elastic scattering

The magnitude of proton polarization can be measured by NMR technique with a reference. Because the NMR gives the free-induction decay signal, which is a voltage or current. For Boltzmann polarization using strong magnetic field and low temperature, the polarization can be calculated. However, when a reference point is not available, the absolute magnitude of proton polarization can be measured using proton-proton elastic scattering. The principle is the nuclear spin-orbital coupling. That creates left-right asymmetry on the scattering cross section.

Because of spin-orbital interaction:

$V_{ls}(r) = f(r) \vec{l} \cdot \vec{s} ,$

where $f(r)$ is the distance function, $\vec{l}$ is the relative angular momentum, $\vec{s}$ is the spin of the incident proton. In the following picture, the spin of the incident proton can be either out of the plane ($\uparrow$ ) or into the plan ($\downarrow$). When the proton coming above, the angular momentum is into the plane ($\downarrow$). The 4 possible sign of the spin-orbital interaction is shown. We can see, when the spin is up, the spin-orbital force repulses the proton above and attracts the proton below. That creates an asymmetry in the scattering cross section.

The cross section is distorted and characterized using analysing power $A_y$. Analyzing power is proportional to the difference between left-right cross-section. By symmetry (parity, time-reversal) consideration, $A_y = 1 + P sin(2\theta)$ (why?), in center of mass frame. In past post, the transformation between difference Lorentz frame. The angle in the $A_y$ has to be expressed in lab angle. The cross section and $A_y$ can be obtained from http://gwdac.phys.gwu.edu/ .

In scattering experiment, the number of proton (yield) is counted in left and right detectors. The yield should be difference when either proton is polarized. The yield is

$Y(\theta, \phi) = L \epsilon \sigma_0 (1 + cos(\phi)A_y(\theta) P) ,$

where $L$ is the luminosity, $\epsilon$ is the detector efficiency, $\sigma_0$ is the integrated cross-section of un-polarized beam and target of the detector, $P$ is the polarization of either the target or beam. When both target and the beam are polarized, the cross section is

$\sigma = \sigma_0 (1 + (P + P_T)A_y + P P_T C_yy),$

where $C_yy$ is spin-spin correlation due to spin-spin interaction of nuclear force.

Using the left-right yield difference, the absolute polarization of the target or the beam can be found using,

$\displaystyle A_y P = \frac{Y_L - Y_R}{Y_L + Y_R} ,$

where $Y_L = Y(\phi =0)$ and $Y_R = Y(\phi=\pi)$.

## Hard & soft thresholding

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

or

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

## Compiling Fortran-77 code in Ubuntu-16

Fortran-77 is a very old code, who lives in 32-bit computer.

In Ubuntu-16, the compiler g++, gcc, or gfortran are “basically the same” (as far as I understand, correct me if I am wrong.) that they only support fortran-95.

In order to compile Fortran-77 code, I tried many way, but the only way is install g77 from external source, and add -m32 for the compiling flag.

or, people can search in google by

 g77_x64_debian_and_ubuntu.tar.gz

people need to extract, change the mod of install.sh to be executable.

 tar -xzvf g77_x64_debian_and_ubuntu.tar.gz
cd g77_x64_debian_and_ubuntu
chmod +x ./install.sh
./install.sh

Somehow, you may face an error in apt-get, saying

Errors were encountered while processing:
g77-3.4-doc 

you can remove that by

cd /var/lib/dpkg/info
sudo rm g77-3.4-doc*
sudo dpkg --remove --force-remove-reinstreq g77-3.4-doc

Thanks (here)

Hope it help. :)

## Algorithm of Wavelet Transform (with Qt class)

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

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

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.

## Coulomb Energy in Nuclear Unit

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