Obtaining radar-variable estimates with low errors is important to correctly characterize the underlying meteorological phenomena.
The estimators presented in Chapter \ref{chap_2} yield estimates of sufficient quality provided that the number of independent measurements is large \citep{Melnikov2004, Yu2012}.
In practice, however, it may be extremely difficult to obtain the required number of measurements due to numerous reasons, some of which were mentioned in the previous chapter.
In this chapter, the Bootstrap DPSD estimator is presented as an alternative method that overcomes some of the limitations of the conventional PSD and DPSD estimators.
The bootstrap concept and definitions of non-dependent and dependent data bootstrap are presented in Section \ref{sec:bs_concepts}.
Section \ref{sec:bootwx} presents applications of these methods to weather radar I/Q time-series signals and explores the utility of each method for this application.
In Section \ref{sec:considerations}, the considerations taken into account in the design of the Bootstrap DPSD estimator are presented.
The criteria to evaluate the performance of the methods and a few signal conditioning strategies that improve the quality of the estimates are explained.
Section \ref{sec:analysis} presents the rationale behind the selection of the best strategies and parameters for the Bootstrap DPSD estimator.
The algorithm of the Bootstrap DPSD estimator is described in Section \ref{sec:algorithm}.
And finally, a summary of this chapter is presented in Section \ref{sec:summary}.

\section{Bootstrap Resampling} \label{sec:bs_concepts}
In a statistical investigation, samples of an unknown distribution can be obtained to make inferences about the actual distribution.
In other words, the sample distribution is an estimate of the real distribution.
It is explained by the law of large numbers that the accuracy of an estimate improves as the number of samples increases.
Logically, as the sample size increases, the sample distribution will converge to the real distribution.
However, obtaining a large number of samples may not always be feasible or practical, in which case there is no choice other than to work with the available data.
In such situations, resampling methods can provide a better estimate of the sample distribution because the distribution, albeit small and discrete, is known and thus can be treated as a parent distribution \citep{Hall1992}. 
In this sense, the statistical characteristics of the resample distribution inferred from the sample distribution can be compared to those of the sample distribution inferred from the unknown, real distribution \citep{Efron1994}.
Consequently, more accurate inferences about the real distribution can be made, provided that the sample distribution is reasonably representative of the underlying distribution \citep{Varian2005}.
According to \citet{Efron1994}: ``The bootstrap is a computer-based method for assigning measures of accuracy to statistical estimates".
This concept was first introduced in \citet{Efron1979} as the random resampling with replacement of a given set of samples, and it is typically used to estimate statistical properties of estimators (e.g., variance, confidence intervals). 
The bootstrapping procedure is similar to a Monte Carlo simulation, where independent samples are drawn from a known distribution, and by drawing a large number of them, an approximation to the real distribution is obtained. 
The difference is that the bootstrap generates a large number of resamples from the sample distribution to obtain a bootstrap distribution with which the statistical quality of an estimator can be assessed \citep{Varian2005}.
Intuitively, the bootstrap simulates the probabilistic aspect of the real world by substituting unknowns with estimates from the available data \citep{Zoubir2004}.
Figure \ref{fig:bootdiag} shows a diagram of the bootstrap concept, where the mechanisms of the real world are compared to the bootstrap world.
A proper description will be presented later.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/bootdiag.pdf}
\caption{Bootstrap resampling. Adapted from \citet{Zoubir2004}.}
\label{fig:bootdiag}
\end{figure}

Although the bootstrap in its simplest form can fail under certain conditions (e.g., correlated data), with proper modification it can be widely extended to many applications \citep{Hall1992}.
The bootstrap has been used in many different areas of engineering, including radar signal processing \citep{Zoubir2004} and some applications related to radar meteorology include rainfall modeling through spatial bootstrapping \citep{Overeem2009}, cloud radar observations and tests for similarities between cloud occurrence profiles \citep{Marchand2006}, the determination of the likelihood of TDS parameters observed in two different wavelengths having the same parent distribution \cite{Bodine2014a}, and the analysis of tornado strength through skill scores \citep{Kingfield2015}.
However, to our knowledge, the bootstrap has not been applied to weather radar signal processing.
%As previously mentioned, the bootstrap's forte is to produce confidence intervals of estimators.

There are bootstrap methods for non-dependent and dependent data, a brief description of the methods is provided hereafter.
In a later section, the bootstrap methods are applied to I/Q signals, and it will be clear that only the dependent data methods can yield meaningful results for spectral estimation.
%The results obtained from bootstrapping weather radar signals will be assessed in Section \ref{sec:bootwx}.

\subsection{Non-Dependent Data Bootstrap}
For non-dependent data, i.e., data without structure, the independent and identically distributed (IID) bootstrap can be formulated as follows \citep{Efron1994}.
Let 
\be \mathbf{X} =\{ X_0, X_1, \ldots, X_{M-1} \} \label{eq:iid1} \ee 
be the measured $M$-point sample set from a given distribution $F$. 
A resample 
\be \mathbf{X}' = \{ X'_0, X'_1, \ldots X'_{M-1} \} \label{eq:iid2} \ee 
of  $\mathbf{X}$ is obtained by reassigning the samples of $\mathbf{X}$ randomly and with replacement with probability
\be P(X'_i = X_j | \mathbf{X}) = \frac{1}{M}  \label{eq:iid3} ,\ee 
for $0 \leq i,j \leq M-1$. 
Then, a number $K'$ of resamples is generated, $\mathbf{X}'_0, \ldots, \mathbf{X}'_{K'-1}$, of which the distribution is $\hat{F}$, an estimate of the distribution $F$. 
Hereafter, the terms resample and bootstrap sample will be used interchangeably.
It is shown in \citet{Hall1992} that as $M$ goes to infinity, $\hat{F}$ approaches $F$. 
In other words, the sample distribution converges to the real distribution as the number of samples increases.
The number $K'$ of resamples needed varies depending on various factors which can be application- or estimator-specific \citep{Zoubir2004}.
From these resamples, measures of the accuracy of a statistical property of interest $\hat{\theta}'= s(\mathbf{X}')$ can be obtained.
For example, if the sample mean estimator is used with a single set, only one estimate of the mean could be computed.
However, by bootstrapping the sample mean estimator, the distribution of the mean can be obtained to assess the errors of the measurement.
%The IID bootstrap assumes every sample from the set is independent and identically distributed. 



\subsection{Dependent Data Bootstrap}
For dependent data, such as the I/Q time-series signals from weather radars, the former definition is inadequate. 
However, it can be extended to fit the definition if the resampling is done by taking blocks of consecutive data instead of single points to try to maintain the sample-time correlation of the data.
There are many different block bootstrapping methods, each with its advantages and selection of parameters \citep{Lahiri2003}. 

The general block bootstrap definition is derived from the definition of the IID bootstrap.
By resampling blocks of data instead of individual samples, that is, the original sequence $X_0,\ldots,X_{M-1}$ is divided into $N$ blocks of consecutive samples as
\be \mathbf{X} =\{ \mathbf{B}_0, \mathbf{B}_1, \ldots, \mathbf{B}_{N-1} \} ,\ee
and
%\be \mathbf{B}_j = \{ X_{jL} , \ldots ,  X_{(j+1)L-1} \} , \qquad 0 \leq j < N .\ee
\be \mathbf{B}_j = \{ X_{jL} , \ldots ,  X_{(j+1)L-1} \} .\ee
for $0 \leq j < N$, where $N$ is the number of blocks, and $L$ is the block length.
The different dependent data bootstrap methods are similar in principle, though the block selection scheme varies for each method.
For now, it suffices to say that these blocks of data can be arranged in different ways yielding different results.
Following this, a bootstrap sample becomes 
\be \mathbf{X}' = \{ \mathbf{B}'_0, \mathbf{B}'_1, \ldots \mathbf{B}'_{N-1} \} , \ee 
where $\mathbf{X}'$, in this case, is obtained by reassigning blocks of consecutive samples in $\mathbf{X}$ randomly with replacement, i.e.,  with probability 
\be P(\mathbf{B}'_i = \mathbf{B}_j | \mathbf{X}) = \frac{1}{N} ,\ee for $0 \leq i,j \leq N-1$.
It should be noted that for larger block lengths, these methods are less computationally expensive because of the reduced number of blocks per bootstrap sample.
However, dividing the data into blocks reduces the degrees of freedom to do resampling.

The concept of dividing the time-series signals into blocks may be similar to Bartlett's and Welch's method of PSD estimation.
In fact, two of the block bootstrap methods that will be described hereafter employ a non-overlapping as well as an overlapping block scheme like the aforementioned PSD estimators.
However, neither of them resample the blocks in the way bootstrap does.
Additionally, a circular and a stationary block bootstrap method will be presented.
The former method allows the samples in the blocks to wrap around the data circularly, while the latter uses a dynamic block length.

\subsubsection{Non-Overlapping Block Bootstrap (NBB)} \label{sec:nbb_def}
This method was introduced by \citet{Carlstein1986} and consists in dividing the signal into blocks that do not overlap. These blocks are then resampled with equal probability \citep{Lahiri2003}.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/nbb_diag.pdf}
\caption{Non-overlapping block bootstrap block scheme. Adapted from \citet{Lahiri2003}.}
\label{fig:nbb}
\end{figure}

