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

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

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

\begin{document}
\chapter{Waveform Analysis Techniques}
\label{sec:waveform}
%Electric fields,
%Antenna Polarizations,
%Frequency Bandwidth,

Radio antennas are sensitive to changes in their surrounding electric fields.
The polarisation of the electric field that a single antenna can record is dependent on the geometry of this antenna.
Therefore, in experiments such as \gls{Auger} or \gls{GRAND}, multiple antennas are incorporated into a single unit to obtain complementary polarisation recordings.
Additionally, the shape and size of antennas affect how well the antenna responds to certain frequency ranges, resulting in different designs meeting different criteria.
\\

%Waveform time series data
%Sampling,
%Waveform + Time vector,
In each radio detector, the antenna presents its signals to an \gls{ADC} as fluctuating voltages.
In turn, the \gls{ADC} records the analog signals with a specified samplerate $f_s$ resulting in a sequence of digitised voltages or waveform.
The $n$-th sample in this waveform is then associated with a time $t[n] = t[0] + n/f_s = t[0] + n\cdot\Delta t$ after the initial sample at $t[0]$.
%In reality, the sampling rate will vary by very small amounts resulting in timestamp inaccuracies called jitter.
\\

% Filtering before ADC
The sampling is limited by the \gls{ADC}'s Nyquist frequency at half its sampling rate.
In addition, various frequency-dependent backgrounds can be reduced by applying a band-pass filter before digitisation.
For example, in \gls{AERA} and in AugerPrime's radio detector \cite{Huege:2023pfb}, the filter attenuates all of the signal except for the frequency interval between $30 \text{--} 80\MHz$.
\\
In addition to a band-pass filter, more complex filter setups are used to remove unwanted components or introduce attenuation at specific frequencies.
For example, in \gls{GRAND} \cite{GRAND:2018iaj}, the total frequency band ranges from $20\MHz$ to $200\MHz$.
such that the FM broadcasting band ($87.5\MHz \text{--} 108\MHz$) falls within this range.
Therefore, notch filters have been introduced to suppress signals in this band.
\\

% Filter and Antenna response
From the above it is clear that the digitised waveform is a measurement of the electric field that is sequentially convolved with the antenna's and filter's response.
Thus to reconstruct properties of the electric field signal from the waveform, both responses must be known.
\\

% Analysis, properties, frequencies, pulse detection, shape matching,
Different methods are available for the analysis of the waveform, and the antenna and filter responses.
A key aspect is determining the frequency-dependent amplitudes (and phases) in the measurements to characterise the responses and, more importantly, select signals from background.
With \glspl{FT}, these frequency spectra can be produced.
This technique is especially important for the sine wave beacon of Section~\ref{sec:beacon:sine}, as it forms the basis of the phase measurement.
\\
The detection and identification of more complex time-domain signals can be achieved using the cross correlation,
which is the basis for the pulsed beacon method of Section~\ref{sec:beacon:pulse}.
\\

\section{Fourier Transforms}% <<<<
\label{sec:fourier}
\glspl{FT} allow for a frequency-domain representation of a time-domain signal.
In the case of radio antennas, it converts a time-ordered sequence of voltages into a set of complex amplitudes that depend on frequency.
By evaluating the \gls{FT} at appropriate frequencies, the frequency spectrum of a waveform is calculated.
This method then allows to modify a signal by operating on its frequency components, i.e.~removing a narrow frequency band contamination within the signal.
\\

% DTFT from CTFT
The continuous \acrlong{FT} takes the 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) \in \mathcal{R}$ into plane waves with complex-valued amplitude $X(f)$ at frequency $f$.
\\
From the complex amplitude $X(f)$, the phase $\pTrue(f)$ and amplitude $A(f)$ are calculated as
\begin{equation*}
	\begin{aligned}
		\phantom{.}
		\pTrue(f) = \arg\left( X(f) \right), && \text{and} && A(f) = 2 \left| X(f) \right|
		.
	\end{aligned}
\end{equation*}
Note the factor $2$ in this definition of the amplitude.
It is introduced to compensate for expecting a real valued input signal $x(t) \in \mathcal{R}$ and mapping negative frequencies to their positive equivalents.
\\

When $x(t)$ is sampled at discrete times, the integral of \eqref{eq:fourier} is discretized in time to result in the \gls{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)$ is sampled a finite number of times $N$ at times $t[n]$.
Note that the amplitude $A(f)$ will now scale with the number of samples~$N$, and thus should be normalised to $A(f) = 2 \left| X(f) \right| / N$.
\\

