Resolution-Adaptive Bayesian Wavelet Denoising

Fits the resolution-adaptive Wendland-Semicircle Bayesian wavelet denoiser, shrinks detail coefficients by posterior means, and reconstructs the signal with the inverse discrete wavelet transform.

wswavelet(
  y,
  filter.number = 10L,
  family = "DaubExPhase",
  bc = "periodic",
  j0 = 0L,
  likelihood = c("gaussian", "laplace"),
  supplied_sd = NULL,
  spike_offset = 2,
  spike_gamma = 2.4,
  beta_quantile = 0.99,
  beta_quantile_type = 8L,
  beta_floor = 1e-8,
  omega_model = c("adaptive", "constant", "fixed"),
  fixed_omega = NULL,
  eta_start = c(0, 0),
  eta_bounds = c(-12, 12),
  laplace_rate = NULL,
  quadrature_n = 48L,
  use_exact_wendland = TRUE,
  mad_levels = 1L,
  return_wavelet = FALSE
)

Arguments

y: A finite numeric signal vector. Its length must be a dyadic integer, such as 128, 256, or 1024.
filter.number: The wavelet filter order passed to wd. The default is 10.
family: The wavelet family passed to wd. Common choices include "DaubExPhase" and "DaubLeAsymm".
bc: The boundary convention passed to wd. The default is "periodic".
j0: The coarsest detail level to shrink. Detail levels below j0 are not shrunk. The default is 0.
likelihood: The coefficientwise working likelihood. Use "gaussian" for the primary version or "laplace" for the robustness-oriented alternative.
supplied_sd: Optional known noise standard deviation. If NULL, the function estimates it from high-frequency wavelet coefficients using a robust MAD rule.
spike_offset: The offset $\ell$ in the level-dependent spike probability. The implementation permits values at least 1; the default is 2.
spike_gamma: The positive resolution-sparsity exponent $\gamma$. Larger values make the spike probability increase more rapidly with resolution.
beta_quantile: The probability $\kappa$ used in the empirical upper-quantile rule for the support scale $\beta_j$. It must lie strictly between 0 and 1; the default is 0.99.
beta_quantile_type: The integer quantile convention passed to quantile. The default is type 8.
beta_floor: A positive numerical floor used to prevent a zero support scale when the noise estimate and coefficients are numerically zero.
omega_model: How the Wendland weight is modeled across resolution levels. "adaptive" estimates an intercept and a resolution trend, "constant" estimates one common weight, and "fixed" uses fixed_omega.
fixed_omega: A fixed Wendland weight in $[0,1]$. Supplying this argument forces omega_model = "fixed"; 0 gives the semicircle-only endpoint and 1 gives the Wendland-only endpoint.
eta_start: Starting values for the empirical-Bayes logistic trend parameters $(\eta_0,\eta_1)$. For a constant mixture, only the first value is used.
eta_bounds: Common finite lower and upper bounds for the optimized logistic parameters. The default is c(-12, 12).
laplace_rate: Optional positive rate parameter for the Laplace working likelihood. If NULL, the function sets $\widehat{\lambda}=1/(2\widehat{\sigma}^2)$ once before optimization and holds it fixed. This argument has no effect for the Gaussian likelihood.
quadrature_n: Number of Gauss-Legendre nodes used for numerical integration. The default is 48. Larger values may improve accuracy but increase computation time.
use_exact_wendland: Logical value. For the Laplace likelihood, TRUE uses the finite-sum Wendland moment formula with a numerical fallback. If FALSE, numerical quadrature is used. For the Gaussian likelihood, numerical quadrature is used.
mad_levels: Number of finest detail levels pooled for the robust noise-scale estimate. The default is 1.
return_wavelet: Logical value. If TRUE, also returns the original and thresholded wavethresh wavelet objects.

Details

The prior for a detail coefficient at resolution level $j$ is $$ \theta_{j,k} \sim \pi_j\delta_0 + (1-\pi_j)\{\omega_j g_W(\theta;\beta_j) +(1-\omega_j)g_S(\theta;\beta_j)\}. $$ The spike probability is $$ \pi_j = 1-(j-j_0+\ell)^{-\gamma}. $$ For the adaptive model, the Wendland weight is $$ \omega_j(\eta)=\operatorname{logit}^{-1}(\eta_0+\eta_1 t_j), \qquad t_j=(j-j_0)/(J-1-j_0). $$ The support scale is selected from the observed detail coefficients using $$ \widehat{\beta}_j = \max\{\widehat{\sigma},Q_{\kappa}(|d_{j,k}|)\}, $$ followed by the positive numerical floor $$ \beta_{\min}=\epsilon_\beta s_\beta,\qquad s_\beta=\max\{1,\widehat{\sigma},\max_{j,k}|d_{j,k}|\},\qquad \widehat{\beta}_j\leftarrow\max\{\widehat{\beta}_j,\beta_{\min}\}, $$ where $\epsilon_\beta$ is beta_floor.

When omega_model = "constant", the fitted Wendland weight is $\operatorname{logit}^{-1}(\eta_0)$ at every level. When omega_model = "fixed", no weight optimization is performed.

For Gaussian errors, the coefficientwise likelihood is $$ L_N(d\mid\theta,\sigma^2)= (2\pi\sigma^2)^{-1/2} \exp\{-(d-\theta)^2/(2\sigma^2)\}. $$ For the Laplace working likelihood, the implementation uses $$ L_L(d\mid\theta,\lambda)=(a/2)\exp\{-a|d-\theta|\}, \qquad a=\sqrt{2\lambda}. $$ The empirical-Bayes parameters are fitted by bounded L-BFGS-B optimization of the joint log marginal likelihood. The posterior mean replaces each detail coefficient, while the scaling coefficients are retained without shrinkage.

