The context of the problem, as established in Chapter \ref{chap_1}, would be incomplete without a detailed description of the concepts.
This chapter synthesizes the fundamental theory required to properly understand the problem. 
In Section \ref{polarimetric_radar}, the basics of the polarimetric weather radar time-series signal processing are presented. 
A description of the TDS and a hypothesis about their spectral polarimetric characteristics are presented in Section \ref{ch2tds}.
Section \ref{psd_estimation} presents methods to estimate the PSD and assesses the advantages and limitations of each method. 
Section \ref{dpsd_estimation} focuses on DPSD estimators in a similar manner.
This chapter concludes with a summary highlighting key aspects of each section.

\section{Polarimetric Radar Signal Processing}
\label{polarimetric_radar}
A complete review on the principles of Doppler weather radar can be found in \citet{Doviak2006}, and with additional details on polarimetric Doppler weather radars in \citet{Bringi2001}.
This section will focus on the signal processing aspects of a polarimetric radar.

\subsection{I/Q Time-Series Signals}
As the radar sends out a pulse of electromagnetic waves, energy is radiated to the scatterers along the path which reflect a portion of this energy back to the radar.
A weather signal is composed of the returns of a large number of hydrometeors in a sub-volume of the atmosphere.
Since the individual returns cannot be resolved, the weather signal is sampled at discrete range-time delays with the use of gating circuits which convert the analog backscattered returns to digital numbers.
The range gates define the approximate range of a resolution volume in space where the scatterers have the largest contribution to the weather signal sample.
The resolution volume is determined by a range weighting function related to the radar filter's amplitude response and the transmitted pulse envelope, and an angular weighting function related to the antenna's radiation pattern.
The echoes of each scatterer in the resolution volume add constructively or destructively (depending on their position within a wavelength) producing a composite complex voltage sample (Figure \ref{fig:iqsigs}).
For a particular range time delay $\tau_s$ and the $m$-th transmitted pulse this is
\be V(\tau_s, mT_s) = \sum \limits ^{N-1} _{i=0} A_i e^{-j \gamma_i} ,\ee
\be \gamma_i = \frac{4 \pi r_i}{\lambda} + \frac{4\pi v_i T_s}{\lambda} - \psi_{si} ,\ee
where $A_i$ is the amplitude of the $i$-th scatterer, $r_i$ is the range to the scatterer, $\lambda$ is the radar wavelength, $v_i$ is the radial velocity of the scatterer, $T_s$ is the pulse repetition time (PRT), $\psi_{si}$ is the phase shift imposed by the $i$-th scatterer, and $N$ is the total number of scatterers.
The amplitude of each scatterer is related to the radar constants and the backscattering cross section ($\sigma_\mathrm{b}$) by:
\be A_i \propto \sqrt{\frac{P_\mathrm{t} \sigma_\mathrm{b} }{ (4 \pi)^3 }} \frac{\lambda g }{ r^2} W_\mathrm{a} W_\mathrm{r} ,\ee
where $P_\mathrm{t}$ is the peak transmitted power, $\sigma_\mathrm{b}$ is the backscattering cross section of the scatterer, $g$ is the antenna gain, $W_\mathrm{a}$ is a weighting function due to the antenna pattern, and $W_\mathrm{r}$ is a range-weighting function.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap2/iq_gates_new.eps}
\caption{I/Q signals for different range gates as a function of sample time $T_s$, with $I(m)$ shown in blue and and $Q(m)$ shown in cyan.}
\label{fig:iqsigs}
\end{figure}
Alternatively, for a given range resolution volume,
\be V(m) = I(m) +jQ(m) ,\ee
\be I(m) = \sum \limits ^{N-1} _{i=0}  |A_i| \cos \gamma_i ,\ee
\be Q(m) = -\sum \limits ^{N-1} _{i=0}  |A_i| \sin \gamma_i ,\ee
represent the I/Q signal, for a dwell defined by $M$ consecutive pulses as $m = 0,\ldots,M-1$.
Hereafter, the range time delay and the pulse repetition time shall be ignored for the sake of simplicity.