Considering a finite sampling size $N$ and periodicity of the signal, the bounds of the integral in \eqref{eq:fourier} have collapsed to $t[0]$ up to $t_{N-1}$.
It follows that the lowest resolvable frequency is $f_\mathrm{lower} = 1/T = 1/(t_{N-1} - 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 = n/f_s$, with $f_s$ the sampling frequency.
Here the highest resolvable frequency is limited by the Nyquist~frequency.
\\

% 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 of the sampling frequency $f = k f_s/N$, becoming the \gls{DFT}
\begin{equation*}
	\label{eq:fourier:dft}
	\phantom{.}
	X(k)
	% = \sum_{n=0}^{N-1}   x(t[n])\, e^{ -i 2 \pi f (t[0] + n/f_s)}
	% = \sum_{n=0}^{N-1}   x(t[n])\, e^{ -i 2 \pi f t[0]}\, e^{ -i 2 \pi f n/f_s}
	% = e^{ -i 2 \pi f t[0]}\, \sum_{n=0}^{N-1}   x(t[n])\, e^{ -i 2 \pi k n}
	= e^{ -i 2 \pi f t[0]} \sum_{n=0}^{N-1}   x(t[n])\, \cdot e^{ -i 2 \pi \frac{k n}{N} }
	.
\end{equation*}

The direct computation of this transform takes $2N$ complex multiplications and $2(N-1)$ complex additions for a single frequency $k$.
When computing this transform for all integer $0 \leq k < N$, this amounts to $\mathcal{O}(N^2)$ complex computations.
\acrlong{FFT}s (\acrshort{FFT}s) are efficient algorithms that derive all $X( 0 \leq k < N)$ in $\mathcal{O}( N \log N)$ calculations.

\begin{figure}
	\begin{subfigure}{0.49\textwidth}
		\includegraphics[width=\textwidth]{methods/fourier/waveforms.pdf}%
		%\caption{}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.49\textwidth}
		\includegraphics[width=\textwidth]{methods/fourier/spectrum.pdf}%
		\label{fig:fourier:dtft_dft}
		%\caption{}
	\end{subfigure}
	\caption{
		\textit{Left:} A waveform sampling a sine wave with white noise.
		\textit{Right:}
			The frequency spectrum of the waveform.
			Comparison of the \gls{DTFT} and \gls{DFT} of the same waveform.
			The \gls{DFT} can be interpreted as sampling the \gls{DTFT} at integer multiple of the waveform's sampling rate $f_s$.
	}
	\label{fig:fourier}
\end{figure}


% Linearity fourier for real/imag
In the previous equations, the resultant quantity $X(f)$ is a complex amplitude.
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 as well,
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(t[n]) \, \cos( 2\pi f t[n] )
		  - i \sum_{n=0}^{N-1} \, x(t[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 \left| X(f) \right| }{N}
		= \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 \arg( X(f) )
		= \arctantwo\left( X_I(f), X_R(f) \right)
	.
\end{equation}

% 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 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.
\\

Figure~\ref{fig:fourier} shows the frequency spectrum of a simulated waveform that is obtained using either a \gls{DFT} or a \gls{DTFT}.
It shows that the \gls{DFT} evaluates the \gls{DTFT} only at certain frequencies.
By missing the correct frequency bin for the sine wave, it estimates both a too low amplitude and the wrong phase for the input function.
\\


% % 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$ up to an overall phase which is dependent on $t[0]$.
Therefore, at the cost of an increased memory allocation, these terms can be precomputed, reducing the number of real multiplications to $2N+1$.

% .. 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 amplitude and phase in realtime.


% >>>>

\section{Cross-Correlation}% <<<<
\label{sec:correlation}
The cross-correlation is a measure of how similar two waveforms $u(t)$ and $v(t)$ are.
By introducing a time delay $\tau$ in one of the waveforms it turns into a function of this time delay,
\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.
% Figure example of correlation and argmax
Figure~\ref{fig:correlation} illustrates how the best time delay $\tau$ between two waveforms can thus be found by finding the maximum cross-correlation.
\\

% Discrete \tau because of sampling
In reality, both waveforms have a finite size, also reducing the time delay $\tau$ resolution to the highest sampling rate of the two waveforms.
When the sampling rates are equal, the time delay variable is effectively shifting one waveform by a number of samples.
\\
% Upsampling? No
Techniques such as upsampling or interpolation can be used to effectively change the sampling rate of a waveform up to a certain degree.
\\

% Approaching analog \tau; or zero-stuffing
Since zero-valued samples do not contribute to the integral of \eqref{eq:correlation_cont}, they can be freely added (or ignored) to a waveform when performing the calculations.
This means two waveforms of different sampling rates can be correlated when the sampling rates are integer multiples of each other, simply by zero-stuffing the slowly sampled waveform.
This allows to approximate an analog time delay between two waveforms when one waveform is sampled at a very high rate as compared to the other.

\begin{figure}
	\centering
	\begin{subfigure}{0.48\textwidth}
		\includegraphics[width=\textwidth]{methods/correlation/waveforms.pdf}
		%\caption{
		%	Two waveforms.
		%}%
		\label{subfig:correlation:waveforms}
	\end{subfigure}
	\hfill
	\begin{subfigure}{0.48\textwidth}
		\includegraphics[width=\textwidth]{methods/correlation/correlation.pdf}
		%\caption{
		%	The correlation of two Waveforms as a function of time.
		%}%
		\label{subfig:correlation}
	\end{subfigure}%
	\caption{
		\textit{Left:} Two waveforms to be correlated with the second waveform delayed by $5$.
		\textit{Right:} The correlation of both waveforms as a function of the time delay $\tau$.
		Here the best time delay (red dashed line) is found at $5$, which would align the maximum amplitudes of both waveforms in the left pane.
	}
	\label{fig:correlation}
\end{figure}

% >>>
\end{document}