mirror of
https://gitlab.science.ru.nl/mthesis-edeboone/m.internship-documentation.git
synced 2024-11-14 03:13:37 +01:00
320 lines
14 KiB
TeX
320 lines
14 KiB
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 (called channels) 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.
|
|
For frequencies at or above this Nyquist frequency, the \gls{ADC} records aliases that appear as lower frequencies in the waveform.
|
|
To prevent such aliases, these frequencies must be removed by a filter before sampling.
|
|
\\
|
|
For air shower radio detection, very low frequencies are also not of interest.
|
|
Therefore, this filter is generally a bandpass filter.
|
|
For example, in \gls{AERA} and AugerPrime's RD\Todo{RD name} the filter attenuates all of the signal except for the frequency interval between $30 \text{--} 80\MHz$.\Todo{citation?}
|
|
\\
|
|
In addition to a bandpass filter, more complex filter setups are used to remove unwanted components or introduce attenuation at specific frequencies.
|
|
For example, in \gls{GRAND}, the total frequency band ranges from $20\MHz$ to $200\MHz$
|
|
such that the FM broadcast band ($87.5\MHz \text{--} 108\MHz$) falls within this range.
|
|
Therefore, notch filters have been introduced to suppress signals in this band.
|
|
\Todo{citation?}
|
|
\\
|
|
|
|
% 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 \acrlong{FT}s these frequency spectra can be produced.
|
|
This technique is especially important for the sinewave 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{Analysis Methods}% <<<
|
|
%\label{sec:waveform:analysis}
|
|
|
|
|
|
\section{Fourier Transforms}% <<<<
|
|
\label{sec:fourier}
|
|
The \gls{FT} allows 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.
|
|
\Todo{citation?}
|
|
%For integer $0 \leq k < N $, efficient algorithms exist that derive all $X( 0 \leq k < N )$ in $\mathcal{O}( N \log N )$ calculations instead of $\mathcal{O}(kcalled \acrlong{FFT}s, sampling a subset of the frequencies.\Todo{citation?}
|
|
|
|
\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 aswell,
|
|
i.e.,~\eqref{eq:fourier:dtft} becomes
|
|
\begin{equation}
|
|
\phantom{.}
|
|
\label{eq:fourier:dtft_decomposed}
|
|
\begin{aligned}
|
|
X(f)
|
|
&
|
|
= X_R(f) + i X_I(f)
|
|
%\\ &
|
|
\equiv \Re(X(f)) + i \Im(X(f))
|
|
\\ &
|
|
= \sum_{n=0}^{N-1} \, x(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 a too low amplitude for the sine wave.
|
|
\\
|
|
|
|
|
|
% % 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$ upto 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 phases in realtime.
|
|
|
|
|
|
% >>>>
|
|
%\section{Pulse Detection}
|
|
|
|
\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.
|
|
|
|
It is defined as
|
|
\begin{equation}
|
|
\label{eq:correlation_cont}
|
|
\phantom{,}
|
|
\Corr(\tau; u, v) = \int_{-\infty}^{\infty} \dif t \, u(t)\, v^*(t-\tau)
|
|
,
|
|
\end{equation}
|
|
where the integral reduces to a sum for a finite amount of samples in either $u(t)$ or $v(t)$.
|
|
Still, $\tau$ remains a continuous variable.
|
|
\\
|
|
% 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.
|
|
However, for the purposes in this document, these methods are not used.
|
|
\\
|
|
|
|
% 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.45\textwidth}
|
|
\includegraphics[width=\textwidth]{methods/correlation/waveforms.pdf}
|
|
%\caption{
|
|
% Two waveforms.
|
|
%}%
|
|
\label{subfig:correlation:waveforms}
|
|
\end{subfigure}
|
|
\hfill
|
|
\begin{subfigure}{0.45\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}
|
|
|
|
% >>>
|
|
\section{Hilbert Transform}% <<<<
|
|
The analytic signal $s_a(t)$ of a waveform $x(t)$ can be obtained using the Hilbert Transform through
|
|
\begin{equation}
|
|
\label{eq:analytic_signal}
|
|
\phantom{,}
|
|
s_a(t) = x(t) + \hat{x}(t)
|
|
,
|
|
\end{equation}
|
|
where $\hat{x}(t)$ is the Hilbert Transformed waveform.
|
|
The Hilbert Transform corresponds to a \gls{FT} where positive frequencies are phase-shifted by $-\pi/2$ and negative frequencies are phase-shifted by $+\pi/2$.
|
|
\\
|
|
|
|
The analytic signal allows to estimate the overall maximum amplitude of a signal irrespective of sign by determining its envelope.
|
|
In Figure~\ref{fig:hilbert_transform}, the envelope of a signal is used to find the time of the maximum amplitude.
|
|
Such a mechanism might be used for timing instead of the cross-correlation described in the previous Section.
|
|
\\
|
|
|
|
\begin{figure}
|
|
\centering
|
|
\includegraphics[width=0.5\textwidth]{pulse/hilbert_timing_interpolation_template.pdf}
|
|
\caption{
|
|
Timing information from the maximum amplitude of the envelope.
|
|
\protect \Todo{noisy trace figure}
|
|
}
|
|
\label{fig:hilbert_transform}
|
|
\end{figure}
|
|
% >>>>
|
|
% >>>
|
|
\end{document}
|