By definition \citep{Lahiri2003},
\be \mathbf{X} =\{ \mathbf{B}_0, \mathbf{B}_1, \ldots, \mathbf{B}_{N-1} , \ldots, X_{M-1}\} \ee
\be \mathbf{B}_j = \{ X_{jL} , \ldots ,  X_{(j+1)L-1} \} \ee
\be \mathcal{B} = \{ \mathbf{B}_0, \mathbf{B}_1, \ldots, \mathbf{B}_{N-1} \} ,\ee
where $\mathcal{B}$ defines all the possible blocks available for resampling (Figure \ref{fig:nbb}).
A bootstrap sample $\mathbf{X}'$ is obtained as
\be \mathbf{X}' =\{ \mathbf{B}'_0, \mathbf{B}'_1, \ldots \mathbf{B}'_{N-1} \} \ee 
%\be \mathbf{X}' = \begin{cases} \{ \mathbf{B}'_1, \mathbf{B}'_2, \ldots \mathbf{B}'_N \} & \text{if } \langle M \rangle _L = 0 \\ 
%\{ \mathbf{B}'_1, \mathbf{B}'_2, \ldots \mathbf{B}'_{N+1} \}  & \text{if } \langle M \rangle _L \neq 0 \end{cases} \ee 
with
\be P(\mathbf{B}'_i = \mathbf{B}_j | \mathcal{B}) = \frac{1}{N} \ee 
\be N= \left \lfloor \frac{M}{L} \right \rfloor  ,\ee
for $0 \leq i,j \leq N-1$, where $M$ is the length of the signal, $L$ is the block length, $N$ is the number of blocks, and $\lfloor\cdot\rfloor$ is the floor function.
If $L$ is an integer multiple of $M$ (i.e. $NL = M$), there will be exactly $N$ blocks in $\mathbf{X}$.
Otherwise, there will be $N$ blocks of length $L$, and the remaining $M-NL$ samples of $\mathbf{X}$ will not have a block assigned.
In this case, an extra block at index $N$ is drawn and the resulting sequence is truncated to have exactly $M$ data points (as the original sequence).

\subsubsection{Moving Block Bootstrap (MBB)} \label{sec:mbb_def}
In the moving block bootstrap (MBB) \citep{Kunsch1989, Liu1992}, unlike the NBB, the blocks are allowed to overlap as depicted in Figure \ref{fig:mbb}. The blocks are then resampled with equal probability \citep{Lahiri2003}. 
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/mbb_diag.pdf}
\caption{Moving block bootstrap block scheme. Adapted from \citet{Lahiri2003}.}
\label{fig:mbb}
\end{figure}

By definition,
\be \mathbf{B}_j = \{ X_{j} , \ldots ,  X_{j+L-1} \} \ee
\be \mathcal{B} = \{ \mathbf{B}_0, \mathbf{B}_1, \ldots, \mathbf{B}_{N'-1} \} ,\ee
for $0\leq j \leq N'-1$, where $N' = M-L+1$ is the number of available blocks (Figure \ref{fig:mbb}).
A bootstrap sample $\mathbf{X}'$ is obtained as
\be \mathbf{X}' =\{ \mathbf{B}'_0, \mathbf{B}'_1, \ldots \mathbf{B}'_{N-1} \} \ee 
\be P(\mathbf{B}'_i = \mathbf{B}_j | \mathcal{B}) = \frac{1}{N'} \ee 
\be N = \left\lceil \frac{M}{L} \right\rceil ,\ee
for $0 \leq i \leq N-1$, where $N$ is the number of blocks to be drawn from $\mathcal{B}$ for each bootstrap sample, $L$ is the block length, and $\lceil\cdot\rceil$ is the ceiling function.
As with the previous scheme, the bootstrap samples exceeding the length of the original sequence will be truncated to $M$.

\subsubsection{Circular Block Bootstrap (CBB)} \label{sec:cbb_def}
A variant of the MBB, which additionally allows the blocks to wrap around the data circularly, is the circular block bootstrap (CBB) \citep{Lahiri2003} depicted in Figure \ref{fig:cbb}.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/cbb_diag.pdf}
\caption{Circular block bootstrap block scheme.}
\label{fig:cbb}
\end{figure}

By definition, 
\be \mathbf{B}_j = \{ X_{\langle j \rangle _M} , \ldots ,  X_{\langle j+L-1 \rangle _M} \} \ee
\be \mathcal{B} = \{ \mathbf{B}_0, \mathbf{B}_1, \ldots, \mathbf{B}_{M-1} \} .\ee
A bootstrap resample $\mathbf{X}'$ is obtained as
\be \mathbf{X}' = \{ \mathbf{B}'_0, \mathbf{B}'_1, \ldots, \mathbf{B}'_{M-1} \}  \ee
\be P(\mathbf{B}'_i = \mathbf{B}_j | \mathcal{B}) = \frac{1}{M} \ee 
\be N = \left\lceil \frac{M}{L} \right\rceil ,\ee
for $0 \leq i \leq N-1$ and $0 \leq j \leq M-1$, where $N$ is the number of blocks per bootstrap resample.
The bootstrap samples exceeding the length of the original sequence will be truncated to $M$.

\subsubsection{Stationary Block Bootstrap (SBB)} \label{sec:sbb_def}
The stationary block bootstrap (SBB) \citep{Politis1994} makes the assumption of stationarity by dynamically varying the block lengths for each bootstrap sample \citep{Lahiri2003}.
The blocks, however, are not allowed to wrap around circularly.
The length of each block is a random number with geometric distribution with parameter $p$.
This allows for even higher variability and randomness than the previous methods.
%but the samples of consecutive blocks lack coherence for dependent data.
%\begin{figure}[!htbp]
%\centering
%\includegraphics[width=0.5\linewidth]{figures/placeholder.png}
%\caption{Stationary block bootstrap block scheme.}
%\label{fig:sbb}
%\end{figure}
Let the bootstrap resample be drawn as
\be \mathbf{X}' = \{ \mathbf{B}(i_0,j_0),  \mathbf{B}(i_1,j_1), \ldots \} \ee
\be \mathbf{B}(i,j) = \{ X_i, \ldots, X_{i+j-1} \} ,\ee
where $i$ defines the starting index of the block, and $j$ defines the length of the block.
Furthermore, $i$ has an uniform distribution:
\be i \sim U(0,M-1) ,\ee
and $j$ has a geometric distribution with probability parameter $p$:
\be P(j = k) = p(1-p)^{k-1} .\ee
For each bootstrap sample, blocks are continuously drawn until the length of the resample is greater than the original sequence.
Then, the resulting sequence is truncated to have a length equal to $M$.

\subsubsection{Other Methods}
Bootstrapping can also be performed in the frequency domain rather than in the time domain.
There are a variety of methods within this group of which the general procedure is to assume a spectral model and then to bootstrap the residuals, i.e., the difference between the model and the spectral estimate.
This procedure yields bootstrapped PSDs directly, as opposed to bootstrapped time-series signals.
Additionally, bootstrapped time-series signals may be obtained with the inverse Fourier transform of the bootstrapped PSD.
Although they make no assumption on data dependence in the time domain, these methods rely on a good spectral model in order to be effective.
Obviously, this is unsuitable for DPSD estimation since the model of the underlying phenomena is unknown, and hence will not be discussed in this work.

\section{Bootstrapping I/Q Time-Series Signals} \label{sec:bootwx}
The diagram in Figure \ref{fig:iqbootdiag1} illustrates the bootstrap process applied to I/Q time-series signals.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/iqbootdiag.pdf}
\caption{Diagram of the bootstrap. Adapted from \citet{Zoubir2004} for I/Q time-series signals.}
\label{fig:iqbootdiag1}
\end{figure}
After gathering the measurements for a dwell, instead of estimating the PSD as per the classical approach, bootstrap samples of the I/Q signals are generated and the PSD of the resamples are estimated.
This yields $K'$ spectra per dwell, which are averaged to compute the DPSD estimates.
The process of averaging bootstrapped estimates is known as bootstrap aggregating \citep{Breiman1996}.
Hereafter, the bootstrap resamples of the I/Q signals shall be referred to as pseudo-realizations.
For weather radar signals, although the samples of the I and Q components are uncorrelated with each other, the random processes controlling the changes in I and Q are correlated \citep{Doviak2006}.
Therefore, the structure of the pseudo-realizations must be similar to the original sequence in order to correctly convey the sample-time coherence.
For this reason, it is evident that the IID bootstrap assumption is inadequate and cannot be applied to estimate the spectra of weather signals.
The different block bootstrapping methods have a different selection of blocks available for resampling, which is why a preliminary assessment is necessary to determine the method that yields the best results for DPSD estimation.
This section illustrates the block bootstrap methods applied to polarimetric weather radar I/Q time-series signals and the DPSD estimates obtained with each method.
At the end of this section a summary of methods is presented.

\subsection{IID Bootstrap} \label{iid_boot_iq}
Given the I/Q signal for the H- and V- channels from Equation (\ref{eq:iqsigs}), 
\be \mathbf{V}_\mathrm{H,V} = [ V_\mathrm{H,V}(0), \ldots, V_\mathrm{H,V}(M-1) ] \ee
and the definitions in Equations (\ref{eq:iid1}-\ref{eq:iid3}), a pseudo-realization is obtained as
\be \mathbf{V}'_\mathrm{H,V} = [ V'_\mathrm{H,V}(0), \ldots, V'_\mathrm{H,V}(M-1) ] , \ee
where 
\be P(V'_\mathrm{H,V}(m') = V_\mathrm{H,V}(m) | \mathbf{V}_\mathrm{H,V}) = \frac{1}{M} ,\ee
and $0 \leq m, m' \leq M-1$. 
%$\mathbf{V}'_\mathrm{H,V}$ is the $i$-th bootstrap pseudo-realization for either the H- or the V-channel accordingly, though the index has been omitted for simplicity.
In this sense, a pseudo-realization is obtained from $\mathbf{V}$ by resampling randomly with replacement, which means the data points are drawn with equal probability, and as such, they can be repeated in the sequence (Figure \ref{fig:iid_iq}).
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/bs_iq_example_iid.eps}
\caption{Obtaining I/Q pseudo-realizations with the IID bootstrap.}
\label{fig:iid_iq}
\end{figure}
The IID bootstrap is a particular case of a block bootstrap with unit block length.

