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

ex-6: went on writing FFT, convolution concluded

Browse Source
This commit is contained in:
Giù Marcer 2020-03-23 22:49:47 +01:00 committed by rnhmjoj
parent 0a31d3b39b
commit 0a8351ea68
1 changed files with 80 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

107
notes/sections/6.md
View File

@ -196,27 +196,26 @@ one.
## 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)$:
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.
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:
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.
$$
\hat{F}[s*k] = \hat{F}[s] \cdot \hat{F}[k] \thus
\hat{F} [s] = \frac{\hat{F}[s*k]}{\hat{F}[k]}
$$
FFT are efficient algorithms for calculating the DFT. Namely, given a set of $n$
values {$z_i$}, each one is transformed into:
Being the histogram a discrete set of data, the Discrete Fourier Transform (DFT)
was emploied. In particular, 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)
@ -239,31 +238,85 @@ 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:
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
floating point numbers. In a packed array the real and imaginary parts of
each complex number are placed in alternate neighboring elements.
In this special case, the sequence of values 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
where $z^*$ is the conjugate of $z$. A sequence with this symmetry is called
'half-complex'. This structure requires 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?}
is real. Thus, only $n$ real numbers are required to store the half-complex
sequence (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 GSL functions aforementioned
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}
\hypertarget{fig:reorder}{%
\centering
\begin{tikzpicture}
\definecolor{cyclamen}{RGB}{146, 24, 43}
% standard histogram
\begin{scope}[shift={(7,0)}]
\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.5,0) rectangle (2.0,1.2);
\draw [thick, cyclamen, fill=cyclamen!25!white] (2.0,0) rectangle (2.5,1.4);
\draw [thick, cyclamen, fill=cyclamen!25!white] (2.5,0) rectangle (3.0,1.4);
\draw [thick, cyclamen, fill=cyclamen!25!white] (3.0,0) rectangle (3.5,1.2);
\draw [thick, cyclamen, fill=cyclamen!25!white] (3.5,0) rectangle (4.0,0.6);
\draw [thick, cyclamen] (4.5,0) -- (4.5,0.2);
\draw [thick, ->] (0,0) -- (5,0);
\draw [thick, ->] (2.5,0) -- (2.5,2);
\end{scope}
% shifted histogram
\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.5,0) rectangle (2.0,0.6);
\draw [thick, cyclamen, fill=cyclamen!25!white] (3.0,0) rectangle (3.5,0.6);
\draw [thick, cyclamen, fill=cyclamen!25!white] (3.5,0) rectangle (4.0,1.2);
\draw [thick, cyclamen, fill=cyclamen!25!white] (4.0,0) rectangle (4.5,1.4);
\draw [thick, ->] (0,0) -- (5,0);
\draw [thick, ->] (2.5,0) -- (2.5,2);
\end{tikzpicture}
\caption{On the left, an example of the DFT as it is given by the gsl function
and the same dataset, on the right, with the rearranged "intuitive"
order of the sequence.}\label{fig:reorder}
}
\end{figure}
When $\hat{F}[s*k]$ and $\hat{F}[k]$ are computed, their normal format must be
restored in order to 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 application.
GSL provides the function `gsl_fft_halfcomplex_unpack` which passes the vectors
from half-complex format to standard complex format. The inverse procedure,
required to compute the inverse transformation of $\hat{F}[s]$, which is not
provided by GSL, was implemented in the code.
The fact that the gaussian kernel is centerd in the middle of the vector and
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
are cut away in order to restore the original number of bins $n$.
---