m-thesis-documentation/documents/thesis/chapters/beacon_discipline.tex

% vim: fdm=marker fmr=<<<,>>>
\documentclass[../thesis.tex]{subfiles}

\graphicspath{
	{.}
	{../../figures/}
	{../../../figures/}
}

\begin{document}
\chapter{Disciplining by Beacon} %<<<
\label{sec:disciplining}
Time synchronisation for autonomous stations is typically performed with a \gls{GNSS} clock in each station.
The time accuracy supplied by the \gls{GNSS} clock ($\sim 10 \ns$) is not enough to do effective interferometry.\Todo{citation?}
To cross the $1 \ns$ accuracy threshold an additional timing mechanism is required.
\\

% High sample rate -> additional clock
For radio antennas, an in-band solution can be created using the antennas themselves together with a transmitter.
This is directly dependent on the sampling rate of the detectors.
With the position of the transmitter known, time delays can be inferred and thus the arrival times at each station individually.
Such a mechanism has been previously employed in \gls{AERA} reaching an accuracy better than $2 \ns$ \cite{PierreAuger:2015aqe}.
\\

% Discrete vs Continuous
The nature of the transmitted radio signal, hereafter beacon, affects both the mechanism of reconstructing the timing information and the measurement of the radio signal for which the antennas have been designed..
Depending on the stability of the station clock, one can choose for employing a continous or an intermittent beacon.
This influences the tradeoff between methods.
\\

% outline of chapter
In the following, the synchronisation scheme for both the continuous and intermittent beacon are elaborated upon.
\Todo{further outline}

\section{Timing Problem} %<<<

\begin{figure}
	\centering
	\includegraphics[width=0.6\textwidth,height=0.7\textheight,keepaspectratio]{beacon/antenna_setup_two.pdf}
	\caption{
		An example setup of two antennas ($A_i$) at different distances from a transmitter ($T$).
	}
	\label{fig:beacon_spatial_setup}
\end{figure}

The setup of an additional in-band synchronisation mechanism using a transmitter reverses the method of interferometry.\todo{Requires part in intro about IF}
\\

% time delay
The distance between the transmitter $T$ and the antenna $A_i$ incur a time delay $(\tProp)_i$ caused by the finite propagation speed of the radio signal over these distances.
Since the signal is an electromagnetic wave, its instantanuous velocity $v$ depends solely on the refractive index~$n$ of the medium as $v = \frac{c}{n}$.
In general, the refractive index of air is dependent on factors such as the pressure and temperature of the air the signal is passing through and the frequencies of the signal.
However, in many cases, the refractive index can be taken constant over the trajectory to simplify models.

As such, the time delay due to propagation can be written as
\begin{equation}
	\label{eq:propagation_delay}
	\phantom{,}
	(\tProp)_i = \frac{ \left|{ \vec{\small T} - \vec{ \small A_i} }\right| }{c} n_{eff}
	,
\end{equation}
where $n_{eff}$ is the effective refractive index over the trajectory of the signal.
\\

If the time of emitting the signal at the transmitter $\tTrueEmit$ is known, this allows to directly synchronise the transmitter and an antenna since
\begin{equation}
	\label{eq:transmitter2antenna_t0}
	\phantom{,}
	%$
	(\tTrueArriv)_i
	=
	\tTrueEmit + (\tProp)_i
	=
	(\tMeasArriv)_i - (\tClock)_i
	%$
	,
\end{equation}
where $(\tTrueArriv)_i$ and $(\tMeasArriv)_i$ are respectively the true and measured arrival time of the signal at antenna $A_i$.
The difference between these two terms gives the clock deviation term $(\tClock)_i$.
\\

% relative timing; synchronising without t0 information
As \eqref{eq:transmitter2antenna_t0} applies for each antenna, two antennas recording the same signal from a transmitter will share the $\tTrueEmit$ term.
In that case, the differences between the true arrival times $(\tTrueArriv)_i$ and propagation delays $(\tProp)_i$ of the antennas can be related as
\begin{equation}
	\label{eq:interantenna_t0}
	\phantom{.}
	\begin{aligned}
		(\Delta \tTrueArriv)_{ij}
			&\equiv (\tTrueArriv)_i - (\tTrueArriv)_j \\
			&= \left[ \tTrueEmit + (\tProp)_i \right] - \left[ \tTrueEmit + (\tProp)_j \right] \\
			%&= \left[ \tTrueEmit - \tTrueEmit \right] + \left[ (\tProp)_i - (\tProp)_j \right] \\
			&= (\tProp)_i - (\tProp)_j
			%\\
			%&
			\equiv (\Delta \tProp)_{ij}
	\end{aligned}
	.
\end{equation}



% mismatch into clock deviation
Combining \eqref{eq:interantenna_t0} and \eqref{eq:transmitter2antenna_t0} then gives the relative clock mismatch $\Delta (\tClock)_{ij}$ as
\begin{equation}
	\label{eq:synchro_mismatch_clocks}
	\phantom{.}
	\begin{aligned}
		(\Delta \tClock)_{ij}
		&\equiv (\tClock)_i - (\tClock)_j \\
		&= \left[ (\tMeasArriv)_i - (\tTrueArriv)_i \right] - \left[ (\tMeasArriv)_j - (\tTrueArriv)_j \right] \\
		&= \left[ (\tMeasArriv)_i - (\tMeasArriv)_j \right] - \left[ (\tTrueArriv)_i - (\tTrueArriv)_j \right] \\
		&= (\Delta \tMeasArriv)_{ij} - (\Delta \tTrueArriv)_{ij} \\
		&= (\Delta \tMeasArriv)_{ij} - (\Delta \tProp)_{ij} \\
	\end{aligned}
	.
