diff --git a/notes/sections/6.md b/notes/sections/6.md index 351627f..f14d7a5 100644 --- a/notes/sections/6.md +++ b/notes/sections/6.md @@ -15,7 +15,7 @@ where: - $E$ is the electric field amplitude, default set $E = \SI{1e4}{V/m}$; - $a$ is the radius of the slit aperture, default set $a = \SI{0.01}{m}$; -- $\theta$ is the angle specified in @fig:fenditure; +- $\theta$ is the angle specified in @fig:slit; - $J_1$ is the Bessel function of first order; - $k$ is the wavenumber, default set $k = \SI{1e-4}{m^{-1}}$; - $L$ default set $L = \SI{1}{m}$. @@ -46,7 +46,7 @@ where: \node [cyclamen] at (5.5,-0.4) {$\theta$}; \node [rotate=-90] at (10.2,0) {screen}; \end{tikzpicture} -\caption{Fraunhofer diffraction.} +\caption{Fraunhofer diffraction.}\label{fig:slit} } \end{figure} @@ -101,11 +101,8 @@ same bin width of the previous one, but a smaller number of them ($m < n$), was filled with $m$ points according to a Gaussian distribution with mean $\mu$, corresponding to the central bin, and variance $\sigma$. Then, the original histogram was convolved with the kernel in order to obtain -the smeared signal. The result is shown in @fig:convolved. - -![Same sample of @fig:original convolved with the - kernel.](images/6_convolved.pdf){#fig:convolved} - +the smeared signal. Some results in terms of various $\sigma$ are shown in +@fig:convolved. The convolution was implemented as follow. Consider the definition of convolution of two functions $f(x)$ and $g(x)$: @@ -124,8 +121,9 @@ $$ $$ where $j$ runs over the bins of the kernel. -For a better understanding, see @fig:dot_conv. Thus, the third histogram was -obtained with $n + m - 1$ bins, a number greater than the initial one. +For a better understanding, see @fig:dot_conv. As can be seen, the third +histogram was obtained with $n + m - 1$ bins, a number greater than the initial +one. \begin{figure} \hypertarget{fig:dot_conv}{% @@ -194,6 +192,104 @@ obtained with $n + m - 1$ bins, a number greater than the initial one. } \end{figure} -\textcolor{red}{Missing various $\sigma$ comparison.} -## Unfolding +## Unfolding with FFT + +Two different unfolding routines were implemented, one of which exploiting the +Fast Fourier Transform. +This method is based on the property of the Fourier transform according to +which, given two functions $f(x)$ and $g(x)$: + +$$ + \hat{F}[f*g] = \hat{F}[f] \cdot \hat{F}[g] +$$ + +where $\hat{F}[\quad]$ stands for the Fourier transform of its argument. +Thus, the implementation of this tecnique lies in the computation of the +Fourier trasform of the signal and the kernel, the product of their transforms +and the anti-transformation of the result. Being the histogram a discrete set of +data, the Discrete Fourier Transform (DFT) was emploied. + +In order to accomplish this procedure, every histogram was transformed into a +vector. The kernel vector was 0-padded in order to make its length the same as +the one of the signal, making it feasable to implement the dot product between +the vectors. + +FFT are efficient algorithms for calculating the DFT. Namely, given a set of $n$ +values {$z_i$}, each one is transformed into: + +$$ + x_j = \sum_{k=0}^{n-1} z_k \exp \left( - \frac{2 \pi i j k}{n} \right) +$$ + +The evaluation of the DFT is a matrix-vector multiplication $W \vec{z}$. A +general matrix-vector multiplication takes $O(n^2)$ operations. FFT algorithms, +instad, use a divide-and-conquer strategy to factorize the matrix into smaller +sub-matrices. If $n$ can be factorized into a product of integers $n_1$, $n_2 +\ldots n_m$, then the DFT can be computed in $O(n \sum n_i) < O(n^2)$ +operations, hence the name. +The inverse Fourier transform is thereby defined as: + +$$ + z_j = \frac{1}{n} + \sum_{k=0}^{n-1} x_k \exp \left( \frac{2 \pi i j k}{n} \right) +$$ + +In GSL, `gsl_fft_complex_forward()` and `gsl_fft_complex_inverse()` are +functions which allow to compute the foreward and inverse transform, +respectively. + +In this special case, the sequence which must be transformed is made of real +numbers, but the Fourier transform is not real: it is a complex sequence wich +satisfies: + +$$ + z_k = z^*_{n-k} +$$ + +Where $z^*$ is the conjugate of $z$. A sequence with this symmetry is called +'half-complex'. This structure requires a particular storage layouts for the +forward transform (from real to half-complex) and inverse transform (from +half-complex to real). As a consequence, the routines are divided into two sets: +`gsl_fft_real` and `gsl_fft_halfcomplex`. The symmetry of the half-complex +sequence implies that only half of the complex numbers in the output need to be +stored. This works for all lengths: when the length is even, the middle value +where $k = n/2$ is real. Thus, only $n$ real numbers are required to store the +half-complex sequence (real and imaginary parts). +The ratio between the kernel and the sigmal transformations where therefore +made between complex numbers. +If the time-step of the DFT is $\Delta$, then the frequency-domain includes +both positive and negative frequencies, ranging from $-1 / (2 \Delta)$ to +$+1 / (2 \Delta$). The positive frequencies are stored from the beginning of +the array up to the middle, and the negative frequencies are stored backwards +from the end of the array. \textcolor{red}{Useful?} + + +--- + +
+![Convolved. $\sigma = 0.05 \Delta \theta$](images/noise-0.05.pdf){width=7cm} +![Deconvolved. $\sigma = 0.05 \Delta \theta$](images/deco-0.05.pdf){width=7cm} + +![Convolved. $\sigma = 0.1 \Delta \theta$](images/noise-0.1.pdf){width=7cm} +![Deconvolved. $\sigma = 0.1 \Delta \theta$](images/deco-0.1.pdf){width=7cm} + +![Convolved. $\sigma = 0.5 \Delta \theta$](images/noise-0.5.pdf){width=7cm} +![Deconvolved. $\sigma = 0.5 \Delta \theta$](images/deco-0.5.pdf){width=7cm} + +![Convolved. $\sigma = 1 \Delta \theta$](images/noise-1.pdf){width=7cm} +![Deconvolved. $\sigma = 1 \Delta \theta$](images/deco-1.pdf){width=7cm} + +Signal convolved with kernel on the left and deconvolved signal on the right. +Increasing values of $\sigma$ from top to bottom. $\Delta \theta$ is the bin +width. +
+ +As can be seen from the plots on the left in @fig:convolved, increasig the +value of $\sigma$ implies a stronger smoothing of the curve, while the +deconvolution process seems not to be affected by $\sigma$ amplitude changes: +it always gives the same outcome. + +It was also implemented the possibility to add a Poisson noise to the +distribution to check weather the deconvolution is affected or not by this +kind of noise.