To assess the quality of the spectral estimates obtained from each bootstrap method, simple simulations for a single signal are run.
The simulation parameters of the bulk radar variables of the signal are summarized in Table \ref{tab:simpar1}. 
\begin{table}
\centering
\begin{tabular}{|c|c|}
\hline
Parameter & Value \\
\hline
SNR & $30$ dB \\
$\bar{v}_r$ & $5$ m/s \\
$\sigma_v$ & $2$ m/s\\
$Z_\mathrm{DR}$ & $1.5$ dB \\
$\rho_\mathrm{HV}$ & $0.9$\\
$M$ & $100$ \\
$K$ & $1$ independent spectrum \\
$K'$ & $100$ pseudo-realizations \\
$N_i$ & $1000$ iterations \\
$v_a$ & $15.7$ m/s \\
\hline
\end{tabular}
\caption{Simulation parameters for a simple single signal case.}
\label{tab:simpar1}
\end{table}
The parameters are chosen to resemble a realistic weather signal, while $\rho_\mathrm{HV}$ is arbitrarily set to a low value so that potential biases become more evident. 
$N_i$ independent realizations of a signal are simulated and, for each of them, $K'$ pseudo-realizations are bootstrapped. 
The PSD of the pseudo-realizations are estimated using the periodogram with a Blackman-Nuttall window, and then averaged to compute the DPSD assuming only one independent spectrum is available.
Simply put, for each independent realization, the PSD is estimated using the average of $K'$ bootstrapped spectra, and $K$ independent spectra (in this case, $K=1$) may be averaged to estimate the DPSD.
The $K$ independent spectra may be obtained from different sources (i.e., different azimuth, range or dwell, or generated from simulations), whereas the $K'$ bootstrapped spectra are obtained from the same source.
The $N_i$ independent realizations are performed to assess the statistical characteristics of the DPSD.
The weather-like I/Q time-series signals are simulated using the procedure described in \citet{Zrnic1975a}, extended to dual-polarized signals with the procedure described by \citet{Galati1995}, as presented in \citet{Torres2003a}.
This procedure is the same for all the methods evaluated hereafter.
%A complete description of the simulation methodology is presented in Chapter \ref{chap_4}.
The estimates using the IID bootstrap are shown in Figure \ref{fig:iid}.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/iid_boot.eps}
\caption{Mean of $N_i=1000$ DPSD estimates obtained with IID bootstrap (blue) and the ideal estimates (black): $s\hat{\mathrm{SNR}}_\mathrm{H}$ (top), $\szdrh$ (middle), and $\srhvh$ (bottom).}
\label{fig:iid}
\end{figure}
To assess the quality of the estimates, the bootstrap DPSD estimates are compared to ideal DPSD estimates obtained with a large number of independent spectra ($K=100$).
It can be seen that the mean of the spectral estimates of the IID bootstrap tend to the average radar variables.
For this single signal case, the values of the DPSD for frequency components with high enough SNR should be approximately equivalent to the bulk radar variables.
However, by assuming no data dependence in the time domain, the information about the distribution of radial velocities cannot be captured.
As previously mentioned, the underlying phenomena controlling the I/Q time-series signals are correlated, thus a dependent data bootstrap must be used.

\subsection{Non-Overlapping Block Bootstrap (NBB)}
The NBB was defined in Section \ref{sec:nbb_def}, where the I/Q time-series signal can be expressed as,
\be \mathbf{V}_\mathrm{H,V} =\{ \mathbf{B}_0, \mathbf{B}_1, \ldots, \mathbf{B}_{N-1} , \ldots, V_\mathrm{H,V}(M-1)\} ,\ee
where
\be \mathbf{B}_j = \{ V_\mathrm{H,V}( jL ) , \ldots ,  V_\mathrm{H,V}[ (j+1)L-1 ] \}, \ee
%\be \mathcal{B} = \{ \mathbf{B}_0, \mathbf{B}_1, \ldots, \mathbf{B}_{N-1} \} \ee
%$\mathcal{B}$ defines all the possible blocks for resampling (Figure \ref{fig:nbb}).
and an I/Q pseudo-realization $\mathbf{V}'_\mathrm{H,V}$ is obtained as
\be \mathbf{V}'_\mathrm{H,V} =\{ \mathbf{B}'_0, \mathbf{B}'_1, \ldots \mathbf{B}'_{N-1} \} ,\ee 
with
\be P(\mathbf{B}'_i = \mathbf{B}_j | \mathcal{B}) = \frac{1}{N} \ee
for $0 \leq i,j \leq N-1$.
$\mathcal{B}$ is the set of all available blocks, $M$ is the length of the signal, $L$ is the block length, and $N= \left \lfloor \frac{M}{L} \right \rfloor$ is the number of blocks.
Figure \ref{fig:nbb_iq} illustrates this process.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/bs_iq_example_nbb.eps}
\caption{Obtaining I/Q pseudo-realizations with the NBB.}
\label{fig:nbb_iq}
\end{figure}
If $N=ML$, there will be exactly $M$ samples in the pseudo-realization.
Otherwise, the last $M-NL$ samples from $\mathbf{V}$ will not have a block assigned and therefore will not be part of any pseudo-realization.
In this case, an additional block is drawn from $\mathcal{B}$ and the pseudo-realization is truncated to $M$ samples.

The DPSD is estimated using the same procedure mentioned in the previous section, but using the NBB scheme for different block lengths (Figure \ref{fig:nbb_dpsd}).
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/nbb_dpsd.eps}
\caption{Mean of $N_i=1000$ DPSD estimates obtained with NBB for $L=M/2$ (blue), $L=M/3$ (green), $L=M/4$ (red), and the ideal estimates (black): $s\hat{\mathrm{SNR}}_\mathrm{H}$ (top), $\szdrh$ (middle), and $\srhvh$ (bottom).}
\label{fig:nbb_dpsd}
\end{figure}
The effect of spectral leakage can be seen in all estimates and it is such that even a very aggressive data window function as the Blackman-Nuttall (SLL -98 dB) is not enough to mitigate it.
This example shows the tradeoff between the selection of the bootstrap block length and the spectral leakage introduced by disrupting the sample-time coherence.
Methods to mitigate this effect will be discussed later, but for now, it is assumed that the discontinuities between blocks are not accounted for in any of the methods.

\subsection{Moving Block Bootstrap (MBB)}
The definition of the MBB was presented in Section \ref{sec:mbb_def}, and the implementation on the I/Q time-series signal is as follows.
Let the set of available blocks for resampling be
\be \mathcal{B} = \{ \mathbf{B}_0, \mathbf{B}_1, \ldots, \mathbf{B}_{N'-1} \} ,\ee
where each block is given by 
\be \mathbf{B}_j = \{ V_\mathrm{H,V}(j) , \ldots ,  V_\mathrm{H,V}(j + L - 1) \} \ee
for $0\leq j \leq N'-1$, where $N' = M-L+1$ is the number of available blocks.
An I/Q pseudo-realization $\mathbf{V}'_\mathrm{H,V}$ is obtained as
\be \mathbf{V}'_\mathrm{H,V}  =\{ \mathbf{B}'_0, \mathbf{B}'_1, \ldots \mathbf{B}'_{N-1} \} ,\ee  
where
\be P(\mathbf{B}'_i = \mathbf{B}_j | \mathcal{B}) = \frac{1}{N'} ,\ee 
for $0 \leq i \leq N-1$, where $N = \left\lceil \frac{M}{L} \right\rceil$ is the number of blocks to be drawn from $\mathcal{B}$ for each pseudo-realization, and $L$ is the block length.
Figure \ref{fig:nbb_iq} illustrates this process.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/bs_iq_example_mbb.eps}
\caption{Obtaining I/Q pseudo-realizations with the MBB.}
\label{fig:mbb_iq}
\end{figure}

The DPSD estimated using MBB for different block lengths is shown in Figure \ref{fig:mbb_dpsd}
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/mbb_dpsd.eps}
\caption{Mean of $N_i=1000$ DPSD estimates obtained with MBB for $L=0.9M$ (blue), $L=0.8M$ (green), $L=0.7M$ (red), $L=0.6M$ (teal) and the ideal estimates (black): $s\hat{\mathrm{SNR}}_\mathrm{H}$ (top), $\szdrh$ (middle), and $\srhvh$ (bottom).}
\label{fig:mbb_dpsd}
\end{figure}
For block lengths of $0.9M$ and $0.8M$ there is almost no spectral leakage, while it is noticeable for $0.7M$ and even higher for $0.6M$.
However, for larger block lengths the DPSD estimates are significantly biased.
The MBB estimates for block lengths of $0.8M$ seem to follow the shape of $\srhv$ better, but $\szdrh$ still has large biases.
As the block length is reduced, it can be seen that the error is also reduced.
Moreover, in the limit where the estimates are unbiased, the block length is $1$, which is the case of the IID bootstrap. 
The tradeoff between the block length and the quality of the estimates is thus evident.
%A proper discussion will be presented in another section.