Value

A list with the following principal components:

estimate: The reconstructed denoised signal.
sigma_hat: The estimated or supplied noise standard deviation.
laplace_rate_hat: The rate used for the Laplace working likelihood.
likelihood: The selected likelihood.
omega_model: The selected mixture-weight model.
eta_hat: The fitted logistic-trend parameters, padded to length two for the constant and fixed cases.
optimization: The optimization result or fixed-weight record.
level_summary: A data frame containing level-specific support scales, prior probabilities, fitted weights, average posterior component probabilities, and numerical diagnostics.
detail: A list containing the observed and fitted coefficients and posterior quantities for each detail level.
original: The input signal.
diagnostics: Quadrature, fallback, non-finite-integral, quantile-type, and MAD-level information.
settings: The principal fitting settings.

If return_wavelet = TRUE, the list also contains wavelet_object and thresholded_wavelet.

References

Sanyal, N. (2026). Resolution-Adaptive Compact-Support Priors for Bayesian Wavelet Denoising: A Wendland-Semicircle Slab Mixture for Low-SNR Signal Recovery. Axioms, 15(9), 678. <DOI:10.3390/axioms15090678>

Author

Nilotpal Sanyal

Examples

# \donttest{
set.seed(1)
n <- 128
x <- seq(0, 1, length.out = n)
y <- sin(6 * pi * x) + rnorm(n, sd = 0.5)

fit <- wswavelet(
  y,
  likelihood = "gaussian",
  filter.number = 6L,
  quadrature_n = 24L
)
fit$estimate
#>   [1] -0.071386952  0.105517006  0.281882808  0.456518416  0.613050068
#>   [6]  0.763681779  0.904313908  1.015850262  1.103742420  1.166099340
#>  [11]  1.214487222  1.239122970  1.236521005  1.186389290  1.119103299
#>  [16]  0.985804528  0.822982477  0.643915311  0.450070575  0.252023775
#>  [21]  0.058895198 -0.128441776 -0.298084132 -0.439455043 -0.533349635
#>  [26] -0.627812771 -0.704314551 -0.780392576 -0.838014732 -0.860555822
#>  [31] -0.852925694 -0.821255685 -0.758325791 -0.682517413 -0.593394681
#>  [36] -0.486688383 -0.376836129 -0.269047302 -0.173840212 -0.110127251
#>  [41] -0.063730924 -0.011039427  0.049511997  0.124583635  0.208469172
#>  [46]  0.286880683  0.366819982  0.449127838  0.526669310  0.605009665
#>  [51]  0.671887783  0.718521252  0.757675031  0.790745007  0.826555253
#>  [56]  0.848403985  0.852906121  0.846357570  0.828118479  0.783500824
#>  [61]  0.724042706  0.630895782  0.518065137  0.370356682  0.201317840
#>  [66]  0.026995247 -0.157277951 -0.315842070 -0.481337076 -0.632646877
#>  [71] -0.780317845 -0.906819713 -1.005251867 -1.088078016 -1.143256109
#>  [76] -1.172464232 -1.174937923 -1.132238979 -1.045170425 -0.908555214
#>  [81] -0.730202253 -0.529696625 -0.316974379 -0.118642830  0.092135258
#>  [86]  0.284003104  0.457451320  0.585058205  0.686252637  0.770669254
#>  [91]  0.842288654  0.920560883  0.983772803  1.006351234  0.996786974
#>  [96]  0.946806507  0.863880713  0.773276346  0.667828520  0.553913592
#> [101]  0.437399645  0.325209991  0.225468459  0.159442304  0.108297620
#> [106]  0.062416024  0.004037752 -0.069158110 -0.153250005 -0.221404659
#> [111] -0.297651789 -0.367571643 -0.429623169 -0.498453995 -0.554225418
#> [116] -0.593600630 -0.621349050 -0.646290775 -0.666426468 -0.688795831
#> [121] -0.704270188 -0.697255823 -0.681305185 -0.639787872 -0.576172469
#> [126] -0.489399576 -0.379048894 -0.235643167
fit$level_summary
#>   level coordinate coefficients spike_probability      beta        omega
#> 0     0  0.0000000            1         0.8105354 0.4566733 6.144175e-06
#> 1     1  0.1666667            2         0.9284007 0.4566733 8.315280e-07
#> 2     2  0.3333333            4         0.9641032 4.4558058 1.125352e-07
#> 3     3  0.5000000            8         0.9789878 1.3257960 1.522998e-08
#> 4     4  0.6666667           16         0.9864345 1.0412432 2.061154e-09
#> 5     5  0.8333333           32         0.9906295 0.9173786 2.789468e-10
#> 6     6  1.0000000           64         0.9931988 0.9969668 3.775135e-11
#>        mean_p0      mean_pW     mean_pS mean_log_m0_W mean_log_m0_S
#> 0 8.279244e-01 1.148495e-06 0.172074482    -0.1699530    -0.2527230
#> 1 9.329233e-01 5.862605e-08 0.067076661    -0.3707288    -0.4208667
#> 2 9.437820e-15 2.347335e-09 0.999999998    -7.0661761    -3.1011819
#> 3 9.379320e-01 5.164522e-10 0.062068014    -1.4312545    -1.2031533
#> 4 9.848497e-01 2.928933e-11 0.015150343    -0.7787522    -0.8478558
#> 5 9.921376e-01 2.486652e-12 0.007862381    -0.5105236    -0.6621466
#> 6 9.936713e-01 2.534684e-13 0.006328689    -0.6418063    -0.7602679
#>   quadrature_fallbacks
#> 0                    0
#> 1                    0
#> 2                    0
#> 3                    0
#> 4                    0
#> 5                    0
#> 6                    0
# }