\subsection{Polarimetric Radar Variables}
Typical radar moments that can be estimated from the I/Q signals are: reflectivity factor ($Z$), mean radial velocity ($\bar{v}_r$), and spectrum width ($\sigma_v$). 
Also, with the simultaneous transmission and reception of H- and V-channel polarized waves (with the WSR-88D), polarimetric radar variables can be estimated, providing more information about the scatterers in a radar resolution volume; these are: differential reflectivity ($Z_\mathrm{DR}$), magnitude of the lag-0 cross-correlation coefficient ($\rho_\mathrm{HV}$, herein referred to as the correlation coefficient), and differential phase ($\phi_\mathrm{DP}$).
The polarimetric I/Q signals shall be represented as
\be V_\mathrm{H,V} (m) = I_\mathrm{H,V} (m) + jQ_\mathrm{H,V} (m) \label{eq:iqsigs} .\ee
Weather radar signals can be considered to be stationary processes because no significant changes in the statistical properties occur during the short observation time of tens of milliseconds; also, they are assumed to be ergodic, i.e., the statistical properties of the ensemble can be inferred from sample-time averages.
Furthermore, since they are random in nature, to make inferences about the properties of the received signals, we must rely on estimating average statistical parameters.
We recall the autocorrelation function of a stationary process $V(m)$ which is defined as: \citep{Doviak2006}
\be R_\mathrm{H,V}(l) = E[ V_\mathrm{H,V}^*(m) V_\mathrm{H,V}(m+l) ] = \lim \limits _{M \rightarrow \infty} \frac{1}{M} \sum \limits ^{M-|l|-1} _{m=0} V_\mathrm{H,V}^*(m) V_\mathrm{H,V}(m+l), \ee
where $l$ is the time lag, with H or V indicating the polarization channel.
Because there is a finite number of samples, only an estimate of the autocorrelation function can be obtained as 
\be \hat{R}_\mathrm{H,V}(l) = \frac{1}{M-|l|} \sum \limits ^{M-|l|-1} _{m=0} V_\mathrm{H,V}^*(m) V_\mathrm{H,V}(m+l) ,\ee
for $|l| \leq M-1$, which is the unbiased autocorrelation function estimator.
An estimate of the average returned power can be obtained from the autocorrelation function at lag 0 as
\be \hat{R}_\mathrm{H,V}(0) = \frac{1}{M} \sum \limits ^{M-1} _{m=0} V_\mathrm{H,V}^*(m)V_\mathrm{H,V}(m) = \hat{P}_\mathrm{H,V} \label{eq:R0} .\ee 
By subtracting the noise power $N$, an estimate of the signal power can be obtained as
\be \hat{S}_\mathrm{H,V} = \hat{P}_\mathrm{H,V} - N_\mathrm{H,V} . \label{eq:Sigpow} \ee 
Also, an estimate of the autocorrelation function at lag 1 ($l=1$) is used for calculating estimates of the mean velocity and spectrum width using the autocovariance estimator (as will be shown later), which can be computed as
\be \hat{R}_\mathrm{H}(1) = \frac{1}{M-1} \sum \limits ^{M-2} _{m=0} V_\mathrm{H}^*(m)V_\mathrm{H}(m+1). \ee
Additionally, the cross-correlation between the H- and V-channel signals is represented as
\be \hat{R}_\mathrm{X}(l) = \frac{1}{M-|l|} \sum \limits ^{M-|l|-1} _{m=0} V_\mathrm{H}^*(m) V_\mathrm{V}(m+l) . \label{eq:RX}\ee
Brief descriptions of the polarimetric radar variables and their estimators are given next.

The general radar terminology for the backscattering cross section per unit volume is reflectivity ($\eta$), and it is important to relate this measure to factors that have meteorological significance.
If the raindrops are spherical and have small diameters compared to the radar wavelength, the radar reflectivity is related to the reflectivity factor ($Z$) by \citep{Doviak2006}
\be \eta = \frac{\pi^5}{\lambda^4}|K_\mathrm{w}|^2 Z ,\ee
where $K_\mathrm{w}$ is the dielectric factor for water.
The reflectivity factor measures of the power returned from the scatterers in a radar resolution volume for a given dwell and it is typically used, along with an assumed DSD, to estimate the rainfall rate and the liquid water content.
The reflectivity factor can be estimated as 
\be \hat{Z}_\mathrm{H} = 10 \log \hat{S}_\mathrm{H} + 20 \log r + 10 \log C ,\ee 
in dBZ, ignoring any losses, where $r$ is the range to the radar volume and $C$ is a constant that depends on characteristics of the radar system such as the wavelength, transmitted power, pulse-width, antenna gain, and antenna weighting patterns, as well as the dielectric factor of the hydrometeor, and other scaling factors.
The reflectivity factor is related to the physical properties of the scatterers, so, high returns are expected for higher concentration of particles and larger electrical sizes.

The mean radial velocity ($\bar{v}_r$) is the power-weighted average of the radial velocities of the scatterers in the radar resolution volume.
It can be estimated using the autocovariance estimator \citep{Doviak2006} as
\be \hat{v}_\mathrm{r} = -\frac{v_\mathrm{a}}{\pi} \arg (\hat{R}_\mathrm{H}(1)), \ee
where $v_a$ is the radar's maximum unambiguous velocity, the argument of $\hat{R}_\mathrm{H}(1)$ is in radians and $\hat{v}_r$ is in m/s.
The radar's maximum unambiguous velocity is the maximum range of radial velocity that can be observed by a Doppler radar and is related to the radar's wavelength and PRT by 
\be v_a = \pm \frac{\lambda}{4T_\mathrm{s}} .\ee

The spectrum width ($\sigma_v$) is a measure of the spread of the radial velocities about the mean, and it is a function of radar system parameters such as beamwidth, bandwidth, and pulse width, and the meteorological parameters that describe the distribution of hydrometeor density and velocity within the resolution volume \citep{Doviak2006}.
There are many mechanisms that can broaden the spectrum width.
Meteorological factors include turbulence and wind shear, which can cause relative motion of the scatterers in the volume.
Antenna rotation can additionally broaden the spectrum width as the location of the resolution volume changes from pulse to pulse causing a decorrelation of signal samples and subsequent increase in spectrum width.
Additionally, depending on the data window used in spectral processing, the spectrum width may also vary.
More aggressive windows have wider main lobe width, resulting in broader spectrum widths than with less aggressive windows.
The spectrum width can be simply estimated as
\be \hat{\sigma}_v = \frac{\sqrt{2} v_\mathrm{a} }{\pi} \sqrt{\left| \ln{ \left( \frac{ \hat{S} }{ | \hat{R}(1) | } \right) } \right| }  .\ee
However, NEXRAD uses a hybrid spectrum width estimator \citep{Meymaris2009} that uses the lags 0/1, 1/2, and 1/3 estimators for wide, medium, and narrow spectrum widths, respectively.