\subsection{Circular Block Bootstrap (CBB)}
The CBB was defined in Section \ref{sec:cbb_def}.
The implementation on the I/Q time-series signal is similar to the MBB, though the blocks are now allowed to wrap data points circularly.
With the CBB, an I/Q pseudo-realization $\mathbf{V}'_\mathrm{H,V}$ is obtained as
\be \mathbf{V}'_\mathrm{H,V}  =\{ \mathbf{B}'_0, \mathbf{B}'_1, \ldots \mathbf{B}'_{M-1} \}, \ee  
where
\be P(\mathbf{B}'_i = \mathbf{B}_j | \mathcal{B}) = \frac{1}{M}, \ee 
for $0 \leq i \leq N-1$, where $N = \left\lceil \frac{M}{L} \right\rceil$ is the number of blocks to be drawn from $\mathcal{B}$ for each pseudo-realization, and $L$ is the block length (Figure \ref{fig:cbb_iq}).
The set of available blocks $\mathcal{B}$ is 
\be \mathcal{B} = \{ \mathbf{B}_0, \mathbf{B}_1, \ldots, \mathbf{B}_{M-1} \} ,\ee
with
\be \mathbf{B}_j = \{ V_\mathrm{H,V}(\langle j \rangle_M) , \ldots ,  V_\mathrm{H,V}( \langle j + L - 1 \rangle_M ) \} , \ee
for $0\leq j \leq M-1$, and there are as many blocks available as $M$ samples in the I/Q signal.
The $\langle \cdot \rangle_x$ is the modulo-$x$ operator.
\begin{figure}[!htbp]
\centering
\includegraphics[width=1\linewidth]{figures/chap3/bs_iq_example_cbb.eps}
\caption{Obtaining I/Q pseudo-realizations with the CBB.}
\label{fig:cbb_iq}
\end{figure}

With the CBB, results similar to those of the MBB are expected.
However, because circular wrapping is allowed, higher spectral leakage is also introduced by the lack of coherence between the samples in blocks closer to the ends of the original sample set.
The DPSD estimated using CBB for different block lengths is shown in Figure \ref{fig:cbb_dpsd}.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/cbb_dpsd.eps}
\caption{Mean of $N_i=1000$ DPSD estimates obtained with CBB for $L=0.9M$ (blue), $L=0.8M$ (green), $L=0.7M$ (red), $L=0.6M$ (teal) and the ideal estimates (black): $s\hat{\mathrm{SNR}}_\mathrm{H}$ (top), $\szdrh$ (middle), and $\srhvh$ (bottom).}
\label{fig:cbb_dpsd}
\end{figure}
Once again the counteracting effects of block length and spectral leakage can be seen.
An advantage of the CBB is that it has a larger number of different blocks available to resample, which yields higher variability or ``randomness" of the pseudo-realizations.
Additionally, it is observed that when compared to the MBB for a same block length, the biases for frequency components with high $s\mathrm{SNR}$ are lower for the CBB although the spectral leakage is higher, as expected.
It is important to note that a very aggressive window is used for the analysis, yet the spectral leakage is significant.
This leakage contaminates the frequency components of noise, and as such, there is a possibility to wrongfully conclude that there is significant signal where there is none.

\subsection{Stationary Block Bootstrap (SBB)}
Following the SBB scheme defined in Section \ref{sec:sbb_def}, the pseudo-realizations are obtained as
\be \mathbf{V}'_\mathrm{H,V} = \{ \mathbf{B}(i_0,j_0),  \mathbf{B}(i_1,j_1), \ldots \}, \ee
where
\be \mathbf{B}(i,j) = \{ V_\mathrm{H,V}(i), \ldots, V_\mathrm{H,V}(i+j-1) \}, \ee
$i$ is uniformly distributed in $U(0, j-1)$ and defines the starting index of the block, and $j$ has a geometric distribution with parameter $p$ and defines the length of the block.
For each pseudo-realization, blocks are continuously drawn until the length of the pseudo-realization is equal to or greater than the original sequence.
Then, it is truncated to have a length equal to $M$.
The estimates for values of $p=0.1$, $0.01$ and $0.01$ are shown in Figure \ref{fig:sbb_dpsd}.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/sbb_dpsd.eps}
\caption{Mean of $N_i=1000$ DPSD estimates obtained with CBB for $p=0.1$ (blue), $p=0.01$ (green), $p=0.001$ (red), and the ideal estimates (black): $s\hat{\mathrm{SNR}}_\mathrm{H}$ (top), $\szdrh$ (middle), and $\srhvh$ (bottom).}
\label{fig:sbb_dpsd}
\end{figure}
The estimates show results similar to those obtained with the MBB and CBB.
The biases for frequency components with high $s\mathrm{SNR}$ are lower compared to the MBB and CBB.
However, due to the randomness of the block length of each pseudo-realization, it is harder to control the spectral leakage on the estimates.

\subsection{Summary of Block Bootstrap Methods}
It has been observed that smaller block lengths can help increase the variability which in turn decreases the errors for the spectral coefficients corresponding to the significant components of the signal, but also decreases the sample-time coherence, which in turn increases the spectral leakage.
Herein, the significant components of the signal are the spectral coefficients for which the $s\mathrm{SNR}$ is higher than $20$ dB.
%This is a general rule that applies to all block bootstrap methods.
The NBB has the lowest variability for any block length because the blocks are not allowed to overlap.
This means that individual blocks from the NBB block population are less correlated.
The MBB has a higher variability than the NBB, but lower than both the CBB and SBB.
Due to this, the errors in the MBB estimates are higher, but the main advantage of the MBB is that its block population best preserves the signal coherence and thus has the lowest spectral leakage.
The CBB has lower errors than the MBB, but higher spectral leakage due to the characteristics of its blocking scheme.
The SBB has the lowest errors because of the high variability in the pseudo-realizations.
However, it is also the worst at preserving signal coherence and the has highest spectral leakage.
With this exercise, it has been understood that the spectral leakage and the error of the block bootstrap methods depend mainly on the block length $L$ and the block arrangement.
For the NBB, MBB and CBB methods, the block length $L$ can be adjusted to obtain better estimator quality, while for the SBB this cannot be done.
The NBB exhibits side-lobe contamination in the spectra even with an aggressive window such as the Blackman-Nuttall window; therefore, it can be concluded that this method is not appropriate for spectral estimation of weather signals.
The SBB performs poorly for detecting features in the spectra, so it is also discarded from the analysis.
Ideally, the spectral leakage characteristics of the MBB and the error characteristics of the CBB would be desired characteristics of a useful block bootstrap method.
Since the CBB is an extension of the MBB, the advantages of both methods can be exploited to obtain optimal results both in terms of  spectral leakage and statistical errors.
Such strategies and considerations are discussed in the following section.

\section{Considerations in the Design of the Bootstrap DPSD Estimator} \label{sec:considerations}
The selection of the best block bootstrap method is based on parameters that result in an ``optimal" quality of the spectral estimates. 
The choice of such parameters is a compromise between a number of factors, such as spectral leakage, statistical errors of estimates, and computational cost.
Before selecting the block bootstrap method and analyzing the performance with different parameters, the criteria to determine the quality of the DPSD estimates must be established.

\subsection{Indicators of the Quality of Spectral Estimates}
\subsubsection{Spectral Leakage}
As mentioned before, the random processes controlling the changes in the I and Q components are correlated.
Provided that the pulse repetition time is small enough, the phasor of the I/Q signal follows a somewhat ordered path \citep{Doviak2006}.
%By bootstrapping the I/Q time-series signals the confidence interval of the spectral estimates can be retrieved, which can be used to obtain a better approximation of the real spectra.
Block bootstrap methods randomly resample blocks of data to generate a pseudo-realization, but by doing so the sample-time coherence is disrupted (i.e., discontinuities are added to the phasor's path) increasing spectral leakage.
As mentioned previously, this leakage contaminates the frequency components of noise or other weaker signals, which can lead to the wrong conclusion when performing spectral analysis.
Therefore, mitigating spectral leakage is important in terms of the quality of the spectral estimates.

%\subsubsubsection{Discontinuities}
\paragraph{Discontinuities}
The two main sources of spectral leakage for block bootstrapping are block discontinuities and circular discontinuities.
A block discontinuity is a transition that occurs between the last (or first) sample of a block and the first (or last) sample of an adjacent block.
Clearly, with an increasing number of blocks, there is also an increasing number of discontinuities.
A circular discontinuity is a transition that occurs between the last (or first) sample of the original sequence and the first (or last) sample of the same sequence when the signal is extended periodically.
Although these incoherencies added by the random resampling of blocks can provide additional variability, they should be kept to a minimum to reduce spectral leakage.
To achieve this, it is important to keep the bootstrap block length as large as possible to have the least amount of discontinuities in the signal while also keeping the population of different blocks as large as possible to have more variability.
This tradeoff must be considered in the design to ensure estimates provide useful information with minimal errors.

\paragraph{Data Windowing}
The choice of the data windowing function to be applied to the original sequence as well as to the pseudo-realizations has an impact on spectral leakage.
As shown by \citet{Harris1978}, the data windowing function smoothes a time-series signal as samples near the edges are attenuated, which improves the spectral estimates because the extraneous frequency components due to the discontinuity at the edges are attenuated.
However, doing so decreases the ability to resolve different peaks in the spectra due to the ``smearing" of the power of frequency components into adjacent components, i.e., the main-lobe of the spectral window is wider.
If a proper windowing function is not applied to the signal, additional spectral leakage may occur.
By de-emphasizing the incoherent pulses, the spectral leakage can be reduced.

\subsubsection{Statistical Errors}
In addition to the other indicators of quality of spectral estimates, the accuracy and precision of the estimates obtained with the block bootstrap methods must be assessed to select the method that produces the best estimates.
As it was observed from the preliminary results of implementing the block bootstrap methods to estimate the DPSD, a major tradeoff between the statistical errors and the spectral leakage is expected.
Clearly, the aim is to develop an estimator with small errors and little spectral leakage.
However, that is not feasible and a compromise between both must be reached.
A reasonable compromise is to determine a block bootstrap method that results in minimal leakage while still providing acceptable error levels for practical purposes.
For this, the quality of the estimates of the different block bootstrap methods is analyzed in terms of spectral leakage and statistical errors.

\subsubsection{Computational Cost}
Because the bootstrap methods inherently depend on generating a large number of pseudo-realizations from one measurement, their computational cost is also of interest since we would like the Bootstrap DPSD estimator to be viable in an operational environment (i.e., a real-time implementation should be feasible).
This is especially important because there are typically hundreds of thousands of I/Q signals per radar scan, which can quickly increase the computational cost by a few orders of magnitude depending on the complexity of the blocking scheme and the number of pseudo-realizations generated.
While this is noted, a proper analysis is beyond the scope of this work, though a discussion about the resampling scheme will be presented later.

\subsection{Strategies to Improve the Quality of the Spectral Estimates}
Although bootstrapping weather radar I/Q signals can enable spectral estimation when the available number of independent spectra is limited, improvements must be made to obtain the best possible results.
The strategies described in this section work independently to optimize the quality of the spectral estimates; thus, the core of the Bootstrap DPSD estimator is a combination of the best block bootstrap method and these pre-processing strategies.

\subsubsection{Signal Coherency Correction}
To mitigate the effects of discontinuities, coherency corrections can be applied to the blocks of a pseudo-realization.
These corrections consist in multiplying an adjacent block in a pseudo-realization by a complex scalar.
\citet{Borowska2016} presents a phase correction technique for phased array radar signals similar to the one presented in this work, though the strategy presented here is better suited for polarimetric radar I/Q time-series signals.
In general, these correction factors have three degrees of freedom: the block used as reference, the polarimetric channel used as reference, and the type of correction.

The propagation of the correction factor can be either forward or backward depending on block used as a reference to compute the correction factor.
That is
\be C_i = \begin{cases} 
\frac{ B'_{i-1}(L-1) }{ B'_{i}(0) } = C_i ^{+} &\text{for forward correction} \\ 
 \frac{ B'_{i+1}(0) }{ B'_{i}(L-1) } = C_i ^{-} &\text{for backward correction} , 
\end{cases} \ee
where $i$ indicates the index of the block (within the pseudo-realization) to be corrected. 
A superscript by the correction factor indicates whether it is a forward or backward correction, and the lack of a superscript indicates that the propagation direction is irrelevant.
The first element of the block has an index of $0$, while the last element of the block has an index of $L-1$.
The propagation of the correction indicates whether the $i$-th block is corrected using the last element of the previous block (forward propagation), or the first element of the next block (backward propagation). 
When applying a forward (or backward) correction, the blocks must be drawn with an additional data point (i.e., a block length of $L+1$) because the first (or last) pulse in the corrected block will be equal to the last (or first) pulse of the contiguous block and, as such, should be eliminated.
Because the I/Q signals are random processes, the choice of forward or backward correction should provide similar results since the correction factors are random.
However, from our experience, it is best to keep an uncorrected block around the center of the signal as the reference to correct the other blocks because the data windowing functions have maximum weight around $M/2$.
In such case, if the center block is the $i$-th block, the blocks with indices greater than $i$ should be corrected using forward correction, and those with indices lesser than $i$ should be corrected using backward correction.

The type of correction could be based on magnitude, phase, or both.
The corrected pseudo-realization then becomes
\be \mathbf{V}' = \left \{ C'_0\mathbf{B}'_0, C'_1 \mathbf{B}'_1, \ldots, C'_{N-1} \mathbf{B}'_{N-1} \right \}  = \{ C' _i \mathbf{B}'_i \} \ee
\be C'_i = \begin{cases} |C_i| & \text{for magnitude correction}\\ e^{j\arg (C_i)}& \text{for phase correction} \\C_i&  \text{for magnitude and phase correction}, \end{cases} \ee
%for $0 \leq i \leq N-1$, where $N$ is the number of blocks in the pseudo-realization.
Figure \ref{fig:mp_correction} shows the propagation and the complex constant correction methods that can be applied to a signal.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/magpha_correction.eps}
\caption{Correcting sample-time coherency on a periodically extended complex sinusoid using different methods.}
\label{fig:mp_correction}
\end{figure}
Here, the signal used to illustrate this concept is a sum of complex sinusoids of different frequencies.
The block of the uncorrected signal is in the center, while the signal is periodically extended to the left and to the right.
The first row shows the extended signal with no corrections, the second row is corrected by magnitude, the third row is corrected by phase, and the last row is corrected by magnitude and phase.
Additionally, the blocks on the left side are backward-corrected, whereas the blocks on the right side are forward-corrected.
Clearly, the magnitude-and-phase-corrected signal appears to have a smoother transition for both its real and imaginary parts while the magnitude- and phase-corrected signals are inferior in correcting the discontinuities.
It can be inferred that the magnitude and phase correction provides the best discontinuity mitigation, and thus, is the preferred type of correction.