\end{equation}
Thus, measuring $(\tMeasArriv)_i$ and determining $(\tProp)_i$ for two antennas provides the synchronisation mismatch between them.
\\

% is relative
As the mismatch is the difference between the antenna clock deviations, this scheme does not allow to uniquely attribute the mismatch to one of the clock deviations $(\tClock)_i$.
Instead, it only gives a relative synchronisation between the antennas.
\\
This can be resolved by knowledge on the $\tTrueEmit$ of the transmitter.


\bigskip


% extending to array
In general, we are interested in synchronising an array of antennas.
As \eqref{eq:synchro_mismatch_clocks} applies for any two antennas in the array, all the antennas that record the signal can determine the synchronisation mismatches simultaneously.
\\
The mismatch terms for any two pairs of antennas sharing a single antenna $\{ (i,j), (j,k) \}$ allows to find the closing mismatch term for $(i,k)$ since
\begin{equation*}
	\label{eq:synchro_closing}
	(\Delta \tClock)_{ij} + (\Delta \tClock)_{jk} + (\Delta \tClock)_{ki} = 0
\end{equation*}
Taking one antenna as the reference antenna with $(\tClock)_r = 0$, the mismatches across the array can be determined by applying \eqref{eq:synchro_mismatch_clocks} over consecutive pairs of antennas and thus all clock deviations $(\tClock)_i$.
\\

% floating offset, minimising total
\Todo{floating offset, matrix minimisation?}



% signals to send, and measure, (\tTrueArriv)_i.
In the former, the mechanism of measuring $(\tMeasArriv)_i$ from the signal has been deliberately left out.
The nature of the beacon allows for different methods to determine $(\tMeasArriv)_i$.
In the following, two approaches for measuring $(\tMeasArriv)_i$ are examined.
\Todo{reword towards next sections?}

%%%% >>>
%%%% Pulse
%%%%
\section{Pulse Beacon}% <<<
\label{sec:beacon:pulse}
If the stability of the clock allows for it, the synchronisation can be performed during a discrete period.
The tradeoff between the gained accuracy and the timescale between synchronisation periods allows for a dead time of the detectors during synchronisation.
The dead time in turn, allows to emit and receive strong signals such as a single pulse.
\\
Schemes using such a ``ping'' can be employed between the antennas themselves.
Appointing the transmitter role to differing antennas additionally opens the way to calibrating the antennas in the array.
\\
Note the following method works fully within the time-domain.

% conceptually simple + filterchain response
The detection of a pulse is conceptually simple.
Before recording a signal at a detector, it is typically put through a filterchain which acts as a bandpass filter.
This causes the sampled pulse to be stretched in time (see Figure~\ref{fig:pulse:filter_response}).
\\

The response of a filter is characterised by the response to an impulse.
In Figure~\ref{fig:pulse:filter_response}, an impulsive signal is filtered using a Butterworth filter which bandpasses the signal between $30\MHz$ and $80\MHz$.
The resulting signal can be used as a template to match against a measured waveform.
\\

A measured waveform will consist of the filtered signal in combination with noise.
Due to the linearity of filters, a noisy waveform can be simulated by summing the components after separately filtering them.
Figure~\ref{fig:pulse:simulated_waveform} shows an example of the waveform obtained when summing these components with a considerable noise component.
\\

\begin{figure}
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\textwidth]{pulse/filter_response.pdf}
		\caption{
			The filter response.
			The amplitudes are not to scale.
			}
		\label{fig:pulse:filter_response}
	\end{subfigure}
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\textwidth]{pulse/antenna_signal_to_noise_6.pdf}
		\caption{
			A simulated waveform with noise.
			Dashed lines indicate signal and noise level.
		}
		\label{fig:pulse:simulated_waveform}
	\end{subfigure}
	\caption{
		Left: A single impulse and the Butterworth filtered signal available to the digitiser in a detector.
		Right: A noisy sampling of the filtered signal. It is derived from the filtered signal by adding filtered gaussian noise.
	}
	\label{fig:pulse:waveforms}
\end{figure}

% pulse finding: signal to noise definition
The impulse response spreads the power of the signal over time.
The peak amplitude gives a measure of this power without needing to integrate the signal.
\\
Since the noise is gaussian distributed in the time domain, it is natural to use the root mean square of its amplitude.
\\
Therefore, in the following, the signal-to-noise ratio will be defined as the maximum amplitude of the filtered signal versus the root-mean-square of the noise amplitudes.

\bigskip
% pulse finding: template correlation: correlation
Detecting the modeled signal from Figure~\ref{fig:pulse:filter_response} in a waveform can be achieved by finding the correlation~\eqref{eq:correlation_cont} between the two signals.
This is a measure of how similar two signals $u(t)$ and $v(t)$ are as a function of the time delay $\tau$.
The maximum is attained when $u(t)$ and $v(t)$  are most similar to each other.
This then gives a measure of the best time delay $\tau$ between the two signals.
\\

The correlation is defined as
\begin{equation}
	\label{eq:correlation_cont}
	\phantom{,}
	\Corr(\tau; u,v) = \int_{-\infty}^{\infty} \dif t \, u(t)\, v^*(t-\tau)
	,
