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}


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


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.





Discrete Fourier Transform (shift of the spectrum)

Leave a comment

as the last post shown the formula of discrete Fourier transform:

a_n = \sum_{i=0}^m { x_i e^{i (n-1) 2\pi (i-1) /m}}

where (m+1) is the number of data-point, n is from 1 to (m+1). if n is bigger then m+1, it is a_n is start to repeat. for example, n = m+1+r

a_{m+1+r} = a_{r}


a_{m+1-r} = a_{-r}

Thus, we can rearrange the a_n from -(m+1)/2  to (m+1)/2. This, set the center be zero-frequency. for example, if m+1 = 1000, then we can have n from 1 to 1000, or -499 to 0 and 1 to 500. since n=1 mean frequency = 0, we can shift n by 1 unit and get n from -500 to 499.

this it very useful if out signal is actually from 2 signals, 1 is raw signal and the other is a reference signal. the signal we get is by subtracting these 2. thus, when the raw signal is very close to the reference signal, the frequency is very small and the peak of Fourier spectrum will go to the edge and hardly recognized. However, by shift the spectrum and set the a_0 at the center. we can see How the raw signal is different from the reference signal.


Discrete Fourier Transform

Leave a comment

we usually record a set of data in time-domain. the signal or data may have some oscillation and we want to to which frequency is inside. thus, we will looking for Fourier Transform, which change from time-domain to frequency domain.

the continuous Fourier Transform is

F(k) = \frac{1}{\sqrt{2\pi}} \int_{\infty}^{\infty}{ f(x) e^{-i k x} dx}

if we see the exponential as cosine and sine function, that we understand it is finding a suitable wave-vector “k”, such that the function f(x) is “parallel” or “part-parallel” with it.

However, we cannot have infinite long data, but just some interval. Thus, the interval will bounded the available frequency.

a_n = \frac{1}{\sqrt{2\pi}} \int_{-L}^{L}{ f(x) e^{-i n \pi x / L} dx}

we can understand the reason is by invert thinking about how a signal is formed by difference frequency. say, we have a time interval of  “2L” unit (if the 2π is absorbed in the interval, the total interval is 2π larger), if a mono-frequency signal has wavelength longer then “2L” unit, or the k-vector is smaller then 1/2L, then, this signal is not going to change within the time-interval and do no contribution. Thus, the interval set the resolution of the k-vector.

when the data is discrete with m+1 data point in interval 2L, we have discrete Fourier transform. the integration becomes:

a_n = \sum_{i=0}^{m} { f(x_i) e^{-i n \pi x_i /L}} \Delta x

Set the \Delta x = 2L/m = 1 , and x_i = (i-1)2L/m , and since the data separation on the data list is defined by ourselves, so, the date does not have this information. And we want the a_n start at 1, not zero.

a_n = \sum_{i=0}^m { f(x_i) e^{-i (n-1) 2\pi (i-1) /m}}

This is the formula for Discrete Fourier Transform without a constant factor.

Remark on mathematica and FFTW user: The function “Fourier”‘s default setting is using – sign in the exponential, while FFTW is using – sign. and also, the amplitude will be divided by √n in mathematica, and FFTW will not. in this post, we use FFTW transform. For mathematica user, the “Fourier” has a setting called “FourierParameters”, by setting it to (1,-1), will get the result same as FFTW.

Differential Cross Section

Leave a comment

In nuclear physics, cross section is a raw data from experiment. Or more precisely differential cross section, which is some angle of the cross section, coz we cannot measure every scatter angle and the differential cross section gives us more detail on how the scattering going on.

The differential cross section (d.s.c.) is the square of the scattering amplitude of the scatter spherical wave, which is the Fourier transform of the density.

d.s.c = |f(\theta)|^2 = Fourier ( \rho (r), \Delta p , r )

Where the angle θ come from the momentum change. So, sometime we will see the graph is plotted against momentum change instead of angle.

By measuring the yield of different angle. Yield is the intensity of scattered particle. We can plot a graph of the Form factor, and then find out the density of the nuclear or particle.

However, the density is not in usual meaning, it depends on what kind of particle we are using as detector. For example, if we use electron, which is carry elected charge, than it can feel the coulomb potential by the proton and it reflected on the “density”, so we can think it is kind of charge density.

Another cross section is the total cross section, which is sum over the d.s.c. in all angle. Thus, the plot always is against energy. This plot give us the spectrum of the particle, like excitation energy, different energy levels.

decay time constant and line width

Leave a comment

the spectrum of energy always has a peak and a line width.

the reason for the line width is, this is decay.

i give 2 explanations, once is from classical point of view and i skipped the explanation for the imaginary part. so, i am not fully understand. the 2nd explanation is look better, but it is from QM. however, there is one hide question for that explanation is, why the imaginary energy is negative?

the simplest understanding of the relation is using fourier transform. (i think)

Fourier transform is changing the time-frame into the frequency frame. i.e, i have a wave, propagating with frequency w. we can see a wave shape when plot with time. and we only see a line, when we plot with frequency, since there is only 1 single frequency. however, for a general wave, it is composite of many different frequencies, using fourier transform can tell us which frequency are involved. And energy is proportional to frequency.

when the particle or state under decay. the function is like

f(t) = Exp(-R t) Exp ( i \omega_0 t)

where the R is decay constant, and ω0 is the wave frequency.

after fourier transform, assume there is nothing for t < 0

F(t) = \frac {1} { R + i ( \omega_0 - \omega )}

the real part is

Re(F(t)) = \frac {R} { R^2 + ( \omega_0 - \omega )^2}

which is a Lorentzian shape and have Full-Width-Half-Maximum (FWHM) is 2R. it comes from the cosine part of the fourier transform. thus, the real part.

and the imaginary part is

Im(F(t)) = \frac {\omega_0 - \omega}{R^2 + ( \omega_0 -\omega )^2 }

the imaginary part is corresponding to the since part, so, we can neglect it. (how exactly why we can neglect it? )

Thus, we can see, if there is no decay, R → 0, thus, there is no line width.

therefore, we can see the line width in atomic transition, say, 2p to 1s. but there are many other mechanism to the line width, like Doppler broadening, or power broadening. So, Decay will product line width, but not every line width is from decay.


another view of this relation is from the quantum mechanics.

the solution of Schroedinger equation is

\Psi (x,t) = \phi(x) Exp \left( - i \frac {E}{\hbar} t \right)

so, the probability conserved with time, i.e.:

|\Psi(x,t)|^2 = |\Psi (x,0)|^2

if we assume the energy has small imaginary part

E = E_0 - \frac {i} {2} R \hbar

( why the imaginary energy is nagative?)

|\Psi(x,t)|^2 = |\Psi (x,0)|^2 Exp ( - R t)

that make the wavefunction be :

\Psi (x,t) = \phi (x) Exp( - i \frac {E}{\hbar} t ) Exp( - \frac {R}{2} t )

what is the meaning of the imaginary energy?

the wave function is on time-domain, but what is “physical”, or observable is in Energy -domain. so, we want Psi[x,E] rather then Psi[x,t], the way to do the transform is by fourier transform.

and after the transform, the probability of finding particle at energy E is given by

|\Psi(x,E)|^2 = \frac {Const.}{R^2 +(\omega_0 - \omega )^2}

which give out the line width in energy.

and the relation between the FWHM(line width) and the decay time is

mean life time ≥ hbar / FWHM

which once again verify the uncertainty principle.