Although it is not explicitly stated in the formulas, the correction factors depend on other correction factors applied to other blocks prior to the correction of the $i$-th block. 
For example, consider the case where all blocks of a particular pseudo-realization are corrected with forward correction in magnitude and phase.
If $\tilde{\mathbf{B}}'_1 = C'_1 \mathbf{B}'_1$ is the first block of the pseudo-realization after a coherency correction, the second block should be computed as $\tilde{\mathbf{B}}'_2 = C'_1 C'_2 \mathbf{B}'_2$.
This is done to consistently propagate the coherence throughout the blocks, otherwise there would still exist discontinuities in the sequence.
From this, it should be clear that the corrections are propagated along the pseudo-realization, which is one reason why it is important to keep the number of blocks as small as possible when corrections are applied.

So far, two degrees of freedom of the correction factor have been presented.
The last one deals with which polarimetric channels are used as a reference to apply the corrections.
That is, the correction factor may be applied to both channels based on the H- or V-channel, to each channel individually, or to both channels simultaneously using a combined correction factor.
A correction based on the H-channel is
\be \begin{aligned}
 \mathbf{V}'_\mathrm{H} = \left \{ C'_{\mathrm{H},0} \mathbf{B}'_{\mathrm{H},0}, C'_{\mathrm{H},1} \mathbf{B}'_{\mathrm{H},1}, \ldots, C'_{\mathrm{H},N-1} \mathbf{B}'_{\mathrm{H},N-1} \right \} = \{ C'_{\mathrm{H},i} \mathbf{B}'_{\mathrm{H},i} \} \\
 \mathbf{V}'_\mathrm{V} = \left \{ C'_{\mathrm{H},0} \mathbf{B}'_{\mathrm{V},0}, C'_{\mathrm{H},1} \mathbf{B}'_{\mathrm{V},1}, \ldots, C'_{\mathrm{H},N-1} \mathbf{B}'_{\mathrm{V},N-1} \right \}  = \{ C'_{\mathrm{H},i} \mathbf{B}'_{\mathrm{V},i} \} , \end{aligned}
\ee
where the $C'_{\mathrm{H},i}$ are correction factors based on blocks of the H-channel.
Similarly, a correction based on the V-channel is
\be \begin{aligned} \mathbf{V}'_\mathrm{H} = \left \{ C'_{\mathrm{V},0} \mathbf{B}'_{\mathrm{H},0}, C'_{\mathrm{V},1} \mathbf{B}'_{\mathrm{H},1}, \ldots, C'_{\mathrm{V},N-1} \mathbf{B}'_{\mathrm{H},N-1} \right \}   = \{ C'_{\mathrm{V},i} \mathbf{B}'_{\mathrm{H},i} \} \\
 \mathbf{V}'_\mathrm{V} = \left \{ C'_{\mathrm{V},0} \mathbf{B}'_{\mathrm{V},0}, C'_{\mathrm{V},1} \mathbf{B}'_{\mathrm{V},1}, \ldots, C'_{\mathrm{V},N-1} \mathbf{B}'_{\mathrm{V},N-1} \right \}   = \{ C'_{\mathrm{V},i} \mathbf{B}'_{\mathrm{V},i} \}  , \end{aligned} \ee
where the $C'_{\mathrm{V},i}$ are correction factors based on blocks of the V-channel.
The choice of either the H- or V-channel based correction should yield similar results for the DPSD, although the spectral leakage will be lower for H-channel if using the H-channel based correction, and vice versa.
The individual correction for each channel is
\be \begin{aligned} \mathbf{V}'_\mathrm{H} = \left \{ C'_{\mathrm{H},0} \mathbf{B}'_{\mathrm{H},0}, C'_{\mathrm{H},1} \mathbf{B}'_{\mathrm{H},1}, \ldots, C'_{\mathrm{H},N-1} \mathbf{B}'_{\mathrm{H},N-1} \right \} = \{ C'_{\mathrm{H},i} \mathbf{B}'_{\mathrm{H},i} \} \\
 \mathbf{V}'_\mathrm{V} = \left \{ C'_{\mathrm{V},0} \mathbf{B}'_{\mathrm{V},0}, C'_{\mathrm{V},1} \mathbf{B}'_{\mathrm{V},1}, \ldots, C'_{\mathrm{V},N-1} \mathbf{B}'_{\mathrm{V},N-1} \right \}   = \{ C'_{\mathrm{V},i} \mathbf{B}'_{\mathrm{V},i} \} .\end{aligned} \ee
And the combined correction is
\be \begin{aligned} \mathbf{V}'_\mathrm{H} = \left \{ C'_{\mathrm{X},0} \mathbf{B}'_{\mathrm{H},0}, C'_{\mathrm{X},1} \mathbf{B}'_{\mathrm{H},1}, \ldots, C'_{\mathrm{X},N-1} \mathbf{B}'_{\mathrm{H},N-1} \right \}   = \{ C'_{\mathrm{X},i} \mathbf{B}'_{\mathrm{H},i} \} \\
 \mathbf{V}'_\mathrm{V} = \left \{ C'_{\mathrm{X},0} \mathbf{B}'_{\mathrm{V},0}, C'_{\mathrm{X},1} \mathbf{B}'_{\mathrm{V},1}, \ldots, C'_{\mathrm{X},N-1} \mathbf{B}'_{\mathrm{V},N-1} \right \}   = \{ C'_{\mathrm{X},i} \mathbf{B}'_{\mathrm{V},i} \} \end{aligned}  \ee
