# pdSpecEst1D: Intrinsic wavelet HPD spectral estimation In pdSpecEst: An Analysis Toolbox for Hermitian Positive Definite Matrices

## Description

`pdSpecEst1D` calculates a (d,d)-dimensional HPD wavelet-denoised spectral matrix estimator by applying the following steps to an initial noisy HPD spectral estimate (obtained with e.g., `pdPgram`):

1. a forward intrinsic AI wavelet transform, with `WavTransf1D`,

2. (tree-structured) thresholding of the wavelet coefficients, with `pdCART`,

3. an inverse intrinsic AI wavelet transform, with `InvWavTransf1D`.

The complete estimation procedure is described in more detail in \insertCiteCvS17pdSpecEst or Chapter 3 of \insertCiteC18pdSpecEst.

## Usage

 ```1 2``` ```pdSpecEst1D(P, order = 5, metric = "Riemannian", alpha = 1, return_val = "f", ...) ```

## Arguments

 `P` a (d,d,m)-dimensional array of HPD matrices, corresponding to a sequence of (d,d)-dimensional HPD matrices of length m, with m = 2^J for some J > 0. `order` an odd integer larger or equal to 1 corresponding to the order of the intrinsic AI refinement scheme, defaults to `order = 5`. Note that if `order > 9`, the computational cost significantly increases as the wavelet transform no longer uses a fast wavelet refinement scheme based on pre-determined weights. `metric` the metric that the space of HPD matrices is equipped with. The default choice is `"Riemannian"`, but this can also be one of: `"logEuclidean"`, `"Cholesky"`, `"rootEuclidean"`, `"Euclidean"` or `"Riemannian-Rahman"`. See also the Details section below. `alpha` an optional tuning parameter in the wavelet thresholding procedure. The penalty (or sparsity) parameter in the tree-structured wavelet thresholding procedure in `pdCART` is set to `alpha` times the estimated universal threshold, defaults to `alpha = 1`. `return_val` an optional argument that specifies whether the denoised spectral estimator is returned or not. See the Details section below. `...` additional arguments for internal use.

## Details

The input array `P` corresponds to an initial noisy HPD spectral estimate of the (d,d)-dimensional spectral matrix at `m` different frequencies, with m = 2^J for some J > 0. This can be e.g., a multitaper HPD periodogram given as output by the function `pdPgram`.
`P` is transformed to the wavelet domain by the function `WavTransf1D`, which applies an intrinsic 1D AI wavelet transform based on a metric specified by the user. The noise is removed by tree-structured thresholding of the wavelet coefficients based on the trace of the whitened coefficients with `pdCART` by minimization of a complexity penalized residual sum of squares (CPRESS) criterion via the fast tree-pruning algorithm in \insertCiteD97pdSpecEst. The penalty or sparsity parameter in the optimization procedure is set equal to `alpha` times the universal threshold, where the noise variance of the traces of the whitened wavelet coefficients are determined from the finest wavelet scale. See \insertCiteCvS17pdSpecEst and Chapter 3 of \insertCiteC18pdSpecEst for further details.
The function computes the forward and inverse intrinsic AI wavelet transform in the space of HPD matrices equipped with one of the following metrics: (i) the affine-invariant Riemannian metric (default) as detailed in e.g., \insertCiteB09pdSpecEst[Chapter 6] or \insertCitePFA05pdSpecEst; (ii) the log-Euclidean metric, the Euclidean inner product between matrix logarithms; (iii) the Cholesky metric, the Euclidean inner product between Cholesky decompositions; (iv) the Euclidean metric; or (v) the root-Euclidean metric. The default choice of metric (affine-invariant Riemannian) satisfies several useful properties not shared by the other metrics, see \insertCiteCvS17pdSpecEst or \insertCiteC18pdSpecEst for more details. Note that this comes at the cost of increased computation time in comparison to one of the other metrics.
If `return_val = 'f'` the thresholded wavelet coefficients are transformed back to the frequency domain by the inverse intrinsic 1D AI wavelet transform via `InvWavTransf1D`, returning the wavelet-denoised HPD spectral estimate.

## Value

The function returns a list with the following five components:

 `f ` a (d,d,m)-dimensional array of HPD matrices, corresponding to the HPD wavelet-denoised estimate of the same resolution as the input array `P`. If `return_val != 'f'`, the inverse wavelet transform of the thresholded wavelet coefficients is not computed and `f` is set equal to `NULL`. `D ` the pyramid of threshold wavelet coefficients. This is a list of arrays, where each array contains the (d,d)-dimensional thresholded wavelet coefficients from the coarsest wavelet scale `j = 0` up to the finest wavelet scale `j = jmax`. `M0 ` a numeric array containing the midpoint(s) at the coarsest scale `j = 0` in the midpoint pyramid. `tree.weights ` a list of logical values specifying which coefficients to keep, with each list component corresponding to an individual wavelet scale starting from the coarsest wavelet scale `j = 0`. `D.raw ` the pyramid of non-thresholded wavelet coefficients in the same format as the component `\$D`.

## Note

The function does not check for positive definiteness of the input matrices, and (depending on the specified metric) may fail if matrices are close to being singular.

## References

\insertAllCited

`pdPgram`, `WavTransf1D`, `InvWavTransf1D`, `pdCART`
 ```1 2``` ```P <- rExamples1D(2^8, example = "bumps")\$P f <- pdSpecEst1D(P) ```