Thesis: First part Beacon Disciplining going to Harm

This commit is contained in:
Eric Teunis de Boone 2023-06-30 00:40:18 +02:00
parent 68dffa658d
commit cb3ad3c85c
2 changed files with 132 additions and 120 deletions

View File

@ -10,56 +10,68 @@
\begin{document}
\chapter{Disciplining with a 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.
The detection of extensive air showers uses detectors distributed over large areas.%<<<
Solutions for precise timing over large distances exist for wire setups, i.e.~White Rabbit.
However, the combination of large distances and the number of detectors make it prohibitively expensive to realise such a setup.
For this reason, the time synchronisation of these autonomous stations is typically performed with a \gls{GNSS} clock in each station.
\\
While obtaining 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} was found to range between a few nanoseconds upto 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 employ radio measurements for \Xmax~determination in these experiments.
\\
% 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 a 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}.
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 succesfully employed in \gls{AERA} reaching an accuracy better than $2 \ns$ \cite{PierreAuger:2015aqe}.
\\
Apart from introducing (and controlling) the transmitter ourselves, other sources could introduce similar usable signals.
However, for such signals to work, they must have a well-determined and stable origin.
\\
% 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..
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 beacon (sine) or one that is emitted at some interval (pulse).
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}
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 timing problem\Todo{rephrase} common to both scenarios is worked out.%>>>
\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}
\\
\section{The Timing Problem} %<<<
% 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.
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 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{figure}
\centering
\includegraphics[width=0.6\textwidth,height=0.4\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).
\Todo{use `real' transmitter and radio for schematic}
}
\label{fig:beacon_spatial_setup}
\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{\small T} - \vec{ \small A_i} }\right| }{c} n_{eff}
(\tProp)_i = \frac{ \left|{ \vec{x}_{T} - \vec{x}_{A_i} }\right| }{c} n_\mathrm{eff}
,
\end{equation}
where $n_{eff}$ is the effective refractive index over the trajectory of the signal.
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
@ -122,34 +134,32 @@ Thus, measuring $(\tMeasArriv)_i$ and determining $(\tProp)_i$ for two antennas
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.
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.
\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
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 clock deviations $(\tClock)_i$.
\\
% floating offset, minimising total
\Todo{floating offset, matrix minimisation?}
%\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?}
In the following sections, two approaches for measuring $(\tMeasArriv)_i$ are examined.
%%%% >>>
%%%% Pulse
@ -160,20 +170,20 @@ If the stability of the clock allows for it, the synchronisation can be performe
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.
The detection of a pulse is conceptually simple, and can be accomplished while working fully in the time-domain.
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.
We can characterise the response of a filter as the response to an impulse.
This impulse response can then be used as template to match against a measured waveform.
In Figure~\ref{fig:pulse:filter_response}, the impulse and the filter's response are shown, where the Butterworth filter bandpasses the signal between $30\MHz$ and $80\MHz$.
\\
A measured waveform will consist of the filtered signal in combination with noise.
@ -182,84 +192,119 @@ Figure~\ref{fig:pulse:simulated_waveform} shows an example of the waveform obtai
\\
\begin{figure}
\begin{subfigure}{0.45\textwidth}
\centering
\begin{subfigure}{0.48\textwidth}
\includegraphics[width=\textwidth]{pulse/filter_response.pdf}
\caption{
The filter response.
The amplitudes are not to scale.
The impulse response of the used filter.
Amplitudes are not to scale.
}
\label{fig:pulse:filter_response}
\end{subfigure}
\hfill
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{pulse/antenna_signal_to_noise_6.pdf}
\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}
\\
\begin{subfigure}{\textwidth}
\includegraphics[width=\textwidth]{pulse/filter_signal_correlation.pdf}
\caption{
}
\label{fig:pulse_correlation}
\end{subfigure}
\caption{
Left: A single impulse and the Butterworth filtered signal available to the digitiser in a detector.
Left: A single impulse and a simulated filtered signal, using a Butterworth filter, 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 (see Section~\ref{sec:correlation}) between the two signals.
Detecting the modeled 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}).
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.
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 smaller than the sampling period cannot be resolved.
Since the filtered signal is sampled discretely, this means the start of the
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 modeled 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.
\\
% pulse finding: time accuracies
\begin{figure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{pulse/correlation_tdt0.2_zoom.pdf}
\caption{
Top: The measured waveform and templated filter response from Figure~\ref{fig:pulse:filter_response}.
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 root mean square 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.
\\
\subsubsection{Simulation}
\Todo{remove heading?}
% simulation
From the above, it is clear that both the \gls{SNR} aswell 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\mathrm{ns}$ and a randomised time-offset $t_\mathrm{true}$.
\\
Second, the matching template is created by sampling the ``analog'' template at the specified sampling rate.
\\
% pulse finding: time accuracies
Afterwards, simulated waveforms are correlated against the matching template obtaining a best time delay $\tau$ per waveform.
Comparing the best time delay $\tau$ with the randomised time-offset $t_\mathrm{true}$, we get a time residual $t_\mathrm{res} = t_\mathrm{true} - \tau$ per waveform.
\\
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 the quantisation and the noise, we fit a gaussian to the histograms.
The width of each gaussian gives us an accuracy on how well the time offset is determined from 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.pdf}
\caption{}
\caption{\gls{SNR} = 5}
\label{}
\end{subfigure}
\hfill
\begin{subfigure}{0.45\textwidth}
\begin{subfigure}{0.47\textwidth}
\includegraphics[width=\textwidth]{pulse/time_residuals/time_residual_hist_tdt1.0e-02_n5.0e+01.pdf}
\caption{}
\caption{\gls{SNR} = 50}
\label{}
\end{subfigure}
\caption{
Time residuals histogram
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}
\label{fig:pulse:snr_histograms}
\end{figure}%>>>
By evaluating the time residuals 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 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.
Pulse timing accuracy obtained by matching a templated pulse for multiple template sampling rates to $N=500$ waveforms sampled at $2\ns$.
Dashed lines indicate the asymptotic best time accuracy ($\Delta t/\sqrt{12}$) per template sampling rate.
\Todo{fit curves?, remove dashed line at 1ns}
}
\label{fig:pulse_snr_time_resolution}
\label{fig:pulse:snr_time_resolution}
\end{figure}
% dead time
@ -365,7 +410,7 @@ This relies on the ability of counting how many beacon periods have passed since
\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 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).
}
\label{fig:beacon:pa}
@ -375,7 +420,7 @@ This relies on the ability of counting how many beacon periods have passed since
\bigskip
% Yay for the sine wave
In the following section, the latter scenario of a (single) sine wave as a beacon is worked out.
In the following section, the latter scenario of sine wave beacons 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.
@ -393,43 +438,10 @@ The digital measurement of the beacon phase is dependent on at least two factors
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}% <<<
%\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$.
@ -437,19 +449,12 @@ Such an algorithm efficiently finds the amplitudes and phases within a trace $x$
% .. 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?}
However, if the frequency of interest is not covered in the specific frequencies $k f_s$, 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 the \acrlong{DTFT} for this frequency directly.
\\
\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?}
@ -467,7 +472,7 @@ opening the way to efficiently measuring the phases in realtime.\Todo{figure?}
\subsubsection{Signal to Noise}% <<<
% Gaussian noise
The phase measurement by employing \eqref{eq:fourier:dtft} is influenced by noise in the detector traces.
The phase measurement 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.

View File

@ -120,6 +120,8 @@
\newcommand{\beaconfreq}{\ensuremath{f_\mathrm{beacon}}}
\newcommand{\Xmax}{\ensuremath{X_\mathrm{max}}}
%% time variables
\newcommand{\tTrue}{t}
\newcommand{\tMeas}{\tau}
@ -154,8 +156,10 @@
\newacronym{LOFAR}{LOFAR}{Low-Frequency Array}
\newacronym{PA}{PA}{Pierre~Auger}
\newacronym{PAObs}{PAO}{Pierre~Auger Observatory}
\newacronym{Auger}{Auger}{\acrlong{PAObs}}
\newacronym{AERA}{AERA}{Auger Engineering Radio~Array}
\newacronym{ADC}{ADC}{Analog-to-Digital~Converter}
%% >>>>
%% <<<< Math
\newacronym{DTFT}{DTFT}{Discrete Time Fourier Transform}
@ -163,5 +167,8 @@
\newacronym{FFT}{FFT}{Fast Fourier Transform}
\newacronym{FT}{FT}{Fourier Transform}
\newacronym{RMS}{RMS}{root-mean-square}
\newacronym{SNR}{SNR}{signal-to-noise ratio}
%% >>>>
%% >>>