Two variants of the combined correction were considered:
\be C_{\mathrm{X},i} = \begin{cases} \frac{B'_{\mathrm{H}, i-1}(L-1) + B'_{\mathrm{V}, i-1}(L-1) }{ B'_{\mathrm{H}, i}(0) + B'_{\mathrm{V}, i}(0)  } = C_{\mathrm{X},i} ^+ & \text{for forward correction} \\
\frac{B'_{\mathrm{H}, i+1}(0) + B'_{\mathrm{V}, i+1}(0) }{ B'_{\mathrm{H}, i}(L-1) + B'_{\mathrm{V}, i}(L-1)  } = C_{\mathrm{X},i} ^- & \text{for backward correction} , \end{cases} \label{eq:cxbad} \ee
and
\be C_{\mathrm{X},i} = \begin{cases} \frac{1}{2} \left[ \frac{B'_{\mathrm{H}, i-1}(L-1)}{B'_{\mathrm{H}, i}(0)} + \frac{B'_{\mathrm{V}, i-1}(L-1)}{B'_{\mathrm{V}, i}(0) } \right] = C_{\mathrm{X},i} ^+ & \text{for forward correction} \\
\frac{1}{2} \left[ \frac{B'_{\mathrm{H}, i+1}(0)}{ B'_{\mathrm{H}, i}(L-1) } + \frac{B'_{\mathrm{V}, i+1}(0) }{ B'_{\mathrm{V}, i}(L-1)  } \right] = C_{\mathrm{X},i} ^- & \text{for backward correction} . \end{cases} \ee
However, the latter variant is selected because it is more robust for a wide range of $\phi_\mathrm{DP}$ values (e.g., the denominator in \ref{eq:cxbad} would become small or tend to zero as the value of $\phi_\mathrm{DP}$ gets close to $180^\circ$, which can produce unexpected results).
It should be noted that the selection of the channel used for the correction is independent of the choice of the type or propagation.
Figure \ref{fig:chan_correction} is provided to illustrate how the different correction factors can be applied to the H- and V-channel signals.
\begin{figure}
\centering
\includegraphics[width=0.9\linewidth]{figures/chap3/chan_correction.eps}
\caption{Methods to apply the correction factors on the H- and V-channel signals.}
\label{fig:chan_correction}
\end{figure}

To summarize the correction methods, the following set of expressions are provided for a clearer interpretation on how these corrections are applied. 
The general form of the corrected signal becomes
\be \begin{aligned} 
\mathbf{V}'_\mathrm{H} = \{ C'' _{\mathrm{H},i} \mathbf{B}'_{\mathrm{H},i} \} \\
\mathbf{V}'_\mathrm{V} = \{ C'' _{\mathrm{V},i} \mathbf{B}'_{\mathrm{V},i} \} 
, \end{aligned} \ee
for $0\leq i \leq N-1$, where $i$ is the block index within the pseudo-realization, $N$ is the number of blocks in the pseudo-realization, and the correction factors $C'' _{\mathrm{H},i}$ and $C'' _{\mathrm{V},i}$ depend on how the channels are corrected as
\be C'' _{\mathrm{H},i},C'' _{\mathrm{V},i} = \begin{cases} 
	C' _{\mathrm{H},i},C' _{\mathrm{H},i} &  \text{for H-channel based correction} \\
	C' _{\mathrm{V},i},C' _{\mathrm{V},i} & \text{for V-channel based correction} \\
	C' _{\mathrm{H},i},C' _{\mathrm{V},i} & \text{for individual correction} \\
	C' _{\mathrm{X},i},C' _{\mathrm{X},i} &  \text{for combined correction}
. \end{cases} \label{eq:corfac_dp} \ee
Then, the factors $C' _{\mathrm{H},i}$, $C' _{\mathrm{V},i}$, and $C' _{\mathrm{X},i}$ depend on the complex constant as
\be C' _{x,i} = \begin{cases} |C_{x,i}| & \text{for magnitude correction}\\ e^{j\arg (C_{x,i})}& \text{for phase correction} \\C_{x,i}&  \text{for magnitude and phase correction} ,\end{cases} \ee
for $x=$ H,V or X, and where the $C_{x,i}$ factors depend on whether the previous or next block is used as a reference, as
\be  \begin{aligned}
C_{\mathrm{H},i} = \begin{cases} 
\frac{ B'_{\mathrm{H},i-1}(L-1) }{ B'_{\mathrm{H},i}(0) } = C_{\mathrm{H},i} ^{+} &\text{for forward correction} \\ 
 \frac{ B'_{\mathrm{H},i+1}(0) }{ B'_{\mathrm{H},i}(L-1) } = C_{\mathrm{H},i} ^{-} &\text{for backward correction} 
, \end{cases} \\
C_{\mathrm{V},i} = \begin{cases} 
\frac{ B'_{\mathrm{V},i-1}(L-1) }{ B'_{\mathrm{V},i}(0) } = C_{\mathrm{V},i} ^{+} &\text{for forward correction} \\ 
 \frac{ B'_{\mathrm{V},i+1}(0) }{ B'_{\mathrm{V},i}(L-1) } = C_{\mathrm{V},i} ^{-} &\text{for backward correction} 
, \end{cases} \\
C_{\mathrm{X},i} = \begin{cases} \frac{1}{2} \left[ \frac{B'_{\mathrm{H}, i-1}(L-1)}{B'_{\mathrm{H}, i}(0)} + \frac{B'_{\mathrm{V}, i-1}(L-1)}{B'_{\mathrm{V}, i}(0) } \right] = C_{\mathrm{X},i} ^+ & \text{for forward correction} \\
\frac{1}{2} \left[ \frac{B'_{\mathrm{H}, i+1}(0)}{ B'_{\mathrm{H}, i}(L-1) } + \frac{B'_{\mathrm{V}, i+1}(0) }{ B'_{\mathrm{V}, i}(L-1)  } \right] = C_{\mathrm{X},i} ^- & \text{for backward correction} . \end{cases}
\end{aligned} \label{eq:corfac_0} \ee
Because some concepts need to be introduced in order to simplify the analysis, the selection of the best channel correction method will be deferred to later.

\subsubsection{Extended Signal}
Although the corrected data provides an added variability to the pseudo-realizations, it is desired to keep the number of original (uncorrected) samples as high as possible, regardless of the coherence correction method used, so that the estimated spectra are as close to the true spectra as possible.
The MBB is suited for this task because it introduces fewer discontinuities when the block length is large.
However, by doing so, the block population available for resampling gets reduced significantly.
The CBB can overcome this limitation since it allows large blocks and has good variability, but the spectral leakage is very high.
In the limit when the block length is $M$ (i.e., same as the original sequence), the signal is as coherent as possible.
Though by assigning blocks of length $M$, the MBB would only have one block in the population.
However, the CBB has exactly $M$ blocks of length $M$ available because the signal is extended periodically.
If the samples of the periodically extended signal were coherent, the CBB would have the best performance of all methods.
Thus, the coherence correction can be applied to the periodic extensions of the original sequence to improve the quality of the estimates.
Therefore, to benefit from both methods, an extended I/Q signal with coherency corrections is proposed.
With this idea, the extended signal is kept as coherent as possible by applying the correction to the periodic extension, while also keeping the original signal unaltered.
To create this extended signal (Figure \ref{fig:extended_signal}), three instances of the signal $\mathbf{V}$ are concatenated.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/extended_signal.eps}
\caption{Generation of the extended signal $X_\mathrm{H,V}(m)$. The left and right extensions of the original sequence are denoted by the superscripts L and R, respectively; while the channel of polarization is indicated by the H and V subscripts, for horizontal and vertical channels, respectively.}
\label{fig:extended_signal}
\end{figure}
Then, the left extension is backward-corrected and the right extension is forward-corrected.
That is,
\be \mathbf{X}_\mathrm{H,V} = \{ \mathbf{V}^{\mathrm{L}}_\mathrm{H,V}, \mathbf{V}_\mathrm{H,V}, \mathbf{V}^{\mathrm{R}}_\mathrm{H,V} \}, \label{eq:Xm}\ee
where
\be \begin{aligned} 
\mathbf{V}^{\mathrm{L}}_\mathrm{H,V} = C''^-_\mathrm{H,V}  \{ V_\mathrm{H,V}(0), \ldots, V_\mathrm{H,V}(M-2) \},  \\
\mathbf{V}^{\mathrm{R}}_\mathrm{H,V} = C''^+_\mathrm{H,V}  \{ V_\mathrm{H,V}(1), \ldots, V_\mathrm{H,V}(M-1) \},
 \end{aligned}  \label{eq:vlr}\ee
