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

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

\begin{document}
\chapter{Synchronising Detectors with a Beacon Signal}
\label{sec:disciplining}
The detection of extensive air showers uses detectors distributed over large areas. %<<<
Solutions for precise timing ($< x\ns$\Todo{fill x}) over large distances exist for cabled setups, e.g.~White~Rabbit~\cite{Serrano:2009wrp}.
However, the combination of large distances and the number of detectors make it prohibitively expensive to realise such a setup for \gls{UHECR} detection.
For this reason, the time synchronisation of these autonomous stations is typically performed with a \gls{GNSS} clock in each station.
\\

To obtain a competitive resolution of the atmospheric shower depth \Xmax with radio interferometry requires an inter-detector synchronisation of better than a few nanoseconds (see Figure~\ref{fig:xmax_synchronise}).
The synchronisation defect in \gls{AERA} using a \gls{GNSS} was found to range between a few nanoseconds up to multiple tens of nanoseconds over the course of a single day (see~\cite[Figure~3]{PierreAuger:2015aqe}).\Todo{copy figure?}
Therefore, an extra timing mechanism must be provided to enable interferometric reconstruction of \gls{EAS}.
\\


% High sample rate -> additional clock
For radio antennas, an in-band solution can be created using the antennas themselves by emitting a radio signal from a transmitter.
With the position of the transmitter known, the time delays can be inferred and thus the arrival times at each station individually.
Such a mechanism has been successfully employed in \gls{AERA} reaching an accuracy better than $2 \ns$ \cite{PierreAuger:2015aqe}.
\\
% Active vs Parasitic
For this section, it is assumed that the transmitter is actively introduced to the array and is therefore fully controlled in terms of produced signals and the transmitting power.
It is foreseeable that ``parasitic'' setups, where sources that are not under control of the experiment introduce signals, can be analysed in a similar manner.
However, for such signals to work, they must have a well-determined and stable origin.\Todo{mention next chapter for auger tv transmitter}
\\

% Impulsive vs Continuous
The nature of the transmitted radio signal, hereafter beacon signal, 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 continuous beacon (e.g.~a~sine~wave) or one that is emitted at some interval (e.g.~a~pulse).
\\

% outline of chapter
In the following, the synchronisation scheme for both the continuous and the recurrent beacon are elaborated upon.
Before going in-depth on the synchronisation using either of such beacons, the synchronisation problem is worked out. %>>>

\section{The Synchronisation Problem} %<<<

% time delay
An in-band solution for synchronising the detectors is effectively a reversal of the method of interferometry in Section~\ref{sec:interferometry}.
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 (see Figure~\ref{fig:beacon_spatial_setup}).
\\

Since the signal is an electromagnetic wave, its instantaneous 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.

