path: root/doc/thesis/chapters/implementation.tex

% vim: set ts=2 sw=2 noet spell:

\chapter{Implementation}

\section{Overview}

First of all the tools which were used in this theses are interdused; software and hardware. Then the sender and receiver chain are explained in detail and how they are implemented. This both chains are the same for the simulation itself and the hardware implementation. Also the channel in details are explained which is inportent for the simulation to illustrate the multipath fading effect. After some measurements are tried, with the help of the bit error rate, to compare it with the simulation.

\section{Software Stack}

\subsection{GNU Radio}
For both the signal processing and the simulations the GNU Radio (GR) toolkit was chosen, as it already had drivers for the USRP hardware. GR is an open-source free software framework that can be used to build signal processing chains and software defined radios (SDR). GR is composed of two parts: a C\texttt{++} library with Python bindings to write signal processing code, and GNU Radio Companion (GRC), a graphical user interface to more easily construct signal processing chains by representing signal processing algorithms as ``blocks'' that are chained together with arrows, essentially drawing a diagram called ``flow graph''. An example of a flow graph is shown in \figref{fig:flowgraph}.

Internally GR works by keeping multiple memory buffers of samples, that are passed as pointers to the signal processing algorithms' ``work functions''. When the signal processing is complete, the output buffer of one block is given to the next block as input according to how they were connected in the flow graph. The structure of a block is shown in the Python listing \ref{lst:gr-block-py}. To improve performance GR creates a thread for each work function to parallelize the workload of the concurrently running signal processing blocks. For more details see the GNU Radio Wiki and User Manual in \cite{GRWiki}.