\end{equation}
where the integral reduces to a sum for a finite amount of samples in either $u(t)$ or $v(t)$.
Still, $\tau$ remains a continuous variable.
\\

% pulse finding: template correlation: template and sampling frequency/sqrt(12)
When the digitiser samples the filtered signal, time offsets smaller than the sampling period that cannot be resolved.
Since the filtered signal is sampled discretely, this means the start of the

\begin{figure}
	\includegraphics[width=\textwidth]{pulse/waveform+correlation.pdf}
	\caption{
		}
	\label{fig:pulse_correlation}
\end{figure}

% pulse finding: time accuracies

\begin{figure}
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{pulse/time_residuals/time_residual_hist_tdt1.0e-02_n5.0e+00.pdf}
		\caption{}
		\label{}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{pulse/time_residuals/time_residual_hist_tdt1.0e-02_n5.0e+01.pdf}
		\caption{}
		\label{}
	\end{subfigure}
	\caption{
		Time residuals histogram
		}
	\label{fig:pulse_snr_histograms}
\end{figure}

\begin{figure}
	\includegraphics[width=\textwidth]{pulse/time_res_vs_snr_multiple_dt.pdf}
	\caption{
		Pulse timing accuracy obtained by correlating a template pulse for multiple template sampling rates.
		Dashed lines indicate the asymptotic best time accuracy ($\tfrac{1}{f\sqrt{12}}$) per template sampling rate.
		}
	\label{fig:pulse_snr_time_resolution}
\end{figure}

% dead time

%%%% >>>
%%%% Sine
%%%%
\section{Sine Beacon}% <<<
\label{sec:beacon:sine}
% continuous -> can be discrete
In the case that the stations need continuous synchronisation, a different route must be taken.
Still, the following method could be applied as an intermittent beacon if required.
\\
% continuous -> affect airshower
If the beacon is emitted continuously, it will be recorded simultaneously with the signals from airshowers.
The strength of the beacon at each antenna must therefore be tuned such to both be prominent enough to be able to synchronise,
 and only affect the airshower signals recording upto a certain degree\Todo{reword}, much less saturating the detector.
\\


% continuous -> period multiplicity
The continuity of the beacon poses a different issue.
Because the beacon must be periodic, differentiating between consecutive periods is not possible using the beacon alone.
The $\tTrueEmit$ term describing the transmitter time in \eqref{eq:transmitter2antenna_t0} is no longer uniquely defined,
\begin{equation}
	\phantom{,}
	\label{eq:period_multiplicity}
	\tTrueEmit = \left[ \frac{\pTrueEmit}{2\pi} + k\right] T
	,
\end{equation}
with $-\pi < \pTrueEmit < \pi$ the phase of the beacon at time $\tTrueEmit$, $T$ the period of the beacon and $k \in \mathbb{Z}$.
\\
This changes the synchronisation mismatches in \eqref{eq:synchro_mismatch_clocks} to
\begin{equation}
	\label{eq:synchro_mismatch_clocks_periodic}
	\phantom{.}
	\begin{aligned}
		(\Delta \tClock)_{ij}
		&\equiv (\tClock)_i - (\tClock)_j \\
		&= (\Delta \tMeasArriv)_{ij} - (\Delta \tTrueArriv)_{ij} \\
		&= (\Delta \tMeasArriv)_{ij} - (\Delta \tProp)_{ij}  - \Delta k_{ij} T\\
		&= \left[ \frac{ (\Delta \pMeasArriv)_{ij}}{2\pi} - \frac{(\Delta \pProp)_{ij} }{2\pi} - \Delta k_{ij} \right] T\\
		&\equiv \left[ \frac{ (\Delta \pClock)_{ij} }{2\pi} - \Delta k_{ij} \right] T\\
	\end{aligned}
	.
\end{equation}

\begin{figure}
	\begin{subfigure}{\textwidth}
		\includegraphics[width=\textwidth]{beacon/08_beacon_sync_timing_outline.pdf}
		\caption{
			Measure two waveforms at different antennas at approximately the same local time (clocks are not synchronised).
		}
		\label{fig:beacon_sync:timing_outline}
	\end{subfigure}
	\begin{subfigure}{\textwidth}
		\includegraphics[width=\textwidth]{beacon/08_beacon_sync_synchronised_outline.pdf}
		\caption{
			Phase alignment syntonising the antennas using the beacon.
		}
		\label{fig:beacon_sync:syntonised}
	\end{subfigure}
	\begin{subfigure}{\textwidth}
		\includegraphics[width=\textwidth]{beacon/08_beacon_sync_synchronised_period_alignment.pdf}
		\caption{
			Lifting period degeneracy ($k=n-m=7$ periods) using the optimal overlap between impulsive signals.
		}
		\label{fig:beacon_sync:period_alignment}
	\end{subfigure}
	\caption{
		Synchronisation scheme for two antennas using a continuous beacon and an impulsive signal, each emitted from a separate transmitter.
		Grey dashed lines indicate periods of the beacon (orange),
		full lines indicate the time of the impulsive signal (blue).
		\\
		Middle panel: The beacon allows to resolve a small timing delay ($\Delta t_\phase$).
		\\
		Lower panel: Expecting the impulsive signals to come from the same source, the overlap between the two impulsive signals is used to lift the period degeneracy ($k=n-m$).
	}
	\label{fig:beacon_sync:sine}
	\todo{
		Redo figure without xticks and spines,
		rename $\Delta t_\phase$,
		also remove impuls time diff?
	}
