Denoise a polarized spectrum via trend filtering

Source: R/denoise_spectrum.R

The denoise_spectrum() function uses quadratic trend filtering to optimally denoise each of the \(IQU\) Stokes parameters of a polarized spectrum, in turn also leading to denoised estimates for the normalized Stokes parameters \(Q/I\) and \(U/I\). Setting compute_uncertainties = TRUE generates a bootstrap ensemble of denoised spectra for each Stokes parameter, which allows variability bands to be computed for each denoised spectrum by then calling bands() on the denoise_spectrum() output.

denoise_spectrum(
  wavelength,
  flux,
  variance,
  mask,
  lambda = c("lambda_min", "lambda_1se"),
  compute_uncertainties = FALSE,
  break_at = 10L,
  min_pix_segment = 10L,
  mc_cores = parallel::detectCores() - 4,
  ...
)

Arguments

wavelength

Vector of wavelength measurements.
flux

Spectropolarimetric measurements, passed as a 3-column tibble, data frame, or matrix, with the columns corresponding to the \(IQU\) Stokes parameters, respectively.
mask

Pixel masks, in a tibble, data frame, or matrix with dimensions matching those of flux. Nonzero elements flag bad pixels.
compute_uncertainties

(Boolean) If TRUE, then bootstrap ensembles are created for each denoised spectrum via trendfiltering::bootstrap_trendfilter(). The ensembles are stored in the function output so that variability bands of any level can quickly be computed by calls to bands() without redundant overhead calculations (see examples). Defaults to compute_uncertainties = FALSE.
break_at

The minimum number of consecutively-masked spectral pixels that will trigger a break in the spectrum. Defaults to break_at = 10.
min_pix_segment

After the segmentation procedure is complete, the resulting segments are examined to ensure that each is sufficiently long for a non-trivial/well-defined denoising analysis. In particular, any segment that has fewer than min_pix_segment unmasked spectral pixels is discarded. Defaults to min_pix_segment = 10.
mc_cores

Multi-core computing using the parallel package: The number of cores to use. Defaults to the number of cores detected on the machine, minus 4.
...

(Optional) Named arguments to be passed to trendfiltering::sure_trendfilter() or trendfiltering::bootstrap_trendfilter(). The arguments x, y, weights, k, and algorithm cannot be overridden.
variances

Measurement variances, in a tibble, data frame, or matrix with dimensions matching those of flux.

Value

An object of class 'polarized_spectrum'. This is a list with the following elements:

n_segments

Number of segments the spectrum was broken into by break_spectrum().

data

The original data set, as a list of n_segments tibbles. Each tibble contains the observed wavelengths (with the union set of masked wavelengths removed), the Stokes \(IQU\) flux measurements, and the measurement variances.

denoised_signals

The set of denoised spectra, as a list of n_segments tibbles. Each tibble contains the wavelength evaluation grid for its respective segment and all of the denoised signal estimates: \(I\), \(Q\), \(U\), \(Q/I\), \(U/I\).

ensembles

If compute_uncertainties = TRUE, a list of bootstrap ensembles for the denoised \(I\), \(Q\), \(U\) parameters, respectively. Each ensemble is returned as an \(n \times B\) matrix. If compute_uncertainties = FALSE, this will return NULL.

I_analysis_summary

Technical summary of the denoising analysis of the \(I\) Stokes parameter.

Q_analysis_summary

Same as above, but for the \(Q\) Stokes parameter.

U_analysis_summary

Same as above, but for the \(U\) Stokes parameter.

References

  1. Politsch et al. (2020a). Trend filtering – I. A modern statistical tool for time-domain astronomy and astronomical spectroscopy. MNRAS, 492(3), p. 4005-4018.

  2. Politsch et al. (2020b). Trend Filtering – II. Denoising astronomical signals with varying degrees of smoothness. MNRAS, 492(3), p. 4019-4032.

See also

bands(), break_spectrum()

Examples

# Any SALT Observatory spectrum can be read into R from its FITS file
# using the "FITSio" R package, as below. In the example below, we will
# analyze the polarized spectrum of a SALT Wolf-Rayet star. For convenience,
# here we've stored the header and data units (HDU) from the FITS file into
# three variables (`sci`, `var`, and `bpm`) that are stored in an R data file
# within this package, so we can source them using a simple call to `data()`.

if (FALSE) {
install.packages("FITSio")
path_to_FITS_file <- "<your_local_path_to_FITS_file>"
file.name <- "WR006_c1_12345678_stokes.fits"
sci <- FITSio::readFITS(paste0(path_to_FITS_file, file.name), hdu = 1)
var <- FITSio::readFITS(paste0(path_to_FITS_file, file.name), hdu = 2)
bpm <- FITSio::readFITS(paste0(path_to_FITS_file, file.name), hdu = 4)
}

####

library(dplyr)
data(polarized_spectrum_WR_star)

wavelength <- seq(
  from = sci$axDat$crval[1],
  by = sci$axDat$cdelt[1],
  length = sci$axDat$len[1]
)

flux <- as_tibble(sci$imDat)
variance <- as_tibble(var$imDat) %>% select(1:3)
mask <- as_tibble(bpm$imDat)

spec_denoised <- denoise_spectrum(
  wavelength,
  flux,
  variance,
  mask,
  compute_uncertainties = TRUE
)