Giu/analistica
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

ex-6: Written about the RL algorithm and generated some plots

Browse Source
This commit is contained in:
Giù Marcer 2020-05-20 16:28:19 +02:00 committed by rnhmjoj
parent 7b8c61e3e5
commit 16c2511289
7 changed files with 169 additions and 169 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
notes/docs/bibliography.bib
View File

@ -118,7 +118,6 @@
publisher={Elsevier} publisher={Elsevier}
} }
@article{ridder17, @article{ridder17,
title={Variance reduction}, title={Variance reduction},
author={Ridder, AAN and Botev, ZI}, author={Ridder, AAN and Botev, ZI},
@ -138,3 +137,20 @@
year={1978}, year={1978},
publisher={Elsevier} publisher={Elsevier}
} }
@book{hecht02,
title={Optics},
year={2002},
publisher={Pearson},
author={Eugene Hecht}
}
@article{lucy74,
title={An iterative technique for the rectification of observed
distributions},
author={Lucy, Leon B},
journal={The astronomical journal},
volume={79},
pages={745},
year={1974}
}

BIN
notes/images/nome.pdf Normal file
View File

Binary file not shown.

BIN
notes/images/original.pdf Normal file
View File

Binary file not shown.

BIN
notes/images/prova.pdf Normal file
View File

Binary file not shown.

BIN
notes/images/smoothed.pdf Normal file
View File

Binary file not shown.

317
notes/sections/6.md
View File