The differential reflectivity $Z_\mathrm{DR}$ is the ratio of the power in the H-channel to the V-channel and it provides an insight about the shape of the scatterers. 
It can be estimated as
\be \hat{Z}_\mathrm{DR} = 10 \log{\left( \frac{\hat{S}_\mathrm{H}}{\hat{S}_\mathrm{V}} \right)} .\ee
Small raindrops are typically spherical in shape and have a $Z_\mathrm{DR}$ of zero due to the power returns in both channels being equal.
Larger raindrops tend to become more oblate in shape, corresponding to the horizontal backscattered power being larger than the backscattered power in the vertical channel, and result in a positive $Z_\mathrm{DR}$.
Other hydrometeors, such as ice crystals, can have a prolate shape, producing the opposite effect and thus having a negative $Z_\mathrm{DR}$ \citep{Bringi2001}.
The $Z_\mathrm{DR}$ of debris can largely depend on its electrical properties, orientation, size and concentration and it can vary in a wide range \citep{Bringi1977}.

The correlation coefficient $\rho_\mathrm{HV}$ is a measure of how similar the H-channel voltage return is, in average, to the V-channel return. 
Therefore, it is a measure of the homogeneity of the returns from scatterers in the radar volume.
To estimate the correlation coefficient, the cross-correlation function between the H- and V-channels at lag $0$ is obtained from (\ref{eq:RX}), then
\be \hat{\rho}_\mathrm{HV} = \frac{|\hat{R}_\mathrm{X} (0)|}{\sqrt{\hat{S}_\mathrm{H} \hat{S}_\mathrm{V} } } . \ee
The returns in the H- and V-channels for any set of homogeneous scatterers will have a high $\rho_\mathrm{HV}$ (e.g., only raindrops).
When the voltage returns from different types of scatterers are mixed (e.g., hydrometeors and debris), the scattering can become less homogeneous, resulting in a lower $\rho_\mathrm{HV}$.

The differential phase $\phi_\mathrm{DP}$ is defined as the backscattered differential phase between the voltage returns from the H- and V-channels (from the scatterers in the resolution volume) and propagation differential phase (from the path the pulse travelled). 
It can be estimated using the cross-correlation function between the H- and V-channels at lag $0$ as
\be \hat{\phi}_\mathrm{DP} = \arg{ \left( \frac{\hat{R}_\mathrm{X} (0)}{\sqrt{\hat{S}_\mathrm{H} \hat{S}_\mathrm{V} } } \right) }.\ee
The differential phase is a measure of the resistance experienced by the electromagnetic waves as they travel through the atmosphere \citep{Kumjian2013}.
For example, waves traveling through a spherical scatterer would experience the same amount of resistance in both polarization channels ($\phi_\mathrm{DP} = 0^\circ$), while non-spherical scatterers would experience a positive difference if they are oblate ($\phi_\mathrm{DP} > 0 ^\circ$) and a negative difference if they are prolate ($\phi_\mathrm{DP} < 0 ^\circ$).
It also provides an insight about the concentration of the scatterers, as it is proportional to the number concentration of particles and tends to increase with increasing particle size \citep{Kumjian2013}.

\subsection{Spectral Analysis}
In addition to the radar moments and polarimetric radar variables, there are tools which provide information about the Doppler velocities of the scatterers within a radar resolution volume.
Hereafter, the term 'bulk' will be used to simply refer to both radar moments and polarimetric variables, and to distinguish them from the spectral variables.
The power spectral density (PSD) represents the power-weighted distribution of radial velocities of the scatterers within a resolution volume \citep{Doviak2006}.
The Wiener-Khinchin theorem states that if a process is wide-sense stationary, the autocorrelation function and the PSD form a Fourier transform pair \citep{Papoulis2002}.
The PSD is defined as
\be sS_\mathrm{H,V,X}(kf_0) \equiv \lim \limits _{M\rightarrow\infty} T_s \sum \limits ^{M-1} _{l=-(M-1)} R_\mathrm{H,V,X}(l) e^{-j2\pi k f_0 T_s l} ,\ee
where $f_0 = 1/ MT_s$, and $k = 0, \ldots, M-1$ is the spectral component index.
The preceding $s$ denotes spectral variables as opposed to bulk variables, and hereafter, $f_0$ will be omitted for simplicity.
Methods to estimate the PSD will be presented in Section \ref{psd_estimation}.

In a similar way, the DPSD represents the distribution of the polarimetric variables as a function of the radial velocities of the scatterers within a resolution volume.
For example, the spectral differential reflectivity is the extension of the differential reflectivity in the frequency domain.
It is defined as the ratio between the horizontal PSD $sS_\mathrm{H}$ and the vertical PSD $sS_\mathrm{V}$:
\be \szdr (k) = \frac{sS_\mathrm{H}(k)}{sS_\mathrm{V}(k)} .\ee
The spectral correlation coefficient and the spectral differential phase are likewise defined using the corresponding PSDs instead of the polarimetric radar variables.
This is,
\be \srhv(k) = \frac{| sS_\mathrm{X} (k) |}{ \sqrt{sS_\mathrm{H}(k) sS_\mathrm{V}(k) } } ,\ee
\be \sphidp(k) = \arg{ \left( \frac{ sS_\mathrm{X} (k) }{ \sqrt{sS_\mathrm{H}(k) sS_\mathrm{V}(k) } } \right)}  .\ee
The DPSD estimators will be presented in Section \ref{dpsd_estimation}.

