Thesis: Beacon: small work

This commit is contained in:
Eric Teunis de Boone 2023-05-24 15:57:36 +02:00
parent fcedcaba82
commit d3f659108f
1 changed files with 207 additions and 130 deletions

View File

@ -11,7 +11,7 @@
\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.
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.
\\
@ -246,10 +246,19 @@ Since the filtered signal is sampled discretely, this means the start of the
% pulse finding: time accuracies
\begin{figure}
\includegraphics[width=0.45\textwidth]{pulse/time_accuracy_histogram_snr5.pdf}
\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
\includegraphics[width=0.45\textwidth]{pulse/time_accuracy_histogram_snr50.pdf}
\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}
@ -275,7 +284,7 @@ In the case that the stations need continuous synchronisation, a different route
Still, the following method could be applied as an intermittent beacon if required.
\\
% continuous -> affect airshower
If the beacon must be emitted continuously to be able to synchronise, it will be recorded simultaneously with the signals from airshowers.
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.
\\
@ -326,7 +335,7 @@ This changes the synchronisation mismatches in \eqref{eq:synchro_mismatch_clocks
\begin{subfigure}{\textwidth}
\includegraphics[width=\textwidth]{beacon/08_beacon_sync_synchronised_period_alignment.pdf}
\caption{
Lifting period degeneracy ($k=m-n=7$ periods) using the optimal overlap between impulsive signals.
Lifting period degeneracy ($k=n-m=7$ periods) using the optimal overlap between impulsive signals.
}
\label{fig:beacon_sync:period_alignment}
\end{subfigure}
@ -337,7 +346,7 @@ This changes the synchronisation mismatches in \eqref{eq:synchro_mismatch_clocks
\\
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=m-n$).
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{
@ -355,7 +364,7 @@ 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}.
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 +
@ -376,7 +385,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 scenario of a (single) sine wave as a beacon is worked out.
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.
@ -433,12 +442,12 @@ These aspects are examined in the following section.
\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 = f_s \tfrac{k}{N}$ determined solely by the number of samples $N$ ($0 \leq k < N$) and the sampling frequency $f_s$.
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).
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.
\\
@ -455,7 +464,7 @@ It decomposes the signal $x(t)$ into complex-valued plane waves $X(f)$ of freque
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}
%\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}
@ -465,7 +474,7 @@ Considering a finite sampling size $N$ and periodicity of the signal, the bounds
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 Nyqvist frequency, is limited by this sampling frequency as $f_\mathrm{nyqvist} = \tfrac{f_s}{2}$.
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
@ -529,6 +538,8 @@ Likewise, the complex phase at a given frequency $\pTrue(f)$ is obtained by
\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
@ -561,9 +572,11 @@ opening the way to efficiently measuring the phases in realtime.\Todo{figure?}
\subsubsection{Signal to Noise}% <<<
% Gaussian noise
The traces will contain noise from various sources, both internal (e.g.~LNA~noise) and external (e.g.~radio~communications) to the detector.
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$,
@ -600,7 +613,7 @@ Integrating \eqref{eq:noise:pdf:joint} over the amplitude $a$, it follows that t
Likewise, the amplitude follows a Rayleigh distribution
\begin{equation}
\label{eq:noise:pdf:amplitude}
\label{eq:pdf:rayleigh}
%\label{eq:pdf:rayleigh}
\phantom{,}
p_A(a; \sigma)
%= p^{\mathrm{RICE}}_A(a; \nu = 0, \sigma)
@ -627,8 +640,8 @@ for which the mean is $\bar{a} = \sigma \sqrt{\frac{\pi}{2}}$ and the standard~d
\end{subfigure}
\caption{
Marginal distribution functions of the noise phasor.
Rayleigh and Rice distributions.
\Todo{expand captions}
Rayleigh and Rice distributions.
}
\label{fig:noise:pdf}
\end{figure}
@ -664,7 +677,7 @@ 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}
%\label{eq:pdf:rice}
\phantom{,}
p_A(a; s, \sigma)
= \frac{a}{\sigma^2}
@ -735,6 +748,50 @@ where
\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{
@ -749,10 +806,142 @@ is the error function.
%
\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 >>>
@ -761,27 +950,6 @@ is the error function.
\bigskip
\chapter{Old work on Sine Beacon}% <<<
\Todo{fully rewrite}
The idea of a sine beacon is semi-analogous to an oscillator in electronic circuits.
A periodic signal is sent out from a transmitter (the oscillator), and captured by an antenna (the chip the oscillator drives).
In a digital circuit, the oscillator often emits a discrete (square wave) signal (see Figure~\ref{fig:beacon:ttl}).
A tick is then defined as the moment that the signal changes from high to low or vice versa.
In this scheme, synchronising requires latching on the change very precisely.
As between the ticks, there is no time information in the signal.
\\
\todo{Possibly Invert story from short->long to long->short}
Instead of introducing more ticks in the same time, and thus a higher frequency of the oscillator, a smooth continous signal can also be used.
This enables the opportunity to determine the phase of the signal by measuring the signal at some time interval.
This time interval has an upper limit on its size depending on the properties of the signal, such as its frequency, but also on the length of the recording.
In Figure~\ref{fig:beacon:sine}, both sampling~1~and~2 can reconstruct the sine wave from the measurements.
Meanwhile, the square wave has some leeway on the precise timing.\todo{reword sentence}
\\
\begin{figure}[h]
\begin{subfigure}{0.45\textwidth}
\includegraphics[width=\textwidth]{beacon/ttl_beacon.pdf}
@ -808,53 +976,10 @@ Meanwhile, the square wave has some leeway on the precise timing.\todo{reword se
\todo{Add fourier spectra?}
\end{figure}
%% Second timescale needed
Instead of driving the antenna, the beacon is meant to synchronise the clock of the antenna with the clock of the transmitter.
With one oscillator, the antenna can work in phase with the transmitter, but the actual synchronization can be off by a multiple of periods.
To be able to determine this offset, a second timescale needs to be introduced in the signal.
\\
This slower timescale allows to count the ticks of the quicker signal.\todo{Extend paragraph}
\begin{figure}
\begin{subfigure}{0.45\textwidth}
% \includegraphics[width=0.5\textwidth]{beacon/sine_beacon_multiple_periods_off.pdf}
\caption{
Two syntonised beacons.
The actual synchronization is off by a multiple of periods.
}
\label{fig:second_timescale:off}
\end{subfigure}
\hfill
\begin{subfigure}{0.45\textwidth}
% \includegraphics[width=0.5\textwidth]{beacon/sine_beacon_multiple_periods_off.pdf}
\caption{
Two syntonised beacons, the actual synchronization is off by a multiple of periods.
}
\label{fig:second_timescale:on}
\end{subfigure}
\caption{
}
\label{fig:second_timescale}
\todo{Fill figure and caption}
\end{figure}
\subsection{Beacons in Airshower timing}% <<<
To setup a time synchronising system for airshower measurements, actually only the high frequency part of the beacon must be employed.
The low frequency part, from which the number of oscillations of the high frequency part are counted, is supplied by the very airshower that is measured.
% >>>
\section{Beacon synchronisation}% <<<
As outlined in Section~\ref{sec:time:beacon}, a beacon can also be employed to synchronise the stations.
\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
@ -867,25 +992,6 @@ This means the local clock difference of the two antennas can be corrected upto
By finding a suitably long timescale signal in addition to the sine wave, the amount of periods $k$ can be determined.
\\
\begin{figure}
\centering
\includegraphics[width=\textwidth]{beacon/08_beacon_sync_timing_outline.pdf}
\caption{
Waveforms of a beacon at two antennas, where the clocks have not been synchronised.
Grey dotted lines indicate periods of the sine wave (orange),
full lines indicate the time of the impulsive signal (blue).
Both are sent out from the same transmitter.
The sine wave allows to resolve a small timing delay ($\Delta t_\phase$),
while the impulsive signal allows to calibrate the amount of cycles ($m$,~$n$) the two clocks are separated.
}
\label{fig:beacon_outline}
\todo{
Redo figure without xticks and spines,
rename $\Delta t_\phase$,
also remove impuls time diff
}
\end{figure}
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.
@ -925,35 +1031,6 @@ However, while in a static setup the value of $k$ can be estimated from the dist
\\
\subsection{Lifting period degeneracy}% <<<
\begin{figure}
\begin{subfigure}[t]{0.5\textwidth}
\includegraphics[width=\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.no_offset.scale4d.pdf}
\label{fig:grid_power:no_offset}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.5\textwidth}
\includegraphics[width=\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.repair_none.scale4d.pdf}
\label{fig:grid_power:repair_none}
\end{subfigure}
\\
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.repair_phases.scale4d.pdf}
\label{fig:grid_power:repair_phases}
\end{subfigure}
\hfill
\begin{subfigure}[b]{0.5\textwidth}
\includegraphics[width=\textwidth]{radio_interferometry/dc_grid_power_time_fixes.py.X400.repair_all.scale4d.pdf}
\label{fig:grid_power:repair_all}
\end{subfigure}
\caption{
}
\label{fig:grid_power_time_fixes}
\end{figure}
% >>>
%>>>