\end{figure}

% lifting period multiplicity -> long timescale
Synchronisation is possible with the caveat of being off by an unknown integer amount of periods $\Delta k_{ij}$.
In phase-locked systems this is called syntonisation.
There are two ways to lift this period degeneracy.
\\
First, if the timescale of the beacon is much longer than the estimated accuracy of another timing mechanism (such as a \gls{GNSS}),
 one can be confident to have the correct period.
In AERA \cite{PierreAuger:2015aqe} for example, the total beacon repeats only after $\sim 1 \us$ (see Figure~\ref{fig:beacon:pa}).
With an estimated accuracy of the \gls{GNSS} below $50 \ns$ the correct beacon period can be determined, resulting in a unique $\tTrueEmit$ transmit time\Todo{reword}.
\\

% lifing period multiplicity -> short timescale counting +
Another scheme is using an additional discrete signal to declare a unique $\tTrueEmit$ (see Figure~\ref{fig:beacon_sync:sine}).
This relies on the ability of counting how many beacon periods have passed since the discrete signal has been recorded.

\begin{figure}
		\includegraphics[width=0.5\textwidth]{beacon/auger/1512.02216.figure2.beacon_beat.png}
		\caption{
			From Ref~\cite{PierreAuger:2015aqe}.
			The beacon signal that the \acrlong*{PAObs} has employed in \gls{AERA}.
			The beating between 4 frequencies gives a total period of $1.1\us$ (indicated by the arrows).
		}
		\label{fig:beacon:pa}
\end{figure}


\bigskip

% Yay for the sine wave
In the following section, the latter scenario of a (single) sine wave as a beacon is worked out.
It involves the tuning of the signal strength to attain the required accuracy.
Later, a mechanism to lift the period degeneracy using an airshower as discrete signal is presented.


%%
%% Phase measurement
\subsection{Phase measurement} % <<<
% <<<
A continuous beacon can syntonise an array of antennas by correcting for the measured difference in beacon phases $(\Delta \pMeasArriv)_{ij}$.
They are derived by applying a \gls{FT} to the traces of each antenna.

The digital measurement of the beacon phase is dependent on at least two factors:
 the strength of the beacon in comparison to other signals (such as noise) and the length of the traces.

Additionally, the \gls{FT} can be performed in a number of ways.

These aspects are examined in the following section.

\begin{figure}[h]
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{beacon/sine_beacon.pdf}
		\caption{
			A waveform of a strong sine wave with gaussian noise.\Todo{Add noise}
		}
		\label{fig:beacon:sine}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{fourier/noisy_sine.pdf}
		\caption{
			Fourier Spectrum of the signals.
			\Todo{Add fourier spectra?}
			}
		\label{fig:beacon:spectrum}
	\end{subfigure}
	\\
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{beacon/ttl_beacon.pdf}
		\caption{
			TTL
		}
		\label{fig:beacon:ttl}
	\end{subfigure}

	\caption{
		Both show two samplings with a small offset in time.
		Reconstructing the signal is easier to do for the sine wave with the same samplelength and number of samples.
	}
	\label{fig:beacon:ttl_sine_beacon}
\end{figure}
% >>>
%
% DTFT
\subsubsection{Discrete Time Fourier Transform}% <<<
% FFT common knowledge ..
The typical method to obtain spectral information from periodic data is the \gls{FFT} (a fast implementation of the \gls{DFT} \eqref{eq:fourier:dft}).
Such an algorithm efficiently finds the amplitudes and phases within a trace $x$ at specific frequencies $f_k = f_s \tfrac{k}{N}$ determined solely by the number of samples $N$ ($0 \leq k < N$) and the sampling frequency $f_s$.
\\

% .. but we require a DTFT
Depending on the frequency of the beacon, the sampling frequency and the number of samples, one can resort to use such a \gls{DFT}.
However, if the frequency of interest is not covered in the specific frequencies, the approach must be modified (e.g. zero-padding or interpolation).\Todo{extend?}

Especially when a single frequency is of interest, a shorter route can be taken by evaluating a discretized \gls{FT} directly.
\\

% DTFT from CTFT
The continuous formulation of the \acrlong{FT} takes the following form,
\begin{equation}
	\label{eq:fourier}
	\phantom{.}
	X(f) = \int_\infty^\infty \dif{t}\, x(t)\, e^{-i 2 \pi f t}
	.
\end{equation}
It decomposes the signal $x(t)$ into complex-valued plane waves $X(f)$ of frequency $f$.

When $x(t)$ is sampled at discrete times, the integral of \eqref{eq:fourier} is discretized in time to result in the \acrlong{DTFT}:
\begin{equation}
	%\tag{DTFT}
	\label{eq:fourier:dtft}
	X(f) = \sum_{n=0}^{N-1}   x(t[n])\, e^{ -i 2 \pi f t[n]}