\begin{figure}%<<<
	\centering
	\begin{subfigure}{0.49\textwidth}%<<<
		%\centering
		\includegraphics[width=\textwidth,height=\textheight,keepaspectratio]{beacon/antenna_setup_two.pdf}
		\caption{
			Schematic of two antennas ($A_i$) at different distances from a transmitter ($T$).
			Each distance incurs a specific time delay $(\tProp)_i$.
			The maximum time delay difference for these antennas is proportional to the baseline distance (green line).
			\protect \Todo{use `real' transmitter and radio for schematic}
		}
		\label{fig:beacon_spatial_setup}
	\end{subfigure}%>>>
	\begin{subfigure}{0.49\textwidth}%<<<
		%\centering
		\includegraphics[width=\textwidth]{beacon/auger/1512.02216.figure2.beacon_beat.png}
		\caption{
			From Ref~\cite{PierreAuger:2015aqe}.
			The beacon signal that the \gls{Auger} has employed in \gls{AERA}.
			The beating between 4 frequencies gives a total period of $1.1\us$ (indicated by the arrows).
			With a synchronisation uncertainty below $100\ns$ from the \gls{GNSS}, it fully resolves the period degeneracy.
			\protect \Todo{incorporate into text}
		}
		\label{fig:beacon:pa}
	\end{subfigure}%>>>
\end{figure}%>>>

As such, the time delay due to the propagation from the transmitter to an antenna can be written as
\begin{equation}
	\label{eq:propagation_delay}
	\phantom{,}
	(\tProp)_i = \frac{ \left|{ \vec{x}_{T} - \vec{x}_{A_i} }\right| }{c} n_\mathrm{eff}
	,
\end{equation}
where $n_\mathrm{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$.\Todo{different symbols math}
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.
\Todo{text continuity}
\\

\subsection{Phase}
% continuous -> period multiplicity% <<<
In the case of a sine beacon, its periodicty poses an issue.
Differentiating between consecutive periods is not possible using the beacon alone.
This effect is observable in the $\tMeasArriv$ term in \eqref{eq:transmitter2antenna_t0}, being no longer uniquely defined, since
\begin{equation}\label{eq:period_multiplicity}%<<<
	\phantom{,}
	f(\tMeasArriv)
	%= \tTrueArriv + kT\\
	 = f(\frac{\pMeasArriv}{2\pi}T)\\
	 = f(\left[ \frac{\pMeasArriv}{2\pi} + k\right] T)\\
	,
\end{equation}%>>>
with $-\pi < \pMeasArriv < \pi$ the phase of the beacon at time $\tMeasArriv$, $T$ the period of the beacon and the unknown period counter $k \in \mathbb{Z}$.
\\
Of course, this means that the clock defects $\tClock$ can only be resolved up to $\tClock < T$.
this period counter $k$,\Todo{complete}
changing the synchronisation mismatches in \eqref{eq:synchro_mismatch_clocks} to
\begin{equation}\label{eq:synchro_mismatch_clocks_periodic}%<<<
	\begin{aligned}
	\phantom{.}
		(\Delta \tClock)_{ij}
		&\equiv (\tClock)_i - (\tClock)_j \\
		&= (\Delta \tMeasArriv)_{ij} - (\Delta \tTrueArriv)_{ij} \\
		&= (\Delta \tMeasArriv)_{ij} - (\Delta \tProp)_{ij} \\
		&= \left[ \frac{ (\Delta \pMeasArriv)_{ij}}{2\pi} - \Delta k'_{ij} \right] T - (\Delta \tProp)_{ij} \\
		&= \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}%>>>

% lifting period multiplicity
Synchronisation is thus possible with the caveat of being off by an unknown integer amount of periods $\Delta k_{ij}$.
In general, using station $j$ as reference, this amount of periods will be written as $k_i$ for the $i$th station.
%In phase-locked systems this is called onisation.
There are at least two ways to lift this period degeneracy.
\\

% lifting period multiplicity -> long timescale
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 repeated only after $\sim 1 \us$ (see Figure~\ref{fig:beacon:pa}).
With an estimated timing accuracy of the \gls{GNSS} under $50 \ns$ the correct beacon period can be determined, resulting in a unique $\tTrueEmit$ transmit time\Todo{reword}.
\\
% lifing period multiplicity -> short timescale counting +
A second method consists of 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.
A special case of this last scenario where the period counters are approximated from an extensive air shower is worked out in Chapter~\ref{sec:single_sine_sync}.
\\%>>>

\subsection{Array synchronisation}
\Todo{text continuity}

% 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 and exploiting \eqref{eq:transmitter2antenna_t0}.
However, for our purposes relative synchronisation is enough.

% extending to array
In general, we are interested in synchronising an array of antennas.
As \eqref{eq:synchro_mismatch_clocks} applies for each pair of antennas in the array, all the antennas that record the signal can determine the synchronisation mismatches simultaneously.
\footnote{
The mismatch terms for any two pairs of antennas sharing one 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 relative clock deviations $(\Delta \tClock)_{ir}$.
\\
% floating offset, minimising total
%\Todo{floating offset, matrix minimisation?}

\begin{figure}%<<<
	\centering
	\begin{subfigure}{0.49\textwidth}
		\centering
		\includegraphics[width=\textwidth]{beacon/field/field_three_left_time_nomax.pdf}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.49\textwidth}
		\centering
		\includegraphics[width=\textwidth]{beacon/field/field_three_left_phase_nomax.pdf}
	\end{subfigure}
	\caption{
		Reconstruction of a signal's origin (\textit{tx}) or direction using three antennas~($a$,~$b$,~$c$).
		For each location, the colour indicates the total deviation from the measured time or phase differences in the array, such that $0$ (blue) is considered a valid location of \textit{tx}.
		The different baselines allow to reconstruct the direction of an impulsive signal (\textit{left pane}) while a periodic signal (\textit{right pane}) gives rise to a complex pattern (see Appendix~\ref{fig:dynamic-resolve:phase:large} for enhanced size).
	}
	\label{fig:dynamic-resolve}
\end{figure}%>>>




% 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, being impulsive or continuous, requires different methods to determine this quantity.
In the following sections, two separate approaches for measuring the arrival time $(\tMeasArriv)_i$ are examined.
\\
%%%% >>>
%%%% Pulse
%%%%
\section{Pulse Beacon}% <<< Impulsive
\label{sec:beacon:pulse}
If the stability of the clock allows for it, the synchronisation can be performed during a discrete period.
The trade-off 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 very strong signals.
\Todo{rephrase p, order of magnitudes}
\\
Schemes using such a ``ping'' might be employed between the antennas themselves.
Appointing the transmitter role to differing antennas additionally opens the way to (self-)calibrating the antennas in the array.
\\
In this section, the idea of using a single pulse as beacon signal is explored.
\\

% conceptually simple + filterchain response
The detection of a (strong) pulse in a waveform is conceptually simple, and can be accomplished while working fully in the time-domain.
Before recording the signal at a detector, the signal at the antenna is typically put through a filter-chain which acts as a band-pass filter.
This causes the sampled pulse to be stretched in time (see Figure~\ref{fig:pulse:filter_response}).
\\

We can characterise the response of a filter as the response to an impulse.
This impulse response can then be used as a template to match against measured waveforms.
In Figure~\ref{fig:pulse:filter_response}, the impulse and the filter's response are shown, where the Butterworth filter band-passes the signal between $30\MHz$ and $80\MHz$.
\\

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}
	\centering
	\begin{subfigure}{0.48\textwidth}
		\includegraphics[width=\textwidth]{pulse/filter_response.pdf}
		\caption{
			The impulse response of the used filter.
			Amplitudes are not to scale.
			}
		\label{fig:pulse:filter_response}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.48\textwidth}
		\includegraphics[width=\textwidth]{pulse/antenna_signals_tdt0.2.pdf}
		\caption{
			A simulated waveform with noise.
			Dashed lines indicate signal and noise level.
		}
		\label{fig:pulse:simulated_waveform}
	\end{subfigure}
	\caption{
		\textit{Left:} A single impulse and a simulated filtered signal, using a Butterworth filter,  available to the digitiser in a detector.
		\textit{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: template correlation: correlation
Detecting the modelled signal from Figure~\ref{fig:pulse:filter_response} in a waveform can be achieved by finding the correlation (see Section~\ref{sec:correlation}) between the two signals (see Figure~\ref{fig:pulse_correlation}).
The correlation 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.
Therefore, this gives a measure of the best time delay $\tau$ between the two signals.
\\

% pulse finding: template correlation: template and sampling frequency/sqrt(12)
When the digitiser samples the filtered signal, time offsets $\tau$ smaller than the sampling period $\Delta t = 1/f_s$ cannot be resolved.
Still, for many measurements under ideal conditions, one can show that the resolution of the timing asymptotically approaches $\Delta t/\sqrt{12}$.
\\
This is an effect of the quantisation of the sampling period, where the time offsets $\tau$ are modelled as a uniform distribution in time bins the size of $\Delta t$.
In that case, the variance of a uniform distribution applies, obtaining this limit.
\\


\begin{figure}
	\centering
	\includegraphics[width=\textwidth]{pulse/correlation_tdt0.2_zoom.pdf}
	\caption{
		\textit{Top:} The measured waveform and templated filter response from Figure~\ref{fig:pulse:filter_response}.
		\textit{Bottom:} The (normalised) correlation between the waveform and template as a function of time delay $\tau$.
		The template is shifted by the time delay found at the maximum correlation (green dashed line), aligning the template and waveform in the top figure.
	}
	\label{fig:pulse_correlation}
\end{figure}

% pulse finding: signal to noise definition
As can be seen in Figure~\ref{fig:pulse:filter_response}, 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.
\\
Expecting the noise to be gaussian distributed in the time domain, it is natural to use the \gls{RMS} of its amplitude as a quantity representing the strength of the noise.
\\
Therefore, the \gls{SNR} will be defined as the maximum amplitude of the filtered signal versus the \gls{RMS} of the noise amplitudes.
\\


\subsection{Timing accuracy}
% simulation
From the above, it is clear that both the \gls{SNR} as well as the sampling rate of the template have an effect on the ability to resolve small time offsets.
To further investigate this, we set up a simulation\footnote{\Todo{Url to repository}} where templates with different sampling rates are matched to simulated waveforms for multiple \glspl{SNR}.

First, an ``analog'' template is rendered at $\Delta t = 10\mathrm{\,fs}$ to be able to simulate small time-offsets.
Each simulated waveform samples this ``analog'' template with $\Delta t = 2\ns$ and a randomised time-offset $\tTrueTrue$.
\\
Second, the matching template is created by sampling the ``analog'' template at the specified sampling rate (here considered are $0.5\ns$, $0.1\ns$ and $0.01\ns$).
\\
% pulse finding: time accuracies
Afterwards, simulated waveforms are correlated (see \eqref{eq:correlation_cont} in Chapter~\ref{sec:correlation}) against the matching template, this obtains a best time delay $\tau$ per waveform by finding the maximum correlation (see Figure~\ref{fig:pulse_correlation}).
%\\
%Finding the best time delay $\tau$ using \eqref{eq:correlation_cont} corresponds to solving
%\begin{equation}\label{eq:argmax_correlation}
%	\begin{array}
%		\phantom{,}
%		\tau
%		&= \argmax \Corr(\tau; u, v)\\
%		&= \argmax\left( \sum u(t)\, v^*(t-\tau) \right)
%		,
%	\end{array}
%\end{equation}
%where we take the discrete version of \eqref{correlation_cont}.
\\
Comparing the best time delay $\tau$ with the randomised time-offset $\tTrueTrue$, we get a time residual $\tResidual = \tTrueTrue - \tau$ per waveform.
\\
For weak signals ($\mathrm{\gls{SNR}} \lesssim 2$), the correlation method will often select wrong peaks.
Therefore a selection criteria is applied on $\tResidual < 2 \Delta t$ to filter such waveforms and low \glspl{SNR} are not considered here.
\\

Figure~\ref{fig:pulse:snr_histograms} shows two histograms ($N=500$) of the time residuals for two \glspl{SNR}.
Expecting the time residual to be affected by both the quantisation and the noise, we fit a gaussian to the histograms.
The width of each such gaussian gives an accuracy on the time offset $\sigma_t$ that is recovered using the correlation method.
\\

\begin{figure}%<<<
	\centering
	\begin{subfigure}{0.47\textwidth}
		\includegraphics[width=\textwidth]{pulse/time_residuals/time_residual_hist_tdt1.0e-02_n5.0e+00.small.pdf}
		\caption{\gls{SNR} = 5}
		\label{}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.47\textwidth}
		\includegraphics[width=\textwidth]{pulse/time_residuals/time_residual_hist_tdt1.0e-02_n5.0e+01.small.pdf}
		\caption{\gls{SNR} = 50}
		\label{}
	\end{subfigure}
	\caption{
		Time residuals histograms ($N=500$) for $\mathrm{\gls{SNR}} = (5, 50)$ at a template sampling rate of $10 \mathrm{\,ps}$.
		}
	\label{fig:pulse:snr_histograms}
\end{figure}%>>>

By evaluating the timing accuracies $\sigma_t$ for some combinations of \glspl{SNR} and template sampling rates, Figure~\ref{fig:pulse:snr_time_resolution} is produced.
It shows that, as long as the pulse is (much) stronger than the noise ($\mathrm{\gls{SNR}} \gtrsim 5$), template matching could achieve a sub-nanosecond timing accuracy even if the measured waveform is sampled at a lower rate (here $\Delta t = 2\ns$).

\begin{figure}
	\centering
	\includegraphics[width=\textwidth]{pulse/time_res_vs_snr_multiple_dt.pdf}
	\caption{
		Pulse timing accuracy obtained by matching $N=500$ waveforms, sampled at $2\ns$, to a templated pulse, sampled at $\Delta t = 0.1\ns$ (blue), $0.05\ns$ (yellow) and $0.001\ns$ (green).
		Dashed lines indicate the asymptotic best time accuracy ($\Delta t/\sqrt{12}$) per template sampling rate.
		\Todo{points in legend}
		}
	\label{fig:pulse:snr_time_resolution}
\end{figure}

%\begin{figure}%<<<
%		\centering
%		\includegraphics[width=0.5\textwidth]{beacon/auger/1512.02216.figure2.beacon_beat.png}
%		\caption{
%			From Ref~\cite{PierreAuger:2015aqe}.
%			The beacon signal that the \gls{Auger} has employed in \gls{AERA}.
%			The beating between 4 frequencies gives a total period of $1.1\us$ (indicated by the arrows).
%			With a synchronisation uncertainty below $100\ns$ from the \gls{GNSS}, it fully resolves the period degeneracy.
%		}
%		\label{fig:beacon:pa}
%\end{figure}%>>>


% dead time

%%%% >>>
%%%% Sine
%%%%
\section{Sine Beacon}% <<< Continuous
\label{sec:beacon:sine}
% continuous -> can be discrete
In the case the stations need continuous synchronisation, a different approach can be taken.
Still, the following method can be applied as a non-continuous beacon if required.
\\
% continuous -> affect air shower
A continuously emitted beacon will be recorded simultaneously with the signals from air showers.
It is therefore important that the beacon does not fully perturb the recording of the air shower signals, but still be prominent enough for synchronising the antennas.
\\
% Use sine wave to filter using frequency
By implementing the beacon signal as one or more sine waves, the beacon can be recovered from the waveform using Fourier Transforms (see Section~\ref{sec:fourier}).
It is then straightforward to discriminate a strong beacon from the air shower signals, resulting in a relatively unperturbed air shower recording for analysis after synchronisation.
\\
Note that for simplicity, the beacon in this section will consist of a single sine wave at $f_\mathrm{beacon} = 51.53\MHz$ corresponding to a period of roughly $20\ns$.
\\

\Todo{text continuity}
By implementing the beacon signal as one or more sine waves, the beacon signal can be recovered from the waveform using Fourier Transforms (see Section~\ref{sec:fourier}).
\\
% FFT common knowledge ..
The typical Fourier Transform implementation, the \gls{FFT}, finds the amplitudes and phases at frequencies $f_m = m \Delta f$ determined solely by properties of the waveform, i.e.~the~sampling frequency $f_s$ and the number of samples $N$ in the waveform ($0 \leq m < N$ such that $\Delta f = f_s / (2N)$).
\\
% .. but we require a DTFT
Depending on the frequency content of the beacon, the sampling frequency and the number of samples, one can resort to use such a \gls{DFT} \eqref{eq:fourier:dft}.
However, if the frequency of interest is not covered in the specific frequencies $f_m$, the approach must be modified (e.g.~by~zero-padding or interpolation).\Todo{extend?}
Especially when only a single frequency is of interest, a simpler and shorter route can be taken by evaluating the \gls{DTFT} \eqref{eq:fourier:dtft} for this frequency directly.
\\

% Signal to Noise
% frequency domain
%A strong beacon consisting of sine waves will show up as peaks in the frequency spectrum.
%An example spectrum is shown in Figure~\ref{fig:sine:snr_definition}, where
% large amplitudes
Of course, like the pulse method, the ability to measure the beacon's sine waves is dependent on the amplitude of the beacon in comparison to noise.
To quantify this comparison in terms of signal to noise ratio,
 we define the signal level to be the amplitude of the frequency spectrum at the beacon's frequency determined by \gls{DTFT} (the orange line in Figure~\ref{fig:sine:snr_definition}),
 and the noise level as the scaled \gls{RMS} of all amplitudes in the noise band determined by \gls{FFT} (blue line in Figure~\ref{fig:sine:snr_definition}).
Since gaussian noise has Rayleigh distributed amplitudes (see Figure~\ref{fig:noise:pdf:amplitude} in Appendix~\ref{sec:phasor_distributions}), this \gls{RMS} is scaled by $1/\sqrt{2\pi}$.
\\
% longer traces
However, for sine waves, an additional method to increase the \gls{SNR} is available.
In the frequency spectrum, the amplitude with respect to gaussian noise also increases with more samples $N$ in a waveform.
Thus, by recording more samples in a waveform, the sine wave is recovered better.
This effect can be seen in Figure~\ref{fig:sine:snr_vs_n_samples} where the signal to noise ratio increases as $\sqrt{N}$.
\\

% spectral leaking, need strong(?) signals
Note that the \gls{DTFT}, as a finite \gls{FT}, suffers from spectral leakage, where signals at adjacent frequencies influence the ability to resolve the signals separately.
Depending on the signal to be recovered, different windowing functions (e.g.~Hann, Hamming, etc.) can be applied to a waveform.
For simplicity, in this document, no special windowing functions are applied to waveforms.

\begin{figure}%<<<
	\centering
	%\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=0.7\textwidth]{fourier/signal_to_noise_definition.pdf}
		\caption{
			Signal to Noise definition in the frequency domain.
			Solid lines are the noise and beacon's frequency spectra obtained with a \gls{FFT}.
			The noise level (blue dashed line) is the $\mathrm{\gls{RMS}}/\sqrt{2 \pi}$ over all frequencies (blue-shaded area).
			The signal level (orange dashed line) is the amplitude calculated from the \gls{DTFT} at $51.53\MHz$ (orange star).
		}
		\label{fig:sine:snr_definition}
\end{figure}
\begin{figure}
		\centering
	%\end{subfigure}
	%\hfill
	%\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=0.7\textwidth]{fourier/signal_to_noise_vs_timelength.pdf}
		\caption{
			Signal to Noise ratio (SNR) as a function of time for waveforms containing only a sine wave and gaussian noise.
			Note that there is little dependence on the sine wave frequency.
			The two branches (up and down triangles) differ by a factor of two in SNR due to their sampling rate.
		}
		\label{fig:sine:snr_vs_n_samples}
	%\end{subfigure}
	%\\
	%\caption{}
\end{figure}%>>>

\subsection{Timing accuracy}
% simulation

% Gaussian noise
The phase measurement of a sine beacon is influenced by other signals in the recorded waveforms.
They can come from various sources, both internal (e.g.~LNA~noise) and external (e.g.~galactic~background) 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, without introducing a frequency dependence,~i.e.~ white noise.
\\

% simulation waveform
To investigate the resolution of the phase measurement, we generate waveforms of a sine wave with known, but differing, phases $\pTrueTrue$.
Gaussian noise is added on top of the waveform in the time-domain, after which the waveform is band-pass filtered\Todo{list frequencies?} .
The phase measurement of the band-passed waveform is then performed by employing a \gls{DTFT}.
We can compare this measured phase $\pMeas$ with the initial known phase $\pTrueTrue$ to obtain a phase residual $\pResidual = \pTrueTrue - \pMeas$.
\\
In Figure~\ref{fig:sine:trace_phase_measure}, the band-passed waveform and the measured sine wave are shown.
Note that the \gls{DTFT} allows for an implementation where samples are missing by explicitly using the samples' timestamps.
This is illustrated in Figure~\ref{fig:sine:trace_phase_measure} by the cut-out of the waveform.
\\

\begin{figure}
	\centering
	%\begin{subfigure}{0.8\textwidth}
		\includegraphics[width=\textwidth]{fourier/analysed_waveform.zoomed.pdf}
		\caption{
			Band-passed waveform containing a sine wave and gaussian time domain noise and the recovered sine wave at $51.53\MHz$.
			Part of the waveform is removed to verify the implementation of the \gls{DTFT} allowing cut-out samples.
		}
		\label{fig:sine:trace_phase_measure}
	%\end{subfigure}
\end{figure}



Figure~\ref{fig:sine:snr_histograms} shows two histograms ($N=100$) of the phase residuals for a medium and a high \gls{SNR}, respectively.
It can be shown that for medium and strong signals, the phase residual will be gaussian distributed (see below).
The width of each fitted gaussian in Figure~\ref{fig:sine:snr_histograms} gives an accuracy on the phase offset that is recovered using the \gls{DTFT}.
\\

% Signal to Noise definition
\begin{figure}
	\begin{subfigure}{0.47\textwidth}
		%\includegraphics[width=\textwidth]{ZH_simulation/bd_antenna_phase_deltas.py.phase.residuals.c5_b_N4096_noise1e4.pdf}
		\includegraphics[width=\textwidth]{fourier/time_residuals/time_residuals_hist_n7.0e+0.small.pdf}
		\caption{$\mathrm{\gls{SNR}} \sim 7$}
		\label{fig:sine:snr_histograms:medium_snr}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.47\textwidth}
		%\includegraphics[width=\textwidth]{ZH_simulation/bd_antenna_phase_deltas.py.phase.residuals.c5_b_N4096_noise1e3.pdf}
		\includegraphics[width=\textwidth]{fourier/time_residuals/time_residuals_hist_n7.0e+1.small.pdf}
		\caption{$\mathrm{\gls{SNR}} \sim 70$}
		\label{fig:sine:snr_histograms:strong_snr}
	\end{subfigure}
	\caption{
		Phase residuals histograms ($N=100$) for $\mathrm{\gls{SNR}} \sim (7, 70)$.
		For medium to strong signals the phase residuals sample a gaussian distribution.
		\Todo{means not zero}
	}
	\label{fig:sine:snr_histograms}
\end{figure}

% Random phasor sum
For gaussian noise, the resolution of the phase measurement can be shown to be distributed by the following equation
(see Appendix~\ref{sec:phasor_distributions} or \cite[Chapter 2.9]{goodman1985:2.9} for derivation),
\begin{equation}\label{eq:random_phasor_sum:phase:sine}
	\phantom{,}
    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 $s$ is the amplitude of the beacon, $\sigma$ the noise amplitude and $\erf{z}$ the error function.
\cite{goodman1985:2.9} names this equation as ``Constant Phasor plus a Random Phasor Sum''.
For sake of brevity, it will be referred to as ``Random Phasor Sum''.
\Todo{use Phasor Sum instead}
\\
This Random Phasor Sum distribution approaches a gaussian distribution when the beacon amplitude is (much) larger than the noise amplitude.
This can be seen in Figure~\ref{fig:sine:snr_time_resolution} where both distributions are shown for a range of \glspl{SNR}.
There, the phase residuals of the simulated waveforms closely follow the distribution.
\\

From Figure~\ref{fig:sine:snr_time_resolution} we can conclude that depending on the \gls{SNR}, the timing accuracy of the beacon is below $1\ns$ for our beacon at $51.53\MHz$.
Since the time accuracy is derived from the phase accuracy, slightly lower frequencies could be used, but they would require a stronger signal to resolve to the same degree.
Likewise, higher frequencies are an available method of linearly improving the time accuracy.
\\

However, as mentioned before, the period duplicity restricts an arbitrary high frequency to be used for the beacon.
For the $51.53\MHz$ beacon, the next Chapter~\ref{sec:single_sine_sync} shows a method of using an additional signal to counter the period degeneracy of a single sine wave.

\begin{figure}
	\includegraphics[width=\textwidth]{beacon/time_res_vs_snr_large.pdf}
	\caption{
		Timing accuracy for a sine beacon as a function of signal to noise ratio for waveforms of $10240$ samples containing a sine wave at $51.53\MHz$ and white noise.
		It can be shown that the phase accuracies (right y-axis) follow a special distribution~\eqref{eq:random_phasor_sum:phase:sine} that is well approximated by a gaussian distribution for $\mathrm{\gls{SNR}} \gtrsim 3$.
		The green dashed line indicates the $1\ns$ level.
		Thus, for a beacon at $51.53\MHz$ and a $\mathrm{\gls{SNR}} \gtrsim 3$, the time accuracy is better than $1\ns$.
		\Todo{remove title}
	}
	\label{fig:sine:snr_time_resolution}
\end{figure}

% Sine Beacon >>>
\end{document}