{\newcommand{\placeholder}[1]{\textit{\(\langle\)\,\textrm{#1}\,\(\rangle\)}}
\begin{lstlisting}[
	language = python, escapechar = {`}, float,
	caption = {A minimal GNU Radio block in Python.},
	captionpos = b, label = {lst:gr-block-py}
]
class myblock(gr.sync_block):
	# there are also other types of blocks such as interpolators 
	# (more outputs that inputs), decimators (more inputs than 
	# outputs) sync blocks have a 1-to-1 input to output ratio
	def __init__(self, `\placeholder{parameters}`):
		gr.sync_block.__init__(self, name="My Block",
			# this block as one input port and one output port
			# with samples that are 64 bit complex numbers
			in_sig=[np.complex64], out_sig=[np.complex64]
		)

	def work(self, inputs, outputs):
		# signal processing goes here
		# inputs and outputs are k x n arrays, where each
		# of the k rows is a port that contains n samples
		return `\placeholder{N of inputs that were processed}`
\end{lstlisting}}

\subsection{Dear PyGUI}

To construct a graphical interface for a demonstration platform the Dear IMGUI (immediate graphical user interface) library was chosen, mainly for its ease of use, wide range of techincal capabilites and high refresh rate. Dear PyGUI (DPG) are the Python bindings for the Dear IMGUI library.

The DPG GUI communicates with the GR flow graphs using the IP/UDP protocol. This decision to separate the project into two parts that communicate over the IP network was made because it is not very easy to extend the graphical interface of GRC without interfering with the sophisticated multi-threaded architecture of GR. Furthermore, this allows to have multiple correctly configure flow graph on disk and to choose which one to run and display on the graphical interface, instead of having a single flow graph whose parameters need to be changed each time.

Also as a side effect, in theory this setup allows to have one computer running the graphical interface, and another remote machine running just the flow graph.

\section{Hardware}

\begin{table}[b]
	%To DO sepzifikationen ampssen / genauer? https://www.ettus.com/wp-content/uploads/2019/01/b200-b210_spec_sheet.pdf
	% https://kb.ettus.com/B200/B210/B200mini/B205mini#FAQ
	\centering
	\begin{tabular}{ll}
		\toprule
		Dimensions               & \(9.7 \times 15.5 \times 1.5\) cm \\
		Ports                    & 2 TX, 2 RX, Half or Full Duplex   \\
		RF frequencies           & 70MHz to 6GHz                     \\
		Bandwidth                & 200kHz -- 56MHz                   \\
		External reference input & 10 MHz                            \\
		\bottomrule
	\end{tabular}
	\caption{USRP B210 specifications.\label{tab:usrp-specs}}
\end{table}

As receivers and transmitter devices for the SDR setup two USRP B210 devices from Ettus Research were used. Some technical specifications are shown in \tabref{tab:usrp-specs}. GR provides off the shelf blocks that interface with the official API provided from Ettus Research.

\section{Sender chain}
\subsection{Data frame}

\begin{figure}
\centering
\input{figures/tikz/packet-frame}
	\caption{
		Structure of framed data packets used in the implementation.
		\label{fig:dataframe}
	}
\end{figure}

To compute the empirical \emph{bit error rate} (BER) of the setup, the data has to be framed by the sender and the bitstream synchronized on the receiver side. The structure of a data packet used in the implementation is shown in \figref{fig:dataframe}. A frame begins with an user specified \(k\)-byte preamble, that in the current implementation serves as synchronization pattern. Another use case for the preamble sequence could be to introduce channel estimation pilot symbols. Following the preamble are 4 bytes encoded using a (31, 26) Hamming code (plus 1 padding bit), that contain metadata about the packet, namely payload ID and payload length. Because the payload length in bytes is encoded in 21 bits, the maximum payload size is 2 MiB, which together with 32 possible unique IDs gives a maximum data transfer with unique frame headers of 64 MiB. These constraints are a result of decisions made to keep the implementation simple.


\subsection{Modulation}

GR provides a constellation modulator block, that already implements several standard constellations (QPSK and 16ary-QAM being of interest for us). The block also already integrates a root raised cosine filter, whose phase bandwidth (roll-off factor) can be given as parameter.

\section{Receiver chain}

\subsection{Envelope detector}

What is here referred to as envelope detector has the purpose of synchronizing the symbols and equalizing the input signal amplitude. This is accomplished in GRC using two blocks: a polyphase clock sync and a CMA equalizer. The input signal for the envelope detector has 4 samples per symbol, while the output has only one sample per symbol.

\paragraph{Polyphase Clock Sync}
%% To Do : nochmals anschauen ob dieese erklärung verständlich ist und richtig interpretiert wurde.
With the the polyphase clock sync the symbols can be synchronized by preforming a time synchronization with the help of multiple filterbanks. For that the derivative of the filtered signal should be minimized which turns to a better SNR. 
% This works with the help of two filterbanks, one of them contains the filters of the signal adapted to the pulse shaping with several phases. The other contains its derivative. So in the time domain it has a sinc shape, for the output Signal the sinc peak should be on a sample, with the fact that sinc(0) = 1 and sinc(0)' = 0 an error signal can be generated which tells how far away from the peak it is. This error Signal should be zero this is possible with the help of a loop second order whish constants the number of the filterbank and the rate. This rate is generated because of the clock difference between the transmitter and receiver to synchronized the receiver the filter goes through the phases. For the output one sample per symbol is enough.

\paragraph{Equalizer}

\skelpar[2]{CMA equalizer.}

\subsection{Frame synchronization and phase correction} \label{sec:phasecorr}

Once the envelope's clock is synchronized in the processing chain the data stream has one sample per symbol. At this point it is necessary to find where each data frame starts or end in order to correctly decode their payloads. For such purpose a special sequence called \emph{access code} is put in front of each frame. Access codes are sequences of samples that are carefully constructed to have an autocorrelation with a high peak at zero, and that rapidly decreases for increasing shifts. In other words, the autocorrelation of an access code high only when the sequence is perfectly aligned. Thus by cross correlating an envelope signal \(r(t)\), that periodically contains an access code \(a(t)\) with the access code itself, and looking for peaks in the result, it is possible to determine where each frame begin. Furthermore by analyzing the values of the peaks it is possible to extract information about the phase and frequency offsets.

To understand how correlation peaks allow for fine phase correction, recall that mathematically the cross correlation (denoted here by \(\star\)) of two complex valued signals is
\begin{equation}
	R_{ra}
	= (r \star a)(t)
	= \int_\mathbb{R} r(\tau) a^*(\tau - t) \,d\tau
	= r(t) * a^*(-t),
\end{equation}
which is equivalent to a convolution, with the left term being time-reversed complex conjugated \cite{Gallager}. This last property is especially useful because it makes possible to implement cross correlation using FIR filters. Some interesting properties of the cross correlation are that correlation with itself (autocorrelation) at \(t = 0\) is
\begin{equation}
	R_{aa} = (a \star a)(0)
	= \int_\mathbb{R} a(\tau) a^*(\tau - 0) \,d\tau
	= \int_\mathbb{R} |a(\tau)|^2 \,d\tau \in \mathbb{R},
\end{equation}
which is a real number. And more importantly the correlation with an out of phase copy \(a'(t) = a(t) e^{j\varphi}\) at 0 is
\begin{equation} \label{eqn:xc-oop-copy}
	% R_{a'a} =
	(a' \star a)(0) 
	= \int_\mathbb{R} a(\tau)e^{j\varphi}  a^*(\tau) \,d\tau
	= R_{aa} e^{j\varphi}.
\end{equation}
The relevant observation to make in \eqref{eqn:xc-oop-copy} that since \(R_{aa}\) is a real number, the phase of the cross correlation at \(t = 0\) is the phase of \(a'(t)\). This fact can be exploited to implement fine phase correction for the received envelope in relatively few steps as follows:
\begin{enumerate}
	\item Compute the cross correlation \(R_{ra}\) of the envelope \(r(t)\) with the access code \(a(t)\),
	\item Find the maximum value of \(\hat{R}_{ra} = \max R_{ra}(t)\) (correlation peak),
	\item Extract the phase offset \(\varphi = \arg \hat{R}_{ra}\),
	\item Remove the phase offset in the envelope by multiplying it with the complex conjugate of the offset, that is \(\hat{r}(t) = r(t) e^{-j\varphi}\).
\end{enumerate}

\subsection{Implementing fine phase and frequency correction}

To implement in GR what was discussed in \S\ref{sec:phasecorr} two blocks shown in \figref{fig:phasecorr-blocks} were used: a correlator estimator block, and a custom block. The former essentially implements the first 3 of the steps discussed at the end of \S\ref{sec:phasecorr}. The correlator estimator block is given a sequence of samples, and when the cross correlation between them and the input stream is higher than a certain threshold (90\% of the amplitude of a perfect autocorrelation), it produces a ``tag'' in the output stream, that contains the phase estimate.

Tags are GR's way of working with metadata that is attached to a sample. Internally tags are just polymorphic data structures containing a number indicating the absolute offset (in samples), and a pair of arbitrary values called ``key'' and ``value''. Tags are passed on from one block to the next like sample streams (unless the block specifies to do otherwise).

Thus the tagged stream is processed with a custom block, of which a simplified version of its work function shown in listing \ref{lst:phasecorr-work}. The custom block also implements fine frequency correction by linearly interpolating the phase estimates between each pair of tags (called chunk). Mathematically this can be rather trivially be formulated for a chunk of \(N\) samples with the
\begin{subequations}
	\begin{align}
		k\text{-th chunk digital frequency} \quad  & \omega_k = (\varphi_{k+1} - \varphi_k) / N, \text{ and the }\\
		k\text{-th chunk phase estimate} \quad & \phi_k(n) = \varphi_k - \omega_k n/N.
	\end{align}
\end{subequations}

% TODO: discuss performance of implementation

\begin{lstlisting}[
	texcl = true, language = python, escapechar = {`},
	float, captionpos = b, label = {lst:phasecorr-work},
	caption = {
		Simplified work function of fine phase correction block that corrects only samples `in the middle'. The version that is actually used does handle edge cases that have been removed here for readability. See also \figref{fig:phasecorr-chunks} for a graphical representation of the inputs and listing \ref{lst:phasecorr-blockphase} for the definition of the \texttt{block\_phase} function.
	},
]
def work(self, inputs, outputs):
	# alias for inputs of the first port
	inp = inputs[0]
	# read phase tags from stream
	is_phase = lambda tag: pmt.to_python(tag.key) == "phase_est"
	tags = filter(is_phase, self.get_tags_in_window(0, 0, len(inp)))
	# create a list of pairs \(((\varphi_0,\varphi_1), (\varphi_1, \varphi_2), \ldots, (\varphi_{k-1}, \varphi_k)))\)
	pairs = zip(tags, tags[1:])
	# compute phase correction between each pair of tags
	chunks = [self.block_phase(start, end) for (start, end) in pairs]
	# flatten array to get \(\phi(n)\) and compute the correction
	phases = np.concatenate(chunks)
	correction = np.exp(-1j * phases) 
	# write to the first output port
	left = tags[0].offset - self.nitems_written(0)
	right = tags[-1].offset - self.nitems_written(0)
	outputs[0][left:right] = inp * correction
	# return how many samples were processed
	return len(outputs[0])
\end{lstlisting}

\begin{figure}
	\centering
	%% TODO: move code into separate file
	\begin{tikzpicture}[
			blk/.style = {
				draw, rectangle, thick, black,
				minimum width = 15mm,
				minimum height = 3mm,
				outer sep = 1mm,
				pattern = vertical lines,
				pattern color = lightgray,
			},
		]

		\foreach \i in {0,1,...,4}{
			\coordinate (blkC\i) at (15mm*\i,0);
			\node[blk] (blk\i) at (blkC\i) {};
			\node[below] (phi\i) at (blk\i.south west) {\(\varphi_{\i}\)};
		}
		% last phase 
		\node[below] (phi5) at (blk4.south east) {\(\varphi_{5}\)};

		% first block
		\node[blk, minimum width = 10mm, xshift = 2.5mm, fill = red!30] (S) at (-15mm,0) {};
		\node[anchor = east] at (S.west) {Input};

		% last block
		\node[blk, minimum width = 9mm, xshift = -3mm, fill = blue!30] (E) at ($(blk4)+(15mm,0)$) {};

		% labels
		\draw[thick, latex-] (blk3.north) to[out = 90, in = 180] ++(5mm,6mm)
			node[right] {Chunk of \(N\) samples};

		\draw[thick, latex-] (blk0.north east) ++(-1mm,0) to[out = 90, in = 0] ++(-5mm,6mm)
			node[left] (tags) {Phase tags};
	\end{tikzpicture}
	\caption{
		Graphical representation of the input samples for the work function of the fine phase and frequency correction block (shown in listing \ref{lst:phasecorr-work}). Roughly every \(N\) samples there is a tag containing the information of the phase error (computed using the cross correlation peak). The white `chunks' of samples can be corrected using their respective left and right tag values. The samples in the red chunk need phase information from the previous block processing. The samples in the blue chunk need a phase information from the future, which is not attainable. Thus for the blue chunk the frequency estimate of the previous chunk is used.
		\label{fig:phasecorr-chunks}
	}
\end{figure}

\begin{lstlisting}[
	texcl = true, language = python, escapechar = {`},
	float, captionpos = b, label = {lst:phasecorr-blockphase},
	caption = {
		Block phase function referenced in listing \ref{lst:phasecorr-work}.
	},
]
def block_phase(self, start, end):
	# compute number of samples between tags
	nsamples = end.offset - start.offset
	# unpack pmt values into start and end phase
	sphase = pmt.to_python(start.value)
	ephase = pmt.to_python(end.value)
	# compute frequency offset between start and end
	phasediff = (ephase - sphase) % (2 * np.pi)
	freq = phasediff / nsamples
	# compute chunk values
	return sphase * np.ones(nsamples) + freq * np.arange(0, nsamples)
\end{lstlisting}


\section{Channel simulations}

Here its possible to add some AWGN noise in the channel line. Different parameters can be changed like the noise voltage, time or the frequency offset.

In this part the fading blocks for the simulation are added. Tow different types are implemented a discrete-time one and a statistical one. First to get a better understanding of the material the discrete-time modeled one and after the more realistic one, based on the Rayleight and Rician distribution.

\skelpar[5]{
	Discuss the multitap FIR model we used. How it is possible to set the delay etc. Also mathematics for the interpolation.
}

%To get a basic line for further simulations a 16QAM has been made. The results of this simulation are shown in \figref{fig:simul16QAM} and \figref{fig:simul16QAM_1} as the red Signal. In \tabref{tab:modulation_settings} some importer Parameter settings for a different modulation scheme are mentioned.
%
%A FIR-Filter was added in the Channel to create a time delay between tow paths. In \figref{fig:simul16QAM} the result includes a direct path and a delayed one. In the plot of \figref{fig:simul16QAM_1} the transmission line dosn't include a direct path. %It's impotent to mention that the delay should be smaller than the symbol rate or a multiple of it. (Stimmt dies , not sure any more)
%
%For the a first simulation with some fading the 16QAM simulation model has been extended with a FIR-Filter in the Chanel. The results of this simulation are shown in \figref{fig:simul16QAM} and \figref{fig:simul16QAM_1} as the blue Signal.

\subsection{Fading with Discrete-time model}

For the statical version according to \ref{sec:discrete-time-model} to implement and illustrat the fading effect, a separate block was created and implemented in the channel. Nearer shown in \ref{lst:fir-block}. This block is based on a FIR filter. It can be displayed with a direct path or without one. With the help of this filter, the delay of the line of side paths are illustrated. In this block it is possible to simulate any number of these paths with different strengths, as long as there is an associated amplitude specified for each delayed path. Unfortunately, these simulation values do not correspond to the realety, because too many incalculable side effects occur, which aren't possiple to ilustrate in this simulation.
This block was additionally implemented with the method described in \ref{sec:fractional-delay} to allow non-integer delay values compared to the samples shown in \figref{fig:fractional-delay-sinc-plot}. Where the sinc function does not select an integer sample. Which in turn means that the other sampled values do not add up to zero.
Thus, they will be distributed among the other whole numbers. A window function could also be implemented to limit these values. Here none was implemented because the sinc function is restricted.

\skelpar[5]{
	Discrabe a perfect plot
}

%\begin{german}
%Für die statische gemäss \ref{sec:discrete-time-model} implementation und veranschaulichen des Fading effekts wurde ein eigener Block kreaiert und in den Channel implementiert \ref{lst:fir-block}. Dieser Block basiert auf einem FIR Filter. Es kann mit direcktem Pfad oder ohne dargestellt werden (Line of Sight). Mit Hilfe dieses Filters wird die Verspätung, anhand Zeitlicher delay, der nebenpfaden dargestellt. Es ist möglich beliebig viele dieser Pfade mit unterschiedlicher stärke zu simulieren, solange das für jeden verspäteten Pfad eine dazugehörige Amplute angegeben wird. Leider werden diese simulationswerte nicht mit der Wircklichkeit übereinstien, da zuviele unkalkulierbaren nebeneffekte auftretten.
%Dieser Block wurde zusätzlich mit der methode in \ref{sec:fractional-delay} beschriben implementiert um nichtganzahlige delay Werte zuzulassen, im vergleich zu den sampels zu erlauben, welche in \figref{fig:fractional-delay-sinc-plot} dargestellt sind. Bei dem die sinc-funktion auf kein ganzzahliges sample fählt. Was widerum bedeutet das bei den anderen abgetasteten Werten nicht null ergeben.
%\end{german}

\begin{lstlisting}[
	texcl = true, language = python, escapechar = {`},
	float, captionpos = b, label = {lst:fir-block},
	caption = {
		Block FIR Filter function referenced in listing \ref{lst:phasecorr-work}.
	},
	]
	def work(self, input_items, output_items):
		inp = input_items[0]
		oup = output_items[0]
		# reads the moduled signal 
		if len(self.amplitudes) != len(self.delays):  
			raise Exception("Amplitudes and Delay length dont match")
		if np.min(self.delays)<0:  
			raise Exception("Delay can't be negativ")
		# Some test that there are as match delays as amplitudes and no negative filter values
		max_order = 2 * np.floor(np.max(self.delays)) + 1 
		# find the maximal order of the Filter
	  	max_samples = np.arange(0, max_order +1) 
		max_len = len(max_samples) 
		sum_x = np.zeros(int(max_len))
	
		for (a,d) in zip(self.amplitudes,self.delays):
			order = 2 * np.floor(d) + 1
			samples = np.arange(0, order +1)  
			h = a*(np.sinc(samples-d)) 
			# impuls respond 
			h_len = np.concatenate([h, np.zeros(max_len-len(h))])
			sum_x += h_len
			# add them al together to have one array with al amplitudes one the position of the delay (in samples)
		sum_x[0] = self.los
		# if there is no direct line ad the first delayed amplitude to a zero else to a one.
		y = np.convolve(inp, sum_x)
		y += np.concatenate([self.temp,np.zeros(len(y)-len(self.temp))])
		oup[:] = y[:len(inp)]
		self.temp = y[len(inp):]   
		# The output need to have the same lenght as the input.     		
		return len(oup)
	
\end{lstlisting}

% Bild einfügen 
\subsection{Fading with Statistical model}

	

% TODO: Quelle https://ch.mathworks.com/help/comm/ug/fading-channels.html?searchHighlight=rician%20fading&s_tid=srchtitle_rician%2520fading_2#a1070327427b1


\begin{german}
	Um den effect des multipaht fadinngs nicht nur statisch darzu stellen, wurde ein zweites model kreiert mit hilfe des Frequency Selective Fading Models von Gnu Radio, gemäss \ref{statistical_model}.
	Welcher nach dem Algorthmus aud dem paper \cite{Alimohammad2009} implementiert wurde. 
	Um die resultate einigermassse nach vollziehen zu können wurde ein MATLAP model zur veranschaulichung erstelle.
	Um ein realistisches beispiel zu haben wurden werte aus dem Skript \cite{Mathis} genomen
	
\end{german}

Some realistic value for this block are:
The first delay when theirs non line of side is zero. The second delayed path depend on the environment of measurement. In an indoor enviroment it is usaely between \(1\cdot10^{-9}\) to \(1\cdot10^{-7}\) and in an outdoor environment between \(1\cdot10^{-7}\) to \(1\cdot10^{-5}\). The rest depends on on the bandwith 

\skelpar[5]{
	Simulation mit Werten aus dem Skript
}

\subsection{Measurements}

\subsection{Empirical BER} \label{sec:ber}

%
%
%

\begin{figure}
	\includegraphics[width=\linewidth]{./figures/pdfs/qam_nogui.pdf}
	\caption{GNU Radio Blocks}
	\label{fig:simul16QAM_block}	
\end{figure}

\begin{figure}
	\includegraphics[width=\linewidth]{./figures/screenshots/QAM16_Fading_2.png}
	\caption{Simulation results}
	\label{fig:simul16QAM}	
\end{figure}

\begin{figure}
	\includegraphics[width=\linewidth]{./figures/screenshots/QAM16_Fading_2_.png}
	\caption{Simulation results}
	\label{fig:simul16QAM_1}	
\end{figure}

\begin{figure}
	\includegraphics[width=\linewidth]{./figures/pdfs/qam_Hardware_1711.pdf}
	\caption{GNU Radio Blocks Hardware}
	\label{fig:simul16QAM_Hardware_Aufbau}	
\end{figure}

\begin{figure}
	\includegraphics[width=\linewidth]{./figures/screenshots/QAM16_Hardware_1711.png}
	\caption{Hardware results}
	\label{fig:simul16QAM__Hardware}	
\end{figure}

% To Do: Picture of the setup