\end{equation}
where $x(t) \in \mathcal{R} $ is sampled at times $t[n]$.
Considering a finite sampling size $N$ and periodicity of the signal, the bounds of the integral in \eqref{eq:fourier} collapse to $t[0]$ up to $t[N]$.
\\
From this it follows that the lowest resolvable frequency is $f_\mathrm{lower} = \tfrac{1}{T} = \tfrac{1}{t[N] - t[0]}$.
\\
Additionally, when the sampling of $x(t)$ is equally spaced, the $t[n]$ terms can be written as a sequence, $t[n] - t[0] = n \Delta t = \tfrac{n}{f_s}$, with $f_s$ the sampling frequency.
The highest resolvable frequency, known as the Nyquist frequency, is limited by this sampling frequency as $f_\mathrm{nyquist} = \tfrac{f_s}{2}$.
\\

% DFT sampling of DTFT / efficient multifrequency FFT
Implementing the above decomposition of $t[n]$, \eqref{eq:fourier:dtft} can be rewritten in terms of multiples $k$ of the sampling frequency, becoming the \acrlong{DFT}
\begin{equation*}
	\label{eq:fourier:dft}
	\phantom{,}
	X(k) = \sum_{n=0}^{N-1}   x[n]\, \cdot e^{ -i 2 \pi {\frac{k n}N} }
	,
\end{equation*}
with $k = \tfrac{f}{f_s}$.
For integer $0 \leq k < N $, efficient algorithms exist that derive all $X( 0 \leq k < N )$ in $\mathcal{O}( N \log N )$ calculations, a~\acrlong{FFT}, sampling a subset of the frequencies.\Todo{citation?}

\bigskip

% Linearity fourier for real/imag
In the previous equations, the resultant quantity $X(f)$ is a complex value.
Since a complex plane wave can be linearly decomposed as
\begin{equation*}
	\phantom{,}
	\label{eq:complex_wave_decomposition}
	\begin{aligned}
		e^{-i x}
		&
		= \cos(x) + i\sin(-x)
		%\\ &
		= \Re\left(e^{-i x}\right) + i \Im\left( e^{-i x} \right)
		,
	\end{aligned}
\end{equation*}
the above transforms can be decomposed into explicit real and imaginary parts aswell,
i.e.,~\eqref{eq:fourier:dtft} becomes
\begin{equation}
	\phantom{.}
	\label{eq:fourier:dtft_decomposed}
	\begin{aligned}
		X(f)
		&
		= X_R(f) + i X_I(f)
		%\\ &
		\equiv \Re(X(f)) + i \Im(X(f))
		\\ &
		=     \sum_{n=0}^{N-1} \, x[n] \, \cos( 2\pi f t[n] )
		  - i \sum_{n=0}^{N-1} \, x[n] \, \sin( 2\pi f t[n] )
		.
	\end{aligned}
\end{equation}

% FT term to phase and magnitude
The normalised amplitude at a given frequency $A(f)$ is calculated from \eqref{eq:fourier:dtft} as
\begin{equation}
	\label{eq:complex_magnitude}
	\phantom{.}
	A(f) \equiv \frac{ 2 \sqrt{ X_R(f)^2 + X_I(f)^2 } }{N}
	.
\end{equation}
Likewise, the complex phase at a given frequency $\pTrue(f)$ is obtained by
\begin{equation}
	\label{eq:complex_phase}
	\phantom{.}
	\pTrue(f) \equiv \arctantwo\left( X_I(f), X_R(f) \right)
	.
\end{equation}
The definition of the amplitude in \eqref{eq:complex_magnitude} contains a factor $2$.
It is introduced to compensate for expecting a real input signal $x(t)$ and mapping negative frequencies to their positive equivalents.
\\

% Recover A\cos(2\pi t[n] f + \phi) using above definitions
Applying \eqref{eq:fourier:dtft_decomposed} to a signal $x(t) = A\cos(2\pi t[n] f + \pTrue)$ with the above definitions obtains
an amplitude $A$ and phase $\pTrue$ at frequency $f$.
When the minus sign in the exponent of \eqref{eq:fourier} is not taken into account, the calculated phase in \eqref{eq:complex_phase} will have an extra minus sign.

\bigskip
% Static sin/cos terms if f_s, f and N static ..
When calculating the \gls{DTFT} for multiple inputs which share both an equal number of samples $N$ and equal sampling frequencies $f_s$, the $\sin$ and $\cos$ terms in \eqref{eq:fourier:dtft_decomposed} are the same for a single frequency $f$.
Therefore, these can be precomputed ahead of time, reducing the number of calculations to $2N$ multiplications.

% .. relevance to hardware if static frequency
Thus, for static frequencies in a continuous beacon, the coefficients for evaluating the \gls{DTFT} can be put into the hardware of the detectors,
opening the way to efficiently measuring the phases in realtime.\Todo{figure?}



% Beacon frequency known -> single DTFT run
% Beacon frequency unknown -> either zero-padding FFT or, DTFT grid search
%When the beacon frequency is known, a single \gls{DTFT} needs to be evaluated.


% Removing the beacon from the signal trace

% >>>
%
% >>>
% Signal to noise
\subsubsection{Signal to Noise}% <<<

% Gaussian noise
The phase measurement by employing \eqref{eq:fourier:dtft} is influenced by noise in the detector traces.
It can come from various sources, both internal (e.g.~LNA~noise) and external (e.g.~radio~communications) to the detector.
A simple noise model is given by gaussian noise in the time-domain, associated to many independent random noise sources.
Especially important is that this simple noise model will affect the phase measurement depending on the strength of the beacon with respect to the noise level.
\\

