From 0a8351ea682e7aa7bf100c1807277e4d1f43e6da Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gi=C3=B9=20Marcer?= Date: Mon, 23 Mar 2020 22:49:47 +0100 Subject: [PATCH] ex-6: went on writing FFT, convolution concluded --- notes/sections/6.md | 107 +++++++++++++++++++++++++++++++++----------- 1 file changed, 80 insertions(+), 27 deletions(-) diff --git a/notes/sections/6.md b/notes/sections/6.md index f14d7a5..a57b461 100644 --- a/notes/sections/6.md +++ b/notes/sections/6.md @@ -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$. ---