To highlight the importance of spectral analysis in weather radar signal processing, typical signatures found in observations are briefly discussed.
% Spectra of weather and clutter, wx+gc
%\subsection{Spectra of weather and ground clutter signals}
Referring to Figure \ref{fig:typspectra}, the power spectrum of weather is Gaussian-like with relatively wide spectrum width and power that depends on the strength of the echoes and the radar sensitivity \citep{Zrnic1975a}.
For ground clutter, the power spectrum is also Gaussian-like with narrow spectrum width, relatively high power, and centered at $v_\mathrm{r} = 0$ m/s \citep{Doviak2006}.
If the resolution volume contains weather but is contaminated by ground clutter, the spectrum will be a combination of both.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap2/wx_gc.eps}
\caption{Power spectral densities of typical signals. Weather (WX, left), ground clutter (GC, center), weather and ground clutter (right). The PSDs correspond to returns from the H channel and are normalized to the noise power. Hence, the $y$-axis represents signal-to-noise ratio in the H channel ($s\mathrm{SNR}_\mathrm{H}$).}
\label{fig:typspectra}
\end{figure}

For a single signal in the resolution volume (i.e., either weather only or ground clutter only) and with the assumption of a Gaussian model, spectral analysis does not provide additional information about the scatterers since the estimates of the radar variables already reflect this.
However, for a composite signal, the radar variables are biased (e.g., they can be skewed towards the dominant signal).
%This means that if the ground clutter signal is equal to or higher than the weather signal in power, the radar measurement will be biased.
In the case where the ground clutter signal is completely dominant over the weather signal, information about the hydrometeors will be lost.
With spectral analysis, the signals can be separated with minimal biases introduced in the radar estimates.
Even in the case where the weather signal has a zero mean radial velocity and it is mixed with a ground clutter signal (Figure \ref{fig:0wxgc}), information about the weather signal can still be extracted because the ground clutter signal is spectrally narrower.
\begin{figure}
\centering
\includegraphics[width=0.6\linewidth]{figures/chap2/zero_wx_gc.eps}
\caption{Power spectrum of zero mean radial velocity weather mixed with ground clutter.}
\label{fig:0wxgc}
\end{figure}
The separation of a ground clutter signal from others is relatively easy because of the distinct characteristics seen in the power spectra, namely, the high power return, the narrow spectrum width and a mean radial velocity about zero.
However, when dealing with groups of scatterers move independently with different radial velocities, it becomes harder to separate them in the PSD alone because of the lack of information to distinguish them from one another, as is the case in resolution volumes with hydrometeors and debris.
In such cases, the DPSDs provide additional spectral information that could be used to identify groups of scatterers moving with different radial velocities that have distinct polarimetric characteristics.
It is hypothesized that these type of spectra could be found near TDSs, which is discussed next.

\section{Tornadic Debris Signatures}\label{ch2tds}
As mentioned previously, TDS are tornado-scale signatures with distinctive polarimetric characteristics collocated near a tornado vortex that are visible in radar observations (Figure \ref{fig:tds}) after a tornado touches ground \citep{Ryzhkov2002}. 
It is hypothesized that they are related to the lofting, ejection and centrifuging of debris by the tornado vortex \citep{Dowell2005}.
Since debris have polarimetric signatures that depend on their electrical size, shape, orientation and concentration \citep{Bodine2014b}, these signatures can be noticeably different from those of hydrometeors and can be used for polarimetric tornado detection \citep{Ryzhkov2005}.
In addition, they improve the detection of tornadoes \citep{Wang2015} when ground observation is limited or impossible (e.g., at night or during heavy rainfall) \citep{Kumjian2008}.
\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{figures/chap2/ppi05202013_200403.pdf}
\caption{Tornadic debris signatures in observations of 20 May 2013 Moore, OK tornado at 20:04:03 UTC with the PX-1000 radar. Reflectivity factor $Z_\mathrm{H}$ (top left), mean radial velocity $\bar{v}_\mathrm{r}$ (top right), differential reflectivity $Z_\mathrm{DR}$ (bottom left), correlation coefficient $\rho_\mathrm{HV}$ (bottom right).}
\label{fig:tds}
\end{figure}

Due to centrifuging in a tornado vortex, debris move at different velocities than hydrometeors.
Since Doppler radars measure the mean scatterers motion rather than the actual wind speed, significant biases can be introduced in the radar measurements \citep{Dowell2005}.
This can occur whether the debris signal is stronger than the weather signal or not.
To effectively assess the potential damage \citep{Bodine2013}, it is important to have accurate measurements to estimate the intensity of the tornado.
\citet{Ryzhkov2005} catalogued typical values of the polarimetric radar variables associated to the TDS as $Z > 45$ dBZ, $Z_\mathrm{DR} < 0.5$ dB, and $\rho_\mathrm{HV} < 0.8$, while \citet{Bodine2013} use more stringent sets of values for $Z$ and $\rho_\mathrm{HV}$: either $Z < 43$ dBZ and $\rho_\mathrm{HV} < 0.82$, or $Z < 51$ dBZ and $\rho_\mathrm{HV}< 0.72$.
With these criteria, it can be inferred that TDS are more likely to be seen in observations due to the debris signal dominating in a particular radar resolution volume \citep{Ryzhkov2005}, thus the mean radial velocities in volumes with TDS are more prone to centrifuging biases.
Furthermore, it is hypothesized that this error in measurement can be corrected if the velocity of the debris is known and can be separated from the velocity of hydrometeors, which are passive tracers of the wind. 
It is argued in this work that such information could be retrieved through spectral analysis.