In the following, this aspect is shortly described in terms of two frequency-domain phasors;
the noise phasor written as $\vec{m} = a \, e^{i\pTrue}$ with phase $-\pi < \pTrue \leq \pi$ and amplitude $a \geq 0$,
and the signal phasor written as $\vec{s} = s \, e^{i\pTrue_s}$, but rotated such that its phase $\pTrue_s = 0$.
\Todo{reword; phasor vs plane wave}
Further reading can be found in Ref.~\cite{goodman1985:2.9}.
\\

% Phasor concept
\begin{figure}
	\label{fig:phasor}
	\caption{
		Phasors picture
	}
\end{figure}

\bigskip

% Noise phasor description
The noise phasor is fully described by the joint probability density function
\begin{equation}
	\label{eq:noise:pdf:joint}
	\phantom{,}
	p_{A\PTrue}(a, \pTrue; \sigma)
	=
	\frac{a}{s\pi\sigma^2} e^{-\frac{a^2}{2\sigma^2}}
	,
\end{equation}
for $-\pi < \pTrue \leq \pi$ and $a \geq 0$.
\\

Integrating \eqref{eq:noise:pdf:joint} over the amplitude $a$, it follows that the phase is uniformly distributed.

Likewise, the amplitude follows a Rayleigh distribution
\begin{equation}
	\label{eq:noise:pdf:amplitude}
	%\label{eq:pdf:rayleigh}
	\phantom{,}
	p_A(a; \sigma)
	%= p^{\mathrm{RICE}}_A(a; \nu = 0, \sigma)
	= \frac{a}{\sigma^2} e^{-\frac{a^2}{2\sigma^2}}
	,
\end{equation}
for which the mean is $\bar{a} = \sigma \sqrt{\frac{\pi}{2}}$ and the standard~deviation is given by $\sigma_{a} = \sigma \sqrt{ 2 - \tfrac{\pi}{2} }$.

\begin{figure}
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{beacon/pdf_noise_phase.pdf}
		\caption{
			The phase of the noise is uniformly distributed.
		}
		\label{fig:noise:pdf:phase}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{beacon/pdf_noise_amplitude.pdf}
		\caption{
			The amplitude of the noise is Rayleigh distribution \eqref{eq:noise:pdf:amplitude}.
		}
		\label{fig:noise:pdf:amplitude}
	\end{subfigure}
	\caption{
		Marginal distribution functions of the noise phasor.
		\Todo{expand captions}
		Rayleigh and Rice distributions.
	}
	\label{fig:noise:pdf}
\end{figure}

\bigskip

% Random phasor sum

In this work, the addition of the signal phasor to the noise phasor will be named ``Random Phasor Sum''.
The addition shifts the mean in \eqref{eq:noise:pdf:joint}
from $\vec{a}^2 = a^2 {\left( \cos \pTrue + \sin \pTrue \right)}^2$
to ${\left(\vec{a} - \vec{s}\right)}^2 = {\left( a \cos \pTrue -s \right)}^2 + {\left(\sin \pTrue \right)}^2$
,
resulting in a new joint distribution
\begin{equation}
	\label{eq:phasor_sum:pdf:joint}
	\phantom{.}
	p_{A\PTrue}(a, \pTrue; s, \sigma)
	= \frac{a}{2\pi\sigma^2}
		\exp[ -
			\frac{
				{\left( a \cos \pTrue - s \right)}^2
				+ {\left( a \sin \pTrue \right)}^2
			}{
				2 \sigma^2
			}
		]
	.
\end{equation}
\\

Integrating \eqref{eq:phasor_sum:pdf:joint} over $\pTrue$ one finds
a Rice (or Rician) distribution for the amplitude,
\begin{equation}
	\label{eq:phasor_sum:pdf:amplitude}
	%\label{eq:pdf:rice}
	\phantom{,}
	p_A(a; s, \sigma)
	= \frac{a}{\sigma^2}
		\exp[-\frac{a^2 + s^2}{2\sigma^2}]
		\;
		I_0\left( \frac{a s}{\sigma^2} \right)
	,
\end{equation}
where $I_0(z)$ is the modified Bessel function of the first kind with order zero.

For the Rician distribution, two extreme cases can be highlighted (as can be seen in Figure~\ref{fig:phasor_sum:pdf:amplitude}).
In the case of a weak signal ($s \ll a$), \eqref{eq:phasor_sum:pdf:amplitude} behaves as a Rayleigh distribution~\eqref{eq:noise:pdf:amplitude}.
Meanwhile, it approaches a gaussian distribution around $s$ when a strong signal ($s \gg a$) is presented.

\begin{equation}
	\label{eq:strong_phasor_sum:pdf:amplitude}
	p_A(a; \sigma) = \frac{1}{\sqrt{2\pi}} \exp[-\frac{{\left(a - s\right)}^2}{2\sigma^2}]
\end{equation}

\begin{figure}
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{beacon/pdf_phasor_sum_phase.pdf}
		\caption{
			The Random Phasor Sum phase distribution \eqref{eq:phasor_sum:pdf:phase}.
		}
		\label{fig:phasor_sum:pdf:phase}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{beacon/pdf_phasor_sum_amplitude.pdf}
		\caption{
			The Random Phasor Sum amplitude distribution \eqref{eq:phasor_sum:pdf:amplitude}.
		}
		\label{fig:phasor_sum:pdf:amplitude}
	\end{subfigure}
	\caption{
			A signal phasor's amplitude in the presence of noise will follow a Rician distribution.
			For strong signals, this approximates a gaussian distribution, while for weak signals, this approaches a Rayleigh distribution.
		\Todo{expand captions}
	}
	\label{fig:phasor_sum:pdf}