@ -2,10 +2,10 @@
## Generating points according to Fraunhöfer diffraction ## Generating points according to Fraunhöfer diffraction
The diffraction of a plane wave through a round slit must be simulated by The diffraction of a plane wave through a round slit is to be simulated by
generating $N =$ 50'000 points according to the intensity distribution generating $N =$ 50'000 points according to the intensity distribution
$I(\theta)$ on a screen at a great distance $L$ from the slit itself: $I(\theta)$ [@hecht02] on a screen at a great distance $L$ from the slit itself
(see @fig:slit):
$$ $$
I(\theta) = \frac{E^2}{2} \left( \frac{2 \pi a^2 \cos{\theta}}{L} I(\theta) = \frac{E^2}{2} \left( \frac{2 \pi a^2 \cos{\theta}}{L}
\frac{J_1(x)}{x} \right)^2 \with x = k a \sin{\theta} \frac{J_1(x)}{x} \right)^2 \with x = k a \sin{\theta}
@ -16,7 +16,7 @@ where:
- $E$ is the electric field amplitude, default set $E = \SI{1e4}{V/m}$; - $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}$; - $a$ is the radius of the slit aperture, default set $a = \SI{0.01}{m}$;
- $\theta$ is the angle specified in @fig:slit; - $\theta$ is the angle specified in @fig:slit;
- $J_1$ is a Bessel function of first kind; - $J_1$ is the Bessel function of first kind;
- $k$ is the wavenumber, default set $k = \SI{1e-4}{m^{-1}}$; - $k$ is the wavenumber, default set $k = \SI{1e-4}{m^{-1}}$;
- $L$ default set $L = \SI{1}{m}$. - $L$ default set $L = \SI{1}{m}$.
@ -51,9 +51,8 @@ where:
\end{figure} \end{figure}
Once again, the *try and catch* method described in @sec:3 was implemented and Once again, the *try and catch* method described in @sec:3 was implemented and
the same procedure about the generation of $\theta$ was employed. This time, the same procedure about the generation of $\theta$ was applied. This time,
though, $\theta$ must be evenly distributed on half sphere: though, $\theta$ must be evenly distributed on half sphere, hence:
\begin{align*} \begin{align*}
\frac{d^2 P}{d\omega^2} = const = \frac{1}{2 \pi} \frac{d^2 P}{d\omega^2} = const = \frac{1}{2 \pi}
&\thus d^2 P = \frac{1}{2 \pi} d\omega^2 = &\thus d^2 P = \frac{1}{2 \pi} d\omega^2 =
@ -73,7 +72,6 @@ though, $\theta$ must be evenly distributed on half sphere:
If $\theta$ is chosen to grew together with $x$, then the absolute value can be If $\theta$ is chosen to grew together with $x$, then the absolute value can be
omitted: omitted:
\begin{align*} \begin{align*}
\frac{d\theta}{dx} = \frac{1}{\sin{\theta}} \frac{d\theta}{dx} = \frac{1}{\sin{\theta}}
&\thus d\theta \sin(\theta) = dx &\thus d\theta \sin(\theta) = dx
@ -85,44 +83,50 @@ omitted:
&\thus \theta = \text{acos} (1 -x) &\thus \theta = \text{acos} (1 -x)
\end{align*} \end{align*}
The sample was binned and stored in a histogram with a customizable number $n$ The so obtained sample was binned and stored in a histogram with a customizable
of bins default set $n = 150$. In @fig:original an example is shown. number $n$ of bins (default set $n = 150$) ranging from $\theta = 0$ to $\theta
= \pi/2$ bacause of the system symmetry. In @fig:original an example is shown.
![Example of an intensity histogram.](images/fraun-original.pdf){#fig:original} ![Example of intensity histogram.](images/original.pdf){#fig:original}
## Gaussian convolution {#sec:convolution} ## Gaussian convolution {#sec:convolution}
The sample must then be smeared with a Gaussian function with the aim to recover The sample has then to be smeared with a Gaussian function with the aim to
the original sample afterwards, implementing a deconvolution routine. recover the original sample afterwards, implementing a deconvolution routine.
For this purpose, a 'kernel' histogram with a odd number $m$ of bins and the For this purpose, a 'kernel' histogram with an even number $m$ of bins and the
same bin width of the previous one, but a smaller number of them ($m < n$), was same bin width of the previous one, but a smaller number of them ($m \sim 6\%
filled with $m$ points according to a Gaussian distribution with mean $\mu$, \, n$), was created according to a Gaussian distribution with mean $\mu$,
corresponding to the central bin, and variance $\sigma$. and variance $\sigma$. The reason why the kernel was set this way will be
descussed lately.
Then, the original histogram was convolved with the kernel in order to obtain Then, the original histogram was convolved with the kernel in order to obtain
the smeared signal. Some results in terms of various $\sigma$ are shown in the smeared signal. As an example, the result obtained for $\sigma = \Delta
[@fig:results1; @fig:results2; @fig:results3]. \theta$, where $\Delta \theta$ is the bin width, is shown in @fig:convolved.
As expected, the smeared signal looks smoother with respect to the original
one.
![Convolved signal.](images/smoothed.pdf){#fig:convolved}
The convolution was implemented as follow. Consider the definition of The convolution was implemented as follow. Consider the definition of
convolution of two functions $f(x)$ and $g(x)$: convolution of two functions $f(x)$ and $g(x)$:
$$ $$
f \otimes g (x) = \int \limits_{- \infty}^{+ \infty} dy f(y) g(x - y) f * g (x) = \int \limits_{- \infty}^{+ \infty} dy f(y) g(x - y)
$$ $$
Since a histogram is made of discrete values, a discrete convolution of the Since a histogram is made of discrete values, a discrete convolution of the
signal $s$ and the kernel $k$ must be computed. Hence, the procedure boils signal ($s$) and the kernel ($k$) must be computed. Hence, the procedure boils
down to an element wise product between $s$ and the reverse histogram of $k$ down to an element wise product between $s$ and the flipped histogram of $k$
for each relative position of the two histograms. Namely, if $c_i$ is the (from the last bin to the first one) for each relative position of the two
$i^{\text{th}}$ bin of the convoluted histogram: histograms. Namely, if $c_i$ is the $i^{\text{th}}$ bin of the convolved
histogram:
$$ $$
c_i = \sum_j k_j s_{i - j} c_i = \sum_j k_j s_{i - j} = \sum_{j'} k_{m - j'} s_{i - m + j'}
\with j' = m - j
$$ $$
where $j$ runs over the bins of the kernel. where $j$ runs over the bins of the kernel.
For a better understanding, see @fig:dot_conv. As can be seen, the third For a better understanding, see @fig:dot_conv: the third histogram turns out
histogram was obtained with $n + m - 1$ bins, a number greater than the initial with $n + m - 1$ bins, a number greater than the original one.
one.
\begin{figure} \begin{figure}
\hypertarget{fig:dot_conv}{% \hypertarget{fig:dot_conv}{%
@ -195,27 +199,32 @@ one.
## Unfolding with FFT ## Unfolding with FFT
Two different unfolding routines were implemented, one of which exploiting the Two different unfolding routines were implemented, one of which exploiting the
Fast Fourier Transform. This method is based on the property of the Fourier Fast Fourier Transform (FFT).
transform according to which, given two functions $f(x)$ and $g(x)$:
This method is based on the convolution theorem, according to which, given two
functions $f(x)$ and $g(x)$:
$$ $$
\hat{F}[f \otimes g] = \hat{F}[f] \cdot \hat{F}[g] \hat{F}[f * g] = \hat{F}[f] \cdot \hat{F}[g]
$$ $$
where $\hat{F}[\quad]$ stands for the Fourier transform of its argument. 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 smeared signal and the kernel, the ratio between their
transforms and the anti-transformation of the result:
$$
\hat{F}[s \otimes k] = \hat{F}[s] \cdot \hat{F}[k] \thus
\hat{F} [s] = \frac{\hat{F}[s \otimes k]}{\hat{F}[k]}
$$
Being the histogram a discrete set of data, the Discrete Fourier Transform (DFT) Being the histogram a discrete set of data, the Discrete Fourier Transform (DFT)
was emploied. In particular, the FFT are efficient algorithms for calculating was appied. When dealing with arrays of discrete values, the theorems still
the DFT. Given a set of $n$ values {$z_i$}, each one is transformed into: holds if the two arrays have the same lenght and a cyclical convolution is
applied. For this reason the kernel was 0-padded in order to make it the same
lenght of the original signal. Besides, the 0-padding allows to avoid unpleasant
side effects due to the cyclical convolution.
In order to accomplish this procedure, every histogram was transformed into a
vector. The implementation lies in the computation of the Fourier trasform of
the smeared signal and the kernel, the ratio between their transforms and the
anti-transformation of the result:
$$
\hat{F}[s * k] = \hat{F}[s] \cdot \hat{F}[k] \thus
\hat{F} [s] = \frac{\hat{F}[s * k]}{\hat{F}[k]}
$$
The FFT are efficient algorithms for calculating the DFT. 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) x_j = \sum_{k=0}^{n-1} z_k \exp \left( - \frac{2 \pi i j k}{n} \right)
$$ $$
@ -227,7 +236,6 @@ 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)$ \ldots n_m$, then the DFT can be computed in $O(n \sum n_i) < O(n^2)$
operations, hence the name. operations, hence the name.
The inverse Fourier transform is thereby defined as: The inverse Fourier transform is thereby defined as:
$$ $$
z_j = \frac{1}{n} z_j = \frac{1}{n}
\sum_{k=0}^{n-1} x_k \exp \left( \frac{2 \pi i j k}{n} \right) \sum_{k=0}^{n-1} x_k \exp \left( \frac{2 \pi i j k}{n} \right)
@ -235,20 +243,12 @@ $$
In GSL, `gsl_fft_complex_forward()` and `gsl_fft_complex_inverse()` are In GSL, `gsl_fft_complex_forward()` and `gsl_fft_complex_inverse()` are
functions which allow to compute the foreward and inverse transform, functions which allow to compute the foreward and inverse transform,
respectively. respectively.
In order to accomplish this procedure, every histogram was transformed into a
vector. The kernel vector was 0-padded and centred in the middle to make its
length the same as that of the signal, making it feasable to implement the
division between the entries of the vectors one by one.
The inputs and outputs for the complex FFT routines are packed arrays of The inputs and outputs for the complex FFT routines are packed arrays of
floating point numbers. In a packed array the real and imaginary parts of floating point numbers. In a packed array the real and imaginary parts of
each complex number are placed in alternate neighboring elements. each complex number are placed in alternate neighboring elements. In this
In this special case, the sequence of values which must be transformed is made special case, the sequence of values to be transformed is made of real numbers,
of real numbers, but the Fourier transform is not real: it is a complex sequence but the Fourier transform is not real: it is a complex sequence wich satisfies:
wich satisfies:
$$ $$
z_k = z^*_{n-k} z_k = z^*_{n-k}
$$ $$
@ -258,14 +258,10 @@ where $z^*$ is the conjugate of $z$. A sequence with this symmetry is called
forward transform (from real to half-complex) and inverse transform (from 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: 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 `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 sequence requires only half of the complex numbers in the output to be stored.
stored. This works for all lengths: when the length is even, the middle value This works for all lengths: when the length is even, the middle value is real.
is real. Thus, only $n$ real numbers are required to store the half-complex Thus, only $n$ real numbers are required to store the half-complex sequence
sequence (half for the real part and half for the imaginary). (half for the real part and half for the imaginary).
If the bin width is $\Delta \theta$, then the DFT domain ranges from $-1 / (2
\Delta \theta)$ to $+1 / (2 \Delta \theta$). The aforementioned GSL functions
store the positive values from the beginning of the array up to the middle and
the negative backwards from the end of the array (see @fig:reorder).
\begin{figure} \begin{figure}
\hypertarget{fig:reorder}{% \hypertarget{fig:reorder}{%
@ -273,7 +269,6 @@ the negative backwards from the end of the array (see @fig:reorder).
\begin{tikzpicture} \begin{tikzpicture}
\definecolor{cyclamen}{RGB}{146, 24, 43} \definecolor{cyclamen}{RGB}{146, 24, 43}
% standard histogram % standard histogram
\begin{scope}[shift={(7,0)}]
\draw [thick, cyclamen] (0.5,0) -- (0.5,0.2); \draw [thick, cyclamen] (0.5,0) -- (0.5,0.2);
\draw [thick, cyclamen, fill=cyclamen!25!white] (1.0,0) rectangle (1.5,0.6); \draw [thick, cyclamen, fill=cyclamen!25!white] (1.0,0) rectangle (1.5,0.6);
\draw [thick, cyclamen, fill=cyclamen!25!white] (1.5,0) rectangle (2.0,1.2); \draw [thick, cyclamen, fill=cyclamen!25!white] (1.5,0) rectangle (2.0,1.2);
@ -284,8 +279,8 @@ the negative backwards from the end of the array (see @fig:reorder).
\draw [thick, cyclamen] (4.5,0) -- (4.5,0.2); \draw [thick, cyclamen] (4.5,0) -- (4.5,0.2);
\draw [thick, ->] (0,0) -- (5,0); \draw [thick, ->] (0,0) -- (5,0);
\draw [thick, ->] (2.5,0) -- (2.5,2); \draw [thick, ->] (2.5,0) -- (2.5,2);
\end{scope}
% shifted histogram % shifted histogram
\begin{scope}[shift={(7,0)}]
\draw [thick, cyclamen, fill=cyclamen!25!white] (0.5,0) rectangle (1.0,1.4); \draw [thick, cyclamen, fill=cyclamen!25!white] (0.5,0) rectangle (1.0,1.4);
\draw [thick, cyclamen, fill=cyclamen!25!white] (1.0,0) rectangle (1.5,1.2); \draw [thick, cyclamen, fill=cyclamen!25!white] (1.0,0) rectangle (1.5,1.2);
\draw [thick, cyclamen, fill=cyclamen!25!white] (1.5,0) rectangle (2.0,0.6); \draw [thick, cyclamen, fill=cyclamen!25!white] (1.5,0) rectangle (2.0,0.6);
@ -294,94 +289,126 @@ the negative backwards from the end of the array (see @fig:reorder).
\draw [thick, cyclamen, fill=cyclamen!25!white] (4.0,0) rectangle (4.5,1.4); \draw [thick, cyclamen, fill=cyclamen!25!white] (4.0,0) rectangle (4.5,1.4);
\draw [thick, ->] (0,0) -- (5,0); \draw [thick, ->] (0,0) -- (5,0);
\draw [thick, ->] (2.5,0) -- (2.5,2); \draw [thick, ->] (2.5,0) -- (2.5,2);
\end{scope}
\end{tikzpicture} \end{tikzpicture}
\caption{On the left, an example of the DFT as it is given by the gsl function \caption{The histogram on the right shows how the real numbers histogram on the
and the same dataset, on the right, with the rearranged "intuitive" left is handled by the dedicated GSL functions`.}\label{fig:reorder}
order of the sequence.}\label{fig:reorder}
} }
\end{figure} \end{figure}
When $\hat{F}[s \otimes k]$ and $\hat{F}[k]$ are computed, their normal format If the bin width is $\Delta \theta$, then the DFT domain ranges from $-1 / (2
must be restored in order to use them as standard complex numbers and compute \Delta \theta)$ to $+1 / (2 \Delta \theta$). As regards the real values, the
the ratio between them. Then, the result must return in the half-complex format aforementioned GSL functions store the positive values from the beginning of
for the inverse DFT application. the array up to the middle and the negative backwards from the end of the array
GSL provides the function `gsl_fft_halfcomplex_unpack()` which passes the (see @fig:reorder).
vectors from half-complex format to standard complex format. The inverse Whilst do not matters if the convolved histogram has positive or negative
procedure, required to compute the inverse transformation of $\hat{F}[s]$, which values, the kernel must be centered in zero in order to compute a correct
is not provided by GSL, was implemented in the code. convolution. This requires the kernel to be made of an ever number of bins
The fact that the gaussian kernel is centerd in the middle of the vector and in order to be possible to cut it into two same-lenght halves.
not in the $\text{zero}^{th}$ bin causes the final result to be shifted of half
the leght of the vector the same as it was produced by a DFT. This makes it
necessary to rearrange the two halfs of the final result.
At the end, the external bins which exceed with respect to the original signal When $\hat{F}[s * k]$ and $\hat{F}[k]$ are computed, they are given in the
are cut away in order to restore the original number of bins $n$. Results are half-complex GSL format and their normal format must be restored in order to
shown in [@fig:results1; @fig:results2; @fig:results3]. use them as standard complex numbers and compute the ratio between them. Then,
the result must return in the half-complex format for the inverse DFT
computation. GSL provides the function `gsl_fft_halfcomplex_unpack()` which
convert the vectors from half-complex format to standard complex format but the
inverse procedure is not provided by GSL and was hence implemented in the
code.
At the end, the external bins which exceed the original signal size are cut
away in order to restore the original number of bins $n$. Results will be
discussed in @sec:conv_Results.
## Unfolding with Richardson-Lucy ## Unfolding with Richardson-Lucy
The Richardson–Lucy (RL) deconvolution is an iterative procedure tipically used The Richardson–Lucy (RL) deconvolution is an iterative procedure tipically used
for recovering an image that has been blurred by a known point spread function. for recovering an image that has been blurred by a known 'point spread
function'.
It is based on the fact that an ideal point source does not appear as a point
but is spread out into the so-called point spread function, thus the observed
image can be represented in terms of a transition matrix
$P$ operating on an underlying image:
Consider the problem of estimating the frequeny distribution $f(\xi)$ of a
variable $\xi$ when the available measure is a sample {$x_i$} of points
dronwn not by $f(x)$ but by an other function $\phi(x)$ such that:
$$ $$
d_i = \sum_{j} u_j \, P_{i, j} \phi(x) = \int d\xi \, f(\xi) P(x | \xi)
$$ $$
where $u_j$ is the intensity of the underlying image at pixel $j$ and $d_i$ is where $P(x | \xi) \, d\xi$ is the probability (presumed known) that $x$ falls
the detected intensity at pixel $i$. Hence, the matrix describes the portion of in the interval $(x, x + dx)$ when $\xi = \xi$. An example of this problem is
signal from the source pixel $j$ that is detected in pixel $i$. precisely that of correcting an observed distribution $\phi(x)$ for the effect
In one dimension, the transfer function can be expressed in terms of the of observational errors, which are represented by the function $P (x | \xi)$,
distance between the source pixel $j$ and the observed $i$: called point spread function.
Let $Q(\xi | x) d\xi$ be the probability that $\xi$ comes from the interval
$(\xi, \xi + d\xi)$ when the measured quantity is $x = x$. The probability that
both $x \in (x, x + dx)$ and $(\xi, \xi + d\xi)$ is therefore given by $\phi(x)
dx \cdot Q(\xi | x) d\xi$ which is identical to $f(\xi) d\xi \cdot P(x | \xi)
dx$, hence:
$$ $$
P_{i, j} = \widetilde{P}(i-j) = P_{i - j} \phi(x) dx \cdot Q(\xi | x) d\xi = f(\xi) d\xi \cdot P(x | \xi) dx
\thus Q(\xi | x) = \frac{f(\xi) \cdot P(x | \xi)}{\phi(x)}
$$ $$
In order to estimate $u_j$ given {$d_i$} and $\widetilde{P}$, the following
iterative procedure can be applied for the estimate $\hat{u}^t_j$ of $u_j$,
where $t$ stands for the iteration number. The $t^{\text{th}}$ step is updated
as follows:
$$ $$
\hat{u}^{t+1}_j = \hat{u}^t_j \sum_i \frac{d_i}{c_i} \, P_{i - j} \thus Q(\xi | x) = \frac{f(\xi) \cdot P(x | \xi)}
\with c_i = \sum_j \hat{u}^t_j \, P_{i - j} {\int d\xi \, f(\xi) P(x | \xi)}
$$ {#eq:first}
which is the Bayes theorem for conditional probability. From the normalization
of $P(x | \xi)$, it follows also that:
$$
f(\xi) = \int dx \, \phi(x) Q(\xi | x)
$$ {#eq:second}
Since $Q (\xi | x)$ depends on $f(\xi)$, @eq:second suggests a reiterative
procedure for generating estimates of $f(\xi)$. With a guess for $f(\xi)$ and
a known $P(x | \xi)$, @eq:first can be used to calculate and estimate for
$Q (\xi | x)$. Then, taking the hint provided by @eq:second, an improved
estimate for $f (\xi)$ is generated, using the observed sample {$x_i$} to give
an approximation for $\phi$.
Thus, if $f^t$ is the $t^{\text{th}}$ estimate, the $t^{\text{th + 1}}$ is:
$$
f^{t + 1}(\xi) = \int dx \, \phi(x) Q^t(\xi | x)
\with
Q^t(\xi | x) = \frac{f^t(\xi) \cdot P(x | \xi)}
{\int d\xi \, f^t(\xi) P(x | \xi)}
$$ $$
where $c_i$ is thereby an estimation of the blurred signal obtained with the from which:
previous estimation of the clean signal.
It has been shown empirically that if this iteration converges, it converges to
the maximum likelihood solution for $u_j$. Writing it in terms of convolution,
it becomes:
$$ $$
\hat{u}^{t+1} = \hat {u}^{t} \cdot \left( \frac{d}{{\hat{u}^{t}} \otimes P} f^{t + 1}(\xi) = f^t(\xi)
\otimes P^{\star} \right) \int dx \, \frac{\phi(x)}{\int d\xi \, f^t(\xi) P(x | \xi)}
P(x | \xi)
$$ {#eq:solution}
If the spread function $P(x | \xi)$ follows a normal distribution with variance
$\sigma$, namely:
$$
P(x | \xi) = \frac{1}{\sqrt{2 \pi} \sigma}
\exp \left( - \frac{(x - \xi)^2}{2 \sigma^2} \right)
$$ $$
where the division and multiplication are element wise, and then, @eq:solution can be rewritten in terms of convolutions:
$P^{\star}$ is the flipped point spread function. $$
f^{t + 1} = f^{t}\left( \frac{\phi}{{f^{t}} * P} * P^{\star} \right)
$$
where $P^{\star}$ is the flipped point spread function [@lucy74].
In this special case, the Gaussian kernel stands for the point spread function
and, dealing with discrete values, the division and multiplication are element
wise and the convolution is to be carried out as described in @sec:convolution.
When implemented, this method results in an easy step-wise routine: When implemented, this method results in an easy step-wise routine:
- create a flipped copy of the kernel; - create a flipped copy of the kernel;
- choose a zero-order estimate for {$c_i$}; - choose a zero-order estimate for {$f(\xi)$};
- compute the convolutions with the method described in @sec:convolution, the - compute the convolutions, the product and the division at each step;
product and the division at each step;
- proceed until a given number of reiterations is achieved. - proceed until a given number of reiterations is achieved.
In this case, the zero-order was set $c_i = 0.5 \, \forall i$ and it was In this case, the zero-order was set $f(\xi) = 0.5 \, \forall \, \xi$.
empirically shown that the better result is given with a number of three steps,
otherwise it starts returnig fanciful histograms. Results are shown in
[@fig:results1; @fig:results2; @fig:results3].
## Results comparison ## Results comparison {#sec:conv_Results}
In [@fig:results1; @fig:results2; @fig:results3] the results obtained for three In [@fig:results1; @fig:results2; @fig:results3] the results obtained for three
different $\sigma$s are shown. The tested values are $\Delta \theta$, $0.5 \, different $\sigma$s are shown. The tested values are $\Delta \theta$, $0.5 \,
@ -399,46 +426,14 @@ In the real world, it is unpratical, since signals are inevitably blurred by
noise. noise.
The same can't be said about the RL deconvolution, which, on the other hand, The same can't be said about the RL deconvolution, which, on the other hand,
looks heavily influenced by the variance magnitude: the greater $\sigma$, the looks heavily influenced by the variance magnitude: the greater $\sigma$, the
worse the deconvoluted result. In fact, given the same number of steps, the worse the deconvolved result. In fact, given the same number of steps, the
deconvolved signal is always the same 'distance' far form the convolved one: deconvolved signal is always the same 'distance' far form the convolved one:
if it very smooth, the deconvolved signal is very smooth too and if the if it very smooth, the deconvolved signal is very smooth too and if the
convolved is less smooth, it is less smooth too. convolved is less smooth, it is less smooth too.
The original signal is shown below for convenience. The original signal is shown below for convenience.
![Example of an intensity histogram.](images/fraun-original.pdf)
<div id="fig:results1">
![Convolved signal.](images/fraun-conv-0.05.pdf){width=12cm}
![Deconvolved signal with FFT.](images/fraun-fft-0.05.pdf){width=12cm}
![Deconvolved signal with RL.](images/fraun-rl-0.05.pdf){width=12cm}
Results for $\sigma = 0.05 \Delta \theta$, where $\Delta \theta$ is the bin
width.
</div>
<div id="fig:results2">
![Convolved signal.](images/fraun-conv-0.5.pdf){width=12cm}
![Deconvolved signal with FFT.](images/fraun-fft-0.5.pdf){width=12cm}
![Deconvolved signal with RL.](images/fraun-rl-0.5.pdf){width=12cm}
Results for $\sigma = 0.5 \Delta \theta$, where $\Delta \theta$ is the bin
width.
</div>
<div id="fig:results3">
![Convolved signal.](images/fraun-conv-1.pdf){width=12cm}
![Deconvolved signal with FFT.](images/fraun-fft-1.pdf){width=12cm}
![Deconvolved signal with RL.](images/fraun-rl-1.pdf){width=12cm}
Results for $\sigma = \Delta \theta$, where $\Delta \theta$ is the bin width.
</div>
It was also implemented the possibility to add a Poisson noise to the It was also implemented the possibility to add a Poisson noise to the
convolved histogram to check weather the deconvolution is affected or not by convolved histogram to check weather the deconvolution is affected or not by
@ -451,14 +446,6 @@ everywhere on the curve and it is particularly evident on the edges, where the
expected data are very small. On the other hand, the Richardson-Lucy routine is expected data are very small. On the other hand, the Richardson-Lucy routine is
less affected by this further complication. less affected by this further complication.
<div id="fig:poisson">
![Deconvolved signal with FFT.](images/fraun-noise-fft.pdf){width=12cm}
![Deconvolved signal withh RL.](images/fraun-noise-rl.pdf){width=12cm}
Results for $\sigma = \Delta \theta$, with Poisson noise.
</div>
In order to quantify the similarity of a deconvolution outcome with the original In order to quantify the similarity of a deconvolution outcome with the original
signal, a null hypotesis test was made up. signal, a null hypotesis test was made up.
Likewise in @sec:Landau, the original sample was treated as a population from Likewise in @sec:Landau, the original sample was treated as a population from
@ -519,7 +506,7 @@ $$
{\sum_{i = 1}^m \sum_{j=1}^n f_{ij}} {\sum_{i = 1}^m \sum_{j=1}^n f_{ij}}
$$ $$
In this case, where the EMD must be applied to two same-lenght histograms, the In this case, where the EMD has to be applied to two same-lenght histograms, the
procedure simplifies a lot. By representing both histograms with two vectors $u$ procedure simplifies a lot. By representing both histograms with two vectors $u$
and $v$, the equation above boils down to [@ramdas17]: and $v$, the equation above boils down to [@ramdas17]:

3
notes/todo
View File

@ -1,3 +0,0 @@
- cambiare simbolo convoluzione
- aggiungere citazioni e referenze
- leggere l'articolo di Lucy