and the correction factors $C''^-_\mathrm{H,V}$ and $C''^+_\mathrm{H,V}$ can be obtained with (\ref{eq:corfac_dp})-(\ref{eq:corfac_0}) depending on the type of correction. 
Here, the superscripts indicate backward and forward corrections specifically, while the type of correction and the polarimetric channel(s) used as reference(s) are unspecified.
It should also be noted that the extensions are of length $M-1$ due to the fact that, with the correction, the first (or last) element of a block is matched to the last (or first) element of the adjacent block.
The length of the extended signal then becomes $M' = 3M - 2$, and the block population for the extended signal is defined as
\be \mathcal{B} = \{ \mathbf{B}_0, \mathbf{B}_1, \ldots, \mathbf{B}_{N'-1} \}, \ee
where
\be \mathbf{B}_j = \{ X_\mathrm{H,V}( j ) , \ldots ,  X_\mathrm{H,V}( j + M - 1 ) \}. \ee
With the use of an extended signal, the block population increases to $N' = 2M - 2$, where the block length is now $M$.
However, a pseudo-realization obtained by resampling from this set of blocks is still prone to spectral leakage, albeit lesser than the case with no coherency correction applied.

\subsubsection{Ratio of Original Samples}
To further improve the quality of the bootstrapped estimates, a careful selection of the blocks can be made.
As mentioned previously, by keeping a larger number of original samples in a pseudo-realization, it is guaranteed that the spectral estimates will be less biased.
Logically, pseudo-realizations with more discontinuities (corrected or not) are more likely to contain more spectral leakage.
Hence, instead of using all the blocks of the extended signal, it is possible to select the set of blocks that contain more original samples than corrected samples.
The concept of the ratio of corrected samples in a block (or pseudo-realization) is introduced as
\be r = \frac{N_\mathrm{corrected}}{M} ,\ee
where $N_\mathrm{corrected}$ is the number of corrected samples in the block (or pseudo-realization).
In this sense, the original sequence has a ratio $r=0$, while the $(M/2)$ and $(3M/2)$-th blocks have a ratio close to $r=0.5$, and the first and the last blocks have a ratio close to $r=1$.
For example, by establishing a threshold for blocks that have a maximum ratio of $r=0.5$, it can be guaranteed that all the blocks in the reduced population will have at least 50\% of the original samples (Figure \ref{fig:block_ratio}).
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/block_ratio.eps}
\caption{Applying data windows on blocks with different ratio.}
\label{fig:block_ratio}
\end{figure}
Clearly, this strategy can additionally reduce the spectral leakage due to the incoherencies remaining in the periodic extension of the original signal.

While the selection of $r$ can be made arbitrarily with tradeoffs in spectral leakage and error statistics, its automatic determination is more convenient.
Due to how the blocks are arranged, $r$ depends on the data windowing function used.
The spectral leakage is reduced with this strategy mainly because the data points near the discontinuities are de-emphasized, while more weight is placed on the original samples.
Essentially, it is comparable to applying a sliding data window on a signal, while allowing the signal to be periodically extended.
A similar approach is presented by \citet{Barbe2010}, where Welch's method is extended by allowing segments to overlap circularly.
The segments that overlap circularly use windows of different length to deal with the circular discontinuities.
However, this method still loses frequency resolution if the length of the segments is not equal to the number of samples.
Moreover, the use of different windows for the segments results in poor SLL for spectral estimation.
Another approach where different data windows are used is with the Thomson's multi-taper method (MTM) for PSD estimation \citep{Thomson1982}. 
In this method, a discrete number of orthogonal windows are used to capture different characteristics of the spectra, then the partial PSD estimates are averaged together.
However, this method is not suitable for the spectral analysis of weather signals because there is no control over the SLL, and therefore estimates show can be severely biased.
The use of an extended signal along with a maximum ratio of corrected samples allows extracting spectral information from data points that are attenuated by the data windowing function without introducing additional biases in terms of both spectral leakage and statistical errors.
As mentioned before, the ratio that maximizes the amount of information depends on the data windowing function used in the analysis.
A good compromise between spectral leakage and statistical errors is obtained by limiting the selection to blocks with a maximum ratio of:
\be r = \frac{1-\sqrt{\frac{1}{M}\sum\limits^{M-1}_{m=0} |d(m)|^2}}{2} = \frac{1-\sqrt{\alpha}}{2} \label{eq:minratio} ,\ee
where $\alpha$ is the mean power of the data windowing function, and the amplitude $\sqrt{\alpha}/2$ is related to the number of points in the data window that contribute the most to the spectral estimate.
It is assumed that the data windows are normalized to a maximum amplitude of 1.
With this definition, $r$ is bounded between $0.5$, where half of the samples are corrected, and $0$, where none of the samples are corrected.
It can be easily seen that more aggressive windows have a lower window power, thus leading to higher maximum ratio and vice versa.
For example, if only the original sequence is used to estimate the PSD, with the use of a data window, part of the information is lost due to the end samples being severely attenuated.
The setting of a maximum ratio threshold for the resampling of blocks can be thought of as if a data window is shifted and applied on the extended signal to recover information that was attenuated by the conventional application of the data window.
Less aggressive windows allow more power from the corrected samples, therefore their corresponding maximum ratio is lower than for more aggressive windows.
In other words, the amplitude of the data windowing function determines the number of corrected samples, as weighted by the window,  that can be present on either end of a sequence before the spectral leakage becomes significant.

\subsubsection{Power Correction}
An additional step is taken to correct the power of the I/Q pseudo-realizations. 
Since the coherence correction factors scale both magnitude and phase of corrected samples, a power correction must be applied to preserve the power of the original signal. 
Let $\hat{P}_\mathrm{H,V}$ be the estimated average power of the H- and V-channels of the original signal, and $\hat{P}'_\mathrm{H,V}$ be the average power of the H- and V-channels of a I/Q pseudo-realization. 
The average power of each pseudo-realization is matched to the average power of the original signal;
mathematically,
\begin{equation}\label{eq:PowCor}
V''_\mathrm{H,V}(m) = \sqrt{\frac{\hat{P}_\mathrm{H,V}}{\hat{P}'_\mathrm{H,V}}} V'_\mathrm{H,V}(m),
\end{equation}
so that
\be \hat{P}''_\mathrm{H,V} = \hat{P}_\mathrm{H,V}. \ee
%This ensures that the amplitudes of the samples of the pseudo-realizations are scaled proportionally to the original signal power.

\subsubsection{Bias Correction}
In general, the expected value of the DPSD estimate and the true value are related by
\be E[s\hat{Z}_\mathrm{dr} (k)] = s{Z}_\mathrm{dr}(k) + \mathrm{bias}[s\hat{Z}_\mathrm{dr}](k), \ee
and
\be E[\srhvh (k)] = \srhv(k) + \mathrm{bias}[\srhvh](k) \label{eq:bias}, \ee
where $sZ_\mathrm{dr}$ is the linear spectral differential reflectivity (i.e., $sZ_\mathrm{DR} = 10 \log sZ_\mathrm{dr}$).
The analytical expressions of the statistical errors of $\szdrh$ and $\srhvh$ were determined in \citet{Yu2012}.
The biases of the DPSD depend on the number of independent spectra $K$, the spectral SNR of the H- and V-channels, $s\mathrm{SNR}_\mathrm{H,V}$, and the real spectral correlation coefficient, $\srhv$.
The bias expressions derived from \citet{Yu2012} are: 
\be \frac{ \mathrm{bias}[s\hat{Z}_\mathrm{dr}](k) }{ s\hat{Z}_\mathrm{dr}(k) } = \frac{1}{\beta K }\left[1 - \srhvh^2(k) \right] \label{eq:bzdr} \ee
\be \frac{ \mathrm{bias}[\srhvh](k) }{ \srhvh(k) } = \frac{1}{\beta K }\left\{ \frac{ [1 - \srhvh^2(k) ]^2}{4\srhvh^2(k)} \right\}  \label{eq:brhv} ,\ee
where $\beta$ is introduced in this work as a error scaling factor accounting for the changes in the quality of bootstrapped estimates prior to the bias correction, and the estimates of the DPSDs are used rather than the true DPSDs. The $s\mathrm{SNR}$ terms are neglected because it is assumed that the SNR of the signal of interest is high (more than 20 dB).
The factor $\beta$ is determined by fitting different values of $r$ such that the error is minimized for all $K$, as
\be \beta = \begin{cases}
(1-r)^{-3.3} - 2(1-r)^{1.1} & \text{for }K=1 \\
(1-r)^{-4.5} - (1-r)^{-2.1} & \text{for }K>1
, \end{cases} \ee
where $r$ is the maximum ratio defined by (\ref{eq:minratio}), and $K$ the number of independent spectra.
By replacing (\ref{eq:bzdr}) and (\ref{eq:brhv}) in (\ref{eq:bias}), the following expressions are obtained
\be s\tilde{Z}_\mathrm{dr}(k) =  s\hat{Z}_\mathrm{dr} (k) \left\{ 1 - \frac{1}{\beta K }[1 - \srhvh^2(k) ] \right\} \label{eq:cor_szdr} \ee
\be s\tilde{\rho}_\mathrm{HV}(k) =  \srhvh(k) \left\{ 1 - \frac{1}{\beta K }\left\{ \frac{[1 - \srhvh^2(k) ]^2}{4\srhvh^2(k)} \right\}  \right\} ,\label{eq:cor_srhv}\ee
where the tilde denotes the bias corrected estimate to distinguish it from the true DPSD.

\section{Analysis and Optimal Selection of Strategies} \label{sec:analysis}
While some parameters have already been selected, there are others that still need further analysis before a final selection can be done.
So far it has been established that neither the MBB nor the CBB were by themselves the best block bootstrap methods, which is why the periodic extension of the signal and the maximum ratio strategies have been proposed.
In this sense, the block bootstrap method of choice is essentially the MBB because the blocks are not allowed to wrap around circularly on the extended signal.
However, the periodic extension idea was inspired by the CBB, even though the method is not a CBB under the formal definition because of the coherency corrections.
A more appropriate description is as a hybrid between the MBB and the CBB.
Hereafter, the proposed method is referred to as the hybrid block bootstrap (or HBB).
The parameters that give useful DPSD estimates under this method have been found to be $L=M$ (i.e., blocks of the same size as the original I/Q signal), while an optimal value for $K'$ is yet to be determined (and will be analyzed later in this section).
Moreover, with $L=M$ and also the limitation set by a maximum value of $r$, the block population $\mathcal{B}$ is further reduced to a point in which it is useful to compare a deterministic resampling scheme, as opposed to a random resampling scheme, i.e., using all or a subset of the possible pseudo-realizations in a systematic order, versus randomly selecting the pseudo-realizations.

The maximum ratio of corrected samples establishes a limit on the block population to control the leakage that is introduced in the spectral coefficients.
The selection of the maximum ratio can be made arbitrarily, although the value that maximizes the information in the set of pseudo-realizations can be computed and it depends on the data windowing function.

The parameter $\beta$ for correcting the biases of the DPSD estimates is dependent on the maximum ratio.
While simpler expressions for $\beta$ have been studied, the current form, which was obtained empirically (as explained in the previous section), gives satisfactory results.
Additionally, it can be automatically obtained for different data windows and number of independent spectra, which could further improve the quality of the Bootstrap DPSD estimates.

By having selected most parameters, the dimension of the problem is simplified.
A statistical assessment of the quality of the estimates using different polarimetric channel correction methods was mentioned previously and will be analyzed next, followed by a comparison of bootstrap and deterministic resampling schemes.

\subsection{Performance of Channel Correction Methods}
Two of the three aspects of the correction factor have been analyzed in previous sections.
An analysis of the polarization channel correction must be made to determine which one yields the best results.
Since it is difficult to objectively assess the performance of each channel correction method, the statistical quality of the estimates obtained with each variant must be analyzed (Figure \ref{fig:chan_correction_stats}).
To estimate the DPSD, the HBB is performed, as indicated earlier.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/chan_correction_stats_new.eps}
\caption{Statistics of the estimator quality for different channel correction methods as a function of the correlation coefficient ($\srhv$). The average noise floor as a measure of spectral leakage (left), the bias and stabdard deviation (SD) of $\szdrh$ (center), and the bias and standard deviation of $\srhvh$ (right) are shown. The simulation parameters are shown in Table \ref{tab:simpar1}, but with $\mathrm{SNR}=40$ dB and a von Hann window.}
\label{fig:chan_correction_stats}
\end{figure}
The simulation parameters used for this analysis are similar to those in Table \ref{tab:simpar1}, though the SNR is set to $40$ dB, and the data window used in the analysis is a von Hann window.
This is done in order to emphasize the spectral leakage as a quality indicator, since the errors (bias and standard deviation) are expected to be similar for all cases.
The spectral leakage is measured by the average noise floor in the PSD of both channels, while the accuracy and precision are measured by the statistical errors of the DPSD.
A proper description on how these statistical analyses are performed will be presented in Chapter \ref{chap_4}.
For now, it suffices to say that the average noise floor is computed as the average SNR of frequency components with $s\hat{\mathrm{SNR}} < 20$ dB over all independent realizations, while the bias and standard deviation are computed as the averages for frequency components with $s\hat{\mathrm{SNR}} \geq 20$ dB based on to the modeled value of $Z_\mathrm{DR}$ and $\rho_\mathrm{HV}$, and over all independent realizations.
The $\rho_\mathrm{HV}$ of the signal is varied to highlight the effects of the different polarization channel correction methods.
For high $\rho_\mathrm{HV}$, all methods should have comparable performance given that the H- and V-channel signals are very similar, and thus the correction factors based on any method should be roughly the same.
For lower $\rho_\mathrm{HV}$, the correction factors can no longer appropriately correct both channels simultaneously.
If the correction is based on either the H- or V-channel, either the V- or H-channel, respectively, will not be completely matched and will be affected with higher spectral leakage, as can be seen in the green curve of the top-left panel and red curve of the bottom-left panel of Figure \ref{fig:chan_correction_stats}.
If both channels are corrected individually, the lowest joint effect of leakage is observed.
However, this negatively affects $\szdrh$ because the correction factors are not necessarily equivalent, and also affects $\srhvh$ because it adds an artificial modification to the correlation coefficient.
Therefore, the method that reasonably corrects the spectral leakage for both channels while keeping good error statistics is the combined correction method.

\subsection{Bootstrap versus Deterministic Resampling}
Since many constraints have been put on the bootstrap resampling mechanism and the block population in order to obtain estimates of good quality, an analysis of the performance of bootstrap as compared to a deterministic resampling scheme is appropriate.
A deterministic resampling scheme consists in selecting all or a subset of the available pseudo-realizations in a systematic order in order to obtain averaged estimate. 
If the block population is small, it would be beneficial to implement this type of scheme because all of the available information would be used.
Figure \ref{fig:rand_v_det1} shows the quality indicators as a function of the number of pseudo-realizations generated.
The same simulation procedure as in the previous analysis is maintained.
\begin{figure}
\includegraphics[width=1\linewidth]{figures/chap3/rand_v_det_1_new.eps}
\caption{Estimator statistics for bootstrap and deterministic resampling schemes as a function of the number of pseudo-realizations ($K'$). The average noise floor as a measure of spectral leakage (left), the bias and SD of $\szdrh$ (center), and the bias and SD of $\srhvh$ (right) are shown. The deterministic resampling scheme in this scenario takes contiguous blocks with increasing ratio as the number of pseudo-realizations increases.}
\label{fig:rand_v_det1}
\end{figure}
The difference between the two estimates is the way the blocks are resampled.
The bootstrapped estimates are shown in blue, while estimates with deterministic sampling are shown in green.
The deterministic resampling scheme used for this case takes contiguous blocks with increasing ratio as the number of pseudo-realizations increases.
It can be observed that with this scheme, while the deterministic sampling starts with a lower leakage, its errors are very high.
In the limit where the maximum number of different blocks has been sampled (i.e., around $K' = 40$ for this particular data window), the deterministic resampling scheme converges with the bootstrap in error statistics.
However, it should be noted that bootstrap resampling achieves these error levels for smaller values of $K'$.
Assuming $10$ pseudo-realizations are good enough, the difference with a deterministic sampling scheme is four-fold.
The number of computations can rapidly increase for a large number of radar bins in an operational setting, though a more thorough analysis is required to estimate the actual difference in computational cost.

A different deterministic sampling scheme is analyzed next.
In this scheme, the blocks with  $r=0$ and $r=0.5$ are sampled first, and successive blocks are taken as the mid-point ratio of the previously sampled blocks, i.e. $r=0.25$ then $r=0.125$ and $r=0.375$, and so on.
With this modification, faster convergence is expected for the deterministic scheme.
\begin{figure}
\includegraphics[width=1\linewidth]{figures/chap3/rand_v_det_2_new.eps}
\caption{Estimator statistics for bootstrap and deterministic resampling as a function of the number of pseudo-realizations ($K'$). A modified deterministic resampling scheme is used. The blocks with $r=0.5$ are sampled first, and successive blocks are taken as the mid-point ratio of the previously sampled blocks.}
\label{fig:rand_v_det2}
\end{figure}
In this case, it is observed that the modified deterministic resampling scheme converges faster to the bootstrap, though a higher noise floor is present because the pseudo-realizations with $r\simeq 0.5$ have the least sample-time coherence from the population.
The standard deviation of estimates is lower, but the spectral leakage is higher (indicated by the average noise floor of $s\mathrm{SNR}_\mathrm{H}$ and $s\mathrm{SNR}_\mathrm{V}$ in Figure \ref{fig:rand_v_det2}).
Moreover, the $\szdrh$ bias is negative and higher than the bootstrap in absolute value.
It may be possible that a deterministic resampling scheme may be vulnerable to artifacts in the signal or sampling biases, and also depending on the scheme, systematic biases may be introduced (e.g., by choosing pseudo-realizations with $r\simeq 0.5$, as mentioned previously).
The bootstrap appears to be superior in both cases.
The statistical errors have a marginal improvement with increasing number of pseudo-realizations starting from approximately $K'=10$.
At this point, it would be a reasonable compromise to keep a lower number of pseudo-realizations to save computational cost.

\section{The Bootstrap DPSD Estimator} \label{sec:algorithm}
All of the strategies and design considerations as well as the selection of parameters have been taken into account to develop the Bootstrap DPSD estimator.
The concepts employed in the design of the Bootstrap DPSD estimator are explained with the necessary details in Sections  \ref{sec:bootwx}, \ref{sec:considerations}, and \ref{sec:analysis}.
Here, a concise description of the estimator is given.

The Bootstrap DPSD estimator is the result of the combining the bootstrap method with the DPSD estimator using averaged periodogram estimates. 
The basic idea is to generate bootstrap samples of the weather radar I/Q time-series signals, in order to construct a bootstrap aggregate of the PSD of the I/Q pseudo-realizations from which a DPSD estimate can be obtained. 
The signals of the H- and V-channels are bootstrapped as a pair in order to preserve the correlation between channels. 
Since the processes controlling the I/Q time series signals are correlated, the signals must be conditioned prior to the generation of bootstrapped pseudo-realizations, such that additional information can be extracted from each signal without degrading the quality of the estimates.
A block bootstrap method is then employed to generate a suitable number of pseudo-realizations from which the PSD are computed.
These are then input to the DPSD estimator which averages the PSDs to obtain a DPSD estimate.
Finally, a bias correction technique is applied to the DPSD estimate to obtain the Bootstrap DPSD estimate.
A diagram of the Bootstrap DPSD estimator is shown in Figure \ref{fig:iqbootdiag}.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap3/BDPSD_diagram.pdf}
\caption{Diagram of the Bootstrap DPSD estimator.}
\label{fig:iqbootdiag}
\end{figure}

The steps taken by the Bootstrap DPSD estimator are summarized as follows:
\begin{enumerate}
\item Construct the periodic extension of the time series signal $X_\mathrm{H,V}(m)$ as in (\ref{eq:Xm}) and (\ref{eq:vlr}), using the magnitude and phase correction, and the combined correction.
\item Compute the maximum ratio of corrected samples $r_\mathrm{max}$ as in (\ref{eq:minratio}) for the selected data windowing function.
\item Generate $K'=20$ pseudo-realizations of length $M$ from $X_\mathrm{H,V}(m)$ limiting the blocks to those with $r\leq r_\mathrm{max}$.
\item Apply the power correction to each pseudo-realization as in (\ref{eq:PowCor}).
\item Compute the PSD estimates using (\ref{eq:sh}) and (\ref{eq:sx}).
\item Average the PSD estimates obtained and compute the biased DPSD estimate using (\ref{eq:szdr}) and (\ref{eq:srhv}).
\item Compute the corrected DPSD estimates with (\ref{eq:cor_szdr}) and (\ref{eq:cor_srhv}) and the results from the previous step.
\end{enumerate}

The statistical performance of the Bootstrap DPSD estimator in various scenarios will be compared to conventional DPSD estimators and to the ideal DPSD in the next chapter.

\section{Summary} \label{sec:summary}
The bootstrap as a tool to measure the accuracy of statistical estimates was introduced.
Such methods are useful to make inferences about the unknown underlying distribution when the sample size is relatively small.
Different methods to bootstrap dependent data (e.g., correlated data) have been presented.
These methods are known as block bootstrap methods, and they may be applied to weather radar I/Q time-series signals, provided that the signals are conditioned in a way that additional information can be extracted from each pseudo-realization without degrading the quality of the estimates.
These considerations and strategies are at the core of the Bootstrap DPSD estimator.
Moreover, the appropriate selection of the parameters and methods plays a big role in the quality of the estimates.
With objective analyses of the most important parameters, the optimal set of values was defined for the Bootstrap DPSD estimator, which was summarized in algorithmic terms.
In the following chapter, the statistical performance of the proposed Bootstrap DPSD estimator will be analyzed under different scenarios.