% Spectra of Debris
%\subsection{TDS spectra}
With a Gaussian assumption, the PSD does not provide additional information if there is only one type of scatterer in the radar resolution volume.
In the case where there is more than one type of scatterer, discrimination might be possible provided that neither signal is completely obscuring the other signal.
However, the additional information provided by the DPSD is needed to appropriately classify the signals.
Because the DPSD shows the distribution of polarimetric variables as a function of the radial velocities, it may be possible to determine the type of the scatterers by correlating the DPSD values to known polarimetric characteristics of the different groups of scatterers.
However, identification of all scatterer types would be impossible.

As an example, a sketch of a bimodal DPSD for a combination of weather and debris signals is shown in Figure \ref{fig:sketch} \citep{Cheong2015}.
The representation of the spectral variables of debris may be exaggerated, but it illustrates a possible way to retrieve the actual velocity of the raindrops from DPSD estimates.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap2/debris_dpsd.eps}
\caption{Sketch of spectral polarimetric variables for a combination of weather and debris. Spectral power $sS_\mathrm{H}$ (left), spectral differential reflectivity $\szdr$ (center), spectral correlation coefficient $\srhv$ (right). Adapted from \citet{Cheong2015}.}
\label{fig:sketch}
\end{figure}
In this scenario, raindrops are moving faster, with a wider spectrum width, positive $\szdr$, and $\srhv$ high and close to $1$.
Conversely, debris is slower, with a narrower spectrum width, negative $\szdr$ and lower $\srhv$.
It is expected that the spectral variables of the scatterers would follow the polarimetric radar variables somewhat closely in value.
Without the DPSD, there would be a higher uncertainty in determining which spectral components of the PSD correspond to weather and which to debris.
However, by inspecting the spectral polarimetric variables it should be possible to separate the different scatterer types.
Although the debris signal is depicted as a Gaussian-shaped signal, it should be noted that such assumption is not realistic.
Nonetheless, it is a simple model that successfully illustrates the idea of the presence of signals with different polarimetric characteristics in the spectra.

%The estimation of DPSD will be presented in a further section.

\section{Estimation of Power Spectral Densities} \label{psd_estimation}
For PSD estimation, there are non-parametric and parametric methods. The former type makes no assumptions about the structure of the underlying phenomena, while the latter assumes their structure can be described with certain models by using different parameters.
The atmosphere is constantly changing, and only over 75\% of the observed power spectra are Gaussian shaped \citep{Janssen1985}, in addition to non-hydrometeor scatterers being present in the radar volume.
Thus, the uncertainty in determining the nature of the event would be high and extremely difficult to parameterize a model for each case. 
Therefore, non-parametric methods are better suited for this work.

Two common non-parametric methods to compute the PSD are the periodogram and the correlogram methods.
The periodogram consists in applying the discrete Fourier transform (DFT) to the I/Q time-series signal to obtain 
\be Z(k) = \sum \limits _{m=0} ^{M-1} V(m)e^{-j2\pi mk /M} ,\ee
where the signal $V$ corresponds to either the H- or V-channel. 
Then an estimate of the PSD can be obtained as
\be s\hat{S}(k) = \frac{|Z(k)|^2}{M} ,\ee
where $M$ is the total number of samples in the dwell, and $k$ is the frequency component index.
The correlogram method \citep{Stoica1997} consists in estimating the autocorrelation function of the I/Q time series and then applying the DFT.
An estimate of the PSD with the correlogram can be computed as
\be s\hat{S}(k) = \sum \limits _{l=-(M-1)} ^{(M-1)} \hat{R}(l) e^{-j2\pi lk /M} .\ee
%\begin{figure}[!htbp]
%\centering
%\includegraphics[width=0.5\linewidth]{figures/placeholder.png}
%\caption{Diagram of methods to obtain the PSD estimate.}
%\label{fig:psddiagram}
%\end{figure}
It can be shown that the periodogram and the correlogram estimators are equivalent under certain conditions \citep{Doviak2006}.
The method preferred herein is the periodogram because of its relative simplicity and low computational cost. 
In spite of its advantages, the periodogram comes with limitations regarding the accuracy and precision of the PSD estimates. 
Many variants of the periodogram \citep{Hayes1996, Stoica1997} have been proposed in the literature to address these issues, some of which will be described later in this chapter.

\subsection{Modified Periodogram}
Since the DFT assumes the underlying signal to be periodic, the discontinuities associated with the periodic extension of an I/Q dwell introduce frequency components that are not related to the underlying phenomenon. 
This effect is known as spectral leakage \citep{Harris1978}.
The modified periodogram is a generalized variation of the periodogram which accounts for this effect by applying a data window or tapering function to smooth out the discontinuities associated with the periodic extension of the time-series signal \citep{Hayes1996}. 
\be Z(k) = \sum \limits _{m=0} ^{M-1} d(m)V(m)e^{-j2\pi mk /M} \label{mperiodogram_1} ,\ee
where $d(m)$ is the data window function.
\begin{figure}
\begin{minipage}[c]{0.49\linewidth}
\centering
\includegraphics[width=1\linewidth]{figures/chap2/iq_period.eps}
\\(a)
\end{minipage}
\begin{minipage}[c]{0.49\linewidth}
\centering
\includegraphics[width=1\linewidth]{figures/chap2/iq_modperiod.eps}
\\(b)
\end{minipage}
\caption{Effects of data windowing in (a) the periodogram, and (b) the modified periodogram. Whereas the former uses a rectangular window (i.e., no window), the latter uses a tapered window (the von Hann window in this example). The in-phase component $I(m)$ shown in blue and the quadrature component $Q(m)$ shown in cyan.}
\label{fig:per_mper}
\end{figure}
Figure \ref{fig:per_mper} shows how the data windowing function is applied for the periodogram and for the modified periodogram.
It is important to note that by applying the taper the signal power is attenuated.
To compensate for this loss, the PSD is normalized by:
\be \alpha = \frac{1}{M} \sum \limits _{m=0} ^{M-1} |d(m)|^2 \label{mperiodogram_2}. \ee
Then, the modified periodogram estimate becomes
\be s\hat{S} (k) = \frac{1}{\alpha} \frac{|Z(k)|^2 }{ M } \label{mperiodogram_3} .\ee