\end{figure}

\bigskip
Like the amplitude distribution \eqref{eq:phasor_sum:pdf:amplitude}, the marginal phase distribution of \eqref{eq:phasor_sum:pdf:joint} results in two extremes cases;
weak signals correspond to the uniform distribution for \eqref{eq:noise:pdf:joint}, while strong signals are well approximated by a gaussian distribution.

The analytic form takes the following complex expression,
\begin{equation}
	\label{eq:phase_pdf:random_phasor_sum}
	p_\PTrue(\pTrue; s, \sigma) =
		\frac{ e^{-\left(\frac{s^2}{2\sigma^2}\right)} }{ 2 \pi }
		+
		\sqrt{\frac{1}{2\pi}}
		\frac{s}{\sigma}
		e^{-\left( \frac{s^2}{2\sigma^2}\sin^2{\pTrue} \right)}
		\frac{\left(
			1 + \erf{ \frac{s \cos{\pTrue}}{\sqrt{2} \sigma }}
		\right)}{2}
		\cos{\pTrue}
\end{equation}
where
\begin{equation}
	\label{eq:erf}
	\phantom{,}
	\erf{\left(z\right)} = \frac{2}{\sqrt{\pi}} \int_0^z \dif{t} e^{-t^2}
	,
\end{equation}
is the error function.

\bigskip

\hrule
% Signal to Noise definition
SNR definition

\begin{figure}
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=0.5\textwidth]{ZH_simulation/signal_to_noise_definition.pdf}
		\caption{
			Signal to Noise definition.
		}
		\label{fig:simu:sine:snr_definition}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{ZH_simulation/ba_measure_beacon_phase.py.A74.masked.pdf}
		\caption{
			Phase measurement in a trace with the pulse at $t=$ removed.\Todo{fill t=}
		}
		\label{fig:simu:sine:trace_phase_measure}
	\end{subfigure}
	\caption{}
	\label{fig:simu:sine}
\end{figure}

\begin{figure}
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{ZH_simulation/bd_antenna_phase_deltas.py.phase.residuals.c5_b_N4096_noise1e4.pdf}
		\caption{}
		\label{fig:simu:sine:phase_residuals:medium_snr}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{ZH_simulation/bd_antenna_phase_deltas.py.phase.residuals.c5_b_N4096_noise1e3.pdf}
		\caption{}
		\label{fig:simu:sine:phase_residuals:strong_snr}
	\end{subfigure}
	\caption{
		Phase residuals between the resolved and the true clock phases.
	}
	\label{fig:simu:sine:phase_residuals}
\end{figure}

\begin{figure}
	\includegraphics[width=0.5\textwidth]{beacon/time_res_vs_snr.pdf}
	\caption{
		Measured Time residuals vs Signal to Noise ratio
	}
	\label{fig:time_res_vs_snr}
\end{figure}

% Signal to Noise >>>

% Phase measurement >>>
%
\subsection{Period degeneracy}% <<<
% period multiplicity/degeneracy
A problem with a continuous beacon is resolving the period multiplicity $\Delta k_{ij}$ in \eqref{eq:synchro_mismatch_clocks_periodic}.
It can be resolved by declaring a shared time $\tTrueEmit$ common to the stations in some fashion, and counting the cycles since $\tTrueEmit$ per station.
\\

\bigskip
% Same transmitter 
When the signal defining $\tTrueEmit$ is emitted from the same transmitter that sends out the beacon signal, the number of periods $k$ can be obtained directly from the signal.
If, however, this signal is sent from a different location, the different distances incur different time delays.
In a static setup, these distances should be measured to have a time delay accuracy of less than one period of the beacon signal.\todo{reword sentence}
\\

\bigskip
% airshower gives t0
If measuring the distances to the required accuracy is not possible, a different method must be found to obtain the correct number of periods.
The total time delay in \eqref{eq:phase_diff_to_time_diff} contains a continuous term $\Delta t_\phase$ that can be determined from the beacon signal, and a discrete term $k T$ where $k$ is the unknown discrete quantity.
\\


\begin{figure}
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.maxima.run0.pdf}
		\caption{
			Combined amplitude maxima near shower axis
		}
		\label{fig:findks:maxima}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.reconstruction.run0.power.pdf}
		\caption{
			Power measurement near shower axis with the $k$s belonging to the maximum in the amplitude maxima.
			\Todo{indicate maximum in plot, square figure}
		}
		\label{fig:findks:reconstruction}
	\end{subfigure}
	\\
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.maxima.run1.pdf}
		\caption{
			Maxima near shower axis, zoomed to the location in \ref{fig:findks:reconstruction} with the highest amplitude.
		}
		\label{fig:findks:maxima:zoomed}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{ZH_simulation/findks/ca_period_from_shower.py.reconstruction.run1.power.pdf}
		\caption{
			Power measurement of new 
		}
		\label{}
	\end{subfigure}
	\caption{
		Iterative $k$-finding algorithm:
		First, in the upper left pane, find the set of period shifts $k$ per point that returns the highest maximum amplitude.
	}
	\label{fig:findks}
\end{figure}