The choice of the data windowing function can reduce the spectral leakage in exchange for frequency resolution \citep{Harris1978}.
Typical choice of windows include the Bartlett window (SLL -25 dB), the von Hann window (SLL -32 dB), the Hamming window (SLL -43 dB), and the Blackman window (SLL -58 dB), among others.
%Alternatively, this data window function may be applied in the frequency domain by convolving the DFT of the time-series signal with the DFT of the data window function as:
%\be Z(k) = \sum \limits _{m=0} ^{M-1} V(m)e^{-j2\pi mk /M} \ee
%\be D(k) = \sum \limits _{m=0} ^{M-1} d(m) e^{-j2\pi mk / M} \ee
%\be \hat{S}(k) = \frac{1}{\alpha} \frac{|D(k) \ast Z(k)|^2}{M} \ee
%though only the data windowing in time domain will be considered.
From here on, this method will be simply referred to as the periodogram.

\subsection{Bartlett's Method}
The method proposed by \citet{Bartlett1950}, consists in dividing the time-series signal into non-overlapping segments. 
The PSD of each segment is computed with the periodogram using a rectangular data window (Figure \ref{fig:iq_bart}), and then these partial estimates are averaged to reduce the variance of the frequency component estimates \citep{Stoica1997}.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap2/iq_bart.eps}
\caption{Bartlett's method makes non-overlapping segments of the I/Q signal, and computes the PSD for each segment using a rectangular data window. In this case, the length of each segment is $L=M/4$, and each segment is enumerated.}
\label{fig:iq_bart}
\end{figure}
The PSD estimate of a single dwell with Bartlett's method is computed as
\be Z_i(k) = \sum _{m=0} ^{L-1} V(m+iL) e^{-j2\pi mk/L} ,\ee
\be s\hat{S}_i(k) = \frac{|Z_i(k)|^2}{L} ,\ee
\be s\hat{S}(k) = \frac{1}{N} \sum \limits _{i=0} ^{N-1} s\hat{S}_i(k) ,\ee
for $ m = 0, \ldots, L-1$, $k = 0, \ldots, L-1$ and  $ i = 0, \ldots, N-1$, where $N= \lfloor M/L\rfloor$ is the number of segments, $L$ is the segment length, and $\lfloor \cdot \rfloor$ is the floor function.

Although Bartlett's method computes a smoother PSD estimate because the partial estimates are independent, it comes with its limitations and tradeoffs. 
The method, by definition, does not make use of data windowing functions and thus is less flexible compared to other methods. 
Since it reduces the number of samples to compute the partial PSD estimates, the frequency resolution gets degraded significantly, although the variance is also reduced by a factor of $N$.
This tradeoff between frequency resolution and variance is clear \citep{Stoica1997}.
Moreover the method has no control over the spectral leakage which makes it unsuitable for the dynamic range of weather signals.

\subsection{Welch's Method}
A variant of Bartlett's method was proposed by \citet{Welch1967} and it allows the segments to overlap (Figure \ref{fig:iq_welch}), as well as introducing an extra degree of freedom when selecting the data windowing function for each segment \citep{Stoica1997}. 
Bartlett's method is reduced to a specific case of the Welch's method when the data window function used is the rectangular window and the segments do not overlap.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap2/iq_welch.eps}
\caption{Welch's method makes overlapping segments of the I/Q signal, and computes the PSD for each segment using a particular data window. In this case, the length of each segment is $L=M/2$ with a $50\%$ overlap and the von Hann window. Each segment is enumerated.}
\label{fig:iq_welch}
\end{figure}
The PSD estimate of a single dwell with Welch's method is obtained as follows:
\be Z_i(k) = \sum _{m=0} ^{L-1} d(m) V[m+i(L-O)] e^{-j2\pi mk/L} ,\ee
\be \alpha = \frac{1}{L} \sum \limits _{m=0} ^{L-1} |d(m)|^2 ,\ee
\be s\hat{S}_i(k) = \frac{1}{\alpha} \frac{|Z_i(k)|^2}{L} ,\ee
\be s\hat{S}(k) = \frac{1}{N} \sum \limits _{i=0} ^{N-1} s\hat{S}_i(k) ,\ee
for $ m = 0, \ldots, L-1$, $k = 0, \ldots, L-1$ and  $ i = 0, \ldots, N-1$, where $N= \lfloor \frac{M-O}{L-O}\rfloor$ is the number of segments, $O$ the amount of segment overlap, and $L$ is the segment length.

Welch's method accounts for the limitations of Bartlett's method by adding more flexibility in the selection of the appropriate parameters that improves the statistical performance of the estimator. 
However, by allowing the overlap of the time-series signal segments, the partial PSD estimates are no longer independent even though the windowing function works to ``decorrelate" these estimates.
The resolution and the variance of Welch's method are better than Bartlett's method \citep{Stoica1997}.  
The performance of the method is analyzed in Chapter \ref{chap_4}.

\subsection{Daniell's Method}
The smoothed periodogram method proposed by \citet{Daniell1946} reduces the variance of the estimate in a different manner than in the methods described above. 
In this method, a moving average filter is applied to the ``raw" periodogram estimate \citep{Stoica1997} 
(See Figure \ref{fig:iq_daniell}).
To obtain the PSD estimate with Daniell's method, we first get an estimate using Equations (\ref{mperiodogram_1})-(\ref{mperiodogram_3}). Then, a moving average filter is applied as
\be s\hat{S}(k) = \frac{1}{2p+1} \sum \limits ^{\langle k+p\rangle_M} _{k' = \langle k-p \rangle_M} s\hat{S}(k') ,\ee
where $2p+1$ is the length of the filter, and $\langle \cdot \rangle_M$ is the modulo $M$ operator.
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap2/iq_daniell.eps}
\caption{Daniell's method smooths the raw PSD estimate. The length of the moving average filter used here is $2p+1 = 3$.}
\label{fig:iq_daniell}
\end{figure}
Clearly this method trades spectral resolution for reduction of variance, and it may lead to higher bias as it smoothes out the raw PSD. Additionally, when the frequency components of the signals of interest are too close, the ability to resolve them individually may be lost due to this smearing effect. 
The method can be extended to a more general version which reduces the spectral leakage by using a weighted moving average filter. 
However, such implementation is not analyzed in this study.
The performance of the method is analyzed in Chapter \ref{chap_4}.

\section{Estimation of Dual-Polarimetric Spectral Densities} \label{dpsd_estimation}
To estimate the DPSDs, the co- and cross-polar power spectra must be obtained first as
\be Z_\HV (k) =  \sum \limits _{m=0} ^{M-1} d(m)V_\HV(m) e^{-j2\pi mk/M } ,\ee
\be \label{eq:sh} s\hat{S}_\HV (k) = \frac{1}{\alpha} \frac{|Z_\HV (k) |^2}{M} ,\ee
\be \label{eq:sx} s\hat{S}_\mathrm{X} (k) = \frac{1}{\alpha} \frac{Z_\mathrm{H} (k) Z^*_\mathrm{V}}{M} ,\ee
To illustrate this, the expressions for the modified periodogram are shown, but it should be noted that the estimates $s\hat{S}_\mathrm{H}$, $s\hat{S}_\mathrm{V}$ and $s\hat{S}_\mathrm{X}$ may be obtained with any PSD estimator.

The polarimetric spectral estimates are obtained in a similar way to the polarimetric radar variables but using the PSD estimates as:
\be \label{eq:szdr}\szdrh(k) = \frac{\sumK s\hat{S}_\mathrm{H} ^{(i)} (k) }{ \sumK s\hat{S}_\mathrm{V}^{(i)} (k) } ,\ee
\be \label{eq:srhv} \srhvh(k) = \frac{\left| \sumK s\hat{S}_\mathrm{X}^{(i)} (k) \right|}{ \sqrt{\sumK s\hat{S}_\mathrm{H}^{(i)} (k) \sumK s\hat{S}_\mathrm{V}^{(i)} (k)} } ,\ee
\be \sphidph(k) = \arg{\left( \frac{\sumK s\hat{S}_\mathrm{X}^{(i)} (k)}{ \sqrt{\sumK s\hat{S}_\mathrm{H} ^{(i)}(k) \sumK s\hat{S}_\mathrm{V}^{(i)} (k)} } \right) } ,\ee
where $K$ is the number of independent spectra that are averaged to obtain useful DPSD estimates.
As previously stated, the spectral differential phase will not be part of the analyses because it does not convey as much information to discriminate hydrometeors from debris.

Operational weather radars perform a scan every few minutes, which yields one independent spectra ($K=1$) for each radar resolution volume. 
Using the periodogram estimator, for the spectral differential reflectivity, the frequency components of the PSD have significantly large variance, providing a poor $\szdr$ estimate. 
Additionally, it can be shown that the spectral correlation coefficient estimate fails to produce any useful results. 
By combining (\ref{eq:srhv}) and (\ref{eq:sh})-(\ref{eq:sx}) for $K=1$:
\begin{eqnarray}
s\hat{\rho}_\mathrm{HV}(k) = \frac{|s\hat{S}_\mathrm{X}(k)|}{\sqrt{s\hat{S}_\mathrm{H}(k)s\hat{S}_\mathrm{V}(k)}} = \frac{|Z_\mathrm{H}(k)Z^*_\mathrm{V}(k)|}{\sqrt{|Z_\mathrm{H}(k)|^2 |Z_\mathrm{V}(k)|^2}} \nonumber\\
=  \frac{|Z_\mathrm{H}(k)||Z^*_\mathrm{V}(k)|}{|Z_\mathrm{H}(k)||Z_\mathrm{V}(k)|} =  1
,\end{eqnarray}
which shows how the $\srhv$ estimate always equals 1.

This limitation can be overcome in different ways.
Perhaps the most obvious solution is to average more independent spectra to reduce the statistical errors of the estimates and to produce meaningful $s\hat{\rho}_\mathrm{HV}$.
To this end, data can be obtained from different sources such as adjacent radar resolution volumes in azimuth, range, or scans. 
Alternatively, the methods described in Section \ref{psd_estimation} provide estimates with better statistical performance than the periodogram estimator, and these could be used to obtain more reliable DPSD estimates with data from one dwell.
Previous works \citep{Bachmann2007a, Moisseev2009, Alku2015, Unal2009, Yanovsky2011, Dufournet2010, Spek2008, Moisseev2006, Moisseev2007, Yu2013, Yu2012} have used one or more of the previously mentioned methods. 

The DPSD can be used as tool for polarimetric spectral analysis, and it has been applied in the retrieval of microphysical parameters \citep{Yanovsky2011, Dufournet2010, Spek2008, Moisseev2006, Moisseev2007, Yu2014}, as well as for detecting and suppressing clutter signals \citep{Bachmann2007a, Moisseev2009, Alku2015, Unal2009}.
Ideally, it is desired to obtain reliable DPSD estimates out of an individual dwell for polarimetric spectral analysis to be both useful for interpretation and operationally feasible. 
Furthermore, it is preferred that resolution loss in any dimension is kept to the minimum in order to retrieve meaningful information for tornado observations.
The advantages and disadvantages of averaging in different dimensions is briefly discussed next.

\subsection{Range Averaging}
More spectra may be acquired from spatially correlated radar resolution volumes. 
If the range gates of a particular radial are chosen to be averaged (Figure \ref{fig:rg_averaging}), the range resolution is degraded by at least a factor of two (e.g., two adjacent range gates).
\begin{figure}
\centering
\includegraphics[width=0.45\linewidth]{figures/chap2/rg_average.eps}
\caption{Range averaging scheme. Radar bins of three consecutive resolution volumes along range are averaged in this example.}
\label{fig:rg_averaging}
\end{figure}
Depending on the range resolution of the radar, this could be quite significant as important spatial features of the weather event may be masked by averaging.
Additionally, since the range dimension of the resolution volume remains constant although the resolution volume increases in size with range due to the antenna's beamwidth and range may be oversampled, range averaging may be a less compromising option in regard to resolution loss. 

\subsection{Azimuth Averaging}
A different spatial averaging can be performed using independent spectra from adjacent radar resolution volumes in azimuth (Figure \ref{fig:az_averaging}).
Similar to range averaging, the azimuthal resolution is degraded at the very least by a factor of two (e.g., resolution volumes from two adjacent rays).
\begin{figure}
\centering
\includegraphics[width=0.45\linewidth]{figures/chap2/az_average.eps}
\caption{Azimuthal averaging scheme. Radar bins of three consecutive azimuths are averaged in this example.}
\label{fig:az_averaging}
\end{figure}
At farther ranges, the radar resolution volume gets wider in azimuth resulting in degraded spatial resolution, while for a constant range, resolution volumes will be similar in size.
Thus, averaging in azimuth may be favored over range averaging at closer ranges.
Additionally, if the azimuth is oversampled or if signal processing techniques such as the super-resolution \citep{Torres2006} are used, the impact of azimuth averaging may be acceptable.

\subsection{Scan-to-Scan Averaging}
Averaging spectra from different scans (Figure \ref{fig:dw_averaging}) can also be performed to obtain better DPSD estimates, provided the spectra are correlated in time (i.e., slow moving phenomena or short scan times).
\begin{figure}
\centering
\includegraphics[width=1\linewidth]{figures/chap2/dw_average.eps}
\caption{Dwell averaging scheme. Radar bins of three consecutive dwells are averaged in this example.}
\label{fig:dw_averaging}
\end{figure}
By doing so, it must be ensured that the observations are based on the same location in space for a given event.
However, to correctly capture the evolution of certain weather events, such as tornadoes, the time between consecutive scans for a given radar resolution volume must be considerably short.

\subsection{Frequency Averaging}
Some of the PSD estimation methods described in the previous section can provide decent DPSD estimates using data from a single dwell by reducing the statistical errors through averaging in the frequency domain.
Albeit useful in cases where no other type of averaging can be performed, these methods degrade the frequency resolution (Bartlett, Welch) or add additional spectral ``smearing" (Daniell).
As a consequence, the quality of these estimates are usually not enough to be used in any quantitative analysis.


% Add a section on differences between Radar Variable estimation and DPSD estimation --
% How to interpret DPSD and error statistics?

\section{Summary}
While the radar variables provide the ``big picture" of a weather event, measurements can be biased by the presence of dominant clutter (i.e., non-meteorological targets).
In such cases, a spectral analysis may be useful to infer the characteristics for groups of scatterers moving at different Doppler velocities.
A potential application for polarimetric spectral analysis is in radar resolution volumes with TDS, where the measurements are more prone to errors due to the weather signal being dominated by the debris signal.
The correction of these biases is important when estimating tornado intensity and damage potential through radar observations.
However, the quality of the DPSD estimates is limited by the number of independent observations that are needed to obtain accurate estimates.
Polarimetric spectral analysis of tornado observations is additionally hindered by the fact that tornadoes evolve fast in time in a relatively small spatial extent, and because debris are centrifuged in the tornadic vortex.
This means that good temporal and spatial resolutions are required to capture important features of the tornado, as well as a good frequency resolution to discriminate the velocities of hydrometeors and debris.
Currently available techniques for DPSD estimation are not sufficient to overcome these constraints.
In the following chapter, the Bootstrap DPSD estimator is introduced which is capable of estimating reliable DPSDs with minimal resolution loss. 
%In this study, an alternative way to estimate the DPSD from one dwell with minimal resolution loss is presented.