\begin{figure}
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{radio_interferometry/trace_overlap/on-axis/dc_grid_power_time_fixes.py.repair_none.axis.trace_overlap.repair_none.pdf}
		\caption{
			Randomised clocks
		}
		\label{fig:simu:sine:period:repair_none}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{radio_interferometry/trace_overlap/on-axis/dc_grid_power_time_fixes.py.repair_phases.axis.trace_overlap.repair_phases.pdf}
		\caption{
			Clock syntonisation
		}
		\label{fig:simu:sine:period:repair_phases}
	\end{subfigure}
	\\
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{radio_interferometry/trace_overlap/on-axis/dc_grid_power_time_fixes.py.no_offset.axis.trace_overlap.no_offset.pdf}
		\caption{
			True clocks
		}
		\label{fig:simu:sine:periods:no_offset}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{radio_interferometry/trace_overlap/on-axis/dc_grid_power_time_fixes.py.repair_full.axis.trace_overlap.repair_full.pdf}
		\caption{
			Full resolved clocks
		}
		\label{fig:simu:sine:periods:repair_full}
	\end{subfigure}
	\caption{
		Trace overlap for a position on the true shower axis.
	}
	\label{fig:simu:sine:periods}
\end{figure}


\begin{figure}
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.repair_none.scale4d.pdf}
		\caption{
			Randomised clocks
		}
		\label{fig:grid_power:repair_none}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.repair_phases.scale4d.pdf}
		\caption{
			Clock syntonisation
		}
		\label{fig:grid_power:repair_phases}
	\end{subfigure}
	\\
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.no_offset.scale4d.pdf}
		\caption{
			True clocks
		}
		\label{fig:grid_power:no_offset}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.repair_all.scale4d.pdf}
		\caption{
			Full resolved clocks
		}
		\label{fig:grid_power:repair_full}
	\end{subfigure}
	\caption{
		Power measurements near the simulation axis with varying degrees of clock deviations.
	}
	\label{fig:grid_power_time_fixes}
\end{figure}


% Period Degeneracy >>>

% Continuous Sine Beacon >>>
% >>>

\bigskip
\chapter{Old work on Sine Beacon}% <<<
\begin{figure}[h]
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{beacon/ttl_beacon.pdf}
		\caption{
			Discrete (square wave) clocks are commonly found in digital circuits.
		}
		\label{fig:beacon:ttl}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\textwidth]{beacon/sine_beacon.pdf}
		\caption{
			A sine wave clock, as will be employed throughout this document.
		}
		\label{fig:beacon:sine}
	\end{subfigure}

	\caption{
		Two different beacon signals with the same frequency.
		Both show two samplings with a small offset in time.
		Reconstructing the signal is easier to do for the sine wave with the same samplelength and number of samples.
	}
	\label{fig:beacon:ttl_sine_beacon}
	\todo{Add fourier spectra?}
\end{figure}


\clearpage
\section{Beacon synchronisation}% <<<

% \delta \phase
As mentioned in Section~\ref{sec:time:beacon}, a beacon consisting of a single sine wave allows to syntonise two antennas by measuring the phase difference of the beacon at both antennas $\Delta \phase = \phase_1 - \phase_2$.
This means the local clock difference of the two antennas can be corrected upto an unknown multiple $k$ of its period, with
\begin{equation}
	\label{eq:phase_diff_to_time_diff}
	\phantom{.}
	\Delta t = \Delta t_\phase + kT = \left(\frac{\Delta \phase}{2\pi} + k\right) T
	.
\end{equation}
By finding a suitably long timescale signal in addition to the sine wave, the amount of periods $k$ can be determined.
\\

In Figure~\ref{fig:beacon_outline}, both such a signal and a sine wave beacon are shown as received at two desynchronised antennas.
The total time delay $\Delta t$ is indicated by the location of the peak of the slow signal.
Part of this delay can be observed as a phase difference $\Delta \phase$ between the two beacons.


% k from coherent sum
\bigskip
The phase difference of the beacon signal obtained in Figure~\ref{fig:beacon_outline} allows to correct small (with respect to the beacon frequency) time delays.
The total time delay may, however, be much larger than one such period.
As shown in \eqref{eq:phase_diff_to_time_diff}, after correcting for the time delay proportional to the phase difference $\Delta t_\phase$, the left-over time delay should be a multiple of the beacon period $kT$.

\bigskip
When the slower signal is transmitted from the transmitter that sent out the beacon signal, then the number of periods $k$ can be obtained directly from the signal.
If, however, the slow signal is sent from a different transmitter, the different distances incur different time delays.
In a static setup, these distance should be measured to such a degree to have a time delay accuracy of about one period of the beacon signal.\todo{reword sentence}
\\

\bigskip
If measuring the distances to the required accuracy is not possible, a different method must be found to obtain the correct number of periods.
The total time delay in \eqref{eq:phase_diff_to_time_diff} contains a continuous term $\Delta t_\phase$ that can be determined from the beacon signal, and a discrete term $k T$ where $k$ is the unknown discrete quantity.
\\

Since $k$ is discrete, the best time delay might be determined from the calibration signal by using a coherent sum
\begin{equation}
	\label{eq:coherent_sum}
	\phantom{,}
	%\chi( t; k) = \sum
	,
\end{equation}
where .., finding the best time delay at the maximum of the sum.
The time delay obtained from the coherent sum

\bigskip
When measuring airshowers, the very signal of the airshower can be used as the calibration signal.
This falls into the dynamic setup described above.
However, while in a static setup the value of $k$ can be estimated from the distances, the distances for each airshower will differ.
\\




%>>>
\end{document}